lme: random effects for replicated growth curves
Hi Ben, Thanks for your suggestion, I tried to play around with the gls function, but didn't manage get an adequate model so far. I am looking now at alternatives to mixed-effect models, such as growth mixture models. I'll keep you posted if I find a proper solution. -----Original Message----- From: Ben Bolker [mailto:bbolker at gmail.com] Sent: Tuesday, August 20, 2013 4:54 PM To: Adrien Combaz; r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] lme: random effects for replicated growth curves [cc'ing back to r-sig-mixed-models ] I'm still not 100% sure what the problem is, and I don't know if I have time to dig through in detail, but here's an example that _might_ help -- if you want to use a grouped correlation structure without incorporating random effects ... d <- data.frame(x=1:1000, f=rep(1:50,each=20), y=rnorm(1000),time=rep(1:20,50)) lme(y~x,correlation=corAR1(form=~time|f),random=~1|f,data=d) gls(y~x,correlation=corAR1(form=~time|f),data=d) ## without RE
On 13-08-20 07:48 AM, Adrien Combaz wrote:
Hi Ben, I realize that my question was quite long and I appreciate that you took the time to go through it and answer it. You were right about my problem with the variance function, when observing the residuals of the new model, I should explicitly mentioned type="pearson" or type="normalized" (both giving the same result). I should have seen it on Pinheiro and Bates book (pp. 218). I had a slightly outdated version of the nlme package where the help of the residuals.lme function mentioned that the "pearson" was the default type. This solved my heteroscedasticity issue. However I still have a correlation problem. My model is now lm0 <- lme( values ~ time*condition, data = dataGroup, random = ~time|subject/condition, weights = varPower(form = ~time), control = list(opt="optim"), method = "REML") The model gives a straight line along time for each subject and condition that is supposed to represent all 12 trials for this subject/condition. Now it appears that a trial deviating from the model at time t will keep on deviating at time t+1 (this is related the increase of variability in time that I mention in my initial question). The variance function mentioned earlier corrects the residuals in a way to keep this variation constant along time, but it still have to account somehow for the correlation along time of the residuals within each trials. What I would like, is to be able to specify that the correlation structure is along time within each trial independently, without having to specifity the trials as a random effect of my model. This the same issue as posted here https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q1/003311.html, unfortunately no answer was provided. So I am still stuck - for the moment - with my other model, which specify a random slope and intercept for each subject and for each trial within subject. lm1 <- lme( values ~ time*condition, data = dataGroup, random = ~time|subject/trialInSub, control = list(opt="optim"), method = "REML"). Unfortunately, as mentioned in the original question, this model leads also to non-normal and correlated residuals. I could not correct for this using variance and/or correlation structures. I am not sure I can do much more with this type of model, and I am actually wondering if those patterns observed in the residuals could point out to the fact that I should use a non-linear model. Adrien Combaz -----Original Message----- From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Ben Bolker Sent: Monday, August 19, 2013 9:39 PM To: r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] lme: random effects for replicated growth curves Adrien Combaz <Adrien.Combaz at ...> writes:
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.
Part of the reason this question may not have gotten the attention it deserves is that it's rather long. I know I often have the reaction "hmm, that looks complicated, I don't have time to look at it right now" to long/involved questions. It is certainly sensible that you want to explain the context of your problem and what you've already tried, but it may be that some problems are just a bit too involved to be easily tractable in a mailing list/forum format. [snip] An increase in variability over time would seem to be more sensibly modeled by your second attempt (varStruct/heteroscedasticity) than by your first (temporal autocorrelation); sometimes it can also be handled by transformation (although that is likely to affect/mess up linearity to some extent at the same time)
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".
Did you make sure to use residuals(.,type="pearson") ? Otherwise, the heteroscedasticity may be correctly modeled but the (raw) residuals won't reflect that the model has taken the heteroscedasticity into account (see also type="normalized", for models with correlation structure: more generally, ?residuals.lme Good luck Ben Bolker
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models