######################################### ### TV field results 2018 ### ### Healthy Prairies II – Activity 2 ### ### Jan. 2019 ### ### Adrien Pozzi ### ######################################### ################################## ### Setting R and loading data ### ################################## ### R session info ### sessionInfo() # R version 3.5.1 (2018-07-02) # Platform: x86_64-apple-darwin15.6.0 (64-bit) # Running under: macOS Sierra 10.12.6 RStudio.Version() ### setting working directory ### setwd("~/Documents/Adrien Pozzi/pro/UMN postdoc 2018-2020/experiments/TV and LB field xp 2018") ### loading the data environment ### load("~/Documents/Adrien Pozzi/pro/UMN postdoc 2018-2020/experiments/TV and LB field xp 2018/TV field xp 2018 workspace.RData") ### save image of R environment ### save.image("~/Documents/pro/UMN postdoc 2018-2020/experiments/TV and LB field xp 2018/TV field xp 2018 workspace.RData") ### packages/library used ### pack<-c("corrplot","ggfortify","vegan","lme4","lmerTest","brms","spdep") #pack1<-c("car","glm2", "lme4","nlme","lmerTest","MASS","fitdistrplus") # create vector of packages #pack2<-c("brms", "vcd", "gamlss","rstan") #pack3<-c("caret","e1071","grid","plot3D","plot3Drgl","beanplot","corrplot") #pack4<-c("ggplot2","ggjoy","ggfortify","ggExtra","ggridges") #pack5<-c("psych","betareg","raster","ade4","vegan","pvclust","gclus","cluster","FD") # install packages lapply(pack, install.packages, character.only = TRUE) #lapply(pack1, install.packages, character.only = TRUE) #lapply(pack2, install.packages, character.only = TRUE) #lapply(pack3, install.packages, character.only = TRUE) #lapply(pack4, install.packages, character.only = TRUE) #lapply(pack5, install.packages, character.only = TRUE) ### loading data ### df<-read.table(file = "TV field xp 2018 data.txt", sep = "\t", header = T, dec = ".") str(df) # structure of the dataframe summary(df) # descriptive statistics for each variable # 394 observations over 443 (89 %) where the sown Dalea species matches the harvested one ### filtering / subsetting data ### df2<-df[which(df$filter=='Dalea species match'),] # filters the 394 right observations df2<-df2[which(df2$source!='COMM'),] # excludes the COMM (commercial source, not sown) individuals df2<-df2[which(df2$source!='indigeneous'),] # excludes the indigeneous individuals df2$source<-df2$source[, drop=T] # drops unused levels df2$species<-df2$species[, drop=T] # drops unused levels df2$species.source<-df2$species.source[, drop=T] # drops unused levels df2$testspecies<-df2$testspecies[, drop=T] # drops unused levels df2$filter<-df2$filter[, drop=T] # drops unused levels df2$observedspecies<-df2$observedspecies[, drop=T] # drops unused levels str(df2) summary(df2) ### reordering levels ### print(levels(df2$disturbance_score)) df2$disturbance_score = factor(df2$disturbance_score, levels(df2$disturbance_score)[c(4,2,3,1)]) print(levels(df2$disturbance_score)) print(levels(df2$disturbance_type)) df2$disturbance_type = factor(df2$disturbance_type, levels(df2$disturbance_type)[c(4,1,2,3)]) print(levels(df2$disturbance_type)) print(levels(df2$source)) df2$source = factor(df2$source, levels(df2$source)[c(1,5,2,3,4)]) print(levels(df2$source)) ### creating new local distance variable (distance from lower left quadrant) df2$coldist.quadrant<-round(df2$coldist, 0) df2$rowdist.quadrant<-round(df2$rowdist, 0) df2$radist.quadrant<-sqrt((df2$coldist.quadrant^2)+(df2$rowdist.quadrant^2)) plot(df2$radist.quadrant, df2$radist) # perfect correlation str(as.factor(df2$radist.quadrant)) # 15 levels str(as.factor(df2$radist)) # 135 levels ### 1st visualisation of data ### # histograms par(mar = c(2, 2, 2, 2)) par(mfrow=c(4,4)) hist(df2$SL, prob=T, main="SL", xlab="Shoot length (cm)") # seems normal lines(density(df2[df2$species=="DACA",]$SL), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$SL), col="purple",lwd=2) #lines(seq(0, max(df2$SL)+1, by=.5), dnorm(seq(0, max(df2$SL)+1, by=.5), mean(df2$SL), regate(df2$SL)), col="blue", lty = 2) # fits to a normal distribution #abline(v = summary(df2[df2$species=="DACA",]$SL)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$SL)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) hist(df2$SM, prob=T, main="SM", xlab="Shoot length (cm)") # seems normal lines(density(df2[df2$species=="DACA",]$SM), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$SM), col="purple",lwd=2) #lines(seq(0, max(df2$SM)+1, by=.5), dnorm(seq(0, max(df2$SM)+1, by=.5), mean(df2$SM), regate(df2$SM)), col="blue", lty = 2) # fits to a normal distribution #abline(v = summary(df2[df2$species=="DACA",]$SM)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$SM)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) barplot(prop.table(table(df2$FN)), main="FN", xlab="Number of fresh nodules") #abline(v = summary(df2$FN)[4], col = "red", lwd=2) # adds mean #legend("topright", "Mean", fill="red", cex=0.75) hist(df2$FNRL, prob=T, main="FNRL", xlab="Shoot length (cm)") # seems normal lines(density(df2[df2$species=="DACA",]$FNRL), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$FNRL), col="purple",lwd=2) #lines(seq(0, max(df2$FNRL)+1, by=.5), dnorm(seq(0, max(df2$FNRL)+1, by=.5), mean(df2$FNRL), regate(df2$FNRL)), col="blue", lty = 2) # fits to a normal distribution #abline(v = summary(df2[df2$species=="DACA",]$FNRL)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$FNRL)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) hist(df2$RL, prob=T, main="RL", xlab="Shoot length (cm)") # seems normal lines(density(df2[df2$species=="DACA",]$RL), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$RL), col="purple",lwd=2) #lines(seq(0, max(df2$RL)+1, by=.5), dnorm(seq(0, max(df2$RL)+1, by=.5), mean(df2$RL), regate(df2$RL)), col="blue", lty = 2) # fits to a normal distribution #abline(v = summary(df2[df2$species=="DACA",]$RL)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$RL)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) hist(df2$RM, prob=T, main="RM", xlab="Shoot length (cm)") # seems normal lines(density(df2[df2$species=="DACA",]$RM), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$RM), col="purple",lwd=2) #lines(seq(0, max(df2$RM)+1, by=.5), dnorm(seq(0, max(df2$RM)+1, by=.5), mean(df2$RM), regate(df2$RM)), col="blue", lty = 2) # fits to a normal distribution #abline(v = summary(df2[df2$species=="DACA",]$RM)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$RM)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) barplot(prop.table(table(df2$DN)), main="DN", xlab="Number of decayed nodules") #abline(v = summary(df2$DN)[4], col = "red", lwd=2) # adds mean #legend("topright", "Mean", fill="red", cex=0.75) hist(df2$FNRM, prob=T, main="FNRM", xlab="Shoot length (cm)") # seems noFNRMal lines(density(df2[df2$species=="DACA",]$FNRM), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$FNRM), col="purple",lwd=2) #lines(seq(0, max(df2$FNRM)+1, by=.5), dnoFNRM(seq(0, max(df2$FNRM)+1, by=.5), mean(df2$FNRM), regate(df2$FNRM)), col="blue", lty = 2) # fits to a noFNRMal distribution #abline(v = summary(df2[df2$species=="DACA",]$FNRM)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$FNRM)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) hist(df2$TL, prob=T, main="TL", xlab="Shoot length (cm)") # seems noTLal lines(density(df2[df2$species=="DACA",]$TL), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$TL), col="purple",lwd=2) #lines(seq(0, max(df2$TL)+1, by=.5), dnoTL(seq(0, max(df2$TL)+1, by=.5), mean(df2$TL), regate(df2$TL)), col="blue", lty = 2) # fits to a noTLal distribution #abline(v = summary(df2[df2$species=="DACA",]$TL)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$TL)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) hist(df2$TM, prob=T, main="TM", xlab="Shoot length (cm)") # seems noTMal lines(density(df2[df2$species=="DACA",]$TM), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$TM), col="purple",lwd=2) #lines(seq(0, max(df2$TM)+1, by=.5), dnoTM(seq(0, max(df2$TM)+1, by=.5), mean(df2$TM), regate(df2$TM)), col="blue", lty = 2) # fits to a noTMal distribution #abline(v = summary(df2[df2$species=="DACA",]$TM)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$TM)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) barplot(prop.table(table(df2$TN)), main="TN", xlab="Number of decayed nodules") #abline(v = summary(df2$TN)[4], col = "red", lwd=2) # adds mean #legend("topright", "Mean", fill="red", cex=0.75) hist(df2$DNRL, prob=T, main="DNRL", xlab="Shoot length (cm)") # seems noDNRLal lines(density(df2[df2$species=="DACA",]$DNRL), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$DNRL), col="purple",lwd=2) #lines(seq(0, max(df2$DNRL)+1, by=.5), dnoDNRL(seq(0, max(df2$DNRL)+1, by=.5), mean(df2$DNRL), regate(df2$DNRL)), col="blue", lty = 2) # fits to a noDNRLal distribution #abline(v = summary(df2[df2$species=="DACA",]$DNRL)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$DNRL)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) hist(df2$SRLR, prob=T, main="SRLR", xlab="Shoot length (cm)") # seems noSRLRal lines(density(df2[df2$species=="DACA",]$SRLR), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$SRLR), col="purple",lwd=2) #lines(seq(0, max(df2$SRLR)+1, by=.5), dnoSRLR(seq(0, max(df2$SRLR)+1, by=.5), mean(df2$SRLR), regate(df2$SRLR)), col="blue", lty = 2) # fits to a noSRLRal distribution #abline(v = summary(df2[df2$species=="DACA",]$SRLR)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$SRLR)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) hist(df2$SRMR, prob=T, main="SRMR", xlab="Shoot length (cm)") # seems noSRMRal lines(density(df2[df2$species=="DACA",]$SRMR), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$SRMR), col="purple",lwd=2) #lines(seq(0, max(df2$SRMR)+1, by=.5), dnoSRMR(seq(0, max(df2$SRMR)+1, by=.5), mean(df2$SRMR), regate(df2$SRMR)), col="blue", lty = 2) # fits to a noSRMRal distribution #abline(v = summary(df2[df2$species=="DACA",]$SRMR)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$SRMR)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) barplot(prop.table(table(df2$NS)), main="NS", xlab="Success of nodulation") hist(df2$DNRM, prob=T, main="DNRM", xlab="Shoot length (cm)") # seems noDNRMal lines(density(df2[df2$species=="DACA",]$DNRM), col="black",lwd=2) lines(density(df2[df2$species=="DAPU",]$DNRM), col="purple",lwd=2) #lines(seq(0, max(df2$DNRM)+1, by=.5), dnoDNRM(seq(0, max(df2$DNRM)+1, by=.5), mean(df2$DNRM), regate(df2$DNRM)), col="blue", lty = 2) # fits to a noDNRMal distribution #abline(v = summary(df2[df2$species=="DACA",]$DNRM)[4], col = "black", lwd=2, lty=2) # adds mean #abline(v = summary(df2[df2$species=="DACA",]$DNRM)[4], col = "purple", lwd=2, lty=2) # adds mean #legend("topright", c("Density estimate", "Gaussian estimate", "Mean"), fill=c("black","blue","red"), cex=0.75) #################################################################################### ### 1. Broad question: Plant and rhizobia trait variation - how is it explained? ### #################################################################################### # 1.1 correlation between all traits using the mean by family for each trait (or raw data) # 1.2 use the PCA approach for visual representation and hypothesis generation. # 1.3 use RDA (constrained PCA) ################### # 1.1 correlation # ################### library(corrplot) corrplot(cor(df2[,c(12:14,22:26,28:39)])) # correlation plot all plants corrplot(cor(df2[df2$NS>0,c(12:14,22:26,28:38,40:43)])) # correlation plot nodulated plants corall<-cor(df2[,c(12:14,22:26,28:39)], method="s") # correlation matrix all plants cornodulated<-cor(df2[df2$NS>0,c(12:14,22:26,28:38,40:43)], method="s") # correlation matrix nodulated plants cortraits<-cor(df2[,c(14,22:24,28:39,48,50,56:57)], method="s") # compute correlation test # mat : is a matrix of data # ... : further arguments to pass to the native R cor.test function #cor.mtest <- function(mat, ...) { mat <- as.matrix(mat) n <- ncol(mat) p.mat<- matrix(NA, n, n) diag(p.mat) <- 0 for (i in 1:(n - 1)) { for (j in (i + 1):n) { tmp <- cor.test(mat[, i], mat[, j], ...) p.mat[i, j] <- p.mat[j, i] <- tmp$p.value } } colnames(p.mat) <- rownames(p.mat) <- colnames(mat) p.mat } pmatall <- cor.mtest(df2[,c(12:14,22:26,28:39)], method="s", alternative="two.sided", exact=F, continuity=T) # spearman correlations in matrix pmatnod <- cor.mtest(df2[df2$NS>0,c(12:14,22:26,28:38,40:43)], method="s", alternative="two.sided", exact=F, continuity=T) pmatraits <- cor.mtest(df2[,c(14,22:24,28:39,48,50,56:57)], method="s", alternative="two.sided", exact=F, continuity=T) # spearman correlations in matrix corrplot(corall, method="circle", type="upper", main="Correlation of measured trait on all plants", p.mat = pmatall$p, sig.level = 0.05, mar = c(2,2,4,2)) corrplot(cornodulated, method="circle", type="upper", main="Correlation of measured trait on nodulated plants", p.mat = pmatnod, sig.level = 0.05, mar = c(2,2,4,2)) corrplot(cortraits, number.cex = .5, main="Correlogram", p.mat = pmatraits$p, insig = "pch", sig.level = .05, pch.cex = 1, pch.col = "black", mar = c(1,1,2,1), order = "hclust", hclust.method="ward.D2", addrect = 2) cortraits[c(5,6,9,10,13,15),c(5,6,9,10,13,15)] pmatraits$p[c(5,6,9,10,13,15),c(5,6,9,10,13,15)] detach("package:corrplot", unload=TRUE) corprimary<-cor(df2[,c(28:29,32:33,36:37)], method="s") # correlation matrix all plants pmatprimary<-cor.mtest(df2[,c(28:29,32:33,36:37)], method="s", alternative="two.sided", exact=F, continuity=T, conf.level = .95) # spearman correlations in matrix corrplot(corprimary, method="circle", type="upper", main="Correlogram", p.mat = pmatprimary$p, sig.level = 0.05, mar = c(2,2,4,2)) corprimary pmatprimary # correlation results # on all plants # radial distance positively correlated with column and row disances: use radial distance only # latitude negatively correlated with longitude and distance as the crow flies: use distance as the crow flies only # SL and RL not correlated: use both # SM and RM positively correlated: could use TM instead # FN and DN not correlated: use both # misc.: longitude correlated with SRMR, FN and TN # on nodulated plants # almost same as above # FS, FP, DS and DP correlated together and with FN and DN: do no add any information # misc.: FN correlated with column distance (E-W gradient in field area) # misc.: longitude correlated withSRMR and TN ################### # 1.2 ordinations # ################### library(ggfortify) library(vegan) # 1.2.1 # unconstrained ordination (PCA) # principle: ordinate traits (SL, RL, SM, RM, FN and DN) without a priori hypotheses # then test the ordinations scores on each axis with species*source, plus distance to source and local distance # then test/map the other traits and environment variables # PCA on both species, using SL, RL, SM, RM, FN and DN variables pca.all<-prcomp(df2[,c(28:29,32:33,36:37)], center=T, scale=T, rank.=2) pca.all<-prcomp(df2[,c(28:29,32:33,36:37)], center=T, scale=T) summary(pca.all) # PCA axis 1 explains 42.68% of variance, axis 2 17.42% (total 60.10%) pca.all<-(autoplot(pca.all, data=df2, colour="species", shape="source", size=2, loadings=T, loadings.colour='Black', loadings.label.colour='black', loadings.label.vjust = -1.25, loadings.label=T, loadings.label.size=3, frame=F, frame.type="convex")+ theme_bw() + theme(legend.justification=c(1,1), legend.position=c(.95,.95))) # plant traits strongly (RM, SM) or rather (RL, SL) positively correlate together and to axis 1 # FN and DN correlate positively together and with axis 2 # root and shoot traits negatively correlate (they are separated by axis 2): nodules and investment in roots are correlated summary(aov(scores(pca.all)[,1]~species*source+radist, data=df2)) # effect on axis 1 of PCA (PC1): # species ** (1 df; F=7.299; p=0.00720) 30% of variation (Mean Sq) # source ns (4 df; F=1.184; p=0.317470) 5% of variation # species:source * (4 Df; F=2.874; p=0.022824) 12% of variation # radial distance *** (1 df; F=12.350; p=0.000494) 50% of variation # residual (383 df) 3% summary(aov(scores(pca.all)[,1]~species*crowdist+radist, data=df2)) # species:source no longer significant summary(aov(scores(pca.all)[,2]~species*source+radist, data=df2)) # effect on axis 2 of PCA (PC2): # species ns (1 df; F=0.222; p=0.638045) 2% of variation (Mean Sq) # source *** (4 df; F=5.641; p=0.000202) 61% of variation # radial distance ns (1 df; F=0.005; p=0.944191) 1% of variation # species:source . (4 Df; F=2.272; p=0.060929) 25% of variation # residual (383 df) 11% summary(aov(scores(pca.all)[,2]~species*crowdist+radist, data=df2)) # no changes summary(aov(scores(pca.all)[,1]~species*source+radist.quadrant, data=df2)) summary(aov(scores(pca.all)[,1]~species*crowdist+radist.quadrant, data=df2)) summary(aov(scores(pca.all)[,2]~species*source+radist.quadrant, data=df2)) summary(aov(scores(pca.all)[,2]~species*crowdist+radist.quadrant, data=df2)) # accurate results of anova on PCA scores pca.1.id<-anova(lm(scores(pca.all)[,1]~species*source+radist, data=df2)) afss.pca.1.id<-pca.1.id$"Sum Sq" round(print(cbind(pca.1.id,PctExp=afss.pca.1.id/sum(afss.pca.1.id)*100)), digits = 4) pca.1.dts<-anova(lm(scores(pca.all)[,1]~species*crowdist+radist, data=df2)) afss.pca.1.dts<-pca.1.dts$"Sum Sq" round(print(cbind(pca.1.dts,PctExp=afss.pca.1.dts/sum(afss.pca.1.dts)*100)), digits = 4) pca.2.id<-anova(lm(scores(pca.all)[,2]~species*source+radist, data=df2)) afss.pca.2.id<-pca.2.id$"Sum Sq" round(print(cbind(pca.2.id,PctExp=afss.pca.2.id/sum(afss.pca.2.id)*100)), digits = 4) pca.2.dts<-anova(lm(scores(pca.all)[,2]~species*crowdist+radist, data=df2)) afss.pca.2.dts<-pca.2.dts$"Sum Sq" round(print(cbind(pca.2.dts,PctExp=afss.pca.2.dts/sum(afss.pca.2.dts)*100)), digits = 4) # fitting environmental variables to unconstrained ordination fit.all<-envfit(ord = pca.all, env = df2[,c(12:14,18,21,22:24,30:31,34:35,38:39,48,50,56,57,58)], permutations = 1000) fit.all vecfit.all<-envfit(ord = pca.all, env = df2[,c(12:14,22:24,30:31,34:35,38:39,48,50,56,57,58)], permutations = 1000) # radist **, while coldist and rowdist both ns: there is a spatial structure but not by columns or rows of cells or quadrants # crowdist **, with latitude ns and longitude *** : the source to introduction distance is significant because of a longitutinal effect # all additional plant and nodule traits significant: expected since they are correlated to primary traits, they should map together biplot(pca.all, cex=1, xlabs=rep("°",394)) plot(vecfit.all, p.max=0.05, cex=0.5) # local spatial structure along PC1 (reminder: species** and species:source *) # source to introduction site distance and longitude correlate positively together along PC2 (reminder: source ***) and negatively correlate with number of nodules # Shoot:Root ratios (length and mass) correlate together along PC2 and negatively correlate with number of nodules # to better see the difference between sources fit.source<-envfit(ord = pca.all, env = df2[,c(21:24)], permutations = 1000) fit.source biplot(cex=0.25, pca.all) plot(fit.source, p.max=0.05, cex=0.5) # ranks of centroids on PC1 (longer/heavier first): CHIP, KWD, PEMB, BLUE, GLYNN # ranks of centroids on PC2 (higher number of nodules firts): BLUE, GLYNN, CHIP, PEMB, KWD cor.test(fit.source$factors$centroids[,1],fit.source$factors$centroids[,2], alternative="t", method="s") # no significant rank correlation for source between plant lenght/mass and nodule number # to better see the difference between species fit.species<-envfit(ord = pca.all, env = df2[,c(1,18)], permutations = 1000) fit.species biplot(cex=0.5, pca.all) plot(fit.species, p.max=0.05, cex=0.5) # ranks of centroids on PC1 (longer/heavier first): DACA, DAPU # ranks of centroids on PC2 (higher number of nodules firts): DAPU, DACA cor.test(fit.species$factors$centroids[,1],fit.species$factors$centroids[,2], alternative="t", method="s") # no significant rank correlation for species between plant lenght/mass and nodule number # to better see the difference between species:source interaction fit.species.source<-envfit(ord = pca.all, env = df2[,c(18,21,27)], permutations = 1000) fit.species.source biplot(cex=0.5, pca.all) plot(fit.species.source, p.max=0.05, cex=0.5) cor.test(fit.species.source$factors$centroids[,1],fit.species.source$factors$centroids[,2], alternative="t", method="s") # no significant rank correlation for species:source between plant lenght/mass and nodule number # correlation between PCA scores of population (species*source) and distance to source cor.test(fit.species.source$factors$centroids[-(1:7),1], crowdist2) # no correlation between scores of populations on PC1 and distance to source cor.test(fit.species.source$factors$centroids[-(1:7),2], crowdist2) # significant negative correlation between scores of populations on PC2 and distance to source # the closer from the source, the higher the scores on PC2 # 1.2.2 # PCA on each species, using SL, RL, SM, RM, FN and DN variables pca.dapu<-prcomp(df2[df2$species=='DAPU',c(28:29,32:33,36:37)], center=T, scale=T, rank.=2) summary(pca.dapu) # need to invert PC1 axis on pca pca.dapu$rotation[,1]<- -pca.dapu$rotation[,1] # PCA axis 1 explains 44.6% of variance, axis 2 20.28% (total 64.8%) autoplot(pca.dapu, data=df2[df2$species=='DAPU',], colour="species", shape="source", size=1, loadings=T, loadings.colour='red', loadings.label.colour='red', loadings.label=T, loadings.label.size=3, frame=F, frame.type="convex")+ ggtitle("PCA on DAPU only")+ theme_bw() # exactly same as when both species analyzed together (see above) summary(aov(scores(pca.dapu)[,1]~source+radist, data=df2[df2$species=='DAPU',])) # effect on axis 1 of PCA (PC1): # source ns # radial distance ** summary(aov(scores(pca.dapu)[,1]~crowdist+radist, data=df2[df2$species=='DAPU',])) # same as above summary(aov(scores(pca.dapu)[,2]~source+radist, data=df2[df2$species=='DAPU',])) # effect on axis 2 of PCA (PC2): # source ns # radial distance ns summary(aov(scores(pca.dapu)[,2]~crowdist+radist, data=df2[df2$species=='DAPU',])) # same as above # accurate PCA results for DAPU pca.1.dapu.id<-anova(lm(scores(pca.dapu)[,1]~source+radist, data=df2[df2$species=='DAPU',])) afss.pca.1.dapu.id<-pca.1.dapu.id$"Sum Sq" round(print(cbind(pca.1.dapu.id,PctExp=afss.pca.1.dapu.id/sum(afss.pca.1.dapu.id)*100)), digits = 4) pca.2.dapu.id<-anova(lm(scores(pca.dapu)[,2]~source+radist, data=df2[df2$species=='DAPU',])) afss.pca.2.dapu.id<-pca.2.dapu.id$"Sum Sq" round(print(cbind(pca.2.dapu.id,PctExp=afss.pca.2.dapu.id/sum(afss.pca.2.dapu.id)*100)), digits = 4) pca.1.dapu.dts<-anova(lm(scores(pca.dapu)[,1]~crowdist+radist, data=df2[df2$species=='DAPU',])) afss.pca.1.dapu.dts<-pca.1.dapu.dts$"Sum Sq" round(print(cbind(pca.1.dapu.dts,PctExp=afss.pca.1.dapu.dts/sum(afss.pca.1.dapu.dts)*100)), digits = 4) pca.2.dapu.dts<-anova(lm(scores(pca.dapu)[,2]~crowdist+radist, data=df2[df2$species=='DAPU',])) afss.pca.2.dapu.dts<-pca.2.dapu.dts$"Sum Sq" round(print(cbind(pca.2.dapu.dts,PctExp=afss.pca.2.dapu.dts/sum(afss.pca.2.dapu.dts)*100)), digits = 4) # accurate PCA results for DACA pca.1.daca.id<-anova(lm(scores(pca.daca)[,1]~source+radist, data=df2[df2$species=='DACA',])) afss.pca.1.daca.id<-pca.1.daca.id$"Sum Sq" round(print(cbind(pca.1.daca.id,PctExp=afss.pca.1.daca.id/sum(afss.pca.1.daca.id)*100)), digits = 4) pca.2.daca.id<-anova(lm(scores(pca.daca)[,2]~source+radist, data=df2[df2$species=='DACA',])) afss.pca.2.daca.id<-pca.2.daca.id$"Sum Sq" round(print(cbind(pca.2.daca.id,PctExp=afss.pca.2.daca.id/sum(afss.pca.2.daca.id)*100)), digits = 4) pca.1.daca.dts<-anova(lm(scores(pca.daca)[,1]~crowdist+radist, data=df2[df2$species=='DACA',])) afss.pca.1.daca.dts<-pca.1.daca.dts$"Sum Sq" round(print(cbind(pca.1.daca.dts,PctExp=afss.pca.1.daca.dts/sum(afss.pca.1.daca.dts)*100)), digits = 4) pca.2.daca.dts<-anova(lm(scores(pca.daca)[,2]~crowdist+radist, data=df2[df2$species=='DACA',])) afss.pca.2.daca.dts<-pca.2.daca.dts$"Sum Sq" round(print(cbind(pca.2.daca.dts,PctExp=afss.pca.2.daca.dts/sum(afss.pca.2.daca.dts)*100)), digits = 4) # fitting environmental variables to unconstrained ordination fit.dapu<-envfit(ord = pca.dapu, env = df2[df2$species=='DAPU',c(12:14,18,21,22:24,30:31,34:35,38:39,48,50,56,57)], permutations = 1000) fit.dapu vecfit.dapu<-envfit(ord = pca.dapu, env = df2[df2$species=='DAPU',c(12:14,22:24,30:31,34:35,38:39,48,50,56,57)], permutations = 1000) # need to invert PC1 vecfit.dapu$vectors$arrows[,1]<- -vecfit.dapu$vectors$arrows[,1] # radist * and rowdist*, while coldist ns: there is a spatial structure by row of cells/quadrants (N-S gradient) # no effect of source to introduction distance # all additional plant and nodule traits significant: expected since they are correlated to primary traits, they should map together biplot(cex=0.5, pca.dapu, xlabs=rep("º",213), main="PCA and variable fitting on DAPU") plot(vecfit.dapu, p.max=0.05, cex=0.5) # factors omitted to improve readability # local spatial structure (radist and rowdist) along PC1 (reminder: radial distance **) # Shoot:Root ratios (length and mass) correlate positively together along PC2 and negatively correlate with number of nodules vecfit.dapu2<-envfit(ord = pca.dapu, env = df2[df2$species=='DAPU',c(12:14,22:24,30:31,34:35,38:39,48:51)], permutations = 1000) vecfit.dapu2$vectors$arrows[,1]<- -vecfit.dapu2$vectors$arrows[,1] biplot(cex=1, pca.dapu, xlabs=rep("°",213), main="PCA and variable fitting DAPU") plot(vecfit.dapu2, p.max=0.05, cex=.75) # with number of nodules standardized by plant trait pca.daca<-prcomp(df2[df2$species=='DACA',c(28:29,32:33,36:37)], center=T, scale=T, rank.=2) summary(pca.daca) # need to invert both PC1 and PC2 pca.daca$rotation[,1]<- -pca.daca$rotation[,1] pca.daca$rotation[,2]<- -pca.daca$rotation[,2] # PCA axis 1 explains 41.9% of variance, axis 2 18.26% (total 60.10%) autoplot(pca.daca, data=df2[df2$species=='DACA',], colour="species", shape="source", size=1, loadings=T, loadings.colour='red', loadings.label.colour='red', loadings.label=T, loadings.label.size=3, frame=F, frame.type="convex")+ ggtitle("PCA of measured trait on daca")+ theme_bw() # differs from when both species analyzed together or from dapu analyzed alone # fresh nodules and decayed nodules negatively correlate along axis 2 summary(aov(scores(pca.daca)[,1]~source+radist, data=df2[df2$species=='DACA',])) # effect on axis 1 of PCA (PC1): # source * # radial distance * summary(aov(scores(pca.daca)[,1]~crowdist+radist, data=df2[df2$species=='DACA',])) # all ns summary(aov(scores(pca.daca)[,2]~source+radist, data=df2[df2$species=='DACA',])) # effect on axis 2 of PCA (PC2): # source ns # radial distance ns summary(aov(scores(pca.daca)[,2]~crowdist+radist, data=df2[df2$species=='DACA',])) # all ns # fitting environmental variables to unconstrained ordination fit.daca<-envfit(ord = pca.daca, env = df2[df2$species=='DACA',c(12:14,18,21,22:24,30:31,34:35,38:39,48,50,56,57)], permutations = 1000) fit.daca vecfit.daca<-envfit(ord = pca.daca, env = df2[df2$species=='DACA',c(12:14,22:24,30:31,34:35,38:39,48,50,56,57)], permutations = 1000) # need to invert PC1 and PC2 vecfit.daca$vectors$arrows[,1]<- -vecfit.daca$vectors$arrows[,1] vecfit.daca$vectors$arrows[,2]<- -vecfit.daca$vectors$arrows[,2] # no significant effect of distance to source or local distance # all additional plant and nodule traits significant: expected since they are correlated to primary traits, they should map together biplot(cex=0.5, pca.daca, xlabs=rep("º",181), main="PCA and variable fitting on DACA") plot(vecfit.daca, p.max=0.05, cex=0.5) # factors omitted to improve readability # local spatial structure (radist and rowdist) along PC1 (reminder: radial distance **) # Shoot:Root ratios (length and mass) correlate together along PC2 and negatively correlate with number of nodules vecfit.daca2<-envfit(ord = pca.daca, env = df2[df2$species=='DACA',c(12:14,22:24,30:31,34:35,38:39,48:51)], permutations = 1000) vecfit.daca2$vectors$arrows[,1]<- -vecfit.daca2$vectors$arrows[,1] vecfit.daca2$vectors$arrows[,2]<- -vecfit.daca2$vectors$arrows[,2] biplot(cex=0.5, pca.daca, xlabs=rep("º",181), main="PCA and variable fitting #2 on DACA") plot(vecfit.daca2, p.max=0.3, cex=0.5) # with number of nodules standardized by plant trait names(df2[,c(12:14,22:24,30:31,34:35,38:39,48:51)]) names(df2) # biplots with ggplot-like package, for Fig 2C-D for Georgiana biplot(cex=1, pca.dapu, xlabs=rep("°",213), main="PCA and variable fitting DAPU") plot(vecfit.dapu2, p.max=0.05, cex=.75) # with number of nodules standardized by plant trait vecfit.dapu3<-envfit(ord = pca.dapu, env = df2[df2$species=='DAPU',c(8:14,21:24,30:31,34:35,38:39,48:51)], permutations = 1000) vecfit.dapu3 biplot(cex=1, pca.daca, xlabs=rep("°",181), main="PCA and variable fitting DACA") plot(vecfit.daca2, p.max=0.05, cex=.75) # with number of nodules standardized by plant trait vecfit.daca3<-envfit(ord = pca.daca, env = df2[df2$species=='DACA',c(8:14,21:24,30:31,34:35,38:39,48:51)], permutations = 1000) vecfit.daca3 str(df2) # 1.3 # constrained ordination (CCA) # principle: ordinate traits (SL, RL, SM, RM, FN and DN) with a priori hypotheses # based on CCA, use species*source, radist and crowdist to constrain # cca.all<-cca(df2[,c(28:29,32:33,36:37)]~species*source+crowdist+radist, data=df2) # anova(cca.all, by="axis", perm=1000) # 2 first axes are significant # anova(cca.all, by="term", perm=1000) # confirms species, source, crowdist and radist are significant # cca.all2<-cca(df2[,c(28:39)]~species*source+crowdist+radist, data=df2) # anova(cca.all2, by="axis", perm=1000) # first axis is significant, second almost # anova(cca.all2, by="term", perm=1000) # confirms species, source, crowdist and radist are significant # cca.dapu<-cca(df2[df2$species=='DAPU',c(28:39)]~species*source+crowdist+radist, data=df2[df2$species=='DAPU',]) # anova(cca.dapu, by="axis", perm=1000) # first axis is significant, second ns # anova(cca.dapu, by="term", perm=1000) # confirms radist is significant # cca.daca<-cca(df2[df2$species=='DACA',c(28:39)]~species*source+crowdist+radist, data=df2[df2$species=='DACA',]) # anova(cca.daca, by="axis", perm=1000) # first axis is significant, second almost # anova(cca.daca, by="term", perm=1000) # confirms source and radist are significant detach("package:ggfortify", unload=TRUE) detach("package:vegan", unload=TRUE) ###################################################### # 2. Evidence for local adaptation? -models by trait # ###################################################### library(lme4) library(lmerTest) library(brms) ############## # 2.1 Models # ############## ########################################### # 2.1.1 lmer/glmer, with source as factor # ########################################### # having both source and distance from source as fixed effects makes the matrix rank deficient # how to test for source effect: factor (source) or variable (distance from source)? # I would opt for using the source in the model, and the distance from source in regressions for each species # m.SL, shoot length model m.SL.rand <- lmer(SL~species*source+radist+(1|quadrant), data=df2, REML=F) m.SL <- lm(SL~species*source+radist, data=df2) m.SL.rq <- lm(SL~species*source+radist.quadrant, data=df2) plot(m.SL, main="Residual vs. fitted") summary(m.SL) anova(m.SL) # signifiance (NumDf, DenDf, F, p-value) # species * # source ns # radist * # species:source ** ranova(m.SL) # signifiance (Df, p-value) # quadrants *** (1 0.0001994) m.SLbayes<-brm(SL~species*source+radist+(1|quadrant), gaussian(link = "identity"), data=df2, seed=1234) summary(m.SLbayes) # bayesian estimations of effects round(ranef(m.SLbayes)$quadrant,2) # effect of each quadrants on shoot length # m.RL, root length model m.RL.rand <- lmer(RL~species*source+radist+(1|quadrant), data=df2, REML=F) m.RL <- lm(RL~species*source+radist, data=df2) m.RL.rq <- lm(RL~species*source+radist.quadrant, data=df2) plot(m.RL, main="Residual vs. fitted") summary(m.RL) anova(m.RL) # signifiance (NumDf, DenDf, F, p-value) # species ns # source ns # radist ns # species:source ns ranova(m.RL) # signifiance (Df, p-value) # quadrants *** (1 3.09e-06) m.RLbayes<-brm(RL~species*source+radist+(1|quadrant), gaussian(link = "identity"), data=df2, seed=1234) summary(m.RLbayes) # bayesian estimations of effects round(ranef(m.RLbayes)$quadrant,2) # effect of each quadrants on root length # m.TL, total length model m.TL.rand <- lmer(TL~species*source+radist+(1|quadrant), data=df2, REML=F) m.TL <- lm(TL~species*source+radist, data=df2) m.TL.rq <- lm(TL~species*source+radist.quadrant, data=df2) plot(m.TL, main="Residual vs. fitted") summary(m.TL) anova(m.TL.rq) # signifiance (NumDf, DenDf, F, p-value) # species . # source ns # radist . # species:source ns ranova(m.TL) # signifiance (Df, p-value) # quadrants *** (1 0.0008408) m.TLbayes<-brm(TL~species*source+radist+(1|quadrant), gaussian(link = "identity"), data=df2, seed=1234) summary(m.TLbayes) # bayesian estimations of effects round(ranef(m.TLbayes)$quadrant,2) # effect of each quadrants on total length # m.SRLR S:R length ratio m.SRLR.rand <- glmer(SRLR~species*source+radist+(1|quadrant), data=df2, family = Gamma(link = 'log')) m.SRLR <- glm(SRLR~species*source+radist, data=df2, family = Gamma(link = 'log')) m.SRLR.rq <- glm(SRLR~species*source+radist.quadrant, data=df2, family = Gamma(link = 'log')) # model fails to converge, max|grad|>0.001, check relgrad.m.SRLR <- with(m.SRLR@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.SRLR)) # <0.001 threshold, should be fine plot(m.SRLR, main="Residual vs. fitted") summary(m.SRLR) anova(m.SRLR.rq, test='F') # all ns m.SRLRbayes<-brm(SRLR~species*source+(1|quadrant), family=c("gamma","log"), data=df2, seed=1234) summary(m.SRLRbayes) # bayesian estimations of effects round(ranef(m.SRLRbayes)$quadrant,2) # effect of each quadrants on s:r length ratio # m.SM shoot dry mass m.SM.rand <- glmer(SM~species*source+radist+(1|quadrant), data=df2, family = Gamma(link = 'log')) # model fails to converge, max|grad|>0.001, check relgrad.m.SM <- with(m.SM@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.SM)) # <0.001 threshold, should be fine m.SM <- glm(SM~species*source+radist, data=df2, family = Gamma(link = 'log')) m.SM.rq <- glm(SM~species*source+radist.quadrant, data=df2, family = Gamma(link = 'log')) plot(m.SM, main="Residual vs. fitted") summary(m.SM) anova(m.SM, test='F') # species *** # source . # radist *** # species:source ** m.SMbayes<-brm(SM~species*source+(1|quadrant), family=c("gamma","log"), data=df2, seed=1234) summary(m.SMbayes) # bayesian estimations of effects round(ranef(m.SMbayes)$quadrant,2) # effect of each quadrants on shoot dry mass # m.RM root dry mass m.RM.rand <- glmer(RM~species*source+radist+(1|quadrant), data=df2, family = Gamma(link = 'log')) # model fails to converge, max|grad|>0.001, check relgrad.m.RM <- with(m.RM@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.RM)) # <0.001 threshold, should be fine m.RM <- glm(RM~species*source+radist, data=df2, family = Gamma(link = 'log')) m.RM.rq <- glm(RM~species*source+radist.quadrant, data=df2, family = Gamma(link = 'log')) plot(m.RM, main="Residual vs. fitted") summary(m.RM) anova(m.RM.rq, test='F') m.RMbayes<-brm(RM~species*source+(1|quadrant), family=c("gamma","log"), data=df2, seed=1234) summary(m.RMbayes) # bayesian estimations of effects round(ranef(m.RMbayes)$quadrant,2) # effect of each quadrants on root dry mass # m.TM total dry mass m.TM.rand <- glmer(TM~species*source+radist+(1|quadrant), data=df2, family = Gamma(link = 'log')) # model fails to converge, max|grad|>0.001, check relgrad.m.TM <- with(m.TM@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.TM)) # <0.001 threshold, should be fine m.TM <- glm(TM~species*source+radist, data=df2, family = Gamma(link = 'log')) m.TM.rq <- glm(TM~species*source+radist.quadrant, data=df2, family = Gamma(link = 'log')) plot(m.TM, main="Residual vs. fitted") summary(m.TM) anova(m.TM, test='F') m.TMbayes<-brm(TM~species*source+(1|quadrant), family=c("gamma","log"), data=df2, seed=1234) summary(m.TMbayes) # bayesian estimations of effects round(ranef(m.TMbayes)$quadrant,2) # effect of each quadrants on total dry mass # m.SRMR S:R dry mass ratio m.SRMR.rand <- glmer(SRMR~species*source+radist+(1|quadrant), data=df2, family = Gamma(link = 'log')) # model fails to converge, max|grad|>0.001, check relgrad.m.SRMR <- with(m.SRMR@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.SRMR)) # <0.001 threshold, should be fine m.SRMR <- glm(SRMR~species*source+radist, data=df2, family = Gamma(link = 'log')) m.SRMR.rq <- glm(SRMR~species*source+radist.quadrant, data=df2, family = Gamma(link = 'log')) plot(m.SRMR, main="Residual vs. fitted") summary(m.SRMR) anova(m.SRMR.rq, test='F') m.SRMRbayes<-brm(SRMR~species*source+(1|quadrant), family=c("gamma","log"), data=df2, seed=1234) summary(m.SRMRbayes) # bayesian estimations of effects round(ranef(m.SRMRbayes)$quadrant,2) # effect of each quadrants on s:r dry mass ratio # m.FN number of fresh nodules m.FN.rand<-glmer(df2$FN~species*source+radist+(1|quadrant), family="poisson", data=df2) # model fails to converge, max|grad|>0.001, check relgrad.m.FN <- with(m.FN@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.FN)) # <0.001 threshold, should be fine m.FN<-glm(df2$FN~species*source+radist, family="poisson", data=df2) m.FN.rq<-glm(df2$FN~species*source+radist.quadrant, family="poisson", data=df2) plot(m.FN) summary(m.FN) anova(m.FN, test='Chisq') m.FNbayes<-brm(FN~species*source+radist+(1|quadrant), family=c("poisson","log"), data=df2, seed=1234) summary(m.FNbayes) # bayesian estimations of effects round(ranef(m.FNbayes)$quadrant,2) # effect of each quadrants of nb fresh nod # m.DN number of decayed nodules m.DN.rand<-glmer(df2$DN~species*source+radist+(1|quadrant), family="poisson", data=df2) # model fails to converge, max|grad|>0.001, check relgrad.m.DN <- with(m.DN@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.DN)) # <0.001 threshold, should be fine m.DN<-glm(df2$DN~species*source+radist, family="poisson", data=df2) plot(m.DN) summary(m.DN) anova(m.DN, test='Chisq') m.DNbayes<-brm(DN~species*source+radist+(1|quadrant), family=c("poisson","log"), data=df2, seed=1234) summary(m.DNbayes) # bayesian estimations of effects round(ranef(m.DNbayes)$quadrant,2) # effect of each quadrants of nb decayed nod # m.TN number of total nodules m.TN.rand<-glmer(df2$TN~species*source+radist+(1|quadrant), family="poisson", data=df2) # model fails to converge, max|grad|>0.001, check relgrad.m.TN <- with(m.TN@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.TN)) # <0.001 threshold, should be fine m.TN<-glm(df2$TN~species*source+radist, family="poisson", data=df2) plot(m.TN) summary(m.TN) anova(m.TN, test='Chisq') m.TNbayes<-brm(TN~species*source+radist+(1|quadrant), family=c("poisson","log"), data=df2, seed=1234) summary(m.TNbayes) # bayesian estimations of effects round(ranef(m.TNbayes)$quadrant,2) # effect of each quadrants of nb total nod # m.NS nodulation success NS m.NS.rand<-glmer(NS~species*source+radist+(1|quadrant), family=binomial(link="logit"), data=df2) # model fails to converge, max|grad|>0.001, check relgrad.m.NS <- with(m.NS@optinfo$derivs,solve(Hessian,gradient)) max(abs(relgrad.m.NS)) # <0.001 threshold, should be fine m.NS<-glm(NS~species*source+radist, family=binomial(link="logit"), data=df2) summary(m.NS) anova(m.NS, test='Chisq') m.NSbayes<-brm(NS~species*source+(1|quadrant), family=c("binomial","logit"), data=df2, seed=1234) summary(m.NSbayes) # bayesian estimations of effects round(ranef(m.NSbayes)$quadrant,2) # effect of each quadrants of nod success # m.FS fresh nodule succes model # only 284 individuals on the 394 total m.FS.rand<-glmer(FS~species*source+radist+(1|quadrant), family=binomial(link="logit"), data=df2[df2$NS > 0,]) # singularity, trying to drop the quadrants since there's less data m.FS<-glm(FS~species*source+radist, family=binomial(link="logit"), data=df2[df2$NS > 0,]) plot(m.FS) summary(m.FS) anova(m.FS, test='Chisq') m.FSbayes<-brm(FS~species*source+radist+(1|quadrant), family=c("binomial","logit"), data=df2[df2$NS > 0,], seed=1234) summary(m.FSbayes) # bayesian estimations of effects round(ranef(m.FSbayes)$quadrant,2) # effect of each quadrants # m.FP fresh nodule proportion model # only 284 individuals on the 394 total # dropping the quadrants since there's less data m.FP<-glm(FP~species*source+radist, family=quasibinomial(link="logit"), data=df2[df2$NS > 0,]) plot(m.FP) summary(m.FP) anova(m.FP, test='Chisq') m.FPbayes<-brm(FN|trials(TN)~species*source+radist+(1|quadrant), family=c("binomial","logit"), data=df2[df2$NS > 0,], seed=1234) summary(m.FPbayes) # bayesian estimations of effects round(ranef(m.FPbayes)$quadrant,2) # effect of each quadrants # m.DS decayed nodule succes model # only 284 individuals on the 394 total m.DS<-glmer(DS~species*source+radist+(1|quadrant), family=binomial(link="logit"), data=df2[df2$NS > 0,]) # fails to converge, trying to drop the quadrants since there's less data m.DS<-glm(DS~species*source+radist, family=binomial(link="logit"), data=df2[df2$NS > 0,]) plot(m.DS) summary(m.DS) anova(m.DS, test='Chisq') m.DSbayes<-brm(DS~species*source+radist+(1|quadrant), family=c("binomial","logit"), data=df2[df2$NS > 0,], seed=1234) summary(m.DSbayes) # bayesian estimations of effects round(ranef(m.DSbayes)$quadrant,2) # effect of each quadrants # m.DP decayed nodule proportion model # only 284 individuals on the 394 total # dropping the quadrants since there's less data m.DP<-glm(DP~species*source+radist, family=quasibinomial(link="logit"), data=df2[df2$NS > 0,]) plot(m.DP) summary(m.DP) anova(m.DP, test='Chisq') m.DPbayes<-brm(DN|trials(TN)~species*source+radist+(1|quadrant), family=c("binomial","logit"), data=df2[df2$NS > 0,], seed=1234) summary(m.DPbayes) # bayesian estimations of effects round(ranef(m.DPbayes)$quadrant,2) # effect of each quadrants # FNRL df2$FNRL<-df2$FN/df2$RL m.FNRL <- glm(FNRL+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.FNRL, main="Residual vs. fitted") summary(m.FNRL) anova(m.FNRL, test='F') # TNRL df2$TNRL<-df2$TN/df2$RL m.TNRL <- glm(TNRL+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.TNRL, main="Residual vs. fitted") summary(m.TNRL) anova(m.TNRL, test='F') # FNRM df2$FNRM<-df2$FN/df2$RM m.FNRM <- glm(FNRM+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.FNRM, main="Residual vs. fitted") summary(m.FNRM) anova(m.FNRM, test='F') # TNRM df2$TNRM<-df2$TN/df2$RM m.TNRM <- glm(TNRM+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.TNRM, main="Residual vs. fitted") summary(m.TNRM) anova(m.TNRM, test='F') # FNTL df2$FNTL<-df2$FN/df2$RL m.FNTL <- glm(FNTL+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.FNTL, main="Residual vs. fitted") summary(m.FNTL) anova(m.FNTL, test='F') # TNTL df2$TNTL<-df2$TN/df2$RL m.TNTL <- glm(TNTL+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.TNTL, main="Residual vs. fitted") summary(m.TNTL) anova(m.TNTL, test='F') # FNTM df2$FNTM<-df2$FN/df2$RM m.FNTM <- glm(FNTM+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.FNTM, main="Residual vs. fitted") summary(m.FNTM) anova(m.FNTM, test='F') # TNTM df2$TNTM<-df2$TN/df2$RM m.TNTM <- glm(TNTM+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.TNTM, main="Residual vs. fitted") summary(m.TNTM) anova(m.TNTM, test='F') # DNRL df2$DNRL<-df2$DN/df2$RL m.DNRL <- glm(DNRL+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.DNRL, main="Residual vs. fitted") summary(m.DNRL) anova(m.DNRL, test='F') # DNRM df2$DNRM<-df2$DN/df2$RM m.DNRM <- glm(DNRM+0.001~species*source+radist, data=df2, family = Gamma(link = 'log')) plot(m.DNRM, main="Residual vs. fitted") summary(m.DNRM) anova(m.DNRM, test='F') detach("package:lme4", unload=TRUE) detach("package:lmerTest", unload=TRUE) detach("package:brms", unload=TRUE) ############################################# # 2.1.2 lmer/glmer, with distance to source # ############################################# m.SL.dts <- lm(SL~species*crowdist+radist, data=df2) m.RL.dts <- lm(RL~species*crowdist+radist, data=df2) m.TL.dts <- lm(TL~species*crowdist+radist, data=df2) m.SRLR.dts <- glm(SRLR~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) m.SM.dts <- glm(SM~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) m.RM.dts <- glm(RM~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) m.TM.dts <- glm(TM~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) m.SRMR.dts <- glm(SRMR~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) m.FN.dts <- glm(FN~species*crowdist+radist, family="poisson", data=df2) m.DN.dts <- glm(DN~species*crowdist+radist, family="poisson", data=df2) m.TN.dts <- glm(TN~species*crowdist+radist, family="poisson", data=df2) m.NS.dts <- glm(NS~species*crowdist+radist, family=binomial(link="logit"), data=df2) m.FS.dts <- glm(FS~species*crowdist+radist, family=binomial(link="logit"), data=df2[df2$NS > 0,]) m.FP.dts <- glm(FP~species*crowdist+radist, family=quasibinomial(link="logit"), data=df2[df2$NS > 0,]) m.DS.dts <- glm(DS~species*crowdist+radist, family=binomial(link="logit"), data=df2[df2$NS > 0,]) m.DP.dts <- glm(DP~species*crowdist+radist, family=quasibinomial(link="logit"), data=df2[df2$NS > 0,]) # FNRL df2$FNRL<-df2$FN/df2$RL m.FNRL.dts <- glm(FNRL+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.FNRL.dts, main="Residual vs. fitted") summary(m.FNRL.dts) anova(m.FNRL.dts, test='F') # TNRL df2$TNRL<-df2$TN/df2$RL m.TNRL.dts <- glm(TNRL+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.TNRL.dts, main="Residual vs. fitted") summary(m.TNRL.dts) anova(m.TNRL.dts, test='F') # FNRM df2$FNRM<-df2$FN/df2$RM m.FNRM.dts <- glm(FNRM+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.FNRM.dts, main="Residual vs. fitted") summary(m.FNRM.dts) anova(m.FNRM.dts, test='F') # TNRM df2$TNRM<-df2$TN/df2$RM m.TNRM.dts <- glm(TNRM+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.TNRM.dts, main="Residual vs. fitted") summary(m.TNRM.dts) anova(m.TNRM.dts, test='F') # FNTL df2$FNTL<-df2$FN/df2$RL m.FNTL.dts <- glm(FNTL+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.FNTL.dts, main="Residual vs. fitted") summary(m.FNTL.dts) anova(m.FNTL.dts, test='F') # TNTL df2$TNTL<-df2$TN/df2$RL m.TNTL.dts <- glm(TNTL+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.TNTL.dts, main="Residual vs. fitted") summary(m.TNTL.dts) anova(m.TNTL.dts, test='F') # FNTM df2$FNTM<-df2$FN/df2$RM m.FNTM.dts <- glm(FNTM+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.FNTM.dts, main="Residual vs. fitted") summary(m.FNTM.dts) anova(m.FNTM.dts, test='F') # TNTM df2$TNTM<-df2$TN/df2$RM m.TNTM.dts <- glm(TNTM+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.TNTM.dts, main="Residual vs. fitted") summary(m.TNTM.dts) anova(m.TNTM.dts, test='F') # DNRL m.DNRL.dts <- glm(DNRL+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.DNRL.dts, main="Residual vs. fitted") summary(m.DNRL.dts) anova(m.DNRL.dts, test='F') # DNRM m.DNRM.dts <- glm(DNRM+0.001~species*crowdist+radist, data=df2, family = Gamma(link = 'log')) plot(m.DNRM.dts, main="Residual vs. fitted") summary(m.DNRM.dts) anova(m.DNRM.dts, test='F') anova(m.SL.dts) # species * # crowdist ns # radist * # spe:crow ns anova(m.RL.dts) # all ns anova(m.TL.dts) # species . # other ns anova(m.SRLR.dts, test='F') anova(m.SM.dts, test='F') anova(m.RM.dts, test='F') anova(m.TM.dts, test='F') anova(m.SRMR.dts, test='F') anova(m.FN.dts, test='Chisq') anova(m.DN.dts, test='Chisq') anova(m.TN.dts, test='Chisq') anova(m.NS.dts, test='Chisq') anova(m.FS.dts, test='Chisq') anova(m.FP.dts, test='Chisq') anova(m.DS.dts, test='Chisq') anova(m.DP.dts, test='Chisq') ### comparing models with source as identity or as disatnce ### require(MASS) require(stargazer) stargazer(m.SL, m.SL.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.SL<-anova(m.SL, test='F') afss.SL <- af.SL$"Sum Sq" round(print(cbind(af.SL,PctExp=afss.SL/sum(afss.SL)*100)), digits = 4) af.SL.dts<-anova(m.SL.dts, test='F') afss.SL.dts <- af.SL.dts$"Sum Sq" round(print(cbind(af.SL.dts,PctExp=afss.SL.dts/sum(afss.SL.dts)*100)), digits = 4) round(AIC(m.SL, m.SL.dts), 0) summary(m.SL) summary(m.SL.dts) aggregate(df2$SL, list(Name=df2$species), mean) aggregate(df2[df2$species=="DAPU",]$SL, list(Name=df2[df2$species=="DAPU",]$source), mean) TukeyHSD(aov(df2$SL~df2$species*df2$source)) # DAPU-DACA -0.5352061 -0.9513841 -0.119028 0.0118553 # DAPU:PEMB-DACA:PEMB -1.46486486 -2.89802837 -0.03170136 0.0405001 # DAPU:KWD-DACA:PEMB -1.70097598 -3.26093284 -0.14101911 0.0203304 # DAPU:KWD-DACA:GLYN -1.60111111 -3.13194766 -0.07027456 0.0320650 stargazer(m.RL, m.RL.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.RL<-anova(m.RL, test='F') afss.RL <- af.RL$"Sum Sq" round(print(cbind(af.RL,PctExp=afss.RL/sum(afss.RL)*100)), digits = 4) af.RL.dts<-anova(m.RL.dts, test='F') afss.RL.dts <- af.RL.dts$"Sum Sq" round(print(cbind(af.RL.dts,PctExp=afss.RL.dts/sum(afss.RL.dts)*100)), digits = 4) round(AIC(m.RL, m.RL.dts), 0) stargazer(m.TL, m.TL.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.TL<-anova(m.TL, test='F') afss.TL <- af.TL$"Sum Sq" round(print(cbind(af.TL,PctExp=afss.TL/sum(afss.TL)*100)), digits = 4) af.TL.dts<-anova(m.TL.dts, test='F') afss.TL.dts <- af.TL.dts$"Sum Sq" round(print(cbind(af.TL.dts,PctExp=afss.TL.dts/sum(afss.TL.dts)*100)), digits = 4) round(AIC(m.TL, m.TL.dts), 0) stargazer(m.SRLR, m.SRLR.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.SRLR<-anova(m.SRLR, test='F') afss.SRLR <- af.SRLR$`Resid. Dev` round(print(cbind(af.SRLR,PctExp=(max(afss.SRLR)-afss.SRLR)/max(afss.SRLR)*100)), digits = 4) af.SRLR.dts<-anova(m.SRLR.dts, test='F') afss.SRLR.dts <- af.SRLR.dts$`Resid. Dev` round(print(cbind(af.SRLR.dts,PctExp=(max(afss.SRLR.dts)-afss.SRLR.dts)/max(afss.SRLR.dts)*100)), digits = 4) round(AIC(m.SRLR, m.SRLR.dts), 0) stargazer(m.SM, m.SM.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.SM<-anova(m.SM, test='F') afss.SM <- af.SM$`Resid. Dev` round(print(cbind(af.SM,PctExp=(max(afss.SM)-afss.SM)/max(afss.SM)*100)), digits = 4) af.SM.dts<-anova(m.SM.dts, test='F') afss.SM.dts <- af.SM.dts$`Resid. Dev` round(print(cbind(af.SM.dts,PctExp=(max(afss.SM.dts)-afss.SM.dts)/max(afss.SM.dts)*100)), digits = 4) round(AIC(m.SM, m.SM.dts), 0) summary(m.SM) summary(m.SM.dts) aggregate(df2$SM, list(Name=df2$species), mean) aggregate(df2[df2$species=="DAPU",]$SM, list(Name=df2[df2$species=="DAPU",]$source), mean) TukeyHSD(aov(df2$SM~df2$species*df2$source)) # DAPU-DACA -6.433593 -9.632728 -3.234457 9.15e-05 # CHIP-PEMB -6.70635362 -13.357111 -0.05559648 0.0469833 # DACA:PEMB-DACA:BLUE 13.8920000 0.630833 27.153167 0.0315370 # DACA:PEMB-DAPU:BLUE 15.1100000 3.426546 26.793454 0.0019091 # DAPU:PEMB-DACA:PEMB -18.6019231 -29.618565 -7.585281 0.0000059 # DACA:CHIP-DACA:PEMB -16.8750000 -28.866295 -4.883705 0.0004185 # DAPU:CHIP-DACA:PEMB -18.0901961 -29.151648 -7.028744 0.0000140 # DAPU:GLYN-DACA:PEMB -15.7852941 -27.954034 -3.616554 0.0018160 # DACA:KWD-DACA:PEMB -12.7651163 -24.251060 -1.279173 0.0162664 # DAPU:KWD-DACA:PEMB -14.5805556 -26.571850 -2.589261 0.0049614 stargazer(m.RM, m.RM.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.RM<-anova(m.RM, test='F') afss.RM <- af.RM$`Resid. Dev` round(print(cbind(af.RM,PctExp=(max(afss.RM)-afss.RM)/max(afss.RM)*100)), digits = 4) af.RM.dts<-anova(m.RM.dts, test='F') afss.RM.dts <- af.RM.dts$`Resid. Dev` round(print(cbind(af.RM.dts,PctExp=(max(afss.RM.dts)-afss.RM.dts)/max(afss.RM.dts)*100)), digits = 4) round(AIC(m.RM, m.RM.dts), 0) stargazer(m.TM, m.TM.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.TM<-anova(m.TM, test='F') afss.TM <- af.TM$`Resid. Dev` round(print(cbind(af.TM,PctExp=(max(afss.TM)-afss.TM)/max(afss.TM)*100)), digits = 4) af.TM.dts<-anova(m.TM.dts, test='F') afss.TM.dts <- af.TM.dts$`Resid. Dev` round(print(cbind(af.TM.dts,PctExp=(max(afss.TM.dts)-afss.TM.dts)/max(afss.TM.dts)*100)), digits = 4) round(AIC(m.TM, m.TM.dts), 0) summary(m.TM) summary(m.TM.dts) aggregate(df2$TM, list(Name=df2$species), mean) aggregate(df2[df2$species=="DAPU",]$TM, list(Name=df2[df2$species=="DAPU",]$source), mean) TukeyHSD(aov(df2$TM~df2$species*df2$source)) # DAPU-DACA -9.162478 -17.80924 -0.5157123 0.0378732 # DAPU:PEMB-DACA:PEMB -33.2566528 -63.0329199 -3.4803857 0.0153263 # DACA:CHIP-DACA:PEMB -36.7519520 -69.1625541 -4.3413498 0.0127200 # DAPU:CHIP-DACA:PEMB -32.1611023 -62.0584848 -2.2637198 0.0236923 # DAPU:GLYN-DACA:PEMB -32.9973768 -65.8875860 -0.1071676 0.0484910 stargazer(m.SRMR, m.SRMR.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.SRMR<-anova(m.SRMR, test='F') afss.SRMR <- af.SRMR$`Resid. Dev` round(print(cbind(af.SRMR,PctExp=(max(afss.SRMR)-afss.SRMR)/max(afss.SRMR)*100)), digits = 4) af.SRMR.dts<-anova(m.SRMR.dts, test='F') afss.SRMR.dts <- af.SRMR.dts$`Resid. Dev` round(print(cbind(af.SRMR.dts,PctExp=(max(afss.SRMR.dts)-afss.SRMR.dts)/max(afss.SRMR.dts)*100)), digits = 4) round(AIC(m.SRMR, m.SRMR.dts), 0) summary(m.SRMR) summary(m.SRMR.dts) aggregate(df2$SRMR, list(Name=df2$species), mean) aggregate(df2[df2$species=="DAPU",]$SRMR, list(Name=df2[df2$species=="DAPU",]$source), mean) TukeyHSD(aov(df2$SRMR~df2$species*df2$source)) # DAPU-DACA -0.1740521 -0.2185865 -0.1295177 0 # GLYN-BLUE 0.13288585 0.02849378 0.23727792 0.0048892 # GLYN-PEMB 0.15498773 0.05837918 0.25159628 0.0001382 # GLYN-CHIP 0.14936836 0.05225699 0.24647973 0.0002992 # KWD-GLYN -0.16929339 -0.26863906 -0.06994772 0.0000408 # DACA:GLYN-DACA:BLUE 0.235450000 0.053657195 0.41724280 0.0018640 # DACA:GLYN-DAPU:BLUE 0.245750000 0.086307163 0.40519284 0.0000600 # DAPU:PEMB-DACA:PEMB -0.205306653 -0.358666788 -0.05194652 0.0010586 # DAPU:CHIP-DACA:PEMB -0.185808161 -0.339792091 -0.03182423 0.0055384 # DAPU:KWD-DACA:PEMB -0.237507508 -0.404435561 -0.07057945 0.0003370 # DACA:GLYN-DAPU:PEMB 0.343826923 0.193864566 0.49378928 0.0000000 # DACA:KWD-DAPU:PEMB 0.198600179 0.051624351 0.34557601 0.0008966 # DAPU:CHIP-DACA:CHIP -0.171911765 -0.327130200 -0.01669333 0.0169476 # DAPU:KWD-DACA:CHIP -0.223611111 -0.391678619 -0.05554360 0.0011799 # DACA:GLYN-DAPU:CHIP 0.324328431 0.173728205 0.47492866 0.0000000 # DACA:KWD-DAPU:CHIP 0.179101687 0.031475085 0.32672829 0.0051249 # DAPU:GLYN-DACA:GLYN -0.177955882 -0.344284285 -0.01162748 0.0251530 # DAPU:KWD-DACA:GLYN -0.376027778 -0.539839721 -0.21221583 0.0000000 # DAPU:KWD-DAPU:GLYN -0.198071895 -0.368593073 -0.02755072 0.0093407 # DAPU:KWD-DACA:KWD -0.230801034 -0.391883429 -0.06971864 0.0002928 stargazer(m.FN, m.FN.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.FN<-anova(m.FN, test='Chisq') afss.FN <- af.FN$`Resid. Dev` round(print(cbind(af.FN,PctExp=(max(afss.FN)-afss.FN)/max(afss.FN)*100)), digits = 4) af.FN.dts<-anova(m.FN.dts, test='Chisq') afss.FN.dts <- af.FN.dts$`Resid. Dev` round(print(cbind(af.FN.dts,PctExp=(max(afss.FN.dts)-afss.FN.dts)/max(afss.FN.dts)*100)), digits = 4) round(AIC(m.FN, m.FN.dts), 0) summary(m.FN) summary(m.FN.dts) aggregate(df2$FN, list(Name=df2$species), mean) aggregate(df2[df2$species=="DAPU",]$FN, list(Name=df2[df2$species=="DAPU",]$source), mean) TukeyHSD(aov(df2$FN~df2$species*df2$source)) # GLYN-BLUE -0.99067829 -1.7215189 -0.2598376 0.0021612 # DACA:GLYN-DACA:BLUE -1.54000000 -2.8127171 -0.26728287 0.0053235 # DAPU:GLYN-DACA:BLUE -1.47529412 -2.7904940 -0.16009425 0.0145418 stargazer(m.DN, m.DN.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.DN<-anova(m.DN, test='Chisq') afss.DN <- af.DN$`Resid. Dev` round(print(cbind(af.DN,PctExp=(max(afss.DN)-afss.DN)/max(afss.DN)*100)), digits = 4) af.DN.dts<-anova(m.DN.dts, test='Chisq') afss.DN.dts <- af.DN.dts$`Resid. Dev` round(print(cbind(af.DN.dts,PctExp=(max(afss.DN.dts)-afss.DN.dts)/max(afss.DN.dts)*100)), digits = 4) round(AIC(m.DN, m.DN.dts), 0) stargazer(m.TN, m.TN.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.TN<-anova(m.TN, test='Chisq') afss.TN <- af.TN$`Resid. Dev` round(print(cbind(af.TN,PctExp=(max(afss.TN)-afss.TN)/max(afss.TN)*100)), digits = 4) af.TN.dts<-anova(m.TN.dts, test='Chisq') afss.TN.dts <- af.TN.dts$`Resid. Dev` round(print(cbind(af.TN.dts,PctExp=(max(afss.TN.dts)-afss.TN.dts)/max(afss.TN.dts)*100)), digits = 4) round(AIC(m.TN, m.TN.dts), 0) summary(m.TN) summary(m.TN.dts) TukeyHSD(aov(df2$FN~df2$species*df2$source)) # GLYN-BLUE -0.99067829 -1.7215189 -0.2598376 0.0021612 # DACA:GLYN-DACA:BLUE -1.54000000 -2.8127171 -0.26728287 0.0053235 # DAPU:GLYN-DACA:BLUE -1.47529412 -2.7904940 -0.16009425 0.0145418 stargazer(m.NS, m.NS.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.NS<-anova(m.NS, test='Chisq') afss.NS <- af.NS$`Resid. Dev` round(print(cbind(af.NS,PctExp=(max(afss.NS)-afss.NS)/max(afss.NS)*100)), digits = 4) af.NS.dts<-anova(m.NS.dts, test='Chisq') afss.NS.dts <- af.NS.dts$`Resid. Dev` round(print(cbind(af.NS.dts,PctExp=(max(afss.NS.dts)-afss.NS.dts)/max(afss.NS.dts)*100)), digits = 4) round(AIC(m.NS, m.NS.dts), 0) stargazer(m.FNRL, m.FNRL.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.FNRL<-anova(m.FNRL, test='F') afss.FNRL <- af.FNRL$`Resid. Dev` round(print(cbind(af.FNRL,PctExp=(max(afss.FNRL)-afss.FNRL)/max(afss.FNRL)*100)), digits = 4) af.FNRL.dts<-anova(m.FNRL.dts, test='F') afss.FNRL.dts <- af.FNRL.dts$`Resid. Dev` round(print(cbind(af.FNRL.dts,PctExp=(max(afss.FNRL.dts)-afss.FNRL.dts)/max(afss.FNRL.dts)*100)), digits = 4) round(AIC(m.FNRL, m.FNRL.dts), 0) stargazer(m.FNRM, m.FNRM.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.FNRM<-anova(m.FNRM, test='F') afss.FNRM <- af.FNRM$`Resid. Dev` round(print(cbind(af.FNRM,PctExp=(max(afss.FNRM)-afss.FNRM)/max(afss.FNRM)*100)), digits = 4) af.FNRM.dts<-anova(m.FNRM.dts, test='F') afss.FNRM.dts <- af.FNRM.dts$`Resid. Dev` round(print(cbind(af.FNRM.dts,PctExp=(max(afss.FNRM.dts)-afss.FNRM.dts)/max(afss.FNRM.dts)*100)), digits = 4) round(AIC(m.FNRM, m.FNRM.dts), 0) stargazer(m.DNRL, m.DNRL.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.DNRL<-anova(m.DNRL, test='F') afss.DNRL <- af.DNRL$`Resid. Dev` round(print(cbind(af.DNRL,PctExp=(max(afss.DNRL)-afss.DNRL)/max(afss.DNRL)*100)), digits = 4) af.DNRL.dts<-anova(m.DNRL.dts, test='F') afss.DNRL.dts <- af.DNRL.dts$`Resid. Dev` round(print(cbind(af.DNRL.dts,PctExp=(max(afss.DNRL.dts)-afss.DNRL.dts)/max(afss.DNRL.dts)*100)), digits = 4) round(AIC(m.DNRL, m.DNRL.dts), 0) stargazer(m.DNRM, m.DNRM.dts, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") af.DNRM<-anova(m.DNRM, test='F') afss.DNRM <- af.DNRM$`Resid. Dev` round(print(cbind(af.DNRM,PctExp=(max(afss.DNRM)-afss.DNRM)/max(afss.DNRM)*100)), digits = 4) af.DNRM.dts<-anova(m.DNRM.dts, test='F') afss.DNRM.dts <- af.DNRM.dts$`Resid. Dev` round(print(cbind(af.DNRM.dts,PctExp=(max(afss.DNRM.dts)-afss.DNRM.dts)/max(afss.DNRM.dts)*100)), digits = 4) round(AIC(m.DNRM, m.DNRM.dts), 0) summary(m.DNRM) summary(m.DNRM.dts) TukeyHSD(aov(df2$DNRM~df2$species*df2$source)) # no significant differences ############################ # 2.1.3 bayesian estimates # ############################ # compute estimates by populations and by quadrants for each trait # for trait estimates # carefull: success and probability are given as logit -> need to transform # carefull: values from models with link = log need to be exponentiated # function to convert logit in prob logit2prob <- function(logit){ odds <- exp(logit) prob <- odds / (1 + odds) return(prob) } # function to reconstruct estimates for each pop extractlogit<-function(model,column){ DACABLUE=fixef(model)[1,column] DACACHIP=fixef(model)[1,column]+fixef(model)[3,column] DACAGLYNN=fixef(model)[1,column]+fixef(model)[4,column] DACAKWD=fixef(model)[1,column]+fixef(model)[5,column] DACAPEMB=fixef(model)[1,column]+fixef(model)[6,column] DAPUBLUE=fixef(model)[1,column]+fixef(model)[2,column] DAPUCHIP=fixef(model)[1,column]+fixef(model)[2,column]+fixef(model)[3,column]+fixef(model)[7,column] DAPUGLYNN=fixef(model)[1,column]+fixef(model)[2,column]+fixef(model)[4,column]+fixef(model)[8,column] DAPUKWD=fixef(model)[1,column]+fixef(model)[2,column]+fixef(model)[5,column]+fixef(model)[9,column] DAPUPEMB=fixef(model)[1,column]+fixef(model)[2,column]+fixef(model)[6,column]+fixef(model)[10,column] lo<-c(DACABLUE,DACACHIP,DACAGLYNN,DACAKWD,DACAPEMB,DAPUBLUE,DAPUCHIP,DAPUGLYNN,DAPUKWD,DAPUPEMB) names(lo)<-c("DACABLUE","DACACHIP","DACAGLYNN","DACAKWD","DACAPEMB","DAPUBLUE","DAPUCHIP","DAPUGLYNN","DAPUKWD","DAPUPEMB") pro<-logit2prob(lo) return(pro) } # function to compile probs by population + 95% CI compileprobs<-function(model){ prob<-extractlogit(model,1) error<-extractlogit(model,2) lowci<-extractlogit(model,3) hici<-extractlogit(model,4) tot<-cbind(prob,error,lowci,hici) return(tot) } fNS<-round(compileprobs(m.NSbayes),2) fNS aggregate(df2$NS, list(Name=df2$species.source), mean) # check OK fFS<-round(compileprobs(m.FSbayes),2) fFS aggregate(df2[df2$NS>0,]$FS, list(Name=df2[df2$NS>0,]$species.source), mean) # check OK fFP<-round(compileprobs(m.FPbayes),2) fFP aggregate(df2[df2$NS>0,]$FP, list(Name=df2[df2$NS>0,]$species.source), mean) # check OK fDS<-round(compileprobs(m.DSbayes),2) fDS aggregate(df2[df2$NS>0,]$DS, list(Name=df2[df2$NS>0,]$species.source), mean) # check OK fDP<-round(compileprobs(m.DPbayes),2) fDP aggregate(df2[df2$NS>0,]$DP, list(Name=df2[df2$NS>0,]$species.source), mean) # check OK # reconstruct estimates for all other traits extractest<-function(model,column){ DACABLUE=fixef(model)[1,column] DACACHIP=fixef(model)[1,column]+fixef(model)[3,column] DACAGLYNN=fixef(model)[1,column]+fixef(model)[4,column] DACAKWD=fixef(model)[1,column]+fixef(model)[5,column] DACAPEMB=fixef(model)[1,column]+fixef(model)[6,column] DAPUBLUE=fixef(model)[1,column]+fixef(model)[2,column] DAPUCHIP=fixef(model)[1,column]+fixef(model)[2,column]+fixef(model)[3,column]+fixef(model)[7,column] DAPUGLYNN=fixef(model)[1,column]+fixef(model)[2,column]+fixef(model)[4,column]+fixef(model)[8,column] DAPUKWD=fixef(model)[1,column]+fixef(model)[2,column]+fixef(model)[5,column]+fixef(model)[9,column] DAPUPEMB=fixef(model)[1,column]+fixef(model)[2,column]+fixef(model)[6,column]+fixef(model)[10,column] lo<-c(DACABLUE,DACACHIP,DACAGLYNN,DACAKWD,DACAPEMB,DAPUBLUE,DAPUCHIP,DAPUGLYNN,DAPUKWD,DAPUPEMB) names(lo)<-c("DACABLUE","DACACHIP","DACAGLYNN","DACAKWD","DACAPEMB","DAPUBLUE","DAPUCHIP","DAPUGLYNN","DAPUKWD","DAPUPEMB") est<-lo return(est) } # compile estimates by pop + 95% CI compilestim<-function(model){ est<-extractest(model,1) error<-extractest(model,2) lowci<-extractest(model,3) hici<-extractest(model,4) tot<-cbind(est,error,lowci,hici) return(tot) } fSL<-round(compilestim(m.SLbayes),2) fSL aggregate(df2$SL, list(Name=df2$species.source), mean) # check OK fRL<-round(compilestim(m.RLbayes),2) fRL aggregate(df2$RL, list(Name=df2$species.source), mean) # check OK fTL<-round(compilestim(m.TLbayes),2) fTL aggregate(df2$TL, list(Name=df2$species.source), mean) # check OK fSRLR<-round(exp(compilestim(m.SRLRbayes)),2) fSRLR aggregate(df2$SRLR, list(Name=df2$species.source), mean) # check OK fSM<-round(exp(compilestim(m.SMbayes)),2) fSM aggregate(df2$SM, list(Name=df2$species.source), mean) # check OK fRM<-round(exp(compilestim(m.RMbayes)),2) fRM aggregate(df2$RM, list(Name=df2$species.source), mean) # check OK fTM<-round(exp(compilestim(m.TMbayes)),2) fTM aggregate(df2$TM, list(Name=df2$species.source), mean) # check OK fSRMR<-round(exp(compilestim(m.SRMRbayes)),2) fSRMR aggregate(df2$SRMR, list(Name=df2$species.source), mean) # check OK fFN<-round(exp(compilestim(m.FNbayes)),2) fFN aggregate(df2$FN, list(Name=df2$species.source), mean) # check OK fDN<-round(exp(compilestim(m.DNbayes)),2) fDN aggregate(df2$DN, list(Name=df2$species.source), mean) # check OK fTN<-round(exp(compilestim(m.TNbayes)),2) fTN aggregate(df2$TN, list(Name=df2$species.source), mean) # check OK fix<-cbind(fSL,fRL,fTL,fSRLR,fSM,fRM,fTM,fSRMR,fFN,fDN,fTN,fNS,fFS,fFP,fDS,fDP)[,c(1,5,9,13,17,21,25,29,33,37,41,45,49,53,57,61)] fix<-as.data.frame(fix) colnames(fix)<-c("fSL","fRL","fTL","fSRLR","fSM","fRM","fTM","fSRMR","fFN","fDN","fTN","fNS","fFS","fFP","fDS","fDP") fix # for random effect (quadrant) estimates (marginal conditions) rSL<-round(as.data.frame(ranef(m.SLbayes)$quadrant),2) rRL<-round(as.data.frame(ranef(m.RLbayes)$quadrant),2) rTL<-round(as.data.frame(ranef(m.TLbayes)$quadrant),2) rSRLR<-round(exp(as.data.frame(coef(m.SRLRbayes))[1])-exp(fixef(m.SRLRbayes)[1]),2) rSM<-round(exp(as.data.frame(coef(m.SMbayes))[1])-exp(fixef(m.SMbayes)[1]),2) rRM<-round(exp(as.data.frame(coef(m.RMbayes))[1])-exp(fixef(m.RMbayes)[1]),2) rTM<-round(exp(as.data.frame(coef(m.TMbayes))[1])-exp(fixef(m.TMbayes)[1]),2) rSRMR<-round(exp(as.data.frame(coef(m.SRMRbayes))[1])-exp(fixef(m.SRMRbayes)[1]),2) rFN<-round(exp(as.data.frame(coef(m.FNbayes))[1])-exp(fixef(m.FNbayes)[1]),2) rDN<-round(exp(as.data.frame(coef(m.DNbayes))[1])-exp(fixef(m.DNbayes)[1]),2) rTN<-round(exp(as.data.frame(coef(m.TNbayes))[1])-exp(fixef(m.TNbayes)[1]),2) rNS<-round(logit2prob(as.data.frame(coef(m.NSbayes))[1])-logit2prob(fixef(m.NSbayes)[1]),2) rFS<-round(logit2prob(as.data.frame(coef(m.FSbayes))[1])-logit2prob(fixef(m.FSbayes)[1]),2) rFP<-round(logit2prob(as.data.frame(coef(m.FPbayes))[1])-logit2prob(fixef(m.FPbayes)[1]),2) rDS<-round(logit2prob(as.data.frame(coef(m.DSbayes))[1])-logit2prob(fixef(m.DSbayes)[1]),2) rDP<-round(logit2prob(as.data.frame(coef(m.DPbayes))[1])-logit2prob(fixef(m.DPbayes)[1]),2) rand<-cbind(rSL[1],rRL[1],rTL[1],rSRLR[1],rSM[1],rRM[1],rTM[1],rSRMR[1],rFN[1],rDN[1],rTN[1],rNS[1],rFS[1],rFP[1],rDS[1],rDP[1]) colnames(rand)<-c("fSL","fRL","fTL","fSRLR","fSM","fRM","fTM","fSRMR","fFN","fDN","fTN","fNS","fFS","fFP","fDS","fDP") rand ######################################## # 2.2 regression by distance to source # ######################################## m.SL.dapu <- lm(SL~crowdist, data=df2[df2$species=="DAPU",]) plot(m.SL.dapu, 1) summary(m.SL.dapu) anova(m.SL.dapu, test='F') # crowdist ns m.SL.daca <- lm(SL~crowdist, data=df2[df2$species=="DACA",]) plot(m.SL.daca, 1) summary(m.SL.daca) anova(m.SL.daca, test='F') # crowdist ns m.SL.all <- lmer(SL~crowdist+(1|species), data=df2) plot(m.SL.all) summary(m.SL.all) anova(m.SL.all, test='F') ranef(m.SL.all) fixef(m.SL.all) stargazer(m.SL.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist ns when controlling for species and radist # plot and regression by species x<-seq(0,500,10) plot(SL~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="SL") y.srmr.dapu<-predict(m.SL.dapu, list(crowdist=x), type="response") lines(x, y.srmr.dapu, col="blue", lwd=2) points(SL~crowdist, data=df2[df2$species=="DACA",], col="red") y.srmr.daca<-predict(m.SL.daca, list(crowdist=x), type="response") lines(x, y.srmr.daca, col="red", lwd=2) m.RL.dapu <- lm(RL~crowdist, data=df2[df2$species=="DAPU",]) plot(m.RL.dapu, 1) summary(m.RL.dapu) anova(m.RL.dapu, test='F') # crowdist ns m.RL.daca <- lm(RL~crowdist, data=df2[df2$species=="DACA",]) plot(m.RL.daca, 1) summary(m.RL.daca) anova(m.RL.daca, test='F') # crowdist ns m.RL.all <- lmer(RL~crowdist+(1|species), data=df2) plot(m.RL.all) summary(m.RL.all) anova(m.RL.all, test='F') stargazer(m.RL.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist ns when controlling for species and radist # plot and regression by species x<-seq(0,500,10) plot(RL~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="RL") y.srmr.dapu<-predict(m.RL.dapu, list(crowdist=x), type="response") lines(x, y.srmr.dapu, col="blue", lwd=2) points(RL~crowdist, data=df2[df2$species=="DACA",], col="red") y.srmr.daca<-predict(m.RL.daca, list(crowdist=x), type="response") lines(x, y.srmr.daca, col="red", lwd=2) m.TL.dapu <- lm(TL~crowdist, data=df2[df2$species=="DAPU",]) plot(m.TL.dapu, 1) summary(m.TL.dapu) anova(m.TL.dapu, test='F') # crowdist ns m.TL.daca <- lm(TL~crowdist, data=df2[df2$species=="DACA",]) plot(m.TL.daca, 1) summary(m.TL.daca) anova(m.TL.daca, test='F') # crowdist ns m.TL.all <- lmer(TL~crowdist+(1|species), data=df2) plot(m.TL.all) summary(m.TL.all) anova(m.TL.all, test='F') stargazer(m.TL.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist ns when controlling for species and radist # plot and regression by species x<-seq(0,500,10) plot(TL~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="TL") y.srmr.dapu<-predict(m.TL.dapu, list(crowdist=x), type="response") lines(x, y.srmr.dapu, col="blue", lwd=2) points(TL~crowdist, data=df2[df2$species=="DACA",], col="red") y.srmr.daca<-predict(m.TL.daca, list(crowdist=x), type="response") lines(x, y.srmr.daca, col="red", lwd=2) m.SRLR.dapu <- glm(SRLR~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.SRLR.dapu, 1) summary(m.SRLR.dapu) anova(m.SRLR.dapu, test='F') # crowdist ns m.SRLR.daca <- glm(SRLR~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.SRLR.daca, 1) summary(m.SRLR.daca) anova(m.SRLR.daca, test='F') # crowdist ns m.SRLR.all <- glmer(SRLR~crowdist+(1|species), data=df2, family = Gamma(link = 'log')) plot(m.SRLR.all) summary(m.SRLR.all) anova(m.SRLR.all, test='F') stargazer(m.SRLR.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist ns when controlling for species and radist # plot and regression by species x<-seq(0,500,10) plot(SRLR~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="SRLR") y.srmr.dapu<-predict(m.SRLR.dapu, list(crowdist=x), type="response") lines(x, y.srmr.dapu, col="blue", lwd=2) points(SRLR~crowdist, data=df2[df2$species=="DACA",], col="red") y.srmr.daca<-predict(m.SRLR.daca, list(crowdist=x), type="response") lines(x, y.srmr.daca, col="red", lwd=2) m.SM.dapu <- glm(SM~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.SM.dapu, 1) summary(m.SM.dapu) anova(m.SM.dapu, test='F') # crowdist ns m.SM.daca <- glm(SM~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.SM.daca, 1) summary(m.SM.daca) anova(m.SM.daca, test='F') # crowdist ns m.SM.all <- glmer(SM~crowdist+(1|species), data=df2, family = Gamma(link = 'log')) plot(m.SM.all) summary(m.SM.all) anova(m.SM.all, test='F') stargazer(m.SM.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist ns when controlling for species and radist # plot and regression by species x<-seq(0,500,10) plot(SM~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="SM") y.srmr.dapu<-predict(m.SM.dapu, list(crowdist=x), type="response") lines(x, y.srmr.dapu, col="blue", lwd=2) points(SM~crowdist, data=df2[df2$species=="DACA",], col="red") y.srmr.daca<-predict(m.SM.daca, list(crowdist=x), type="response") lines(x, y.srmr.daca, col="red", lwd=2) m.RM.dapu <- glm(RM~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.RM.dapu, 1) summary(m.RM.dapu) anova(m.RM.dapu, test='F') # crowdist ns m.RM.daca <- glm(RM~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.RM.daca, 1) summary(m.RM.daca) anova(m.RM.daca, test='F') # crowdist ns m.RM.all <- glmer(RM~crowdist+(1|species), data=df2, family = Gamma(link = 'log')) plot(m.RM.all) summary(m.RM.all) anova(m.RM.all, test='F') stargazer(m.RM.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist ns when controlling for species and radist # plot and regression by species x<-seq(0,500,10) plot(RM~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="RM") y.srmr.dapu<-predict(m.RM.dapu, list(crowdist=x), type="response") lines(x, y.srmr.dapu, col="blue", lwd=2) points(RM~crowdist, data=df2[df2$species=="DACA",], col="red") y.srmr.daca<-predict(m.RM.daca, list(crowdist=x), type="response") lines(x, y.srmr.daca, col="red", lwd=2) m.TM.dapu <- glm(TM~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.TM.dapu, 1) summary(m.TM.dapu) anova(m.TM.dapu, test='F') # crowdist ns m.TM.daca <- glm(TM~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.TM.daca, 1) summary(m.TM.daca) anova(m.TM.daca, test='F') # crowdist ns m.TM.all <- glmer(TM~crowdist+(1|species), data=df2, family = Gamma(link = 'log')) plot(m.TM.all) summary(m.TM.all) anova(m.TM.all, test='F') stargazer(m.TM.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist ns when controlling for species and radist # plot and regression by species x<-seq(0,500,10) plot(TM~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="TM") y.srmr.dapu<-predict(m.TM.dapu, list(crowdist=x), type="response") lines(x, y.srmr.dapu, col="blue", lwd=2) points(TM~crowdist, data=df2[df2$species=="DACA",], col="red") y.srmr.daca<-predict(m.TM.daca, list(crowdist=x), type="response") lines(x, y.srmr.daca, col="red", lwd=2) m.SRMR.dapu <- glm(SRMR~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.SRMR.dapu, 1) summary(m.SRMR.dapu) anova(m.SRMR.dapu, test='F') # crowdist * m.SRMR.daca <- glm(SRMR~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.SRMR.daca, 1) summary(m.SRMR.daca) anova(m.SRMR.daca, test='F') # crowdist ** m.SRMR.all <- glmer(SRMR~crowdist+(1|species), data=df2, family = Gamma(link = 'log')) plot(m.SRMR.all) summary(m.SRMR.all) anova(m.SRMR.all, test='F') stargazer(m.SRMR.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist ns when controlling for species and radist # plot and regression by species x<-seq(0,500,10) plot(SRMR~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="SRMR") y.srmr.dapu<-predict(m.SRMR.dapu, list(crowdist=x), type="response") lines(x, y.srmr.dapu, col="blue", lwd=2) points(SRMR~crowdist, data=df2[df2$species=="DACA",], col="red") y.srmr.daca<-predict(m.SRMR.daca, list(crowdist=x), type="response") lines(x, y.srmr.daca, col="red", lwd=2) m.FN.dapu<-glm(FN~crowdist, family="poisson", data=df2[df2$species=="DAPU",]) plot(m.FN.dapu) summary(m.FN.dapu) anova(m.FN.dapu,test='Chisq') # crowdist ns m.FN.daca<-glm(FN~crowdist, family="poisson", data=df2[df2$species=="DACA",]) plot(m.FN.daca) summary(m.FN.daca) anova(m.FN.daca,test='Chisq') # crowdist *** m.FN.all<-glmer(FN~crowdist+(1|species), family="poisson", data=df2) plot(m.FN.all) summary(m.FN.all) anova(m.FN,test='Chisq') stargazer(m.FN.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist *** # plot and regression by species x<-seq(0,500,10) plot(FN~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="FN") y.FN.dapu<-predict(m.FN.dapu, list(crowdist=x), type="response") lines(x, y.FN.dapu, col="blue", lwd=2) points(FN~crowdist, data=df2[df2$species=="DACA",], col="red") y.FN.daca<-predict(m.FN.daca, list(crowdist=x), type="response") lines(x, y.FN.daca, col="red", lwd=2) m.DN.dapu<-glm(DN~crowdist, family="poisson", data=df2[df2$species=="DAPU",]) plot(m.DN.dapu) summary(m.DN.dapu) anova(m.DN.dapu,test='Chisq') # crowdist ns m.DN.daca<-glm(DN~crowdist, family="poisson", data=df2[df2$species=="DACA",]) plot(m.DN.daca) summary(m.DN.daca) anova(m.DN.daca,test='Chisq') # crowdist . m.DN.all<-glmer(DN~crowdist+(1|species), family="poisson", data=df2) plot(m.DN.all) summary(m.DN.all) anova(m.DN.all,test='Chisq') stargazer(m.DN.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist * # plot and regression by species x<-seq(0,500,10) plot(DN~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="DN") y.DN.dapu<-predict(m.DN.dapu, list(crowdist=x), type="response") lines(x, y.DN.dapu, col="blue", lwd=2) points(DN~crowdist, data=df2[df2$species=="DACA",], col="red") y.DN.daca<-predict(m.DN.daca, list(crowdist=x), type="response") lines(x, y.DN.daca, col="red", lwd=2) m.TN.dapu<-glm(TN~crowdist, family="poisson", data=df2[df2$species=="DAPU",]) plot(m.TN.dapu) summary(m.TN.dapu) anova(m.TN.dapu,test='Chisq') # crowdist ns m.TN.daca<-glm(TN~crowdist, family="poisson", data=df2[df2$species=="DACA",]) plot(m.TN.daca) summary(m.TN.daca) anova(m.TN.daca,test='Chisq') # crowdist * m.TN.all<-glmer(TN~crowdist+(1|species), family="poisson", data=df2) plot(m.TN.all) summary(m.TN.all) anova(m.TN.all,test='Chisq') stargazer(m.TN.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # crowdist *** # plot and regression by species x<-seq(0,500,10) plot(TN~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="TN") y.TN.dapu<-predict(m.TN.dapu, list(crowdist=x), type="response") lines(x, y.TN.dapu, col="blue", lwd=2) points(TN~crowdist, data=df2[df2$species=="DACA",], col="red") y.TN.daca<-predict(m.TN.daca, list(crowdist=x), type="response") lines(x, y.TN.daca, col="red", lwd=2) m.NS.dapu<-glm(NS~crowdist, family=binomial(link="logit"), data=df2[df2$species=="DAPU",]) plot(m.NS.dapu) summary(m.NS.dapu) anova(m.NS.dapu,test='Chisq') # crowdist ns m.NS.daca<-glm(NS~crowdist, family=binomial(link="logit"), data=df2[df2$species=="DACA",]) plot(m.NS.daca) summary(m.NS.daca) anova(m.NS.daca,test='Chisq') # crowdist ns m.NS.all<-glmer(NS~crowdist+(1|species), family=binomial(link="logit"), data=df2) plot(m.NS.all) summary(m.NS.all) anova(m.NS.all,test='Chisq') stargazer(m.NS.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") # plot and regression by species x<-seq(0,500,10) plot(NS~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="NS") y.NS.dapu<-predict(m.NS.dapu, list(crowdist=x), type="response") lines(x, y.NS.dapu, col="blue", lwd=2) points(NS~crowdist, data=df2[df2$species=="DACA",], col="red",) y.NS.daca<-predict(m.NS.daca, list(crowdist=x), type="response") lines(x, y.NS.daca, col="red", lwd=2) m.FS.dapu<-glm(FS~crowdist, family=binomial(link="logit"), data=df2[df2$species=="DAPU",]) plot(m.FS.dapu, 1) summary(m.FS.dapu) anova(m.FS.dapu,test='Chisq') # crowdist ns m.FS.daca<-glm(FS~crowdist, family=binomial(link="logit"), data=df2[df2$species=="DACA",]) plot(m.FS.daca, 1) summary(m.FS.daca) anova(m.FS.daca,test='Chisq') # crowdist ns # plot and regression by species x<-seq(0,500,10) plot(FS~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="FS") y.FS.dapu<-predict(m.FS.dapu, list(crowdist=x), type="response") lines(x, y.FS.dapu, col="blue", lwd=2) points(FS~crowdist, data=df2[df2$species=="DACA",], col="red",) y.FS.daca<-predict(m.FS.daca, list(crowdist=x), type="response") lines(x, y.FS.daca, col="red", lwd=2) m.DS.dapu<-glm(DS~crowdist, family=binomial(link="logit"), data=df2[df2$species=="DAPU",]) plot(m.DS.dapu, 1) summary(m.DS.dapu) anova(m.DS.dapu,test='Chisq') # crowdist ns m.DS.daca<-glm(DS~crowdist, family=binomial(link="logit"), data=df2[df2$species=="DACA",]) plot(m.DS.daca, 1) summary(m.DS.daca) anova(m.DS.daca,test='Chisq') # crowdist ns # plot and regression by species x<-seq(0,500,10) plot(DS~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="DS") y.DS.dapu<-predict(m.DS.dapu, list(crowdist=x), type="response") lines(x, y.DS.dapu, col="blue", lwd=2) points(DS~crowdist, data=df2[df2$species=="DACA",], col="red") y.DS.daca<-predict(m.DS.daca, list(crowdist=x), type="response") lines(x, y.DS.daca, col="red", lwd=2) m.FP.dapu<-glm(FP~crowdist, family=quasibinomial(link="logit"), data=df2[df2$species=="DAPU",]) plot(m.FP.dapu, 1) summary(m.FP.dapu) anova(m.FP.dapu,test='Chisq') # crowdist ns m.FP.daca<-glm(FP~crowdist, family=quasibinomial(link="logit"), data=df2[df2$species=="DACA",]) plot(m.FP.daca, 1) summary(m.FP.daca) anova(m.FP.daca,test='Chisq') # crowdist ns # plot and regression by species x<-seq(0,500,10) plot(FP~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="FP") y.FP.dapu<-predict(m.FP.dapu, list(crowdist=x), type="response") lines(x, y.FP.dapu, col="blue", lwd=2) points(FP~crowdist, data=df2[df2$species=="DACA",], col="red") y.FP.daca<-predict(m.FP.daca, list(crowdist=x), type="response") lines(x, y.FP.daca, col="red", lwd=2) m.DP.dapu<-glm(DP~crowdist, family=quasibinomial(link="logit"), data=df2[df2$species=="DAPU",]) plot(m.DP.dapu, 1) summary(m.DP.dapu) anova(m.DP.dapu,test='Chisq') # crowdist ns m.DP.daca<-glm(DP~crowdist, family=quasibinomial(link="logit"), data=df2[df2$species=="DACA",]) plot(m.DP.daca, 1) summary(m.DP.daca) anova(m.DP.daca,test='Chisq') # crowdist ns # plot and regression by species x<-seq(0,500,10) plot(DP~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="DP") y.DP.dapu<-predict(m.DP.dapu, list(crowdist=x), type="response") lines(x, y.DP.dapu, col="blue", lwd=2) points(DP~crowdist, data=df2[df2$species=="DACA",], col="red") y.DP.daca<-predict(m.DP.daca, list(crowdist=x), type="response") lines(x, y.DP.daca, col="red", lwd=2) # normalizing number of nodules per root unit m.FNRL.dapu <- glm(FNRL+0.001~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.FNRL.dapu, 1) summary(m.FNRL.dapu) anova(m.FNRL.dapu, test='F') # crowdist ns m.FNRL.daca <- glm(FNRL+0.001~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.FNRL.daca, 1) summary(m.FNRL.daca) anova(m.FNRL.daca, test='F') # crowdist *** m.FNRL.all<-glmer(FNRL+0.001~crowdist+(1|species), family=Gamma(link = 'log'), data=df2) plot(m.FNRL.all) summary(m.FNRL.all) anova(m.FNRL.all,test='Chisq') stargazer(m.FNRL.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") m.FNRM.dapu <- glm(FNRM+0.001~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.FNRM.dapu, 1) summary(m.FNRM.dapu) anova(m.FNRM.dapu, test='F') # crowdist ns m.FNRM.daca <- glm(FNRM+0.001~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.FNRM.daca, 1) summary(m.FNRM.daca) anova(m.FNRM.daca, test='F') # crowdist ns m.FNRM.all<-glmer(FNRM+0.001~crowdist+(1|species), family=Gamma(link = 'log'), data=df2) plot(m.FNRM.all) summary(m.FNRM.all) anova(m.FNRM.all,test='Chisq') stargazer(m.FNRM.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") m.DNRL.dapu <- glm(DNRL+0.001~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.DNRL.dapu, 1) summary(m.DNRL.dapu) anova(m.DNRL.dapu, test='F') # crowdist ns m.DNRL.daca <- glm(DNRL+0.001~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.DNRL.daca, 1) summary(m.DNRL.daca) anova(m.DNRL.daca, test='F') # crowdist almost m.DNRL.all<-glmer(DNRL+0.001~crowdist+(1|species), family=Gamma(link = 'log'), data=df2) plot(m.DNRL.all) summary(m.DNRL.all) anova(m.DNRL.all,test='Chisq') stargazer(m.DNRL.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") m.DNRM.dapu <- glm(DNRM+0.001~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.DNRM.dapu, 1) summary(m.DNRM.dapu) anova(m.DNRM.dapu, test='F') # crowdist ns m.DNRM.daca <- glm(DNRM+0.001~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.DNRM.daca, 1) summary(m.DNRM.daca) anova(m.DNRM.daca, test='F') # crowdist ns m.DNRM.all<-glmer(DNRM+0.001~crowdist+(1|species), family=Gamma(link = 'log'), data=df2) plot(m.DNRM.all) summary(m.DNRM.all) anova(m.DNRM.all,test='Chisq') stargazer(m.DNRM.all, type = "text", digits = 3, star.cutoffs = c(0.05, 0.01, 0.001), digit.separator = "") m.FNTL.dapu <- glm(((FN/TL)+0.001)~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.FNTL.dapu, 1) summary(m.FNTL.dapu) anova(m.FNTL.dapu, test='F') # crowdist ns m.FNTL.daca <- glm(((FN/TL)+0.001)~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.FNTL.daca, 1) summary(m.FNTL.daca) anova(m.FNTL.daca, test='F') # crowdist *** m.FNTM.dapu <- glm(((FN/TM)+0.001)~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.FNTM.dapu, 1) summary(m.FNTM.dapu) anova(m.FNTM.dapu, test='F') # crowdist ns m.FNTM.daca <- glm(((FN/TM)+0.001)~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.FNTM.daca, 1) summary(m.FNTM.daca) anova(m.FNTM.daca, test='F') # crowdist ns m.TNRL.dapu <- glm(((TN/RL)+0.001)~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.TNRL.dapu, 1) summary(m.TNRL.dapu) anova(m.TNRL.dapu, test='F') # crowdist ns m.TNRL.daca <- glm(((TN/RL)+0.001)~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.TNRL.daca, 1) summary(m.TNRL.daca) anova(m.TNRL.daca, test='F') # crowdist *** m.TNRM.dapu <- glm(((TN/RM)+0.001)~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.TNRM.dapu, 1) summary(m.TNRM.dapu) anova(m.TNRM.dapu, test='F') # crowdist ns m.TNRM.daca <- glm(((TN/RM)+0.001)~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.TNRM.daca, 1) summary(m.TNRM.daca) anova(m.TNRM.daca, test='F') # crowdist ns m.TNTL.dapu <- glm(((TN/TL)+0.001)~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.TNTL.dapu, 1) summary(m.TNTL.dapu) anova(m.TNTL.dapu, test='F') # crowdist ns m.TNTL.daca <- glm(((TN/TL)+0.001)~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.TNTL.daca, 1) summary(m.TNTL.daca) anova(m.TNTL.daca, test='F') # crowdist *** m.TNTM.dapu <- glm(((TN/TM)+0.001)~crowdist, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) plot(m.TNTM.dapu, 1) summary(m.TNTM.dapu) anova(m.TNTM.dapu, test='F') # crowdist ns m.TNTM.daca <- glm(((TN/TM)+0.001)~crowdist, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) plot(m.TNTM.daca, 1) summary(m.TNTM.daca) anova(m.TNTM.daca, test='F') # crowdist ns # plot and regression by species x<-seq(0,500,10) plot(FN/TM~crowdist, data=df2[df2$species=="DAPU",], col="blue", main="FN/TM") y.sTMr.dapu<-predict(m.FNTM.dapu, list(crowdist=x), type="response") lines(x, y.sTMr.dapu, col="blue", lwd=2) points(FN/TM~crowdist, data=df2[df2$species=="DACA",], col="red") y.sTMr.daca<-predict(m.FNTM.daca, list(crowdist=x), type="response") lines(x, y.sTMr.daca, col="red", lwd=2) @@@ # other plotting method install.packages("gdtools") install.packages("devtools") devtools::install_github("davidgohel/ggiraph") install.packages("ggiraph") install.packages("ggiraphExtra") library(ggiraphExtra) ggPredict(m.NS.dapu,se=TRUE,interactive=TRUE,digits=3) @@@ ################################################################# # 2.3 spatial autocorrelation of trait variation across TV plot # ################################################################# # local distance (radial distance from lower left corner of TV plot, in m) was found significant for SL, SM, RM, TM, DN, and NS (and DS on nodulated plants only) # variables stored in columns 28, 32, 33, 34, 37 and 39 # for each quadrant, we have: # geographic coordinates (in cm) fom lower left corner of experimental plot # estimates of each plant and rhizobia traits, # records of disturbance by deer, ants or gophers # records of competition with weeds # Do the geographic distance (in the experiment) and/or the disturbance/competition explains distances in plant traits among quadrants? # computing dissimilarity matrices # at cell resolution # for euclidian distance between cells, euc.dist.cell euc.dist.cell<-dist(df2[,12:13]) summary(euc.dist.cell) # each trait separately, for both species val.dist.cell.SL<-dist(df2[,28]) val.dist.cell.RL<-dist(df2[,29]) val.dist.cell.TL<-dist(df2[,30]) val.dist.cell.SRLR<-dist(df2[,31]) val.dist.cell.SM<-dist(df2[,32]) val.dist.cell.RM<-dist(df2[,33]) val.dist.cell.TM<-dist(df2[,34]) val.dist.cell.SRMR<-dist(df2[,35]) val.dist.cell.FN<-dist(df2[,36]) val.dist.cell.DN<-dist(df2[,37]) val.dist.cell.TN<-dist(df2[,38]) val.dist.cell.NS<-dist(df2[,39]) val.dist.cell.FNRL<-dist(df2[,48]) val.dist.cell.TNRL<-dist(df2[,49]) val.dist.cell.FNRM<-dist(df2[,50]) val.dist.cell.TNRM<-dist(df2[,51]) val.dist.cell.FNTL<-dist(df2[,52]) val.dist.cell.TNTL<-dist(df2[,53]) val.dist.cell.FNTM<-dist(df2[,54]) val.dist.cell.TNTM<-dist(df2[,55]) val.dist.cell.DNRL<-dist(df2[,56]) val.dist.cell.DNRM<-dist(df2[,57]) # mantel tests to test for sparial autocorrelation # at the cell resolution on actual trait values # whithout affording for species and source # SL mantel.euc.val.cell.SL<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.SL, method="s", permutations=999) mantel.euc.val.cell.SL # overall **, r=0.0401, p=0.003 mant.cor.SL<-mantel.correlog(val.dist.cell.SL, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.SL, main="SL") # # RL mantel.euc.val.cell.RL<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.RL, method="s", permutations=999) mantel.euc.val.cell.RL # overall ns, r=0.0142, p=0.175 mant.cor.RL<-mantel.correlog(val.dist.cell.RL, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.RL, main="RL") # # TL mantel.euc.val.cell.TL<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.TL, method="s", permutations=999) mantel.euc.val.cell.TL # overall ., r=0.0239, p=0.051 mant.cor.TL<-mantel.correlog(val.dist.cell.TL, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.TL, main="TL") # # SRLR mantel.euc.val.cell.SRLR<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.SRLR, method="s", permutations=999) mantel.euc.val.cell.SRLR # overall ns, r=0.0052, p=0.371 mant.cor.SRLR<-mantel.correlog(val.dist.cell.SRLR, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.SRLR, main="SRLR") # # SM mantel.euc.val.cell.SM<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.SM, method="s", permutations=999) mantel.euc.val.cell.SM # overall **, r=0.0424, p=0.006 mant.cor.SM<-mantel.correlog(val.dist.cell.SM, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.SM, main="SM") # # RM mantel.euc.val.cell.RM<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.RM, method="s", permutations=999) mantel.euc.val.cell.RM # overall *, r=0.0369, p=0.015 mant.cor.RM<-mantel.correlog(val.dist.cell.RM, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.RM, main="RM") # # TM mantel.euc.val.cell.TM<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.TM, method="s", permutations=999) mantel.euc.val.cell.TM # overall *, r=0.0383, p=0.019 mant.cor.TM<-mantel.correlog(val.dist.cell.TM, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.TM, main="TM") # # SRMR mantel.euc.val.cell.SRMR<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.SRMR, method="s", permutations=999) mantel.euc.val.cell.SRMR # overall ns, r=0.0091, p=0.273 mant.cor.SRMR<-mantel.correlog(val.dist.cell.SRMR, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.SRMR, main="SRMR") # # FN mantel.euc.val.cell.FN<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.FN, method="s", permutations=999) mantel.euc.val.cell.FN # overall ns, r=0.0113, p=0.246 mant.cor.FN<-mantel.correlog(val.dist.cell.FN, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.FN, main="FN") # # DN mantel.euc.val.cell.DN<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.DN, method="s", permutations=999) mantel.euc.val.cell.DN # overall *, r=0.0334, p=0.039 mant.cor.DN<-mantel.correlog(val.dist.cell.DN, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.DN, main="DN") # # TN mantel.euc.val.cell.TN<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.TN, method="s", permutations=999) mantel.euc.val.cell.TN # overall ns, r=0.0081, p=0.283 mant.cor.TN<-mantel.correlog(val.dist.cell.TN, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.TN, main="TN") # # NS mantel.euc.val.cell.NS<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.NS, method="s", permutations=999) mantel.euc.val.cell.NS # overall ns, r=-0.0269, p=0.98 mant.cor.NS<-mantel.correlog(val.dist.cell.NS, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.NS, main="NS") # # FS mantel.euc.val.cell.FS<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.FS, method="s", permutations=999) mantel.euc.val.cell.FS # nd mant.cor.FS<-mantel.correlog(val.dist.cell.FS, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.FS, main="FS") # nd # FP mantel.euc.val.cell.FP<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.FP, method="s", permutations=999) mantel.euc.val.cell.FP # nd mant.cor.FP<-mantel.correlog(val.dist.cell.FP, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.FP, main="FP") # nd # DS mantel.euc.val.cell.DS<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.DS, method="s", permutations=999) mantel.euc.val.cell.DS # nd mant.cor.DS<-mantel.correlog(val.dist.cell.DS, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.DS, main="DS") # nd # DP mantel.euc.val.cell.DP<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.DP, method="s", permutations=999) mantel.euc.val.cell.DP # nd mant.cor.DP<-mantel.correlog(val.dist.cell.DP, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.DP, main="DP") # nd # FNRL mantel.euc.val.cell.FNRL<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.FNRL, method="s", permutations=999) mantel.euc.val.cell.FNRL # overall ns, r=0.0040, p=0.379 mant.cor.FNRL<-mantel.correlog(val.dist.cell.FNRL, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.FNRL, main="FNRL") # # FNRM mantel.euc.val.cell.FNRM<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.FNRM, method="s", permutations=999) mantel.euc.val.cell.FNRM # overall nd, r=0.0002, p=0.487 mant.cor.FNRM<-mantel.correlog(val.dist.cell.FNRM, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.FNRM, main="FNRM") # # DNRL mantel.euc.val.cell.DNRL<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.DNRL, method="s", permutations=999) mantel.euc.val.cell.DNRL # mant.cor.DNRL<-mantel.correlog(val.dist.cell.DNRL, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.DNRL, main="DNRL") # # DNRM mantel.euc.val.cell.DNRM<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.DNRM, method="s", permutations=999) mantel.euc.val.cell.DNRM # mant.cor.DNRM<-mantel.correlog(val.dist.cell.DNRM, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.DNRM, main="DNRM") # # FNTL mantel.euc.val.cell.FNTL<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.FNTL, method="s", permutations=999) mantel.euc.val.cell.FNTL # overall ns, r=0.003995, p=0.385 mant.cor.FNTL<-mantel.correlog(val.dist.cell.FNTL, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.FNTL, main="FNTL") # # FNTM mantel.euc.val.cell.FNTM<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.FNTM, method="s", permutations=999) mantel.euc.val.cell.FNTM # overall ns, r=0.0002159, p=0.485 mant.cor.FNTM<-mantel.correlog(val.dist.cell.FNTM, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.FNTM, main="FNTM") # # TNRL mantel.euc.val.cell.TNRL<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.TNRL, method="s", permutations=999) mantel.euc.val.cell.TNRL # overall ns, r=-0.002241, p=0.541 mant.cor.TNRL<-mantel.correlog(val.dist.cell.TNRL, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.TNRL, main="TNRL") # # TNRM mantel.euc.val.cell.TNRM<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.TNRM, method="s", permutations=999) mantel.euc.val.cell.TNRM # overall ns, r=0.001076, p=0.431 mant.cor.TNRM<-mantel.correlog(val.dist.cell.TNRM, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.TNRM, main="TNRM") # # TNTL mantel.euc.val.cell.TNTL<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.TNTL, method="s", permutations=999) mantel.euc.val.cell.TNTL # overall ns, r=-0.002241, p=0.528 mant.cor.TNTL<-mantel.correlog(val.dist.cell.TNTL, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.TNTL, main="TNTL") # # TNTM mantel.euc.val.cell.TNTM<-mantel(xdis=euc.dist.cell, ydis=val.dist.cell.TNTM, method="s", permutations=999) mantel.euc.val.cell.TNTM # overall ns, r=0.001076, p=0.447 mant.cor.TNTM<-mantel.correlog(val.dist.cell.TNTM, euc.dist.cell, r.type="s", cutoff=F) plot(mant.cor.TNTM, main="TNTM") # # at the cell resolution on actual trait values # separately for each species # DACA euc.dist.cell.daca<-dist(df2[df2$species=="DACA",12:13]) summary(euc.dist.cell.daca) val.dist.cell.daca.SL<-dist(df2[df2$species=="DACA",28]) val.dist.cell.daca.RL<-dist(df2[df2$species=="DACA",29]) val.dist.cell.daca.TL<-dist(df2[df2$species=="DACA",30]) val.dist.cell.daca.SRLR<-dist(df2[df2$species=="DACA",31]) val.dist.cell.daca.SM<-dist(df2[df2$species=="DACA",32]) val.dist.cell.daca.RM<-dist(df2[df2$species=="DACA",33]) val.dist.cell.daca.TM<-dist(df2[df2$species=="DACA",34]) val.dist.cell.daca.SRMR<-dist(df2[df2$species=="DACA",35]) val.dist.cell.daca.FN<-dist(df2[df2$species=="DACA",36]) val.dist.cell.daca.DN<-dist(df2[df2$species=="DACA",37]) val.dist.cell.daca.TN<-dist(df2[df2$species=="DACA",38]) val.dist.cell.daca.NS<-dist(df2[df2$species=="DACA",39]) val.dist.cell.daca.FNRL<-dist(df2[df2$species=="DACA",48]) val.dist.cell.daca.TNRL<-dist(df2[df2$species=="DACA",49]) val.dist.cell.daca.FNRM<-dist(df2[df2$species=="DACA",50]) val.dist.cell.daca.TNRM<-dist(df2[df2$species=="DACA",51]) val.dist.cell.daca.FNTL<-dist(df2[df2$species=="DACA",52]) val.dist.cell.daca.TNTL<-dist(df2[df2$species=="DACA",53]) val.dist.cell.daca.FNTM<-dist(df2[df2$species=="DACA",54]) val.dist.cell.daca.TNTM<-dist(df2[df2$species=="DACA",55]) val.dist.cell.daca.DNRL<-dist(df2[df2$species=="DACA",56]) val.dist.cell.daca.DNRM<-dist(df2[df2$species=="DACA",57]) # SL mantel.euc.val.cell.daca.SL<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.SL, method="s", permutations=999) mantel.euc.val.cell.daca.SL # overall ns, r=0.02795, p=0.078 mant.cor.daca.SL<-mantel.correlog(val.dist.cell.daca.SL, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.SL, main="SL") # significant spatial autocorrelation at 1.6 m # RL mantel.euc.val.cell.daca.RL<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.RL, method="s", permutations=999) mantel.euc.val.cell.daca.RL # overall ns, r=0.02178, p=0.17 mant.cor.daca.RL<-mantel.correlog(val.dist.cell.daca.RL, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.RL, main="RL") # significant spatial autocorrelation at 11, 17, 27 and 40 m, but no clear pattern # TL mantel.euc.val.cell.daca.TL<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.TL, method="s", permutations=999) mantel.euc.val.cell.daca.TL # overall ns, r=0.01306, p=0.247 mant.cor.daca.TL<-mantel.correlog(val.dist.cell.daca.TL, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.TL, main="TL") # significant spatial autocorrelation at 1.6, 4.8, 11 and 27 m, but no clear pattern # SRLR mantel.euc.val.cell.daca.SRLR<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.SRLR, method="s", permutations=999) mantel.euc.val.cell.daca.SRLR # overall ns, r=-0.003824, p=0.532 mant.cor.daca.SRLR<-mantel.correlog(val.dist.cell.daca.SRLR, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.SRLR, main="SRLR") # significant spatial autocorrelation at 1.6, and 11 m # SM mantel.euc.val.cell.daca.SM<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.SM, method="s", permutations=999) mantel.euc.val.cell.daca.SM # overall ns, r=0.03232, p=0.119 mant.cor.daca.SM<-mantel.correlog(val.dist.cell.daca.SM, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.SM, main="SM") # significant spatial autocorrelation at 1.6, 4.8, 8, 14 and 17 m # RM mantel.euc.val.cell.daca.RM<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.RM, method="s", permutations=999) mantel.euc.val.cell.daca.RM # overall *, r=0.04494, p=0.046 mant.cor.daca.RM<-mantel.correlog(val.dist.cell.daca.RM, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.RM, main="RM") # significant spatial autocorrelation at 1.6, 4.8, 8 and 14 m # TM mantel.euc.val.cell.daca.TM<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.TM, method="s", permutations=999) mantel.euc.val.cell.daca.TM # overall ., r=0.04048, p=0.056 mant.cor.daca.TM<-mantel.correlog(val.dist.cell.daca.TM, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.TM, main="TM") # # SRMR mantel.euc.val.cell.daca.SRMR<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.SRMR, method="s", permutations=999) mantel.euc.val.cell.daca.SRMR # ns, r=0.02698, p=0.126 mant.cor.daca.SRMR<-mantel.correlog(val.dist.cell.daca.SRMR, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.SRMR, main="SRMR") # significant spatial autocorrelation at 1.6, 4.8, 8 and 14 m # FN mantel.euc.val.cell.daca.FN<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.FN, method="s", permutations=999) mantel.euc.val.cell.daca.FN # overall ., r=0.03766, p=0.059 mant.cor.daca.FN<-mantel.correlog(val.dist.cell.daca.FN, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.FN, main="FN") # no significant spatial autocorrelation at any scale # DN mantel.euc.val.cell.daca.DN<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.DN, method="s", permutations=999) mantel.euc.val.cell.daca.DN # overall ns, r=-5.922e-05, p=0.491 mant.cor.daca.DN<-mantel.correlog(val.dist.cell.daca.DN, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.DN, main="DN") # significant spatial autocorrelation at 1.6 m # TN mantel.euc.val.cell.daca.TN<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.TN, method="s", permutations=999) mantel.euc.val.cell.daca.TN # ns, r=0.02144, p=0.18 mant.cor.daca.TN<-mantel.correlog(val.dist.cell.daca.TN, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.TN, main="TN") # no significant spatial autocorrelation at any scale # NS mantel.euc.val.cell.daca.NS<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.NS, method="s", permutations=999) mantel.euc.val.cell.daca.NS # overall ns, r=-0.04164, p=0.987 mant.cor.daca.NS<-mantel.correlog(val.dist.cell.daca.NS, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.NS, main="NS") # significant spatial autocorrelation at 4.8 and 30 m # FS mantel.euc.val.cell.daca.FS<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.FS, method="s", permutations=999) mantel.euc.val.cell.daca.FS # mant.cor.daca.FS<-mantel.correlog(val.dist.cell.daca.FS, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.FS, main="FS") # # FP mantel.euc.val.cell.daca.FP<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.FP, method="s", permutations=999) mantel.euc.val.cell.daca.FP # mant.cor.daca.FP<-mantel.correlog(val.dist.cell.daca.FP, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.FP, main="FP") # # DS mantel.euc.val.cell.daca.DS<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.DS, method="s", permutations=999) mantel.euc.val.cell.daca.DS # mant.cor.daca.DS<-mantel.correlog(val.dist.cell.daca.DS, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.DS, main="DS") # # DP mantel.euc.val.cell.daca.DP<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.DP, method="s", permutations=999) mantel.euc.val.cell.daca.DP # mant.cor.daca.DP<-mantel.correlog(val.dist.cell.daca.DP, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.DP, main="DP") # # FNRL mantel.euc.val.cell.daca.FNRL<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.FNRL, method="s", permutations=999) mantel.euc.val.cell.daca.FNRL # overall ., r=0.03167, p=0.094 mant.cor.daca.FNRL<-mantel.correlog(val.dist.cell.daca.FNRL, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.FNRL, main="FNRL") # no significant spatial autocorrelation at any scale # FNRM mantel.euc.val.cell.daca.FNRM<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.FNRM, method="s", permutations=999) mantel.euc.val.cell.daca.FNRM # overall *, r=0.0518, p=0.028 mant.cor.daca.FNRM<-mantel.correlog(val.dist.cell.daca.FNRM, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.FNRM, main="FNRM") # no significant spatial autocorrelation at any scale # DNRL mantel.euc.val.cell.daca.DNRL<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.DNRL, method="s", permutations=999) mantel.euc.val.cell.daca.DNRL # overall ns mant.cor.daca.DNRL<-mantel.correlog(val.dist.cell.daca.DNRL, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.DNRL, main="DNRL") # no significant spatial autocorrelation at any scale # DNRM mantel.euc.val.cell.daca.DNRM<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.DNRM, method="s", permutations=999) mantel.euc.val.cell.daca.DNRM # overall ns mant.cor.daca.DNRM<-mantel.correlog(val.dist.cell.daca.DNRM, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.DNRM, main="DNRM") # no significant spatial autocorrelation at any scale # FNTL mantel.euc.val.cell.daca.FNTL<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.FNTL, method="s", permutations=999) mantel.euc.val.cell.daca.FNTL # overall ., r=0.03167, p=0.095 mant.cor.daca.FNTL<-mantel.correlog(val.dist.cell.daca.FNTL, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.FNTL, main="FNTL") # no significant spatial autocorrelation at any scale # FNTM mantel.euc.val.cell.daca.FNTM<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.FNTM, method="s", permutations=999) mantel.euc.val.cell.daca.FNTM # overall *, r=0.0518, p=0.028 mant.cor.daca.FNTM<-mantel.correlog(val.dist.cell.daca.FNTM, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.FNTM, main="FNTM") # no significant spatial autocorrelation at any scale # TNRL mantel.euc.val.cell.daca.TNRL<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.TNRL, method="s", permutations=999) mantel.euc.val.cell.daca.TNRL # overall ns, r=0.006372, p=0.392 mant.cor.daca.TNRL<-mantel.correlog(val.dist.cell.daca.TNRL, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.TNRL, main="TNRL") # no significant spatial autocorrelation at any scale # TNRM mantel.euc.val.cell.daca.TNRM<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.TNRM, method="s", permutations=999) mantel.euc.val.cell.daca.TNRM # overall ns, r=0.02376, p=0.172 mant.cor.daca.TNRM<-mantel.correlog(val.dist.cell.daca.TNRM, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.TNRM, main="TNRM") # significant spatial autocorrelation at 20 m # TNTL mantel.euc.val.cell.daca.TNTL<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.TNTL, method="s", permutations=999) mantel.euc.val.cell.daca.TNTL # overall ns, r=0.006372, p=0.417 mant.cor.daca.TNTL<-mantel.correlog(val.dist.cell.daca.TNTL, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.TNTL, main="TNTL") # no significant spatial autocorrelation at any scale # TNTM mantel.euc.val.cell.daca.TNTM<-mantel(xdis=euc.dist.cell.daca, ydis=val.dist.cell.daca.TNTM, method="s", permutations=999) mantel.euc.val.cell.daca.TNTM # overall ns, r=0.02376, p=0.172 mant.cor.daca.TNTM<-mantel.correlog(val.dist.cell.daca.TNTM, euc.dist.cell.daca, r.type="s", cutoff=F) plot(mant.cor.daca.TNTM, main="TNTM") # # DAPU euc.dist.cell.dapu<-dist(df2[df2$species=="DAPU",12:13]) summary(euc.dist.cell.dapu) val.dist.cell.dapu.SL<-dist(df2[df2$species=="DAPU",28]) val.dist.cell.dapu.RL<-dist(df2[df2$species=="DAPU",29]) val.dist.cell.dapu.TL<-dist(df2[df2$species=="DAPU",30]) val.dist.cell.dapu.SRLR<-dist(df2[df2$species=="DAPU",31]) val.dist.cell.dapu.SM<-dist(df2[df2$species=="DAPU",32]) val.dist.cell.dapu.RM<-dist(df2[df2$species=="DAPU",33]) val.dist.cell.dapu.TM<-dist(df2[df2$species=="DAPU",34]) val.dist.cell.dapu.SRMR<-dist(df2[df2$species=="DAPU",35]) val.dist.cell.dapu.FN<-dist(df2[df2$species=="DAPU",36]) val.dist.cell.dapu.DN<-dist(df2[df2$species=="DAPU",37]) val.dist.cell.dapu.TN<-dist(df2[df2$species=="DAPU",38]) val.dist.cell.dapu.NS<-dist(df2[df2$species=="DAPU",39]) val.dist.cell.dapu.FNRL<-dist(df2[df2$species=="DAPU",48]) val.dist.cell.dapu.TNRL<-dist(df2[df2$species=="DAPU",49]) val.dist.cell.dapu.FNRM<-dist(df2[df2$species=="DAPU",50]) val.dist.cell.dapu.TNRM<-dist(df2[df2$species=="DAPU",51]) val.dist.cell.dapu.FNTL<-dist(df2[df2$species=="DAPU",52]) val.dist.cell.dapu.TNTL<-dist(df2[df2$species=="DAPU",53]) val.dist.cell.dapu.FNTM<-dist(df2[df2$species=="DAPU",54]) val.dist.cell.dapu.TNTM<-dist(df2[df2$species=="DAPU",55]) val.dist.cell.dapu.DNRL<-dist(df2[df2$species=="DAPU",56]) val.dist.cell.dapu.DNRM<-dist(df2[df2$species=="DAPU",57]) # SL mantel.euc.val.cell.dapu.SL<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.SL, method="s", permutations=999) mantel.euc.val.cell.dapu.SL # overall **, r=0.05428, p=0.008 mant.cor.dapu.SL<-mantel.correlog(val.dist.cell.dapu.SL, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.SL, main="SL") # significant spatial autocorrelation at 1.6, 4.8 and 8 m # RL mantel.euc.val.cell.dapu.RL<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.RL, method="s", permutations=999) mantel.euc.val.cell.dapu.RL # overall ns, r=0.0009191, p=0.468 mant.cor.dapu.RL<-mantel.correlog(val.dist.cell.dapu.RL, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.RL, main="RL") # significant spatial autocorrelation at 11 m # TL mantel.euc.val.cell.dapu.TL<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.TL, method="s", permutations=999) mantel.euc.val.cell.dapu.TL # overall ., r=0.03037, p=0.077 mant.cor.dapu.TL<-mantel.correlog(val.dist.cell.dapu.TL, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.TL, main="TL") # significant spatial autocorrelation at 4.8 m # SRLR mantel.euc.val.cell.dapu.SRLR<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.SRLR, method="s", permutations=999) mantel.euc.val.cell.dapu.SRLR # overall ns, r=0.0104, p=0.276 mant.cor.dapu.SRLR<-mantel.correlog(val.dist.cell.dapu.SRLR, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.SRLR, main="SRLR") # significant spatial autocorrelation at 11 m # SM mantel.euc.val.cell.dapu.SM<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.SM, method="s", permutations=999) mantel.euc.val.cell.dapu.SM # overall *, r=0.05068, p=0.026 mant.cor.dapu.SM<-mantel.correlog(val.dist.cell.dapu.SM, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.SM, main="SM") # significant spatial autocorrelation at 1.6, 4.8, 8, 14 and 17 m # RM mantel.euc.val.cell.dapu.RM<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.RM, method="s", permutations=999) mantel.euc.val.cell.dapu.RM # overall ., r=0.03126, p=0.077 mant.cor.dapu.RM<-mantel.correlog(val.dist.cell.dapu.RM, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.RM, main="RM") # significant spatial autocorrelation at 1.6 and 4.8 m # TM mantel.euc.val.cell.dapu.TM<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.TM, method="s", permutations=999) mantel.euc.val.cell.dapu.TM # overall ., r=0.03699, p=0.059 mant.cor.dapu.TM<-mantel.correlog(val.dist.cell.dapu.TM, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.TM, main="TM") # significant spatial autocorrelation at 1.6, 4.8, 8 and 17 m # SRMR mantel.euc.val.cell.dapu.SRMR<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.SRMR, method="s", permutations=999) mantel.euc.val.cell.dapu.SRMR # overall ns, r=-0.0137, p=0.737 mant.cor.dapu.SRMR<-mantel.correlog(val.dist.cell.dapu.SRMR, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.SRMR, main="SRMR") # significant spatial autocorrelation at 30 m # FN mantel.euc.val.cell.dapu.FN<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.FN, method="s", permutations=999) mantel.euc.val.cell.dapu.FN # overall ns, r=-0.009593, p=0.632 mant.cor.dapu.FN<-mantel.correlog(val.dist.cell.dapu.FN, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.FN, main="FN") # no significant spatial autocorrelation at any scale # DN mantel.euc.val.cell.dapu.DN<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.DN, method="s", permutations=999) mantel.euc.val.cell.dapu.DN # overall *, r=0.06056, p=0.012 mant.cor.dapu.DN<-mantel.correlog(val.dist.cell.dapu.DN, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.DN, main="DN") # significant spatial autocorrelation at 1.6, 4.8 and 8 m # TN mantel.euc.val.cell.dapu.TN<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.TN, method="s", permutations=999) mantel.euc.val.cell.dapu.TN # overall ns, r=-0.00653, p=0.626 mant.cor.dapu.TN<-mantel.correlog(val.dist.cell.dapu.TN, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.TN, main="TN") # no significant spatial autocorrelation at any scale # NS mantel.euc.val.cell.dapu.NS<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.NS, method="s", permutations=999) mantel.euc.val.cell.dapu.NS # overall ns, r=0., p=0.889 mant.cor.dapu.NS<-mantel.correlog(val.dist.cell.dapu.NS, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.NS, main="NS") # significant spatial autocorrelation at 1.6 m # FS mantel.euc.val.cell.dapu.FS<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.FS, method="s", permutations=999) mantel.euc.val.cell.dapu.FS # nd mant.cor.dapu.FS<-mantel.correlog(val.dist.cell.dapu.FS, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.FS, main="FS") # nd # FP mantel.euc.val.cell.dapu.FP<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.FP, method="s", permutations=999) mantel.euc.val.cell.dapu.FP # nd mant.cor.dapu.FP<-mantel.correlog(val.dist.cell.dapu.FP, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.FP, main="FP") # nd # DS mantel.euc.val.cell.dapu.DS<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.DS, method="s", permutations=999) mantel.euc.val.cell.dapu.DS # nd mant.cor.dapu.DS<-mantel.correlog(val.dist.cell.dapu.DS, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.DS, main="DS") # nd # DP mantel.euc.val.cell.dapu.DP<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.DP, method="s", permutations=999) mantel.euc.val.cell.dapu.DP # nd mant.cor.dapu.DP<-mantel.correlog(val.dist.cell.dapu.DP, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.DP, main="DP") # nd # FNRL mantel.euc.val.cell.dapu.FNRL<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.FNRL, method="s", permutations=999) mantel.euc.val.cell.dapu.FNRL # overall ns, r=-0.02373, p=0.858 mant.cor.dapu.FNRL<-mantel.correlog(val.dist.cell.dapu.FNRL, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.FNRL, main="FNRL") # no significant spatial autocorrelation at any scale # FNRM mantel.euc.val.cell.dapu.FNRM<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.FNRM, method="s", permutations=999) mantel.euc.val.cell.dapu.FNRM # overall ns, r=-0.04383, p=0.977 mant.cor.dapu.FNRM<-mantel.correlog(val.dist.cell.dapu.FNRM, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.FNRM, main="FNRM") # no significant spatial autocorrelation at any scale # DNRL mantel.euc.val.cell.dapu.DNRL<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.DNRL, method="s", permutations=999) mantel.euc.val.cell.dapu.DNRL # overall * mant.cor.dapu.DNRL<-mantel.correlog(val.dist.cell.dapu.DNRL, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.DNRL, main="DNRL") # no significant spatial autocorrelation at any scale # DNRM mantel.euc.val.cell.dapu.DNRM<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.DNRM, method="s", permutations=999) mantel.euc.val.cell.dapu.DNRM # overall * mant.cor.dapu.DNRM<-mantel.correlog(val.dist.cell.dapu.DNRM, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.DNRM, main="DNRM") # no significant spatial autocorrelation at any scale # FNTL mantel.euc.val.cell.dapu.FNTL<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.FNTL, method="s", permutations=999) mantel.euc.val.cell.dapu.FNTL # overall ns, r=-0.02373, p=0.852 mant.cor.dapu.FNTL<-mantel.correlog(val.dist.cell.dapu.FNTL, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.FNTL, main="FNTL") # no significant spatial autocorrelation at any scale # FNTM mantel.euc.val.cell.dapu.FNTM<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.FNTM, method="s", permutations=999) mantel.euc.val.cell.dapu.FNTM # overall ns, r=-0.04383, p=0.972 mant.cor.dapu.FNTM<-mantel.correlog(val.dist.cell.dapu.FNTM, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.FNTM, main="FNTM") # no significant spatial autocorrelation at any scale # TNRL mantel.euc.val.cell.dapu.TNRL<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.TNRL, method="s", permutations=999) mantel.euc.val.cell.dapu.TNRL # overall ns, r=-0.0153, p=0.734 mant.cor.dapu.TNRL<-mantel.correlog(val.dist.cell.dapu.TNRL, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.TNRL, main="TNRL") # no significant spatial autocorrelation at any scale # TNRM mantel.euc.val.cell.dapu.TNRM<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.TNRM, method="s", permutations=999) mantel.euc.val.cell.dapu.TNRM # overall ns, r=-0.01639, p=0.73 mant.cor.dapu.TNRM<-mantel.correlog(val.dist.cell.dapu.TNRM, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.TNRM, main="TNRM") # no significant spatial autocorrelation at any scale # TNTL mantel.euc.val.cell.dapu.TNTL<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.TNTL, method="s", permutations=999) mantel.euc.val.cell.dapu.TNTL # overall ns, r=-0.0153, p=0.754 mant.cor.dapu.TNTL<-mantel.correlog(val.dist.cell.dapu.TNTL, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.TNTL, main="TNTL") # no significant spatial autocorrelation at any scale # TNTM mantel.euc.val.cell.dapu.TNTM<-mantel(xdis=euc.dist.cell.dapu, ydis=val.dist.cell.dapu.TNTM, method="s", permutations=999) mantel.euc.val.cell.dapu.TNTM # overall ns, r=-0.01639, p=0.745 mant.cor.dapu.TNTM<-mantel.correlog(val.dist.cell.dapu.TNTM, euc.dist.cell.dapu, r.type="s", cutoff=F) plot(mant.cor.dapu.TNTM, main="TNTM") # # at quadrant resolution # make sure vector 'rand'is complete and accurate # adding quadrant center coordinates in a 2-dimensional orthogonal (Cartesian) system (in meters) rand$y<-rep(c(0.5,2,7,20,34),5) rand$x<-c(rep(0.5,5),rep(2,5),rep(7,5),rep(20,5),rep(34,5)) rand$z<-sqrt((rand$x^2)+(rand$y^2)) plot(rand$y, rand$x, pch=0, cex=2, main="Spatial distribution of quadrants", xlab="Distance from lower left (SW) corner (m)", ylab="Distance from lower left (SW) corner (m)") arrows(x0=25, y0=27.5, x1=25, y1=30, length=0.1, lwd=2) text(x=26, y = 28.5, labels = "N", cex=2) text(x=rand$x, y=rand$y, labels=LETTERS[1:25], cex=0.75) # for euclidian distance between quadrants, euc.dist.quad euc.dist.quad<-dist(scale(rand[,17:18])) # for distance between quadrant estimates, all 6 traits, est.dist.quad.all6 est.dist.quad.all6<-dist(scale(rand[,c(1,5:7,10,12)])) # each trait separately est.dist.quad.SL<-dist(scale(rand[,1])) est.dist.quad.SM<-dist(scale(rand[,5])) est.dist.quad.RM<-dist(scale(rand[,6])) est.dist.quad.TM<-dist(scale(rand[,7])) est.dist.quad.DN<-dist(scale(rand[,10])) est.dist.quad.NS<-dist(scale(rand[,12])) #est.dist.quad.DS<-dist(scale(rand[,9])) # at the quadrant resolution on trait estimates, whithout affording for species and source # for 6 traits together mantel.euc.est.quad.all6<-mantel(xdis=euc.dist.quad, ydis=est.dist.quad.all6, method="s", permutations=999) mantel.euc.est.quad.all6 # ns, r=0.01657, p=0.393 # SL mantel.euc.est.quad.SL<-mantel(xdis=euc.dist.quad, ydis=est.dist.quad.SL, method="s", permutations=999) mantel.euc.est.quad.SL # ns, r=0.04049, p=0.278 # SM mantel.euc.est.quad.SM<-mantel(xdis=euc.dist.quad, ydis=est.dist.quad.SM, method="s", permutations=999) mantel.euc.est.quad.SM # ns, r=0.04328, p=0.25 # RM mantel.euc.est.quad.RM<-mantel(xdis=euc.dist.quad, ydis=est.dist.quad.RM, method="s", permutations=999) mantel.euc.est.quad.RM # ns, r=0.007782, p=0.403 # TM mantel.euc.est.quad.TM<-mantel(xdis=euc.dist.quad, ydis=est.dist.quad.TM, method="s", permutations=999) mantel.euc.est.quad.TM # ns, r=0.01495, p=0.379 # DN mantel.euc.est.quad.DN<-mantel(xdis=euc.dist.quad, ydis=est.dist.quad.DN, method="s", permutations=999) mantel.euc.est.quad.DN # ns, r=-0.0889, p=0.883 # NS mantel.euc.est.quad.NS<-mantel(xdis=euc.dist.quad, ydis=est.dist.quad.NS, method="s", permutations=999) mantel.euc.est.quad.NS # ns, r=-0.001987, p=0.453 ######################################### # 3. Evidences for local co-adaptation? # ######################################### # hyp: plant dry mass correlates with number of nodules cor(df2[,c(32:34,36:38)], method="p") cor.mtest(df2[,c(32:34,36:38)], method="p") cor(df2[df2$species=="DACA",c(32:34,36:38)], method="p") cor.mtest(df2[df2$species=="DACA",c(32:34,36:38)], method="p") cor(df2[df2$species=="DAPU",c(32:34,36:38)], method="p") cor.mtest(df2[df2$species=="DAPU",c(32:34,36:38)], method="p") # for DACA: weak but significant positive correlation between nodule number (FN and TN) and plant mass (SM, RM, TM) trait (p-values=[0.01;0.03]; Pearson rhos=[0.15-0.18]) # for DAPU: weak but significant positive correlation between nodule (FN an TN) number and plant mass (RM, TM) trait (p-values=[0.004;0.01]; Pearson rhos=[0.14-0.19]) # correlation are stronger for total number of nodules and root mass, confirmed by PCA # hyp: plant dry mass associated with success of nodulation m.SM.daca.ns <- glm(SM~NS, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) m.RM.daca.ns <- glm(RM~NS, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) m.TM.daca.ns <- glm(TM~NS, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) summary(m.SM.daca.ns) # * summary(m.RM.daca.ns) # * summary(m.TM.daca.ns) # * # for DACA, significant effect of NS on plant mass m.SM.dapu.ns <- glm(SM~NS, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) m.RM.dapu.ns <- glm(RM~NS, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) m.TM.dapu.ns <- glm(TM~NS, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) summary(m.SM.dapu.ns) # ns summary(m.RM.dapu.ns) # ns summary(m.TM.dapu.ns) # ns # for DAPU, no effect of NS on plant mass m.SM.ns <- glmer(SM~NS+(1|species)+(1|source), data=df2, family = Gamma(link = 'log')) m.RM.ns <- glmer(RM~NS+(1|species)+(1|source), data=df2, family = Gamma(link = 'log')) m.TM.ns <- glmer(TM~NS+(1|species)+(1|source), data=df2, family = Gamma(link = 'log')) summary(m.SM.ns) # . summary(m.RM.ns) # . summary(m.TM.ns) # . ranef(m.SM.ns) fixef(m.SM.ns) ranef(m.RM.ns) fixef(m.RM.ns) ranef(m.TM.ns) fixef(m.TM.ns) # almost significant effect of NS on plant mass when controlling for species an source # hyp: total mass associated with total number of nodules m.TM.daca.tn <- glm(TM~TN, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) summary(m.TM.daca.tn) # for DACA, significant effect of TN on TM m.TM.dapu.tn <- glm(TM~TN, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) summary(m.TM.dapu.tn) # for DAPU, significant effect of TN on TM m.TM.tn <- glmer(TM~TN+(1|species)+(1|source), data=df2, family = Gamma(link = 'log')) summary(m.TM.tn) # significant effect of TN on TM when controlling for species and source ranef(m.TM.tn) fixef(m.TM.tn) # hyp: total mass associated with number of fresh nodules m.TM.daca.fn <- glm(TM~FN, data=df2[df2$species=="DACA",], family = Gamma(link = 'log')) summary(m.TM.daca.fn) # for DACA, no effect of FN on TM m.TM.dapu.fn <- glm(TM~FN, data=df2[df2$species=="DAPU",], family = Gamma(link = 'log')) summary(m.TM.dapu.fn) # for DAPU, almost significant effect of FN on TM m.TM.fn <- glmer(TM~FN+(1|species)+(1|source), data=df2, family = Gamma(link = 'log')) summary(m.TM.fn) # significant effect of FN on TM when controlling for species and source ranef(m.TM.fn) fixef(m.TM.fn) # hyp: plant length associated with success of nodulation m.SL.daca.ns <- lm(SL~NS, data=df2[df2$species=="DACA",]) m.RL.daca.ns <- lm(RL~NS, data=df2[df2$species=="DACA",]) m.TL.daca.ns <- lm(TL~NS, data=df2[df2$species=="DACA",]) summary(m.SL.daca.ns) # ns summary(m.RL.daca.ns) # * summary(m.TL.daca.ns) # * # for DACA, significant effect of NS on RL and TL m.SL.dapu.ns <- lm(SL~NS, data=df2[df2$species=="DAPU",]) m.RL.dapu.ns <- lm(RL~NS, data=df2[df2$species=="DAPU",]) m.TL.dapu.ns <- lm(TL~NS, data=df2[df2$species=="DAPU",]) summary(m.SL.dapu.ns) # ns summary(m.RL.dapu.ns) # * summary(m.TL.dapu.ns) # ns # for DAPU, no effect of NS on RL only m.SL.ns <- lmer(SL~NS+(1|species)+(1|source), data=df2) m.RL.ns <- lmer(RL~NS+(1|species)+(1|source), data=df2) m.TL.ns <- lmer(TL~NS+(1|species)+(1|source), data=df2) summary(m.SL.ns) # ns summary(m.RL.ns) # ns summary(m.TL.ns) # ns ranef(m.SL.ns) fixef(m.SL.ns) ranef(m.RL.ns) fixef(m.RL.ns) ranef(m.TL.ns) fixef(m.TL.ns) # no significant effect of NS on plant length when controlling for species an source # hyp: total length associated with total number of nodules m.TL.daca.tn <- lm(TL~TN, data=df2[df2$species=="DACA",]) summary(m.TL.daca.tn) # for DACA, no significant effect of TN on TL m.TL.dapu.tn <- lm(TL~TN, data=df2[df2$species=="DAPU",]) summary(m.TL.dapu.tn) # for DAPU, almost significant effect of TN on TL # hyp: total length associated with number of fresh nodules m.TL.daca.fn <- lm(TL~FN, data=df2[df2$species=="DACA",]) summary(m.TL.daca.fn) # for DACA, significant effect of FN on TL m.TL.dapu.fn <- lm(TL~FN, data=df2[df2$species=="DAPU",]) summary(m.TL.dapu.fn) # for DAPU, almost significant effect of FN on TL m.TL.fn <- lmer(TL~FN+(1|species)+(1|source), data=df2) summary(m.TL.fn) ranef(m.TL.fn) fixef(m.TL.fn) var.test(TM~NS, data=df2[df2$species=="DACA",], alternative = "two.sided") # equal t.test(TM~NS, data=df2[df2$species=="DACA",], alternative="less", var.equal=T) # *, t = -2.233, df = 179, p-value = 0.0134 var.test(TM~NS, data=df2[df2$species=="DAPU",], alternative = "two.sided") # equal t.test(TM~NS, data=df2[df2$species=="DAPU",], alternative="less", var.equal=T) # ns, t = 0.14311, df = 211, p-value = 0.5568 var.test(TL~NS, data=df2[df2$species=="DACA",], alternative = "two.sided") # equal t.test(TL~NS, data=df2[df2$species=="DACA",], alternative="less", var.equal=T) # *, t = -2.3171, df = 179, p-value = 0.01081 var.test(TL~NS, data=df2[df2$species=="DAPU",], alternative = "two.sided") # equal t.test(TL~NS, data=df2[df2$species=="DAPU",], alternative="less", var.equal=T) # ns, t = -0.82154, df = 211, p-value = 0.2061 anova(lm(NS~species, data=df2), test="Chisq") # ns aggregate(df2$NS, list(Name=df2$species), mean) anova(lm(TN~species, data=df2), test="Chisq") # ns aggregate(df2$TN, list(Name=df2$species), mean) anova(lm(FN~species, data=df2), test="Chisq") # ns aggregate(df2$FN, list(Name=df2$species), mean) aggregate(df2$TM, list(Name=df2$NS), mean) aggregate(df2[df2$species=="DACA",]$TM, list(Name=df2[df2$species=="DACA",]$NS), mean) aggregate(df2[df2$species=="DAPU",]$TM, list(Name=df2[df2$species=="DAPU",]$NS), mean) aggregate(df2$TL, list(Name=df2$NS), mean) aggregate(df2[df2$species=="DACA",]$TL, list(Name=df2[df2$species=="DACA",]$NS), mean) aggregate(df2[df2$species=="DAPU",]$TL, list(Name=df2[df2$species=="DAPU",]$NS), mean) # hyp: DACA and DAPU differ in nodulation and fate (FN or DN) of nodules aggregate(df2$FN, list(Name=df2$species), mean) aggregate(df2$DN, list(Name=df2$species), mean) aggregate(df2$TN, list(Name=df2$species), mean) aggregate(df2$NS, list(Name=df2$species), mean) # are these difference significant? var.test(df2$FN~df2$species, alternative = "two.sided") # NOT equal variances t.test(df2$FN~df2$species, alternative = "two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means var.test(df2$DN~df2$species, alternative = "two.sided") # NOT equal variances t.test(df2$DN~df2$species, alternative = "two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means var.test(df2$TN~df2$species, alternative = "two.sided") # equal variances t.test(df2$TN~df2$species, alternative = "two.sided", mu=0, paired=F, var.equal=T, conf.level=0.95) # equal means var.test(df2$NS~df2$species, alternative = "two.sided") # equal variances t.test(df2$NS~df2$species, alternative = "two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means var.test(df2$FNRL~df2$species, alternative = "two.sided") # NOT equal variances t.test(df2$FNRL~df2$species, alternative = "two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means var.test(df2$FNRM~df2$species, alternative = "two.sided") # NOT equal variances t.test(df2$FNRM~df2$species, alternative = "two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means var.test(df2$DNRL~df2$species, alternative = "two.sided") # NOT equal variances t.test(df2$DNRL~df2$species, alternative = "two.sided", mu=0, paired=F, var.equal=T, conf.level=0.95) # equal means var.test(df2$DNRM~df2$species, alternative = "two.sided") # equal variances t.test(df2$DNRM~df2$species, alternative = "two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means # trying another approach cor.test(df2$FN, df2$DN, method="s") cor.test(df2[df2$species=='DACA',]$FN, df2[df2$species=='DACA',]$DN, method="s") cor.test(df2[df2$species=='DAPU',]$FN, df2[df2$species=='DAPU',]$DN, method="s") cor.nod.all<-cor(df2[df2$NS>0,36:38], method="s") cor.nod.daca<-cor(df2[df2$NS>0 & df2$species=='DACA',36:38], method="s") cor.nod.dapu<-cor(df2[df2$NS>0 & df2$species=='DAPU',36:38], method="s") pmat.nod.all <- cor.mtest(cor.nod.all, method="s") pmat.nod.daca <- cor.mtest(cor.nod.daca, method="s") pmat.nod.dapu <- cor.mtest(cor.nod.dapu, method="s") corrplot(cor.nod.all, method="circle", type="upper", main="Correlation of measured trait on all plants", p.mat = pmat.nod.all, sig.level = 0.05, mar = c(2,2,4,2)) corrplot(cor.nod.daca, method="circle", type="upper", main="Correlation of measured trait on DACA", p.mat = pmat.nod.daca, sig.level = 0.05, mar = c(2,2,4,2)) corrplot(cor.nod.dapu, method="circle", type="upper", main="Correlation of measured trait on DAPU", p.mat = pmat.nod.dapu, sig.level = 0.05, mar = c(2,2,4,2)) # no significant correlations between nodulation traits when species analyzed together # nor when analyzed separately ### LCCMR slides ### plot(df2[df2$species=="DACA",]$TN, df2[df2$species=="DACA",]$SM, col="red") fit.daca <- lm(df2[df2$species=="DACA",]$SM~df2[df2$species=="DACA",]$TN, data=df2[df2$species=="DACA",]) lines(df2[df2$species=="DACA",]$TN, fitted(fit.daca), col="red") points(df2[df2$species=="DAPU",]$TN, df2[df2$species=="DAPU",]$SM, col="blue") fit.dapu <- lm(df2[df2$species=="DAPU",]$SM~df2[df2$species=="DAPU",]$TN, data=df2[df2$species=="DAPU",]) lines(df2[df2$species=="DAPU",]$TN, fitted(fit.dapu), col="blue") frac.w.nod<-aggregate(df2$NS, list(Name=df2$species.source), mean) descdist(frac.w.nod[,2], discrete = F, boot = 500) # uniform dts<-aggregate(df2$crowdist, list(Name=df2$species.source), mean) dts.miles<-dts dts.miles[,2]<-round(dts[,2]/1.609344,0) dts.miles fit.daca<-lm(frac.w.nod[1:5,2]~dts.miles[1:5,2]) fit.dapu<-lm(frac.w.nod[6:10,2]~dts.miles[6:10,2]) fit.all<-lm(frac.w.nod[,2]~dts.miles[,2]) par(mar=c(5,5,2,5)) plot(dts.miles[1:5,2], frac.w.nod[1:5,2]+0.01, cex=2.5, cex.axis=2, col="black", pch=1, lwd=3.25, ylim=c(0.5,1), ylab="Fraction of plants with nodules", cex.lab=2, font.lab=2 , xlim=c(0,300), xlab="Distance to the source (miles)", cex.lab=2, font.lab=2 ) #lines(dts.miles[1:5,2], fitted(fit.daca), col="black", lwd=7.5) points(dts.miles[6:10,2]-0.01, frac.w.nod[6:10,2], cex=2.75, col="purple", pch=21) #lines(dts.miles[6:10,2], fitted(fit.dapu), col="purple", lwd=7.5) legend(x="topright", y=NULL, title="Prairie clover species", legend=c("White","Purple"), fill=c("white","purple"), cex=1.5, inset=0.02, box.lty=0) lines(dts.miles[,2], fitted(fit.all), col="black", lwd=7.5) ### equality of variances and means for nodulation traits between species ### var.test(x=df2[df2$species=="DACA",]$FN, y=df2[df2$species=="DAPU",]$FN, alternative = "two.sided") # NOT equal variances t.test(x=df2[df2$species=="DACA",]$FN, y=df2[df2$species=="DAPU",]$FN, alternative="two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # NOT equal means var.test(x=df2[df2$species=="DACA",]$DN, y=df2[df2$species=="DAPU",]$DN, alternative = "two.sided") # NOT equal variances t.test(x=df2[df2$species=="DACA",]$DN, y=df2[df2$species=="DAPU",]$DN, alternative="two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means var.test(x=df2[df2$species=="DACA",]$TN, y=df2[df2$species=="DAPU",]$TN, alternative = "two.sided") # NOT equal variances t.test(x=df2[df2$species=="DACA",]$TN, y=df2[df2$species=="DAPU",]$TN, alternative="two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means var.test(x=df2[df2$species=="DACA",]$NS, y=df2[df2$species=="DAPU",]$NS, alternative = "two.sided") # equal variances t.test(x=df2[df2$species=="DACA",]$NS, y=df2[df2$species=="DAPU",]$NS, alternative="two.sided", mu=0, paired=F, var.equal=T, conf.level=0.95) # equal means summary(as.factor(df2[df2$species=="DACA",]$NS)) summary(as.factor(df2[df2$species=="DAPU",]$NS)) prop.test(x=c(133, 151), n=c(48+133, 62+151), alternative="two.sided", correct=F) # two-proportions z-test # no differences in proportion of nodulated plants var.test(x=df2[df2$species=="DACA",]$FNRL, y=df2[df2$species=="DAPU",]$FNRL, alternative = "two.sided") # NOT equal variances t.test(x=df2[df2$species=="DACA",]$FNRL, y=df2[df2$species=="DAPU",]$FNRL, alternative="two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means (nearly different) var.test(x=df2[df2$species=="DACA",]$FNRM, y=df2[df2$species=="DAPU",]$FNRM, alternative = "two.sided") # equal variances t.test(x=df2[df2$species=="DACA",]$FNRM, y=df2[df2$species=="DAPU",]$FNRM, alternative="two.sided", mu=0, paired=F, var.equal=T, conf.level=0.95) # equal means var.test(x=df2[df2$species=="DACA",]$DNRL, y=df2[df2$species=="DAPU",]$DNRL, alternative = "two.sided") # NOT equal variances t.test(x=df2[df2$species=="DACA",]$DNRL, y=df2[df2$species=="DAPU",]$DNRL, alternative="two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means var.test(x=df2[df2$species=="DACA",]$DNRM, y=df2[df2$species=="DAPU",]$DNRM, alternative = "two.sided") # NOT equal variances t.test(x=df2[df2$species=="DACA",]$DNRM, y=df2[df2$species=="DAPU",]$DNRM, alternative="two.sided", mu=0, paired=F, var.equal=F, conf.level=0.95) # equal means ### co-evolution = correlation between plant and rhizobia traits ### cor.test(x = df2[df2$species=="DACA",]$TL, y = df2[df2$species=="DACA",]$NS, method = "s", alternative = "two.sided", exact = F, continuity = T) # **, S = 782289, p-value = 0.004869, rho = 0.208418 cor.test(x = df2[df2$species=="DACA",]$TM, y = df2[df2$species=="DACA",]$NS, method = "s", alternative = "two.sided", exact = F, continuity = T) # **, S = 791303, p-value = 0.007151, rho = 0.1992964 cor.test(x = df2[df2$species=="DAPU",]$TL, y = df2[df2$species=="DAPU",]$NS, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 1518646, p-value = 0.4073, rho = 0.05707171 cor.test(x = df2[df2$species=="DAPU",]$TM, y = df2[df2$species=="DAPU",]$NS, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 1672693, p-value = 0.5756, rho = -0.03857612 cor.test(x = df2[df2$species=="DACA",]$TL, y = df2[df2$species=="DACA",]$TN, method = "s", alternative = "two.sided", exact = F, continuity = T) # *, S = 830147, p-value = 0.03144, rho = 0.1599909 cor.test(x = df2[df2$species=="DACA",]$TM, y = df2[df2$species=="DACA",]$TN, method = "s", alternative = "two.sided", exact = F, continuity = T) # **, S = 783726, p-value = 0.005183, rho = 0.2069637 cor.test(x = df2[df2$species=="DAPU",]$TL, y = df2[df2$species=="DAPU",]$TN, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 1398584, p-value = 0.05512, rho = 0.1316187 cor.test(x = df2[df2$species=="DAPU",]$TM, y = df2[df2$species=="DAPU",]$TN, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 1464955, p-value = 0.1887, rho = 0.09040858 cor.test(x = df2[df2$species=="DACA",]$TL, y = df2[df2$species=="DACA",]$FN, method = "s", alternative = "two.sided", exact = F, continuity = T) # **, S = 766973, p-value = 0.002444 cor.test(x = df2[df2$species=="DACA",]$TM, y = df2[df2$species=="DACA",]$FN, method = "s", alternative = "two.sided", exact = F, continuity = T) # *, S = 809602, p-value = 0.01487 cor.test(x = df2[df2$species=="DAPU",]$TL, y = df2[df2$species=="DAPU",]$FN, method = "s", alternative = "two.sided", exact = F, continuity = T) # *, S = 1381871, p-value = 0.03839 cor.test(x = df2[df2$species=="DAPU",]$TM, y = df2[df2$species=="DAPU",]$FN, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 1470433, p-value = 0.206 cor.test(x = df2[df2$species=="DACA",]$TL, y = df2[df2$species=="DACA",]$DN, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 1010608, p-value = 0.7625 cor.test(x = df2[df2$species=="DACA",]$TM, y = df2[df2$species=="DACA",]$DN, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 897405, p-value = 0.2184 cor.test(x = df2[df2$species=="DAPU",]$TL, y = df2[df2$species=="DAPU",]$DN, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 1511888, p-value = 0.3736 cor.test(x = df2[df2$species=="DAPU",]$TM, y = df2[df2$species=="DAPU",]$DN, method = "s", alternative = "two.sided", exact = F, continuity = T) # ns, S = 1506845, p-value = 0.3496