6
$\begingroup$

I have an experiment in which I presented multiple stimuli to participants and wanted to control for the order in which the stimuli for shown. I am curious if it's possible to only account for order in the random effects without including a term in the fixed effects (i.e., since order in this case is a nuisance variable). The question would be if the second formula below is valid or not:

mod1 <- lmer(Response ~ Order + (Order|Subject), data = dat)
mod2 <- lmer(Response ~ 1 + (Order|Subject), data = dat)

In terms of a reproducible example, the sleepstudy data could be used:

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
fm2 <- lmer(Reaction ~ 1 + (Days | Subject), sleepstudy)

The same question would apply: would fm2 be reasonable fit under the assumption that I did not want the random slope for Days within Subject to be contingent on the fixed effect for Days. In other words, for fm2, the random effects would account for all aspects of the data relative to Subject besides the fixed effect intercept. When plotting the random slopes with sjPlot, the results seem reasonable given the expected outcome:

library(sjPlot)
sjp.lmer(fm1, type = "rs.ri")

fm1: lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

sjp.lmer(fm2, type = "rs.ri")

fm2: lmer(Reaction ~ 1 + (Days | Subject), sleepstudy)

library(ggplot2)
ggplot(sleepstudy) + geom_line(aes(x = Days, y = Reaction, color = Subject))

ggplot(sleepstudy) + geom_line(aes(x = Days, y = Reaction, color = Subject))

$\endgroup$
  • $\begingroup$ Are you sure that your two plots are not transposed? I would guess that your first plot, where the mean slope seems to be $0$, is the one without the fixed effect, while the second, where the mean slope seems to be positive, is the model with the fixed effect. $\endgroup$ – Matthew Drury Mar 1 '16 at 18:27
  • 1
    $\begingroup$ The plots are definitely in the correct order. The reason the one higher up (i.e., the plot of fm1) has a mean slope that appears to be 0 is because the fixed effect for Days is essentially "centering" the random effects in the plot (i.e., the random slopes in the first plot must be considered relative to the fixed effect slope for Days, which is positive overall). By contrast, the random slopes in fm2 do not include a fixed effect for the slope (i.e., Days), so here the slopes have an overall positive effect because the fixed main effect is not included, so the random slopes appear positive. $\endgroup$ – David Dobolyi Mar 1 '16 at 20:01
  • $\begingroup$ That makes perfect sense, you are, of course, correct. I think I was misinterpreting your plots as the posterior predictions from the model, though even in that case my conclusion was dubious. Sorry about that. I really like question by the way. $\endgroup$ – Matthew Drury Mar 1 '16 at 20:08
  • $\begingroup$ Thanks Matthew! And definitely no need for any apology--I've spent a ton of time wrapping my head around this, and what these plots show/ought to show is not immediately intuitive. In fact, that's the main reason I'm asking the clarification here about the whole random slopes without fixed slopes question, since I can't get a clear answer on whether or not it's reasonable in general. $\endgroup$ – David Dobolyi Mar 1 '16 at 20:24
  • $\begingroup$ Is your goal inference? $\endgroup$ – Roland Mar 2 '16 at 8:43
3
$\begingroup$

I believe this question to be very similar to the often wondered "must one always include an intercept term in a linear regression", for which the agreed upon answer is "yes, unless you have an extremely good reason not to".

I tried to think through what would happen without the fixed effect term before running any experiment. Let's write your two models out in detail. The first, with the fixed effect slope, is

$$ y \sim N(\mu_{\alpha} +\alpha_{[i]} + (\mu_{\beta} + \beta_{[i]}) x, \sigma) $$ $$ \alpha \sim N(0, \sigma_{\alpha}) $$ $$ \beta \sim N(0, \sigma_{\beta}) $$

where $x$ is the number of days, and we have a random intercept $\alpha_{[i]}$, and a random slope $\beta_{[i]}$ for each subject. In the other case, where there is no fixed slope, the model is

$$ y \sim N(\mu_{\alpha} + \alpha_{[i]} + \beta_{[i]} x, \sigma) $$ $$ \alpha \sim N(0, \sigma_{\alpha}) $$ $$ \beta \sim N(0, \sigma_{\beta}) $$

The difference is that in the second model, we a priori assume that the mean of the random slopes is zero. This means, we expect the slopes associated to the various subjects to distribute evenly around a slope of $0$ (for example, half should be negative and half positive).

Now, in the model on your data this does not seem to be true. In your second plot the estimated slopes within each subject are all positive. It looks like this model is invalid for your data. The inclusion of the fixed slope includes the mean of the subject-wise slopes as a degree of freedom, and in this plot you see the random slopes cluster evenly around zero, as you would like.

