install.packages("nonnest2") library(ggplot2) library(ggpmisc) library(car) library(nonnest2) dat <- read.table("C:/Users/YiMing/Desktop/MantelTest/mantel_data.txt", sep="", header=T) head(dat) # scatter plot for TMRCA my.formula <- y ~ x ggplot(dat, aes(x=distance, y=TMRCA, shape=drainage, col=drainage)) + ggtitle("Coalescence Time") + geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 11550, ymax = 21000), alpha = 0.01, fill="grey", color=NA) + geom_point() + scale_x_continuous(name ="Geographical distance (km)") + scale_y_continuous(name ="TMRCA (yr)") + scale_shape_manual(values=c(16, 1)) + geom_smooth(method=lm, se=FALSE, linetype = "dashed", size=0.5)+ scale_color_manual(values=c("black","grey50"))+ stat_poly_eq(formula = my.formula, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "right", size=3, parse = TRUE) + theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) ########################## model building and comparison ###################### # do lm for each each group newdat <- dat[with(dat, order(drainage)),] #tmrca_within.lm <- lm(TMRCA~distances, data=newdat[1:30,]) #summary(tmrca_within.lm) #tmrca_between.lm <- lm(TMRCA~distances, data=newdat[31:96,]) #summary(tmrca_between.lm) ############### model K=6 ############# # create regression with interaction term to check the interaction drainage.tmrca.lm <- lm(TMRCA~drainage*distance, data=newdat) summary(drainage.tmrca.lm) # the interactoin is not significant, move on drainage.tmrca.lm <- lm(TMRCA~drainage+distance, data=newdat) summary(drainage.tmrca.lm) drainage.tmrca.aov <- Anova(drainage.tmrca.lm, type=2) AIC(drainage.tmrca.lm) BIC(drainage.tmrca.lm) ############### model K=4 ############# ggplot(dat, aes(x=distance, y=TMRCA, shape=K4, col=K4)) + ggtitle("Coalescence Time") + geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 11550, ymax = 21000), alpha = 0.01, fill="grey", color=NA) + geom_point() + scale_x_continuous(name ="Geographical distance (km)") + scale_y_continuous(name ="TMRCA (yr)") + scale_shape_manual(values=c(16, 1)) + geom_smooth(method=lm, se=FALSE, linetype = "dashed", size=0.5)+ scale_color_manual(values=c("black","grey50"))+ stat_poly_eq(formula = my.formula, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "right", size=3, parse = TRUE) + theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) # create regression with interaction term to check the interaction K4.tmrca.lm <- lm(TMRCA~K4*distance, data=newdat) summary(K4.tmrca.lm) # the interactoin is not significant, move on K4.tmrca.lm <- lm(TMRCA~K4+distance, data=newdat) summary(K4.tmrca.lm) K4.tmrca.aov <- Anova(K4.tmrca.lm, type=2) K4.tmrca.aov AIC(K4.tmrca.lm) BIC(K4.tmrca.lm) ################ model K=5 ############ # create regression with interaction term to check the interaction K5.tmrca.lm <- lm(TMRCA~K5*distance, data=newdat) summary(K5.tmrca.lm) # the interactoin is not significant, move on K5.tmrca.lm <- lm(TMRCA~K5+distance, data=newdat) summary(K5.tmrca.lm) K5.tmrca.aov <- Anova(K5.tmrca.lm, type=2) K5.tmrca.aov AIC(K5.tmrca.lm) BIC(K5.tmrca.lm) ################ model morphotype ############# # create regression with interaction term to check the interaction morphotype.tmrca.lm <- lm(TMRCA~morphotype*distance, data=newdat) summary(morphotype.tmrca.lm) # the interactoin is not significant, move on morphotype.tmrca.lm <- lm(TMRCA~morphotype+distance, data=newdat) summary(morphotype.tmrca.lm) morphotype.tmrca.aov <- Anova(morphotype.tmrca.lm, type=2) morphotype.tmrca.aov AIC(morphotype.tmrca.lm) BIC(morphotype.tmrca.lm) ################## model comparisons ############ ## K4 vs K6 icci(morphotype.tmrca.lm, K4.tmrca.lm, conf.level = 0.95, ll1 = llcont, ll2 = llcont) vuongtest(morphotype.tmrca.lm, K4.tmrca.lm) ## K5 vs K6 icci(drainage.tmrca.lm, K5.tmrca.lm, conf.level = 0.95, ll1 = llcont, ll2 = llcont) vuongtest(drainage.tmrca.lm, K5.tmrca.lm) ## morphotype vs K4 icci(morphotype.tmrca.lm, K4.tmrca.lm, conf.level = 0.95, ll1 = llcont, ll2 = llcont) vuongtest(morphotype.tmrca.lm, K4.tmrca.lm) ## morphotype vs K5 icci(morphotype.tmrca.lm, K5.tmrca.lm, conf.level = 0.95, ll1 = llcont, ll2 = llcont) vuongtest(morphotype.tmrca.lm, K5.tmrca.lm) # scatter plot for mean Fst ggplot(dat, aes(x=distance, y=meanFST, shape=drainage, color=drainage)) + ggtitle("mean FST") + geom_point() + scale_x_continuous(name ="Geographical distance (km)") + scale_y_continuous(name ="mean FST") + scale_shape_manual(values=c(16, 1)) + geom_smooth(method=lm, se=FALSE, linetype = "dashed", size=0.5)+ scale_color_manual(values=c("black","grey50"))+ stat_poly_eq(formula = my.formula, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", size=3, parse = TRUE) + theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) # create regression with interaction term to check the interaction meanfst.lm <- lm(meanFST~K4*distance, data=newdat) summary(meanfst.lm) # the interactoin is not significant, move on meanfst.aov <- Anova(meanfst.lm, type=2) meanfst.aov ######################## For LCP ########################## dat <- read.table("C:/Users/wengz/Box/Chapter 1/MantelTest/mantel_data2.txt", sep="", header=T) head(dat) ## scatter plot for FST my.formula <- y ~ x ggplot(dat, aes(x=LCP, y=meanFST, shape=drainage, color=drainage)) + ggtitle("mean FST") + geom_point() + scale_x_continuous(name ="Least Cost Path (km)") + scale_y_continuous(name ="mean FST") + scale_shape_manual(values=c(16, 1)) + geom_smooth(method=lm, se=FALSE, linetype = "dashed", size=0.5)+ scale_color_manual(values=c("black","grey50"))+ stat_poly_eq(formula = my.formula, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", size=3, parse = TRUE) + theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) ## The outlier of x will strongly affect the result of ancova, remove them dat <- read.table("C:/Users/wengz/Box/Chapter 1/MantelTest/mantel_data3.txt", sep="", header=T) head(dat) ggplot(dat, aes(x=LCP, y=meanFST, shape=drainage, color=drainage)) + ggtitle("mean FST") + geom_point() + scale_x_continuous(name ="Least Cost Path (km)") + scale_y_continuous(name ="mean FST") + scale_shape_manual(values=c(16, 1)) + geom_smooth(method=lm, se=FALSE, linetype = "dashed", size=0.5)+ scale_color_manual(values=c("black","grey50"))+ stat_poly_eq(formula = my.formula, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", size=3, parse = TRUE) + theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) # scatter plot for TMRCA my.formula <- y ~ x ggplot(dat, aes(x=distance, y=TMRCA, shape=drainage, col=drainage)) + ggtitle("Coalescence Time") + geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 11550, ymax = 21000), alpha = 0.01, fill="grey", color=NA) + geom_point() + scale_x_continuous(name ="Geographical distance (km)") + scale_y_continuous(name ="TMRCA (yr)") + scale_shape_manual(values=c(16, 1)) + geom_smooth(method=lm, se=FALSE, linetype = "dashed", size=0.5)+ scale_color_manual(values=c("black","grey50"))+ stat_poly_eq(formula = my.formula, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "right", size=3, parse = TRUE) + theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) ########################## model building and comparison ###################### # do lm for each each group newdat <- dat[with(dat, order(drainage)),] tmrca_within.lm <- lm(TMRCA~LCP, data=newdat[1:57,]) summary(tmrca_within.lm) tmrca_between.lm <- lm(TMRCA~LCP, data=newdat[58:87,]) summary(tmrca_between.lm) ############### model K=6 ############# # create regression with interaction term to check the interaction drainage.tmrca.lm <- lm(TMRCA~drainage*LCP, data=newdat) summary(drainage.tmrca.lm) # the interactoin is not significant, move on drainage.tmrca.lm <- lm(TMRCA~drainage+LCP, data=newdat) summary(drainage.tmrca.lm) drainage.tmrca.aov <- Anova(drainage.tmrca.lm, type=2) AIC(drainage.tmrca.lm) BIC(drainage.tmrca.lm) ############### model K=4 ############# ggplot(dat, aes(x=LCP, y=TMRCA, shape=K4, col=K4)) + ggtitle("Coalescence Time") + geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 11550, ymax = 21000), alpha = 0.01, fill="grey", color=NA) + geom_point() + scale_x_continuous(name ="Geographical distance (km)") + scale_y_continuous(name ="TMRCA (yr)") + scale_shape_manual(values=c(16, 1)) + geom_smooth(method=lm, se=FALSE, linetype = "dashed", size=0.5)+ scale_color_manual(values=c("black","grey50"))+ stat_poly_eq(formula = my.formula, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "right", size=3, parse = TRUE) + theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) # create regression with interaction term to check the interaction K4.tmrca.lm <- lm(TMRCA~K4*LCP, data=newdat) summary(K4.tmrca.lm) # the interactoin is not significant, move on K4.tmrca.lm <- lm(TMRCA~K4+LCP, data=newdat) summary(K4.tmrca.lm) K4.tmrca.aov <- Anova(K4.tmrca.lm, type=2) K4.tmrca.aov AIC(K4.tmrca.lm) BIC(K4.tmrca.lm) ################ model K=5 ############ # create regression with interaction term to check the interaction K5.tmrca.lm <- lm(TMRCA~K5*LCP, data=newdat) summary(K5.tmrca.lm) # the interactoin is not significant, move on K5.tmrca.lm <- lm(TMRCA~K5+LCP, data=newdat) summary(K5.tmrca.lm) K5.tmrca.aov <- Anova(K5.tmrca.lm, type=2) K5.tmrca.aov AIC(K5.tmrca.lm) BIC(K5.tmrca.lm) ################ model morphotype ############# # create regression with interaction term to check the interaction morphotype.tmrca.lm <- lm(TMRCA~morphotype*LCP, data=newdat) summary(morphotype.tmrca.lm) # the interactoin is not significant, move on morphotype.tmrca.lm <- lm(TMRCA~morphotype+LCP, data=newdat) summary(morphotype.tmrca.lm) morphotype.tmrca.aov <- Anova(morphotype.tmrca.lm, type=2) morphotype.tmrca.aov AIC(morphotype.tmrca.lm) BIC(morphotype.tmrca.lm) ################## model comparisons ############ ## K4 vs K6 icci(drainage.tmrca.lm, K4.tmrca.lm, conf.level = 0.95, ll1 = llcont, ll2 = llcont) vuongtest(drainage.tmrca.lm, K4.tmrca.lm) ## K5 vs K6 icci(drainage.tmrca.lm, K5.tmrca.lm, conf.level = 0.95, ll1 = llcont, ll2 = llcont) vuongtest(drainage.tmrca.lm, K5.tmrca.lm) ## morphotype vs K4 icci(morphotype.tmrca.lm, K4.tmrca.lm, conf.level = 0.95, ll1 = llcont, ll2 = llcont) vuongtest(morphotype.tmrca.lm, K4.tmrca.lm) ## morphotype vs K5 icci(morphotype.tmrca.lm, K5.tmrca.lm, conf.level = 0.95, ll1 = llcont, ll2 = llcont) vuongtest(morphotype.tmrca.lm, K5.tmrca.lm) # scatter plot for mean Fst ggplot(dat, aes(x=LCD, y=meanFST, shape=drainage, color=drainage)) + ggtitle("mean FST") + geom_point() + scale_x_continuous(name ="Geographical distance (km)") + scale_y_continuous(name ="mean FST") + scale_shape_manual(values=c(16, 1)) + geom_smooth(method=lm, se=FALSE, linetype = "dashed", size=0.5)+ scale_color_manual(values=c("black","grey50"))+ stat_poly_eq(formula = my.formula, aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), label.x = "left", size=3, parse = TRUE) + theme_classic()+ theme(plot.title = element_text(hjust = 0.5)) # create regression with interaction term to check the interaction meanfst.lm <- lm(meanFST~K4*distance, data=newdat) summary(meanfst.lm) # the interactoin is not significant, move on meanfst.aov <- Anova(meanfst.lm, type=2) meanfst.aov