R语言回归建模全流程：从lm到混合效应模型与GAM

之前在处理一批生态学调查数据时，我反复卡在同一个问题上：数据有明显的分组结构，而普通的线性回归默认所有观测相互独立。网上教程要么只讲 lm，要么直接上混合效应模型，却很少有人把从数据清洗、基础回归、广义线性模型、混合效应模型、时间/空间/系统发育结构到 GAM、再到结果绘图输出的整套流程串起来。这次我把整个学习路径整理成六个单元，并配好了可运行的 R 代码。无论是生态学、医学、教育学还是社科领域的从业者，只要你的数据存在分组、重复测量、空间采样或物种亲缘关系，这篇文章都能帮你建立一套完整的回归建模思路。

本文读者画像：会用一点 R，但没系统建过模型；或者已经用 lm 做过回归，但面对嵌套数据、非线性关系时不知道下一步怎么走。学完本文，你可以独立完成“数据清洗 → 模型选择 → 模型拟合 → 模型诊断 → 结果可视化”的一整套流程，并理解 lm、glm、lmm、glmm、GAM 这几类模型之间的区别与联系。

1. 背景与核心概念

1.1 回归分析到底在解决什么问题

回归分析是研究响应变量（也叫因变量、被解释变量）与一个或多个解释变量（也叫自变量、协变量）之间关系的统计方法集合。最简单的线性回归用一条直线描述两者关系：

y = β0 + β1x + ε

其中 β0 是截距，β1 是斜率，ε 是随机误差。这个公式看起来简单，但实际数据分析中真正困难的部分不是拟合公式，而是判断数据是否满足模型假设：独立性、正态性、方差齐性、线性关系。

普通线性回归（lm）适用于独立且同分布的数据，但真实数据往往不是这样。例如：同一所学校里多个学生的成绩彼此相关，同一个样地多次采集的土壤数据彼此相关，同一物种不同个体因为亲缘关系而相似。这些数据的“非独立性”会直接导致普通回归的标准误被低估，从而把不显著的效应误判为显著。

1.2 什么是混合效应模型

混合效应模型（Mixed Effects Model）在固定效应之外引入随机效应（Random Effects），用来刻画组内相关性。固定效应是我们关心的、需要估计和解释的变量；随机效应则描述数据的分层结构或重复测量结构带来的随机波动。

一个最简单的随机截距模型可以写成：

y_ij = β0 + β1x_ij + u_i + ε_ij

这里 i 表示第 i 个组，j 表示组内第 j 个观测。u_i 是每个组自己的随机截距，通常假设服从均值为 0、方差为 σ²_u 的正态分布。这样，同一组内的观测共享同一个 u_i，天然产生了组内相关性。

混合效应模型的优势在于：既能估计我们关心的固定效应，又能把组间差异当作方差来源来处理，而不是逐组单独建模或把组别当成普通分类变量。这样既保留了一般性结论，又避免了假重复（Pseudo-replication）问题。

1.3 为什么需要完整学习六大单元

很多初学者直接从混合效应模型或者 GAM 开始学，结果一遇到报错就卡住。真正高效的路径是：先掌握 R 语言数据操作，再学会基础回归，理解模型假设，然后引入随机效应，最后扩展到非线性平滑项和特殊数据结构。六大单元正好对应这条路径：

1. R 语言基础与数据操作
2. 线性回归与广义线性回归（lm/glm）
3. 混合效应模型（lmm/glmm）
4. 时间、空间与系统发育分析
5. 广义加性模型（GAM）
6. 结果可视化与出图

每个单元不是孤立的。lmm 是 lm 的自然延伸，GAM 又可以看作 glm 的平滑化扩展。理解了它们之间的关系，你在面对新数据时才能快速确定合适的模型。

2. 环境准备与版本说明

2.1 R 与 RStudio 环境

R 语言支持 Windows、macOS 和 Linux。建议安装 R 4.1 以上版本，并配合 RStudio 使用。RStudio 不是必须的，但它自带脚本编辑器、变量查看器、绘图窗口和包管理界面，对模型调试非常友好。

关于系统版本，不需要追求最新，稳定即可。关键是把整个分析放在一个可复现的环境中。建议用项目文件夹管理。

2.2 需要安装的扩展包

本文会用到以下 R 包：

