使用 lmer 进行重复测量线性混合效应模型

机器算法验证 r 方差分析 混合模式 重复测量 lme4-nlme
2022-01-30 08:37:32

编辑2：我最初认为我需要对一个因素进行重复测量的双因素方差分析，但我现在认为线性混合效应模型更适合我的数据。我想我几乎知道需要发生什么，但仍然有几点感到困惑。

我需要分析的实验如下所示：

  • 受试者被分配到几个治疗组之一
  • 多天对每个受试者进行测量
  • 所以：
    • 对象嵌套在治疗中
    • 治疗与日交叉

（每个受试者只被分配一次治疗，并且每天对每个受试者进行测量）

我的数据集包含以下信息：

  • 主题 = 阻塞因子（随机因子）
  • 天 = 受试者内或重复测量因子（固定因子）
  • 治疗=主体因素之间（固定因素）
  • Obs = 测量（因）变量

更新 好的，所以我去和一位统计学家交谈，但他是 SAS 用户。他认为模型应该是：

治疗 + 天 + 受试者（治疗）+ 天*受试者（治疗）

显然他的符号与 R 语法不同，但这个模型应该解释：

  • 治疗（固定）
  • 日（固定）
  • 治疗*日互动
  • 嵌套在治疗中的对象（随机）
  • 天与“治疗中的对象”交叉（随机）

那么，这是正确的语法吗？

m4 <- lmer(Obs~Treatment*Day + (1+Treatment/Subject) + (1+Day*Treatment/Subject), mydata)

我特别关心与“治疗中的对象”部分交叉的日子是否正确。是否有人熟悉 SAS，或者有信心了解他的模型中发生了什么，能够评论我对 R 语法的悲伤尝试是否匹配？

这是我之前在构建模型和编写语法方面的尝试（在答案和评论中讨论）：

m1 <- lmer(Obs ~ Treatment * Day + (1 | Subject), mydata)

我如何处理主题嵌套在治疗中的事实？与m1以下有何不同：

m2 <- lmer(Obs ~ Treatment * Day + (Treatment|Subject), mydata)
m3 <- lmer(Obs ~ Treatment * Day + (Treatment:Subject), mydata)

并且是m2和m3等价的（如果不是，为什么）？

另外，如果我想指定相关结构（如correlation = corAR1），是否需要使用 nlme 而不是 lme4 ？根据重复测量，对于对一个因素进行重复测量的重复测量分析，协方差结构（同一受试者的测量之间的相关性的性质）很重要。

当我尝试进行重复测量方差分析时，我决定使用 II 型 SS；这仍然相关吗？如果是，我该如何指定？

以下是数据的示例：

mydata <- data.frame(
  Subject  = c(13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 
               34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 
               19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
               40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 
               29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65), 
  Day       = c(rep(c("Day1", "Day3", "Day6"), each=28)), 
  Treatment = c(rep(c("B", "A", "C", "B", "C", "A", "A", "B", "A", "C", "B", "C", 
                      "A", "A", "B", "A", "C", "B", "C", "A", "A"), each = 4)), 
  Obs       = c(6.472687, 7.017110, 6.200715, 6.613928, 6.829968, 7.387583, 7.367293, 
                8.018853, 7.527408, 6.746739, 7.296910, 6.983360, 6.816621, 6.571689, 
                5.911261, 6.954988, 7.624122, 7.669865, 7.676225, 7.263593, 7.704737, 
                7.328716, 7.295610, 5.964180, 6.880814, 6.926342, 6.926342, 7.562293, 
                6.677607, 7.023526, 6.441864, 7.020875, 7.478931, 7.495336, 7.427709, 
                7.633020, 7.382091, 7.359731, 7.285889, 7.496863, 6.632403, 6.171196, 
                6.306012, 7.253833, 7.594852, 6.915225, 7.220147, 7.298227, 7.573612, 
                7.366550, 7.560513, 7.289078, 7.287802, 7.155336, 7.394452, 7.465383, 
                6.976048, 7.222966, 6.584153, 7.013223, 7.569905, 7.459185, 7.504068, 
                7.801867, 7.598728, 7.475841, 7.511873, 7.518384, 6.618589, 5.854754, 
                6.125749, 6.962720, 7.540600, 7.379861, 7.344189, 7.362815, 7.805802, 
                7.764172, 7.789844, 7.616437, NA, NA, NA, NA))
