glmmPQL: std.errors in summary and vcov differ

A reproducible example would be nice.  On the other hand, we can
reproduce this with the example in ?glmmPQL:


library(MASS)
g1 <- glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,
                      family = binomial, data = bacteria)
s1 <- summary(g1)$tTab[,"Std.Error"]
s2 <- sqrt(diag(vcov(g1)))
all.equal(s1,s2)
## [1] "Mean relative difference: 0.009132611"

Let's dig in to see what's happening. As a starting point, class(g1)
is c("glmmPQL","lme"), and if we look at methods("vcov"),
methods("summary"), we can see there are no glmmPQL methods, so these
computations are being handled by nlme:::summary.lme and
nlme:::vcov.lme

vcov.lme is

object$varFix

summary.lme has

stdFixed <- sqrt(diag(as.matrix(object$varFix)))
  if (adjustSigma && object$method == "ML")
        stdFixed <- stdFixed * sqrt(object$dims$N/(object$dims$N -
            length(stdFixed)))

So what's happening is that vcov() is giving you the uncorrected
(maximum likelihood) estimate of the variance-covariance matrix, while
summary is giving you the bias-corrected version.  Since your sample
size is large, the difference is tiny.

Dear list,

With glmmPQL from package MASS,  I ran a logistic regression model with
an intercept only, which is random across 33 countries:

m1 <- glmmPQL(MATH_Top1 ~ 1,
               random = list(country = ~ 1),
               family=binomial, data=pisas)
summary(m1)
sqrt(vcov(m1))

There is a total of N=22997 students in data.frame pisas.

The summary(m1) functions shows:

Fixed effects: MATH_Top1 ~ 1
                 Value Std.Error    DF  t-value p-value
(Intercept) -2.176564 0.1140811 22964 -19.0791       0


whereas the sqrt(vcov(m1)) function shows:

             (Intercept)
(Intercept)   0.1140786


My question is why the two std. error estimates 0.1140811 and 0.1140786
of the fixed part of the intercept differ slightly. The same holds for
std. errors of the fixed effects of predictor variables, which I did not
add above for reasons of simplicity. Thanks for any help!

Ben.



        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models