• tidyverse / dplyr / ggplot2：数据清洗与绘图
• readxl / readr：数据导入
• lme4 / glmmTMB / nlme：混合效应模型
• mgcv / gratia / ggeffects：GAM 与模型可视化
• performance / DHARMa：模型诊断
• forecast：时间序列
• spdep：空间统计
• ape / phytools / phylolm：系统发育分析
• emmeans：边际均值比较
• car：共线性诊断等工具

安装命令很简单：

install.packages(c("tidyverse", "readxl", "lme4", "glmmTMB", "nlme", "mgcv", "gratia", "ggeffects", "performance", "DHARMa", "forecast", "spdep", "ape", "phytools", "phylolm", "emmeans", "car"))

如果你在国内网络环境下安装较慢，可以设置镜像：

options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

2.3 验证环境是否正常

安装完成后，加载核心包并查看版本：

library(lme4) library(mgcv) library(ggplot2) packageVersion("lme4") packageVersion("mgcv")

如果输出类似 “1.1-35.1”、“1.9-1” 的版本号，说明环境基本就绪。后续不同机器上的包版本可能有差异，但本文的代码使用的是这些包的基础接口，兼容性较好。

3. 第一单元：R 语言基础与数据操作

3.1 R 语言核心对象

R 中最常用的数据对象是向量、因子、数据框和列表。向量是最基础的单位：

# 数值向量 height <- c(1.62, 1.70, 1.75, 1.68, 1.73) # 字符向量 site <- c("A", "B", "A", "C", "B") # 因子：分类变量在建模中的标准形式 site <- factor(site, levels = c("A", "B", "C")) # 数据框：建模中最常用的二维数据结构 df <- data.frame(height = height, site = site)

因子变量在 R 中非常重要。建模时，字符型的列可能会被误判为连续变量；使用 factor() 显式声明分类变量的水平顺序，可以避免截距水平混乱的问题。

3.2 数据导入与清洗

实际分析的第一步一定是读入数据。CSV 是最常见的格式：

# 读取 CSV data <- read.csv("mydata.csv", stringsAsFactors = FALSE, fileEncoding = "UTF-8") # 读取 Excel library(readxl) data <- read_excel("mydata.xlsx", sheet = 1)

读取之后，建议先看结构和缺失值：

str(data) # 查看每列类型 head(data) # 查看前 6 行 summary(data) # 描述统计

dplyr 是清洗数据的利器。常见的操作包括筛选、选择列、生成新变量、分组汇总：

library(dplyr) data_clean <- data %>% filter(!is.na(weight)) %>% # 去掉缺失体重 select(site, treatment, weight, height) %>% # 只保留所需列 mutate(bmi = weight / height^2) %>% # 生成新变量 arrange(desc(bmi)) # 按 BMI 排序

这里要特别注意：filter(!is.na(weight))会把 weight 为空的记录删除。如果缺失值比例很高，删除可能引入偏差，需要结合业务判断处理方式。

3.3 分组汇总与快速可视化

描述统计是建模前必须做的一步。分组均值、标准差和样本量能为你提供最基本的分布信息：

data_clean %>% group_by(treatment) %>% summarise( mean_bmi = mean(bmi, na.rm = TRUE), sd_bmi = sd(bmi, na.rm = TRUE), n = n() )

快速可视化可以先用 ggplot2：

library(ggplot2) ggplot(data_clean, aes(x = treatment, y = bmi, fill = treatment)) + geom_boxplot() + theme_bw() + labs(title = "不同处理组的 BMI 分布")

可视化不是为了发文章，而是帮你发现异常值和分布形状。如果某个处理组只有 3 个样本，后面建模型时要格外谨慎。

4. 第二单元：线性回归与广义线性回归

4.1 线性回归 lm 的核心用法

线性回归对应的核心函数是lm()。以 R 内置的 mtcars 数据为例，分析油耗 mpg 与重量 wt、马力 hp 的关系：

model_lm <- lm(mpg ~ wt + hp, data = mtcars) summary(model_lm)

输出中最值得关注的是：

• Estimate：回归系数，表示在其他变量不变时，该变量每增加一个单位，响应变量的平均变化量。
• Std. Error：系数的标准误。
• Pr(>|t|)：显著性检验的 p 值。
• R-squared / Adjusted R-squared：模型解释了多少变异。
• F-statistic：整个模型的整体显著性。

运行后可以看到，wt 的系数大约为 -3.87，hp 的系数大约为 -0.03，说明车重对油耗的影响远大于马力。

4.2 回归诊断与共线性检查

