7 交叉随机因子线性混合效应模型

7.1 概括被试和刺激的组合

心理学实验的一个常见目标是检验关于在特定类型的刺激下行为发生的主张(有时是关于该行为神经基础的结论)。这些刺激可能是单词、图像、声音、视频或故事。你可能想检验的主张的一些例子包括：

  • 在听第二语言的单词时，双语者是否会受到母语单词的干扰？
  • 人们在心情好的时候和心情不好的时候对面孔的吸引力评分是否不同？
  • 相对于更中性的图像，观看舒适的图像是否有助于减轻压力?
  • 在阅读模糊描述目标个体的场景时，被潜意识启动后，人们是否更有可能对目标所属的社会群体做出假设?

关于所有这些主张需要注意的一点是，它们的形式都是：“当类型为X的个体遇到类型为Y的刺激时，我们的测量结果会发生什么变化”，其中X是从目标被试的总体中抽取的，Y是从目标刺激的总体中抽取的。换句话说，我们试图对涉及被试和刺激的采样单元相遇产生的特定类别的事件做出可推广的主张(Barr, 2017)。但就像我们不能从目标被试的总体中采样所有可能的被试一样，我们也不能从目标刺激的总体中采样所有可能的刺激。因此，在做推断时，我们需要考虑由两种采样过程引入的估计不确定性(Clark, 1973; Coleman, 1964; Judd et al., 2012; Yarkoni, 2019)。线性混合效应模型通过允许我们的模型公式中有多个随机因子而特别容易做到这一点(Baayen et al., 2008)。

这里有一个简单的例子，你感兴趣的是检验人们是否认为猫、狗或日落的图片看起来更舒适。你想对猫、狗和日落这类图片做出一般性陈述，而不是关于你碰巧采样的具体图片。假设你从Google图片中随机为每类选取4张图片(你绝对需要更多的图片才能做出一般性陈述，但我们选择少量图片以使示例简单)。因此，你的刺激表可能如下所示：

stimulus_id category file
1 cat cat1.jpg
2 cat cat2.jpg
3 cat cat3.jpg
4 cat cat4.jpg
5 dog dog1.jpg
6 dog dog2.jpg
7 dog dog3.jpg
8 dog dog4.jpg
9 sunset sunset1.jpg
10 sunset sunset2.jpg
11 sunset sunset3.jpg
12 sunset sunset4.jpg

然后你抽取了4个被试来进行舒适度评分。再次强调，4个被试对于真实研究来说太少了，但为了便于解释我们保持数据量很小。

subject_id age date
1 39 2020-05-08
2 41 2020-05-08
3 38 2020-05-16
4 38 2020-05-19

现在，因为每个受试者对每张图片都进行了“舒适度”评分，所以你会有一个完整的数据集，包含所有的subject_id和stimulus_id组合。这就是我们所说的“交叉随机因子”。你可以使用tidyr(加载tidyverse时会一起加载)中的crossing()函数来创建包含所有这些组合的表。

crossing(subjects %>% select(subject_id),
         stimuli %>% select(-category))
subject_id stimulus_id file
1 1 cat1.jpg
1 2 cat2.jpg
1 3 cat3.jpg
1 4 cat4.jpg
1 5 dog1.jpg
1 6 dog2.jpg
1 7 dog3.jpg
1 8 dog4.jpg
1 9 sunset1.jpg
1 10 sunset2.jpg
1 11 sunset3.jpg
1 12 sunset4.jpg
2 1 cat1.jpg
2 2 cat2.jpg
2 3 cat3.jpg
2 4 cat4.jpg
2 5 dog1.jpg
2 6 dog2.jpg
2 7 dog3.jpg
2 8 dog4.jpg
2 9 sunset1.jpg
2 10 sunset2.jpg
2 11 sunset3.jpg
2 12 sunset4.jpg
3 1 cat1.jpg
3 2 cat2.jpg
3 3 cat3.jpg
3 4 cat4.jpg
3 5 dog1.jpg
3 6 dog2.jpg
3 7 dog3.jpg
3 8 dog4.jpg
3 9 sunset1.jpg
3 10 sunset2.jpg
3 11 sunset3.jpg
3 12 sunset4.jpg
4 1 cat1.jpg
4 2 cat2.jpg
4 3 cat3.jpg
4 4 cat4.jpg
4 5 dog1.jpg
4 6 dog2.jpg
4 7 dog3.jpg
4 8 dog4.jpg
4 9 sunset1.jpg
4 10 sunset2.jpg
4 11 sunset3.jpg
4 12 sunset4.jpg

因为4个被试对12个刺激做出反应，所以结果表将有48行。

7.2 lme4语法用于交叉随机因子

如何分析这样的数据呢？回忆上一个章节中提到的，在lme4中，表示被试水平上的随机截距和预测变量x的随机斜率的模型的公式语法为y ~ x + (1 + x | subject_id)，其中用包含在括号中的竖线|部分进行了随机效应的指定。竖线右侧的变量subject_id指定了随机因子的水平。括号中竖线左侧的公式1 + x指定了与该因子相关的随机效应，在这种情况下是随机截距和x的随机斜率。最好将这个括号中的部分(1 + x | subject_id)理解为对lme4::lmer()关于如何构建协方差矩阵以捕捉由随机因子被试引入的方差的指示。到目前为止，你应该意识到，这个指示将指示一个二维协方差矩阵的估计，一个维度是截距方差，一个是斜率方差。

