[R-sig-ME] glmmPQL: std.errors in summary and vcov differ
Ben Bolker
bbolker at gmail.com
Wed May 31 13:40:56 CEST 2017
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
More information about the R-sig-mixed-models
mailing list