Jarrett Byrnes writes,

I’m currently attempting to work out how one would look at coefficients in a multilevel model that vary by more than one type of grouping after working through your ARM book. For example, let’s say I’m looking at how much algae a snail eats, and how this changes with snail size. The slope of this line can vary by both the species of snail, and, for example, then by the location from which I collected the snail.

Algae Eaten ~ N(B*Snail Mass, sigma_snail)

B = Mean Slope + Deviation due to Species + Deviation due to Collection Location.

Of course, each deviation would be N(0, sigma_grouping)
I would code this, in lmer, as

eaten.lmer<-lmer(Consume ~ Mass + 0 + (Mass-1|Species)+(Mass-1|Location) Note, I'm constraining the intercept. Biologically, a snail of mass 0 cannot eat any algae. And, I'm assuming linearity just for the sake of argument. So, a few questions: If I wanted to look at the error around each coefficient estimate, would I add errors from each of the random effects together? If I wanted to look at the slopes for each grouping individually, would I just look at the estimate of the fixed effect+deviation due to grouping with the standard error of just that grouping? Also, if you do not have complete crossing of all factors (i.e., let's say you didn't find some species at some locations), how would this effect your estimates- particularly if you were more concerned about determining the how slope varies for sampled snails, rather than using the information to make predictions? I've written a first stab at taking any lmer object and plotting each fixed effect, which I’m attaching. It may be a little clumsy currently, but it gets the job done, and allows one to avoid coding model specific plotting functions – although it currently suffers from my not really having any answers to the above questions. I’ll try and keep puzzling it through!

My comments:

1. Indeed, a species with mass 0 cannot eat any algae, but I’d still include an intercept (and I’d let that intercept vary by species and location). I don’t know about your example, but in general you can get nonzero asymptotes as x approaches zero.

2. Why not use the log scale? I can imagine that you might have more variability on the log scale at low levels, but you can handle that in the variance function.

3. The question about standard errors is a good one. The standard errors that come from lmer() are for individual coefficients, and I don’t think there’s a really easy way to combining.

With full Bayesian inference, standard errors come out automatically from the simulations. With an approximate algorithm such as used in lmer(), you can get reasonable standard errors that ignore uncertainty in the variance parameters. One solution is to use mcsamp() (our wrapper for the mcmcsamp() function that goes with lmer()) but then you’re relying on simulation and possibly not-entirely-debugged code. It is straightforward enough to get standard errors using matrix manipulation and weighted regression on the augmented data matrix. Maybe we should write a program to do this.

4. Finally, yes, I like the idea of plotting these things.