但我们并不仅局限于估计被试的随机效应；我们还可以通过简单地在公式中添加另一个项来指定刺激的随机效应。例如：

y ~ x + (1 + x | subject_id) + (1 + x | stimulus_id)。

这个公式通过x对y进行了回归，包含了被试水平上的随机截距和斜率的随机效应和刺激水平上的随机截距和斜率的随机效应。这样，拟合的模型将捕捉到我们估计中的两种不确定性来源——由被试抽样引入的不确定性以及由(刺激)项目抽样引入的不确定性。我们现在估计的是两个独立的协方差矩阵，一个用于被试，一个用于(刺激)项目。在上述例子中，这两个矩阵都将具有相同的2x2结构，但情况未必总是如此。我们可以通过改变每个竖线|符号左侧的公式灵活地更改随机效应结构。例如，如果我们有另一个预测变量w，我们可以有：

y ~ x + (x | subject_id) + (x + w | stimulus_id)

这将为被试估计相同的2x2矩阵，但(刺激)项目的协方差矩阵现在是一个3x3矩阵(截距、x的斜率和 w的斜率)。虽然这提供了很大的灵活性，但随着随机效应结构变得更加复杂，估计过程变得更加困难，并且不太可能收敛得到结果。

7.3 指定随机效应

随机效应结构的选择并不是随着线性混合效应模型而出现的新问题。在使用t检验和ANOVA的传统方法中，当你选择使用何种流程时，你隐含地选择了随机效应结构。如上一章所讨论的，如果你选择配对样本t检验而不是独立样本t检验，这相当于选择拟合lme4::lmer(y ~ x + (1 | subject))而不是lm(y ~ x)。同样地，你可以将混合模型ANOVA运行为aov(y ~ x + Error(subject_id))，这相当于随机截距模型，或者aov(y ~ x + Error(x / subject_id))，这相当于随机斜率模型。心理学中进行验证性分析(confirmatory analyses)的传统是使用研究设计所证明的最大随机效应结构。所以，如果你有一个带有假重复的单因素数据，你会使用aov(y ~ x + Error(x / subject_id))而不是aov(y ~ x + Error(subject_id))。类似地，如果你使用线性混合效应模型来分析具有假重复的数据，你应该使用lme4::lmer(y ~ x + (1 + x | subject_id))而不是lme4::lmer(y ~ x + (1 | subject_id))。换句话说，你应该考虑由重复抽样相同被试或刺激引入的所有非独立性来源。这种方法被称为最大随机效应法(maximal random effects approach)或设计驱动法(design-driven approach)来指定随机效应结构(Barr et al., 2013)。未能考虑到设计引入的依赖性可能会导致标准差过小，进而导致p值小于实际值，从而提高假阳性概率(I类错误)。在某些情况下，这可能导致效力降低，从而提高假阴性概率(II类错误)。因此，进行分析时必须密切注意随机效应结构。

线性混合效应模型几乎不可避免地包含所有设计中随机因子的随机截距。因此，如果你的随机因子是由subject_id和stimulus_id标识的被试和刺激，那么至少你的模型语法将包含(1 | subject_id) + (1 | stimulus_id)。但你会在模型中有各种预测变量，因此关键问题是：我应该允许哪些预测变量在什么抽样单位上变化？例如，如果你的模型的固定效应部分是一个2x2因子设计，因子A和B，y ~ a * b + ...，你可能会有各种随机效应结构，包括(但不限于)：

  1. 仅随机截距：y ~ a * b + (1 | subject_id) + (1 | stimulus_id)
  2. 被试水平上的a随机截距和(刺激)项目水平上的随机截距：y ~ a * b + (a | subject_id) + (1 | stimulus_id)
  3. 被试水平上的a和b的随机截距和斜率及其交互效应，(刺激)项目水平上的随机截距：y ~ a * b + (a * b | subject_id) + (1 | stimulus_id)
  4. 被试水平上的a和b的随机截距和斜率及其交互效应，(刺激)项目水平上的a和b的随机截距和斜率及其交互效应：y ~ a * b + (a * b | subject_id) + (a * b | stimulus_id)。

有一点需要明确：

“由设计证明的最大随机效应结构”与“最大可能的随机效应结构”不同；也就是说，它并不意味着自动将所有预测变量的所有随机斜率放入所有随机因子中。你必须遵循下一部分中的随机效应指南，以决定是否应包含某个随机斜率。

一些作者建议采用“数据驱动”的随机效应结构替代设计驱动的随机效应结构，建议研究者仅在数据进一步证明的情况下才包含设计所证明的随机斜率(Matuschek et al., 2017)。例如，你可以使用虚无假设检验来确定是否包括一个被试水平上的x随机斜率显著提高了模型拟合，只有在显著提高时才包含该效应。虽然当随机斜率非常小时，这可能会在理论意义上提高检验的效力，但也会增加未知的假阳性风险，因此在检验的情景中是否应该采用这种方法是值得怀疑的。所以，我们不推荐基于数据的方法。

