Below are two long questions about variance components and multilevel models.

Shinichi Nakagawa forwards these questions from Hoger Schielzeth:

(A)
If you got a split-plot design, where treatments A are applied to individuals and your are also interested in the within-individual covariate C (and possibly the AxC-interaction). I think this is a common design. E.g. you are measuring egg sizes in a number of female house sparrows (5 eggs from one clutch). Half of the subjects get supplementary feeding. Then you want to test if the increase (or decreases) in egg sizes depends on the treatment (i.e. if there is an interaction between laying order slope and supplementary feeding). I think most people would use a model like
lmer(eggSize ~ Treatment * LayOrder + (1|FemaleID)
and would consider
lmer(eggSize ~ Treatment * LayOrder + (LayOrder|FemaleID)
an extra sophistication that might reduce residual error. Would you agree on this?

Actually, your confidence intervals are much to narrow (of LayOrder and Treatment x LayOrder), when there is between-individual variation in slopes (which is probably the norm!) and you do not included random slopes. Gelman & Hill are not mentioning this problem, but you can see the same effect in their models (The radon-county-floor examples).

Do you know any treatment of this issue in the literature?

(B)
If you included random slopes: Do you know how to interpret the variance components? The problem is that intercept and residual variation are directly comparable, since they are on the same scale (the scale of the data), but the random slope variance is on a different scale (its just
the variation between slopes of individuals). Scaling the input makes it a bit better, but I think the problem still exists.

# Random effects:
# Groups Name Variance Std.Dev. Corr
# ClutchID (Intercept) 0.0212473 0.145764
# scale(LayOrder) 0.0033531 0.057906 -0.138
# Residual 0.0039384 0.062757

Shinichi writes:

The question (B) I thought is a wonderful point which Gelman & Hill do not really mention – I am interested to hear what would you do?

Actually, he later on rephrased the question (A) when I did not quite understand what he asked.

#########################################################
The easiest example, I came up with is a paired t test situation. What the paired t test actually does, is taking differences within individuals (and importantly individuals might vary in their overall level, i.e. their intercepts) and that it tests, whether the population of slopes (random slopes!) is different from zero. There are basically two mixed model alternatives to this.

(1) response ~ treatment + (1 | individual) – random intercepts
(2) repsonse ~ treatment + (treatment | individual) – random slopes

The latter is equivilant to a pair t test. The former performs about as bad as ignoring individual variation (i.e lm(response ~ treatment)). And still I think, most people would think model (1) does what they want. Check it out, that’s how it is! Actually, it controls for peudoreplication, by not for the pseudoreplication you really want to control for (the pseudoreplication of the slopes). Maybe in this case its easy enough, so people will agree on model (2), but if the design is slightly more complicated…
#####################################################

To my knowledge, we are kind of taught to think the model (1) is the equivalent of paired-t test.

I think that he is writing a short note for Evolutionary Biologists to make his point above (now everyone uses mixed models but probably nobody is aware of this problem).

My reply:

(A) I discuss many of these issues, in particular the difficulties with analyzing split-plot designs using standard software, in my 2005 paper on Anova; see page 6 of that paper. We also discuss split-plot designs in our book. I agree that it’s good to allow slopes to vary–as long as the program does not crash! It’s better to include more in your model, and a challenge with thes programs is to allow this. I’m thinking we should have an improved lmer that puts some weak prior information on the group-level covariance matrix to get stable estimates.

(B) If you scale the input variables, then your intercepts and slopes should be on comparable scales. We discuss this in our book (chapter 4, I believe), also see this paper.

(C) I recommend you display your fitted models using display() rather than print() or summary(). All will be more interpretable. Numbers like Variance = 0.0033531 are horrible. First, it’s way too many digits. Second, it’s on the wrong scale. It’s the sd that you want to look at, not the variance.

(D) In the simple example of the paired t-test, you can work with person-level summaries, which should simplify things a bit.