拟合完模型不能直接下结论。线性回归有几条关键假设需要验证：残差正态性、方差齐性、无强影响点、解释变量之间无严重共线性。

基础诊断图可以直接用 plot()：

plot(model_lm)

这会生成四张图：残差与拟合值图、QQ 图、尺度位置图和残差与杠杆图。你不需要每张都看懂，但至少要关注残差图中是否存在明显的喇叭形或弯曲趋势。

更系统的检查可以用 performance 包：

library(performance) check_model(model_lm)

共线性问题可以用方差膨胀因子（VIF）检查。VIF 超过 5 或 10 通常说明存在较强共线性：

library(car) vif(model_lm)

如果 VIF 过高，可以考虑删除其中一个变量，或使用正则化方法。

4.3 逻辑回归 glm 实战

当响应变量是二分类（0/1）时，线性回归不再适用。此时使用广义线性模型（GLM），通过链接函数把响应变量与线性预测项连接起来。最常用的是逻辑回归：

model_glm <- glm(am ~ mpg + wt, data = mtcars, family = binomial) summary(model_glm)

这里 am 表示汽车变速箱类型（0 = 自动，1 = 手动），属于典型的二分类变量。family = binomial 表示我们假设 am 服从二项分布，并使用 logit 链接函数。

逻辑回归的系数表示 log-odds（对数优势比）的变化。解释时通常会取指数得到优势比：

exp(coef(model_glm))

比如 mpg 的系数为 1.26，说明 mpg 每增加一个单位，车子是手动挡的优势比变为原来的 exp(1.26) ≈ 3.5 倍。

4.4 泊松回归处理计数数据

如果响应变量是非负整数（例如样方中的物种个体数、医院接诊人数），通常使用泊松回归：

count_data <- data.frame( x = rnorm(100, 10, 2), count = rpois(100, lambda = 5) ) model_pois <- glm(count ~ x, data = count_data, family = poisson) summary(model_pois)

泊松回归假设方差等于均值，但实际计数数据常有过度离散（方差大于均值）。如果发现残差偏差与自由度的比例远大于 1，可以考虑负二项回归。R 中可以用MASS::glm.nb()：

library(MASS) model_nb <- glm.nb(count ~ x, data = count_data) summary(model_nb)

4.5 什么时候应该进入混合效应模型

如果你在诊断 lm/glm 时发现残差仍然存在明显的组内聚集，比如同一个地点、同一个个体、同一批样本之间存在相关性，说明数据不是完全独立的，此时就应该进入第三单元，考虑混合效应模型。

判断是否非独立，可以从研究设计入手：数据有没有分组变量？是否对同一对象重复测量？样方是否嵌套在样地中？只要答案是肯定的，你就需要在模型中加入随机效应。

5. 第三单元：混合效应模型

5.1 固定效应与随机效应

混合效应模型的关键是区分固定效应和随机效应。

固定效应指的是我们关心的、希望估计并报告的解释变量，例如不同处理之间的差异。随机效应则代表抽样来源或分组结构，例如学校、样地、个体、年份。随机效应的作用不是得到具体某个组的效应值，而是估计组间方差、刻画组内相关性。

最常见的随机效应写法：

(1 | group) # 随机截距：每个组有不同的基线水平 (1 | group1 / group2) # 嵌套结构：group2 嵌套在 group1 内 (1 + x | group) # 随机斜率：每个组对 x 的响应不同

5.2 随机截距模型 lmer

模拟一份学校教学实验数据。20 所学校，每所 30 名学生，三种教学法 A、B、C：

set.seed(123) school_id <- rep(paste0("S", 1:20), each = 30) teaching <- rep(c("A", "B", "C"), each = 10, times = 20) school_effect <- rep(rnorm(20, 0, 5), each = 30) score <- 60 + ifelse(teaching == "A", 8, ifelse(teaching == "B", 3, 0)) + school_effect + rnorm(600, 0, 8) students <- data.frame(school_id, teaching, score)

将 school_id 转为因子，然后用 lme4 拟合随机截距模型：

library(lme4) students$school_id <- factor(students$school_id) students$teaching <- factor(students$teaching, levels = c("A", "B", "C")) model_lmm <- lmer(score ~ teaching + (1 | school_id), data = students) summary(model_lmm)

