I am just dipping my toe into the ocean that is linear effects models and am working through Barr et al.'s 'Keeping it Maximal' paper, trying to figure out the best way to fit a lmem for my experiment. Say you have three groups given three different types of drug over three days: 100mg on first day, 50mg on second day, 10mg on last day. The outcome measure is how they feel that next day on some scale (e.g. mood), before they are given their daily dose (i.e. so we are measuring the effects of the previous day's dose). However participants don't come in at exactly the same time each day, thus as each time of measurement the drug will have had less time to take effect.
I would like to know how best to include that random effect of 'time elapsed since dose' into this model, and just how best to fit the model really.
This is a toy dataset. I have not built any trends into it.
dose100mg <- c(6,2,9,4,6,5,2,4,6,7,3,2)
dose50mg <- c(1,2,4,3,6,1,3,3,2,1,4,1)
dose10mg <- c(8,9,7,9,6,7,8,9,8,7,1,3)
timeD1 <- c(24.2,20.5,26,30,22,26,19,23,29,30,24,16)
timeD2 <- c(24,16,28,20,19,28,30,20,18,15,27,32)
timeD3 <- c(21,28,29,30,29,17,23,18,24,16,28,21)
subject <- c(1,2,3,4,5,6,7,8,9,10,11,12)
group <- factor(c(0,1,0,2,1,2,0,2,1,2,1,0))
df <- data.frame(subject, group, dose100mg, dose50mg, dose10mg, cov)
Turn it from wide to long
require(tidyr)
df <- gather(df, dose, score, dose100mg:dose50mg:dose10mg)
Now add the 'hours elapsed since last dose' variable to the dataframe (btw: if anyone knows how to build this into the gather function above I'd appreciate it)
df$hrsElapsed <- c(timeD1, timeD2, timeD3)
Now fit a model. First group*dose plus with random intercepts for subject.
require(lme4)
# random intercepts
anDf_randomintercepts <- lmer(score ~ group*dose + (1|subject), data = df)
anova(anDf_randomintercepts)
Next random slopes, and my first question. Is it better to include hrsElasped as a covariate, like this?
anDf_randomSlopes <- lmer(score ~ group*dose + hrsElapsed + (1|subject) + (1+hrsElapsed|subject), data = df)
Or to include it as a random effect? like this
nDf_randomSlopes <- lmer(score ~ group*dose + (1|subject) + (1+hrsElapsed|subject), data = df)
I know it's not the latter because I get an error message.
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0450795 (tol = 0.002, component 1)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
But I don't know WHY this doesn't work. I would have thought time elapsed would be exactly the sort of variable you'd want to assign to random effects.
What am I doing wrong?
An ancillary question pertains to fitting random slopes for the group-by-subject effect
anDf_randomSlopes <- lmer(score ~ group*dose + (1|subject) + (1+group|subject), data = df)
When I run this i get the error message
Error: number of observations (=36) <= number of random effects (=36) for term (1 + group | subject); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
Why doesn't this work? What does it mean?
Best Answer
I the following model might be good, but the toy data set may be too small for this kind of model. There is a general (linear) effect of hrsElapsed, and this effect is allowed to be different between individuals.
However, I would simplify it to this:
...which I believe is exactly the same model. In the first model, the random effect of subject is included twice.
The model below probably fails to converge because you include the slope of hrsElapsed as a random effect but not as a fixed effect:
I think there must be a fixed effect whose slope can vary between individuals.
As for the model below, the problem is that you use group|subject as a random slope, which gives 36 different combinations (1/group 0, subject 1/group 1, subject 1/group 2, subject 2/group 0 etc.) plus 12 random intercepts for the subjects. The model will not run because you need to have more observations than random effects to be estimated even though only 12 of the 36 random slope combinations are present (subject 1 always belongs to group 0, subject 2 always belongs to group 1 etc).