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.