[R-sig-ME] glmmPQL: std.errors in summary and vcov differ

[Ben Bolker]

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.


On Wed, May 31, 2017 at 4:37 AM, Ben Pelzer <b.pelzer at maw.ru.nl> wrote:
> 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