Hi, I tried to do a nested Anova with the attached Data. My response variable is "survivors" and I would like to know the effect of (insect-egg clutch) "size", "position" (of clutch on twig) and "clone" (/plant genotype) on the survival of eggs (due to predation). Each plant was provided with three different sizes of clutches (45,15,5) and had pseudo-replications of size 15 and 5 on it. _Question 1: right code for nested Anova_ Here part of the code I used: clone1<-as.factor(ifelse(clone=="Gudrun",1,ifelse(clone=="Loden",2,ifelse(clone=="021",3,4)))) plant1<-as.factor(ifelse(plant==1,1,ifelse(plant==2,2,ifelse(plant==3,3,ifelse(plant==4,4,ifelse(plant==5,5,ifelse(plant==6,6, ifelse(plant==7,7,ifelse(plant==8,8,ifelse(plant==9,9,ifelse(plant==10,10,ifelse(plant==11,11,ifelse(plant==12,12, ifelse(plant==13,13,ifelse(plant==14,14,ifelse(plant==15,15,ifelse(plant==16,16,ifelse(plant==17,17,ifelse(plant==18,18, ifelse(plant==19,19,ifelse(plant==20,20,ifelse(plant==21,21,ifelse(plant==22,22,ifelse(plant==23,23,ifelse(plant==24,24, ifelse(plant==25,25,ifelse(plant==26,26,ifelse(plant==27,27,ifelse(plant==28,28,ifelse(plant==29,29,ifelse(plant==30,30, ifelse(plant==31,31,ifelse(plant==32,32,ifelse(plant==33,33,ifelse(plant==34,34,ifelse(plant==35,35,ifelse(plant==36,36, ifelse(plant==37,37,ifelse(plant==38,38,ifelse(plant==39,39,ifelse(plant==40,40,ifelse(plant==41,41,ifelse(plant==42,42,43))))))))))))))))))))))))))))))))))))))))))) size1<-as.factor(ifelse(size==5,1,ifelse(size==15,2,3))) position1<-as.factor(ifelse(position==1,1,ifelse(position==2,2,ifelse(position==3,3, ifelse(position==4,4,ifelse(position==5,5,ifelse(position==6,6,ifelse(position==7,7, ifelse(position==8,8,ifelse(position==9,9,ifelse(position==10,10,ifelse(position==11,11,ifelse(position==12,12,13))))))))))))) ANOVA<-aov(survivors~(clone1*size1*position1)+Error(plant1/(size1*position1))) After that command it says: "Error() model is singular". Even after googling (looks like many people had that problem) I am still uncertain how to solve that, mainly because I am not sure if the used code for the Anova meets my set up requirements. If I am using "ANOVA<-aov(survivors~clone1*plant1/(size1*position1)))" R sometimes is even crashing down. _Question 2: normal distribution and homogeneity of variance_ So fare I tested "normal" Anovas by doing a linear model following plotting and KS-test and Bartlett-test including Q-Q-plot. fm<-lm(survivors~clone1*plant1*size1*position1) Residuen<-resid(fm) Here R is crashing again!!! ks.test(Residuen,mean(Residuen),sd(Residuen)) plot(density(Residuen)) qqnorm(Residuen) par(mfrow=c(2,2)) plot(fm) bartlett.test(survivors~clone1*plant1*position1*size1) Is that too much to compute for R? Is it OK to look for the preconditions for that nested Anova in that way, or do I have to use the same code for the Bartlett also (or evaluate only by looking at the Graphs (Zuur et al. 2010))?????? Thank you very much in advance J?rg -------------- next part -------------- An embedded and charset-unspecified text was scrubbed... Name: survival.txt URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20110912/2103af20/attachment.txt>
joerg stephan <stephanj <at> rhrk.uni-kl.de> writes:> > Hi, > > I tried to do a nested Anova with the attached Data. My response > variable is "survivors" and I would like to know the effect of > (insect-egg clutch) "size", "position" (of clutch on twig) and "clone" > (/plant genotype) on the survival of eggs (due to predation). Each plant > was provided with three different sizes of clutches (45,15,5) and had > pseudo-replications of size 15 and 5 on it.This may not be what you wanted at all, but here is my take on what you could try with these data. You have binomial data with lots of zeros and ones (i.e. 0 or 100% survival), so they are unlikely to be transformable to normality however you try. You *might* be able to get away with modeling these as normally distributed and then rely on permutation tests to get your significance levels right, but (tricky as it is) I actually think GLMMs (Zuur et al 2009, Bolker et al 2009) are the best way to go ... ## read in the data x <- read.table("aovmisc.dat",header=TRUE) ## take a look summary(x) ## with(x,table(clone,as.factor(size),as.factor(position))) ## make categorical versions of size & plant variables, ## and reconstitute the "total number surviving" variable x <- transform(x, fsize=factor(size), fplant=factor(plant), tsurv=round(survivors*size) ) library(ggplot2) ## trying to plot everything ... zspace <- opts(panel.margin=unit(0,"lines")) ## squash panels together theme_update(theme_bw()) ## white background ## plot all data: survival vs position, with all plants represented; ## separate subplots for each clone ## overlay GLM fit ggplot(x,aes(x=position,y=survivors,size=size,colour=fplant))+ geom_point(alpha=0.5)+ facet_wrap(~clone)+xlim(0,15)+geom_line(size=0.5,alpha=0.2)+zspace+ geom_smooth(aes(group=1,weight=size), colour="black",method="glm",family="binomial") ## it looks like there is a positive effect of position, and also ## a positive effect of size (large bubbles toward high-survival) ## there are a few outliers in 'position' -- are these typos? subset(x,position>15) ## focus on size and clone instead, suppress position: ggplot(x,aes(x=fsize,y=survivors,colour=clone))+stat_sum(aes(size=..n..))+ facet_grid(.~clone)+ geom_smooth(method="glm",family="binomial",aes(weight=size, x=as.numeric(fsize))) ## OR ## + stat_summary(fun.y="mean",geom="line",aes(x=as.numeric(fsize)))+zspace library(mgcv) ggplot(x,aes(x=position,y=survivors,colour=clone))+stat_sum(aes(size=..n..))+ facet_grid(fsize~clone)+ geom_smooth(method="gam",family="binomial")+xlim(0,15)+zspace ## appears to be an effect of size, and a clone*size interaction. ## possible bimodal distribution (0/1) at size=5? with(x,table(clone,plant)) ## each plant is only in one clone ## (explicit nesting) library(lme4) ## fit with size*clone interaction, main effect of position ## no need to nest plant in clone because they are explicitly nested ## (i.e. unique plant IDs) g1 <- glmer(cbind(tsurv,size)~fsize*clone+position+(1|fplant), data=x, family="binomial") ## doesn't LOOK like overdispersion ... sum(residuals(g1)^2) ## >> residual df nrow(x)-length(fixef(g1)) pchisq(778,df=537) ## ... but try fitting observation-level random effect anyway x <- transform(x,obs=seq(nrow(x))) g2 <- update(g1,.~.+(1|obs)) summary(g2) ## among-observation variance 2x among-plant variance anova(g1,g2) ## looks much better (10 AIC units lower / p <<< 0.05 ## drop1(g2,test="Chisq") ## looks like we can drop the size*clone interaction? g3 <- update(g2,.~.-clone:fsize) drop1(g3,test="Chisq") ## strong effects of size, position; ## weak effects of clone fixef(g3)>