[R-sig-ME] how to compute Confident intervals for BLUPS for lme function in nlme library
Reinhold Kliegl
reinhold.kliegl at gmail.com
Sun Jan 18 18:20:23 CET 2009
On Fri, Jan 16, 2009 at 12:09 PM, Nick Isaac <njbisaac at googlemail.com> wrote:
> 1) Figure 4 in your manuscript is very similar to the plot produced by the
> code in your example. You desribe the bars around the conditional modes as
> '95% prediction intervals'. Are these synonymous with the posterior
> variances of the postVar attribute, or did you apply some kind of
> transformation?
They are the default output of the dotplot() function; no
transformations necessary.
>
> 2) Is it appropriate to cite your manuscript? If so, should the citation
> include information not on the website (e.g. journal title)?
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.
>
> 3) You recently posted a question about sorting the output from ranef() for
> a dotplot (13/10/08). Did you ever figure out how to do this? I guess one
> could always extract the components, sort them directly and plot them in
> another function...
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
>
>
> 2009/1/14 Reinhold Kliegl <reinhold.kliegl at gmail.com>
>
>> Well, intervals() is a pretty good nlme response, I think, even if it
>> was not exactly answering the question. The lme4 equivalent is
>> dotplot() 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
>>
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list