why pdLogChol initialised at a model fit does not match the pdLogChol of the fitted model?

Question

I'm trying to build and initialize a pdLogChol object from a fitted lme model. However, although the pdLogChol is started from the model fit, the VarCorr function gives different values when compared with that of the model fit. I expected to see identical values but that's not the case. Is this a normal behaviour?

library(nlme)
fm2 <- lme(distance ~ 1, data = Orthodont, random = ~ 1|Subject)
coef(fm2$modelStruct$reStruct$Subject, unconstrained = F)
pd <- pdLogChol(coef(fm2$modelStruct$reStruct$Subject, unconstrained = F))

VarCorr(fm2) # this must be equal to
VarCorr(pd) # but it is not. Why?

Answer 1

It appears that the random-effects variance is parameterized relative to the residual variance:

(cc <- coef(fm2$modelStruct$reStruct$Subject, unconstrained=FALSE)) 
## 0.7610834 
(res_var <- sigma(fm2)^2)
## [1] 4.929783
cc*res_var
## var((Intercept)) 
## 3.751976 

Compare with fitted model:

VarCorr(fm2)
## Subject = pdLogChol(1) 
##             Variance StdDev  
## (Intercept) 3.751976 1.937002
## Residual    4.929783 2.220312

Looking at the VarCorr.lme method, it seems that the way to do this scaling more idiomatically is to pass an appropriate sigma argument to VarCorr():

pd <- pdLogChol(coef(fm2$modelStruct$reStruct$Subject, 
   unconstrained = TRUE))
VarCorr(pd,sigma=sigma(fm2))
##    Variance   StdDev
## V1 3.751976 1.937002
## attr(,"formStr")
## [1] "pdLogChol(ULL)"