I've seen questions about fitting a slope in a linear mixed model as a random slope or as both a fixed and random slope come up on this list before (see these two threads:
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q1/017750.html
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q1/017495.html ).
This just came up for me again when I was reviewing an analysis, and I was hoping to get some further clarification so I can make a reasonable argument for changing the analysis (or not).  Maybe I'm over-thinking this and none of this matters if the focus of the analysis is on inference to any fixed effects and that is essentially unchanged between models?
When we write out the definition for our random variables in a mixed model, we specify them as coming from distributions with a mean of 0.  I hadn't ever really thought this through explicitly to the idea that the fitting method relies on that mean of 0 for the random effects, but after reading the threads above and reviewing some theory on REML, that appears to be true.  Is that right?  Theoretically, would it be possible to change the definition of the random effects to a distribution with a mean of the mean random effect instead of 0 (outside using a MCMC fitting method)?  It seems easy to simply shift the location of the distribution of a random variable, which wouldn't change the variances in a normal distribution, but I bet this is oversimplifying things.
I've set up an example with the sleepstudy dataset to explore what happens when I fit a model with only random slopes to one fit "correctly" with both a fixed and random slope.
First fit with an uncorrelated random intercept and slope, and a fixed slope for "Days":
> fit1 = lmer(Reaction ~ 1 + Days + (1|Subject) + (0 + Days|Subject), data = sleepstudy)
The mean of the random slope effects is 0:
> mean(ranef(fit1)[["Subject"]]$Days)
[1] 2.479219e-13
The estimated fixed effects:

> fixef(fit1)

(Intercept)        Days

  251.40510    10.46729
The variances from the model:

Groups   Name        Variance Std.Dev.

 Subject  (Intercept) 627.581  25.0516

 Subject  Days         35.858   5.9882

 Residual             653.583  25.5653

Then I remove the fixed slope.  The main things that stand out to me are that the mean of the random slope effects is not the same as the fixed effect slope, and that the slope variance is very different.

> fit2 = lmer(Reaction ~ 1 + (1|Subject) + (0 + Days|Subject), data = sleepstudy)

> mean(ranef(fit2)[["Subject"]]$Days)

[1] 9.908601

> fixef(fit2)

(Intercept)

   253.9192

> VarCorr(fit2)

 Groups   Name        Variance Std.Dev.

 Subject  (Intercept) 610.37   24.706

 Subject  Days        140.95   11.872

 Residual             655.65   25.606

It was actually seeing the difference in the variances that really surprised me.  I fooled around with calculating variances from the random effects to see if I learned anything.  This is a bit naīve, but it does look like that 0 mean for the random slopes is important here.
With fixed slope:

> var(ranef(fit1)[["Subject"]]$Days)

[1] 30.4038

> sum((ranef(fit1)[["Subject"]]$Days - 0 )^2)/17

[1] 30.4038

Without fixed slope:

> var(ranef(fit2)[["Subject"]]$Days)

[1] 38.70345

> sum((ranef(fit2)[["Subject"]]$Days - 0 )^2)/17

[1] 142.6591

I'd be happy to hear how anyone one else has addressed this (in a review or otherwise),
Ariel


	[[alternative HTML version deleted]]

