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 predictionX01 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. Exponentiating the intercept coefficients from the fitted regression, the level=0 values match when the random effect intercept is included> 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 Thanks! Mike

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 includedDo 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.9828449This 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