Questions on the results from glmmPQL(MASS)

Sat, Dec 6, 2008 12:10 PM

zhijie zhang wrote:

This is a philosophical difference between S-PLUS and R.
Since glmmPQL uses quasi-likelihood, technically there is
no log-likelihood (hence no AIC nor BIC, which are based
on the log-likelihood) for this model -- the argument is that
one is limited to looking at Wald tests (testing the Z- or
t-statistics, i.e. parameter estimates divided by estimated
standard errors) for inference in this case.

zhijie zhang wrote:

The (Intercept) random effect is the variance in intercept
across grouping factors .
The residual (0.78) is (I believe) the individual-level error estimated for
the
underlying linear mixed model -- you can probably ignore this.

zhijie zhang wrote:

Q3:In SAS and other softwares, we can get *p*-values for the random effect
'sigma', but i donot see the *p*-values in the results of R/SPLUS. I have
used attributes() to look for them, but no *p* values. Anybody knows how
to
get *p*-values for the random effect 'sigma',.
  Any suggestions or help are greatly appreciated.
#R Results:MASS' version 7.2-44; R version 2.7.2
library(MASS)
summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,family =
binomial,
data = bacteria))

Linear mixed-effects model fit by maximum likelihood
 Data: bacteria
  AIC BIC logLik
   NA  NA     NA

Random effects:
 Formula: ~1 | ID
        (Intercept)  Residual
StdDev:    1.410637 0.7800511

Variance function:
 Structure: fixed weights
 Formula: ~invwt
Fixed effects: y ~ trt + I(week > 2)
                    Value          Std.Error      DF       t-value
p-value
(Intercept)      3.412014    0.5185033   169       6.580506  0.0000
trtdrug         -1.247355     0.6440635     47      -1.936696  0.0588
trtdrug+        -0.754327    0.6453978     47      -1.168779  0.2484
I(week > 2)TRUE -1.607257 0.3583379 169     -4.485311  0.0000
 Correlation:
                (Intr) trtdrg trtdr+
trtdrug         -0.598
trtdrug+        -0.571  0.460
I(week > 2)TRUE -0.537  0.047 -0.001

#S-PLUS8.0: The results are the same as R except the followings:
AIC      BIC    logLik
1113.622 1133.984 -550.8111

#SAS9.1.3
proc nlmixed data=b;
 parms b0=-1 b1=1 b2=1 sigma=0.4;
 yy=b0+u+b1*trt+b2*week;
 p=1/(1+exp(-yy));
 Model response~binary(p);
 Random u~normal(0,sigma) subject=id;
Run;

 -2 Log Likelihood = 192.2
 AIC (smaller is better)=200.2
AICC (smaller is better) =200.3
BIC (smaller is better)= 207.8


                                      Parameter Estimates

                       Standard
  Parameter  Estimate     Error    DF  t Value  Pr > |t|   Alpha
Lower     Upper  Gradient

  b0             3.4966    0.6512    49     5.37    <.0001    0.05
2.1880    4.8052  -4.69E-6
  trt              -0.6763    0.3352    49    -2.02    0.0491    0.05
-1.3500  -0.00266  -0.00001
  I(week>2)   -1.6132    0.4785    49    -3.37    0.0015    0.05   -2.5747
-0.6516  -9.35E-7
  sigma        1.5301    0.9632    49     1.59    0.1186    0.05
-0.4054    3.4656  -2.42E-6

In general (alas) it is *extremely* difficult to get *correct*
p-values for effects (both fixed and random, although fixed
might be even worse) in GLMMs, despite the fact that SAS
will happily give them to you.  In general you can get
p-values for random effects via a likelihood ratio test on
the difference of nested models with and without 
the relevant effects.  In this case that's a little bit trickier
because (1) glmmPQL won't give you log-likelihoods 
(2) glmmPQL won't fit models without any random effects
at all, and comparing log-likelihoods across different
fitting procedures is always a little risky -- you have to make
sure they are including the same constants in the log-likelihood.
A partial solution is to use the lmer function in the lme4 package:

 lme2 <- lmer(y ~ trt + I(week > 2)+(1|ID), family = binomial, data =
bacteria)

which gives you similar estimates (a good idea in any case, because
in any case PQL is not reliable for binary data -- see Breslow 2003).
This will give you a likelihood etc. for the model.  You still need to
work out whether comparing AIC/BIC log-likelihood between lmer
and glm (which you need to fit the model without random effects)
is sensible.

  I would strongly recommend that you follow up further questions
on this topic to r-sig-mixed-models at r-project.org , which is a 
special mailing list for mixed models.

  good luck,
   Ben Bolker

View this message in context: http://www.nabble.com/Questions-on-the-results-from-glmmPQL%28MASS%29-tp20867935p20874037.html
Sent from the R help mailing list archive at Nabble.com.

Questions on the results from glmmPQL(MASS)

Thread (3 messages)