R-平均模型的绘图模型适合离散变量

Question

我有一组线性混合模型，并创建了一个平均模型。我想绘制模型拟合的两个级别，其中包括平均模型中的两个因素。一个简单的例子：

library(lme4)
library(MuMIn)
mtcars2 <- mtcars
mtcars2$vs <- factor(mtcars2$vs)
gl <- lmer(mpg ~ am + disp + hp + qsec + (1 | cyl), mtcars2, 
           REML = FALSE, na.action = 'na.fail')
d <- dredge(gl)
av <- model.avg(d, subset = cumsum(weight) <= 0.95)
summary(av)

Call:
model.avg(object = d, subset = cumsum(weight) <= 0.95)
Component model call: 
lme4::lmer(formula = mpg ~ <7 unique rhs>, data = mtcars2, REML = FALSE, na.action = na.fail)
Component models: 
     df logLik   AICc delta weight
13    5 -77.81 167.92  0.00   0.37
123   6 -76.34 168.05  0.13   0.35
134   6 -77.54 170.43  2.51   0.11
1234  7 -76.25 171.16  3.24   0.07
23    5 -79.85 172.00  4.08   0.05
2     4 -81.63 172.75  4.83   0.03
124   6 -78.99 173.34  5.42   0.02
Term codes: 
  am disp   hp qsec 
   1    2    3    4 
Model-averaged coefficients:  
(full average) 
             Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept) 25.457505   6.467643    6.648016   3.829 0.000129 ***
am           4.103425   1.861593    1.898182   2.162 0.030636 *  
hp          -0.043829   0.017926    0.018265   2.400 0.016415 *  
disp        -0.009419   0.011834    0.011983   0.786 0.431821    
qsec         0.081973   0.284147    0.292015   0.281 0.778929    
(conditional average) 
            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept) 25.45751    6.46764     6.64802   3.829 0.000129 ***
am           4.46519    1.46823     1.51835   2.941 0.003273 ** 
hp          -0.04651    0.01471     0.01515   3.070 0.002140 ** 
disp        -0.01793    0.01068     0.01099   1.632 0.102634    
qsec         0.40421    0.51757     0.53873   0.750 0.453075    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Relative variable importance: 
                     hp   am   disp qsec
Importance:          0.94 0.92 0.53 0.20
N containing models:    5    5    5    3

我想绘制由完整平均模型估计的am的效果。

通常，我会使用lsmeans::lsmeans(gl, ~am)或lmerTest::lsmeansLT(gl, 'am')，并为两组绘制最小二乘手段及其置信区间。

我该如何为平均模型做同样的事情？

Answer 1

（这是一个修订的答案，经过一些讨论和进一步的发现。请注意，我是emmeans软件包作者。）

这似乎是有效的。

首先，定义 emmeans 的方法：

library(emmeans)
terms.averaging = function(x, ...)
    terms(x$formula)
recover_data.averaging = emmeans:::recover_data.lm
    ### NOTE: still have to provide 'data' argument
emm_basis.averaging = function(object, trms, xlev, grid, ...) {
    bhat = coef(object, full = TRUE)
    m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)
    X = model.matrix(trms, m, contrasts.arg = object$contrasts)
    V = vcov(object, full = TRUE)
    dffun = function(k, dfargs) NA
    dfargs = list()
    list(X=X, bhat=bhat, nbasis=estimability::all.estble, V=V, 
         dffun=dffun, dfargs=dfargs, misc=list())
}

需要terms方法，因为没有一种方法。另外两个是根据lm对象的现有方法改编的。现在有一个捕获：vcov()调用要求对象具有非NULL "modelList"属性。并且您的av对象失败。但是model.avg帮助页面底部的示例显示了要做什么：

cs95 = get.models(d, cumsum(weight) <= .95)
AV = model.avg(cs95)

现在，AV具有所需的属性。现在我们得到：

em = emmeans(AV, ~ am, at = list(am = c("0", "1")), data = mtcars)
em
## am   emmean       SE df asymp.LCL asymp.UCL
##  0 15.42665 2.985460 NA  9.575257  21.27805
##  1 19.53008 1.986149 NA 15.637297  23.42286
pairs(em)
## contrast  estimate       SE df z.ratio p.value
## 0 - 1    -4.103425 1.861593 NA  -2.204  0.0275

请注意，对比结果与模型摘要表中的av匹配估计值和未调整的SE。

注意：使用coef(..., full = FALSE)和vcov(... full = FALSE)产生了一个非阳性非限性协方差矩阵，从而对EMMS进行负方差估计。

我警告说，尽管这似乎在计算上起作用，但这并不意味着答案是正确的！