Dear nlme users,

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 (S1, ..., S8) and 2 experimental conditions (C0, C1). 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 response 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:

lm0 <- lme( response ~ time*condition, random = ~time | subject/condition, data = psdData)
(time being coded as an integer and condition and subject as factors)

this model would allow a random intercept and slope for each subject and an additional random variation in slope and intercept within each subject for each level of condition. The idea here was to do something similar to the Machines example from Pinheiro and Bates book (pp. 23-24), since I have replications it seemed reasonable to add the "condition within subject" random effect (pp. 27).

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, the reason for this was that as the responses appear linear, for a given trial the deviation from the model at time t depends on the deviation from time t-1.

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 = ~ stimDuration))
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.

Greeting,

Adrien Combaz

Katholieke Universiteit Leuven
Laboratorium voor Neuro- en Psychofysiologie
Campus Gasthuisberg, O&N2 Herestraat 49, bus 1021
3000 Leuven, Belgium
Tel.: +32 16 330428
Fax: +32 16 330734
E-mail: Adrien.Combaz@med.kuleuven.be<mailto:Adrien.Combaz@med.kuleuven.be>


	[[alternative HTML version deleted]]

