knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
What factors are there? Fixed or random?
He samples in seven estuaries along the New South Wales coast (three of which are "Pristine", four are "Modified"), and in each estuary, he takes 4-7 samples and counts the creepy crawlies therein.
Modification is a factor, taking levels "Pristine" and "Modified"
Estuary is a factor, taking seven levels. If these were sampled randomly, and we want to make inferences across all estuaries on the New South Wales coast, it oculd be treated as a random factor.
The 4-7 samples at each estuary are the replicates, so they shouldn't be added to the model, variation in replicates will enter through the error term.
library(ecostats) data(estuaries) plot(Total~Estuary,data=estuaries,col=c(4,2,2,4,2,4,2)) legend("bottomleft",legend=c("Modified","Pristine"),col=c(4,2),pch=15,pt.cex=2)
library(ecostats) data(estuaries) library(lme4) ft_estu = lmer(Total~Mod+(1|Estuary),data=estuaries) summary(ft_estu)
There is some evidence of an effect of Mod
, since the estimated coefficient is more than double its standard error (so a 95\% confidence interval would not cover zero). The effect appears to be a decrease in total abundance in pristine estuaries.
par(mfrow=c(1,2),mar=c(3,3,1,1),mgp=c(1.75,0.75,0)) ft_estu = lmer(Total~Mod+(1|Estuary),data=estuaries) scatter.smooth(residuals(ft_estu)~fitted(ft_estu), xlab="Fitted values",ylab="Residuals") abline(h=0,col="red") scatter.smooth(residuals(ft_estu)~predict(ft_estu,re.form=NA), xlab="Fitted values (no random effects)",ylab="Residuals") abline(h=0,col="red")
anova
to compare mixed effects models for the estuary dataft_estu = lmer(Total~Mod+(1|Estuary),data=estuaries,REML=F) ft_estuInt = lmer(Total~(1|Estuary),data=estuaries,REML=F) anova(ft_estuInt,ft_estu)
There is some evidence of an effect of modification.
confint(ft_estu)
rft=ranef(ft_estu,condVar=T) library(lattice) dotplot(rft)
estuaries$isMod = as.numeric(estuaries$Mod=="Modified") estuaries$isPri = as.numeric(estuaries$Mod!="Modified") ft_estuDiff = lmer(Total~Mod+(0+isMod|Estuary)+(0+isPri|Estuary),data=estuaries,REML=F) summary(ft_estuDiff) BIC(ft_estu,ft_estuDiff)
BIC suggests that we didn't need a different variance term for each level of Mod
. (It also estimated the cross-estuary variance to be zero for modified estuaries, leading to a warning in the output.)
data(aphidsBACI) str(aphidsBACI)
OK so we are looking for a Treatment:Time
interaction, but to account for repeated measures of each plot, we want a random effect for Plot
in the model.
ft_aphids=lmer(logcount~Treatment*Time+(1|Plot),data=aphidsBACI) ft_aphidNull=lmer(logcount~Time+(1|Plot),data=aphidsBACI) anova(ft_aphidNull,ft_aphids)
Which gives us marginal evidence of an effect of bird exclusion of aphid counts.
Compare this to what we got when adding Plot
as a fixed effect:
lm_aphids=lm(logcount~Plot+Treatment*Time,data=aphidsBACI) anova(lm_aphids)
Interestingly, this $P$-value is a lot smaller.
summary(lm_aphids) summary(ft_aphids)
We get the same point estimate for the Treatment:Time
effect, but the standard error is slightly smaller in the random effects model.
data(estuaryZone) cols=c("blue","red","lightblue","pink") plot(Total~interaction(Estuary,Zone),data=estuaryZone,col=cols[c(1,2,2,1,2,1,2,3,4,4,3,4,3,4)]) legend("bottomright",legend=c("Mod-Inner","Prist-Inner","Mod-Outer","Pris-Outer"),col=cols,pch=15,pt.cex=2)
It looks like there is an effect of Modification, not sure if there is an interaction (the effect seems more striking in Outer than Inner zones)
par(mfrow=c(1,2),mar=c(3,3,1,1),mgp=c(1.75,0.75,0)) library(lme4) lme_MZ = lmer(Total~Zone*Mod + (Zone|Estuary), data=estuaryZone ) scatter.smooth(residuals(lme_MZ)~fitted(lme_MZ), xlab="Fitted values",ylab="Residuals") abline(h=0,col="red") scatter.smooth(residuals(lme_MZ)~predict(lme_MZ,re.form=NA), xlab="Fitted values (no random effects)",ylab="Residuals") abline(h=0,col="red")
There is a suggestion of less total abundance as fitted values increase, which is super-weird. But it's not too alarming...
lme_MplusZ = lmer(Total~Zone+Mod + (Zone|Estuary), data=estuaryZone ) anova(lme_MplusZ,lme_MZ)
No evidence of an interaction between Zone and Modification. Testing for a Mod
main effect:
lme_Z = lmer(Total~Zone + (Zone|Estuary), data=estuaryZone ) anova(lme_Z,lme_MplusZ)
There is strong evidence that total abundance is different between Modified and Pristine estuaries. The boxplot suggests abundance is higher in Modified estuaries, and looking at the data, this appears to be mostly due to high counts of Balanus.variegatus
, especially in outer modified zones.
Mod
fixed effect in Exercise 6.1.nBoot=500 bStat=rep(NA,nBoot) ft_estu = lmer(Total~Mod+(1|Estuary),data=estuaries) for(iBoot in 1:nBoot) { estuaries$TotalSim=unlist(simulate(ft_estu)) ft_i = lmer(TotalSim~Mod+(1|Estuary),data=estuaries) bStat[iBoot] = fixef(ft_i)[2] } sd(bStat) #standard error of Mod effect
And if we compare this to the standard error from summary
:
summary(ft_estu)
we see the estimated standard error is r summary(ft_estu)$coef[2,2]
, which is pretty close to the value we got by simulation.
ft_noestu = lm(Total~Mod,data=estuaries) library(ecostats) anovaPB(ft_noestu,ft_estu,n.sim=99)
So we have no evidence of an Estuary
effect.
Use the parametric bootstrap to get a formal test for a zone:mod
interaction.
We can just run the old analysis and change from anova
to anovaPB
:
lme_MZ = lmer(Total~Zone*Mod + (Zone|Estuary), data=estuaryZone, REML=FALSE ) lme_MplusZ = lmer(Total~Zone+Mod + (Zone|Estuary), data=estuaryZone, REML=FALSE ) anovaPB(lme_MplusZ,lme_MZ,n.sim=99)
There is no evidence of an interaction.
(Ignore the warnings in the output -- this is random stuff that was thrown up in bootstrap resamples that didn't get a good fit.)
How do results compare to those from when you were using the anova
function?
Results are similar to what we saw before. The only thing that is different is the $P$-value, but it is very similar (suggesting there was no need for a parametric bootstrap here!).
This would all have been so much easier if there wasn't a random effect in the model... do we really need Estuary
in there?
lme_MZ = lmer(Total~Zone*Mod + (Zone|Estuary), data=estuaryZone, REML=FALSE ) lme_MZnoest = lm(Total~Zone+Mod, data=estuaryZone) anovaPB(lme_MZnoest,lme_MZ,n.sim=99)
We have no evidence of an Estuary
effect either!
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.