# Repeated measures analysis: why nest experimental factors within subject factor?

Consider a pure repeated measures design, with (let's say) 3 experimental within-subject factors A, B, and C, and (for simplicity) 2 levels per factor. So we have 2*2*2 = 8 measurements per subject.

Now I would like to test the fixed effects with a linear mixed effects model. I have read in several sources (for example Andy Field's Book "Discovering Statistics using R", and this site: http://www.jason-french.com/tutorials/repeatedmeasures.html ) that with lme, one should use the following syntax:

model <- lme(dv ~ A*B*C, random = ~1|id/A/B/C)


However, I do not understand why you would "nest" the factors within the subject in the random part of the model, and not just use (1|id). What is the point of this, and what does it do?

Conceptually, I don't understand why one would nest the experimental fixed factors within the random subject factor. The way I understood nesting until now, you would only use it to account for the fact that certain lower factor levels only exist within certain higher factor levels - like pupils within classes within schools within cities, etc. How does this apply to a repeated measures design with fully crossed within-subject factors?

Mathematically, the way I understood this is that such a model would first estimate a random intercept for each subject, capturing random differences in the average values of the dependent variable between subjects. So in the case of, (let's say) 20 subjects, we get 20 different random intercepts. Then, apparently, the model estimates random intercepts for each combination of subject and level of factor A (resulting in 40 random intercepts), then for each combination of subject, factor A and factor B (80 random intercepts), all the way down to the most specific level, where we get as many estimated random intercepts as we have measurements (160). What is the point of this, and why would we not only estimate a random intercept per subject (1|subject)? Also, wouldn't all of these random intercepts together explain the dependent variable nearly perfectly, and leave little to nothing to be explained by the fixed effects?

Lastly, my intuition tells me that these random intercepts should at least partially explain the same information as would be captured by entering random slopes of the experimental factors into the model. Is that correct?

• This is a purely theoretical question - I don't have a specific real situation. What I am interested in is why this is recommended as a standard way to analyze repeated measures data (at least by some authors), and what it means conceptually and mathematically to have a random intercept on each of these levels of "nested" experimental factors. – Lukas McLengersdorff Jun 13 '17 at 9:11
• +1. This is an excellent question. It's very closely related to something I asked here some time ago: stats.stackexchange.com/questions/232109 - but unfortunately there is no satisfactory answer in that thread (existing answers might only add to the confusion, so be careful). I am planning to write an answer there myself at some point, but I still need to do some investigations beforehand... – amoeba Jun 14 '17 at 21:19
• Still, some partial answers. (i) You are right, nothing is really "nested" here. We use the "nesting operator" because we want to have (1|id) + (1|id:A) + (1|id:B) + (1|id:A:B) + ... etc. structure, but it would not be correct to say that A,B,C are nested in subjects. They are not. (ii) We want (1|id) + (1|id:A) + ... structure because it mimics the approach of repeated measures ANOVA (RM-ANOVA). I can't properly explain why RM-ANOVA does that though. (iii) Yes, you can have random slopes (A*B*C - A:B:C | id) instead. But it's very different model when A,B,C have many levels. – amoeba Jun 14 '17 at 21:23
• Thanks! I did some research myself in the meantime, and found that with lme4, you can nearly perfectly replicate the RM ANOVA (with three factors A,B,C, like in my example) with the model 'y~ABC + (1|id) + (1|id:A) + (1|id:B) + (1|id:C) + (1|id:A:B) + (1|id:A:C) + (1|id:B:C)'. Or, put more simply, a random intercept for subject, plus a random intercept for every interaction of subject and every main effect and interaction of the within-subject factors, except for the highest order interaction. – Lukas McLengersdorff Jul 10 '17 at 15:20
• @amoeba no worries. I didn't realise this was an old question until I answered it and got a revival badge :D – Robert Long Jul 15 '20 at 12:24

This is a very interesting question. I have been thinking this and related things for a long time.

For me the key to understanding this is to realise that: Random intercepts for a grouping factor are not always sufficient to capture the random variation in the data that is in excess of the residual variation. Because of this, we sometimes see models with random intecepts for interactions between a fixed factor and a grouping variable (and even sometimes random intercepts just for a fixed factor). Generally we advise that a factor can be fixed or random (intercepts) but not both - however there are important exceptions of which your example here is one.

Another hindrence to understanding this is for people like me who come from a multilevel modelling / mixed model background in observational social and medical sciences, we are often caught up in thinking about repeated measures and nesting vs crossed random effects without an understanding that things are a bit different in experimental analysis. More on this a bit later.

Judging by the comments we have both discovered the same thing. In the context of a repeated measures ANOVA, if you want to obtain the same results with lmer then you fit:

 y ~ A + B + (1|id) + (1|id:A) + (1|id:B)


where I have discarded factor C without loss of generality.

and the reason why some people specify 1|id/A/B is that they are using nmle:lme and not lme4:lmer. I am not sure why this is needed in lme() but I am fairly sure that to replicate a repeated measures anova - where there is variation for each combination of id and the factors - then you fit the model above in lmer(). Note that (1|id/A/B) seems similar, however it is wrong because it would also fit (1|id:A:B) which is indistinguishable from the residual variance (as noted in your comment also).

It is important to note (and therefore worth repeating) that we only fit this type of model where we have reason to believe that there is variation for each combination of id and the factors. Typically with mixed models we would not do this. We need to understand experimental design. One type of experiment where this is common is the so-called split-plot design where blocking has also been used. This type of experimental design employs randomisation at different "levels" - or rather different combinations of factors, and this is why analysis of such experiemnts often includes random intercept terms that at first glance seem odd. However, the random structure is a property of the experimental design and without knowledge of this, it is virtually impossible to select the correct structure.

So, with regards to the your actual question where the experiment has a repeated factorial design, we can use our friend, simulation, to investigate further.

We will simulate data for the models:

 y ~ A + B + (1|id)


and

 y ~ A + B + (1|id) + (1|id:A) + (1|id:B)


and look at what happens when we use both models to analyse both datasets.

set.seed(15)
n <- 100 # number of subjects
K <- 4 # number of measurements per subject

# set up covariates
df <- data.frame(id = rep(seq_len(n), each = K),
A = rep(c(0,0,1,1), times = n),
B = rep(c(0,1), times = n * 2)
)

#
df$$y <- df$$A + 2 * df$$B + 3 * df$$id + rnorm(n * K)

m0 <- lmer(y ~ A + B + (1|id) , data = df)
m1 <- lmer(y ~ A + B + (1|id) + (1|id:A) + (1|id:B), data = df)
summary(m0)


m0 is the "right" model for these data because we simply created y with fixed effects for id (which we capture with random intercepts) and unit variance. This is a bit of an abuse but it is convenient and does what we want:

 Groups   Name        Variance Std.Dev.
id       (Intercept) 842.1869 29.0205
Residual               0.9946  0.9973
Number of obs: 400, groups:  id, 100

Fixed effects:
Estimate Std. Error t value
(Intercept) 50.47508    2.90333   17.39
A            1.01277    0.09973   10.15
B            2.06675    0.09973   20.72


as we can see, we recover unit variance in the residual and good estimates for the fixed effects. However:

> summary(m1)

Random effects:
Groups   Name        Variance  Std.Dev.
id:B     (Intercept) 0.000e+00  0.0000
id:A     (Intercept) 8.724e-03  0.0934
id       (Intercept) 8.422e+02 29.0204
Residual             9.888e-01  0.9944
Number of obs: 400, groups:  id:B, 200; id:A, 200; id, 100

Fixed effects:
Estimate Std. Error t value
(Intercept) 50.47508    2.90334   17.39
A            1.01277    0.10031   10.10
B            2.06675    0.09944   20.78


This is a singular fit - zero estimates for the variance of the id:B term and close to zero for id:A - which we happen to know is correct here because we didn't simulate any variance for those "levels". Also we find:

> anova(m0, m1)

m0: y ~ A + B + (1 | id)
m1: y ~ A + B + (1 | id) + (1 | id:A) + (1 | id:B)
npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
m0    5 1952.8 1972.7 -971.39   1942.8
m1    7 1956.8 1984.7 -971.39   1942.8 0.0052  2     0.9974


meaning that we strongly prefer the (correct) model m0

So now we simulate data with variation at these "levels" too. Since m1 is the model we want to simlulate for, we an use it's design matrix for the random effects:

# design matrix for the random effects
Z <- as.matrix(getME(m1, "Z"))

# design matrix for the fixed effects
X <-  model.matrix(~ A + B, data = df)

betas <- c(10, 2, 3) # fixed effects coefficients
D1 <- 1 # SD of random intercepts for id
D2 <- 2 # SD of random intercepts for id:A
D3 <- 3 # SD of random intercepts for id:B

# we simulate random effects
b <- c(rnorm(n*2, sd = D3), rnorm(n*2, sd = D2), rnorm(n, sd = D1))
# the order here is goverened by the order that lme4 creates the Z matrix

# linear predictor
lp <- X %*% betas + Z %*% b

# add residual variance of 1
df\$y <- lp + rnorm(n * K)

m2 <- lmer(y ~ A + B + (1|id) + (1|id:A) + (1|id:B), data = df)
m3 <- lmer(y ~ A + B + (1|id), data = df)
summary(m2)


'm2 is the corect model here and we obtain:

Random effects:
Groups   Name        Variance Std.Dev.
id:B     (Intercept) 6.9061   2.6279
id:A     (Intercept) 4.4766   2.1158
id       (Intercept) 2.9117   1.7064
Residual             0.8704   0.9329
Number of obs: 400, groups:  id:B, 200; id:A, 200; id, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  10.3870     0.3866  26.867
A             1.8123     0.3134   5.782
B             3.0242     0.3832   7.892



the SD for the id intercept is a little high, but otherwise we have good estimates for the random and fixed effects. On the other hand:

> summary(m3)

Random effects:
Groups   Name        Variance Std.Dev.
id       (Intercept) 6.712    2.591
Residual             8.433    2.904
Number of obs: 400, groups:  id, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  10.3870     0.3611  28.767
A             1.8123     0.2904   6.241
B             3.0242     0.2904  10.414



although the fixed effects point estimates are OK, their standard errors are larger. The random structure is, of course, completely wrong. And finally:

> anova(m2, m3)

Models:
m3: y ~ A + B + (1 | id)
m2: y ~ A + B + (1 | id) + (1 | id:A) + (1 | id:B)
npar    AIC    BIC   logLik deviance Chisq Df Pr(>Chisq)
m3    5 2138.1 2158.1 -1064.07   2128.1
m2    7 1985.7 2013.7  -985.87   1971.7 156.4  2  < 2.2e-16 ***


showing that we strongly prefer m2`

• +1 and a very good answer. But I wish you wrote a bit more in the paragrah where you mention split-plots, ideally giving examples. What is the specific randomization that is used there and why is it different from what you are used to (like you say earlier). Perhaps then you could come back to it in the end of your answer and explain what would the two simulations correspond to in real world. When would the data be expected to be more like in simulation 2. vs simulation 1? – amoeba Jul 15 '20 at 12:21
• @amoeba Thanks ! I didn't go into detail about split plot because the question was specifically about RM ANOVA. However, the similarities are striking. I'm still doing my research on split plot and other designs so I'm not yet comfortable to go into much more detail. I probaly won't add it to this answer - I will look for a question about split plots specifically and let you know if I answer it. I might also have a go at a new answer to your old question - it's been bugging me since we talked at length in comments about it 4 years ago and never found a satisfactory solution/explanation :/ – Robert Long Jul 15 '20 at 12:36
• Okay - fair enough not to add details about split-plot specifically, but I still think that having some discussion on when the repeated-measures data are expected to be more like sim1 or more like sim2 would be important. Currently you are simulating data for sim2 using the random effect matrix which is quite opaque and so the reader is left wondering when would the real-world data be like that. (And yes, absolutely please do eventually post a new answer to that old Q of mine!) – amoeba Jul 15 '20 at 13:36
• +1 and accept! Very positively surprised that my old question still got an answer :) – Lukas McLengersdorff Jul 18 '20 at 13:44
• Lukas, I was on the same chapter of Andy Field's book that you were on 3 years ago and ended up having the same question. I set up a bounty and fortunately Robert answered it. Thank you Robert! – Neal Aug 2 '20 at 21:46