Random effect with zero SD in LMM

In my mixed-effects model there are one fixed effect and two random effects (subject and time of measurement).

fit <- lmer(DV ~ group + (1|subject) + (1|time))

Now for some DVs variance and standard deviation for random effect time are 0. On pages 10/11 of his book (lme4: Mixed-effects modeling with R) Douglas Bates wrote:

An estimate of 0 for σ does not mean that there is no variation between the groups. (...) The estimate, σ1 = 0, simply indicates that the level of “between-group” variability is not sufficient to warrant incorporating random effects in the model. The important point to take away from this example is that we must allow for the estimates of variance components to be zero. We describe such a model as being degenerate, in the sense that it corresponds to a linear model in which we have removed the random effects associated with Batch. Degenerate models can and do occur in practice. Even when the final fitted model is not degenerate, we must allow for such models when determining the parameter estimates through numerical optimization.

So I do not entirely understand what that means for me. On the hand it sounds like it's not justified to include time as a random effect ("..not sufficient to warrant incorporating random effects in the model"), on the other hand it sounds as if it's OK for variance/SD to be zero (as long as interpreted correctly).

Can you help me with the correct interpretation of the paragraph?

This issue was recently discussed on the r-sig-mixed-models mailing list: See this initial posting and the replies (including one by Douglas Bates).

See also the r-sig-mixed-models FAQ -- in particular, the section: Why is my random effect variance estimated as zero, or correlations estimated as +/- 1? What should I do about it?

But as a more direct answer, consider a simple example: $n$ subjects measured $p$ times and we want to account for the multilevel structure of the data by adding a random intercept at the subject level.

library(nlme)

n <- 10 ### number of subjects
p <- 5  ### number of measurements per subject

### subject id variabble
subj <- rep(1:n, each=p)

### simulate some data (with very little subject-level variance)
set.seed(137814)
yi <- 5 + rep(rnorm(n,0,.01), each=p) + rnorm(n*p,0,1)

### means of the subjects
aggregate(yi, by=list(subj), FUN=mean)

This yields:

Group.1        x
1        1 4.573781
2        2 5.277764
3        3 5.147564
4        4 4.848544
5        5 4.784689
6        6 4.439724
7        7 5.132603
8        8 4.522901
9        9 4.692272
10      10 5.330356

Note that they are not all the same. So, it appears as if there is variability at the subject level. Let's fit the multilevel model:

### fit multilevel model
res <- lme(yi ~ 1, random = ~ 1 | subj)
summary(res)

The results:

Linear mixed-effects model fit by REML
Data: NULL
AIC      BIC    logLik
153.3184 158.9939 -73.65922

Random effects:
Formula: ~1 | subj
(Intercept) Residual
StdDev: 3.568663e-05 1.045392

Fixed effects: yi ~ 1
Value Std.Error DF  t-value p-value
(Intercept) 4.87502 0.1478408 40 32.97479       0

Standardized Within-Group Residuals:
Min         Q1        Med         Q3        Max
-2.5476373 -0.6569860  0.2263918  0.7942828  2.0033502

Number of Observations: 50
Number of Groups: 10

Note that the subject-level variance is estimated to be essentially zero. In essence, this means that the variability in the means at the subject level is not really greater than what we would expect based on the estimated amount of error variance for these data. So, going back to your quote: "An estimate of 0 for $\sigma$ does not mean that there is no variation between the groups." But it suggests for the data at hand that adding this variance component doesn't improve the fit of the model.

In fact, let's try this by just leaving out the random effect at the subject level (so we end up with just a very simple regression model with an intercept):

res <- lm(yi ~ 1)
summary(res)

The results:

Call:
lm(formula = yi ~ 1)

Residuals:
Min      1Q  Median      3Q     Max
-2.6633 -0.6868  0.2367  0.8303  2.0943

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   4.8750     0.1478   32.98   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.045 on 49 degrees of freedom

Essentially the same results as obtained earlier based on the multilevel model. So, you can fit the multilevel model (nothing wrong with that), but in essence, if a particular variance component is estimated to be (essentially) zero, the results are practically the same as if you had left out that random effect.

• Thanks for the great answer! So I'll leave in the random effect (for the sake of consistency across DVs). – user45065 Aug 19 '14 at 9:58