Skip to content

how to compute Confident intervals for BLUPS for lme function in nlme library

9 messages · Greg Lee, Reinhold Kliegl, Douglas Bates +3 more

Original
Reinhold Kliegl
# 
Well, intervals() is a pretty good nlme response, I think, even if it
was not exactly answering the question. The lme4 equivalent is
doplot() which produces a so-called "caterpillar plot" of conditional
means.

For example:
?ranef

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
rr1 <- ranef(fm1, postVar = TRUE)
dotplot(rr1,scales = list(x = list(relation = 'free')))[["Subject"]]

It is important to be clear about the distinction between random
effect estimates of variances and correlations, provided as part of
the model output and the conditional means for groups/units based on
("predicted with") them. At least, correlations computed from
conditional means can differ strongly from the correlation estimates.
The short story: You should only trust the model estimates. As Greg
Lee said, conditional means are to be handled with care; you can not
apply any of the usual inferential statistics to them. Nevertheless,
the caterpillar plots are highly informative and diagnostic about, for
example, whether you need a random effect to account for unit-by-unit
variability. They may also lead you to  discover distinct subgroups
suggestive of a fixed effect (for a future model). There is a
manuscript at the top of my publications page for download (and
constructive feedback) walking through an example.

Reinhold Kliegl
On Wed, Jan 14, 2009 at 2:03 AM, Greg Lee <sp8ial at gmail.com> wrote:
Douglas Bates
# 
On Wed, Jan 14, 2009 at 12:50 AM, Reinhold Kliegl
<reinhold.kliegl at gmail.com> wrote:
There's a typo there.  The function is dotplot(), not doplot().
I agree with all of the above.  As Reinhold says, the BLUPs are the
conditional mean vector of the distribution of the random effects
given Y = y, the observed data and evaluated at the parameter
estimates.  Sometimes I use the term "conditional modes".  For a
linear mixed model the conditional distribution is multivariate
Gaussian so the mean and the mode coincide.  In generalized linear
mixed models or nonlinear mixed models they do not necessarily
coincide and what is calculated is the conditional mode.

The intervals shown in the dotplot come from the evaluation of the
conditional means and the conditional variances of the random effects
given the observed data.  You can obtain the conditional variances by
including the argument postVar = TRUE in a call to ranef. (The
argument is "postVar" and not "condVar" because the person who
requested it used the term "posterior variance", which is what this is
if you take a Bayesian or empirical Bayes approach.)
Nathan Leon Pace, MD, MStat
# 
I ask the indulgence of the group for some assistance in a short
lecture/seminar.

I'm suppose to give a 30 minute overview of mixed models (GL, NL, and GNL)
to a small group of mostly physician researchers - some of whom have
experience using the NONMEM package to do pharmacokinetic/pharmacodynamic
population modeling.

I looked at some my favorite archives/blogs/academic sites, but can't find a
existing outline, syllabus, .ppt file, etc, for such a purpose.

If someone can point me toward an example, I'd appreciate seeing how someone
else did a similar presentation.

Nathan
1 day later
2 days later
Reinhold Kliegl
#
On Fri, Jan 16, 2009 at 12:09 PM, Nick Isaac <njbisaac at googlemail.com> wrote:

            

      
      
      
    

        
They are the default output of the dotplot() function; no
transformations necessary.

            

      
      
      
    

        
Actually, I am require to include "This is a preprint of an article
submitted for consideration in Visual Cognition (c) 2008 Taylor &
Francis" on the title page.

            

      
      
      
    

        
Actually, I was able to help myself. I copied the original dotplot()
function to a private dotplot.RK() function.  There,  I added the
argument "refvar" (reference variable) to the dotplot function and
replaced:

ss$.nn <- rep.int(reorder(factor(rownames(x)), x[[1]]), ncol(x))

with

ss$.nn <- rep.int(reorder(factor(rownames(x)), x[[refvar]]), ncol(x))

This worked for my models, but the solution will not hold for the
general form of the model, as detailed in the following comment by
Douglas Bates in an off-line exchange:

"The tricky part, as always, is whether that argument can be applied
to general forms of the model.  In general the ranef.mer class is a
list of arrays in which the numbers of rows and columns do not need to
be, and usually aren't, consistent.  The names of the columns don't
need to be consistent either.

One could specify the refvar as a number or as a name and get it to
work there but there should be a failsafe line that specifies what to
do if the number is greater than the number of columns or if the name
is not in the column names."

So some caution is in order.

Reinhold Kliegl