使用R中的lme4进行重复测量数据的多层模型

Question

我有一个数据框架，其中包含大约200个人的帖子和后续测量。在这项研究中，我们试图找出运动参与与心理困扰症状之间是否存在相关性。我们有两个测量时段（后续和后续），这些测量是在一次关于健康和运动的研讨会后进行的。后续在研讨会后6个月进行，后续一年后进行。我们得出以下假设：“肥胖人士在研讨会后一年内参加运动与后续心理困扰症状显著正相关。” 我认为，因变量是心理困扰，而自变量是参加体育活动。 数据结构如下：

Df
$ measurement_period : Factor w/ 2 levels "0","1": 1 1 1 1
$   psychological_distress ; int 12 45 32 85
$  participation : Factor w/ 2 levels "0","1": 1 1 1 1
$  id : num 1 2 3 4

在阅读了这里的一些帖子后，我们认为模型中有两个层次：1）测量期（帖子和跟进）2）id

首先，我们使用以下代码进行了无条件模型（仅拦截模型以确认多级模型是否合适，希望这是正确的）：

test <-lmer(psychological_distress ~1+(1|id),data=Df

但我们不确定模型是否适合给定的数据结构，以及一级和二级分类是否正确。

非常感谢！

- Timon Doo

Answer 1

您的模型：

lmer(psychological_distress ~ 1 + (1|id) , data = Df)

这是一个方差分量模型。它会告诉您在心理困扰的变异中，有多少归因于id级别，以及有多少归因于单位/剩余级别。这并不能回答您的研究问题：

我们试图找出参与运动和痛苦症状之间是否存在相关性。

为了深入研究，您需要将参与变量作为固定效应，并且还要包括时间变量及其交互作用。因此，首先我建议您考虑以下内容：

lmer(psychological_distress ~ measurement_period*participation  + (1|id) , data = Df)

- Robert Long

---

谢谢您的回答。我有一个问题，不确定为什么您插入了参与和测量期之间的交互作用，这是必要的吗？提前感谢您。 - Timon Doo

不客气。这是一个很好的问题。在您的特定情况下，交互可能并不相关。我倾向于默认包含它，以防群组之间的关联存在差异。 - Robert Long

Answer 2

一个关于如何使用lme4拟合纵向和增长模型的好网站是https://rpsychologist.com/r-guide-longitudinal-lme-lmer。正如罗伯特所指出的，并在该网站上展示的那样，通常有必要适合“时间”和“组”（例如，治疗组与对照组）之间的交互作用，以查看每个组的结果如何随时间而变化。您可以通过查看系数来看到这种变化，但通常更容易绘制（调整后的）预测。以下是一个玩具示例：

library(parameters)
library(datawizard)
library(lme4)
library(ggeffects)
data("qol_cancer")

# filter two time points
qol_cancer <- data_filter(qol_cancer, time %in% c(1, 2))

# create fake treatment/control variable
set.seed(123)
treatment <- sample(unique(qol_cancer$ID), size = length(unique(qol_cancer$ID)) / 2, replace = FALSE)
qol_cancer$treatment <- 0
qol_cancer$treatment[qol_cancer$ID %in% treatment] <- 1

qol_cancer$time <- as.factor(qol_cancer$time)
qol_cancer$treatment <- factor(qol_cancer$treatment, labels = c("control", "treatment"))

m <- lmer(QoL ~ time * treatment + (1 + time | ID), 
          data = qol_cancer,
          control = lmerControl(check.nobs.vs.nRE = "ignore"))

model_parameters(m)
#> # Fixed Effects
#> 
#> Parameter                        | Coefficient |   SE |         95% CI | t(368) |      p
#> ----------------------------------------------------------------------------------------
#> (Intercept)                      |       70.74 | 2.15 | [66.52, 74.97] |  32.90 | < .001
#> time [2]                         |        0.27 | 2.22 | [-4.10,  4.64] |   0.12 | 0.905 
#> treatment [treatment]            |        4.88 | 3.04 | [-1.10, 10.86] |   1.60 | 0.110 
#> time [2] * treatment [treatment] |        1.95 | 3.14 | [-4.23,  8.13] |   0.62 | 0.535 
#> 
#> # Random Effects
#> 
#> Parameter                 | Coefficient
#> ---------------------------------------
#> SD (Intercept: ID)        |       15.14
#> SD (time2: ID)            |        7.33
#> Cor (Intercept~time2: ID) |       -0.62
#> SD (Residual)             |       14.33
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald t-distribution approximation.

ggpredict(m, c("time", "treatment")) |> plot()

关于交互项的统计显著性：总结中的p值可能会误导。如果您真正关心时间点之间或组别（治疗组与对照组）之间的统计学差异，建议计算包括p值的成对比较。您可以使用emmeans软件包来完成此操作。

library(emmeans)
emmeans(m, c("time", "treatment")) |> contrast(method = "pairwise", adjust = "none")
#>  contrast                          estimate   SE  df t.ratio p.value
#>  time1 control - time2 control       -0.266 2.22 186  -0.120  0.9049
#>  time1 control - time1 treatment     -4.876 3.04 186  -1.604  0.1105
#>  time1 control - time2 treatment     -7.092 2.89 316  -2.453  0.0147
#>  time2 control - time1 treatment     -4.610 2.89 316  -1.594  0.1118
#>  time2 control - time2 treatment     -6.826 2.73 186  -2.497  0.0134
#>  time1 treatment - time2 treatment   -2.216 2.22 186  -0.997  0.3199
#> 
#> Degrees-of-freedom method: kenward-roger

^{此内容由reprex软件包（v2.0.1）于2022-05-22创建。}

在这里，您可以看到，在时间点1，治疗组和对照组在其生活质量方面没有差异，但是在时间点2有所不同。

- Daniel

---

很好的例子（+1）:) - Robert Long

非常感谢你们的建议和示例，这对我理解多层次模型和我的数据有很大帮助。 - Timon Doo

所以，丹尼尔，在你的情况下，你的输出意味着时间作为 QoL 预测因素的影响不显著（p > .05），表明在两个时间段内 QoL 没有显著变化。关于时间和治疗之间的互动作用，你可以说它也不显著（p > .05），表明结果 QoL 在每组（治疗组和对照组）中随时间没有显著变化。我是对的吗？ - Timon Doo

我仍然没有理解条件生长模型中删除随机斜率（lmer(QoL ~ time * treatment + (1 | ID), data=data)）和添加随机斜率（lmer(QoL ~ time * treatment + (time | ID), data=data)）之间的区别。您如何知道在您的数据中可以去除随机斜率？或者您是通过指定（1 + time | ID）将其包含在代码中的？ - Timon Doo

我添加了一个关于统计显著差异的例子。我希望这能澄清我们在哪里找到QoL的显著变化，无论是随时间还是在不同组之间。 - Daniel

是的，这意味着我将"time"作为随机斜率添加了进去。对于纵向数据，假设结果可能发生变化的可能性因人而异是可行的。这时，您需要将"time"作为随机斜率包含在内。我认为一个很好的解释这个问题的网站是http://mfviz.com/hierarchical-models/。 - Daniel