As for inference from the parameters in your model, I believe this misstatement of the model will cause the following parameter estimates to be bias

  • The subject-wise slopes will be biased towards zero, because the assumption of mean zero in the likelihood will pull them towards zero.
  • The estimated standard deviation of the random slopes will be too large, because inflating this parameter lets the slopes cluster around their true, non-zero mean without being penalized so severely.

Here I'll create some simulated data where the true subject-wise mean slope is non-zero

library("lme4")
library("arm")
set.seed(154)

N_classes = 50
N_obs <- 10000

random_intercepts <- structure(
  rnorm(N_classes), names = as.character(1:N_classes)
)

random_slopes <- structure(
  rnorm(N_classes, mean = 1), names = as.character(1:N_classes)
)

classes <- sample(as.character(1:N_classes), size = N_obs, replace = TRUE)
x <- runif(N_obs)
y <- random_intercepts[classes] + random_slopes[classes] * x + rnorm(N_obs)

df <- data.frame(class = factor(classes), x = x, y = y)

The first model estiamtes all true parameters well

> M <- lmer(y ~ x + (x | class), data = df)
> display(M)
lmer(formula = y ~ x + (x | class), data = df)
        coef.est coef.se
(Intercept) 0.01     0.15   
x           1.02     0.15   

Error terms:
 Groups   Name        Std.Dev. Corr 
 class    (Intercept) 1.03          
          x           1.01     0.19 
 Residual             1.00       

Look's like here all the parameters are estimated well, including the standard deviation of the random slopes.

Here's the model without the fixed slope

> N <- lmer(y ~ (x | class), data = df)
> display(N)
lmer(formula = y ~ (x | class), data = df)
coef.est  coef.se 
   -0.14     0.15 

Error terms:
 Groups   Name        Std.Dev. Corr 
 class    (Intercept) 1.04          
          x           1.43     0.24 
 Residual             1.00     

The estimate of the random slope standard deviation is 1.43, confirming my intuition that it would be biased to be too large.

The mean of the subject-wise slopes in the model M comes out well

> mean(fixef(M)["x"] + ranef(M)$class$x)
[1] 1.015418

It doesn't seem like my intuition was quite correct on the other model

> mean(ranef(N)$class$x)
[1] 0.9858566

It looks like the model took fitting the data a bit more seriously than making sure the normality of random slope assumption was totally met. Altogether, it looks like the inflation of the random slope standard deviation is the most serious issue.

$\endgroup$
  • $\begingroup$ This is very interesting, and thank you for taking the time to follow up. I do still wonder about aspects of your answer though, and also wanted to clarify a few things: 1) you mentioned that all subject slopes are positive in the second model (fm2), but in fact, one is negative (the teal one), and the random effect plot closely matches a quick spaghetti plot of the raw data (as expected--I've added this plot to the first post); 2) I need to add this to a separate comment (next up)... $\endgroup$ – David Dobolyi Mar 7 '16 at 22:15
  • $\begingroup$ 2) everything else you say makes sense, and I see your point regarding the ranef SD; however, I still wonder if I had another fixed effect term, e.g., Manipulation (a 2-level factor), would this fixed effect be interpretable regardless of how the random slope for Days is modeled (i.e., with vs. without the fixed effect included). Stated another way, I wonder how reasonable inference is for a given fixed effect parameter in a case like this: mod1 <- lmer(Response ~ Order + Manipulation + (Order|Subject), data = dat) vs. mod2 <- lmer(Response ~ Manipulation + (Order|Subject), data = dat) $\endgroup$ – David Dobolyi Mar 7 '16 at 22:19
  • $\begingroup$ I think I can respond to 1) at the moment. I'm not surprised that the slopes match a quick plot of your raw data given what happened at the end of my answer. Looks like the multilevel model puts more weight into fitting the data than sticking to the distributional assumptions underlying the model. You're right that one is negative, so i'll amend to "almost all" : ) $\endgroup$ – Matthew Drury Mar 7 '16 at 22:22
  • $\begingroup$ Thanks Matthew! =) Indeed, it makes sense given your very helpful/thorough reply. I suppose that's why I'm still curious about a term like Manipulation (the made up fixef) and its interpretation independent of Order being modeled as fixef/ranef or just ranef. Inclusion of the Order fixef obviously has repercussions on how we think of/interpret the ranef Order slope term (obviously if that interpretation was important, including the fixef makes much more sense as you point out), yet given both models fit well regardless of how we handle Order perhaps Manipulation can be interpreted either way. $\endgroup$ – David Dobolyi Mar 7 '16 at 22:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.