7.3.1 为分类因子选择随机效应的规则

线性混合效应模型的随机效应结构——即关于哪些效应在什么抽样单位上变化的假设——对于确保你的参数反映抽样引入的不确定性至关重要(Barr et al., 2013)。首先，请注意我们专注于代表设计变量(design variable)的预测变量，这些变量具有理论意义，并且你将对其进行推断检验。如果你有代表控制变量(control variable)的预测变量，而你不打算对其进行统计检验，那么通常不需要随机斜率。

以下规则源自 Barr et al. (2013) 和 Barr (2013) 。如果你想了解这些指南的更多信息，请参考这些论文。请记住，只有当你有重复测量数据时，你才能使用混合效应模型，这可能是由于假重复和/或存在被试内(或刺激内)因子。在交叉随机因子的情况下，你不可避免地会有假重复——由于多个被试的多个刺激，每个被试的多个观测值，以及多个被试的每个刺激的多个观测值。确定随机效应结构的关键是弄清楚哪些因子是被试内或刺激内，以及设计中的假重复位于何处。你对被试使用一次规则以确定公式(1 + ... | subject_id)部分的形式，并对刺激使用一次规则以确定公式(1 + ... | stimulus_id)部分的形式。以下规则中的“单位”或“抽样单位”一词，根据需要替换为“被试”或“刺激”。

以下是规则：

  1. 如果对抽样单位进行了重复测量，则需要该随机因子的随机截距：(1 | unit_id)；
  2. 如果一个因子x是单位间的，则不需要该因子的随机斜率；
  3. 确定所考虑的单位的被试内因子的最高阶交互效应。如果在这些组合定义的每个单元内有假重复(即每个单元的多个观测值)，那么对于该单位，你将需要该交互效应的斜率以及所有低阶效应的斜率。如果没有假重复，则不需要任何随机斜率。

前两条规则很简单，但第三条需要一些解释。首先问自己：我们怎么知道某个因子是单位间还是单位内的？

一个简单的方法是使用dplyr::count()函数(加载tidyverse时同时加载)，它会提供频率计数。假设你对因子A是否在被试内还是被试者间感兴趣，对下面虚拟的2x2x2因子数据abc_data，其中 \(A\)、\(B\)和\(C\)是你的设计因子。

library(tidyverse)

## 运行代码创建"abc_data"表
abc_subj <- tibble(subject_id = 1:4,
                   B = rep(c("B1", "B2"), times = 2))

abc_item  <- tibble(stimulus_id = 1:4,
                    A = rep(c("A1", "A2"), each = 2),
                    C = rep(c("C1", "C2"), times = 2))

abc_data <- crossing(abc_subj, abc_item) %>%
  select(subject_id, stimulus_id, everything())

查看\(A\)是被试内还是被试间，使用：

abc_data %>%
  count(subject_id, A)
## # A tibble: 8 × 3
##   subject_id A         n
##        <int> <chr> <int>
## 1          1 A1        2
## 2          1 A2        2
## 3          2 A1        2
## 4          2 A2        2
## 5          3 A1        2
## 6          3 A2        2
## 7          4 A1        2
## 8          4 A2        2

在结果表中，你可以看到每个被试都接受了\(A\)的所有水平，说明这是个被试内因子。那么对于\(B\)和\(C\)呢？

abc_data %>%
  count(subject_id, B)
## # A tibble: 4 × 3
##   subject_id B         n
##        <int> <chr> <int>
## 1          1 B1        4
## 2          2 B2        4
## 3          3 B1        4
## 4          4 B2        4
abc_data %>%
  count(subject_id, C)
## # A tibble: 8 × 3
##   subject_id C         n
##        <int> <chr> <int>
## 1          1 C1        2
## 2          1 C2        2
## 3          2 C1        2
## 4          2 C2        2
## 5          3 C1        2
## 6          3 C2        2
## 7          4 C1        2
## 8          4 C2        2

\(B\)是被试间(每个被试只接受1个水平)，\(C\)是被试内(每个被试接受所有水平)。

练习

回答关于abc_data的问题。

  • 因子\(A\)的水平是在刺激间还是刺激内? 
abc_data %>%
  count(stimulus_id, A)
## # A tibble: 4 × 3
##   stimulus_id A         n
##         <int> <chr> <int>
## 1           1 A1        4
## 2           2 A1        4
## 3           3 A2        4
## 4           4 A2        4
# between
  • 因子\(B\)的水平是在刺激间还是刺激内? 
abc_data %>%
  count(stimulus_id, B)
## # A tibble: 8 × 3
##   stimulus_id B         n
##         <int> <chr> <int>
## 1           1 B1        2
## 2           1 B2        2
## 3           2 B1        2
## 4           2 B2        2
## 5           3 B1        2
## 6           3 B2        2
## 7           4 B1        2
## 8           4 B2        2
# within
  • 因子\(C\)的水平是在刺激间还是刺激内? 
