[R-sig-ME] Modelling random effects for only part of the observations (in lme4)
Hufthammer, Karl Ove
karl.ove.hufthammer at helse-bergen.no
Fri Jan 16 09:36:24 CET 2015
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
More information about the R-sig-mixed-models
mailing list