图模型适合离散变量，从平均模型

I have a set of linear mixed models, and have created an average model. I'd like to plot the model fits for two levels of a factor, included in the average model. A simple example:

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

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

I want to plot the effect of am as estimated by the full averaged model.

我想用全平均模型估计am的影响。

Normally I would use lsmeans::lsmeans(gl, ~am) or lmerTest::lsmeansLT(gl, 'am') and plot the least squares means for the two groups and their confidence intervals.

通常我会使用lsmeans::lsmeans(gl， ~am)或lmerTest::lsmeansLT(gl, am)，并对两个组和它们的置信区间绘制最小二乘法。

How can I do the same for the average model?

如何对平均模型做同样的处理呢?

1 个解决方案

#1

(This is a revised answer, after some discussion and further findings. Note that I'm the emmeans package author.)

(经过讨论和进一步的发现，这是一个修正的答案。请注意，我是emmeans包的作者。

Here is something that appears to work.

这似乎是可行的。

First, define methods needed by the emmeans package:

首先，定义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())
}

The terms method is needed because there isn't one. The other two are adapted from the existing methods for lm objects. Now there is one catch: the vcov() call requires the object to have a non-NULL "modelList" attribute. And your av object fails. But the examples at the bottom of the help page for model.avg shows what to do:

因为没有terms方法，所以需要terms method。另外两个是根据lm对象的现有方法改编的。现在有一个问题:vcov()调用要求对象具有非空的“modelList”属性。你的av对象失败了。但是在模型的帮助页面底部的例子。avg显示了要做什么:

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

Now, AV has the required attribute. Now we get:

现在，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

Note that the contrast result matches the estimate and unadjusted SE for av in the model summary table.

注意，对比结果与模型汇总表中对av的估计和未调整的SE相符。

Note: Using coef(..., full = FALSE) and vcov(... full = FALSE) yielded a non-positive-definite covariance matrix, resulting in negative variance estimates for the EMMs.

注意:使用系数(…， full = FALSE, vcov(…full = FALSE)产生了一个非正定协方差矩阵，导致EMMs的方差估计为负。

And I caution that while this seems to work computationally, that does not imply that the answers are right!

我要提醒大家的是，尽管这似乎在计算上行得通，但这并不意味着答案是正确的!

#1