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: [EMAIL PROTECTED]
> Ph. +30 281 0337849, +30 281 0337806
> FAX +30 281 0337822
> Web site: http://www.hcmr.gr/
>
> ______________________________________________
> R-help@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.
______________________________________________
R-help@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.
