This seems to be a growth model scenario. Suppose we had the following variables:
occasion: taking values
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
score ~ occasion + (1|ID)
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:
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) + fixef(m1)*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) + fixef(m2)*x + fixef(m2)*(x^2)
m3 <- lmer(score~occasion + I(occasion^2) + I(occasion^3) + (1|ID),data=dt2)
fun3 <- function(x) fixef(m3) + fixef(m3)*x + fixef(m3)*(x^2) + fixef(m3)*(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) + fixef(m0)*x)/(1+exp(fixef(m0) + fixef(m0)*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.