On Nov 13, 2007 9:55 AM, Jimenez-Leal William
<william.jimenezleal@lsc.gov.uk> wrote:
> Dear S plus users,
> I am fitting a bunch of mixed effects logistic regressions using glmmPQL and
> I have a little problem extracting their random coefficients. I am modeling
> the probability of passing an exam against student's prior attainment,
> allowing the latter to vary randomly across schools. Looks something like
> this
> a<-glmmPQL(POINTS~PA +PA2+PA3,random=~PA|INST, family=binomial(link=logit),
> data=x, niter=1000)
> ('pa' is the prior attainment, 'inst' is school)
> Normally, I just produce a summary of this object, which gives you the
> corresponding standard deviations and correlations, and work out the
> variance-covariance matrix from there. This is the summary:
> summary(a)
> Linear mixed-effects model fit by maximum likelihood
> Data: k
> AIC BIC logLik
> 18142.24 18189.96 -9063.118
> Random effects:
> Formula: ~ PA | INST
> Structure: General positive-definite
> StdDev Corr
> (Intercept) 3.68174225 (Intr)
> PA 0.09902275 0.253
> Residual 0.60168538
> Variance function:
> Structure: fixed weights
> Formula: ~ invwt
> (Etc, etc)
> Number of Observations: 2881
> Number of Groups: 142
> The problem is that when I extract this matrix from the summary object
> directly, the random effects numbers do not match.
> summary(a)$modelStruct$reStruct
> Random effects:
> Formula: ~ PA | INST
> Structure: General positive-definite
> StdDev Corr
> (Intercept) 6.1190489 (Intr)
> PRIORC 0.1645756 0.253
> Residual 1.0000000
> Would anyone tell me what I am missing here? Clearly, the matrix is being
> transformed when summarized, but why and how?
The actual optimization in the lme function is done on the scale of
the relative variance-covariance of the random effects. That is, the
variance-covariance of the random effects is expressed as a multiple
of the "sigma^2" in the variance of the observations given the random
effects. If you extract only the reStruct it will show the residual
variance and standard deviation as unity.
For this model you may question which of the variance-covariance
estimates you want to use. In the case of a binomial distribution for
the responses there is no separate variance parameter. The
conditional mean determines the conditional variance. The estimate of
a "Residual" variance occurs only because of the PQL approximation.
According to the model the second variance covariance matrix, where
the residual variance and standard deviation has been set to unity,
would be the correct one to use.
|