lme: random effects for replicated growth curves

Question

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.

Answer 1

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.