通过线性模型和广义线性模型(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
包)来近似计算这些区间。不过,这些函数的使用通常需要模型对象以及可能的其他参数,并且需要仔细考虑随机效应的影响。
这看起来相当熟悉,预测区间总是比置信区间大。
在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的运行可能需要一些时间。