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的运行可能需要一些时间。

 

相关文章
|
3月前
|
数据可视化 数据挖掘 图形学
R语言基础可视化:使用ggplot2构建精美图形的探索
【8月更文挑战第29天】 `ggplot2`是R语言中一个非常强大的图形构建工具,它基于图形语法提供了一种灵活且直观的方式来创建各种统计图形。通过掌握`ggplot2`的基本用法和美化技巧,你可以轻松地将复杂的数据转化为直观易懂的图形,从而更好地理解和展示你的数据分析结果。希望本文能够为你探索`ggplot2`的世界提供一些帮助和启发。
|
3月前
|
数据可视化 数据挖掘 数据处理
R语言高级可视化技巧:使用Plotly与Shiny制作互动图表
【8月更文挑战第30天】通过使用`plotly`和`shiny`,我们可以轻松地创建高度互动的数据可视化图表。这不仅增强了图表的表现力,还提高了用户与数据的交互性,使得数据探索变得更加直观和高效。本文仅介绍了基本的使用方法,`plotly`和`shiny`还提供了更多高级功能和自定义选项,等待你去探索和发现。希望这篇文章能帮助你掌握使用`plotly`和`shiny`制作互动图表的技巧,并在你的数据分析和可视化工作中发挥更大的作用。
|
6月前
|
数据可视化 数据挖掘 API
【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
|
3月前
|
数据可视化
R语言可视化设计原则:打造吸引力十足的数据可视化
【8月更文挑战第30天】R语言可视化设计是一个综合性的过程,需要综合运用多个设计原则来创作出吸引力十足的作品。通过明确目标、选择合适的图表类型、合理运用色彩与视觉层次、明确标注与引导视线以及引入互动性与动态效果等原则的应用,你可以显著提升你的数据可视化作品的吸引力和实用性。希望本文能为你提供一些有益的启示和帮助。
|
6月前
|
数据采集 数据可视化
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
|
6月前
|
Web App开发 数据可视化 数据挖掘
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
|
6月前
|
移动开发 数据可视化
广义线性模型beta二项分布的淋巴结疾病风险预测可视化R语言2实例合集|附数据代码
广义线性模型beta二项分布的淋巴结疾病风险预测可视化R语言2实例合集|附数据代码
|
2月前
|
数据采集 机器学习/深度学习 数据可视化
R语言从数据到决策:R语言在商业分析中的实践
【9月更文挑战第1天】R语言在商业分析中的应用广泛而深入,从数据收集、预处理、分析到预测模型构建和决策支持,R语言都提供了强大的工具和功能。通过学习和掌握R语言在商业分析中的实践应用,我们可以更好地利用数据驱动企业决策,提升企业的竞争力和盈利能力。未来,随着大数据和人工智能技术的不断发展,R语言在商业分析领域的应用将更加广泛和深入,为企业带来更多的机遇和挑战。
|
8天前
|
数据挖掘 C语言 C++
R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。
【10月更文挑战第21天】时间序列分析是一种重要的数据分析方法,广泛应用于经济学、金融学、气象学、生态学等领域。R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。本文将介绍使用R语言进行时间序列分析的基本概念、方法和实例,帮助读者掌握R语言在时间序列分析中的应用。
27 3
|
6月前
|
机器学习/深度学习 数据可视化
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为2
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为