如何从 logit 模型（glmer）中获取两组成功概率差异的概况置信区间？

Question

我正在努力将从 logit 模型获得的对数优势比配置文件置信区间转换为概率。我想知道如何计算两组之间差异的置信区间。

如果 p 值 > 0.05，则 95% CI 的差异应介于零以下和零以上之间。但是，我不知道当必须对对数比率取幂时如何获得负值。因此，我尝试计算其中一组 (B) 的 CI，看看 CI 的下限和上限与 A 组的估计值有何不同。我认为这不是计算差值 CI 的正确方法，因为 A 的估计值也是不确定的。

如果有人能帮助我，我会很高兴。

library(lme4)    
# Example data: 
set.seed(11)
treatment = c(rep("A",30), rep("B", 40))
site = rep(1:14, each = 5)
presence = c(rbinom(30, 1, 0.6),rbinom(40, 1, 0.8))
df = data.frame(presence, treatment, site)

# Likelihood ratio test 
M0 = glmer(presence ~ 1 + (1|site), family = "binomial", data = df)
M1 = glmer(presence ~ treatment + (1|site), family = "binomial", data = df)
anova(M1, M0)

# Calculating confidence intervals
cc <- confint(M1, parm = "beta_")
ctab <- cbind(est = fixef(M1), cc)
cdat = as.data.frame(ctab)

# Function to back-transform to probability (0-1)
unlogit = function(y){
    y_retransfromed = exp(y)/(1+exp(y))
    y_retransfromed
}

# Getting estimates
A_est = unlogit(cdat$est[1]) 
B_est = unlogit(cdat$est[1] + cdat$est[2])
B_lwr = unlogit(cdat$est[1] + cdat[2,2])
B_upr = unlogit(cdat$est[1] + cdat[2,3])

Difference_est = B_est - A_est

# This is how I tried to calculate the CI of the difference
Difference_lwr = B_lwr - A_est
Difference_upr = B_upr - A_est

# However, I believe this is wrong because A_est is also “uncertain”

如何得到存在概率之差的置信区间？

Answer 1

我们可以通过以下方式计算平均治疗效果。根据原始数据，创建两个新数据集，其中一个所有单元都接受治疗 A，另一个所有单元都接受治疗 B。现在，根据您的模型估计（在您的情况下，M1），我们计算这两个数据集中每个单元的预测结果。然后我们计算两个数据集之间结果的平均差异，以获得我们估计的平均治疗效果。在这里，我们可以编写一个函数，它接受一个 glmer 对象并计算平均治疗效果：

ate <- function(.) {
  treat_A <- treat_B <- df
  treat_A$treatment <- "A"
  treat_B$treatment <- "B"
  c("ate" = mean(predict(., newdata = treat_B, type = "response") -
    predict(., newdata = treat_A, type = "response")))
}
ate(M1)
#        ate 
# 0.09478276

我们如何得到不确定区间？我们可以使用 bootstrap，即 re-estimate 模型多次使用从您的原始数据中随机生成的样本，每次计算平均处理效果。然后我们可以使用 bootstrapped 平均治疗效果的分布来计算我们的不确定性区间。这里我们使用 bootMer 函数

生成 100 个模拟

out <- bootMer(M1, ate, seed = 1234, nsim = 100)

并检查效果的分布：

quantile(out$t, c(0.025, 0.5, 0.975))
#        2.5%         50%       97.5% 
# -0.06761338  0.10508751  0.26907504

如何从 logit 模型（glmer）中获取两组成功概率差异的概况置信区间？

How to obtain profile confidence intervals of the difference in probability of success between two groups from a logit model (glmer)?

r

lme4

confidence-interval

logistic-regression