STAT315 Week 7 线性混合模型 LMM

STAT315 第七周线性混合模型 LMM 笔记，欢迎交流。接有偿咨询&辅导。

上图是三种不同的线性混合模型：

仅随机截距（Intercept only）

定义：
在这种情况下，每个群组（如每个被试）的随机效应只影响截距，而斜率是固定的。这意味着不同群组的回归线平行但在y轴上的起点不同。

数学表示：
$$y_{ij}={\beta }_0+{\beta }_1{x}_{ij}+u_{0i}+{ϵ}_{ij}$$

其中：

$u_{0i}$ 是第 $i$ 个群组的随机截距。

${ϵ}_{ij}$ 是观测误差。

R代码示例：

model_intercept_only <- lme(response ~ time, random = ~ 1 | subject, data = data)
summary(model_intercept_only)

仅随机斜率（Random Slope only）

定义：
在这种情况下，每个群组的随机效应只影响斜率，而截距是固定的。这意味着不同群组的回归线从同一点开始但斜率不同。

数学表示：
$$y_{ij}={\beta }_0+({\beta }_1+u_{1i})x_{ij}+{ϵ}_{ij}$$

其中：

$u_{1i}$ 是第 $i$ 个群组的随机斜率。

${ϵ}_{ij}$ 是观测误差。

R代码示例：

model_slope_only <- lme(response ~ time, random = ~ time - 1 | subject, data = data)
summary(model_slope_only)

随机截距和随机斜率（Intercept and Random Slope）

定义：
在这种情况下，每个群组的随机效应同时影响截距和斜率。这意味着不同群组的回归线既可以在y轴上不同起点，又可以具有不同的斜率。

数学表示：
$$y_{ij}=({\beta }_0+u_{0i})+({\beta }_1+u_{1i})x_{ij}+{ϵ}_{ij}$$

其中：

$u_{0i}$ 是第 $i$ 个群组的随机截距。

$u_{1i}$ 是第 $i$ 个群组的随机斜率。

${ϵ}_{ij}$ 是观测误差。

R代码示例：

model_intercept_and_slope <- lme(response ~ time, random = ~ time | subject, data = data)
summary(model_intercept_and_slope)

解释图示

图示解释了三种不同随机效应模型的区别：

左图：仅随机截距。不同颜色的点表示不同群组，回归线平行但截距不同。

中图：仅随机斜率。不同颜色的点表示不同群组，回归线从同一点开始但斜率不同。

右图：随机截距和随机斜率。不同颜色的点表示不同群组，回归线既有不同的截距又有不同的斜率。

这张幻灯片介绍了线性混合效应模型（LMM）中的估计方法，特别是最大似然估计（MLE）和限制最大似然估计（REML）。让我们详细讨论这两种估计方法及其在R语言中的应用。

最大似然估计（Maximum Likelihood Estimation, MLE）

定义：
最大似然估计是通过最大化似然函数来估计模型参数的方法。MLE广泛应用于固定效应模型，因为它满足某些统计要求，如适用于非正态分布等。

优点：

MLE是一种通用方法，适用于广泛的模型和分布。

在大样本情况下，MLE估计量具有一致性和渐近正态性。

缺点：

在混合效应模型中，MLE估计量可能有偏，特别是方差成分的估计。这种偏差与因素的水平数量有关（水平数量越少，偏差越大）。

限制最大似然估计（Restricted Maximum Likelihood Estimation, REML）

定义：
限制最大似然估计是一种改进的MLE方法，通过最大化剔除固定效应后的残差的似然函数来估计模型参数。

优点：

REML一般比MLE在估计方差成分时产生更小的偏差。

特别适用于混合效应模型，因为它考虑了固定效应对方差成分估计的影响。

缺点：

REML在估计固定效应时仍可能有偏差。

在R语言中的应用

在R语言中，可以使用nlme包或lme4包来进行MLE和REML估计。默认情况下，这些包使用REML进行估计，但可以通过参数设置来选择MLE。

使用nlme包

# 安装并加载nlme包
install.packages("nlme")
library(nlme)

# 构建数据示例
data <- data.frame(
subject = rep(1:10, each = 5),
time = rep(1:5, times = 10),
response = rnorm(50)
)