abc_data %>%
  count(stimulus_id, C)
## # A tibble: 4 × 3
##   stimulus_id C         n
##         <int> <chr> <int>
## 1           1 C1        4
## 2           2 C2        4
## 3           3 C1        4
## 4           4 C2        4
# between

好的，我们已经确定了哪些因子是被试内和哪些被试间，以及哪些因子是刺激内和哪些刺激间。

第二条规则告诉我们，如果一个因子是单位间的，你不需要该因子的随机斜率。事实上，无法估计单位间因子的随机斜率。如果你考虑到随机斜率捕捉效应在单位之间的变化，那么只有当你在该因子的所有水平上测量了你的响应变量时，才能估计出这种变化。例如，如果你有一个名为处理组(实验组、对照组)的2水平因子，除非某个被试经历了该因子的两个水平(即它必须是被试内的)，否则无法估计“治疗”对该被试的效应。

我们现在如何应用第三条规则来确定我们对单位内因子需要哪些随机斜率呢？

考虑到\(A\)和\(C\)是被试内的，\(B\)是被试间的。所以被试内因子的最高阶交互作用是\(AC\)。如果我们对每个被试在每个\(AC\)组合中有假重复，那么我们将需要\(AC\)交互效应的随机斜率以及主效应\(A\)和\(C\)的随机斜率。我们如何确定这一点呢？

abc_data %>%
  count(subject_id, A, C)
## # A tibble: 16 × 4
##    subject_id A     C         n
##         <int> <chr> <chr> <int>
##  1          1 A1    C1        1
##  2          1 A1    C2        1
##  3          1 A2    C1        1
##  4          1 A2    C2        1
##  5          2 A1    C1        1
##  6          2 A1    C2        1
##  7          2 A2    C1        1
##  8          2 A2    C2        1
##  9          3 A1    C1        1
## 10          3 A1    C2        1
## 11          3 A2    C1        1
## 12          3 A2    C2        1
## 13          4 A1    C1        1
## 14          4 A1    C2        1
## 15          4 A2    C1        1
## 16          4 A2    C2        1

这向我们展示了对于每个\(AC\)组合，只有1个观测值，因此我们不需要\(AC\)的随机斜率，也不需要\(A\)或\(C\)的随机斜率。被试水平上的随机效应公式就是(1 | subject_id)。

刺激水平上的随机因子需要哪些随机斜率?

你有1个刺激内因子\(B\)，它进行了假重复。

abc_data %>%
  count(stimulus_id, B)
## # A tibble: 8 × 3
##   stimulus_id B         n
##         <int> <chr> <int>
## 1           1 B1        2
## 2           1 B2        2
## 3           2 B1        2
## 4           2 B2        2
## 5           3 B1        2
## 6           3 B2        2
## 7           4 B1        2
## 8           4 B2        2

因此对于刺激的公式是(B | stimulus_id), 完整的lme4公式是:

y ~ A * B * C + (1 | subject_id) + (B | stimulus_id).

7.3.2 排除不收敛和“畸形拟合(singular fit)”

当你试图拟合具有最大随机效应的模型时，可能会遇到几个不同的问题。回想一下，线性混合效应模型的估计算法是迭代(iterative)——即拟合算法以逐步的方式搜索使数据最合适的参数值。有时它会一直搜索但找不到，这种情况下，你会收到“收敛警告”。当这种情况发生时，信任未收敛模型的任何估计值都不是一个好主意，你需要简化模型结构然后重试。

另一种情况是你会收到关于“畸形拟合”的信息。当估计过程包含一个或多个随机因子产生的方差-协方差矩阵有(1)完全或近乎完全(1.00，-1.00)的正或负相关，(2)一个或多个方差接近于0，或(3)两者兼有时，会出现这种信息。可能可以忽略此信息，但也有理由简化模型结构直到信息消失。

如何简化模型以处理收敛问题或畸形拟合信息呢？这需要小心操作。我建议以下策略：

  1. 将所有协方差参数约束为0。这可以使用双竖线||语法来完成，例如，将(a * b | subject)更改为(a * b || subject)。如果仍然遇到估计问题，那么：
  2. 检查非收敛或畸形模型的参数估计值。是否有任何斜率变量为0或接近于0？删除这些变量并重新拟合模型，重复此步骤直到收敛警告/畸形拟合信息消失。

有关收敛问题的更多技术细节以及如何处理，参见?lme4::convergence和?lme4::isSingular。

注：推荐阅读混合线性模型的实现加深理解。同时可以考虑使用贝叶斯方法，参见 Fiechter (2024) 。

7.4 模拟交叉随机因子的数据

对于这些练习，我们将生成与实验相对应的模拟数据，该实验具有被试内和项目间的单、双水平因子(自变量)。假设实验涉及一组单词的词汇决策(例如，“PINT”是单词还是非单词？)，因变量是反应时(以毫秒为单位)，自变量是词类(名词与动词)。我们想把被试和单词都视为随机因子(以便我们能够推广到被试遇到单词这一事件的总体)。你可以使用网页app进行操作，这允许你操纵数据生成参数并观察其对数据的影响。

