[R-sig-ME] repeated measures lme question

[ONKELINX, Thierry]

Dear Melissa,

Cc'ing back to r-sig-mixed-models. Please keep the mailing list in cc.

Pinheiro & Bates (2004) and Zuur et al (2009) use LRT's to test the model. Zuur et al (2009) is focused on ecology which might be more relevant for your analysis.

Have a look at this post of Douglas Bates on p-values https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx op inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey


-----Oorspronkelijk bericht-----
Van: Melissa Dawes [mailto:melissa.dawes op slf.ch]
Verzonden: dinsdag 13 augustus 2013 16:01
Aan: ONKELINX, Thierry
Onderwerp: Re: [R-sig-ME] repeated measures lme question

Dear Thierry,

Many thanks for your reply and suggestion. I didn't know it was possible to keep rw from the pre-treatment years as individual observations in this way, and this approach is clearly preferable to using a mean value.
I tried your example and also applied it to my dataset, but I am still a bit confused about the model. If I ask for anova results in your example, denDF seem too high, particularly in m0 where CO2 is given 1385 denDF. Shouldn't CO2 in this model have <20 denDF in order to account for the random effects structure and whole tree as the unit of replication? If it is indeed correct this way, do you suggest using likelihood ratio tests rather than p values from anova tests (e.g. if LR test indicated the interaction is n.s. the model with CO2 + yearCat could be tested against one with only yearCat)?

Thanks in advance for your additional help. I would also be grateful for any suggestions of books / scripts that deal with this specific type of problem, as I often encounter similar datasets in other aspects of my experiment.

Kind regards,
Melissa



On 13/08/2013 11:10, ONKELINX, Thierry wrote:
> Dear Melissa,
>
> I would go for a different model with a different dataset. Instead of lumping the pre-treatment data I would keep them as individual observations, all having the 'A' treatment.
>
> library(nlme)
> set.seed(12345)
> #create a toy dataset
> dataset<- expand.grid(Year = 1995:2009, ID = 1:20, Disc = 1:5)
> dataset$CO2<- factor(ifelse(dataset$Year>= 2001&  dataset$ID>= 11,
> "E", "A"))
> dataset$rw<- rnorm(nrow(dataset), mean = c(2, 4)[dataset$CO2])
> dataset$yearCat<- factor(dataset$Year)
>
> #create a new variable to define the interaction.
> dataset$yearCatCO2<- ifelse(dataset$CO2 == "A", "A",
> paste(dataset$yearCat, "E", sep = ""))
>
> #model with interaction between year and treatment
> m1<- lme(rw ~ yearCat + yearCatCO2, random = ~1|ID/Disc, data =
> dataset, correlation = corAR1(form = ~ Year), method = "ML") #model
> with only main effects of year and treatment
> m0<- lme(rw ~ yearCat + CO2, random = ~1|ID/Disc, data = dataset,
> correlation = corAR1(form = ~ Year), method = "ML") #a likelihoodratio
> test for the interaction between treatment and year anova(m1, m0)
>
>
> Best regards,
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and Forest team Biometrie&  Kwaliteitszorg / team Biometrics&  Quality
> Assurance Kliniekstraat 25
> 1070 Anderlecht
> Belgium
> + 32 2 525 02 51
> + 32 54 43 61 85
> Thierry.Onkelinx op inbo.be
> www.inbo.be
>
> To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
> ~ Sir Ronald Aylmer Fisher
>
> The plural of anecdote is not data.
> ~ Roger Brinner
>
> The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
>
>
> -----Oorspronkelijk bericht-----
> Van: r-sig-mixed-models-bounces op r-project.org
> [mailto:r-sig-mixed-models-bounces op r-project.org] Namens Melissa Dawes
> Verzonden: maandag 12 augustus 2013 10:54
> Aan: r-sig-mixed-models op r-project.org
> Onderwerp: [R-sig-ME] repeated measures lme question
>
> Dear All,
>
> I am trying to fit a repeated-measures linear mixed effects model
> using lme and have some doubts about the ANOVA results. My experiment
> consists of 20 trees, which were assigned to ambient (A) or elevated
> (E) CO2 for
> 9 years. The trees were then harvested and ring width was measured on multiple stem discs per tree, including several pre-treatment years.
> Therefore, my dataset has the following variables:
>
> CO2 (A or E treatment)
> ID (code for individual tree, 20 with E and 20 with A CO2) yearCat
> (year of ring width, 2001-2009; set as categorical variable) Disc
> (individual disc with ring width measurements, a-c except for one tree
> with 4 discs a-d) rw (ring width measurement for individual year
> within disc within tree)
> cov95to00 (mean ringwidth over pre-treatment years 1995 to 2000 for an
> individual disc)
>
> I would like to analyze the effect of CO2 enrichment and whether it varied during the 9 years of treatment, using the pre-treatment ring width measured on individual discs as a covariable. So far, I have set up the following model, including corAR1 to account for auto-correlation of residuals between years:
>
> rw.model<- lme(rw ~ cov95to00 + CO2 * yearCat, random = ~ 1 | ID/Disc,
> data = ringsyr, na.action = na.omit, method = "REML", corAR1())
>
> However, the ANOVA output shows that year and CO2:yearCat are tested at the individual disc level, which means a very large denom DF, and I am not convinced this is correct.
>
> numDF denDFF-value p-value
> (Intercept)14128.025640.0048
> cov95to001403.331490.0754
> CO21182.705690.1173
> yearCat8412 103.19579<.0001
> CO2:yearCat84123.375340.0009
>
> I also get huge confidence intervals (using the function intervals()) for the random effect associated with Disc or even the message "cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance", which makes me think the model structure might be incorrect.
>
> While I am more confident in a model where measurements are averaged at the tree level (mean of all discs per tree for each year), I think a model with individual discs is better, especially because then I can use pre-treatment measurements that are specific to individual discs rather than averaged at the tree level. I get rather different ANOVA results when I average at the tree level.
>
> I would be grateful for any suggestions about how to correctly apply lme in this situation.
>
> Thank you,
> Melissa
>
>
>          [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> * * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * *
> * Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
> The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.

--
----
Melissa Dawes, PhD
Mountain Ecosystems Research Group
WSL Institute for Snow and Avalanche Research SLF Flüelastrasse 11, CH-7260 Davos Dorf, Switzerland
Tel: +41 (0)81 417 0271
email: melissa.dawes op slf.ch
web: http://www.wsl.ch/info/mitarbeitende/martinm/index_EN

* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.