############################################################################ #####Adult C. zealandica abundance and sex ratio through time ###Awatere Valley 2014 rm(list=ls()) X<-read.table('chisq-sexKN14-LM.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) M=P.total.M F=P.total.F d=days detach(X) hist(M,breaks=seq(0,120,10),col="green",main="Histogram") plot(density(M),main="Density function") lm1<-lm(M~d) par(mfrow=c(2,2)) plot(lm1) resid(lm1) names(lm1) par(mfrow=c(1,1)) hist(resid(lm1),data=X) summary(lm1) m=M/100 f=F/100 m5<-c(0.99999, 0.5, 0.99999, 0.683, 0.525, 0.5714286, 0.6, 0.85, 0.85, 0.66667, 0.73333, 0.6833333, 0.6, 0.63333, 0.5166667, 0.75, 0.65, 0.616667, 0.63333, 0.65, 0.466667, 0.5166667, 0.65, 0.6, 0.466667, 0.43333, 0.55, 0.55, 0.483333) library(betareg) glm10<-betareg(m5~d, link="logit") par(mfrow=c(2,2)) plot(glm10) resid(glm10) names(glm10) par(mfrow=c(1,1)) hist(resid(glm10),data=X) summary(glm10) plot(m5~d, pch=16, type="p", ylab="Proportion of males", xlab="Days after flight begun") newX<-expand.grid(d=seq(from=1,to=35,by=0.1)) newY<-predict(glm10,newdata=newX,type="response") newY1<-predict(glm10,newdata=newX, type="variance") sd=sqrt(newY1) se=sd/(sqrt(29)) addThese<-data.frame(newX,newY,se) addThese lines(newY~d,data=addThese,col="red",lwd=2) #Aqui no agregue el intervalo de confianza, pero puede ser agregado lines(newY+1.96*se~d, data=addThese,col="red",lty=2,lwd=2) #al agregar 1.96*se.fit calculo el intervalo de confianza al 95%. Si solo dejo se.fit, me muestra el error estandart. lines(newY-1.96*se~d, data=addThese,col="red",lty=2,lwd=2) #####Awatere Valley 2015 rm(list=ls()) X<-read.table('chisq-sexKN15-LM.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) M1=P.total.M F1=P.total.F d1=days detach(X) hist(M1,breaks=seq(0,120,10),col="green",main="Histogram") plot(density(M1),main="Density function") lm3<-lm(M1~d1) par(mfrow=c(2,2)) plot(lm3) resid(lm3) names(lm3) par(mfrow=c(1,1)) hist(resid(lm3),data=X) summary(lm3) m1=M1/100 f1=F1/100 library(betareg) glm11<-betareg(m1~d1, link="logit") par(mfrow=c(2,2)) plot(glm11) resid(glm11) names(glm11) par(mfrow=c(1,1)) hist(resid(glm11),data=X) summary(glm11) plot(m1~d1, pch=16, type="p", ylab="Proportion of males", xlab="Days after flight begun") newX<-expand.grid(d1=seq(from=1,to=25,by=0.1)) newY<-predict(glm11,newdata=newX,type="response") newY1<-predict(glm11,newdata=newX, type="variance") sd=sqrt(newY1) se=sd/(sqrt(20)) addThese<-data.frame(newX,newY,se) addThese lines(newY~d1,data=addThese,col="red",lwd=2) #Aqui no agregue el intervalo de confianza, pero puede ser agregado lines(newY+1.96*se~d1, data=addThese,col="red",lty=2,lwd=2) #al agregar 1.96*se.fit calculo el intervalo de confianza al 95%. Si solo dejo se.fit, me muestra el error estandart. lines(newY-1.96*se~d1, data=addThese,col="red",lty=2,lwd=2) #####Blenheim 2015 rm(list=ls()) X<-read.table('chisq-sexWH15-LM.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) M2=P.total.M F2=P.total.F d2=days detach(X) hist(M2,breaks=seq(0,120,10),col="green",main="Histogram") plot(density(M2),main="Density function") lm4<-lm(M2~d2) par(mfrow=c(2,2)) plot(lm4) resid(lm4) names(lm4) par(mfrow=c(1,1)) hist(resid(lm4),data=X) summary(lm4) m2=M2/100 f2=F2/100 library(betareg) glm12<-betareg(m2~d2, link="logit") par(mfrow=c(2,2)) plot(glm12) resid(glm12) names(glm12) par(mfrow=c(1,1)) hist(resid(glm12),data=X) summary(glm12) plot(m2~d2, pch=16, type="p", ylab="Proportion of males", xlab="Days after flight begun") newX<-expand.grid(d2=seq(from=1,to=37,by=0.1)) newY<-predict(glm12,newdata=newX,type="response") newY1<-predict(glm12,newdata=newX, type="variance") sd=sqrt(newY1) se=sd/(sqrt(29)) addThese<-data.frame(newX,newY,se) addThese lines(newY~d2,data=addThese,col="red",lwd=2) #Aqui no agregue el intervalo de confianza, pero puede ser agregado lines(newY+1.96*se~d2, data=addThese,col="red",lty=2,lwd=2) #al agregar 1.96*se.fit calculo el intervalo de confianza al 95%. Si solo dejo se.fit, me muestra el error estandart. lines(newY-1.96*se~d2, data=addThese,col="red",lty=2,lwd=2) ########################################################################### #######Adult removal and daily sex-ratio trends ######Awatere Valley 2015 rm(list=ls()) X<-read.table('5minKN15.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) m17=mean t17=c(5,10,15,20) detach(X) library(betareg) glm17<-betareg(m17~t17, link="logit") par(mfrow=c(2,2)) plot(glm17) resid(glm17) names(glm17) par(mfrow=c(1,1)) hist(resid(glm17),data=X) summary(glm17) attach(X) d=mean e=c(5,10,15,20) f=se detach(X) x1<-e y1<-d ybar1<-f par(mar=c(5,6,4,2)) xy.error.bars1<-function (x1,y1,ybar1){ plot(x1, y1, pch=16, col="black", xlab = "Time period (min)", ylab="Proportion of males", ylim=c(min(0),max(0.4))) arrows(x1, y1-ybar1, x1, y1+ybar1, code=3, angle=90, length=0.05)} xy.error.bars1(x1,y1,ybar1) points(x1, y1,col="black",pch=16) newX<-expand.grid(t17=seq(from=5,to=20,by=0.1)) newY<-predict(glm17,newdata=newX,type="response") newY1<-predict(glm17,newdata=newX, type="variance") sd=sqrt(newY1) se=sd/(sqrt(4)) addThese<-data.frame(newX,newY,se) addThese lines(newY~t17,data=addThese,col="red",lwd=2) #Aqui no agregue el intervalo de confianza, pero puede ser agregado lines(newY+1.96*se~t17, data=addThese,col="red",lty=2,lwd=2) #al agregar 1.96*se.fit calculo el intervalo de confianza al 95%. Si solo dejo se.fit, me muestra el error estandart. lines(newY-1.96*se~t17, data=addThese,col="red",lty=2,lwd=2) #####Blenheim 2015 rm(list=ls()) X<-read.table('5minWH15.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) m18=mean t18=c(5,10,15,20,25) detach(X) library(betareg) glm18<-betareg(m18~t18, link="logit") par(mfrow=c(2,2)) plot(glm18) resid(glm18) names(glm18) par(mfrow=c(1,1)) hist(resid(glm18),data=X) summary(glm18) attach(X) d=mean e=c(5,10,15,20,25) f=se detach(X) x1<-e y1<-d ybar1<-f par(mar=c(5,6,4,2)) xy.error.bars1<-function (x1,y1,ybar1){ plot(x1, y1, pch=16, col="black", xlab = "Time period (min)", ylab="Proportion of males", ylim=c(min(0),max(0.4))) arrows(x1, y1-ybar1, x1, y1+ybar1, code=3, angle=90, length=0.05)} xy.error.bars1(x1,y1,ybar1) points(x1, y1,col="black",pch=16) newX<-expand.grid(t18=seq(from=5,to=25,by=0.1)) newY<-predict(glm18,newdata=newX,type="response") newY1<-predict(glm18,newdata=newX, type="variance") sd=sqrt(newY1) se=sd/(sqrt(5)) addThese<-data.frame(newX,newY,se) addThese lines(newY~t18,data=addThese,col="red",lwd=2) #Aqui no agregue el intervalo de confianza, pero puede ser agregado lines(newY+1.96*se~t18, data=addThese,col="red",lty=2,lwd=2) #al agregar 1.96*se.fit calculo el intervalo de confianza al 95%. Si solo dejo se.fit, me muestra el error estandart. lines(newY-1.96*se~t18, data=addThese,col="red",lty=2,lwd=2) ############################################################################ #########Sex-ratio and adult abundance correlation #####Awatere Valley 2014 rm(list=ls()) X<-read.table('corKN14.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) a=ab.KN14 b=prop.KN14/100 detach(X) cor.test(a,b, method="spearm", exact=FALSE) #####Awatere Valley 2015 rm(list=ls()) X<-read.table('corKN15.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) a=ab.KN15 b=prop.KN15/100 detach(X) cor.test(a,b, method="spearm", exact=FALSE) #####Blenheim 2015 rm(list=ls()) X<-read.table('corWH15.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) a=ab.WH15 b=prop.WH15/100 detach(X) cor.test(a,b, method="spearm", exact=FALSE) ########################################################################### ########5-min removal experiment ########abundance and sex-ratio correlation ####Blenheim 2015 rm(list=ls()) X<-read.table('corWH15abun.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) a=abun b=prop.SR detach(X) cor.test(a,b, method="spearm", exact=FALSE) ####Awatere Valley 2015 rm(list=ls()) X<-read.table('corKN15abun.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) c=abun d=prop.SR detach(X) cor.test(c,d, method="spearm", exact=FALSE) ############################################################################ ########5-min removal experiment ########adult abundance at each time period ####Blenheim 2015 rm(list=ls()) X<-read.table('5minabundance.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) a=mean b=c(5,10,15,20,25) c=se detach(X) x<-b y<-a ybar<-c #####Awatere Valley 2015 X<-read.table('5minabunKN15.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) d=mean e=c(5,10,15,20) f=se detach(X) x1<-e y1<-d ybar1<-f par(mar=c(5,6,4,2)) xy.error.bars1<-function (x,y, ybar,x1,y1,ybar1){ plot(x1, y1, pch=16, col="darkorange2") plot(x, y, pch=17, col="dodgerblue3", xlab = "Time period (min)", ylab=expression(Adult~abundance~(beetles~plant^{-1})), ylim=c(min(0),max(30))) arrows(x, y-ybar, x, y+ybar, code=3, angle=90, length=0.05) arrows(x1, y1-ybar1, x1, y1+ybar1, code=3, angle=90, length=0.05)} xy.error.bars1(x,y,ybar,x1,y1,ybar1) points(x1, y1,col="darkorange2",pch=16) points(x, y,col="dodgerblue3",pch=17) sequence<-order(x) lines(x[sequence],y[sequence],lwd=2,col="dodgerblue3") sequence<-order(x1) lines(x1[sequence],y1[sequence],lwd=2,lty=2,col="darkorange2") text(5,5, "a",cex=1) text(10,21, "b",cex=1) text(15,28.5, "b",cex=1) text(15,11.5, "ab",cex=1) text(20,10.5, "ab",cex=1) text(20,4, "a",cex=1) text(25,3, "a",cex=1) ############################################################################## ########Tukey contrasts for the lattest experiments ##########Awatere Valley 2015 rm(list=ls()) X<-read.table('ANOVA5minSR.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) hist(abundance,breaks=seq(0,150,10),col="green",main="Histogram") ### Dice que va de 0 a 7, con el eje cada 0.1 plot(density(abundance),main="Density function") ### densidad de frecuencia detach(X) glm1<-glm(abundance~period, family=poisson, data=X) par(mfrow=c(2,2)) plot(glm1) resid(glm1) names(glm1) par(mfrow=c(1,1)) hist(resid(glm1),data=X) summary(glm1) library(MASS) glm2<-glm.nb(abundance~period, data=X) par(mfrow=c(2,2)) plot(glm2) resid(glm2) names(glm2) par(mfrow=c(1,1)) hist(resid(glm2),data=X) summary(glm2) sum(resid(glm2,type="pearson")^2) 1 - pchisq(summary(glm2)$deviance, summary(glm2)$df.residual) ###To see if the model fits the data ###The result is the p-value. p>0.05 then ###the model fits the data library(multcomp) tukey<-(glht(glm2, mcp(period="Tukey"))) summary(tukey) par(mar=c(5,7,4,2)) plot(tukey) #############Blenheim 2015 rm(list=ls()) X<-read.table('ANOVA5minSRKN15.txt', sep='\t', dec=',',h=T) names(X) str(X) dim(X) attach(X) hist(abundance,breaks=seq(0,150,10),col="green",main="Histogram") ### Dice que va de 0 a 7, con el eje cada 0.1 plot(density(abundance),main="Density function") ### densidad de frecuencia detach(X) glm1<-glm(abundance~period, family=poisson, data=X) par(mfrow=c(2,2)) plot(glm1) resid(glm1) names(glm1) par(mfrow=c(1,1)) hist(resid(glm1),data=X) summary(glm1) library(MASS) glm2<-glm.nb(abundance~period, data=X) par(mfrow=c(2,2)) plot(glm2) resid(glm2) names(glm2) par(mfrow=c(1,1)) hist(resid(glm2),data=X) summary(glm2) sum(resid(glm2,type="pearson")^2) 1 - pchisq(summary(glm2)$deviance, summary(glm2)$df.residual) ###To see if the model fits the data ###The result is the p-value. p>0.05 then ###the model fits the data library(multcomp) tukey<-(glht(glm2, mcp(period="Tukey"))) summary(tukey) par(mar=c(5,7,4,2)) plot(tukey)