到目前为止，你应该已经掌握了模拟具有交叉随机效应的研究数据所需的所有拼图。DeBruine & Barr (2021) 提供了更详细的、逐步的练习演示。

这里是被试\(s\)和项目\(i\)的反应时\(Y_{si}\)的数据生成过程(data generating process，DGP)：

第1层:

\[\begin{equation} Y_{si} = \beta_{0s} + \beta_{1} X_{i} + e_{si} \end{equation}\]

第2层:

\[\begin{equation} \beta_{0s} = \gamma_{00} + S_{0s} + I_{0i} \end{equation}\]

\[\begin{equation} \beta_{1} = \gamma_{10} + S_{1s} \end{equation}\]

方差成分:

\[\begin{equation} \langle S_{0s}, S_{1s} \rangle \sim N\left(\langle 0, 0 \rangle, \mathbf{\Sigma}\right) \end{equation}\]

\[\begin{equation} \mathbf{\Sigma} = \left(\begin{array}{cc}{\tau_{00}}^2 & \rho\tau_{00}\tau_{11} \\ \rho\tau_{00}\tau_{11} & {\tau_{11}}^2 \\ \end{array}\right) \end{equation}\]

\[\begin{equation} I_{0s} \sim N\left(0, {\omega_{00}}^2\right) \end{equation}\]

\[\begin{equation} e_{si} \sim N\left(0, \sigma^2\right) \end{equation}\]

在上述等式中，\(X_i\)是是一个数值型预测变量编码，它决定项目\(i\)所处的条件，如-.5表示名词，.5表示动词。

我们可以化简第1层和第2层为：

\[Y_{si} = \beta_0 + S_{0s} + I_{0i} + (\beta_1 + S_{1s})X_{i} + e_{si}\]

其中：

参数 符号 描述
\(Y_{si}\) Y 被试\(s\)对于项目\(i\)的反应时;
\(\beta_0\) b0 总均值;
\(S_{0s}\) S_0s 被试\(s\)的随机截距;
\(I_{0i}\) I_0i 项目\(i\)的随机截距;
\(\beta_1\) b1 单词类型的固定效应(斜率);
\(S_{1s}\) S_1s 被试水平上的随机斜率;
\(X_{i}\) cond 字型的偏差编码预测变量;
\(\tau_{00}\) tau_00 被试水平上的随机截距标准差
\(\tau_{11}\) tau_11 被试水平上的随机斜率标准差
\(\rho\) rho 随机截距与斜率的相关性
\(\omega_{00}\) omega_00 项目水平上的随机截距标准差
\(e_{si}\) err 被试\(s\)项目\(i\)的误差残差
\(\sigma\) sig 误差残差标准差

7.4.1 设置环境并定义DGP的参数

如果你想在这个练习中得到和其他人一样的结果，我们都应该为随机数生成器提供相同的值。同时，让我们加载我们需要的包。

library(lme4)
library(tidyverse)

set.seed(11709)

现在让我们定义DGP参数。

nsubj <- 100 # 被试数
nitem <- 50  # 必须是偶数

b0 <- 800 # 总均值
b1 <- 80 # 80毫秒差异
effc <- c(-.5, .5) # 偏差编码

omega_00 <- 80 # 项目水平上的随机截距标准差(omega_00)

## 被试水平上的方差-协方差矩阵
tau_00 <- 100 # 被试水平上的随机截距标准差
tau_11 <- 40 # 被试水平上的随机斜率标准差
rho <- .2 # 截距和斜率的相关性

sig <- 200 # 残差(标准差)

你将创建3个表：

名称 描述
subjects 包含subj_id和被试随机效应的被试数据表
items 包含item_id和项目随机效应的刺激数据表
trials 列举被试/刺激之间组合的试验表

然后将三个表中的信息合并在一起，根据上面的模型公式计算响应变量。

7.4.2 生成刺激样本

让我们随机生成50项。像下面一样创建一个叫item的表，其中iri是项目水平上的随机截距(提取于方差为\(\omega_{00}^2\) = 6400的正态分布)。一半的词形是名词(cond = -.5)，一半的词形是动词(cond = .5)。

