3

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

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)orlmerTest::lsmeansLT(gl, 'am')并绘制两组的最小二乘均值及其置信区间。

我怎样才能对普通模型做同样的事情？

4

1 回答 1

3

（这是经过一些讨论和进一步发现后的修订答案。请注意，我是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)产生了一个非正定协方差矩阵，导致 EMM 的方差估计为负。

我警告说，虽然这似乎在计算上有效，但这并不意味着答案是正确的！

于 2017-11-25T04:18:07.883 回答