3
$\begingroup$

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?
$\endgroup$

1 Answer 1

2
$\begingroup$

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)"
$\endgroup$
2
  • $\begingroup$ but in order to use VarCorr(pd,sigma=sigma(fm2)), I need to use pd <- pdLogChol(coef(fm2$modelStruct$reStruct$Subject, unconstrained = T)). I'm I right? $\endgroup$
    – utobi
    Nov 26, 2018 at 22:30
  • 1
    $\begingroup$ that seems right ... $\endgroup$
    – Ben Bolker
    Nov 26, 2018 at 22:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.