summary 输出里有两大部分。固定效应部分显示教学法 A、B、C 之间的差异；随机效应部分给出 school_id 的标准差，大约在 5 左右，这正好对应我们模拟的学校间标准差。如果忽略学校分层直接使用 lm，学校效应会被当作误差的一部分，导致固定效应的标准误偏小。

5.3 随机斜率模型

不同学校对教学法的响应可能不同，也就是说教学法斜率存在校际差异。此时可以加随机斜率：

model_lmm_slope <- lmer(score ~ teaching + (1 + teaching | school_id), data = students) summary(model_lmm_slope)

随机斜率模型的参数更多，也更容易出现收敛问题。一个常见建议是：当随机斜率模型的方差分量很小、几乎为 0 时，优先选择更简单的随机截距模型。统计上叫简约原则。

5.4 广义混合效应模型 glmer

当响应变量是二分类、计数或非正态数据，且又存在分组结构时，需要使用广义线性混合效应模型（GLMM）。lme4 中的对应函数是glmer()。

在上面的模拟数据里生成一个二分类变量 pass：

students$pass <- ifelse(score > 65, 1, 0) model_glmm <- glmer(pass ~ teaching + (1 | school_id), data = students, family = binomial) summary(model_glmm)

glmer 的 family 参数与 glm 一致，支持 binomial、poisson 等。输出中的系数仍是 log-odds。解释方式与 glm 完全一致，只是标准误会因为随机效应的存在而更大、更可信。

如果数据量大、随机效应结构复杂，lme4 可能比较慢，可以尝试 glmmTMB：

library(glmmTMB) model_glmm2 <- glmmTMB(pass ~ teaching + (1 | school_id), data = students, family = binomial) summary(model_glmm2)

5.5 模型比较与预测

嵌套模型之间的比较可以使用似然比检验：

model_lmm_simple <- lmer(score ~ teaching + (1 | school_id), data = students) model_lmm_full <- lmer(score ~ teaching + (1 + teaching | school_id), data = students) anova(model_lmm_simple, model_lmm_full)

p 值不显著说明复杂模型没有显著改善拟合，可以继续使用简单模型。非嵌套模型则使用 AIC 比较，AIC 越小越好。

预测时如果不希望包括随机效应，使用re.form = NA：

students$pred_fixed <- predict(model_lmm, re.form = NA) students$pred_full <- predict(model_lmm, re.form = NULL)

固定效应预测值用于绘制边际效应图，完整预测值则用于查看模型对原始数据的拟合程度。

6. 第四单元：时间、空间与系统发育分析

6.1 时间序列数据

时间序列数据的核心特征是相邻观测之间存在自相关。R 内置数据 AirPassengers 是经典的月度客运量数据，长度为 144，有明显的季节趋势。首先把数据转成 ts 对象，再拟合 ARIMA：

library(forecast) data_ts <- AirPassengers str(data_ts) fit_arima <- auto.arima(data_ts, seasonal = TRUE) summary(fit_arima)

预测未来 12 个月：

forecast_result <- forecast(fit_arima, h = 12) plot(forecast_result)

如果你更希望在混合效应模型框架中处理时间相关性，可以使用 nlme 包的相关结构：

library(nlme) model_gls <- gls(y ~ time + x, data = df, correlation = corAR1(form = ~ time | id)) summary(model_gls)

corAR1 表示一阶自回归相关结构，适合等时间间隔的重复测量数据。这里的 y、time、x、id 需要根据你自己的数据定义。

6.2 空间自相关与空间回归

如果样本点在地理空间上分布，距离近的样点可能更相似，这就是空间自相关。忽略空间自相关同样会使标准误偏小。

最简单的检查方法是计算残差的 Moran's I：

library(spdep) set.seed(123) coords <- cbind(runif(50, 0, 10), runif(50, 0, 10)) y <- 2 + 0.5 * coords[, 1] + rnorm(50, 0, 1) x <- 0.3 * coords[, 2] + rnorm(50, 0, 1) model_sp <- lm(y ~ x) nb <- knn2nb(knearneigh(coords, k = 4)) lw <- nb2listw(nb) moran.test(resid(model_sp), lw)

如果 Moran's I 检验显著，说明残差存在空间结构。此时可以用空间回归模型：

model_sar <- lagsarlm(y ~ x, data = data.frame(y, x), listw = lw) summary(model_sar)

在混合效应模型中加入空间相关结构也是常见做法。glmmTMB 支持多种空间协方差结构，不过设置和收敛诊断比较复杂，建议先从简单结构开始，并用性能检查确认改进效果。

