[R] glmmPQL and predict

Ben Bolker bbolker at gmail.com
Tue Jan 10 17:38:00 CET 2012


Mike Harwood <harwood262 <at> gmail.com> writes:

> Is the labeling/naming of levels in the documentation for the
> predict.glmmPQL function "backwards"?  The documentation states "Level
> values increase from outermost to innermost grouping, with level zero
> corresponding to the population predictions".  Taking the sample in
> the documentation:
> 
> fit <- glmmPQL(y ~ trt + I(week > 2), random = ~1 |  ID,
>                family = binomial, data = bacteria)
> 
> > head(predict(fit, bacteria, level = 0, type="response"))
> [1] 0.9680779 0.9680779 0.8587270 0.8587270 0.9344832 0.9344832
> > head(predict(fit, bacteria, level = 1, type="response"))
>       X01       X01       X01       X01       X02       X02
> 0.9828449 0.9828449 0.9198935 0.9198935 0.9050782 0.9050782
> > head(predict(fit, bacteria, type="response")) ## population prediction
>       X01       X01       X01       X01       X02       X02
> 0.9828449 0.9828449 0.9198935 0.9198935 0.9050782 0.9050782
> 
> The returned values for level=1 and level=<unspecified> match, which
> is not what I expected based upon the documentation.

  Well, the documentation says: "Defaults to the highest or innermost level of
grouping", which is level 1 in this case -- right?

> Exponentiating
> the intercept coefficients from the fitted regression, the level=0
> values match when the random effect intercept is included

  Do you mean "is NOT included" here?

  0.9680779 (no random effect, below) matches the level=0 prediction above
  0.9828449 (include random effect, below) matches the level=1 prediction,
which is also the default prediction, above.

> 
> > 1/(1+exp(-3.412014)) ## only the fixed effect
> [1] 0.9680779
> > 1/(1+exp(-1*(3.412014+0.63614382))) ## fixed and random effect intercepts
> [1] 0.9828449

  This all matches my expectations.  If your expectations still go
in the other direction, could you explain in more detail?

  By the way, I recommend r-sig-mixed-models at r-project.org for
mixed model questions in general ...

  Ben Bolker



More information about the R-help mailing list