Effect of baseline on change over time in mixed models?

Question

I was wondering whether and how it is possible to model, within some sample, change of outcome over time which is dependent on the baseline value of this outcome, using a mixed model?

Imagine, for example, a situation where the same knowledge test is administered 5 times to the same group of people. Since the questions are always the same, the students will learn the correct answers over time and will score higher on each administration. However, in those who scored high in the first place, there will be less change than in those who at first scored poorly. Thus, it is quite apparent that the rate of change depends on the baseline value.

I know that in mixed models, I could include a random slope for time in addition to a random intercept, to account for the fact that in some students there will be more change than in others. However, am I right to assume that it is not possible or meaningful to include the value of the first measurement as a baseline covariate (and its interaction with time)? It doesn't "feel" right to me in any case. But on the other hand, it boggles my mind that it wouldn't be possible to explicitly model the effect of the baseline value using a fixed effect. I must admit I'm somewhat confused over this. Any help would be highly appreciated.

Answer 1

This seems to be a growth model scenario. Suppose we had the following variables:

• occasion: taking values 1,2,3,4,5 to reflect the occasion that test was taken, 1 being the first, or baseline.
• ID: the identifier of the each participant.
• score: the test score for this participant on this test occasion.

Random intercepts for ID will take care of the different baselines (subject to having sufficient participants.

Thus a simple linear mixed effects model for these data is (using lme4 syntax):

score ~ occasion + (1|ID)

or

score ~ occasion + (occasion|ID)

where the latter allows the linear slope of occasion to vary among participants

However, for the particular example in the OP, we have the additional problem that the score variable is bounded above by the maximum score on the test. To allow for this, we need to cater for non-linear growth. This could be achieved in a variety of ways, the simplest being the addition of quadratic and possibly cubic terms to the model:

score ~ occasion + I(occasion^2) + I(occasion^3) + (1|ID)

Let's look at a toy example:

require(lme4)
require(ggplot2)

dt2 <- structure(list(occasion = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4), score = c(55.5, 74.5, 92.5, 97.5, 98.5, 54.5, 81.5, 94.5, 97.5, 98.5, 47.5, 68.5, 86.5, 96.5, 98.5, 56.5, 86.5, 91.5, 97.5, 98.5, 60.5, 84.5, 95.5, 97.5, 99.5, 73.5, 87.5, 96.5, 98.5, 99.5), ID = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L), .Label = c("1", "2", "3", "4", "5", "6"), class = "factor")), .Names = c("occasion", "score", "ID"), row.names = c(25L, 26L, 27L, 28L, 29L, 31L, 32L, 33L, 34L, 35L, 37L, 38L, 39L, 40L, 41L, 43L, 44L, 45L, 46L, 47L, 49L, 50L, 51L, 52L, 53L, 55L, 56L, 57L, 58L, 59L), class = "data.frame")

m1 <- lmer(score~occasion+(1|ID),data=dt2)

fun1 <- function(x) fixef(m1)[1] + fixef(m1)[2]*x

ggplot(dt2,aes(x=occasion,y=score, color=ID)) + geom_line(size=0.65) + geom_point() +
 stat_function(fun=fun1, geom="line", size=1, colour="black")

Here we have plots for 6 participants that have been measured over 5 successive occasions, and we have plotted the fixed effects with the solid black line. Clearly this is not a good model for these data, so we introduce a quadratic term, and then a cubic term, after centering the data to reduce collinearity:

dt2$occasion <- dt2$occasion - mean(dt2$occasion)

m2 <- lmer(score~occasion + I(occasion^2) + (1|ID),data=dt2)
fun2 <- function(x) fixef(m2)[1] + fixef(m2)[2]*x + fixef(m2)[3]*(x^2)

m3 <- lmer(score~occasion + I(occasion^2) + I(occasion^3) + (1|ID),data=dt2)
fun3 <- function(x) fixef(m3)[1] + fixef(m3)[2]*x + fixef(m3)[3]*(x^2) + fixef(m3)[4]*(x^3)


p2 <- ggplot(dt2,aes(x=occasion,y=score, color=ID)) + geom_line(size=0.5) + geom_point()

p2 + stat_function(fun=fun2, geom="line", size=1, colour="black")

Here we see that the quadratic model is an obvious improvement over the linear-only model, but is not ideal because it underestimates the scores for the final measurement, and overestimates it for the previous one.

The cubic model on the other hand appears to work very well:

p2 + stat_function(fun=fun3, geom="line", size=1, colour="black")

A slightly more sophisticated approach is to recognise the upper bound explicity, and use (for example) a logistic growth curve model. One way to accomplish this is to transform the outcome to a proportion (of the upper bound), say $\pi$ and then model the logit of this proportion, $\pi/(1-\pi)$ as the outcome of a linear mixed effects model. In addition to recognising the upper bound, this has the added advantage of modelling heteroscasticity in the residuals of the untransformed data, since it seems likely that over successive tests (assuming that the results get better) there will be less variance.

Putting this into practice, as expected, this also models the overall trend in the data very well:

pi <- dt2$score/100
dt2$logitpi <- log(pi/(1-pi))

m0 <- lmer(logitpi~occasion+(1|ID),data=dt2)
funlogis <- function(x) 100*exp(fixef(m0)[1] + fixef(m0)[2]*x)/(1+exp(fixef(m0)[1] + fixef(m0)[2]*x))
p2 + stat_function(fun=funlogis, geom="line", size=0.5, colour="black")

The following shows the cubic mode and the logistic growth models plotted together, and we see very little difference between them, though as mentioned above we may prefer the logistic growth model due to the heteroscedasticity issue:

p2 +  stat_function(fun=fun3, geom="line", size=1, colour="black")  +
stat_function(fun=funlogis, geom="line", size=1, colour="blue")

A more sophisticated approach still would be to use a nonlinear mixed effects model where the logistic growth curve is modelled explicitly, allowing for random variation in the parameters of the logistic function itself.