R混合效应模型入门：lme4 处理重复测量与嵌套数据

前面几节的所有模型都有一个共同前提：观测之间相互独立。重复测量和嵌套数据直接违反这一点，而这类设计在实验研究里普遍存在——同一只动物测多次、同一个班级的学生、同一批次的样本。

578 行不等于 578 个独立样本

ChickWeight 记录 50 只小鸡在不同天数（第 0 到 21 天）的体重，共 578 条观测。

str(ChickWeight)

Classes ‘nfnGroupedData’, ‘nfGroupedData’, ‘groupedData’ and 'data.frame':  578 obs. of  4 variables:
 $ weight: num  42 51 59 64 76 93 106 125 149 171 ...
 $ Time  : num  0 2 4 6 8 10 12 14 16 18 ...
 $ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
 $ Diet  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...

输出里的 Ord.factor 和 nfnGroupedData 这些额外类来自 ChickWeight 自带的分组数据属性，不影响 lmer() 的使用，只影响 str() 和自动绘图的表现。

开头 12 行全部来自第 1 号小鸡：体重从 42 克一路长到 205 克，是一条平滑的个体生长曲线，相邻两次测量几乎完全可以由前一次推测出来。这 12 条记录携带的信息远少于 12 个独立样本，把 578 行当成 578 份独立证据，等于高估了数据里的信息量。

普通 lm() 错在哪

summary(lm(weight ~ Time, data = ChickWeight))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  27.4674     3.0365   9.046   <2e-16 ***
Time          8.8030     0.2397  36.725   <2e-16 ***
---
Residual standard error: 38.91 on 576 degrees of freedom
Multiple R-squared:  0.7007,  Adjusted R-squared:  0.7002
F-statistic:  1349 on 1 and 576 DF,  p-value: < 2.2e-16

问题出在 Residual standard error: 38.91。这个数字是「小鸡之间的差异」和「同一只小鸡不同天之间的差异」混在一起的结果，哪一个都代表不了。模型只估计了一个误差项，却要拿它去检验所有系数，于是：

• 对个体内变量（Time），标准误偏保守。混合模型给出的 Time 标准误是 0.1755，比这里的 0.2397 小 27%，因为个体内的变异其实没那么大。
• 对个体间变量（比如 Diet，每只小鸡只吃一种饲料），标准误明显偏小。用 lm(weight ~ Diet) 时 Diet2 的标准误是 7.867，混合模型给的是 10.24——普通 lm 把标准误低估了约 30%，p 值因此偏小，假阳性风险直接上升。
• 数据不平衡时（ChickWeight 有 5 只小鸡提前退出）系数本身也会偏。Time 的斜率在 lm 里是 8.803，在混合模型里是 8.726。

三个问题方向不同，但根源相同：独立性的假设不成立，误差结构被写错了。

随机截距：让每只小鸡有自己的起点

lme4 包的 lmer() 用公式里的括号项指定随机效应：

library(lme4)
m1 <- lmer(weight ~ Time + (1 | Chick), data = ChickWeight)
summary(m1)

Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ Time + (1 | Chick)
   Data: ChickWeight

REML criterion at convergence: 5619.4

Scaled residuals:
    Min      1Q  Median      3Q     Max
-3.0086 -0.5528 -0.0888  0.4975  3.4588

Random effects:
 Groups   Name        Variance Std.Dev.
 Chick    (Intercept) 717.9    26.79
 Residual             799.4    28.27
Number of obs: 578, groups:  Chick, 50

Fixed effects:
            Estimate Std. Error t value
(Intercept)  27.8451     4.3877   6.346
Time          8.7261     0.1755  49.716

Correlation of Fixed Effects:
     (Intr)
Time -0.422

(1 | Chick) 读作「为每只小鸡估计一个随机截距」。它的含义是：允许每只小鸡有自己的起始体重，这些起点围绕总体均值波动，波动幅度由方差参数描述。1 表示只有截距随个体变化，斜率对所有人相同。