## # A tibble: 50 × 3
##    item_id  cond    I_0i
##      <int> <dbl>   <dbl>
##  1       1  -0.5   14.9 
##  2       2   0.5  -86.3 
##  3       3  -0.5  -12.8 
##  4       4   0.5  -13.9 
##  5       5  -0.5   55.6 
##  6       6   0.5  -45.9 
##  7       7  -0.5  -42.0 
##  8       8   0.5  -87.6 
##  9       9  -0.5  -97.4 
## 10      10   0.5  -85.2 
## 11      11  -0.5  135.  
## 12      12   0.5   83.2 
## 13      13  -0.5  -44.7 
## 14      14   0.5    8.59
## 15      15  -0.5 -156.  
## 16      16   0.5  -57.6 
## 17      17  -0.5  -38.7 
## 18      18   0.5   39.6 
## 19      19  -0.5  105.  
## 20      20   0.5   30.3 
## 21      21  -0.5 -115.  
## 22      22   0.5   -3.40
## 23      23  -0.5 -218.  
## 24      24   0.5   53.0 
## 25      25  -0.5  -86.9 
## 26      26   0.5  -65.4 
## 27      27  -0.5  172.  
## 28      28   0.5 -152.  
## 29      29  -0.5   25.1 
## 30      30   0.5 -156.  
## 31      31  -0.5   47.7 
## 32      32   0.5  -46.3 
## 33      33  -0.5   48.0 
## 34      34   0.5   62.8 
## 35      35  -0.5  -75.4 
## 36      36   0.5  -35.9 
## 37      37  -0.5  -48.5 
## 38      38   0.5   29.3 
## 39      39  -0.5  -55.5 
## 40      40   0.5   69.5 
## 41      41  -0.5  196.  
## 42      42   0.5   77.6 
## 43      43  -0.5  -45.0 
## 44      44   0.5  204.  
## 45      45  -0.5   32.1 
## 46      46   0.5  -63.9 
## 47      47  -0.5  145.  
## 48      48   0.5   66.2 
## 49      49  -0.5  -23.9 
## 50      50   0.5   97.3

rep()

rnorm()

items <- tibble(item_id = 1:nitem,
                cond = rep(c(-.5, .5), times = nitem / 2),
                I_0i = rnorm(nitem, 0, sd = omega_00))

7.4.3 生成被试样本

要生成被试水平上的随机效应，你需要从双变量正态分布中生成数据。为此，我们将使用函数MASS::mvrnorm。

记住：不要仅仅为了获得这个函数而运行library("MASS")，因为MASS包中有一个函数select()会覆盖tidyverse版本的select()。由于我们只需要MASS包中的mvrnorm()函数，我们可以直接通过pkgname::function语法来访问它，即MASS::mvrnorm()。

你的被试表格应如下所示：

## # A tibble: 100 × 3
##     subj_id      S_0s     S_1s
##       <int>     <dbl>    <dbl>
##   1       1  -80.0      -0.763
##   2       2   44.6      54.5  
##   3       3    8.74    -20.4  
##   4       4  -38.6     -23.8  
##   5       5  -83.3      29.2  
##   6       6  -70.9     -13.8  
##   7       7  -21.4      46.0  
##   8       8    2.33      8.39 
##   9       9   62.3     -58.2  
##  10      10  238.        7.72 
##  11      11  -92.5       2.14 
##  12      12   58.5     -65.8  
##  13      13 -204.      -38.8  
##  14      14  -91.6       5.46 
##  15      15   51.1     -38.8  
##  16      16  142.      -12.9  
##  17      17   46.0       6.60 
##  18      18  -56.7     -54.8  
##  19      19  -10.1      62.1  
##  20      20 -226.      -19.3  
##  21      21 -158.      -18.5  
##  22      22  102.        8.99 
##  23      23  -12.7     -70.6  
##  24      24  135.       -9.50 
##  25      25   62.0     -52.5  
##  26      26    0.0653   32.8  
##  27      27 -117.       70.8  
##  28      28 -232.        3.43 
##  29      29   70.9      50.8  
##  30      30 -123.       22.8  
##  31      31  268.       30.0  
##  32      32  -18.7     -25.0  
##  33      33   50.8     -31.0  
##  34      34  -43.1     -28.9  
##  35      35  -10.1      28.3  
##  36      36   65.6      18.2  
##  37      37 -123.       -4.63 
##  38      38  -94.8      10.3  
##  39      39   77.7     -22.5  
##  40      40  -59.1      52.4  
##  41      41  -91.2    -103.   
##  42      42  -66.6      -2.14 
##  43      43   -4.40      0.305
##  44      44   69.7      10.2  
##  45      45  -77.5     -10.4  
##  46      46  -17.8     -48.2  
##  47      47 -103.       47.0  
##  48      48   22.8     -39.3  
##  49      49  -31.1     -34.9  
##  50      50  -26.4      40.0  
##  51      51   47.8      26.0  
##  52      52  -93.2     -42.7  
##  53      53   28.9      51.4  
##  54      54  -19.3      11.5  
##  55      55   53.6      21.5  
##  56      56  -27.4     -21.4  
##  57      57  -67.7     -32.1  
##  58      58   59.2      13.4  
##  59      59  -53.1       2.44 
##  60      60  104.        7.41 
##  61      61  -20.7     -78.7  
##  62      62   55.9     -15.7  
##  63      63  114.      -29.1  
##  64      64  -57.7     -34.7  
##  65      65  -38.7      -9.14 
##  66      66 -106.      -58.0  
##  67      67   99.1     -37.6  
##  68      68  -56.9      21.0  
##  69      69  -50.4      -0.407
##  70      70   27.5      -2.69 
##  71      71  139.      -32.2  
##  72      72   44.9       8.53 
##  73      73  -14.8      71.7  
##  74      74   33.7     -52.6  
##  75      75    2.03     27.8  
##  76      76 -134.       37.0  
##  77      77   24.4      20.7  
##  78      78  -60.6     -36.7  
##  79      79   31.1      16.9  
##  80      80  -34.9       9.68 
##  81      81  206.       17.3  
##  82      82   -7.19    -25.4  
##  83      83  182.       46.0  
##  84      84   55.7      21.7  
##  85      85 -149.      -44.0  
##  86      86 -193.      -73.2  
##  87      87  167.       13.9  
##  88      88  160.        3.87 
##  89      89   84.1      82.1  
##  90      90   97.2      -6.55 
##  91      91 -205.     -125.   
##  92      92  -75.1       6.76 
##  93      93  -95.3     -46.5  
##  94      94  106.       38.6  
##  95      95  -42.4      11.3  
##  96      96   74.0     -21.1  
##  97      97 -245.      -25.3  
##  98      98 -113.       -1.88 
##  99      99   68.8      30.6  
## 100     100  136.       44.2

