rm(list=ls()) setwd("C:/Users/Natalie/Documents/Research_UMBC/Reinforcement") stream<-read.csv("Stream analysis.csv") size<-read.csv("Size_noEF.csv") rein<-read.csv("Rein data.csv") fem<-rein[c(1:72),] mal<-rein[c(73:142),] fMFR<-fem[c(1:18),] fFCT<-fem[c(19:36),] fEFB<-fem[c(37:54),] fLNC<-fem[c(55:72),] fAllo<-fem[c(1:36),] fSym<-fem[c(37:72),] mMFR<-mal[c(1:18),] mFCT<-mal[c(19:36),] mEFB<-mal[c(37:52),] mLNC<-mal[c(53:70),] mAllo<-mal[c(1:36),] mSym<-mal[c(37:70),] Allo<-rein[c(1:36,73:108),] Sym<-rein[c(37:72,109:142),] #################################### ##Libraries #lsmeans if(!require(FSA)){install.packages("FSA")} if(!require(psych)){install.packages("psych")} if(!require(lsmeans)){install.packages("lsmeans")} if(!require(car)){install.packages("car")} library(FSA) library(lsmeans) #mann-whitney u tests library(coin) #effect sizes library(effsize) #wilcox t tests library(exactRankTests) #gls library(nlme) ####################################################### #Results #Size difference between stimulus pairs head(size) shapiro.test(size$S1) shapiro.test(size$S2) qqnorm(size$S1) qqline(size$S1) qqnorm(size$S2) qqline(size$S2) wilcoxsign_test(size$S1~size$S2,paired=T) #Female sym v. female allo shapiro.test(fAllo$SOP_5cm) qqnorm(fAllo$SOP_5cm) qqline(fAllo$SOP_5cm) shapiro.test(fSym$SOP_5cm) qqnorm(fSym$SOP_5cm) qqline(fSym$SOP_5cm) g=factor(c(rep("fAllo$SOP_5cm",length(fAllo$SOP_5cm)),rep("fSym$SOP_5cm",length(fSym$SOP_5cm)))) v=c(fAllo$SOP_5cm,fSym$SOP_5cm) wilcox_test(v~g, distribution="exact") cohen.d(fSym$SOP_5cm,fAllo$SOP_5cm) mean(fSym$SOP_5cm) f<-sd(fSym$SOP_5cm) f/sqrt(length(fSym$SOP_5cm)) mean(fAllo$SOP_5cm) f<-sd(fAllo$SOP_5cm) f/sqrt(length(fAllo$SOP_5cm)) #Male sym v. male allo shapiro.test(mAllo$SOP_5cm) qqnorm(mAllo$SOP_5cm) qqline(mAllo$SOP_5cm) shapiro.test(mSym$SOP_5cm) qqnorm(mSym$SOP_5cm) qqline(mSym$SOP_5cm) wilcox.test(mSym$SOP_5cm,mAllo$SOP_5cm,paired=F) g=factor(c(rep("mAllo$SOP_5cm",length(mAllo$SOP_5cm)),rep("mSym$SOP_5cm",length(mSym$SOP_5cm)))) v=c(mAllo$SOP_5cm,mSym$SOP_5cm) wilcox_test(v~g, distribution="exact") cohen.d(mSym$SOP_5cm,mAllo$SOP_5cm) mean(mSym$SOP_5cm) f<-sd(mSym$SOP_5cm) f/sqrt(length(mSym$SOP_5cm)) mean(mAllo$SOP_5cm) f<-sd(mAllo$SOP_5cm) f/sqrt(length(mAllo$SOP_5cm)) #SOP EF females summary(fEFB) shapiro.test(fEFB$SOP_5cm) t.test(fEFB$SOP_5cm,mu=0,type="two.tail") mean(fEFB$SOP_5cm) sd(fEFB$SOP_5cm)/sqrt(length(fEFB$SOP_5cm)) #SOP EF males summary(mEFB) shapiro.test(mEFB$SOP_5cm) t.test(mEFB$SOP_5cm,mu=0,type="two.tail") mean(mEFB$SOP_5cm) sd(mEFB$SOP_5cm)/sqrt(length(mEFB$SOP_5cm)) #SOP LC females summary(fLNC) shapiro.test(fLNC$SOP_5cm) t.test(fLNC$SOP_5cm,mu=0,type="two.tail") mean(fLNC$SOP_5cm) sd(fLNC$SOP_5cm)/sqrt(length(fLNC$SOP_5cm)) #SOP LC males summary(mLNC) qqnorm(mLNC$SOP_5cm) qqline(mLNC$SOP_5cm) shapiro.test(mLNC$SOP_5cm) wilcox.test(mLNC$SOP_5cm,mu=0,alternative="two.sided") mean(mLNC$SOP_5cm) sd(mLNC$SOP_5cm)/sqrt(length(mLNC$SOP_5cm)) qnorm(1-(0.003485/2)) #SOP fMR summary(fMFR) shapiro.test(fMFR$SOP_5cm) t.test(fMFR$SOP_5cm,mu=0,type="two.tail") mean(fMFR$SOP_5cm) sd(fMFR$SOP_5cm)/sqrt(length(fMFR$SOP_5cm)) #SOPmMR summary(mMFR) shapiro.test(mMFR$SOP_5cm) t.test(mMFR$SOP_5cm,mu=0,alternative="two.sided") mean(mMFR$SOP_5cm) sd(mMFR$SOP_5cm)/sqrt(length(mMFR$SOP_5cm)) #SOP f FC summary(fFCT) shapiro.test(fFCT$SOP_5cm) t.test(fFCT$SOP_5cm,mu=0,type="two.tail") mean(fFCT$SOP_5cm) sd(fFCT$SOP_5cm)/sqrt(length(fFCT$SOP_5cm)) #SOP mFC summary(mFCT) shapiro.test(mFCT$SOP_5cm) t.test(mFCT$SOP_5cm,mu=0,alternative="two.sided") mean(mFCT$SOP_5cm) sd(mFCT$SOP_5cm)/sqrt(length(mFCT$SOP_5cm)) ############################################################# ##Modeling full<-aov(SOP_5cm ~ Context * Context/Pop * Sex, data=rein) summary(full) AIC(full) r1<-aov(SOP_5cm ~ Context * Context/Pop + Sex, data=rein) summary(r1) AIC(r1) r2<-aov(SOP_5cm ~ Context/Pop * Sex -Context:Sex, data=rein) summary(r2) AIC(r2) r3<-aov(SOP_5cm ~ Context,data=rein) summary(r3) AIC(r3) #final model = r2 #r2 AIC - 225.9 (within 3 from minimal model (r3)) #LS means population comparison lsmeans(r2,pairwise ~ Context, adjust ="bonferroni") #P = 0.006, significant difference between sym and allo t1<-aov(SOP_5cm ~ Context + Context/Pop * Sex, data=rein) summary(t1) lsmeans(t1,pairwise ~ Pop|Context:Sex, adjust ="bonferroni") #separate for sexes #Stream trials stream.1<-stream[c(1:5,18)] zs<-stream.1[c(10:12),] za<-stream.1[c(4:6,16:18),] m1<-aov(Iso_index ~ context/pop *species,data=stream.1) summary(m1) AIC(m1) m1.1<-aov(Iso_index ~ context+ pop +species,data=stream.1) summary(m1.1) AIC(m1.1) m1.2<-aov(Iso_index ~ context+ pop,data=stream.1) summary(m1.2) AIC(m1.2) m1.3<-aov(Iso_index ~ context,data=stream.1) summary(m1.3) AIC(m1.3) #context only significant term, no real difference in AIC #adding in random effects m2<-gls(Iso_index ~ context, method = "REML", data=stream.1) summary(m2) m3<-lme(Iso_index ~ context, random =~1|trial, method= "REML",data=stream.1) m4<-lme(Iso_index ~ context, random =~1|trial/replicate, method= "REML",data=stream.1) summary(m3) summary(m4) AIC(m2) AIC(m3) AIC(m4) plot(stream.1$context,stream.1$Iso_index) #heterogenous data, need to account using GLS vf2 <- varIdent(form= ~ 1 | context) m5<-gls(Iso_index ~ context, method = "REML", data=stream.1,weights=vf2) summary(m5) m6<-lme(Iso_index ~ context, random =~1|trial,weights=vf2, method= "REML",data=stream.1) m7<-lme(Iso_index ~ context, random =~1|trial/replicate, weights=vf2,method= "REML",data=stream.1) AIC(m5) AIC(m6) AIC(m7) AIC(m5) - AIC(m7) summary(m6) anova(m5,m6) anova(m5,m7) # #Do i need a random effect? m2.lm<-lm(Iso_index ~ context, data=stream.1) summary(m2.lm) E<-rstandard(m2.lm) boxplot(E~trial,data=stream.1) boxplot(E~replicate,data=stream.1) boxplot(E~trial/replicate,data=stream.1) final.1<-gls(Iso_index ~ context, method = "REML", data=stream.1,weights=vf2) final.2<-lme(Iso_index ~ context, random =~1|trial/replicate, weights=vf2,method= "REML",data=stream.1) anova(final.1,final.2) final.3<-lme(Iso_index ~ context, random =~1|trial, weights=vf2,method= "REML",data=stream.1) anova(final.1,final.3) anova(final.2,final.3) AIC(final.1) AIC(final.2) AIC(final.3) summary(final.1) E<-resid(final.1,type="normalized") F<-fitted(final.1) plot(x=F,y=E) boxplot(E~context,data=stream.1) #Does solicitaitons have large effect? head(stream) nosol<-stream[-c(10:13,16,17)] nosol[,6]+nosol[,8]+nosol[,10] nosol$totC<-nosol[,6]+nosol[,6]+nosol[,10] nosol$totH<-nosol[,7]+nosol[,9]+nosol[,11] nosol$Iso<-(nosol$totC - nosol$totH)/(nosol$totC + nosol$totH) sol <-gls(Iso ~ context,weights = vf2, method = "REML",data = nosol) AIC(sol) summary(sol) ############################################## #Figures #Context v SOP, combined sexes op<-par(mar=c(5,7,4,2)) boxplot(Sym$SOP_5cm,Allo$SOP_5cm,ylim=c(-1,1),axes=F) title(xlab= "Context",ylab= "Strength of preference",cex.lab=1.5,line=3.5) box() axis(1,at=seq(1,2,by=1),labels=c("Sympatric","Allopatric"),cex.axis=1.25) axis(2, at=seq(-1,1, by=.5),las=2,cex.axis=1.25) par(op) text(1.5,.95, expression("*"),cex=2) #Context v SOP, separate sexes op<-par(mar=c(5,7,4,2)) boxplot(fSym$SOP_5cm,fAllo$SOP_5cm,ylim=c(-1,1),axes=F) title(xlab= "Context",ylab= "Strength of preference",cex.lab=1.5,line=3.5) box() axis(1,at=seq(1,2,by=1),labels=c("Sympatric","Allopatric"),cex.axis=1.25) axis(2, at=seq(-1,1, by=.5),las=2,cex.axis=1.25) par(op) text(.55,.99, expression("(a)"),cex=1.25) text(1.5,.95, expression("*"),cex=2) op<-par(mar=c(5,7,4,2)) boxplot(mSym$SOP_5cm,mAllo$SOP_5cm,ylim=c(-1,1),axes=F) title(xlab= "Context",ylab= "Strength of preference",cex.lab=1.5,line=3.5) box() axis(1,at=seq(1,2,by=1),labels=c("Sympatric","Allopatric"),cex.axis=1.25) axis(2, at=seq(-1,1, by=.5),las=2,cex.axis=1.25) par(op) text(.55,.99, expression("(b)"),cex=1.25) #Sex*Pop summary(rein) fe<-fEFB$SOP_5cm me<-mEFB$SOP_5cm fl<-fLNC$SOP_5cm ml<-mLNC$SOP_5cm fm<-fMFR$SOP_5cm mm<-mMFR$SOP_5cm ff<-fFCT$SOP_5cm mf<-mFCT$SOP_5cm a<-mean(fe) adev<-sd(fe) aerr<-adev/sqrt(length(fe)) b<-mean(me) bdev<-sd(me) berr<-bdev/sqrt(length(me)) c<-mean(fl) cdev<-sd(fl) cerr<-cdev/sqrt(length(fl)) d<-mean(ml) ddev<-sd(ml) derr<-ddev/sqrt(length(ml)) e<-mean(fm) edev<-sd(fm) eerr<-edev/sqrt(length(fm)) f<-mean(mm) fdev<-sd(mm) ferr<-fdev/sqrt(length(mm)) g<-mean(ff) gdev<-sd(ff) gerr<-gdev/sqrt(length(ff)) h<-mean(mf) hdev<-sd(mf) herr<-hdev/sqrt(length(mf)) avg<-cbind(a,b,c,d,e,f,g,h) serr<-cbind(aerr,berr,cerr,derr,eerr,ferr,gerr,herr) x<-c(1,1,2,2,3,3,4,4) op<-par(mar=c(7,7,5,2)) plot(x,avg,ylim = range(c(avg-serr+.3,avg+serr-.25)),pch=19,col=c("black","grey"), xlab = "", ylab="SOP",axes=F,cex.lab=1.5,cex=1.5) box() arrows(x,avg-serr,x,avg+serr,length=0.05,angle=90,code=3,col=c("black","grey")) axis(1,at=seq(1,4,by=1),labels=c("EF","LC","MF","FC"),cex.axis=1.25) axis(2, at=seq(0,2, by=0.25),las=2,cex.axis=1.2) par(op) abline(0,0,lty=3) axis(1,at = seq(1,2,by=1),labels=c("",""),tck=0,line=1,cex=2,lwd=2) mtext("Sympatric",1,line=1.5,at=1.5,cex=1.25) axis(1,at = seq(3,4,by=1),labels=c("",""),tck=0,line=1,cex=2,lwd=2) mtext("Allopatric",1,line=1.5,at=3.5,cex=1.25) mtext("Population",1,line=3,at=2.5,cex=1.5) legend("topright",inset = 0.02, legend=c("Female","Male"),col=c("black","grey"),pch=16) #Total isolation index plot stream$context = factor(stream$context, levels = c("sym","allo")) op<-par(mar=c(5,7,4,2)) plot(stream$context,stream$Iso_index,ylim=c(0.4,1.0),axes=F) title(xlab= "Context",ylab= "Total Isolation Index",cex.lab=1.5,line=3.5) box() axis(1,at=seq(1,2,by=1),labels=c("Sympatric","Allopatric"),cex.axis=1.25) axis(2, at=seq(-1,1.0, by=0.10),las=2,cex.axis=1.25) par(op) text(1.5,.98, expression("*"),cex=2) #Zonale iso index plot streamZON<-stream[c(4:6,10:12,16:18),] symZ<-streamZON[c(4:6),] alloZ<-streamZON[-c(4:6),] fc<-alloZ[c(1:3),] mf<-alloZ[c(4:6),] c1<-mean(symZ$MF_con) c2<-mean(symZ$MF_het) sd1<-sd(symZ$MF_con) sd2<-sd(symZ$MF_het) m1<-sd1/sqrt(3) m2<-sd2/sqrt(3) d1<-mean(symZ$MM_con) d2<-mean(symZ$MM_het) dd1<-sd(symZ$MM_con) dd2<-sd(symZ$MM_het) de1<-dd1/sqrt(3) de2<-dd2/sqrt(3) e1<-mean(symZ$Un_con) e2<-mean(symZ$Un_het) ed1<-sd(symZ$Un_con) ed2<-sd(symZ$Un_het) ee1<-ed1/sqrt(3) ee2<-ed2/sqrt(3) f1<-mean(symZ$Suc_con) f2<-mean(symZ$Suc_het) fd1<-sd(symZ$Suc_con) fd2<-sd(symZ$Suc_het) fe1<-fd1/sqrt(3) fe2<-fd2/sqrt(3) g1<-mean(symZ$Sp_con) g2<-mean(symZ$Sp_het) gd1<-sd(symZ$Sp_con) gd2<-sd(symZ$Sp_het) ge1<-gd1/sqrt(3) ge2<-gd2/sqrt(3) H<-c(c1,c2) I<-c(d1,d2) J<-c(e1,e2) K<-c(f1,f2) L<-c(g1,g2) op <- par(mfrow = c(2,2),mar=c(2,4,2,2),mai = c(.1,.5, 0.2,.5 )) plot<-barplot(H,ylab="",xlab="",names=c("",""),ylim=c(0,100),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F,main="Sympatric") box() #location of err bars x=c(.7,1.9) h=c(c1,c2) lower=c(c1,c2) upper=c(c1+sd1,c2+sd2) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,cex.axis=1.25) text(1.3,95, expression(Male-female ~paste("chase",sep="")~italic("I = 0.76")),cex=1.2) plot<-barplot(Atab,ylab="", xlab="",names=c("",""),ylim=c(0,300),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F,main="Allopatric") box() #location of err bars x=c(.7,1.9) lower=c(con.mf,het.mf) upper=c(con.mf+sd.con,het.mf+sd.het) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,at=seq(0,300,by=100),cex.axis=1.25) text(1.3,285, expression(Male-female ~paste("chase",sep="")~italic("I = 0.69")),cex=1.2) plot<-barplot(I,ylab="", xlab="",names=c("",""),ylim=c(0,100),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F) box() #location of err bars x=c(.7,1.9) h=c(d1,d2) lower=c(d1,d2) upper=c(d1+dd1,d2+dd2) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,cex.axis=1.25) text(1.3,95, expression(Male-male ~paste("chase",sep="")~italic("I = 0.90")),cex=1.2) plot<-barplot(Btab,ylab="", xlab="",names=c("",""),ylim=c(0,275),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F) box() #location of err bars x=c(.7,1.9) lower=c(con.mm,het.mm) upper=c(con.mm+sd.cmm,het.mm+sd.hmm) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,cex.axis=1.25) text(1.3,260, expression(Male-male ~paste("chase",sep="")~italic("I = 0.91")),cex=1.2) op <- par(mfrow = c(2,2),mar=c(2,4,2,2)) plot<-barplot(J,ylab="", xlab="",names=c("",""),ylim=c(0,200),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F) box() #location of err bars x=c(.7,1.9) h=c(e1,e2) lower=c(e1,e2) upper=c(e1+ed1,e2+ed2) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,at=seq(0,200,by=25),cex.axis=1.25) text(1.3,190, expression(Solicit ~paste("(unsuccessful)",sep="")~italic("I = 0.94")),cex=1.2) plot<-barplot(Ctab,ylab="", xlab="",names=c("",""),ylim=c(0,300),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F) box() #location of err bars x=c(.7,1.9) lower=c(con.un,het.un) upper=c(con.un+sd.cun,het.un+sd.hun) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,cex.axis=1.25) text(1.3,285, expression(Solicit ~paste("(unsuccessful)",sep="")~italic("I = 0.64")),cex=1.2) plot<-barplot(K,ylab="", xlab="",names=c("",""),ylim=c(0,200),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F) box() #location of err bars x=c(.7,1.9) h=c(f1,f2) lower=c(f1,f2) upper=c(f1+fd1,f2+fd2) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,at=seq(0,200,by=25),cex.axis=1.25) text(1.3,190, expression(Solicit ~paste("(successful)",sep="")~italic("I = 0.96")),cex=1.2) plot<-barplot(Dtab,ylab="", xlab="",names=c("",""),ylim=c(0,250),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F) box() #location of err bars x=c(.7,1.9) lower=c(con.su,het.su) upper=c(con.su+sd.csu,het.su+sd.hsu) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,cex.axis=1.25) text(1.3,235, expression(Solicit ~paste("(successful)",sep="")~italic("I = 0.85")),cex=1.2) op <- par(mfrow = c(2,2),mar=c(2,4,2,2)) plot<-barplot(L,ylab="", xlab="",names=c("Conspecific","Heterospecific"),ylim=c(0,40),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F) box() axis(2,las=2,cex.axis=1.25) #location of err bars x=c(.7,1.9) h=c(g1,g2) lower=c(g1,g2) upper=c(g1+gd1,g2+gd2) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) text(1.3,38, expression(Spawn~italic("I = 1.00")),cex=1.2) plot<-barplot(Etab,ylab="", xlab="",names=c("Conspecific","Heterospecific"),ylim=c(0,60),cex.lab=1.5, cex.axis=1.25,cex=1.25,axes=F) box() #location of err bars x=c(.7,1.9) lower=c(con.s,het.s) upper=c(con.s+sd.cs,het.s+sd.hs) arrows(x,lower,x,upper,lwd=1.5,angle=90,code=3,length=0.05) axis(2,las=2,cex.axis=1.25) text(1.3,57, expression(Spawn~italic("I = 0.99")),cex=1.2)