如何提取lmer随机效应的置信区间及dotplot(ranef())绘制的截距值
lmer随机效应点估计与置信区间提取方法
你之前在fm1对象中找不到对应值,是因为随机效应的置信区间相关数据没有存储在原始模型对象的主列表中,而是需要通过带condVar=TRUE参数的ranef()返回对象的属性获取,具体操作如下:
步骤1:生成带条件方差的随机效应对象
首先运行ranef()时必须携带condVar=TRUE参数,才会输出计算置信区间所需的后验方差:
library(lme4) data(sleepstudy) fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) # 生成带后验方差的随机效应对象 re_obj <- ranef(fm1, condVar = TRUE)
步骤2:提取随机效应点估计值
re_obj$Subject中直接存储了每个受试者的随机截距和*随机斜率(Days)*的点估计值:
# 查看随机效应点估计 head(re_obj$Subject)
步骤3:提取后验方差并计算95%置信区间
后验方差存储在re_obj$Subject的postVar属性中,是一个三维数组,对角线元素为对应效应的方差,开平方得到标准误后即可计算正态近似的95%置信区间,和dotplot()绘制的区间完全一致:
# 提取后验方差数组 post_var <- attr(re_obj$Subject, "postVar") # 整合所有结果为数据框 ci_result <- data.frame( Subject = rownames(re_obj$Subject), # 随机截距相关指标 intercept_est = re_obj$Subject[["(Intercept)"]], intercept_ci_low = re_obj$Subject[["(Intercept)"]] - 1.96 * sqrt(post_var[1,1,]), intercept_ci_high = re_obj$Subject[["(Intercept)"]] + 1.96 * sqrt(post_var[1,1,]), # 随机斜率(Days)相关指标 days_est = re_obj$Subject[["Days"]], days_ci_low = re_obj$Subject[["Days"]] - 1.96 * sqrt(post_var[2,2,]), days_ci_high = re_obj$Subject[["Days"]] + 1.96 * sqrt(post_var[2,2,]) ) # 查看结果 head(ci_result)
可选简化方案
如果不需要手动计算,可以直接调用broom.mixed包的tidy()函数一键导出:
library(broom.mixed) re_ci <- tidy(fm1, effects = "ran_vals", conf.int = TRUE) head(re_ci)
返回结果自动包含所有随机效应的估计值、标准误、95%置信上下限。
内容的提问来源于stack exchange,提问作者Luis M. García
相关产品推荐
相关产品推荐

