首页 文章

来自glmer输出的优势比和置信区间

提问于
浏览
10

我制作了一个模型,该模型着眼于许多变量及其对妊娠结局的影响 . 结果是分组二进制 . 一群动物将有34个怀孕和3个空,接下来将有20个怀孕和4个空等等 .

我使用 glmer 函数对这些数据进行了建模,其中y是妊娠结果(怀孕或空白) .

mclus5 <- glmer(y~adg + breed + bw_start + year + (1|farm),
                data=dat, family=binomial)

我得到所有通常的系数等输出但是对于解释我想将其转换为每个系数的优势比和置信区间 .

在过去的逻辑回归模型中,我使用了以下代码

round(exp(cbind(OR=coef(mclus5),confint(mclus5))),3)

这将很好地提供我想要的东西,但它似乎不适用于我运行的模型 .

有谁知道我可以通过R为我的模型获得此输出的方式?

3 回答

  • 4

    唯一真正的区别是你必须使用 fixef() 而不是 coef() 来提取固定效应系数( coef() 给出每组的估计系数) .

    我将使用 lme4 包中的内置示例进行说明 .

    library("lme4")
    gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                   data = cbpp, family = binomial)
    

    固定效应系数和置信区间,对数 - 赔率表:

    cc <- confint(gm1,parm="beta_")  ## slow (~ 11 seconds)
    ctab <- cbind(est=fixef(gm1),cc)
    

    (如果你想要更快但更准确的Wald置信区间,你可以使用 confint(gm1,parm="beta_",method="Wald") ;这相当于@Gorka的答案,但稍微方便一些 . )

    Exponentiate得到比值比:

    rtab <- exp(ctab)
    print(rtab,digits=3)
    ##               est 2.5 % 97.5 %
    ## (Intercept) 0.247 0.149  0.388
    ## period2     0.371 0.199  0.665
    ## period3     0.324 0.165  0.600
    ## period4     0.206 0.082  0.449
    

    一个稍微简单/更通用的解决方案:

    library(broom.mixed)
    tidy(gm1,conf.int=TRUE,exponentiate=TRUE,effects="fixed")
    

    对于Wald间隔,或为配置文件置信区间添加 conf.method="profile" .

  • 3

    我相信还有另一种更快的方法(如果你的结果不太准确) .

    来自:http://www.ats.ucla.edu/stat/r/dae/melogit.htm

    首先,我们得到估计的置信区间

    se <- sqrt(diag(vcov(mclus5)))
    # table of estimates with 95% CI
    tab <- cbind(Est = fixef(mclus5), LL = fixef(mclus5) - 1.96 * se, UL = fixef(mclus5) + 1.96 * se)
    

    然后比值比为95%CI

    print(exp(tab), digits=3)
    
  • 14

    我认为其他选项只是使用包 emmeans

    library(emmeans)
    data.frame(confint(pairs(emmeans(fit, ~ factor_name,type="response"))))
    

相关问题