[R] glmmPQL and predict

[Prof Brian Ripley]

The whole of idea of 'level' in mixed models is confusing to some. 
Professor Snijders (who teaches our students) and Professor Bates label 
from opposite ends.

But, assuming this is my work in package MASS (Master Harwood: it is 
childish, to put it mildly, to fail to give due credit), it follows lme 
in package nlme: glmmPQL's predict method just passes this on to nlme, 
and the documentation is identical.

So not only was it unfair to fail to mention which package and whose 
work this was, it was even more unfair to attribute personal lack of 
understanding to my work and not package nlme.

Just because R and many contributed packages are free does not entitle 
you to treat them as zero-value: very much to the contrary.

On 10/01/2012 16:38, Ben Bolker wrote:
> 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
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595