Plot model fit for discrete variable, from average model

Question

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.

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.

How can I do the same for the average model?

Answer 1

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

Here is something that appears to work.

First, define methods needed by the emmeans package:

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:

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

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

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.

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

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