#SDavies_DistanceAnalysis # Oct 2013 -etreml@unimelb.edu.au #Oct 30, SD edits #Jan 30, 2014 #reformatted for submission Oct 2014 library(ecodist) library(ggplot2) setwd("/Users/Ted/Desktop/Mol_ecol_data") #PLD45, 65, 90 modify the data you load in each time. dFst_OW<-read.table(file='Dig_Island_Fst_BelowDiag_2014.csv', row.names=1, header=TRUE, sep=',') hFst_OW<-read.table(file='Hya_Island_Fst_BelowDiag_2014.csv', row.names=1, header=TRUE, sep=',') load("PLD90_M_DispersalDistance.RData") Ddist_init<-SubMat EW<-SubMat[lower.tri(SubMat)] SubMat.t<-t(SubMat) WE<-SubMat.t[lower.tri(SubMat.t)] Ddist.mean<-rowMeans(cbind(EW,WE)) labs<-read.table(file='Island_Labels2014.csv', row.names=1, header=TRUE, sep=',') Labs.l<-labs[lower.tri(labs)] df90=data.frame(ew90=EW, we90=WE) df65=data.frame(ew65=EW, we65=WE) df<-data.frame(ew45=EW, we45=WE) sum=data.frame(df90$ew90, df90$we90, df65$ew65, df65$we65, df$ew45, df$we45) sum$rel90=sum$df90.ew90-sum$df90.we90 sum$rel45=sum$df.ew45-sum$df.we45 sum$rel65=sum$df65.ew65-sum$df65.we65 sum summary(sum) sump=subset(sum[,7:9]) write.csv(sump,file="rel_ew-we.csv") #modified csv in excel sump2<-read.table(file='rel_ew-we.csv',header=TRUE, sep=',') sump2$X=as.factor(sump2$X) summary(sump2) summary(df) #Make the figure for east to west at PLD=45 ggplot(df, aes(-(we45/log(10)), -(ew45/log(10)))) + geom_point(data=df, pch=21, size=2, col="black", fill="grey")+ # geom_smooth(method="lm")+ geom_text(aes(label=Labs.l), size=2.5,vjust = -1)+ ylab("log10(Pmig[East--> West])")+ xlab("log10(Pmig[West--> East])")+ xlim(-11,0)+ ylim(-11,0)+ theme_bw() + geom_abline(intercept=0, slope=1, linetype=3) summary(df65) #Make the figure for east to west at PLD=65 ggplot(df65, aes(-(we65/log(10)), -(ew65/log(10)))) + geom_point(data=df65, pch=21, size=2, col="black", fill="grey")+ # geom_smooth(method="lm")+ geom_text(aes(label=Labs.l), size=2.5,vjust = -1)+ ylab("log10(Pmig[East--> West])")+ xlab("log10(Pmig[West--> East])")+ xlim(-9.5,-0.5)+ ylim(-9.5,-0.5)+ theme_bw() + geom_abline(intercept=0, slope=1, linetype=3) summary(df90) #Make the figure for east to west at PLD=90 --> shows migration direction is equal in both directions for most islands with the exception of Guam and Kosrae. ggplot(df90, aes(-(we90/log(10)), -(ew90/log(10)))) + geom_point(data=df90, pch=21, size=2, col="black", fill="grey")+ # geom_smooth(method="lm")+ geom_text(aes(label=Labs.l), size=2.5,vjust = -1)+ ylab("log10(Pmig[East--> West])")+ xlab("log10(Pmig[West--> East])")+ xlim(-8.5,0)+ ylim(-8.5,0)+ theme_bw() + geom_abline(intercept=0, slope=1, linetype=3) #This figure will go into manuscript. ######Let's see what PLD min, max, mean works best for each species #Now let's look at A. hyacinthus Fst_OW<-read.table(file='Hya_Island_Fst_BelowDiag_2014.csv', row.names=1, header=TRUE, sep=',') load("PLD45_M_DispersalDistance.RData") Ddist_init<-SubMat Labs<-read.table(file='Island_Labels2014.csv', row.names=1, header=TRUE, sep=',') Sp = 'Hya45' #Extract as vectors #Fst/OW data Fst.v<-Fst_OW[lower.tri(Fst_OW)] Rous.Fst<-Fst.v/(1-Fst.v) Fst_OW.t<-t(Fst_OW) OWD.v<-Fst_OW.t[lower.tri(Fst_OW.t)] EW<-SubMat[lower.tri(SubMat)] SubMat.t<-t(SubMat) WE<-SubMat.t[lower.tri(SubMat.t)] H45mean<-rowMeans(cbind(EW,WE)) H45min<-apply(cbind(EW,WE), 1, min) H45max<-apply(cbind(EW,WE), 1, max) mantel.fit<-ecodist::mantel(Rous.Fst~OWD.v, nperm=10000) mod=(lm(Rous.Fst~OWD.v)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) # (Intercept) OWD.v # 1.186599e-02 1.934616e-05 mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # Hya45IBD mantel test: r = 0.74, pval = 0.047" #plot & save Fst vs OW df<-data.frame(ew=EW, we=WE, owd=OWD.v, Rous=Rous.Fst) labs<-read.table(file='Island_Labels2014.csv', row.names=1, header=TRUE, sep=',') Labs.l<-labs[lower.tri(labs)] ggplot(df, aes(owd, Rous)) + geom_point()+ geom_smooth(method="lm")+ geom_text(aes(label=Labs.l), size=3,vjust = -1)+ ylab("All Pair-wise FST/(1-FST)")+ xlab("Over Water Distance (km)")+ theme_bw() + ggtitle(lm.title) #now with dispersal matrix mantel.fit<-ecodist::mantel(Rous.Fst~D45mean, nperm=10000) mod=(lm(Rous.Fst~D45mean)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # "Hya45mean mantel test: r = 0.799, pval = 0.023" mantel.fit<-ecodist::mantel(Rous.Fst~D45min, nperm=10000) mod=(lm(Rous.Fst~D45min)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Hya45min mantel test: r = 0.865, pval = 0.016" mantel.fit<-ecodist::mantel(Rous.Fst~D45max, nperm=10000) mod=(lm(Rous.Fst~D45max)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Hya45 mantel test: r = 0.696, pval = 0.037" ####Now hya for 65 load("PLD65_M_DispersalDistance.RData") Ddist_init<-SubMat Sp = 'Hya65' #Extract as vectors EW<-SubMat[lower.tri(SubMat)] SubMat.t<-t(SubMat) WE<-SubMat.t[lower.tri(SubMat.t)] H65mean<-rowMeans(cbind(EW,WE)) H65min<-apply(cbind(EW,WE), 1, min) H65max<-apply(cbind(EW,WE), 1, max) mantel.fit<-ecodist::mantel(Rous.Fst~H65mean, nperm=10000) mod=(lm(Rous.Fst~H65mean)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Hya65 mantel test: r = 0.843, pval = 0.01" mantel.fit<-ecodist::mantel(Rous.Fst~H65min, nperm=10000) mod=(lm(Rous.Fst~H65min)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Hya65 mantel test: r = 0.886, pval = 0.014" #This one won so plot it df<-data.frame(ew=EW, we=WE, owd=OWD.v, Rous=Rous.Fst, hyamin65=H65min) labs<-read.table(file='Island_Labels2014.csv', row.names=1, header=TRUE, sep=',') Labs.l<-labs[lower.tri(labs)] ggplot(df, aes(hyamin65, Rous)) + geom_point()+ geom_smooth(method="lm")+ geom_text(aes(label=Labs.l), size=3,vjust = -1)+ ylab("All Pair-wise FST/(1-FST)")+ xlab("Minimum Modeled Distance (PLD=65)")+ theme_bw() + ggtitle(lm.title) mantel.fit<-ecodist::mantel(Rous.Fst~H65max, nperm=10000) mod=(lm(Rous.Fst~H65max)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Hya65 mantel test: r = 0.728, pval = 0.032" #####Now for Hya 90 load("PLD90_M_DispersalDistance.RData") Ddist_init<-SubMat Sp = 'Hya90' #Extract as vectors EW<-SubMat[lower.tri(SubMat)] SubMat.t<-t(SubMat) WE<-SubMat.t[lower.tri(SubMat.t)] H90mean<-rowMeans(cbind(EW,WE)) H90min<-apply(cbind(EW,WE), 1, min) H90max<-apply(cbind(EW,WE), 1, max) mantel.fit<-ecodist::mantel(Rous.Fst~H90mean, nperm=10000) mod=(lm(Rous.Fst~D90mean)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Hya90 mantel test: r = 0.859, pval = 0.009" mantel.fit<-ecodist::mantel(Rous.Fst~H90min, nperm=10000) mod=(lm(Rous.Fst~H90min)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Hya90 mantel test: r = 0.867, pval = 0.01" mantel.fit<-ecodist::mantel(Rous.Fst~H90max, nperm=10000) mod=(lm(Rous.Fst~H90max)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Hya90 mantel test: r = 0.753, pval = 0.034" dfhya<-data.frame(ew=EW, we=WE, owd=OWD.v, Rous=Rous.Fst, H45min=H45min, H45mean=H45mean, H45max=H45max, H65min=H65min, H65mean=H65mean, H65max=H65max, H90min=H90min, H90mean=H90mean, H90max=H90max, lab=Labs.l) write.csv(dfhya,file="dfhya.csv") #####NOW A. digitifera Fst_OW<-read.table(file='Dig_Island_Fst_BelowDiag_2014.csv', row.names=1, header=TRUE, sep=',') load("PLD65_M_DispersalDistance.RData") Ddist_init<-SubMat Sp = 'Dig65' Fst.v<-Fst_OW[lower.tri(Fst_OW)] Rous.Fst<-Fst.v/(1-Fst.v) Fst_OW.t<-t(Fst_OW) OWD.v<-Fst_OW.t[lower.tri(Fst_OW.t)] #Extract as vectors EW<-SubMat[lower.tri(SubMat)] SubMat.t<-t(SubMat) WE<-SubMat.t[lower.tri(SubMat.t)] D65mean<-rowMeans(cbind(EW,WE)) D65min<-apply(cbind(EW,WE), 1, min) D65max<-apply(cbind(EW,WE), 1, max) mantel.fit<-ecodist::mantel(Rous.Fst~OWD.v, nperm=10000) mod=(lm(Rous.Fst~OWD.v)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) # (Intercept) OWD.v # 1.153210e-02 7.088037e-06 mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Dig65 mantel test: r = 0.616, pval = 0.007" df<-data.frame(ew=EW, we=WE, owd=OWD.v, Rous=Rous.Fst) labs<-read.table(file='Island_Labels2014.csv', row.names=1, header=TRUE, sep=',') Labs.l<-labs[lower.tri(labs)] ggplot(df, aes(owd, Rous)) + geom_point()+ geom_smooth(method="lm")+ geom_text(aes(label=Labs.l), size=3,vjust = -1)+ ylab("All Pair-wise FST/(1-FST)")+ xlab("Over Water Distance (km)")+ theme_bw() + ggtitle(lm.title) mantel.fit<-ecodist::mantel(Rous.Fst~D65mean, nperm=10000) mod=(lm(Rous.Fst~D65mean)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # Dig65mean mantel test: r = 0.736, pval = 0.001 mantel.fit<-ecodist::mantel(Rous.Fst~D65min, nperm=10000) mod=(lm(Rous.Fst~D65min)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # "Dig65min mantel test: r = 0.479, pval = 0.042" mantel.fit<-ecodist::mantel(Rous.Fst~D65max, nperm=10000) mod=(lm(Rous.Fst~D65max)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # "Dig65 mantel test: r = 0.822, pval = 0" #now for A. digi 45 load("PLD45_M_DispersalDistance.RData") Ddist_init<-SubMat Sp = 'Dig45' EW<-SubMat[lower.tri(SubMat)] SubMat.t<-t(SubMat) WE<-SubMat.t[lower.tri(SubMat.t)] D45mean<-rowMeans(cbind(EW,WE)) D45min<-apply(cbind(EW,WE), 1, min) D45max<-apply(cbind(EW,WE), 1, max) mantel.fit<-ecodist::mantel(Rous.Fst~D45mean, nperm=10000) mod=(lm(Rous.Fst~D45mean)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Dig45mean mantel test: r = 0.712, pval = 0.001" mantel.fit<-ecodist::mantel(Rous.Fst~D45min, nperm=10000) mod=(lm(Rous.Fst~D45min)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Dig45min mantel test: r = 0.486, pval = 0.027" mantel.fit<-ecodist::mantel(Rous.Fst~D45max, nperm=10000) mod=(lm(Rous.Fst~D45max)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Dig45 mantel test: r = 0.803, pval = 0" #now for A. digi 90 load("PLD90_M_DispersalDistance.RData") Ddist_init<-SubMat Sp = 'Dig90' EW<-SubMat[lower.tri(SubMat)] SubMat.t<-t(SubMat) WE<-SubMat.t[lower.tri(SubMat.t)] D90mean<-rowMeans(cbind(EW,WE)) D90min<-apply(cbind(EW,WE), 1, min) D90max<-apply(cbind(EW,WE), 1, max) mantel.fit<-ecodist::mantel(Rous.Fst~D90mean, nperm=10000) mod=(lm(Rous.Fst~D90mean)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Dig90 mantel test: r = 0.708, pval = 0.003" mantel.fit<-ecodist::mantel(Rous.Fst~D90min, nperm=10000) mod=(lm(Rous.Fst~D90min)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Dig90 mantel test: r = 0.41, pval = 0.084" mantel.fit<-ecodist::mantel(Rous.Fst~D90max, nperm=10000) mod=(lm(Rous.Fst~D90max)) #Classic R2, ignoring autocorrelation summary(mod)   coef(mod) mantel.r<- round(mantel.fit[1], 3) mantel.pval<- round(mantel.fit[2], 3) lm.title <- paste(Sp, " mantel test: r = ",mantel.r,", pval = ",mantel.pval, sep="") lm.title # [1] "Dig90max mantel test: r = 0.826, pval = 0" df<-data.frame(ew=EW, we=WE, owd=OWD.v, Rous=Rous.Fst, d90max=D90max) labs<-read.table(file='Island_Labels2014.csv', row.names=1, header=TRUE, sep=',') Labs.l<-labs[lower.tri(labs)] ggplot(df, aes(d90max, Rous)) + geom_point()+ geom_smooth(method="lm")+ geom_text(aes(label=Labs.l), size=3,vjust = -1)+ ylab("All Pair-wise FST/(1-FST)")+ xlab("Maximum Modeled Distance (PLD=90)")+ theme_bw() + ggtitle(lm.title) #A. digi WINNER [1] "Dig90max mantel test: r = 0.826, pval = 0" # A. hya WINNER [1] "Hya65min minmantel test: r = 0.886, pval = 0.014" dfdig<-data.frame(ew=EW, we=WE, owd=OWD.v, Rous=Rous.Fst, D45min=D45min, D45mean=D45mean, D45max=D45max, D65min=D65min, D65mean=D65mean, D65max=D65max, D90min=D90min, D90mean=D90mean, D90max=D90max, lab=Labs.l) write.csv(dfdig,file="dfdig.csv") #making figures for the manuscript tot<-read.table(file='summary_TremlModel.csv', header=TRUE, sep=',') str(tot) #FST~OWD ggplot(tot, aes(owd, Rous, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(colour=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=lab), size=3,vjust = -1)+ ylab("All Pair-wise FST/(1-FST)")+ xlab("Over Water Distance (km)")+ theme_bw() # FST~90max digi's best ggplot(tot, aes(X90max, Rous)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(colour=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=lab), size=3,vjust = -1)+ ylab("All Pair-wise FST/(1-FST)")+ xlab("Maximum Modeled Distance")+ theme_bw() # FST~65min hya's best ggplot(tot, aes(X65min, Rous)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(colour=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=lab), size=3,vjust = -1)+ ylab("All Pair-wise FST/(1-FST)")+ xlab("Minimum Modeled Distance")+ theme_bw() #diversity decreases as you go west #Shannon~OWD div<-read.table(file='Shan_NA.csv', header=TRUE, sep=',') str(div) ggplot(div, aes(OWD, Sha, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(colour=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=pop), size=3,vjust = -1)+ ylab("Mean Shannon (sHua) ")+ xlab("Relative Over Water Distance from Palau (km)")+ xlim(0,4200)+ theme_bw() #Shannon~Ddistbestfits (hya=65min, digi=90max) ggplot(div, aes(DDfromPalau, Sha, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(color=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=pop), size=3,vjust = -1)+ ylab("Mean Shannon (sHua) ")+ xlab("Best Modeled Distance")+ theme_bw() #Ho~OWD ggplot(div, aes(OWD, Ho, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(color=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=pop), size=3,vjust = -1)+ ylab("Observed Heterozygosity")+ xlab("Relative Over Water Distance from Palau (km)")+ xlim(0,4200)+ theme_bw() ggplot(div, aes(DDfromPalau, Ho, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(color=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=pop), size=3,vjust = -1)+ ylab("Observed Heterozygosity")+ xlab("Best Modeled Distance")+ theme_bw() ###Na~OWD and modeled distance ggplot(div, aes(OWD, Na, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(color=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=pop), size=3,vjust = -1)+ ylab("Mean Number of Alleles")+ xlab("Relative Over Water Distance from Palau (km)")+ xlim(0,4200)+ theme_bw() ggplot(div, aes(DDfromPalau, Na, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(color=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=pop), size=3,vjust = -1)+ ylab("Mean Number of Alleles")+ xlab("Best Modeled Distance")+ theme_bw() ###Private~OWD and modeled distance ggplot(div, aes(OWD, priv, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(color=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=pop), size=3,vjust = -1)+ ylab("Mean Number of Private Alleles")+ xlab("Relative Over Water Distance from Palau (km)")+ xlim(0,4200)+ theme_bw() ggplot(div, aes(DDfromPalau, priv, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(color=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=pop), size=3,vjust = -1)+ ylab("Mean Number of Private Alleles")+ xlab("Best Modeled Distance")+ theme_bw() #models for A. hya hya=subset(div, spp=="h") lm1=lm(Ho~OWD, data=hya) summary(lm1) # Multiple R-squared: 0.3876, Adjusted R-squared: 0.2651 # F-statistic: 3.165 on 1 and 5 DF, p-value: 0.1354 lm2=lm(Ho~DDfromPalau, data=hya) summary(lm2) # Multiple R-squared: 0.758, Adjusted R-squared: 0.6975 # F-statistic: 12.53 on 1 and 4 DF, p-value: 0.02402 lm3=lm(Sha~OWD, data=hya) summary(lm3) coef(lm3) # Multiple R-squared: 0.8244, Adjusted R-squared: 0.7892 # F-statistic: 23.47 on 1 and 5 DF, p-value: 0.004696 # (Intercept) OWD # 1.694881e+00 -5.195732e-05 lm4=lm(Sha~DDfromPalau, data=hya) summary(lm4) coef(lm4) # Multiple R-squared: 0.889, Adjusted R-squared: 0.8613 # F-statistic: 32.05 on 1 and 4 DF, p-value: 0.0048 lm5=lm(Na~OWD, data=hya) summary(lm5) # Multiple R-squared: 0.2405, Adjusted R-squared: 0.08864 # F-statistic: 1.584 on 1 and 5 DF, p-value: 0.2638 lm6=lm(Na~DDfromPalau, data=hya) summary(lm6) # Multiple R-squared: 0.01494, Adjusted R-squared: -0.2313 # F-statistic: 0.06066 on 1 and 4 DF, p-value: 0.8176 lm7=lm(priv~OWD, data=hya) summary(lm7) # Multiple R-squared: 0.6065, Adjusted R-squared: 0.5502 # F-statistic: 10.79 on 1 and 7 DF, p-value: 0.01341 lm7=lm(priv~DDfromPalau, data=hya) summary(lm7) # Multiple R-squared: 0.7687, Adjusted R-squared: 0.7224 # F-statistic: 16.62 on 1 and 5 DF, p-value: 0.009576 #models for A. digi dig=subset(div, spp=="d") lm1d=lm(Ho~OWD, data=dig) summary(lm1d) # Multiple R-squared: 0.7131, Adjusted R-squared: 0.6721 # F-statistic: 17.4 on 1 and 7 DF, p-value: 0.004182 lm2d=lm(Ho~DDfromPalau, data=dig) summary(lm2d) # Multiple R-squared: 0.5308, Adjusted R-squared: 0.4369 # F-statistic: 5.656 on 1 and 5 DF, p-value: 0.06331 lm3d=lm(Sha~OWD, data=dig) summary(lm3d) coef(lm3d) # Multiple R-squared: 0.8413, Adjusted R-squared: 0.8186 # F-statistic: 37.11 on 1 and 7 DF, p-value: 0.000495 # (Intercept) OWD # 1.660185e+00 -7.621157e-05 lm4d=lm(Sha~DDfromPalau, data=dig) summary(lm4d) coef(lm4d) # Multiple R-squared: 0.9595, Adjusted R-squared: 0.9513 # F-statistic: 118.3 on 1 and 5 DF, p-value: 0.000114 lm5d=lm(Na~OWD, data=dig) summary(lm5d) # Multiple R-squared: 0.6505, Adjusted R-squared: 0.6006 # F-statistic: 13.03 on 1 and 7 DF, p-value: 0.008627 lm6d=lm(Na~DDfromPalau, data=dig) summary(lm6d) # Multiple R-squared: 0.7057, Adjusted R-squared: 0.6469 # F-statistic: 11.99 on 1 and 5 DF, p-value: 0.01799 lm7d=lm(priv~OWD, data=dig) summary(lm7d) # Multiple R-squared: 0.01066, Adjusted R-squared: -0.1542 # F-statistic: 0.06467 on 1 and 6 DF, p-value: 0.8078 lm7d=lm(priv~DDfromPalau, data=dig) summary(lm7d) # Multiple R-squared: 0.000616, Adjusted R-squared: -0.2492 # F-statistic: 0.002466 on 1 and 4 DF, p-value: 0.9628 #Principle Coordinate Analysis #A. hya library(ggplot2) pc1=c(-0.135952679, -0.085356825, -0.07141665, 0.119450742, 0.134881563, 0.07095831, -0.032564462) pc2=c(0.051826492, 0.029989545, 0.020115169, 0.023476897, 0.01304221, 0.01903228, -0.157482593) names=c("PAL", "NGU", "YAP", "CHU", "POH", "KOS", "MAJ") dfpc<-data.frame(PC1=pc1, PC2=pc2, names=names) dfpc$names <- factor(dfpc$names, levels=c("PAL", "NGU", "YAP", "CHU", "POH", "KOS", "MAJ")) ggplot(dfpc, aes(PC1, PC2, group=names)) + geom_vline(xintercept = 0) + geom_hline(yintercept = 0) + xlim(-0.16,0.15)+ ylim(-0.18,0.08)+ geom_point(aes(colour=factor(names), fill = factor(names)), shape=21, size = 5) + scale_fill_manual(values=c("red", "orange", "yellow", "blue", "purple", "brown", "pink")) + # scale_colour_manual(values=c("red", "orange", "yellow", "blue", "purple", "brown", "pink"))+ scale_colour_manual(values=c("black", "black", "black","black","black","black","black"))+ geom_text(aes(label=names), size=4,vjust = -1)+ ylab("PC1")+ xlab("PC2")+ theme_bw() #A. dig dpc1=c(-0.077343578,-0.040608109, -0.053604025, -0.018665653, -0.002668865, 0.055508183, 0.042744528, 0.042973173, 0.051664345) dpc2=c(-0.036, -0.019, 0.094247947, -0.020069416, -0.040264543, 0.001770318, -0.013004983, 0.0240937, 0.008206543) dnames=c("PAL", "YAP", "GUA", "CHU", "POH", "KOS", "KWA", "MAJ", "ARN") dfpcd<-data.frame(PC1=dpc1, PC2=dpc2, names=dnames) dfpcd$names <- factor(dfpcd$names, levels=c("PAL", "YAP", "GUA", "CHU", "POH", "KOS", "KWA", "MAJ", "ARN")) ggplot(dfpcd, aes(PC1, PC2, group=names)) + geom_vline(xintercept = 0) + geom_hline(yintercept = 0) + xlim(-0.08,0.06)+ ylim(-0.07, 0.11)+ geom_point(aes(colour=factor(names), fill = factor(names)), shape=21, size = 5) + scale_fill_manual(values=c("red", "yellow", "green", "blue", "purple", "brown", "darkturquoise", "pink", "grey")) + scale_colour_manual(values=c("black", "black", "black","black","black","black","black","black","black"))+ geom_text(aes(label=names), size=4,vjust = -1)+ ylab("PC1")+ xlab("PC2")+ theme_bw() #determine the slope and significance across space d=subset(p, spp=="d") lm1=lm(priv~dist, data=d) summary(lm1) coef(lm1) # Residuals: # Min 1Q Median 3Q Max # -0.26744 -0.20061 0.00621 0.08728 0.46519 # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 8.078e-01 1.664e-01 4.854 0.00185 ** # dist -1.985e-04 6.045e-05 -3.284 0.01341 * # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # Residual standard error: 0.2637 on 7 degrees of freedom # (2 observations deleted due to missingness) # Multiple R-squared: 0.6065, Adjusted R-squared: 0.5502 # F-statistic: 10.79 on 1 and 7 DF, p-value: 0.01341 # (Intercept) dist # 0.8078093830 -0.0001985401 h=subset(p, spp=="h") lm2=lm(priv~dist, data=h) summary(lm2) # Residuals: # Min 1Q Median 3Q Max # -0.21746 -0.05821 0.02904 0.10550 0.11976 # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 3.163e-01 7.331e-02 4.314 0.00501 ** # dist -6.083e-06 2.392e-05 -0.254 0.80775 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # Residual standard error: 0.135 on 6 degrees of freedom # (3 observations deleted due to missingness) # Multiple R-squared: 0.01066, Adjusted R-squared: -0.1542 # F-statistic: 0.06467 on 1 and 6 DF, p-value: 0.8078 coef(lm2) # (Intercept) dist # 3.162761e-01 -6.083220e-06 plot(priv~spp, data=p) t.test(priv~spp,equal.var=TRUE, data=p) ####Comparing the slopes of digi hya IBD dat<-read.csv(file='summary_data.csv') summary(dat) aov0=aov(Rous~OWD*spp, data=dat) anova(aov0,test="LRT") aov1=aov(Shan~OWD*spp, data=dat) anova(aov1,test="LRT") ##Adding in Jost's D jost<-read.csv(file='summary_data_October2014.csv') summary(jost) ##Jost's OWD ggplot(jost, aes(OWD, Jost.d.D, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(colour=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=Labs.min), size=3,vjust = -1)+ ylab("All Pair-wise Jost's D")+ xlab("Over Water Distance (km)")+ theme_bw() ##Jost's Dispersal model ggplot(jost, aes(Ddist.mean, Jost.d.D, group=spp)) + geom_point(aes(pch=spp), size=4.5, colour="black")+ geom_point(aes(colour=spp, pch=spp), size=4)+ geom_smooth(aes(fill=factor(spp), linetype=spp),method="lm", color="black")+ geom_text(aes(label=Labs.min), size=3,vjust = -1)+ ylab("All Pair-wise Jost's D")+ xlab("Mean Modeled Distance (90 day)")+ theme_bw() ###Jost differences between species with distance and modelled distance aov2=aov(Jost.d.D~OWD*spp, data=jost) anova(aov2,test="LRT") # Response: Jost.d.D # Df Sum Sq Mean Sq F value Pr(>F) # OWD 1 0.017801 0.017801 23.4187 1.166e-05 *** # spp 1 0.045441 0.045441 59.7812 3.005e-10 *** # OWD:spp 1 0.002123 0.002123 2.7929 0.1006 # Residuals 53 0.040287 0.000760 # A. hya hya=subset(jost, spp=="h") lmhya=lm(OWD~Jost.d.D, data=hya) summary(lmhya) # Residual standard error: 1031 on 19 degrees of freedom # Multiple R-squared: 0.2896, Adjusted R-squared: 0.2523 # F-statistic: 7.747 on 1 and 19 DF, p-value: 0.01185 tremlhya=lm(Ddist.mean~Jost.d.D, data=hya) summary(tremlhya) # Residual standard error: 1.458 on 13 degrees of freedom # (6 observations deleted due to missingness) # Multiple R-squared: 0.7293, Adjusted R-squared: 0.7085 # F-statistic: 35.02 on 1 and 13 DF, p-value: 5.084e-05 # A. digi dig=subset(jost, spp=="d") lmdig=lm(OWD~Jost.d.D, data=dig) summary(lmdig) # Residual standard error: 852.3 on 34 degrees of freedom # Multiple R-squared: 0.4846, Adjusted R-squared: 0.4694 # F-statistic: 31.97 on 1 and 34 DF, p-value: 2.425e-06 tremldig=lm(Ddist.mean~Jost.d.D, data=dig) summary(tremldig) # Residual standard error: 2.595 on 19 degrees of freedom # (15 observations deleted due to missingness) # Multiple R-squared: 0.5003, Adjusted R-squared: 0.474 # F-statistic: 19.02 on 1 and 19 DF, p-value: 0.0003358 aov3=aov(Jost.d.D~Ddist.mean*spp, data=jost) anova(aov3,test="LRT") # Response: Jost.d.D # Df Sum Sq Mean Sq F value Pr(>F) # Ddist.mean 1 0.003838 0.003838 8.1178 0.007603 ** # spp 1 0.039427 0.039427 83.4025 1.989e-10 *** # Ddist.mean:spp 1 0.013828 0.013828 29.2510 6.053e-06 *** # Residuals 32 0.015127 0.000473 #####Comparing frequency of most common allele spp<-read.csv(file='spp_alleles.csv', header=TRUE, sep=',', row.names=1) head(spp) plot(digi~hya, data=spp) lm1=lm(digi~hya, data=spp) summary(lm1) ggplot(spp, aes(digi, hya)) + geom_point(aes(), size=4.5, colour="black")+ geom_smooth(aes(),method="lm", color="black")+ ylab("Freq A. hyacinthus")+ xlab("Freq A. digitifera")+ theme_bw()