Suppose I have a dataset where each individual is measured at three timepoints on a variable (we'll call it disease status). What I want to do is estimate the relationship between time and disease. Obviously a mixed model is appropriate, where each person's individual regression line is modeled then aggregated. However, I know that the relationship across timepoints is curvilinear; people tend to increase in disease activity until treated, then disease activity attenuates. Because of this, I decided to model a polynomial term. I am interested in the average shape of the function, aggregated across individuals. On the individual level, each person's regression fits perfectly (3df and 3 parameters estimated), which is usually a problem. However, across individuals it will not fit perfectly (because each person has a different trajectory). Does anyone see a problem with doing this?

Below is some R code to simulate the basic approach I am taking:

    #### make it reproducible

    #### preallocate matrix
data = data.frame(matrix(nrow=50*3, ncol=2))
names(data) = c("disease", "timepoint")
data$ID = 1:50; data= data[order(data$ID),]

    #### loop through each person and give them a disease score/timepoint
for (i in unique(data$ID)){

    #### create varying timepoints
    data$timepoint[data$ID==i] = sort(runif(3, 1, 6))

    #### create the polynomial function with fixed effects
    data$disease[data$ID==i] = data$timepoint[data$ID==i]*rnorm(1,.5,.5) + data$timepoint[data$ID==i]^2*rnorm(1,.25,.5)


        #### model it
mod = lmer(disease~timepoint + I(timepoint^2) + (timepoint + I(timepoint^2) | ID), data=data)
mod2 = lmer(disease~timepoint + (timepoint| ID), data=data)

        ##### plot the data
x = seq(from=min(data$timepoint), to=max(data$timepoint), length.out=50)
y = fixef(mod)[1] + x*fixef(mod)[2] + x^2*fixef(mod)[3]
plot(data$timepoint, data$disease)

        #### see if polynomial terms are significant
anova(mod, mod2)

        #### look at effects
  • 1
    $\begingroup$ I think the residual error and the rest of the variance components will be jointly unidentifiable. $\endgroup$ – Ben Bolker Sep 25 '13 at 13:33
  • $\begingroup$ In my previous example, I didn't have a "person's" effect, and before the variance estimates were all zero. Now, however, they show as non-zero. Is that what you meant by "jointly identifiable?" Or is there something else I am missing? $\endgroup$ – dfife Sep 25 '13 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.