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