[R-sig-ME] Obtaining the variance-covariance matrix for mixed-model averaged parameters

[Ben Bolker]

Paul Mensink <Paul.Mensink at ...> writes:

>    I'm having trouble obtaining a variance-covariance matrix for
> averaged model parameters from a mixed model using the vcov function
> within the MuMIn package (v. 1.9.5).  I've attached an example
> below.  The vcov function works fine with a non-mixed averaged model
> (Example 1); however, when I attempt to use it for a mixed model in
> returns NaNs (Example 2).  In the MuMIn documentation for model.avg
> it says that fitting a model selection object with fit=TRUE should
> store the data required by vcov, so I also attempted to to it that
> way as well (Example 3).  The actual data I am using includes a
> polynomial term (I am aware of the special considerations needed
> when model averaging models that include polynomials), and I need
> the variance-covariance matrix as I am trying to estimate confidenc
> e intervals on the quadratic extremum (as per Hirschberg and Lye
> 2004).

> My questions are:

> (1) Is this a problem with MuMIn, or is it not possible to obtain a
> variance-covariance matrix for mixed model averaged parameters?

  Since the mixed model average parameters are just a weighted
average of the parameters from the individual models this seems
at least plausible to me.  (I guess there's an interesting question
here about how/whether the estimates from the different models
are themselves correlated -- I suspect they are ...)


 This works for me with lme4 1.0-4 (the about-to-be-released
version, available on R-forge and Github)  and MuMin 1.9.5 ...


> #Load MuMIn package
library(MuMIn)
library(lme4)
sessionInfo()
## other attached packages:
## [1] lme4_1.0-4      Matrix_1.0-13   lattice_0.20-23 MuMIn_1.9.5    


> #load Cement data set included in MuMIn package
data(Cement)
> 
#Add in a column, X5, to use as a random factor
set.seed(101)
X5 <- sample(letters[1:4], 13, replace = TRUE) 
Cement1=cbind(Cement, X5)
 
[non-mixed version deleted to make Gmane happy]
 
> ############EXAMPLE # 2 - MIXED MODEL#################################
> 
> # build a global model, specifying X5 as a random intercept term
> # REML set to false for model selection
> 
glob2= lmer(y~X1+X2+(1|X5), REML=FALSE, data=Cement1)
#dredge global model
model.sel2= dredge(glob2)
#obtain model list from model selection object
model.list2 = get.models(model.sel2)
#average model parameters from model list
model.ave2=model.avg(model.list2)
# obtain variance-covariance matrix
vcov(model.ave2)

##             (Intercept)            X1            X2
## (Intercept)  4.02050573 -0.0373575227 -0.0705880379
## X1          -0.03735752  0.0113184076 -0.0009780072
## X2          -0.07058804 -0.0009780072  0.0016174380
## 

#########EXAMPLE # 3 - Mixed model using model.selection  fit=TRUE##########
# build a mixed global model
# REML set to false for model selection
 
glob3= lmer(y~X1+X2+(1|X5), REML=FALSE, data=Cement1)
 
#dredge global model
model.sel3= dredge(glob3)
#directly use model selection object to average model parameters 
model.ave3=model.avg(model.sel3, fit=TRUE)

# obtain variance-covariance matrix
vcov(model.ave3)

## same as above