[R-sig-ME] Modelling random effects for only part of the observations (in lme4)

[Hufthammer, Karl Ove]

Jarrod Hadfield  wrote:
> Create a binary vector (lets call it b) with 0 for observations where
> you don't want random effects and a 1 where you do. Then fit:
> 
> (b|group)

Thank you for the suggestion. I have now tried it, but can't get it to work. I do get a variety of error messages, depending on how I create the 'group' variable, though. :)

These are the errors/warnings: Either (with group = gr1 in code below)

─────────────
Error: number of observations (=100) <= number of random effects (=134) for term (b | gr); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
─────────────

or (with group = gr2 in code below)

─────────────
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?
─────────────

or (with group = gr3 in code below)

─────────────
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels
─────────────

Here's some example code. The simulation is based on an assumption that belonging to a (treatment) group *adds* a random effect, i.e., it adds to the total subject error (so that the *residual* error is identical in both arms). I guess it could be argued that a model where the total *subject error* was identical in the two arms, and belonging to a group would only add intragroup correlation, would also be reasonable. I'm not sure which model is most realistic.


mu = 10                       # Average level in control group
treat = 5                     # Treatment effect (additive)
n = 100                       # Number of individuals
prop.treat = .5               # Treatment proportion
n.treat = round(n*prop.treat) # Number in control group
n.contr = n - n.treat         # Number in treatment group
grsize = 3                    # Group size (fixed, but could be random)
sigeps = 3                    # Residual standard deviation
siggrp = 1                    # Group standard deviation

# Set up treatment factor
arm = factor(c(rep("treatment", n.treat), rep("control", n.contr)))

# Create observations
y = rnorm(n, mu, sd = sigeps)
y[arm=="treatment"] = y[arm=="treatment"] + treat +
                      rep(rnorm(ceiling(n.treat/3), 0, siggrp),
                          each=grsize)[1:n.treat]

# Create group (with control subjects each having a unique group)
gr = c(rep(1:ceiling(n.treat/3), each=3)[1:n.treat], (n.treat+1):n)

# Create group (with control subjects sharing a group)
gr2 = gr
gr2[(n.treat+1):n] = 0

# Create group (with control subjects having group value NA)
gr3=gr2
gr3[gr3==0]=NA

# Convert group variables to factors
gr = factor(gr)
gr2 = factor(gr2)
gr3 = factor(gr3)

# Correct (?) random effects model
b = as.numeric(arm) - 1 # Equal to 1 for random effect and 0 for no random effect
table(b, arm)

# Try to fit lme4 models
library(lme4)
lmer(y ~ arm + (b|gr))
lmer(y ~ arm + (b|gr2))
lmer(y ~ arm + (b|gr3))

-- 
Karl Ove Hufthammer