读 summary：两部分

Random effects 部分有两行：

• Chick (Intercept) 717.9 / 26.79：小鸡之间起始体重的方差 717.9（标准差 26.79 克）。
• Residual 799.4 / 28.27：同一只小鸡在自己生长曲线周围的波动，即个体内残差。

模型的因变量总方差可以粗略看成 717.9 + 799.4 = 1517.3，其中 717.9 / 1517.3 = 47.3% 来自个体差异。这个比例叫组内相关系数（ICC，intraclass correlation coefficient），0.47 意味着同一只小鸡的两次测量之间的相关性很高——这正是不能用普通 lm() 的量化证据。

Fixed effects 部分和 lm() 的输出格式相同，有 Estimate、Std. Error、t value，但没有 Pr(>|t|) 这一列。这不是输出被截断了，是 lme4 作者有意为之：混合模型里固定效应的分母自由度没有唯一定义（它取决于随机效应的方差估计），硬算出来的 p 值会误导人。

需要 p 值有三个选择：

1. 用 lmerTest 包的 lmer() 替代，它自动做 Satterthwaite 自由度近似：

library(lmerTest)
summary(lmer(weight ~ Time + (1 | Chick), data = ChickWeight))

2. 用似然比检验（likelihood ratio test，LRT）比较嵌套模型：

m1_ml <- lmer(weight ~ Time + (1 | Chick), data = ChickWeight, REML = FALSE)
m0_ml <- lm(weight ~ Time, data = ChickWeight)
anova(m1_ml, m0_ml)

Data: ChickWeight
Models:
m0ml: weight ~ Time
m1: weight ~ Time + (1 | Chick)
     npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
m0ml    3 5876.8 5889.9 -2935.4    5870.8
m1      4 5630.3 5647.8 -2811.2    5622.3 248.45  1  < 2.2e-16 ***

加入随机截距让模型显著变好（χ² = 248.45，p < 2.2e-16）。注意这里必须写 REML = FALSE：不是因为随机效应本身不能做似然比检验（两个模型固定效应相同时，用 REML 比较随机效应是可以的），而是因为 m0_ml 是 lm 对象，它的对数似然按 ML 计算——REML 和 ML 的似然不能混着比。检验方差参数是否为零时，原假设还落在参数空间的边界上，χ² 分布给出的 p 值偏保守，实际 p 值比输出的更小——所以这里的「显著」是稳妥的结论。

第三个选择是用 pbkrtest 包的 KRmodcomp() 做 Kenward-Roger 校正，它的近似在小样本下更精确，代价是计算慢得多。

日常报告里，lmerTest 的 Satterthwaite 近似最省事，审稿人也认。

随机斜率：(1 + Time | Chick)

随机截距假定所有小鸡生长速度相同。但从生物学上想，有的小鸡长得快、有的长得慢，斜率也应该允许变化：

m2 <- lmer(weight ~ Time + (1 + Time | Chick), data = ChickWeight)
summary(m2)

Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ Time + (1 + Time | Chick)
   Data: ChickWeight

REML criterion at convergence: 4827.5

Scaled residuals:
    Min      1Q  Median      3Q     Max
-2.6731 -0.5563 -0.0277  0.5010  3.4939

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Chick    (Intercept) 140.54   11.855
          Time         14.14    3.761   -0.95
 Residual             163.51   12.787
Number of obs: 578, groups:  Chick, 50

Fixed effects:
            Estimate Std. Error t value
(Intercept)  29.1780     1.9573   14.91
Time          8.4531     0.5408   15.63

Correlation of Fixed Effects:
     (Intr)
Time -0.871

这里 Random effects 多了两样东西。Time 那一行是随机斜率的方差（14.14，标准差 3.761 克/天），说明小鸡的生长速度确实有差异。Corr 列的 −0.95 是随机截距和随机斜率之间的相关系数，含义是「起始体重越低的小鸡长得越快」——负相关很强。这条信息本身往往就有生物学意义，别忽略。

