R语言广义线性混合模型（GLMM）bootstrap预测置信区间可视化
R语言广义线性混合模型（GLMM）bootstrap预测置信区间可视化
通过线性模型和广义线性模型（GLM），预测函数可以返回在观测数据或新数据上预测值的标准误差（点击文末“阅读原文”获取完整代码数据）。
相关视频
然后，利用这些标准误差绘制出拟合回归线周围的置信区间或预测区间。置信区间（CI）的重点在于回归线，其可以解释为（假设我们绘制的是95%的置信区间）：“如果我们重复抽样X次，那么回归线将有95%的概率落在这个区间内”。另一方面，预测区间的重点在于单个数据点，其可以解释为（同样假设我们绘制的是95%的置信区间）：“如果我们在这些特定的解释变量值上抽样X次，那么响应值将有95%的概率落在这个区间内”。
对于广义线性混合模型（GLMM），预测函数不允许推导标准误差，原因是：“没有计算预测标准误差的选项，因为很难定义一种有效的方法来将方差参数中的不确定性纳入其中”。这意味着目前没有办法将拟合的随机效应标准差的估计（其估计值可能或多或少准确）纳入预测值标准误差的计算中。不过，我们仍然可以推导置信区间或预测区间，但需要注意，我们可能会低估估计值的不确定性。
library(lme4) # 加载lme4包，用于线性混合效应模型的分析
# 第一个案例：简单的线性混合效应模型，从10个组中模拟100个数据点，具有一个连续的固定效应变量
x <- runif(100, 0, 10) # 生成100个在0到10之间的均匀分布随机数，作为固定效应变量x
# 固定效应系数
fixed <- c(1, 0.5) # 设定固定效应系数，这里假设截距为1，x的系数为0.5
# 随机效应
rnd <- rnorm(10, 0, 0.7) # 生成10个来自均值为0、标准差为0.7的正态分布的随机数，作为随机效应
# 模型
# 使用类似predict的方法计算第一个置信区间和预测区间
newdat <- data.frame(x = seq(0, 10, length = 20)) # 生成一个新的数据框newdat，其中x是从0到10的等差序列，长度为20
这段代码是继续上面的线性混合效应模型（LMM）分析的，它计算了预测值、预测区间和置信区间，并使用bootMer函数进行了自助法（bootstrap）来估计置信区间。接下来，我会逐步解释这段代码的内容：
# 生成新数据框newdat的模型矩阵
mm <- model.matrix(~x, newdat)
# 根据固定效应计算新数据框的预测值
newdat$y <- mm %*% fixef(m)
# 使用vcov函数计算模型协方差矩阵，并使用tcrossprod计算其转置和原始矩阵的乘积
# 然后与模型矩阵mm相乘，得到预测值的方差
# 计算总方差，包括随机效应方差（此处仅考虑随机截距）
# 注意：这里假设随机效应只有一个，即随机截距，对于更复杂的模型需要调整
# 在newdat数据框中添加预测值、预测区间的下限和上限、置信区间的下限和上限
newdat <- data.frame(
newdat,
plo = newdat$y - 1.96 * sqrt(pvar1), # 预测区间的下限
# 第二版：使用bootMer进行自助法估计置信区间
# 定义一个函数，该函数应用于nsim次模拟，返回拟合值
predFun <- function(.) mm %*% fixef(.)
# 将自助法得到的置信区间的下限和上限添加到newdat数据框中
newdat$blo <- bb_se[1,]
# 绘制原始数据、拟合线、预测区间和置信区间
plot(y ~ x, data)
lines(newdat$x, newdat$y, col="red", lty=2, lwd=3)
# 添加图例
legend("topleft", legend=c("Fitted line", "Prediction interval", "Confidence interval", "Bootstrapped CI"),
这段代码的主要功能是：
- 使用模型矩阵和固定效应系数来计算新数据点的预测值。
- 计算预测值的方差（
pvar1），进而得到预测区间。 - 计算包含随机效应方差的总方差（
tvar1），进而得到置信区间。 - 使用
bootMer函数进行自助法抽样，估计置信区间。 - 最后，绘制原始数据、拟合线、预测区间和置信区间。
需要注意的是，这段代码假设随机效应只有一个随机截距。对于包含其他类型随机效应的模型，计算总方差时需要相应地进行调整。此外，bootMer函数可能需要较长时间来执行，特别是当模型复杂或自助法抽样次数较多时。
在上述代码中，模拟数据的生成和模型的拟合都是基于线性混合效应模型（LMM）的。然而，计算置信区间（CI）和预测区间（PI）的部分并没有给出具体的实现，因为对于线性混合效应模型，这些区间的计算通常比线性模型更复杂。通常，我们会使用自助法（bootstrap）或者基于模型的近似方法来估计这些区间。在R中，可以使用bootMer函数（来自lme4包）或predictInterval函数（来自merTools包）来近似计算这些区间。不过，这些函数的使用通常需要模型对象以及可能的其他参数，并且需要仔细考虑随机效应的影响。
这看起来相当熟悉，预测区间总是比置信区间大。
01
02
03
04
在predict.merMod函数的帮助页面中，lme4包的作者写道，bootMer应该是从广义线性混合模型（GLMM）推导置信区间的首选方法。那里的想法是从模型中模拟N次新数据，然后获取一些感兴趣的统计数据。在我们的案例中，我们感兴趣的是通过推导自举拟合值来获取回归线的置信区间。bb$t是一个矩阵，其中列是观测值，行是不同的自举样本。为了得到拟合线的95%置信区间，我们需要获取排序后的自举值的[0.025N,0.975N]范围的值。
即使对每个自举样本都计算了新的随机效应值（因为bootMer中默认use.u=FALSE），自举的置信区间也非常接近“正常”的置信区间。
现在让我们转向一个更复杂的例子，一个具有两个交叉随机效应的泊松广义线性混合模型（Poisson GLMM）：
# 第二个案例，具有两个交叉随机效应和泊松响应的更复杂设计
m <- glmer(y ~ x + (1|f1) + (1|f2), data, family = "poisson")
# 对于广义线性混合模型，我们需要在添加/删除标准误后反变换预测值
newdat <- data.frame(x = seq(0, 10, length = 20))
tvar1 <- pvar1 + VarCorr(m)$f1[1] + VarCorr(m)$f2[1] # 对于更复杂的模型，必须进行调整
# 使用bootMer的第二个版本
bb <- bootMer(m, FUN = predFun, nsim = 200)
# 绘图
plot(y ~ x, data)
lines(newdat$x, newdat$y, col = "red", lty = 2, lwd = 3)
同样，在这种情况下，自助法置信区间（Bootstrapped CI）与“常规”置信区间非常接近。在这里，我们已经看到了三种不同的方法来推导表示回归线（CI）和响应点（PI）周围不确定性的区间。选择哪种方法取决于您想看到什么（我拟合的线的周围不确定性的程度，或者如果我抽样新的观测值，它们会取什么值），以及复杂模型的计算能力，因为对于具有许多观测值和复杂模型结构的广义线性混合模型（GLMM），bootMer的运行可能需要一些时间。