Assessing linearity in a mixed effects model

Question

I have developmental data collected across several grades (1-6), where each child in each grade is measured many times. I would like to be able to assess whether there are any linear or non-linear trends in the response variable across grade. Does it make sense to run a first lmer treating grade as continuous, obtain the residuals, then run a second lmer treating grade as a factor? That is:

fit1 = lmer(
   formula = response ~ (1|child)+grade_as_numeric
   , data = my_data
   , family = gaussian
)
my_data$resid = residuals(fit1)
fit2 = lmer(
   formula = resid ~ (1|child)+grade_as_factor
   , data = my_data
   , family = gaussian
)

As I understand it, fit1 will tell me if there are any linear trends in the data, while fit2 will tell me if there are any non-linear trends in the data in addition to the linear trends obtained in fit1.

If this is sensible, how might I apply it to a second binomial response variable given that the residuals from a binomial model are not 0/1?

Answer 1

Not really an answer, but I was interested in trying it out... I assume that the pattern is not easily recognisable just by plotting? So I tried to make up some data that might behave this way:

set.seed(69)
id<- rep(1:20, each=6)
x<-rep(1:6, 20)
y<-jitter(x+id/5, factor=5) + jitter(sin(x), factor=5)
df1<-data.frame(id, x, y)

plot(y~x)
xyplot(y~x|id, data=df1, type="l")

For me, if I had this data without knowing how it was made, I think I would have trouble picking out the overlaid signal and perhaps assume it was linear. Resid vs fitted of lmer1 (below) doesn't show much, but the resid vs x is more suggestive.

lmer1<-lmer(y~x+(1|id), data=df1, REML=F)
xyplot(resid(lmer1)~fitted(lmer1), type=c("p", "smooth"))
xyplot(resid(lmer1)~x, type=c("p", "smooth"))

Using your suggestion of conducting an ANOVA on the residuals gives a significant effect, indicating perhaps some kind of systematic difference:

lmer2<-lmer(resid(lmer1)~factor(x)+(1|id), data=df1, REML=F) #sd attributed to id is 0
anova(lmer2)

So perhaps this method may be useful to determine whether you need to include higher order terms by using, maybe, increasing order polynomials:

lmer3.1 <- lmer(y~poly(x,2)+(1|id), data=df1, REML=F)
lmer3.2 <- lmer(y~poly(x,3)+(1|id), data=df1, REML=F)
lmer3.3 <- lmer(y~poly(x,4)+(1|id), data=df1, REML=F)
anova(lmer1, lmer3.1, lmer3.2, lmer3.3)

In this method the cubic function 'wins' and might be a useful approximation of what's going on and comes quite close to the generated model:

lmer4<-lmer(y~x+sin(x)+(1|id), data=df1)
anova(lmer3.2, lmer4)

I don't really know if this helps or not, but hopefully this simulated data mirrors your problem and somebody can use this to give a more exact answer to your question. I'm not sure about the binomial part of your question.