如何从一个 lmer 模型中提取所有非参考水平连续变量系数的置信区间？

如何解决如何从一个 lmer 模型中提取所有非参考水平连续变量系数的置信区间？

我正在使用 lmer() 拟合 y ~ v + m + s + m:s + (1|subunit) 模型。 s 是一个与 m 相互作用的连续变量，一个具有 3 个级别的分类因子：A、B 和 C >.

拟合模型使用A作为因子m的参考水平：

fit_ref_A <- lmer(y ~ 1 + v + m*s + (1|subunit),data=df)

B 和 C 的参数估计值可以仅使用 A 的估计值和 B 的估计值来计算strong> 和 C。我对提取置信区间感兴趣。

运行 confint() 给出变量 s 在 A 处的斜率的置信区间。我还对 B 和 C 处 s 的斜率的置信区间感兴趣，而不仅仅是 C 处的斜率差异的置信区间强>A。有没有办法从 fit_ref_A 中提取它？到目前为止，我唯一能弄清楚的是使用 B 作为参考重新调平，适应新的配合，然后使用 C 作为参考重新调平，并且适合第三次。

问题：有没有办法从 fit_ref_A 中提取所有内容（尤其是置信区间）？

代码：

library(lme4)
# create the dataset,unbalanced at the lowest stratum ( 2 repeats for m==A instead of 3)
set.seed(2)
s_levels <- 1:5
m_levels <- c("A","B","C")
v_levels <- c("L2","L3","L4")
reps <- 1:3
df <- expand.grid(rep=reps,s=s_levels,m=m_levels,v=v_levels)
df$subunit <- as.factor(paste(df$v,"-",df$m,df$s,sep=""))
df$y <- rnorm(nrow(df),1)
df <- subset(df,!(rep==3 & m=="A"))  # drop the 3rd repeat for m=="A"
table(df$m)  # shows 30 for A,45 for B,45 for C as expected

# fit 3 different models,with three different reference levels for 'm'
fit_ref_A <- lmer(y ~ 1 + v + m*s + (1|subunit),data=df)

df$m <- relevel(df$m,ref = "B")
fit_ref_B <- lmer(y ~ 1 + v + m*s + (1|subunit),ref = "C")
fit_ref_C <- lmer(y ~ 1 + v + m*s + (1|subunit),data=df)

# Calculate the confidence intervals for the continuous variable s at the three
# different levels for categorical factor m.  Must use 3 separate fits.

cis_at_m_reference_A <- confint(fit_ref_A)
cis_at_m_reference_B <- confint(fit_ref_B)
cis_at_m_reference_C <- confint(fit_ref_C)

cis_at_m_reference_A["s",]
cis_at_m_reference_B["s",]
cis_at_m_reference_C["s",]

# Any way to just extract all three from fit_ref_A?

解决方法

您可以通过高斯误差传播获得近似 CI：

sum(fixef(fit_ref_A)[c("s","mB:s")]) + 
  c(-1.96,1.96) * sqrt(sum(vcov(fit_ref_A)[c("s","mB:s"),c("s","mB:s")]))
#[1] -0.3346310  0.1863014

或者你可以引导：

myboot <- bootMer(fit_ref_A,function(x) {
  cf <- fixef(x)
  c(sA = cf[["s"]],sB = cf[["s"]] + cf[["mB:s"]],sC = cf[["s"]] + cf[["mC:s"]])
},nsim = 1e4,seed = 42)

apply(myboot$t,2,quantile,probs = c(0.025,0.975))
#              sA         sB         sC
#2.5%  -0.4022927 -0.3415690 -0.3969831
#97.5%  0.2041610  0.1858731  0.1266355