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

[Paul Mensink]

Hi there,

   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 confidence 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?
(2) If it is a problem with MuMIn, is there any other way to obtain the variance-covariance matrix for mixed model averaged parameters? 

Thanks in advance,
Paul


#Load MuMIn package
library(MuMIn)

#load Cement data set included in MuMIn package
data(Cement)

#Add in a column, X5, to use as a random factor
X5 <- sample(letters[1:4], 13, replace = TRUE) 
Cement1=cbind(Cement, X5)


#############EXAMPLE 1 - NON-MIXED##############################################

# build a non-mixed global model
# REML set to false for model selection

glob1= lm(y~X1+X2, data=Cement1)

#dredge global model
model.sel1= dredge(glob1)

#obtain model list from model selection object
model.list1 = get.models(model.sel1)

#average model parameters from model list
model.ave1=model.avg(model.list1)

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

            (Intercept)           X1           X2
(Intercept)  5.22662252 -0.048565024 -0.091764378
X1          -0.04856502  0.014713921 -0.001271409
X2          -0.09176438 -0.001271409  0.002102662



############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)         NaN     NaN   NaN
X1                     NaN    NaN    NaN
X2                     NaN    NaN    NaN


#########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)

                             (Intercept)  X1    X2
(Intercept)         NaN     NaN   NaN
X1                     NaN    NaN    NaN
X2                     NaN    NaN    NaN