What is the minimum recommended number of groups for a random effects factor?

Question

I'm using a mixed model in R (lme4) to analyze some repeated measures data. I have a response variable (fiber content of feces) and 3 fixed effects (body mass, etc.). My study only has 6 participants, with 16 repeated measures for each one (though two only have 12 repeats). The subjects are lizards that were given different combinations of food in different 'treatments'.

My question is: can I use subject ID as a random effect?

I know this is the usual course of action in longitudinal mixed effects models, to take account of the randomly sampled nature of the subjects and the fact that observations within subjects will be more closely correlated than those between subjects. But, treating subject ID as a random effect involves estimating a mean and variance for this variable.

• Since I have only 6 subjects (6 levels of this factor), is this enough to get an accurate characterization of the mean and variance?

• Does the fact that I have quite a few repeated measurements for each subject help in this regard (I don't see how it matters)?

• Finally, If I can't use subject ID as a random effect, will including it as a fixed effect allow me to control for the fact that I have repeated measures?

Edit: I'd just like to clarify that when I say "can I" use subject ID as a random effect, I mean "is it a good idea to". I know I can fit the model with a factor with just 2 levels, but surely this would be in-defensible? I'm asking at what point does it become sensible to think about treating subjects as random effects? It seems like the literature advises that 5-6 levels is a lower bound. It seems to me that the estimates of the mean and variance of the random effect would not be very precise until there were 15+ factor levels.

Answer 1

Short answer: Yes, you can use ID as random effect with 6 levels.

Slightly longer answer: The @BenBolker's GLMM FAQ says (among other things) the following under the headline "Should I treat factor xxx as fixed or random?":

One point of particular relevance to 'modern' mixed model estimation (rather than 'classical' method-of-moments estimation) is that, for practical purposes, there must be a reasonable number of random-effects levels (e.g. blocks) — more than 5 or 6 at a minimum.

So you are at the lower bound, but on the right side of it.

Answer 2

In an attempt to figure out the minimum number of groups for a multilevel model I looked at the book Data Analysis Using Regression and Mulitilevel/Hierarchical models by Gelman and Hill (2007).

They appear to address this topic in Chapter 11, Section 5 (page 247) where they write that when there are < 5 groups then multilevel models typically add little over classical models. However, they appear to write that there is little risk to applying a multilevel model.

The same authors appear to return to this topic in Chapter 12, Section 9 (pages 275-276). There they write that advice on the minimum number of groups for a multilevel model is misguided. There they again say that multilevel models often add little over classical models when number of groups is small. However, they also write that multilevel models should do no worse than no-pooling regression (where no-pooling appears to mean that group indicators are used in classical regression).

On pages 275-276 the authors have a specific subsection for the case of one or two groups (e.g., male versus female). Here they write that they typically express the model in classical form. However, they state that multilevel modeling can be useful even with only one or two groups. They write that with one or two groups multilevel modeling reduces to classical regression.

My impression from this is that classical regression is one end of a continuum of models, i.e., a special case of a multilevel model.

Based on the above, my impression is that classical regression and multilevel modeling will return almost identical estimates when there are only two groups and that using multilevel models with only one, two, three, four, five or six groups is okay.

I will attempt to modify this answer in the future with R code and a small data set comparing estimates obtained with both approaches when using two groups.

Answer 3

For what it's worth, I did a bit of a simulation study to look at the stability of the variance estimate for a relatively simple LMM (using the sleepstudy dataset available through lme4). The first way generates all possible subject combinations for ngroups number of subjects, and refits the model for each possible combination. The second takes several random subsets of subjects.

library(lme4)
library(ggplot2)
library(tidyr)

m0 <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy,
           control = lmerControl(optimizer = "nloptwrap"))
# set the number of factor levels
ngroups <- 3:18 
# generate all possible combinations
combos <- lapply(X = ngroups, 
                 FUN = function(x) combn(unique(sleepstudy$Subject), x)) 

# allocate output (sorry, this code is entirely un-optimized)
out <- list(matrix(NA, ncol(combos[[1]]), 1), matrix(NA, ncol(combos[[2]]), 1),
            matrix(NA, ncol(combos[[3]]), 1), matrix(NA, ncol(combos[[4]]), 1),
            matrix(NA, ncol(combos[[5]]), 1), matrix(NA, ncol(combos[[6]]), 1),
            matrix(NA, ncol(combos[[7]]), 1), matrix(NA, ncol(combos[[8]]), 1),
            matrix(NA, ncol(combos[[9]]), 1), matrix(NA, ncol(combos[[10]]), 1),
            matrix(NA, ncol(combos[[11]]), 1), matrix(NA, ncol(combos[[12]]), 1),
            matrix(NA, ncol(combos[[13]]), 1), matrix(NA, ncol(combos[[14]]), 1),
            matrix(NA, ncol(combos[[15]]), 1), matrix(NA, ncol(combos[[16]]), 1))
