[R] before-after control-impact analysis with R

[Mark Difford]

Hi Nikolaos,

>> My question again is: Why can't I reproduce the results? When I try a 
>> simple anova without any random factors:

Lack of a "right" result probably has to do with the type of analysis of
variance that is being done. The default in R is to use so-called Type I
tests, for good reason. SAS, I think, still uses a type of test that many
authorities consider to be meaningless, as a general, first-level, test.

There is a long, long thread on this, going back to about approx April/May
1999, when a suitable Ripplyism was coined, which I believe found its way
into the fortunes package. But

RSiteSearch("type III")

should do for a start. Also see

?drop1
library(car)
?Anova

HTH, Mark.


Nikolaos Lampadariou wrote:
> 
> Hello everybody,
> 
> In am trying to analyse a BACI experiment and I really want to do it 
> with R (which I find really exciting).  So, before moving on I though it 
> would be a good idea to repeat some known experiments which are quite 
> similar to my own. I tried to reproduce 2 published examples but without 
> much success. The first one in particular is a published dataset 
> analysed with SAS by McDonald et al. (2000). They also provide the 
> original data as well as the SAS code. I don't know much about SAS and i 
> really want to stick to R. So here follow 3 questions based on these 2 
> papers.
> 
> 
> Paper 1
> (McDonald et al. 2000. Analysis of count data from before-after 
> control-impact studies. J. Agric. Biol. Envir Stat., 5:262-279)
> 
> Data were collected from 1984, 1989, 1990, 1991, 1993, and 1996. Samples 
> from 1984 are considered as before an impact and the remaining years as 
>   after the impact. Each year, 96 transects were sampled (36 in the 
> oiled area and 60 in the non-oiled area; "0" is for oiled and "1" for 
> non-oiled). The authors compare 3 different ways of analysing the data 
> including glmm. The data can be reproduced with the following commands 
> (density is fake numbers but I can provide the original data since I've 
> typed them in anyway):
> 
>  >x<-c(rep("0",36),rep("1",60))
>  >oiled<-rep(x,6)
>  >year<-c(rep(1984,96), rep(1989,96), rep(1990,96), rep(1991,96), 
> rep(1993,96), rep(1996,96))
>  >density<-runif(576, min=0, max=10)
>  >transect<-c(rep(1:96,6))
>  >oil<-data.frame("oiled"=oiled, "transect"=transect, "year"=year, 
> "density"=density)
> 
> Question 1:
> I can reproduce the results of the repeated measures anova with:
>  >oil.aov1<-aov(density~factor(year)*factor(oiled)+Error(factor(transect))
> 
> But why is the following command not working?
>  >oil.aov2<-aov(density~oiled*year + Error(oiled/transect), data=oil)
> 
> After reading the R-help archive, as well as Chambers and Hasties 
> (Statistical Models in S) and Pinheiro's and Bates (Mixed effects models 
> in S and S-plus) I would expect that the correct model is the oil.aov2. 
> As you might see from the data, transect is nested within oiled and I 
> still get the wrong results when I try +Error(transect) or 
> +Error(oiled:transect) and many other variants.
> 
> Question 2 (on the same paper):
> The authors conclude that it is better to analyse the data with a 
> Generalized Linear (Mixed) Model Technique. I tried lme and after 
> reading Douglas Bates article (May 2005, vol. 5/1) also lmer as follows:
>  >oil.lme<-lme(density~year*oiled, random=~1|oiled/transect)
> or
>  >harvest.lmer<-lmer(abund~H*BA+BA/year+BA/year:H+site:BA+(1|H/site))
> but again no luck. I will get always some error messages or some 
> interactions missing.
> 
> 
> Paper 2
> (Keough & Quinn, 2000. Legislative vs. practical protection of an 
> intertidal shoreline in southeastern Australia. Ecol. Appl. 10: 871-881)
> 
> Data were collected from 1989, 1990, 1991, 1993, 1994, 1995, 1996, 
> 1997). the 1989-1991 are the before impact years and the other 4 years 
> are the after the impact years. Eight sites were samples (2 protected 
> and 6 impacted). In this dataset, site is nested within harvest (H) and 
> year is nested within before-after (BA). Also, site is considered as 
> random by the authors. The data (fake again) can be produced with the 
> following commands:
> 
>  >site<-c(rep(c("A1","A2", "RR1", "RR2", "WT1", "WT2", "WT3", "WT4"),8))
>  >H<-c(rep(c("exp", "exp", "prot", "pro", "exp", "exp", "exp", "exp"), 8))
>  >year<-c(rep(1989,8), rep(1990,8), rep(1991,8), rep(1993,8), 
> rep(1994,8), rep(1995,8), rep(1996,8), rep(1997,8))
>  >BA<-c(rep("bef",24), rep("after",40))
>  >abund<-runif(64, min=0, max=10)
>  >harvest<-data.frame(abund, BA, H, site, year)
> 
> Question 3.
> The authors use a complex anova design and here is part of their anova 
> table which shows the design and the df.
> 
> source.of.var      df
> Harvesting(H)     1, 6
> before-after(BA)  1, 6
> H x BA            1, 6
> Site(H)           6, 31
> Year(BA)          6, 31
> Site x BA         6, 31
> Year x H          6, 31
> Resid.            31
> 
> 
> My question again is: Why can't I reproduce the results? When I try a 
> simple anova without any random factors:
>  >harvest.lm<-lm(abund~H*BA+H/site+BA/year+site:BA+year:H)
> I get close enought but the results are different (at least I think they 
> are different because the nomin. df are different).
> 
> So obviously I need to assign sites as a random factor somehow. So when 
> I try to include site (which is nested in H) as a random factor and also 
> year nested in BA (as the authors did) the best I can come up with is:
>  >harvest.lme<-lme(abund~H*BA+BA/year+BA/year:H, random=~1|H/site)
> But I get a warning message (Warning message:In pf(q, df1, df2, 
> lower.tail, log.p) : NaNs produced) and also I don't know where to put 
> the site x BA term (whatever I try I get error messages).
> 
> I really apologise for the long post but after a week of studying and 
> trying as many ideas and examples I could find and think of I felt that 
> I really need advice from some more experienced users if I really want 
> to use this magnificent tool correctly.
> 
> Thanks in advance
> 
> -- 
> -------------------------------------------------------
> Nikolaos Lampadariou
> Hellenic Centre for Marine Research
> P.O. Box 2214
> 710 03 Heraklion Crete
> GREECE
> 
> e-mail: nlamp at her.hcmr.gr
> Ph. +30 281 0337849, +30 281 0337806
> FAX +30 281 0337822
> Web site: http://www.hcmr.gr/
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: http://www.nabble.com/before-after-control-impact-analysis-with-R-tp19024219p19027892.html
Sent from the R help mailing list archive at Nabble.com.