One observation per level in mixed-effect model

Question

Field explains how to analyse repeated-measures data using linear mixed-effect models (LME). See Field et al., Discovering Statistics Using R, 2012, p. 573.

However, the way he specifies the model, there appears to be only one observation per level of the grouping factors. Is this a mistake in the textbook? If not, why not? It seems to me the random effects specify a full model and fit the data (almost) exactly.

The code is as follows:

library(reshape2)
library(nlme)
# Load dataset:
dat.wide <- read.delim("http://www.sagepub.com/dsur/study/DSUR%20Data%20Files/Chapter%2013/Bushtucker.dat")

dat.wide
#   participant stick_insect kangaroo_testicle fish_eye witchetty_grub
# 1          P1            8                 7        1              6
# 2          P2            9                 5        2              5
# 3          P3            6                 2        3              8
# 4          P4            5                 3        1              9
# 5          P5            8                 4        5              8
# 6          P6            7                 5        6              7
# 7          P7           10                 2        7              2
# 8          P8           12                 6        8              1


dat <- melt(dat.wide, variable.name="animal", value.name="retch")
head(dat)
#   participant       animal retch
# 1          P1 stick_insect     8
# 2          P2 stick_insect     9
# ...

# set contrasts, not relevant to question but keeps example same as in book:
PartvsWhole <- c(1, -1, -1, 1)
TesticlevsEye <- c(0, -1, 1, 0)
StickvsGrub <- c(-1, 0, 0, 1)
contrasts(dat$animal) <- cbind(PartvsWhole, TesticlevsEye, StickvsGrub)

# Fit intercept term, then add "animal" term.
# NOTE: random effects are animal nested within participant, as in textbook.
# This would presumably give one observation per group?    
lme1 <- lme(retch ~ 1, random=~1|participant/animal, data=dat, method="ML")
lme2 <- lme(retch ~ 1 + animal, random=~1|participant/animal, data=dat, method="ML")

# I have checked these are the same results as in textbook
anova(lme1, lme2)    
summary(lme2)

# residuals are near-zero (e-05)
resid(lme2)

ran1 <- random.effects(lme1)
ran2 <- random.effects(lme2)
# same number of random effects as observations:
nrow(ran1$animal)
    nrow(ran2$animal)
nrow(dat)

The concern of one observation per level seems to be confirmed by using package lme4 instead:

library(lme4)
# Using lme4 produces an error:
# "Error in checkNlevels(reTrms$flist, n = n, control) : 
# number of levels of each grouping factor must be < number of observations"    
lmer1 <- lmer(retch ~ 1 + (1|participant/animal), data=dat, REML=F)
lmer2 <- lmer(retch ~ 1 + animal + (1|participant/animal), data=dat, REML=F)
anova(lmer1, lmer2)

# see http://stackoverflow.com/questions/19713228/lmer-returning-error-grouping-factor-must-be-number-of-observations
# can force fit with lme4:
# control=lmerControl(check.nobs.vs.nlev = "ignore",
#                     check.nobs.vs.rankZ = "ignore",
#                     check.nobs.vs.nRE="ignore"))

lmer1 <- lmer(retch ~ 1 + (1|participant/animal), data=dat, REML=F,
              control=lmerControl(check.nobs.vs.nlev="ignore",
                                  check.nobs.vs.rankZ="ignore",
                                  check.nobs.vs.nRE="ignore"))
lmer2 <- lmer(retch ~ 1 + animal + (1|participant/animal), data=dat, REML=F,
              control=lmerControl(check.nobs.vs.nlev="ignore",
                                  check.nobs.vs.rankZ="ignore",
                                  check.nobs.vs.nRE="ignore"))
# ignoring errors, we get the same results as with lme:
anova(lmer1, lmer2)
anova(lme1, lme2) # same

Should the random term should be 1|participant instead?

score 4 · Accepted Answer · answered Oct 19 '14 at 15:43

It certainly seems that animal nested within participant is confounded with the residual variance, although it would help to see a relevant descriptive excerpt from the book (those pages aren't available on Google Books).

with(dat,table(participant,animal))

shows the lack of replication that you refer to.

Another hint that something is wrong with the model:

intervals(lme1)
## Error in intervals.lme(lme1) : 
  cannot get confidence intervals on var-cov components: 
  Non-positive definite approximate variance-covariance

It's also worth pointing out (I don't know if the authors do) that the estimated among-participant variance is zero, which could be due to a lack of replication or to a negative correlation within individuals (e.g. someone who likes grubs less than average would like stick insects more than average); this could in principle be handled by a compound symmetry correlation structure with a negative correlation.

I've posted a similar example.

One observation per level in mixed-effect model

1 Answers1

Linked