# took ~ 2.5 hrs on my laptop, commented out for safety
#system.time(for(ii in 1:length(combos)) {
#    for(jj in 1:ncol(combos[[ii]])) {
#    sls <- sleepstudy[sleepstudy$Subject %in% combos[[ii]][,jj],]
#    out[[ii]][jj] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
#        }
#    })

# pad with zeros, not all were equal
# from http://stackoverflow.com/questions/11148429/r-convert-asymmetric-list-to-matrix-number-of-elements-in-each-sub-list-diffe
max.len <- max(sapply(out, length))
corrected.list <- lapply(out, function(x) {c(x, rep(NA, max.len - length(x)))})
mat <- do.call(rbind, corrected.list)
mat <- data.frame(t(mat))
names(mat) <- paste0('s',3:18)
mat <- gather(mat, run, value)

ggplot(mat, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

The dotted black line is the original point estimate of the variance, and the facets represent different numbers of subjects (s3 being groups of three subjects, s4 being four, etc.).

And the alternative way:

ngroups <- 3:18
reps <- 500
out2<- matrix(NA, length(ngroups), reps)

for (ii in 1:length(ngroups)) {
    for(j in 1:reps) {
        sls <- sleepstudy[sleepstudy$Subject %in% sample(unique(sleepstudy$Subject), ngroups[i], replace = FALSE),]
        out2[i,j] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
    }
}
out2 <- data.frame(t(out2))
names(out2) <- paste0('s',3:18)
out2 <- gather(out2, run, value)

ggplot(out2, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

It appears (for this example, anyway) that the variance doesn't really stabilize until there are at least 14 subjects, if not later.

Answer 4

You could also use a Bayesian mixed model - in that case the uncertainty in the estimation of the random effects is fully taken care of in the calculation of the 95% prediction credible intervals. The new R package brms and function brm, for example, allow for a very easy transition from an lme4 frequentist mixed model to a Bayesian one, as it has near identical syntax.

Answer 5

Angrist and Pischke's "Mostly Harmless Econometrics" has a section titled, "Fewer than 42 clusters", in which they semi-jokingly say,

Therefore, following the... dictum that the answer to life, the universe and everything is 42, we believe the question is: How many clusters are enough for reliable inference using the standard cluster adjustment [akin to the variance estimator in GEE]?

The way my econometrics instructor used to answer questions like yours is, "America is a free country, you can do whatever you like. But if you want to have your paper published, you need to be able to defend what you've done." In other words, you will likely be able to run R or Stata or HLM or Mplus or SAS PROC GLIMMIX code with 6 subjects (and switch to these alternative packages if the one of your choice does not run this), but you will likely have very difficult time defending this approach and justifying asymptotic tests.

I believe that by default, including a variable as a random slope implies including that as a fixed effect, as well, and you need to jump through a lot of syntax hoops if you only want to have this as a random effect with the mean of zero. That's a sensible choice that the software developers have made for you.

Answer 6

I would not use a random effects model with only 6 levels. Models using a 6-level random effect can sometime be run using many statistical programs and sometimes give unbiased estimates, but:

1. I think there is an arbitrary consensus in the statistical community that 10-20 is the minimum number. If you want to have your research published, you'll be advised to look for a journal without statistical review (or be able to justify your decision using fairly sophisticated language).
2. With so few clusters, the between cluster variance is likely to be poorly estimated. Poor estimation of between cluster variance usually translates into poor estimation of the standard error of the coefficients of interest. (random effects models rely on the number of clusters theoretically going to infinity).
3. Often the models simply don't converge. Have you tried running your model? I would surprised with only 12-16 measures per subject that the models converge. When I've managed to get this sort of model to converge I've had hundreds of measurements per cluster.

This issue is addressed in most standard textbooks in the field and you've sort of addressed them in your question. I don't think I'm giving you any new information.

Answer 7

It has been a long time since the original question but I thought I might add a few points pertinent to model selection.

1 - As long as the model is identified (i.e. you have degrees of freedom in the parameter space) you should be able to TRY to fit the model. Depending on the optimization method model may or may not converge. In any case I would not try to include more than 1 or 2 random effects and definitely not more than 1 cross level interaction. In the specific case of the problem presented here if we suspect an interaction between lizard specific characteristics (e.g. age, size, etc.) and treatment/measure characteristics group size 6 may not be enough to make precise enough estimates.

2 - As a couple of answers mention, convergence may be an issue. However my experience is that while social science data has huge convergence problem due to measurement issues, life sciences and especially bio-chemical repeated measures have much smaller standard errors. It all depends on the data generating process. In social and economic data we have to work at various levels of abstraction. In biological and chemical and most certainly astronomical data measurement error is a less of a problem.