lme: random effects for replicated growth curves

I am measuring the evolution of the brain response to a visual stimulation over time. The measures are done every seconds from 1 second to 14 seconds (each measure at time t gives a value summarizing the magnitude of the brain response from time 0 to time t). I have 8 subjects and 2 experimental conditions. For each subject and condition I replicate the measurement 12 times. I obtain therefore for each subject and condition 12 growth curves.

My aim is to assess the influence of the experimental condition on the response profile. The responses appear to be fairly linear in time. Additionally, a noticeable aspect of the obtained results is the increase of the response variability over time (for each subject individually as well as between subjects). I am not sure about how to specify properly my random effect for the linear mixed model (using the nlme package).

My initial idea was to start with the following model (time being coded as an integer and condition and subject as factors):lm0 <- lme( response ~ time*condition, random = ~time | subject/condition, data = psdData).

Now due to the increase of variability of the response over time, I also get an increase of variability of the residuals over time. I initially wanted to try correcting for this using correlation structures, unfortunately, due to the replications, it seems that I cannot use the correlation option; when trying lm1 <- update(lm0, corr = corCAR1(form = ~time)) or lm1 <- update(lm0, corr = corCAR1(form = ~time | subject/condition)), I obtain the error:

Error in Initialize.corCAR1(X[[2L]], ...) : covariate must have unique values within groups for "corCAR1" objects


and when trying lm1 <- update(lm0, corr = corCAR1(form = ~time | replicate/subject/condition)), I obtain the error:

Error in lme.formula(response ~ time * condition, data = psdData,  : incompatible formulas for groups in 'random' and 'correlation'


As I didn't manage to use correlation structures, I tried to use variance functions to model heteroscedasticity: lm1 <- update(lm0, weights = varPower(form = ~ time))

But when looking at the residuals, it still displayed an increased variance over time. Same thing happened when using "varExp" or "varConstPower".

I tried then another model, this time considering trials as units of measurements (and not subject as previously) and including the trial to trial variability within each subject. I created the factor trialInSub <- condition:replicate which list all 24 trials of each subject (12 per condition) and wrote the model:

lm0 <- lme( response ~ time*condition, random = ~time | subject/trialInSub, data = psdData)


giving a random intercept and slope for each subject and for each trial within subject.

As this model fits a line for each replicate, I don't have anymore the issue about increasing variance over time. However when looking at the autocorrelation of the models residuals, I would suspect something is wrong with the model: from lag 2 to 7, the autocorrelation decreases from 0.5 to -0.5 and from lag 7 to 10 increases back to 0.04. Additionally I have a clear linear trend when plotting residuals against fitted values of the model for time = 1 (it looks fairly normally distributed for all other times) I tried different correlation structures but it didn't change much my results.

I am really not sure that I handle the problem in a proper way and that the models I tried really reflect the design of the experiment and the structure of the data. I hope someone can direct me towards a solution that would be statistically sound.

Here are a couple thoughts:

First, the reason that you're getting that error message is that the correlation has to apply to the same grouping as your random effects parameters. The way you have it set up now, the random effect is based on subject and condition, but not replicate. So all 12 replicates have the same random effect level, which means that when you specify a correlation structure, R wants that structure to define correlations between 12*14=168 different observations, which is not what you've got in mind. So one option is you could specify a second random effect to pertain to each replicate. There may also be a way to specify that the data is grouped by replicate without requiring a random effect to be associated, possibly with a groupedData object, but I'm not positive how the syntax would go. But that's why you're getting the error message for the corCAR1 command, it needs that grouping to be the same as the random effects groupings.

Also, I'm not entirely familiar with the varPower function in R but when you look at the residuals, are they standardized to their modeled variance level? That is, if the varPower function specifies a larger variance expected for certain observations, are the residuals from those observations normalized to that higher variance? If they're not, that could explain why you're still observing higher variance in the residuals even when you use the varPower structure. It's because the variance has been accounted for in the model fitting, but not adjusted for when you visualize it. In general, I don't believe that specifying a weight vector will cause the residuals to be normalized in R by default, but that could depend on the package and I could be wrong about that.

Third, regarding your observed autocorrelation pattern, that is a pretty standard pattern for something like an AR(1) setup. For a basic AR(1) model, you would expect it to have absolute value decaying exponentially to zero, and possibly oscillating around zero in the process. So I think it would make sense for you to try and use grouped data somehow to accommodate the AR(1) structure the way you've been trying to do. There should be a way to specify that, and I think it involves a groupedData object, but I don't know the details. You can consult http://users.stat.umn.edu/~sandy/courses/8311/handouts/lme.pdf for some guidance.

Good luck, keep us posted.

• Hi Mike, you were right about the VarPower function, I should explicitly mention type="pearson or type=residual. This solves the heteroscedasticity issue but not the correlation issue (I'll develop as an answer below). Concerning the second model and the observed autocorrelation pattern, neither AR(1) nor other structures really helps. It doesn't really decays exponentially to zero, it is more linearly decreasing from 0.5 to -0.5 and then linearly increasing back to 0. I am wondering if this pattern observed in the residuals could point out to the fact that I should use a non-linear model. Aug 20, 2013 at 10:17
• About the correlation issue, I would like to be able to specify that the correlation structure is along time within each trial independently, without having to specify the trials as a random effect of my model. This the same issue as posted here link. It seems that no matter how the groups are defined in the GroupedData object, the groups for the lme object are fully defined by the random effects of the model. I could not find a way to specify a grouping factor for an lme object without it being a random effect. Aug 20, 2013 at 10:30
• Ya I guess the groupedData object is not going to be effective for that, now that I read more closely. I'm trying to experiment with this kind of issue and I actually am able to create lme models with different groupings for the correlation and the random effect. Have you tried running this without using the update statement? Like just including the "corr=" statement in the initial model command? This has been bugging me for a while now and I had a similar problem earlier this summer. Aug 21, 2013 at 21:11
• Apparently, the problem is not the update statement, but it seems that the group specified in the correlation must be nested within the groups mentioned as random effect. Having lm0 <- lme( response ~ time*condition, random = ~time | subject/condition, data = psdData), I was able to run lm1<-update(lm0, corr = corCAR1(form=~time | subject/condition/trial)). Though i'm not sure how well it works, my residuals still seem correlated, additionally fm1.lme\$apVar gives "Non-positive definite approximate variance-covariance" and phi is estimate to 1 for the CorCAR1. I'll keep on digging... Aug 22, 2013 at 12:57