2个回答

我认为你的方法是正确的。模型m1为每个主题指定一个单独的截距。模型m2为每个主题添加了一个单独的斜率。您的斜率跨越几天，因为受试者只参加一个治疗组。如果您m2按如下方式编写模型，则更明显的是您为每个主题建模了一个单独的截距和斜率

m2 <- lmer(Obs ~ Treatment * Day + (1+Day|Subject), mydata)

这相当于：

m2 <- lmer(Obs ~ Treatment + Day + Treatment:Day + (1+Day|Subject), mydata)

即治疗的主要影响，天和两者之间的相互作用。

我认为只要您不在治疗组中重复受试者 ID，您就不必担心嵌套问题。哪种模型是正确的，实际上取决于您的研究问题。是否有理由相信除了治疗效果之外，受试者的斜率也会有所不同？您可以运行这两种模型并将它们与它们进行比较，anova(m1,m2)以查看数据是否支持其中任何一种。

我不确定你想用 model 表达m3什么？嵌套语法使用 a /，例如(1|group/subgroup)。

我认为您不必担心与这么少的时间点的自相关。

我对您的自相关错误问题（也没有关于 lme4 与 nlme 中的不同实现）发表评论感到不舒服，但我可以与其他人交谈。

您的模型m1是一个随机截距模型，其中您包含了治疗和天之间的跨级别交互（天的影响允许在治疗组之间变化）。为了允许随时间的变化在参与者之间有所不同（即明确地模拟随时间变化的个体差异），您还需要允许 Day 的影响是随机的。为此，您将指定：

m2 <- lmer(Obs ~ Day + Treatment + Day:Treatment + (Day | Subject), mydata)

在这个模型中：

  • 如果在第 0 天治疗参考类别的预测得分 = 0，则截距
  • 天的系数是治疗参考类别的天数每增加 1 个单位的预测随时间变化
  • 治疗组的两个虚拟代码的系数（由 R 自动创建）是每个剩余治疗组与 Day=0 时的参考类别之间的预测差异
  • 两个交互项的系数是时间（天）对参考类别和其余治疗组之间预测分数的影响的差异。

截距和 Day 对分数的影响都是随机的（允许每个受试者在 Day=0 时具有不同的预测分数以及随时间的不同线性变化）。截距和斜率之间的协方差也在建模中（允许它们共同变化）。

如您所见，两个虚拟变量的系数的解释以 Day=0 为条件。他们会告诉您参考类别在 Day=0 时的预测分数是否与其余两个治疗组显着不同。因此，决定将 Day 变量居中的位置很重要。如果您以第 1 天为中心，则系数会告诉您第 1 天参考类别的预测分数是否与其余两个组的预测分数显着不同。这样，您可以查看组之间是否存在预先存在的差异。如果您以第 3 天为中心，则系数会告诉您第 3 天参考类别的预测分数是否与剩余两组的预测分数有显着差异。这样，您可以在干预结束时查看各组之间是否存在差异。

最后，请注意主题不嵌套在治疗中。您的三种治疗方法不是您想要将结果概括到的水平群体中的随机水平——相反，正如您所提到的，您的水平是固定的，并且您只想将您的结果概括到这些水平。（更不用说，如果您只有 3 个上层单元，则不应使用多级建模；参见 Maas & Hox, 2005。）相反，治疗是一个 2 级预测器，即在天数内采用单个值的预测器（一级单位）为每个科目。因此，它仅作为预测变量包含在您的模型中。

参考：
Maas, CJM, & Hox, JJ (2005)。为多级建模提供足够的样本量。方法论：欧洲行为和社会科学研究方法杂志，1，86-92。

其它你可能感兴趣的问题
上一篇什么是准二项分布（在 GLM 的背景下）？ 下一篇样本量为 1 时的统计检验