回忆一下：

  • tau_00: 被试水平上随机截距的标准差
  • tau_11: 被试水平上随机斜率的标准差
  • rho : 截距和斜率的相关系数

covariance = rho * tau_00 * tau_11

matrix(    tau_00^2,            rho * tau_00 * tau_11,
        rho * tau_00 * tau_11,      tau_11^2, ...)
as_tibble(mx) %>%
  mutate(subj_id = ...)
cov <- rho * tau_00 * tau_11

mx <- matrix(c(tau_00^2, cov,
               cov,      tau_11^2),
             nrow = 2)

by_subj_rfx <- MASS::mvrnorm(nsubj, mu = c(S_0s = 0, S_1s = 0), Sigma = mx)

subjects <- as_tibble(by_subj_rfx) %>%
  mutate(subj_id = row_number()) %>%
  select(subj_id, everything())

7.4.4 生产试次样本

每个试次都是独特被试和刺激的相遇，在这个实验中，每个被试都接受了每个刺激。生产一个表trials列出实验中所有的试次。注意：每个被试都会接受每种刺激1次。使用crossing()函数创建所有可能的试次。

现在应用这个示例生成下面的表，其中err是残差项，提取于\(N \sim \left(0, \sigma^2\right)\)，\(\sigma\)是err_sd。

## # A tibble: 5,000 × 3
##    subj_id item_id    err
##      <int>   <int>  <dbl>
##  1       1       1  382. 
##  2       1       2  283. 
##  3       1       3   30.4
##  4       1       4 -282. 
##  5       1       5 -239. 
##  6       1       6   73.4
##  7       1       7  -98.4
##  8       1       8 -189. 
##  9       1       9 -410. 
## 10       1      10  102. 
## # ℹ 4,990 more rows
trials <- crossing(subj_id = subjects %>% pull(subj_id),
                   item_id = items %>% pull(item_id)) %>%
  mutate(err = rnorm(nrow(.), mean = 0, sd = sig))

7.4.5 合并subjects，items和trials

合并subjects，items和trials中的信息，创建完整的数据集dat_sim，看起来像这样：

## # A tibble: 5,000 × 7
##    subj_id item_id  S_0s  I_0i   S_1s  cond    err
##      <int>   <int> <dbl> <dbl>  <dbl> <dbl>  <dbl>
##  1       1       1 -80.0  14.9 -0.763  -0.5  382. 
##  2       1       2 -80.0 -86.3 -0.763   0.5  283. 
##  3       1       3 -80.0 -12.8 -0.763  -0.5   30.4
##  4       1       4 -80.0 -13.9 -0.763   0.5 -282. 
##  5       1       5 -80.0  55.6 -0.763  -0.5 -239. 
##  6       1       6 -80.0 -45.9 -0.763   0.5   73.4
##  7       1       7 -80.0 -42.0 -0.763  -0.5  -98.4
##  8       1       8 -80.0 -87.6 -0.763   0.5 -189. 
##  9       1       9 -80.0 -97.4 -0.763  -0.5 -410. 
## 10       1      10 -80.0 -85.2 -0.763   0.5  102. 
## # ℹ 4,990 more rows

inner_join()

dat_sim <- subjects %>%
  inner_join(trials, "subj_id") %>%
  inner_join(items, "item_id") %>%
  arrange(subj_id, item_id) %>%
  select(subj_id, item_id, S_0s, I_0i, S_1s, cond, err)

7.4.6 创建响应变量

根据下面的公式模型将响应变量Y添加到数据中：

\[Y_{si} = \beta_0 + S_{0s} + I_{0i} + (\beta_1 + S_{1s})X_{i} + e_{si}\]

那么，最终的表(dat_sim2)看起来是这样的：