6.3 系统发育数据分析

物种数据往往共享进化历史，亲缘关系近的物种表型可能更相似。系统发育广义最小二乘（PGLS）是处理此类非独立性的常用方法。核心思想是依据系统发育树构建物种间的协方差矩阵。

以 phylolm 包为例：

library(ape) library(phylolm) set.seed(123) tree <- rtree(50) trait_x <- rTraitCont(tree) trait_y <- 2 + 1.2 * trait_x + rTraitCont(tree) fit_phylo <- phylolm(trait_y ~ trait_x, phy = tree, model = "BM") summary(fit_phylo)

model = "BM" 表示布朗运动模型，是最常用的系统发育相关结构。实际研究中，trait_x 和 trait_y 是你自己测定的物种性状，tree 则来自分子系统学或公开发表的物种树。分析前还要检查树是否是二叉、是否包含所有物种。

7. 第五单元：广义加性模型（GAM）

7.1 GAM 的基本思想

广义加性模型（Generalized Additive Model，GAM）是广义线性模型的平滑化扩展。它允许响应变量与解释变量之间不是直线关系，而是通过一组平滑函数来拟合。

核心公式：

y = β0 + f1(x1) + f2(x2) + ...

f 是平滑函数，由数据自动弯曲，因此你不需要先验地指定是二次函数还是三次函数。这是 GAM 最大的优点，尤其在趋势未知、关系复杂的探索性分析中非常实用。

7.2 mgcv 包基础示例

mgcv 是 R 中最成熟的 GAM 包。看一个例子：分析 mpg 与 wt 的关系，同时调整 hp：

library(mgcv) model_gam <- gam(mpg ~ s(wt) + hp, data = mtcars) summary(model_gam)

summary 中注意两类信息：

• 参数项部分：hp 的系数和 p 值。
• 平滑项部分：s(wt) 的有效自由度（edf）。edf 接近 1 表示近似线性，edf 越大表示曲线弯曲程度越高。

绘制平滑项：

plot(model_gam, pages = 1)

7.3 交互平滑与广义 GAM

如果两个变量共同影响响应，且这种影响是非线性的，可以使用张量积平滑：

model_gam_te <- gam(mpg ~ te(wt, hp), data = mtcars) summary(model_gam_te)

te() 会拟合一个二维交互曲面，适合两个变量联合效应。如果只想加入一个主效应的交互，比如其中一个变量线性、另一个平滑，可以写成：

model_gam_interact <- gam(mpg ~ s(wt) + hp + s(wt, by = hp), data = mtcars)

GAM 同样支持非正态响应：

model_gam_binomial <- gam(am ~ s(mpg) + wt, data = mtcars, family = binomial) summary(model_gam_binomial)

7.4 GAM 诊断与调参

拟合完 GAM 后，使用 gam.check() 检查基函数数量和残差：

gam.check(model_gam)

gam.check 会输出 k-index，如果 k-index 明显小于 1，说明默认的 k 值太小，需要增大：

model_gam_k <- gam(mpg ~ s(wt, k = 15), data = mtcars)

k 不是越大越好，过大的 k 容易导致过度拟合。实际使用时可以结合有效自由度判断：如果 edf 明显低于 k，说明 k 已经足够；如果 edf 接近 k 上限，说明需要增大 k。

整体模型评估可以用 AIC：

AIC(model_lm, model_gam)

如果 GAM 的 AIC 明显小于线性模型，说明数据中确实存在非线性关系。

8. 第六单元：结果可视化与出图

8.1 ggplot2 基础图形

结果可视化的核心工作是把模型结论翻译成读者容易理解的图形。ggplot2 是 R 中最常用的绘图包。先看基础散点与回归线：

library(ggplot2) p <- ggplot(mtcars, aes(x = wt, y = mpg)) + geom_point(size = 3, alpha = 0.7) + geom_smooth(method = "lm", se = TRUE) + theme_bw(base_size = 14) + labs(x = "Weight (1000 lbs)", y = "Miles per gallon") print(p)

8.2 绘制模型预测结果

展示回归模型时，通常绘制预测值和置信区间，而不是只画原始数据。这可以使用 ggeffects 包：

library(ggeffects) pred_lm <- ggpredict(model_lm, terms = "wt") plot(pred_lm) + labs(title = "Linear model prediction")

