本文探讨了混合线性模型在分析纵向数据中的应用，关注时间效应、群体效应以及它们之间的交互作用。模型包括随机效应和随机截距，通过实例展示了如何检验这些效应，以及如何使用参数曲线处理非线性关系和预测个体响应。
目录
Interaction between time and group 交互
Testing the interaction between group and time
Graphical presentation of the fitted model 作图展示
前言
下面正式进入纵向数据的分析，从前文我们知道时间（time）是最重要的影响因素，其次就是分组，这里分组指的性别，推而广之，如果是临床资料，分组可能是不同的研究中心，不同的药物等等。
在正式进入下一站前，我们需要思考几个有意思的问题：
原文：In the next sections, we are going to answer the following three main questions:
- Time effect: What is the shape of the trajectory of the mean response over time?
- Group effect: What is the average difference between groups of individuals?
- Interaction between time and group: How the relationship between the response and time vary according to groups of individuals?
时间效应：随着时间推移，平均响应轨迹的形状是什么？ 群体效应：不同个体群体之间的平均差异是多少？ 时间和群体之间的交互作用：响应和时间之间的关系如何随着个体群体的变化而变化？
也就是时间，分组，时间与分组
混合线性模型LMM拟合
初始
lin_0 <- lmer(distance ~ 1 + (1 | id), data = dental_long)
summary(lin_0)这是第一个模型，模型中id作为随机截距，也就是说我们认为基线存在个体差异影响response。
ci.lin(lin_0) 
牙距的估计边际均值是24.02315mm。随机效应反映受试者变异性的估计方差为3.752；反映受试者内变异性的误差项估计方差为4.930。任意两次重复测量之间的相关性（ICC）为3.76/（3.76+4.93）= 0.43.
可以使用lmer test中的ranova函数得出的类似ANOVA表格来对方差分量进行检验。
ranova(lin_0)
小的 P 值表明存在个体间的异质性，这支持选择混合效应模型而不是仅固定效应模型的证据。
Time effect 时间
Question: Is the mean response varying with time?
一个建模策略是使用指示变量来对不同时间点的平均响应进行对比。在牙齿生长的例子中，我们测量了8岁（基线）、10岁、12岁和14岁时的响应值：
lin_age <- lmer(distance ~ measurement + (1 | id), data = dental_long)
summary(lin_age) 在这个模型里面我们拟合的时间效应作为固定效应，id仍然作为随机截距
8岁儿童的平均牙距为22.18毫米。 10岁和8岁儿童平均响应的差异是β1=0.98mm，12岁和8岁差异是2.46mm，14岁和8岁差异是3.90mm。
emmeans 函数可以用来计算边际平均值及其相应的置信区间，非常有用。
tidy(emmeans(lin_age, "measurement"), conf.int = TRUE)
我们可以通过检验空假设来测试平均响应在不同时间是否保持恒定，即用于建模时间的所有回归系数同时等于零。
Anova(lin_age) 
Wald检验结果表明平均响应随时间变化的轨迹并不是平坦的。我们可以使用predict函数来展示拟合模型的估计轨迹。
dental_fit <- bind_cols(
dental_long, pred_age = predict(lin_age, re.form = ~ 0)
)
ggplot(dental_fit, aes(age, distance)) +
geom_line(aes(group = factor(id))) +
geom_point(aes(y = pred_age), col = "blue", size = 2) +
labs(x = "Age, years", y = "Dental growth, mm") 
Group effect 分组
Question: Is the mean response varying in the two groups of individuals?
这次分组变量：性别作为固定效应。
lin_agesex <- lmer(distance ~ measurement + sex + (1 | id), data = dental_long)
summary(lin_agesex)
ci.lin(lin_agesex)
男孩和女孩的平均牙距差异为2.3毫米（95%置信区间为0.83-3.81毫米），调整了时间因素。emmeans 可以用于估计不同性别和时间测量组合的平均响应值。
tidy(emmeans(lin_agesex, c("measurement", "sex")), conf.int = TRUE) 
Interaction between time and group 交互
Question: Is the change of the mean response over time varying according to group of individuals?
lin_agesexinter <- lmer(distance ~ measurement*sex + (1 | id), data = dental_long)
lin_agesexinter交互作用项作为固定效应，此处不需要单独列出sex和measurement，R默认都有。
女孩在8岁时的平均牙距为21.2毫米 女孩在10岁时的平均牙距比8岁时高出1.04毫米 女孩在14岁时的平均牙距比8岁时高出2.9毫米 男孩在8岁时的平均牙距为21.2+1.7= 22.9 毫米 男孩在10岁时的平均牙距比8岁时高出0.93毫米
Testing the interaction between group and time
Anova(lin_agesexinter, type = 3)
Wald检验检测到年龄和性别之间存在一些交互作用（ = 0.07），这表明平均响应曲线可能不是平行的。 可以得到性别和测量时间组合的边际预测值，如下所示：
tidy(emmeans(lin_agesexinter, c("measurement", "sex")), conf.int = TRUE)
男孩和女孩之间的平均响应差异取决于年龄。 随时间平均响应的变化取决于性别。
图形展示模型
dental_fit$pred_agesexinter <- predict(lin_agesexinter, re.form = ~ 0)
ggplot(dental_fit, aes(age, distance)) +
geom_line(aes(group = factor(id))) +
geom_point(aes(y = pred_agesexinter, col = sex), size = 2) +
labs(x = "Age, years", y = "Dental growth, mm", col = "Sex")
参数曲线 Parametric curve 线性 非线性
在混合线性效应模型中，参数曲线通常用于建模因变量和一个或多个预测变量之间的复杂非线性关系。这在一些研究领域中可能非常重要，特别是当线性模型无法恰当地捕捉到数据的特征时。参数曲线可以更好地适应数据，并提供更准确的预测和解释。
在许多情况下，特别是涉及生物医学、心理学和社会科学等领域的研究中，混合线性效应模型和参数曲线的组合可以帮助建模复杂的数据结构，并更全面地理解数据中的变化。因此，可以确定在混合线性效应模型中，参数曲线在某些情境下可能会非常重要。
换而言之，我们要进行线性及非线性的判断，首先构建模型
lin_agecsexinter <- lmer(distance ~ sex*age + (1 | id), data = dental_long)
summary(lin_agecsexinter) 
我们来看一下模型参数的实际意义：age，sex，age*sex交互项
女孩0岁时的平均距离为17.37mm，0岁时，男孩和女孩牙距相差-1.03mm，对于女孩群体，平均响应每增加一岁会增加0.48。0.3代表着性别回归系数在特定年龄取值时的额外变化，反之亦然。
0.3代表着性别在特定年龄取值时对回归系数的额外影响。这意味着性别对因变量的影响不仅取决于它本身的主效应，还受到年龄的影响而发生变化。例如，如果年龄对性别的影响产生了变化，那么0.3就代表了这种变化的额外影响。换句话说，性别对因变量的影响程度并非固定不变，而是随着年龄变化而有所不同，并且0.3表示了这种差异的额外改变。也就是说性别对因变量的影响需要额外加上0.3来反映在特定年龄取值时的影响。也就是说一个14岁的女孩计算牙间距，需要考虑性别和年龄对牙齿的影响，即性别对牙齿的系数再加上0.3.
Test linear hypotheses 线性假设检验
K <- rbind(
c(0, 1, 0, 8),
c(0, 1, 0, 14),
c(0, 0, 1, 1),
c(0, 0, 2, 2),
c(1, 1, 8, 8),
c(1, 1, 14, 14),
c(0, 0, 14-8, 14-8)
)
rownames(K) <- c("Difference between boys and girls at age 8 years",
"Difference between boys and girls at age 14 years",
"Linear trend for 1 year increase of age among boys",
"Linear trend for 2 years increase of age among boys",
"Mean response at age 8 years among boys",
"Mean response at age 14 years among boys",
"Difference between age 14 and 8 among boys")
tidy(glht(lin_agecsexinter, linfct = K), conf.int = TRUE) %>%
gt() %>%
fmt_number(columns = -1,decimals = 2)这段R代码建立了一个系数矩阵K，然后利用glht函数进行多重假设检验，然后使用tidy函数进行结果整理并使用gt函数对结果进行排版输出。接下来我们来逐步解释：
首先，使用rbind函数创建了一个系数矩阵K，每一行代表一个系数向量，总共有7行。每一行包含4个系数，分别代表不同的假设。使用rownames函数为系数矩阵K添加了行名，分别对应不同的假设。接下来使用glht函数对混合效应模型lin_agecsexinter进行多重比较假设检验。linfct参数指定了要进行的线性假设检验的系数线性组合。紧接着，使用tidy函数对多重比较假设检验的结果进行整理，并设置conf.int参数为TRUE以获得置信区间。最后，使用管道符（%>%）将整理后的结果传递给gt函数进行格式化输出。gt函数用于利用gt包对表格数据进行格式化美化，fmt_number函数调整数字格式为两位小数。
Graphical presentation of the fitted model 作图展示
pred_agecsexinter <- expand.grid(
age = seq(8, 14, .5),
sex = levels(dental_long$sex)
) %>%
bind_cols(pred = predict(lin_agecsexinter, newdata = ., re.form = ~ 0))
ggplot(pred_agecsexinter, aes(age, pred, col = sex)) +
geom_line() +
labs(x = "Age, years", y = "Dental growth, mm", col = "Sex")这段R代码首先使用expand.grid函数创建了一个预测数据框，其中age的取值范围为8到14岁，以0.5为间隔，sex取值为dental_long数据集中sex变量的水平。
接着使用管道符（%>%）将预测数据框传递给bind_cols函数，并通过predict函数基于混合效应模型lin_agecsexinter进行预测，得到的预测值被添加到预测数据框中。
之后，利用ggplot函数创建了一个散点图，其中预测数据框pred_agecsexinter中age作为x轴，pred作为y轴，sex用颜色进行区分。geom_line函数添加了连接预测值的曲线，labs函数用于添加轴标签和图例标题。最终绘制了预测的结果图。
我们可以看到一个线性关系。
Prediction random effects
在许多应用程序中，推断着重于固定效应（随时间的整体响应变化）。然而，我们也可以估计（或预测）随时间变化的个体特定响应轨迹。
这被称为“最佳线性无偏预测器”（或BLUP）。它们也被称为经验贝叶斯预测。其思想是根据响应的实际观察值，预测每个个体的随机效应（与总体均值的偏离）。对于第i个个体，预测的BLUP响应轨迹是估计的总体平均响应轮廓和主体观察到的响应轮廓的加权组合。这解释了为什么经验BLUP通常被说成是将第i个主体的预测响应轮廓“缩小”到总体平均值。
假设一个只有年龄（连续）作为协变量的随机截距模型。
lin_agec <- lmer(distance ~ age + (1 | id), data = dental_long)
lin_agec
特定个体，我们选择id10和id21
sid <- c(10, 21)
expand.grid(
age = seq(8, 14, .5),
id = sid
) %>%
bind_cols(
indiv_pred = predict(lin_agec, newdata = .),
marg_pred = predict(lin_agec, newdata = ., re.form = ~ 0)
) %>%
left_join(
filter(dental_long, id %in% sid), by = c("id", "age")
) %>%
ggplot(aes(age, indiv_pred, group = id, col = factor(id))) +
geom_line() +
geom_point(aes(y = distance)) +
geom_line(aes(y = marg_pred, col = "Marginal"), lwd = 1.5) +
labs(x = "Age, years", y = "Dental growth, mm", col = "Curve")
上图提供了随机截距模型的图形表示。 人口中随时间的边际（整体）平均响应随年龄呈线性增长（以实线粗蓝线表示）。 两个特定个体的条件平均响应以不同颜色的细实线表示。
随机效应和随机截距同时存在
同样的，时间/年龄 都可以作为随机斜率，这增加了模型的复杂性，解释起来，更加复杂。
lin_agecr <- lmer(distance ~ age + (age | id), data = dental_long)
summary(lin_agecr) 
这个模型假设个体在响应的基线水平（截距）上不仅有所不同，而且在随时间变化的平均响应上也有所不同。
VarCorr(lin_agecr)
as.data.frame(VarCorr(lin_agecr))
随机截距和斜率之间存在负相关性（-0.32）：基线（8年）处的平均响应越大，随时间变化的增长速率可能就越低。我们可以得到边际（整体）和个体特定的平均响应轨迹的预测值。
expand.grid(
age = seq(8, 14, .5),
id = unique(dental_long$id)
) %>%
bind_cols(
indiv_pred = predict(lin_agecr, newdata = .),
marg_pred = predict(lin_agecr, newdata = ., re.form = ~ 0)
) %>%
ggplot(aes(age, indiv_pred, group = id)) +
geom_line(col = "grey") +
geom_line(aes(y = marg_pred), col = "blue", lwd = 1.5) +
labs(x = "Age, years", y = "Dental growth, mm")
随机斜率中的分组变量
分组变量：性别；
lin_agecsexinterr <- lmer(distance ~ age*sex + (age | id), data = dental_long)
summary(lin_agecsexinterr) 
在预测牙距时存在性别和时间之间显著的交互作用（交互作用的P值为0.032）。
作图展示模型变化
expand.grid(
age = seq(8, 14, .5),
id = unique(dental_long$id)
) %>%
left_join(select(dental, sex, id), by = "id") %>%
bind_cols(
indiv_pred = predict(lin_agecsexinterr, newdata = .),
marg_pred = predict(lin_agecsexinterr, newdata = ., re.form = ~ 0)
) %>%
ggplot(aes(age, indiv_pred, group = id)) +
geom_line(col = "grey") +
geom_line(aes(y = marg_pred, col = sex), lwd = 1.5) +
labs(x = "Age, measurements", y = "Mean Dental growth, mm")
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
