R语言广义线性混合模型(GLMM)bootstrap预测置信区间可视化

简介: R语言广义线性混合模型(GLMM)bootstrap预测置信区间可视化


通过线性模型和广义线性模型(GLM),预测函数可以返回在观测数据或新数据上预测值的标准误差点击文末“阅读原文”获取完整代码数据

相关视频

image.png

然后,利用这些标准误差绘制出拟合回归线周围的置信区间或预测区间。置信区间(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"),

这段代码的主要功能是:

  1. 使用模型矩阵和固定效应系数来计算新数据点的预测值。
  2. 计算预测值的方差(pvar1),进而得到预测区间。
  3. 计算包含随机效应方差的总方差(tvar1),进而得到置信区间。
  4. 使用bootMer函数进行自助法抽样,估计置信区间。
  5. 最后,绘制原始数据、拟合线、预测区间和置信区间。

需要注意的是,这段代码假设随机效应只有一个随机截距。对于包含其他类型随机效应的模型,计算总方差时需要相应地进行调整。此外,bootMer函数可能需要较长时间来执行,特别是当模型复杂或自助法抽样次数较多时。

在上述代码中,模拟数据的生成和模型的拟合都是基于线性混合效应模型(LMM)的。然而,计算置信区间(CI)和预测区间(PI)的部分并没有给出具体的实现,因为对于线性混合效应模型,这些区间的计算通常比线性模型更复杂。通常,我们会使用自助法(bootstrap)或者基于模型的近似方法来估计这些区间。在R中,可以使用bootMer函数(来自lme4包)或predictInterval函数(来自merTools包)来近似计算这些区间。不过,这些函数的使用通常需要模型对象以及可能的其他参数,并且需要仔细考虑随机效应的影响。

441614909142cc8308af2e6ab68a1144.png

这看起来相当熟悉,预测区间总是比置信区间大。

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)


4d2c19a362553a4cbdcc7623da97e58b.png

同样,在这种情况下,自助法置信区间(Bootstrapped CI)与“常规”置信区间非常接近。在这里,我们已经看到了三种不同的方法来推导表示回归线(CI)和响应点(PI)周围不确定性的区间。选择哪种方法取决于您想看到什么(我拟合的线的周围不确定性的程度,或者如果我抽样新的观测值,它们会取什么值),以及复杂模型的计算能力,因为对于具有许多观测值和复杂模型结构的广义线性混合模型(GLMM),bootMer的运行可能需要一些时间。

 

相关文章
|
1月前
|
数据可视化 数据挖掘 API
【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
|
1月前
|
Web App开发 数据可视化 数据挖掘
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
|
1月前
|
数据采集 数据可视化
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
|
1月前
|
移动开发 数据可视化
广义线性模型beta二项分布的淋巴结疾病风险预测可视化R语言2实例合集|附数据代码
广义线性模型beta二项分布的淋巴结疾病风险预测可视化R语言2实例合集|附数据代码
|
1月前
|
机器学习/深度学习 数据可视化
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为2
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为
|
1月前
利用R语言进行典型相关分析实战
利用R语言进行典型相关分析实战
|
1月前
|
机器学习/深度学习 算法
R语言分类回归分析考研热现象分析与考研意愿价值变现
R语言分类回归分析考研热现象分析与考研意愿价值变现
|
1月前
|
数据可视化 定位技术
R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化
R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化
|
1月前
|
机器学习/深度学习 数据可视化 算法
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为1
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为
|
1月前
|
机器学习/深度学习 数据可视化 算法
R语言聚类分析、因子分析、主成分分析PCA农村农业相关经济指标数据可视化|数据分享
R语言聚类分析、因子分析、主成分分析PCA农村农业相关经济指标数据可视化|数据分享

热门文章

最新文章