Search code examples
rlme4model-fitting

Error in fitting mixed model: number of observations <= number of random effects


I am trying fitting several mixed linear models.

 summary(model3 <- lme4::lmer(LPP2POz ~ COND + (COND|ID), data = dataLPP2POz))

whose output is

Error: number of observations (=75) <= number of random effects (=75) for term (COND | ID); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable

Can anyone can figure out why this model turns back this error?

Here's the dataset:

> dput(head(dataLPP2POz))
structure(list(ID = structure(c(1L, 1L, 1L, 2L, 2L, 2L), .Label = c("01", 
"04", "06", "07", "08", "09", "10", "11", "12", "13", "15", "16", 
"17", "18", "19", "21", "22", "23", "25", "27", "28", "30", "44", 
"46", "49"), class = "factor"), GR = c("RP", "RP", "RP", "RP", 
"RP", "RP"), SES = c("V", "V", "V", "V", "V", "V"), COND = structure(c(1L, 
2L, 3L, 1L, 2L, 3L), .Label = c("NEG-CTR", "NEG-NOC", "NEU-NOC"
), class = "factor"), LPP2POz = c(7.91468942320841, 9.94838815736199, 
10.2186482048953, 1.07455889922813, 1.65917850515029, 3.22422743232682
)), row.names = c(NA, 6L), class = "data.frame")

Solution

  • The specification COND + (COND|ID) indicates that there is a fixed effect of condition and a random effect of condition across individuals; that is, every individual/condition combination gets a random-effect value. Approximately

    y_ij = alpha_i + epsilon_{c,ij} + epsilon_{r,ij}
    

    (although that's not quite how lmer parameterizes it); alpha_i is the fixed effect of condition i, epsilon_{c,ij} is the random effect for condition i for individual j, and epsilon_{r,ij} is the residual variation for condition i/individual j. The problem is that because every individual is only measured once in every condition, the two epsilon terms are confounded (jointly unidentifiable).

    • As @Bloxx says, you can force lmer to fit the model anyway by specifying control = lmerControl(check.nobs.vs.nRE = "ignore") in your lmer() call. If you do this, though, you can't interpret the random effects covariance matrix unless you are an expert/really know what you're doing (the fixed effects should still be OK).
    • If you wanted to fit the same model in glmmTMB you could set dispformula = ~0 to zero out the residual variance (actually it gets set to a very small but non-zero value ...)
    • For some kinds of simpler models (but not this one), as in this example, you could simplify the random effects model so that it no longer clashed with the residual variance term. For example, if your model was (1|ID/COND) instead (i.e. a compound symmetric model/intercept variation of ID and COND within ID, equivalent to (1|ID) + (1|ID:COND)) you could solve the problem by dropping the nested term (which is exactly the same as the residual variance term) and simplify the model to (1|ID).