6 单随机因子线性混合效应模型
6.1 什么时候，为什么，你想用线性混合效应建模取代传统的分析?
我们已经反复强调了，心理学中的许多常用技术可以被视为一般线性模型的特例。这意味着可以用回归来替代这些技术。事实上，你可以使用以下四个函数之一来分析几乎任何在心理学中可以想象得到的数据集。
| 取样设计 | 数据类型 | 函数 | 描述 |
|---|---|---|---|
| 单层 | 连续、正态分布 | base::lm() |
简单线性模型 |
| 单层 | 计数或分类 | base::glm() |
广义线性模型 |
| 多层 | 连续、正态分布 | lme4::lmer() |
线性混合效应模型 |
| 多层 | 计数或分类 | lme4::glmer() |
广义线性混合效应模型 |
在本章结束时，你应该:
- 理解如何使用回归来替代标准的ANOVA和t检验分析，适用于具有单一随机因子和连续因变量的数据；
- 能够进行模型比较(
anova())以测试效果； - 能够使用R回归公式语法来表达各种研究设计。
要决定使用四个回归函数中的哪一个，你需要回答两个问题。
- 因变量代表哪种类型的数据，它是如何分布的？
- 数据是多层(multilevel)的还是单层(single level)的？
这些函数的参数在所有四个版本中都非常相似。我们将在本课程后面部分学习如何分析计数和分类数据。现在，我们将重点放在连续数据上，但原理可以推广到其他类型的数据上。
以下是单层数据(没有重复测量的数据)的对比表：
| 检验类型 | 传统方法 | 回归方法 |
|---|---|---|
| 单样本t检验 | t.test(y, mu = c) |
lm(y ~ 1, offset = c) |
| 独立样本t检验 | t.test(x, y) |
lm(y ~ x) |
| 单因素方差分析 | aov(y ~ x) |
lm(y ~ x) |
| 多因素方差分析 | aov(y ~ a * b) |
lm(y ~ a * b) |
以上所有设计都是被试间设计，没有重复测量(请注意，在多因子情况下，我们可能会将a和b替换为偏差编码的数值型预测变量，原因已在交互章节中讨论过)。
混合效应模型在处理多层数据时可以发挥作用。数据通常是多层的，原因有如下三种(多个原因可能同时适用)：
- 你有一个被试内因子，和/或
- 你有假重复(pseudoreplications)，和/或
- 你有多个刺激项目(我们将在下一章中讨论)。
此时，回顾一下被试间因子与被试内因子的意义是个好主意。在sleepstudy数据中，你有被试内因子Day(实际上它更像是一个数值型变量，而不是一个因子；但它在每个被试内有多个变动值)。
假重复发生在同一条件下进行多次测量的情况下。例如，设想一个研究，你随机分配被试饮用两种饮料中的一种——酒精或水——然后进行一个简单的反应时任务，当灯光闪烁时尽快按下按钮。你可能会对每个被试进行多次反应时测量；假设你在100次试验中测量它。你有一个被试间因子(饮料类型)和每个被试的100次观测，对于每组的20个被试来说。一些新手在分析这些数据时常犯的错误是尝试进行t检验。当你有假重复(或多个刺激)时，你不能直接使用传统的t检验。你必须首先为每个被试计算均值，然后对均值而不是原始数据进行分析。虽然ANOVA的一些版本可以处理假重复，但使用线性混合效应模型可能更好，它可以更好地处理复杂的依赖结构(dependency structure)。
以下是多层数据的对比表：
| 检验类型 | 传统方法 | 回归方法 |
|---|---|---|
| 有假重复的单样本t检验 | 计算平均值并使用t.test(x_mean) |
lmer(y ~ (1 | subject), offset = c) |
| 无假重复的配对样本t检验 | t.test(x, y, paired = TRUE) |
lmer(y ~ x + (1 | subject)) |
| 有假重复的配对样本t检验 | 计算平均值并使用t.test(x_mean, y_mean) |
lmer(y ~ x + (1 + x | subject)) |
| 无假重复的重复测量方差分析 | aov(y ~ x + Error(subject)) |
lmer(y ~ x + (1 | subject)) |
| 有假重复的重复测量方差分析 | aov(y ~ x + Error(subject/x)) |
lmer(y ~ x + (1 + x | subject)) |
| 无假重复的多因子方差分析，包含a和b | aov(y ~ a * b + Error(subject)) |
lmer(y ~ a * b + (1 | subject)) |
| 有假重复的多因子方差分析 | aov(y ~ a * b + Error(subject/(a * b)) |
lmer(y ~ a * b + (1 + a * b | subject) |
广义线性模型/回归框架相比t检验和ANOVA的主要卖点之一就是灵活性。在上一章中，我们使用了sleepstudy数据，发现只有在线性混合效应模型框架才能妥善处理这些数据。尽管回归有许多优点，但如果你的数据是平衡的，且在不违反任何检验假设的情况下可以合理使用t检验或ANOVA，那么这样做是可取的；这些方法在心理学中有悠久的历史，并且更为广泛地被理解。
6.2 示例：多层数据上的独立样本\(t\)检验
让我们考虑一种情况，你在测试酒精摄入对简单反应时的影响(例如，灯亮后尽快按下按钮)。为了简化，假设你收集了14个被试的数据，这些被试被随机分配在两种干预后进行10次简单反应时试验：饮用一品脱酒精(处理条件)或饮用安慰剂饮料(安慰剂条件)。你在每组中有7名被试。请注意，实际研究中需要更多的被试。
这个网页app展示了这种研究的模拟数据。被试P01到P07属于安慰剂条件，而被试T01到T07属于处理条件。请停下并查看！
如果我们要对这些数据进行t检验，首先需要计算被试的均值，否则观测值不是独立的。你可以按如下方法进行操作。(如果你想运行下面的代码，可以从上面的网页app下载示例数据，并将其保存为independent_samples.csv)。
library(tidyverse)
dat <- read_csv("data/independent_samples.csv", col_types = "cci")
subj_means <- dat %>%
group_by(subject, cond) %>%
summarise(mean_rt = mean(RT)) %>%
ungroup()
subj_means## # A tibble: 14 × 3
## subject cond mean_rt
## <chr> <chr> <dbl>
## 1 P01 P 354
## 2 P02 P 384.
## 3 P03 P 391.
## 4 P04 P 404.
## 5 P05 P 421.
## 6 P06 P 392
## 7 P07 P 400.
## 8 T08 T 430.
## 9 T09 T 432.
## 10 T10 T 410.
## 11 T11 T 455.
## 12 T12 T 450.
## 13 T13 T 418.
## 14 T14 T 489.然后可以使用t.test()进行\(t\)检验。
##
## Welch Two Sample t-test
##
## data: mean_rt by cond
## t = -3.7985, df = 11.32, p-value = 0.002807
## alternative hypothesis: true difference in means between group P and group T is not equal to 0
## 95 percent confidence interval:
## -76.32580 -20.44563
## sample estimates:
## mean in group P mean in group T
## 392.3143 440.7000虽然这样分析没有问题，但对数据进行聚合会丢失信息。我们可以在上面的网页app中看到，实际上存在两种不同的变异来源：简单反应时的逐试次变异(trial-by-trial variability，用\(\sigma\)表示)和被试相对于总体均值(\(\gamma_{00}\))的快慢的变异。被试\(s\)在试次\(t\)上的反应时(\(Y_{st}\))的数据生成过程如下所示。
第1层:
\[\begin{equation} Y_{st} = \beta_{0s} + \beta_{1} X_{s} + e_{st} \end{equation}\]
第2层:
\[\begin{equation} \beta_{0s} = \gamma_{00} + S_{0s} \end{equation}\]
\[\begin{equation} \beta_{1} = \gamma_{10} \end{equation}\]
方差成分:
\[\begin{equation} S_{0s} \sim N\left(0, {\tau_{00}}^2\right) \end{equation}\]
\[\begin{equation} e_{st} \sim N\left(0, \sigma^2\right) \end{equation}\]
在上面的方程中，\(X_s\)是一个数值型预测变量，用于编码被试\(s\)的条件；例如：0表示安慰剂，1表示处理。
相比于一个简单模型，多层方程显得有些繁琐，但是当我们遇到更复杂的设计时，熟悉一下多层格式是值得的；我们可以将第1层和第2层简化为
\[\begin{equation} Y_{st} = \gamma_{00} + S_{0s} + \gamma_{10} X_s + e_{st}, \end{equation}\]
与上章节中的sleepstudy数据不同，现在的数据每个被试只有一个随机效应\(S_{0s}\)，没有随机斜率。每个被试只出现在两个处理条件中的一个中，因此不可能估计安慰剂与酒精的效应在被试之间的变化情况。对于这些数据，我们拟合的混合效应模型，具有随机截距但没有随机斜率，被称为随机截距模型(random intercepts model)。
随机截距模型可以充分捕捉上述两种变异来源：参数\({\tau_{00}}^2\)中的总体平均反应时的被试间变异，以及参数\(\sigma^2\)中的逐试次变异。我们可以使用下面的公式计算个体差异在总变异中所占的比例。
\[ICC = \frac{{\tau_{00}}^2}{{\tau_{00}}^2 + \sigma^2}\]
这一数值被称为组内相关系数(intra-class correlation coefficient)，告诉你数据中的聚类程度。它的范围从0到1，0表示所有变异都归因于残差变异，1表示所有变异都归因于被试间差异。
拟合随机截距模型的lmer语法是lmer(RT ~ cond + (1 | subject), dat, REML=FALSE)。首先让我们创建自己的数值型预测变量，以明确我们在使用虚拟编码。
## # A tibble: 2 × 2
## cond cond_d
## <chr> <int>
## 1 P 0
## 2 T 1现在估计模型。
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: RT ~ cond_d + (1 | subject)
## Data: dat2
##
## AIC BIC logLik deviance df.resid
## 1451.8 1463.5 -721.9 1443.8 136
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.67117 -0.66677 0.01656 0.75361 2.58447
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 329.3 18.15
## Residual 1574.7 39.68
## Number of obs: 140, groups: subject, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 392.314 8.339 47.045
## cond_d 48.386 11.793 4.103
##
## Correlation of Fixed Effects:
## (Intr)
## cond_d -0.707使用上面app中的滑块并检查lmer输出面板，直到你了解输出结果如何映射到模型参数中。
6.2.1 什么时候使用随机截距模型是合适的?
当然，混合效应模型仅在你有多层数据时才合适。随机截距模型适用于任何单样本或具有假重复的组间数据(如果在这种情况下你没有假重复，那么你就没有多层数据，可以直接使用普通回归，例如lm())。
“仅随机截距”模型也适用于设计中有被试内因素的情况，但只有当你没有假重复数据时才合适；也就是说，它只适用于每个被试在被试内因素的每个水平上只有一个观测值的情况。如果每个被试在每个水平/单元上有多个观测值，那么你需要用随机斜率来丰富你的随机效应结构，在下一节会进行讲述。如果你有多个观测值的原因是每个被试对相同的刺激集合做出反应，那么你可能需要考虑一个对被试和刺激进行交叉随机效应的混合效应模型，在下一章会进行讲述。
同样的逻辑也适用于有多个被试内因素的因子设计。在因子设计中，如果每个被试在由被试内因素相互组合形成的每个单元中只有一个观测值，那么随机截距模型是合适的。例如，如果\(A\)和\(B\)是2个两水平的被试内因素，你需要检查每个被试在\(A_1B_1\)、\(A_1B_2\)、\(A_2B_1\)和\(A_2B_2\)中是否只有1个观测值。如果有多个观测值，你将需要考虑在模型中加入随机斜率。
6.3 用回归表达研究设计和进行检验
为了在线性混合效应模型中重现所有的t检验/ANOVA风格的分析，你需要更好地理解两件事：(1)如何在回归公式中表达你的研究设计，以及(2)如何为你进行的任何检验获得p值。在线性混合效应方法中，后者并不是理所应当的，因为在许多情况下，lme4::lmer()的输出默认不会给出p值，这反映了存在多种获取p值的方法(Luke, 2017)。或者，你可能会为单个回归系数获得p值，但你想进行的检验是一个组合测试，需要同时测试多个参数。
首先，让我们仔细看看R和lme4中的回归公式语法。对于基于被试的随机效应的回归模型\(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_m x_m + e\)，公式如下：
y ~ 1 + x1 + x2 + ... + xm + (1 + x1 + ... | subject)
残差项是隐含的，所以没有提及；仅列出预测变量。左边部分(波浪号~前面)指定响应变量，右侧是预测变量。括号中的项(... | ...)特定于lme4。竖线|的右侧指定一个变量名称(在本例中为subject)，用于编码随机因子的水平。左侧指定你希望在这些水平上变化的回归系数。公式开头的1指定你需要一个截距，默认情况下已包括在内，因此可以省略。括号内的1指定一个随机截距，这也是默认包含的；括号内提到的预测变量指定随机斜率。因此，你可以将上述公式等效写成：
y ~ x1 + x2 + ... + xm + (x1 + ... | subject)。
R公式中还有另一种重要的简写方式，我们已经在交互效应章节中看到过，即被称为“星号语法”用于指定交互的a * b。如果你的设计中有两个因子a和b，并且希望模型中包含所有主效应和交互效应，可以使用：
y ~ a * b
这等同于y ~ a + b + a:b，其中a和b是编码A和B主效应的预测变量，a:b编码AB交互效应。使用星号语法指定交互比手动拼写它们要节省很多拼写(和错误)。这在下面的例子中更容易看出，该例子编码了一个包括因子A、B和C的2x2x2因子设计：
y ~ a * b * c
这等同于
y ~ a + b + c + a:b + a:c + b:c + a:b:c
其中a:b、a:c和b:c是双因子交互，a:b:c是三因子交互。你也可以在括号内的随机效应项中使用星号语法，例如：
y ~ a * b * c + (a * b * c | subject)
等同于
y ~ a + b + c + a:b + a:c + b:c + a:b:c + (a + b + c + a:b + a:c + b:c + a:b:c | subject)
尽管很复杂，但这种设计在心理学中并不罕见(例如，所有因子都在被试内，且每个条件有多个刺激)！
6.3.1 2水平以上的因子
如上所述，您可以使用公式lm(y ~ x)在回归中进行单因素ANOVA，对于单层数据，或者对于没有假重复的多层数据，使用公式lmer(y ~ x + (1 | subject))。在公式中，预测变量x应为factor()类型，R(默认情况下)会将其转换为\(k-1\)个虚拟编码的数值型预测变量(每个水平1个；有关信息，参见交互效应)。
自己编码数值预测变量通常是明智的，特别是在你的目标是进行ANOVA风格的主效应和交互效应检验的时候。这时如果你有变量是factor类型，可能会很困难。因此，对于一个名为meal(breakfast，lunch，dinner)的3水平因子，你可以创建两个变量，lunch_v_breakfast和dinner_v_breakfast，如下所示。
| factor level | lunch_v_breakfast | dinner_v_breakfast |
|---|---|---|
| breakfast | -1/3 | -1/3 |
| lunch | +2/3 | -1/3 |
| dinner | -1/3 | +2/3 |
如果你的因变量是calories，模型将会是：
calories ~ lunch_v_breakfast + dinner_v_breakfast。
但是如果你想要让meal与另一个两水平因素——time_of_week(分别编码为 -0.5 和 +0.5 的weekday和weekend)进行交互，因为你认为每顿饭消耗的卡路里在这个变量的不同水平上会有所不同。那么你的模型将是：
calories ~ (lunch_v_breakfast + dinner_v_breakfast) * time_of_week。
包含交互效应是我们选择“偏差”编码方案的原因。我们在与2水平变量相关的预测变量周围加上括号，以便每个预测变量都与time_of_week交互。上面的星形语法是以下形式的简写:
calories ~ lunch_v_breakfast + dinner_v_breakfast + time_of_week + lunch_v_breakfast:time_of_week + dinner_v_breakfast:time_of_week.
这是估计3x2因子设计参数的“回归方法”。
6.3.2 多参数检验
当你处理的设计中所有分类因子的水平都不超过两个时，与给定因子相关的回归系数检验将等同于ANOVA框架中的效应检验，前提是你使用求和或偏差编码。但是在上述示例中，我们有一个3x2设计，其中两个预测变量分别编码meal的主效应(lunch_v_breakfast和dinner_v_breakfast)。让我们模拟一些数据并使用aov()运行一个单因素ANOVA，然后我们将使用回归函数lm()复制该分析(请注意，对于多层数据的混合效应模型，使用lme4::lmer()代替base::lm()也能得到相同的结果)。
## 模拟一些数据
set.seed(62)
meals <- tibble(meal = factor(rep(c("breakfast", "lunch", "dinner"),
each = 6)),
time_of_week = factor(rep(rep(c("weekday", "weekend"),
each = 3), 3)),
calories = rnorm(18, 450, 50))
## 使用求和编码代替默认的虚拟(处理)编码
options(contrasts = c(unordered = "contr.sum", ordered = "contr.poly"))
aov(calories ~ meal * time_of_week, data = meals) %>%
summary()## Df Sum Sq Mean Sq F value Pr(>F)
## meal 2 2164 1082 0.380 0.692
## time_of_week 1 5084 5084 1.783 0.207
## meal:time_of_week 2 4767 2384 0.836 0.457
## Residuals 12 34209 2851我们得到3个\(F\)检验，每个主效应各1个(meal和time_of_week)，交互效应1个。如果我们用lm()来拟合模型呢？
## 添加我们自己的数值型变量
meals2 <- meals %>%
mutate(lunch_v_breakfast = if_else(meal == "lunch", 2/3, -1/3),
dinner_v_breakfast = if_else(meal == "dinner", 2/3, -1/3),
time_week = if_else(time_of_week == "weekend", 1/2, -1/2))
## 再此确认我们的编码
distinct(meals2, meal, time_of_week,
lunch_v_breakfast, dinner_v_breakfast, time_week)## # A tibble: 6 × 5
## meal time_of_week lunch_v_breakfast dinner_v_breakfast time_week
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 breakfast weekday -0.333 -0.333 -0.5
## 2 breakfast weekend -0.333 -0.333 0.5
## 3 lunch weekday 0.667 -0.333 -0.5
## 4 lunch weekend 0.667 -0.333 0.5
## 5 dinner weekday -0.333 0.667 -0.5
## 6 dinner weekend -0.333 0.667 0.5## 拟合回归模型
mod <- lm(calories ~ (lunch_v_breakfast + dinner_v_breakfast) *
time_week, data = meals2)
summary(mod)##
## Call:
## lm(formula = calories ~ (lunch_v_breakfast + dinner_v_breakfast) *
## time_week, data = meals2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.522 -35.895 -4.063 42.061 73.081
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 451.15 12.58 35.848 1.42e-13 ***
## lunch_v_breakfast -25.23 30.83 -0.818 0.429
## dinner_v_breakfast -20.59 30.83 -0.668 0.517
## time_week 33.61 25.17 1.335 0.207
## lunch_v_breakfast:time_week -58.31 61.65 -0.946 0.363
## dinner_v_breakfast:time_week 17.93 61.65 0.291 0.776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53.39 on 12 degrees of freedom
## Multiple R-squared: 0.2599, Adjusted R-squared: -0.04843
## F-statistic: 0.843 on 5 and 12 DF, p-value: 0.5447好的，这个输出结果看起来非常不同！在这种情况下，如何执行类似ANOVA的检验呢？你已经得到了lunch_v_breakfast和dinner_v_breakfast的估计值，但如何将其转换为对meal主效应的单一检验呢？同样，你有两个交互项lunch_v_breakfast:tow和dinner_v_breakfast:tow，如何将其转换为对交互效应的单一检验呢？
解决方案是使用模型比较进行多参数检验，在R中由anova()函数实现。为了检验meal的主效应，你需要比较包含编码该因素的两个预测变量(lunch_v_breakfast和dinner_v_breakfast)的模型与一个不包含这两个预测变量但其他部分相同的模型。你可以重新写个模型并排除这些项来重新拟合模型，或者使用update()函数并移除这些项(这是一种简写方法)。我们先写出来。
## 拟合模型
mod_main_eff <- lm(calories ~ time_week +
lunch_v_breakfast:time_week + dinner_v_breakfast:time_week,
meals2)
## 比较模型
anova(mod, mod_main_eff)## Analysis of Variance Table
##
## Model 1: calories ~ (lunch_v_breakfast + dinner_v_breakfast) * time_week
## Model 2: calories ~ time_week + lunch_v_breakfast:time_week + dinner_v_breakfast:time_week
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 12 34209
## 2 14 36373 -2 -2163.9 0.3795 0.6921好的，现在这里是使用简写方法update()函数的等效版本，它将你想要更新的模型作为第一个参数，然后使用特殊语法来进行你的更改。对于公式，我们使用. ~ . -lunch_v_breakfast -dinner_v_breakfast~。虽然这个公式看起来有些奇怪，但点.表示“在模型公式的这一侧(~左侧)保持原模型中的所有内容不变”。因此，公式. ~ .表示使用与原模型相同的公式；也就是说update(mod, . ~ .)将拟合与上述完全相同的模型。相比之下，. ~ . -x -y意味着“左侧保持不变(相同的因变量)，但从右侧移除变量x和y”。
mod_main_eff2 <- update(mod, . ~ . -lunch_v_breakfast -dinner_v_breakfast)
anova(mod, mod_main_eff2)## Analysis of Variance Table
##
## Model 1: calories ~ (lunch_v_breakfast + dinner_v_breakfast) * time_week
## Model 2: calories ~ time_week + lunch_v_breakfast:time_week + dinner_v_breakfast:time_week
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 12 34209
## 2 14 36373 -2 -2163.9 0.3795 0.6921如你所见，这给出了和之前相同的结果。
如果你想检验time_of_week的主效应，去掉这个预测因子。
## Analysis of Variance Table
##
## Model 1: calories ~ (lunch_v_breakfast + dinner_v_breakfast) * time_week
## Model 2: calories ~ lunch_v_breakfast + dinner_v_breakfast + lunch_v_breakfast:time_week +
## dinner_v_breakfast:time_week
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 12 34209
## 2 13 39294 -1 -5084.1 1.7834 0.2065试着弄清楚如何检验交互效应。
mod_interact <- update(mod, . ~ . -lunch_v_breakfast:time_week
-dinner_v_breakfast:time_week)
anova(mod, mod_interact)## Analysis of Variance Table
##
## Model 1: calories ~ (lunch_v_breakfast + dinner_v_breakfast) * time_week
## Model 2: calories ~ lunch_v_breakfast + dinner_v_breakfast + time_week
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 12 34209
## 2 14 38977 -2 -4767.4 0.8362 0.4571我们使用lm()进行模型比较得到了与使用aov()完全相同的结果。尽管这涉及更多的步骤，但学习这种方法是值得的，因为它最终会让你有更多的灵活性。