# 使用REML进行估计（默认）
model_reml <- lme(response ~ time, random = ~ 1 | subject, data = data)
summary(model_reml)

# 使用MLE进行估计
model_mle <- lme(response ~ time, random = ~ 1 | subject, data = data, method = "ML")
summary(model_mle)

使用lme4包

# 安装并加载lme4包
install.packages("lme4")
library(lme4)

# 使用REML进行估计（默认）
model_reml <- lmer(response ~ time + (1 | subject), data = data)
summary(model_reml)

# 使用MLE进行估计
model_mle <- lmer(response ~ time + (1 | subject), data = data, REML = FALSE)
summary(model_mle)

总结

MLE：广泛应用，但在混合效应模型中可能产生偏差，特别是方差成分的估计。

REML：针对MLE的偏差问题进行改进，通常更适用于混合效应模型。

R中的应用：nlme和lme4包都可以进行MLE和REML估计，默认情况下使用REML。

诊断图

1. Q-Q图（Quantile-Quantile Plot）

目的：检查残差的正态性。
方法：将模型的残差与标准正态分布的分位数进行比较。如果残差是正态分布的，点应大致沿对角线排列。

2. 残差图（Residual Plot）

目的：检查残差的均值是否为零，以及残差是否均匀分布（即无模式）。
方法：绘制拟合值与残差的散点图。理想情况下，残差应随机分布在零附近，不显示任何模式。

这张幻灯片讨论了线性混合效应模型（LMM）的推断问题，特别是p值、t统计量和预测随机效应的问题。以下是关键点的详细解释：

1. p值的缺失

在使用lme4包的summary()函数进行LMM分析时，p值通常不会给出。这是因为在混合效应模型中，F统计量存在问题：

F统计量的问题：通常情况下，F统计量用于测试固定效应项的显著性，但在混合效应模型中，使用anova()函数会假设随机效应的估计值是已知的真值，这不合适。

过小的p值：由于随机效应实际上是估计的，F统计量不完全符合F分布，因此计算出的p值往往过小。

因此，lme4包的开发者选择省略p值，认为没有答案比可能的错误答案更好。

2. t统计量

t统计量是根据F统计量的平方根计算的，因此同样存在问题：

不适当的结果：lme4包中提供的t统计量在检测固定效应显著性时结果可能不可靠。

3. 置信区间

置信区间是一种在缺乏p值的情况下进行推断的有用方法：

轮廓似然方法（Profile Likelihood Methods）：建议使用轮廓似然方法来获得系数的置信区间。

4. 随机效应的检验

设计特征：随机效应通常作为设计特征包含在实验中，因此不应像固定效应那样进行显著性检验。

未来实验设计：虽然随机效应不适合显著性检验，但可以用于设计未来实验的信息。

5. 预测随机效应

混合效应模型的主要问题之一是预测随机效应：

随机效应不是参数：因为模型中的随机效应是随机变量而不是固定参数，所以不应直接估计它们。

期望值：可以通过查看期望值来处理随机效应，或者从贝叶斯角度出发，查看后验均值。

R中的函数：在lme4包中，有一个名为ranef(modelname)的函数可以用来获得预测的随机效应。

在R中的应用

让我们看一下在R中如何进行这些推断和诊断。

使用lme4包进行LMM推断

1. 拟合模型

library(lme4)

# 创建示例数据
set.seed(123)
data <- data.frame(
subject = rep(1:10, each = 5),
time = rep(1:5, times = 10),
response = rnorm(50)
)

# 拟合线性混合效应模型
model <- lmer(response ~ time + (1 | subject), data = data)
summary(model)

2. 置信区间

# 使用confint函数计算置信区间
confint(model)

3. 提取随机效应

# 提取随机效应
random_effects <- ranef(model)
print(random_effects)

总结

• p值和t统计量：在LMM中，由于F统计量的问题，p值和t统计量的可靠性受到影响，因此通常省略p值。
• 置信区间：使用轮廓似然方法获得固定效应的置信区间。
• 随机效应的检验：不建议对随机效应进行显著性检验，但可以用于未来实验设计的信息。
• 预测随机效应：使用ranef函数提取随机效应的估计值。