Menu
So, I've fitted a linear mixed model with two random intercepts in R:
Y = X beta + Z b + e_i,
where b ~ MVN (0, Sigma); X and Z are the fixed- and random-effects model matrices respectively, and beta and b are the fixed-effect parameters and random-effects BLUPs/conditional modes.
I would like to get my hands on the underlying covariance matrix of b, which doesn't seem to be a trivial thing in lme4 package. You can get only the variances by VarCorr, not the actual correlation matrix.
According to one of the package vignettes (page 2), you can calculate the covariance of beta: e_i * lambda * t(lambda). And all those components you can extract from the output of lme4.
I was wondering if this is the way to go? Or do you have any other suggestions?
970
The variance-covariance matrix is key: it determines the probability of drawing random effect pairs ⟨S0s,S1s⟩ ⟨ S 0 s , S 1 s ⟩ from the population. You have seen these before, in the chapter on correlation and regression. The covariance matrix is always a square matrix (equal numbers of columns and rows).
To accomplish this in LMER just add the variables for which we want to add random slopes to the random part of the input. This means that (1|class) becomes (1+sex+extrav |class) . We can see that all the fixed regression slopes are still significant.
lme4: Linear Mixed-Effects Models using 'Eigen' and S4 The models and their components are represented using S4 classes and methods. The core computational algorithms are implemented using the 'Eigen' C++ library for numerical linear algebra and 'RcppEigen' "glue".
The random effects: (1 + Time | Chick) which allows individual chicks to vary randomly in terms of their intercept (starting weight) and their effect of Time (weight change over time, also called a “random slope”, but I think that terminology can get confusing when fitting models with nonlinear predictors).
From ?ranef:
If ‘condVar’ is ‘TRUE’ each of the data frames has an attribute called ‘"postVar"’ which is a three-dimensional array with symmetric faces; each face contains the variance-covariance matrix for a particular level of the grouping factor. (The name of this attribute is a historical artifact, and may be changed to ‘condVar’ at some point in the future.)
Set up an example:
library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
rr <- ranef(fm1,condVar=TRUE)
Get the variance-covariance matrix among the b values for the intercept
pv <- attr(rr[[1]],"postVar")
str(pv)
##num [1:2, 1:2, 1:18] 145.71 -21.44 -21.44 5.31 145.71 ...
So this is a 2x2x18 array; each slice is the variance-covariance matrix among the conditional intercept and slope for a particular subject (by definition, the intercepts and slopes for each subject are independent of the intercepts and slopes for all other subjects).
To convert this to a variance-covariance matrix (see getMethod("image",sig="dgTMatrix") ...)
library(Matrix)
vc <- bdiag( ## make a block-diagonal matrix
lapply(
## split 3d array into a list of sub-matrices
split(pv,slice.index(pv,3)),
## ... put them back into 2x2 matrices
matrix,2))
image(vc,sub="",xlab="",ylab="",useRaster=TRUE)
61
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
