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

[Steve Candy]

Hi Ben

 

I just noticed a similar problem with lme( ) when using the method="ML"
option. The summary(lmeObject) parameter SEs and  diag(lmeObject
$varFix)^0.5 do not exactly match whereas they do if the "REML" option is
used. See below. This doesn't entire explain what's going on, since the
default option in lme() is method="REML", so the default option in the call
of lme() from glmmPQL should also be method="REML". Unless somehow its been
set as "ML" by glmmPQL. It might be worth using the control option or in
"further arguments for lme" to specify the method in lme() explicitly
through glmmPQL . Either way it would be good to find out why there is this
mismatch for method="ML".

 

> lme.01r <- lme(fixed=varDV  ~ Condition_f+Time_f+Condition_f:Time_f,
method="REML", random=~1 | PRE_Class/ID_f, 

+                    data=data_CrepA)

> 

> summary(lme.01r)

Linear mixed-effects model fit by REML

Data: data_CrepA 

       AIC      BIC    logLik

  2949.689 3023.109 -1459.845

 

Random effects:

Formula: ~1 | PRE_Class

        (Intercept)

StdDev:   0.2882176

 

Formula: ~1 | ID_f %in% PRE_Class

        (Intercept) Residual

StdDev:   0.7799864 0.848057

 

Fixed effects: varDV ~ Condition_f + Time_f + Condition_f:Time_f 

                            Value Std.Error  DF   t-value p-value

(Intercept)             1.6042292 0.2287389 715  7.013363  0.0000

 

...snip

 

> (summary(lme.01r)$tTable)[1,]

       Value    Std.Error           DF      t-value      p-value 

1.604229e+00 2.287389e-01 7.150000e+02 7.013363e+00 5.400002e-12 

> lme.01r$coefficients$fixed[1]

(Intercept) 

   1.604229 

            0.01312038             0.02010276             0.01336595
0.02042342             0.02986909             0.04595728

> lme.01r$varFix[1,1]^0.5

[1] 0.2287389

> sqrt(vcov(lme.01r)[1,1])

[1] 0.2287389

> ### SE the same as the summary output

 

> Now use ML

> 

> lme.01r <- lme(fixed=varDV  ~ Condition_f+Time_f+Condition_f:Time_f,
method="ML", random=~1 | PRE_Class/ID_f, 

+                    data=data_CrepA)

> summary(lme.01r)

Linear mixed-effects model fit by maximum likelihood

Data: data_CrepA 

       AIC      BIC    logLik

  2922.624 2996.225 -1446.312

 

Random effects:

Formula: ~1 | PRE_Class

        (Intercept)

StdDev:   0.2415063

 

Formula: ~1 | ID_f %in% PRE_Class

        (Intercept)  Residual

StdDev:   0.7810628 0.8428382

 

Fixed effects: varDV ~ Condition_f + Time_f + Condition_f:Time_f 

                            Value Std.Error  DF   t-value p-value

(Intercept)             1.5993361 0.2107460 715  7.588926  0.0000

 

...snip

 

> (summary(lme.01r)$tTable)[1,]

       Value    Std.Error           DF      t-value      p-value 

1.599336e+00 2.107460e-01 7.150000e+02 7.588926e+00 1.008003e-13 

> lme.01r$coefficients$fixed[1]

(Intercept) 

   1.599336 

> lme.01r$varFix[1,1]^0.5

[1] 0.2094765

> sqrt(vcov(lme.01r)[1,1])

[1] 0.2094765

> ### SE NOT the same as the summary output

 

>Dear list,

 

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

 

.snip

 

>>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.

 

Dr Steven G. Candy

Director/Consultant

SCANDY STATISTICAL MODELLING PTY LTD

(ABN: 83 601 268 419)

70 Burwood Drive

Blackmans Bay, TASMANIA, Australia 7052

Mobile: (61) 0439284983

 


	[[alternative HTML version deleted]]