Solved – Variance-Covariance Structure in lme/nlme of random effects

Question

I'm simulating data (only one level of grouping) and then I fit a (non-linear) mixed model.

mod1 <- nlme(conc~SSmicmen(time, Vm, K),data=groupedData(conc~time|subject,data=sample1),fixed=Vm+K~1,start=c(Vm=5,K=10))

I want to extract the variance-covariance structure of the random effects. Therefore I use VarCorr and I get:

    > VarCorr(mod1)
subject = pdLogChol(list(Vm ~ 1,K ~ 1)) 
         Variance   StdDev    Corr  
Vm       0.10380769 0.3221920 Vm    
K        0.04527981 0.2127905 -0.977
Residual 0.54599683 0.7389160  

So, after very long research I still don't know whether this output now gives me the covariance matrix of the random effects, or the precision factor.

See the book of Pinheiro and Bates: in their model assumption (page 311) they assume the random effects to be normally distributed with expectation zero and covariance matrix $\psi$.

Furthermore, there is the precision factor such that
$$
\psi^{-1}=1/\sigma^2\cdot\Delta^T\Delta.
$$
So I am really confused now, what VarCorr is exactly providing, is it $\Delta$, $\psi$ or something else?
Unfortunately none of these possibilities yields a variance-covariance matrix $\psi$ which is similar (not at all!!) to the one I used to simulate the data, but I don't know whether this is a normal thing when dealing with mixed effects models.

Answer 1

I know that lme reports the actual variance-covariance matrix (not the precision factor or the scaled variance-covariance matrix. So, to double-check the results for nlme, I set up the same model in lme and nlme:

library(nlme)
m0 <- lme(distance~age,random=~age|Subject,Orthodont,method="ML")
VarCorr(m0)
##             Variance   StdDev    Corr  
## (Intercept) 4.81407327 2.1940996 (Intr)
## age         0.04619252 0.2149244 -0.581
## Residual    1.71620466 1.3100399       

Compare with nlme results:

m1 <- nlme(distance~a+b*age,
           fixed=a+b~1,
           random=a+b~1|Subject,
           Orthodont,
           start=c(a=1,b=1),method="ML")
##          Variance   StdDev    Corr  
## a        4.81408829 2.1941031 a     
## b        0.04619255 0.2149245 -0.581
## Residual 1.71620375 1.3100396       

I can imagine either that you're mistaken about the parameterization of your simulation (it's very easy to do, e.g. working on the std dev rather than the variance scale), or that the sampling distribution of your variances is very wide (so you can get an estimate far from the true value); in the latter case, try fitting a few different realizations of the simulation to see how much the results vary ...