[R-sig-ME] R-sig-mixed-models Digest, Vol 44, Issue 19

[Dave Schoeman]

Sorry, all, there is a typo in the code submitted below. The SECOND line of code should read:

	H<-c(rep(c("exp", "exp", "pro", "pro", "exp", "exp", "exp", "exp"), 8))

i.e., the first code for protection should be "pro" like the rest, rather than "prot".

Very sorry about that.

- Dave


> Message: 2
> Date: Fri, 6 Aug 2010 13:17:31 +0100
> From: Dave Schoeman <david.schoeman at gmail.com>
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] Multiple before-after control-impact analysis
> Message-ID: <75069CDA-59D7-4A76-BA83-FE89F7044145 at gmail.com>
> Content-Type: text/plain; charset=us-ascii
> 
> I'd like to pick up a thread posted by Nikolaos Lampadariou in 2008. I am trying to analyse the MBACI design of Keough & Quinn (Legislative vs. practical protection of an 
> intertidal shoreline in southeastern Australia. Ecol. Appl. 2000. 10: 871-881). The design is simple: samples were taken at 8 Sites, of which 2 were "protected" and 6 were not; protection ended some time into the study; samples were taken annually at each site for three years prior to the end of protection and for 5 years afterwards. The response variable is a single value for each year-site combination. This leaves us with the fixed factors BA (two levels: before or after), H (two levels: protected or not) and Year (nested within BA). The random factor Site has 8 levels and is nested within H. Keough and Quinn use a repeated-measures ANOVA to fit the following terms: H, BA, H*BA, Site(H), Year(BA), Site(H)*BA and Year(BA)*H. Because there is no replication at the highest level interaction, this is left out of the model.
> 
> Nikolaos kindly provided code to simulate the resulting data:
> 	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)
> 
> He suggested fitting the model using lmer something like this:
> 	harvest.lmer<-lmer(abund~H*BA+BA/year+BA/year:H+H/site:BA+(1|H/site), harvest)
> 
> When I run this model, I get the dreaded "Error in mer_finalize(ans) : Downdated X'X is not positive definite, 12." error. Following recent advice on this thread, I tried dropping terms one at a time, and found that the problem lies with the H/site:BA term; omitting this allows the model to run without problem.
> 
> My questions are:
> 1 - Does the specified model look correct according to the experimental design described (or could the model specification causing the error)?
> 2 - Would I be justified in simply omitting the offending term?
> 
> Any comments appreciated.
> 
> - Dave
> R-sig-mixed-models mailing list
> R-sig-mixed-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 
> End of R-sig-mixed-models Digest, Vol 44, Issue 19
> **************************************************
>