## # A tibble: 5,000 × 8
##    subj_id item_id     Y  S_0s  I_0i   S_1s  cond    err
##      <int>   <int> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>
##  1       1       1 1078. -80.0  14.9 -0.763  -0.5  382. 
##  2       1       2  957. -80.0 -86.3 -0.763   0.5  283. 
##  3       1       3  698. -80.0 -12.8 -0.763  -0.5   30.4
##  4       1       4  464. -80.0 -13.9 -0.763   0.5 -282. 
##  5       1       5  497. -80.0  55.6 -0.763  -0.5 -239. 
##  6       1       6  787. -80.0 -45.9 -0.763   0.5   73.4
##  7       1       7  540. -80.0 -42.0 -0.763  -0.5  -98.4
##  8       1       8  483. -80.0 -87.6 -0.763   0.5 -189. 
##  9       1       9  173. -80.0 -97.4 -0.763  -0.5 -410. 
## 10       1      10  776. -80.0 -85.2 -0.763   0.5  102. 
## # ℹ 4,990 more rows

注意：这是该模型完整的分解表

... %>% 
  mutate(Y = ...) %>%
  select(...)
dat_sim2 <- dat_sim %>%
  mutate(Y = b0 + S_0s + I_0i + (S_1s + b1) * cond + err) %>%
  select(subj_id, item_id, Y, everything())

7.4.7 拟合模型

现在你创建完了模拟数据，使用lme4::lmer()拟合模型，并运行summary()吧。

mod_sim <- lmer(Y ~ cond + (1 + cond | subj_id) + (1 | item_id),
                dat_sim2, REML = FALSE)

summary(mod_sim, corr = FALSE)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Y ~ cond + (1 + cond | subj_id) + (1 | item_id)
##    Data: dat_sim2
## 
##      AIC      BIC   logLik deviance df.resid 
##  67639.4  67685.0 -33812.7  67625.4     4993 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6357 -0.6599 -0.0251  0.6767  3.7685 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  subj_id  (Intercept)  9464.8   97.29       
##           cond          597.7   24.45   0.68
##  item_id  (Intercept)  8087.0   89.93       
##  Residual             40305.0  200.76       
## Number of obs: 5000, groups:  subj_id, 100; item_id, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   793.29      16.26  48.782
## cond           77.65      26.18   2.967

现在看看你能否在summary()的输出结果中识别出数据生成参数。

首先，尝试找到\(\beta_0\)和\(\beta_1\)。

parameter variable input estimate
\(\hat{\beta}_0\) b0 800 793.293
\(\hat{\beta}_1\) b1 80 77.652

尝试找到随机效应的估计参数\(\tau_{00}\)，\(\tau_{11}\)，\(\rho\)，\(\omega_{00}\)和\(\sigma\)。

parameter variable input estimate
\(\hat{\tau}_{00}\) tau_00 100.0 97.287
\(\hat{\tau}_{11}\) tau_11 40.0 24.448
\(\hat{\rho}\) rho 0.2 0.675
\(\hat{\omega}_{00}\) omega_00 80.0 89.928
\(\hat{\sigma}\) sig 200.0 200.761

注：推荐阅读 Brown (2021)

参考文献

Baayen, R. H., Davidson, D. J., & Bates, D. M. (2008). Mixed-effects modeling with crossed random effects for subjects and items. Journal of Memory and Language, 59(4), 390–412. https://doi.org/10.1016/j.jml.2007.12.005
Barr, D. J. (2013). Random effects structure for testing interactions in linear mixed-effects models. Frontiers in Psychology, 4, 328. https://doi.org/10.3389/fpsyg.2013.00328
Barr, D. J. (2017). Generalizing over encounters: Statistical and theoretical considerations. In Oxford handbook of psycholinguistics. Oxford University Press. https://doi.org/10.1093/oxfordhb/9780198786825.013.39
Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255–278. https://doi.org/10.1016/j.jml.2012.11.001
Brown, V. A. (2021). An introduction to linear mixed-effects modeling in r. Advances in Methods and Practices in Psychological Science, 4(1). https://doi.org/10.1177/2515245920960351
Clark, H. H. (1973). The language-as-fixed-effect fallacy: A critique of language statistics in psychological research. Journal of Verbal Learning and Verbal Behavior, 12(4), 335–359. https://doi.org/10.1016/S0022-5371(73)80014-3
Coleman, E. B. (1964). Generalizing to a language population. Psychological Reports, 14(1), 219–226. https://doi.org/10.2466/pr0.1964.14.1.219
DeBruine, L., & Barr, D. J. (2021). Understanding mixed effects models through data simulation. Advances in Methods and Practice in Psychological Science, 4(1). https://doi.org/10.1177/2515245920965119
Fiechter, J. L. (2024). Drawing generalizable conclusions from multilevel models: Commentary on van de calseyde and efendić (2022). Psychological Science, 35(6), 694–699. https://doi.org/10.1177/09567976241245411
Judd, C. M., Westfall, J., & Kenny, D. A. (2012). Treating stimuli as a random factor in social psychology: A new and comprehensive solution to a pervasive but largely ignored problem. Journal of Personality and Social Psychology, 103(1), 54–69. https://doi.org/10.1037/a0028347
Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H., & Bates, D. (2017). Balancing Type I error and power in linear mixed models. Journal of Memory and Language, 94, 305–315. https://doi.org/10.1016/j.jml.2017.01.001
Yarkoni, T. (2019). The generalizability crisis. https://doi.org/10.31234/osf.io/jqw35