[R-sig-ME] need help with predicted values from lmer with newdata

[Paul Johnson]

On Sat, Jun 23, 2012 at 6:27 AM, Ben Bolker <bbolker at gmail.com> wrote:
> Joshua Wiley <jwiley.psych at ...> writes:

>
>  Paul, would you be willing to generate a trimmed-down version of your
> example that simply generates multi-trait data?
>

I noticed Joshua already working on help page. Great. I don't want to
duplicate his effort, i'm going to keep working on example usages
which Ben may be able to use in the documentation.  I think this
example will be good for lmer, once it is tighter.

I'm attaching most recent longish version to this note. I've already
cut out some crap, but it is not tightened down to minimums yet
because I am stuck on 2 issues. After I understand those, I can cut
out some more crap. The surprise I had last friday with the predicted
values has made me feel cautious/worried.  Also, I will supply a
version that does not depend on my other package "rockchalk".

First, the merMod question.

I want to script up a simulation that fits lots of lmer and calculates
the "bias" and variance of the estimates of the random effects.  If
you run the script I attach, you get a rough idea of how one
particular run would be analyzed. But I'll scale this up, eventually.

What is the recommended way to extract estimated standard deviations
of random effects for a model like this:

mm3 <- lmer( y3 ~ x1 + x2 + x3 + (1|Mind) + (-1 + x1 | Mind),
data=dat, verbose=3)

Don't say "allow correlation between random intercept and random
slope" right now. I allow that later. Here we are wondering, "if we
ignore correlation between random effects, what happens?"

Here's what I tried, after reading lme4 source code. Problem, the
names on the first two list elements are "Mind". So I can't get the
std dev of the slope of x1

## The usual spiel goes
summary(mm3)
mm3VarCorr <- VarCorr(mm3)
mm3VarCorr

## Yank out the standard deviation estimates. idiom stolen from
formatVC in lmer code.
mm3reStdDev <- c(lapply(mm3VarCorr, attr, "stddev"), list(Residual =
attr(mm3VarCorr, "sc")))
mm3reStdDev

## Calculate bias of this run's estimate:
mm3reStdDev[["Mind"]]["(Intercept)"] - STDEb0
##mm3reStdDev[["Mind"]]["x1"] - STDEb1
## Can't figure how to get STD(b1)

I can get it if I access the list members with their numbers

attr(mm3VarCorr[[2]], "stddev")

but in a general setting, this is not good because it pre-supposes I
know that mm3VarCorr has 2 elements in that list.


Second, this is a substantive question.  I've run this script 100s of
times, plugging in various values of the variance of the X's, the
variance of the b's, and the residual standard deviation.

In many cases, the estimates from lmer come out almost exactly right,
but sometimes I fiddle the parameters so that some nonsense appears to
happen.  In the attached script, there is some nonsense, I expect you
will see if you run it.

Here there are 2 nonsensical results:

mm2 <- lmer( y2 ~ x1 + x2 + x3 + (1 | Mind), data=dat)

gives variance estimate for Mind of 0.  So, apparently, there's some
ratio of STDEb0/STDE at which STDEb0 becomes undetectable?  Possible.

And in mm4, where it says the correlation of the random effects is 1.00.

It doesn't always happen, depending on the relative sizes of STDE,
STDEb0, and STDEb1, I can produce things I understand.

I started calculating the correlation of the fitted and the observed
values to make sense out of this,

cor(fitted(mm2), dat$y2)

but don't quite get it. Yet.

Come to think of it, I will just attach the transcript so you don't
have to run this.

But you should run it to see my graphs, which are lovely.

-- 
Paul E. Johnson
Professor, Political Science    Assoc. Director
1541 Lilac Lane, Room 504     Center for Research Methods
University of Kansas               University of Kansas
http://pj.freefaculty.org            http://quant.ku.edu