模型比较：两个模型的固定效应完全相同，可以直接比 AIC。

AIC(m1, m2)

   df      AIC
m1  4 5627.398
m2  6 4839.499

AIC 降了将近 800，随机斜率值得加。

符号速记：(1 | g) 随机截距；(1 + x | g) 随机截距加随机斜率，且估计两者的相关；(0 + x | g) 只要随机斜率；(1 + x || g)（双竖线）表示截距和斜率独立、不估计相关。

提取结果

fixef(m1)
ranef(m1)$Chick[1:3, , drop = FALSE]
VarCorr(m1)
sigma(m1)

(Intercept)        Time
  27.844165    8.726255
   (Intercept)
18   0.2739469
16 -26.2298297
15 -25.2208687
 Groups   Name        Std.Dev.
 Chick    (Intercept) 26.793
 Residual             28.274
[1] 28.24714

fixef() 给固定效应（总体平均的截距和斜率），ranef() 给每只小鸡偏离总体均值多少——第 16 号小鸡比平均轻 26.2 克。ranef() 的结果可以直接拿去做后续分析，比如挑出生长异常的个体。sigma() 是残差标准差。

边界奇异拟合警告

如果随机结构超出了数据能支撑的范围，lmer() 会给出这个警告：

lmer(weight ~ Time + (1 + Time | Diet), data = ChickWeight)

boundary (singular) fit: see help('isSingular')

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Diet     (Intercept)   14.913  3.862
          Time           3.499  1.871   -1.00
 Residual             1157.779 34.026
Number of obs: 578, groups:  Diet, 4

Diet 只有 4 个水平，却要估计两个方差加一个相关系数。结果 Corr 顶到了 −1.00 这个理论边界上——这不是「发现了完全负相关」，而是模型在告诉你「这个参数我估不出来」。奇异拟合（singular fit）意味着某个方差成分被估成了 0，或者相关系数顶到了 ±1，模型处于参数空间的边界，此时的方差估计和标准误都不可信。

处理办法有三个，按推荐顺序：

1. 简化随机结构。这里正确的做法是用 (1 | Chick)，因为 Diet 是组间的处理变量，本来就不该做随机斜率。
2. 如果确实需要两者但相关性估不出来，去掉相关：写成 (1 + Time || Chick)。
3. 检查是不是分组水平太少。少于 5 到 6 个水平的分组变量，通常不该放随机效应——方差参数本身估不准。

看到警告不要直接忽略，也不要不看结果就照抄。用 isSingular(m1) 可以程序化判断，方便在批量分析里做检查。

固定效应还是随机效应

判断标准不是「哪个变量」，而是研究问题：

• 固定效应：你关心的、想估计其效应大小并做推断的水平。比如三种饲料处理之间的差异，结论要推广到「这三种饲料」以外，水平是精心选择的。
• 随机效应：你只想控制其变异、不关心具体每个水平。比如 50 只小鸡，没人关心第 16 号小鸡比平均轻多少，只希望个体差异不要污染对 Time 的估计。

同一个变量在不同研究里可以是固定效应也可以是随机效应。把分组变量当固定效应处理（lm(weight ~ Time + Chick)）会消耗 49 个自由度，估计出一堆没人看的系数；当随机效应处理则只花 1 个方差参数。反过来，把只有 3 个水平的处理变量当随机效应，方差估计极不可靠，也是常见错误。

下一步

混合模型是线性模型的扩展：固定效应部分怎么写、系数怎么解读，和 /r/modeling/linear-regression 完全一致；当因变量是二分类时，把 lmer() 换成 glmer(..., family = binomial)，随机效应的写法不变。固定效应与随机效应的概念区别，可以先回看 /r/modeling/linear-regression 里对固定效应的处理，再看本篇怎么把随机效应加进去。