Prediction of random effects in glmer()

Mon, Feb 15, 2021 1:58 PM

Thanks Ben.

To be more specific, I want to infer from your comment that the intervals are not guaranteed to cover the true unknown value 95% of the time under an infinite sequence of experiments. Is there anything else?

Thanks very much for the Cox paper - I hadn't seen that.

Cheers,

Andrew

--
Andrew Robinson
Director, CEBRA and Professor of Biosecurity,
School/s of BioSciences and Mathematics & Statistics
University of Melbourne, VIC 3010 Australia
Tel: (+61) 0403 138 955 Email: apro at unimelb.edu.au
Website: https://researchers.ms.unimelb.edu.au/~apro at unimelb/
I acknowledge the Traditional Owners of the land I inhabit, and pay my respects to their Elders.

On Feb 16, 2021, 8:40 AM +1100, Ben Bolker <bbolker at gmail.com>, wrote:

Good question.
Going back to a previous answer in this thread; these are empirical
Bayesian prediction intervals. That is, they may be well-calibrated but
they don't have the same theoretical status as 'true' frequentist CIs.

To be completely honest, I don't know what is or isn't known about
the calibration properties of intervals based on the conditional SD (I
would guess they're not bad but really don't know the theory).

This looks like a good reference on EB confidence intervals:

Cox, D. R. ?Prediction Intervals and Empirical Bayes Confidence
Intervals.? Journal of Applied Probability 12, no. S1 (ed 1975): 47?55.
https://doi.org/10.1017/S0021900200047550.

Based on a *very* quick read it seems as though constructing CIs for
the conditional modes based on the SDs is OK subject to the assumption
that we can plug in/condition on the 'theta'/variance estimates ...

On 2/15/21 4:25 PM, Andrew Robinson wrote:

Hi Ben,

I note with interest below:

Note that these SDs **are not appropriate for testing hypotheses

about individual levels of the random effect**: one of the tradeoffs of

using REs in the first place is that you trade parsimony for the ability

to do hypothesis testing on individual predictions.


Does this mean that they are also inappropriate for computing interval
predictions for individual predictions?

Where might one read more about this observation?

Many thanks,

Andrew

--
Andrew Robinson
Director, CEBRA and Professor of Biosecurity,
School/s of BioSciences and Mathematics & Statistics
University of Melbourne, VIC 3010 Australia
Tel: (+61) 0403 138 955 Email: apro at unimelb.edu.au
Website: https://researchers.ms.unimelb.edu.au/~apro at unimelb/
I acknowledge the Traditional Owners of the land I inhabit, and pay my
respects to their Elders.

On Feb 16, 2021, 8:16 AM +1100, Ben Bolker <bbolker at gmail.com>, wrote:

The penalty term ||u||^2 appears first in the likelihood expression
(eq 11 in the paper),

int exp( - (sum(d_i(y_obs,u)) + ||u||^2)/2) * C du

where C is a normalization constant I don't feel like writing out.

The sum(d_i(.)) part is the deviance of the GLM part of the model.
The other part would come out to exp(-||u||^2/2) if we factored it out;
this expresses the idea that the u values (the 'spherized' random
effects, i.e. on a scale where they are iid N(0,1) variables) are
assumed to be drawn from a Normal distribution.

This stuff is probably better explained in the original JSS lmer
paper (where the same penalty term appears, for the same reason).

An *analog* of standard errors is indeed available for the
conditional modes (this is also explained in the JSS paper); we call
these "conditional variances" and "conditional standard deviations" - we
can't quite call them "standard errors" because the u values are not
estimates in the technical sense (this is the same song-and-dance as
when we call something a "BLUP" rather than an "estimate").

When you call ranef(.) with condVar=TRUE (which is the default), the
conditional covariance matrices are indeed computed and returned, albeit
in a rather inconvenient form: attr(ranef(my_model)$grpvar,"postVar") is
a 3D array of stacked covariance matrices, or a list of 3D arrays (see
?ranef.merMod). If you just want the conditional standard deviations
for each random effect value, as.data.frame(ranef(my_model)) will
probably be more convenient.

(Again, the details behind this computation are explained in the JSS
lmer paper ...)

Note that these SDs **are not appropriate for testing hypotheses
about individual levels of the random effect**: one of the tradeoffs of
using REs in the first place is that you trade parsimony for the ability
to do hypothesis testing on individual predictions.

For your second question, I'm not sure what you mean by the
"conditional mean and conditional variance of the random effect"; can
you explain further/give an example?

Ben

On 2/15/21 3:24 PM, Ravi Varadhan wrote:

Dear Ben,
Thanks for your response. I went back and looked at the draft JSS
paper you sent me. It does describe how the random effects are
predicted as conditional modes, using a penalized, iteratively
weighted least squares. However, I still have some questions. Why is
the penalty term ||u||^2 added? What does this mean? Does glmer then
provide standard errors for the predicted random effects (I don't
think it does)?

One more question: it would be nice to also have an option for
conditional mean and conditional variance of the random effect,
although conditional variance would underestimate the true variance
of the prediction.

Thank you,
Ravi

Prediction of random effects in glmer()

Thread (12 messages)