[R-sig-ME] Prediction of random effects in glmer()

[Ben Bolker]

    I don't know. I suspect that they *will* have nominal or 
close-to-nominal coverage (a restatement of your comment below) under 
reasonably weak asymptotic assumptions.
   It wouldn't be too hard to do some simulations to check the coverage 
for some simple examples, but I shouldn't succumb to the temptation 
right now.

   cheers
    Ben

On 2/15/21 4:58 PM, Andrew Robinson wrote:
> 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 using unimelb.edu.au
> Website: https://researchers.ms.unimelb.edu.au/~apro@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 using 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 using unimelb.edu.au
>>> Website: https://researchers.ms.unimelb.edu.au/~apro@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 using 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
>>>>> ________________________________
>>>>> From: Ravi Varadhan
>>>>> Sent: Thursday, February 11, 2021 8:30 PM
>>>>> To: r-sig-mixed-models using r-project.org <r-sig-mixed-models using r-project.org>
>>>>> Subject: Prediction of random effects in glmer()
>>>>>
>>>>> Hi,
>>>>> I would like to know how the prediction of random effects is done in
>>>>> the GLMM modeling using the lme4::glmer function, i.e. how the
>>>>> BLUP-like predictions are made in the glmer() function?
>>>>>
>>>>> Does it use frequentist prediction or empirical Bayes or full Bayes
>>>>> posterior? Is there any documentation of the prediction methodology?
>>>>>
>>>>> Thanks in advance.
>>>>>
>>>>> Ravi
>>>>>
>>>>> [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-mixed-models using r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models 
>>>>> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>>>>> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models 
>>>>> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>>
>>>>>
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models using r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models 
>>>> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>>>> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models 
>>>> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>>
>>>>
>>