对于混合效应模型，预测时默认包含随机效应。如果你想展示固定效应的边际效应，需要在 ggpredict 中指定：

pred_lmm <- ggpredict(model_lmm, terms = "teaching") plot(pred_lmm)

这幅图会展示三种教学法的预测均值和置信区间，是论文和报告中非常有用的图。

8.3 绘制随机效应与平滑项

随机截距的分布可以用 dotplot 展示：

library(lme4) dotplot(ranef(model_lmm, condVar = TRUE))

每所学校对应一个点加误差线，可以直观看到哪些学校显著高于或低于平均水平。

GAM 的平滑项可以用 mgcv 自带函数绘制：

plot(model_gam, pages = 1, shade = TRUE)

更美观的方式是使用 gratia 包：

library(gratia) draw(model_gam)

draw 会自动绘制所有平滑项，并附带置信区间，输出风格更适合直接用于报告。

8.4 输出高清图片

RStudio 的绘图窗口适合交互查看，但要用于论文或汇报，需要保存为高清图片：

dir.create("figures", showWarnings = FALSE) p <- ggplot(mtcars, aes(x = wt, y = mpg)) + geom_point(size = 3) + geom_smooth(method = "lm", se = TRUE) + theme_bw(base_size = 14) ggsave("figures/scatter_mpg_wt.png", p, width = 8, height = 6, dpi = 300)

dpi = 300 可以满足大部分期刊要求。保存 PDF 矢量图则适合需要文字编辑的场景：

ggsave("figures/scatter_mpg_wt.pdf", p, width = 8, height = 6)

9. 常见问题与排查思路

9.1 高频问题速查表

9.2 三个典型问题的排查过程

第一个是 lmer 收敛警告。一个常见处理办法是：

model_lmm <- lmer(score ~ teaching + (1 | school_id), data = students, control = lmerControl(optCtrl = list(maxfun = 100000)))

如果这样仍不收敛，通常意味着随机效应结构过于复杂，而不是迭代次数不够。

第二个是 glmer 提示 singularity。这表示模型估计出随机效应方差接近 0，说明该随机效应没有起到应有的作用。此时最简单的做法是去掉这个随机效应，重新拟合后再比较 AIC。

第三个是中文乱码。建议优先把所有数据文件统一保存为 UTF-8 编码，并在读入时明确指定：

data <- read.csv("data.csv", fileEncoding = "UTF-8")

如果文件是从 Excel 另存的 CSV，注意 Excel 默认可能在 Windows 下保存为 GBK，此时需要尝试 fileEncoding="GBK"。

10. 最佳实践与学习路线

10.1 建模工程化流程

一个规范的 R 分析项目应该有清晰的流程和文件结构。我建议每个项目按以下步骤运行：

1. 数据读取后先看缺失值、异常值和变量类型。
2. 用直方图或箱线图检查响应变量分布。
3. 根据研究问题选择模型家族：连续变量用高斯、0/1 用二项、计数用泊松或负二项。
4. 先拟合基线模型，再逐步添加随机效应和复杂结构。
5. 每次修改模型后记录 AIC、系数估计和诊断结果。
6. 用 sessionInfo() 保存环境信息，保证结果可复现。

一个小建议：把数据清洗、模型拟合、结果可视化分成三个脚本文件，这样当数据更新时只需要重新运行第一个脚本，模型和图形自动更新。

10.2 模型报告建议

写论文或技术报告时，建议至少报告以下内容：

• 数据类型和样本量。
• 模型公式。
• 固定效应的估计值、标准误、置信区间和 p 值。
• 随机效应的方差分量。
• 模型总体的 AIC 或 R²。
• R 及关键扩展包的版本。

只有报告了这些信息，别人才可能复现你的分析，这也是统计可复现性的基本要求。

10.3 给初学者的一条建议

如果你刚接触这套流程，不要急着把所有模型都学会。先拿一份自带分组结构的模拟数据，分别用 lm、lmer 拟合同样的固定效应，然后对比两种模型输出的标准误。当你亲眼看到 lm 的标准误明显偏小时，你就真正理解了混合效应模型的价值。

之后再进入 GAM、空间分析和系统发育分析时，你面对的不是一堆孤立的新方法，而是一套统一思想下的不同工具。数据有非线性关系，用 GAM；数据有分组结构，用混合效应模型；数据有空间或进化相关，用相关结构或系统发育模型。把每个工具放到正确的位置，你的分析能力会上一个台阶。