# 0) Loading librairies, functions and data #### library(boral) library(reshape2) library(ggplot2) library(tidyr) library(corrplot) library(igraph) library(RColorBrewer) library(vegan) library(Cairo) library(PBSmapping) library(maps) library(data.table) library(wesanderson) library(plyr) library(gridExtra) library(zoo) library(ade4) library(packfor) library(MASS) library(ellipse) library(FactoMineR) library(adegraphics) library(adespatial) library(SoDA) library(factoextra) library(venneuler) library(ppcor) library(viridis) source("evplot.R") source("hcoplot.R") source("scalog.R") source("quickMEM.R") source("panelutils.R") source("cleanplot.pca.R") source("screestick.R") source("box.cox.chord (Appendix 4).R") source("triplot.rda.R") source("sr.value.R") source("my.circleplot.R") phyto <- read.csv("Phyto.csv") zoop <- read.csv("Zoop.csv") zoop <- zoop[,-c(15,33,38,39)] # remove CACN, CYCN = Stade copépodites (not sp level), remove ROTI = Rotifères indéterminés, remove NAUP morpho <- read.csv("Morphometry.csv") physico <- read.csv("Physicochem.csv") phyto_synt <- read.csv("Phyto_Synt.csv") phyto_synt_sub <- phyto_synt[,c("Station","N","n","Chloro")] zoop_synt <- read.csv("Zoop_Synt.csv") zoop_synt_sub <- zoop_synt[,c("Station","DZOO","DCL","DCA")] fish <- read.csv("Poisson.csv") vars_synt <- read.csv("Vars_synt.csv") # Merge data files by Station ID data_m <- merge(x=zoop, y=phyto, by="Station",all.x=TRUE) # zooplankton + phytoplankton data_m <- merge(x=data_m, y=morpho, by="Station",all.x=TRUE) # + morphometry data_m <- merge(x=data_m, y=physico, by="Station",all.x=TRUE) # + physico-chemistry data_m <- merge(x=data_m, y=phyto_synt_sub,all.x=TRUE) # + phyto synthetic data_m <- merge(x=data_m, y=zoop_synt_sub,all.x=TRUE) # + zoop synthetic data_m <- merge(x=data_m, y=fish, by="Station",all.y=TRUE) # + fish (Note: 5 sites missing in fish data) # Now that rows are re-order by Station, re-subset and transform data zoop <- data_m[,2:35] phyto <- data_m[,36:122] morpho <- data_m[,123:137] physico <- data_m[,138:163] phyto_synt <- data_m[,164:166] zoop_synt <- data_m[,167:169] fish <- data_m[,172:189] # do not include columns "CUE" and "n" # Species data transformations phyto.log.chord <- box.cox.chord(phyto,bc.exp=0); rownames(phyto.log.chord) <- data_m[,1] # log.chord phyto.hel <- box.cox.chord(phyto,bc.exp=0.5);rownames(phyto.hel) <- data_m[,1] # hellinger phyto.chord <- box.cox.chord(phyto,bc.exp=1);rownames(phyto.chord) <- data_m[,1] # chord zoop.log.chord <- box.cox.chord(zoop,bc.exp=0); rownames(zoop.log.chord) <- data_m[,1] # log.chord zoop.hel <- box.cox.chord(zoop,bc.exp=0.5);rownames(zoop.hel) <- data_m[,1] # hellinger zoop.chord <- box.cox.chord(zoop,bc.exp=1);rownames(zoop.chord) <- data_m[,1] # chord fish.log.chord <- box.cox.chord(fish,bc.exp=0); rownames(fish.log.chord) <- data_m[,1] # log.chord fish.hel <- box.cox.chord(fish,bc.exp=0.5);rownames(fish.hel) <- data_m[,1] # hellinger fish.chord <- box.cox.chord(fish,bc.exp=1);rownames(fish.chord) <- data_m[,1] # chord # Environmental data scaling env <- cbind(physico,morpho) env <- scale(env); env <- as.data.frame(env) sapply(env, function(x) any(is.na(x))) env <- env[!sapply(env, function(x) any(is.na(x)))] env <- env[,-c(18)] # Cond appears twice # 1.1) Species Association #### # Phytoplancton # Try different transformations {chord, Hellinger, log-chord} # data.transform <- box.cox.chord(phyto,bc.exp=0) # log.chord # data.transform <- box.cox.chord(phyto,bc.exp=0.5) # hellinger data.transform <- box.cox.chord(phyto,bc.exp=1) # chord # Clustering methods (K-means et Ward) data.norm <- scale(data.transform) datanorm.t <- t(data.norm) # Transpose scaled matrix # K-means partitionning kmin <- 2 # between k = 2 ... p <- ncol(data.norm) kmax <- p/2 # ... and k = p/2 where p is the number of columns in the matrix part_ch <- cascadeKM(datanorm.t, kmin, kmax, iter = 9999) part_ch plot(part_ch) # Calinski-Harabasz is maximized par the partition in 2 groups # W de Kendall - Start with all variables in the same group # H0: variables are not concordant (if p < 0.05 reject H0 and look for groups of concordant variables) out1 <- kendall.global(data.transform, nperm=9999) # The fonction transforms the data into ranks out1 # Significant: at least one species is concordant with another... but which one(s) => a posterior test needed # W de Kendall - now separate the variables into different groups # H0: this variable is not concordant (if p < 0.05, reject H0, and you've identified concordant species in the subgroup) out.post.1 <- kendall.post(data.transform, mult="holm", nperm=9999) out.post.1$A_posteriori_tests data.cor = cor(data.transform , method="spearman") # Spearman r matrix among columns data.D = as.dist(1 - data.cor) # Transform r into Distance data.clust = hclust(data.D, "ward.D2") # Groupement agglomératif de Ward par(mgp=c(2,0.6,0.5),mai=c(0.8,0.75,0.4,0)); plot(data.clust, hang=-1,cex=.5) # 2 groups formed. If it doesn't work try 3 groups, etc. data.2 = cutree(data.clust, k=4) # Cut the dendrogram ; separate species 2, 3, 4 groups etc. data.2 rect.hclust(data.clust, k=4,border="darkgreen") out.2 <- kendall.global(data.transform ,group=data.2,mult="holm") # Global test of each group out.2 out.post.2 <- kendall.post(data.transform , group=data.2, mult="holm", nperm=9999) out.post.2$A_posteriori_tests_Group[[1]] out.post.2$A_posteriori_tests_Group[[2]] out.post.2$A_posteriori_tests_Group[[3]] out.post.2$A_posteriori_tests_Group[[4]] # ACP to determine visually how many groups are needed # Examine whether correlations among these groups are positives or negatives # pal <- colorRampPalette(brewer.pal(11,"Spectral"))(11) # pal <- pal[c(2,4,8,10)] # palette(pal) pca.phyto = rda(data.transform, scale.unit=TRUE, ncp=5, graph=TRUE) summary(pca.phyto) plot(pca.phyto, display=c("sp", "lc", "cn"),type="n",xlab="PC1 (15.5%)",ylab="PC2 (12.6%)") # Scaling 2 phyto2.sc = scores(pca.phyto, choices=1:2, display="sp") arrows(0, 0, phyto2.sc[,1], phyto2.sc[,2], lty=1, lwd=1, length=0.1,col=data.2) text(phyto2.sc[,1], phyto2.sc[,2], rownames(phyto2.sc),cex=0.9,col=data.2) # Create object for groups to calculate RV coefficient and p-values below # Note: RV are R-squared - They cannot be negative phyto1 <- data.transform[,data.2 == 1] phyto2 <- data.transform[,data.2 == 2] phyto3 <- data.transform[,data.2 == 3] phyto4 <- data.transform[,data.2 == 4] # Zooplancton # data.transform <- box.cox.chord(zoop,bc.exp=0) # log.chord # data.transform <- box.cox.chord(zoop,bc.exp=0.5) # hellinger data.transform <- box.cox.chord(zoop,bc.exp=1) # chord data.norm <- scale(data.transform) datanorm.t <- t(data.norm) # K-means kmin <- 2 p <- ncol(data.norm) kmax <- p/2 part_ch <- cascadeKM(datanorm.t, kmin, kmax, iter = 9999) part_ch plot(part_ch) # Maximized at 2 for C-H (and chord transformation) # W Kendall out1 <- kendall.global(data.transform, nperm=9999) out1 # W Kendall, sépare maintenant les variables en différents groupes out.post.1 <- kendall.post(data.transform, mult="holm", nperm=9999) out.post.1$A_posteriori_tests data.cor = cor(data.transform , method="spearman") data.D = as.dist(1 - data.cor) data.clust = hclust(data.D, "ward.D2") par(mgp=c(2,0.6,0.5),mai=c(0.8,0.75,0.4,0)); plot(data.clust, hang=-1,cex=0.8) # 2 groupes formés. Si marche pas va essayer 3, etc. data.2 = cutree(data.clust, k=4) data.2 rect.hclust(data.clust, k=4,border="blue") out.2 <- kendall.global(data.transform ,group=data.2,mult="holm") out.2 out.post.2 <- kendall.post(data.transform , group=data.2, mult="holm", nperm=9999) out.post.2$A_posteriori_tests_Group[[1]] out.post.2$A_posteriori_tests_Group[[2]] out.post.2$A_posteriori_tests_Group[[3]] out.post.2$A_posteriori_tests_Group[[4]] # PCA pca.zoop = rda(data.transform, scale.unit=TRUE, ncp=5, graph=TRUE) summary(pca.zoop) plot(pca.zoop, display=c("sp", "lc", "cn"),type="n",xlab="PC1 (35%)",ylab="PC2 (15%)") # Scaling 2 zoop2.sc = scores(pca.zoop, choices=1:2, display="sp") arrows(0, 0, zoop2.sc[,1], zoop2.sc[,2], lty=1, lwd=1, length=0.1,col=data.2) text(zoop2.sc[,1], zoop2.sc[,2], rownames(zoop2.sc),cex=0.9,col=data.2) # RV Coefficients zoop1 <- data.transform[,data.2 == 1] zoop2 <- data.transform[,data.2 == 2] zoop3 <- data.transform[,data.2 == 3] zoop4 <- data.transform[,data.2 == 4] # Fish # data.transform <- box.cox.chord(fish,bc.exp=0) # log.chord # data.transform <- box.cox.chord(fish,bc.exp=0.5) # hellinger data.transform <- box.cox.chord(fish,bc.exp=1) # chord data.norm <- scale(data.transform) datanorm.t <- t(data.norm) # K-means kmin <- 2 p <- ncol(data.norm) kmax <- p/2 part_ch <- cascadeKM(datanorm.t, kmin, kmax, iter = 9999) part_ch plot(part_ch) # Maximized at 2 for C-H (and chord transformation) # W Kendall out1 <- kendall.global(data.transform, nperm=9999) out1 out.post.1 <- kendall.post(data.transform, mult="holm", nperm=9999) out.post.1$A_posteriori_tests data.cor = cor(data.transform , method="spearman") data.D = as.dist(1 - data.cor) data.clust = hclust(data.D, "ward.D2") par(mgp=c(2,0.6,0.5),mai=c(0.8,0.75,0.4,0)); plot(data.clust, hang=-1,cex=1) data.2 = cutree(data.clust, k=3) data.2 rect.hclust(data.clust, k=3,border="red") out.2 <- kendall.global(data.transform ,group=data.2,mult="holm") out.2 out.post.2 <- kendall.post(data.transform , group=data.2, mult="holm", nperm=9999) out.post.2$A_posteriori_tests_Group[[1]] out.post.2$A_posteriori_tests_Group[[2]] out.post.2$A_posteriori_tests_Group[[3]] # PCA pca.fish = rda(data.transform, scale.unit=TRUE, ncp=5, graph=TRUE) summary(pca.fish) plot(pca.fish, display=c("sp", "lc", "cn"),type="n",xlab="PC1 (35.7%)",ylab="PC2 (20%)") # Scaling 2 fish2.sc = scores(pca.fish, choices=1:2, display="sp") arrows(0, 0, fish2.sc[,1], fish2.sc[,2], lty=1, lwd=1, length=0.1,col=data.2) text(fish2.sc[,1], fish2.sc[,2], rownames(fish2.sc),cex=0.9,col=data.2) # RV Coefficients fish1 <- data.transform[,data.2 == 1] fish2 <- data.transform[,data.2 == 2] fish3 <- data.transform[,data.2 == 3] # 1.2) RV coefficients #### tab1 <- cbind(phyto1,phyto2,phyto3,phyto4, zoop1,zoop2,zoop3,zoop4, fish1,fish2,fish3) grn <- c(ncol(phyto1),ncol(phyto2),ncol(phyto3),ncol(phyto4), ncol(zoop1),ncol(zoop2),ncol(zoop3),ncol(zoop4), ncol(fish1),ncol(fish2),ncol(fish3)) QC.mfa <- MFA(tab1, group=grn, type=c("c","c","c","c", "c","c","c","c", "c","c","c"), ncp=5, name.group=c("Phyto 1","Phyto 2","Phyto 3","Phyto 4", "Zoop 1","Zoop 2","Zoop 3","Zoop 4", "Fish 1","Fish 2","Fish 3")) # Table 1: adding p-values to RV coefficient output rvp <- QC.mfa$group$RV rvp[1, 2] <- coeffRV(phyto1, phyto2)$p.value rvp[1, 3] <- coeffRV(phyto1, phyto3)$p.value rvp[1, 4] <- coeffRV(phyto1, phyto4)$p.value rvp[1, 5] <- coeffRV(phyto1, zoop1)$p.value rvp[1, 6] <- coeffRV(phyto1, zoop2)$p.value rvp[1, 7] <- coeffRV(phyto1, zoop3)$p.value rvp[1, 8] <- coeffRV(phyto1, zoop4)$p.value rvp[1, 9] <- coeffRV(phyto1, fish1)$p.value rvp[1, 10] <- coeffRV(phyto1, fish2)$p.value rvp[1, 11] <- coeffRV(phyto1, fish3)$p.value rvp[2, 3] <- coeffRV(phyto2, phyto3)$p.value rvp[2, 4] <- coeffRV(phyto2, phyto4)$p.value rvp[2, 5] <- coeffRV(phyto2, zoop1)$p.value rvp[2, 6] <- coeffRV(phyto2, zoop2)$p.value rvp[2, 7] <- coeffRV(phyto2, zoop3)$p.value rvp[2, 8] <- coeffRV(phyto2, zoop4)$p.value rvp[2, 9] <- coeffRV(phyto2, fish1)$p.value rvp[2, 10] <- coeffRV(phyto2, fish2)$p.value rvp[2, 11] <- coeffRV(phyto2, fish3)$p.value rvp[3, 4] <- coeffRV(phyto3, phyto4)$p.value rvp[3, 5] <- coeffRV(phyto3, zoop1)$p.value rvp[3, 6] <- coeffRV(phyto3, zoop2)$p.value rvp[3, 7] <- coeffRV(phyto3, zoop3)$p.value rvp[3, 8] <- coeffRV(phyto3, zoop4)$p.value rvp[3, 9] <- coeffRV(phyto3, fish1)$p.value rvp[3, 10] <- coeffRV(phyto3, fish2)$p.value rvp[3, 11] <- coeffRV(phyto3, fish3)$p.value rvp[4, 5] <- coeffRV(phyto4, zoop1)$p.value rvp[4, 6] <- coeffRV(phyto4, zoop2)$p.value rvp[4, 7] <- coeffRV(phyto4, zoop3)$p.value rvp[4, 8] <- coeffRV(phyto4, zoop4)$p.value rvp[4, 9] <- coeffRV(phyto4, fish1)$p.value rvp[4, 10] <- coeffRV(phyto4, fish2)$p.value rvp[4, 11] <- coeffRV(phyto4, fish3)$p.value rvp[5, 6] <- coeffRV(zoop1, zoop2)$p.value rvp[5, 7] <- coeffRV(zoop1, zoop3)$p.value rvp[5, 8] <- coeffRV(zoop1, zoop4)$p.value rvp[5, 9] <- coeffRV(zoop1, fish1)$p.value rvp[5, 10] <- coeffRV(zoop1, fish2)$p.value rvp[5, 11] <- coeffRV(zoop1, fish3)$p.value rvp[6, 7] <- coeffRV(zoop2, zoop3)$p.value rvp[6, 8] <- coeffRV(zoop2, zoop4)$p.value rvp[6, 9] <- coeffRV(zoop2, fish1)$p.value rvp[6, 10] <- coeffRV(zoop2, fish2)$p.value rvp[6, 11] <- coeffRV(zoop2, fish3)$p.value rvp[7, 8] <- coeffRV(zoop3, zoop4)$p.value rvp[7, 9] <- coeffRV(zoop3, fish1)$p.value rvp[7, 10] <- coeffRV(zoop3, fish2)$p.value rvp[7, 11] <- coeffRV(zoop3, fish3)$p.value rvp[8, 9] <- coeffRV(zoop4, fish1)$p.value rvp[8, 10] <- coeffRV(zoop4, fish2)$p.value rvp[8, 11] <- coeffRV(zoop4, fish3)$p.value rvp[9, 10] <- coeffRV(fish1, fish2)$p.value rvp[9, 11] <- coeffRV(fish1, fish3)$p.value rvp[10, 11] <- coeffRV(fish2, fish3)$p.value round(rvp[-12,-12],5) # 2.1) MFA #### # Variable selection via PCA contibution circle # Phyto phyto.pca <- rda(phyto.chord) cleanplot.pca(phyto.pca,ax1=1,ax2=2) cleanplot.pca(phyto.pca,ax1=2,ax2=3) cleanplot.pca(phyto.pca,ax1=3,ax2=4) cleanplot.pca(phyto.pca,ax1=4,ax2=5) # Zoop zoop.pca <- rda(zoop.chord) cleanplot.pca(zoop.pca,ax1=1,ax2=2) cleanplot.pca(zoop.pca,ax1=2,ax2=3) cleanplot.pca(zoop.pca,ax1=3,ax2=4) cleanplot.pca(zoop.pca,ax1=4,ax2=5) # Fish fish.pca <- rda(fish.chord) cleanplot.pca(fish.pca,ax1=1,ax2=2) cleanplot.pca(fish.pca,ax1=2,ax2=3) cleanplot.pca(fish.pca,ax1=3,ax2=4) cleanplot.pca(fish.pca,ax1=4,ax2=5) # Env env <- cbind(physico,morpho) # re-runto bring back variables with NAs that were removed above env <- scale(env); env <- as.data.frame(env) env <- env[,-c(18)] # Cond appears twice env.pca <- rda(na.omit(env)) cleanplot.pca(env.pca,ax1=1,ax2=2) cleanplot.pca(env.pca,ax1=2,ax2=3) cleanplot.pca(env.pca,ax1=3,ax2=4) cleanplot.pca(env.pca,ax1=4,ax2=5) # Running MFA with dominant vatiables (selected using cleanplot.pca() above) phyto.esp <- phyto.chord[,c("MALL","NAVI","MEMI","SEMI","TEMI","METE","BICH","OOCY","APNI","ASFO","SYNA")] zoop.esp <- zoop.chord[,c("COCO","KELO","KETA","POVU","TRCY","HOGI","DAGA","LEMI","MEED")] fish.esp <- fish.chord[,c("SAFO","CACO","ESLU","STIVI","PEFL")] env.dis <- env[,c("H","Ca","Alc","Tr","TOC","DOC","Al","Zmax","Long","Zrel","MEI","Vol")] tab1 <- cbind(phyto.esp,zoop.esp,fish.esp,env.dis) grn <- c(ncol(phyto.esp),ncol(zoop.esp),ncol(fish.esp),ncol(env.dis)) QC.mfa <- MFA(tab1, group=grn, type=c("c","c","c","s"), ncp=5, name.group=c("Phyto","Zoop","Fish","Env")) summary(QC.mfa) plot(QC.mfa,choix = "var",habillage = "group", shadowtext = FALSE) # Alternative MFA plots # Plot only the significant variables (correlations) aa <- dimdesc(QC.mfa, axes = 1:2, proba = 0.05) varsig <- QC.mfa$quanti.var$cor[unique(c(rownames(aa$Dim.1$quanti), rownames(aa$Dim.2$quanti))), ] plot(varsig[, 1:2],asp = 1, type = "n", xlim = c(-1, 1), ylim = c(-1, 1)) abline(h = 0, lty = 3) abline(v = 0, lty = 3) symbols(0, 0, circles = 1, inches = FALSE, add = TRUE) arrows(0, 0, varsig[, 1], varsig[, 2], length = 0.08, angle = 20) for (v in 1 : nrow(varsig)) { if (abs(varsig[v, 1]) > abs(varsig[v, 2])) { if (varsig[v, 1] >= 0) pos <- 4 else pos <- 2 } else { if (varsig[v, 2] >= 0) pos <- 3 else pos <- 1 } text(varsig[v, 1], varsig[v, 2], labels = rownames(varsig)[v], pos = pos) } # Plot of partial axes plot(QC.mfa, choix="axes") fviz_mfa_axes(QC.mfa) # First dimension of env, fish, and zoop are highly correlated # The second dimension of the MFA is correlated to the first dimension of phyto, and the second dimensions of zoop, env. # Numeric variables (correlation circle) plot(QC.mfa, choix="var") # Sites (global and partial viewpoints) plot(QC.mfa, choix="ind", habillage="group", partial="all", invisible="quali") # Proportion of variance explained eig.val <- get_eigenvalue(QC.mfa) fviz_screeplot(QC.mfa) # The proportion of variances retained by the different dimensions (axes) # Eigenvalues, scree plot and broken stick model ev <- QC.mfa$eig[, 1] names(ev) <- paste("MFA", 1 : length(ev)) screestick(ev, las = 2) # Groups of variables group <- get_mfa_var(QC.mfa, "group") fviz_mfa_var(QC.mfa, "group") # red colour - active groups, green colour - suppl groups. # The plot illustrates the correlation between groups and dimensions. # The coordinates of the four active groups vary along the first dimension. # Env and Fish contribute the most to Dim 1, followed by Zoop and then Phyto. # Phyto has the highest coordinates on Dim 2 indicating a highest contribution to # this dimension, whereas the other three groups contribute about the same. # Bar plot of groups contribution to the dimensions fviz_contrib(QC.mfa, "group", axes = 1) # Contribution to the first dimension fviz_contrib(QC.mfa, "group", axes = 2) # Contribution to the second dimension fviz_contrib(QC.mfa, "group", axes = 3) # Contributions to the third-fifth dimensions... fviz_contrib(QC.mfa, "group", axes = 4) fviz_contrib(QC.mfa, "group", axes = 5) # Quantitative variables quanti.var <- get_mfa_var(QC.mfa, "quanti.var") quanti.var fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE) # Dims 1 and 2 fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE,axes=c(2,3)) # And other axes fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE,axes=c(3,4)) fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE,axes=c(4,5)) # Quantitative variables colored by groups plot(QC.mfa,choix="var",habillage="group",cex=0.8,select="contrib 24",axes=c(1,2)) # Plotting without fviz (so that arrows are equal) plot(QC.mfa,choix = "var",habillage = "group", shadowtext = FALSE,axes=c(2,3)) plot(QC.mfa,choix = "var",habillage = "group", shadowtext = FALSE,axes=c(3,4)) plot(QC.mfa,choix = "var",habillage = "group", shadowtext = FALSE,axes=c(4,5)) # same plot but with points instead of arrows fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE, geom = c("point", "text"), legend = "bottom") # Contributions of individual vars fviz_contrib(QC.mfa, choice = "quanti.var", axes = 1, top = 20,palette = "jco") # Contributions to dimension 1 # The red dashed line on the graph above indicates the expected average value, if contributions were uniform. fviz_contrib(QC.mfa, choice = "quanti.var", axes = 2, top = 20, palette = "jco") # Contributions to dimension 2 fviz_contrib(QC.mfa, choice = "quanti.var", axes = 3, top = 20, palette = "jco") # Contributions to dimension 3 fviz_contrib(QC.mfa, choice = "quanti.var", axes = 4, top = 20, palette = "jco") # Contributions to dimension 4 fviz_contrib(QC.mfa, choice = "quanti.var", axes = 5, top = 20, palette = "jco") # Contributions to dimension 5 # Most contributing quantitative variables can be highlighted on scatter plot and gradient of color fviz_mfa_var(QC.mfa, "quanti.var", col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), col.var.sup = "violet", repel = TRUE, geom = c("point", "text")) fviz_mfa_var(QC.mfa, "quanti.var", col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), col.var.sup = "violet", repel = TRUE) # Graph of sites fviz_mfa_ind(QC.mfa, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE) # Graph of sites for other dimensions fviz_mfa_ind(QC.mfa, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE,axes=c(2,3)) fviz_mfa_ind(QC.mfa, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE,axes=c(3,4)) fviz_mfa_ind(QC.mfa, col.ind = "cos2",gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE,axes=c(4,5)) fviz_mfa_ind(QC.mfa, partial = c("231B", "20B", "230","228","229","215B","216","164")) # 2.2) LVM #### # PCA of env variables to extract dominant PCs, which will be used as expl vars in LVMs pca.env <- rda(na.exclude(data_m[,c(123:154,156:162)]),scale=TRUE) # row 7 and 12 have NAs for MEI and DevVol screeplot(pca.env,bstick = TRUE,npcs = length(pca.env$CA$eig)) plot(pca.env,choices=c(1,2)) cleanplot.pca(pca.env) res.PCA.env <- PCA(na.exclude(data_m[,c(123:154,156:162)]),scale.unit=TRUE,ncp=5,graph=T) # Plots by axes PCA_12 <- fviz_pca_var(res.PCA.env, axes=c(1,2),repel=TRUE, col.var = "contrib", # gradient.cols = c("white", "#f47a27", "#ed4040"),title = "PCA - Environmental variables", gradient.cols = c("gray90", "gray70", "gray20"),title = "PCA - Environmental variables", ggtheme = theme_minimal()) + xlim(c(1,-1)) + ylim(c(1,-1)) PCA_23 <- fviz_pca_var(res.PCA.env, axes=c(2,3), repel=TRUE, col.var = "contrib", # gradient.cols = c("white", "#ffd700","#c7c216"),title = "PCA - Environmental variables", gradient.cols = c("gray90", "gray70", "gray20"),title = "PCA - Environmental variables", ggtheme = theme_minimal()) + xlim(c(1,-1)) + ylim(c(1,-1)) PCA_34 <- fviz_pca_var(res.PCA.env, axes=c(3,4), repel=TRUE, col.var = "contrib", # gradient.cols = c("white", "#ffd700","#c7c216"),title = "PCA - Environmental variables", gradient.cols = c("gray90", "gray70", "gray20"),title = "PCA - Environmental variables", ggtheme = theme_minimal()) + xlim(c(1,-1)) + ylim(c(1,-1)) PCA_45 <- fviz_pca_var(res.PCA.env, axes=c(4,5), repel=TRUE, col.var = "contrib", # gradient.cols = c("white","#b5066a","#721f81"),title = "PCA - Environmental variables", gradient.cols = c("gray90", "gray70", "gray20"),title = "PCA - Environmental variables", ggtheme = theme_bw()) + xlim(c(1,-1)) + ylim(c(1,-1)) grid.arrange(PCA_12,PCA_23,PCA_34,PCA_45,ncol=2) Loading1 <- scores(pca.env,choices=1:7)$species[,1] # PCA1 = Cond, Ca, Alc, CIT, pH, DIC, K, SO4, Mg Loading2 <- scores(pca.env,choices=1:7)$species[,2] # PCA2 = Max depth, Mean depth, Vol, Rel depth, Tr Loading3 <- scores(pca.env,choices=1:7)$species[,3] # PCA3 = Order, TOC, DOC, CA, Colour, Cl, Al, Loading4 <- scores(pca.env,choices=1:7)$species[,4] # PCA4 = SA, Wmax, Lmax, Mg, Lat Loading5 <- scores(pca.env,choices=1:7)$species[,5] # PCA5 = Lake dev volume, Altitude Loading6 <- scores(pca.env,choices=1:7)$species[,6] # PCA6 = Na, K, H Loading7 <- scores(pca.env,choices=1:7)$species[,7] # PCA7 = CA, NO3, Order, Cl, Cu sort(Loading1) # and so on... axes1 <- scores(pca.env,choices=1:7)$sites[,1] axes2 <- scores(pca.env,choices=1:7)$sites[,2] axes3 <- scores(pca.env,choices=1:7)$sites[,3] axes4 <- scores(pca.env,choices=1:7)$sites[,4] axes5 <- scores(pca.env,choices=1:7)$sites[,5] axes6 <- scores(pca.env,choices=1:7)$sites[,6] axes7 <- scores(pca.env,choices=1:7)$sites[,7] # Remove lines with missing values for MEI data_mNAr <- data_m[-7,] data_mNAr <- data_mNAr[-11,] # Add PCA axes to merged data data_mNAr <- cbind(data_mNAr,scores(pca.env,choices=1:7)$sites) write.csv(data_mNAr,"data_mNAr") data_mNAr <- read.csv("data_mNAr"); data_mNAr <- data_mNAr[,-1] # Re-separate sp and env (PCs) data phyto <- data_mNAr[,36:122] zoop <- data_mNAr[,2:35] fish <- data_mNAr[,172:189] env <- data_mNAr[,c(190:196)] # select latent variables from PCA # Transform species data (used Hellinger, chord, or log-chord) phyto.hell <- decostand(phyto,"hellinger") zoop.hell <- decostand(zoop,"hellinger") fish.hell <- decostand(fish,"hellinger") names(phyto.hell) <- paste(c(rep("phyto",length(names(phyto.hell)))),names(phyto.hell),sep="_") names(zoop.hell) <- paste(c(rep("zoop",length(names(zoop.hell)))),names(zoop.hell),sep="_") names(fish.hell) <- paste(c(rep("fish",length(names(fish.hell)))),names(fish.hell),sep="_") # Select phyto, zoop or fish data abund <- phyto.hell # abund <- zoop.hell # abund <- fish.hell # abund <- cbind(phyto.hell,zoop.hell,fish.hell) # Use PCs as predictors (selecting one at a time for fit.lvm.PCX below) covX <- cbind(env[,"PC1"],env[,"PC1"]^2) # covX <- cbind(env[,"PC2"],env[,"PC2"]^2) # covX <- cbind(env[,"PC3"],env[,"PC3"]^2) # covX <- cbind(env[,"PC4"],env[,"PC4"]^2) # covX <- cbind(env[,"PC5"],env[,"PC5"]^2) # covX <- cbind(env[,"PC6"],env[,"PC6"]^2) # covX <- cbind(env[,"PC7"],env[,"PC7"]^2) covAll <- cbind(env[,"PC1"],env[,"PC1"]^2,env[,"PC2"],env[,"PC2"]^2, env[,"PC3"],env[,"PC3"]^2,env[,"PC4"],env[,"PC4"]^2, env[,"PC5"],env[,"PC5"]^2,env[,"PC6"],env[,"PC6"]^2, env[,"PC7"],env[,"PC7"]^2) fit.lvm.PCs <- boral(y = abund, X = covAll, lv.control = list(num.lv = 2), family = "normal", save.model = TRUE, prior.control = list(type = c("cauchy","cauchy","cauchy","halfcauchy"), hypparams = c(2.5^2,2.5^2,2.5^2,2.5^2))) fit.lvm.PCX <- boral(y = abund, X = covX, lv.control = list(num.lv = 2), family = "normal", save.model = TRUE, prior.control = list(type = c("cauchy","cauchy","cauchy","halfcauchy"), hypparams = c(2.5^2,2.5^2,2.5^2,2.5^2))) # Testing for different number of latent variables # fit.lvm.PCs.lv2 <- boral(y = abund, X = covAll, lv.control = list(num.lv = 2), family = "normal", save.model = TRUE, prior.control = list(type = c("cauchy","cauchy","cauchy","halfcauchy"), hypparams = c(2.5^2,2.5^2,2.5^2,2.5^2))) # fit.lvm.PCs.lv3 <- boral(y = abund, X = covAll, lv.control = list(num.lv = 3), family = "normal", save.model = TRUE, prior.control = list(type = c("cauchy","cauchy","cauchy","halfcauchy"), hypparams = c(2.5^2,2.5^2,2.5^2,2.5^2))) # fit.lvm.PCs.lv4 <- boral(y = abund, X = covAll, lv.control = list(num.lv = 4), family = "normal", save.model = TRUE, prior.control = list(type = c("cauchy","cauchy","cauchy","halfcauchy"), hypparams = c(2.5^2,2.5^2,2.5^2,2.5^2))) # fit.lvm.PCs.lv5 <- boral(y = abund, X = covAll, lv.control = list(num.lv = 5), family = "normal", save.model = TRUE, prior.control = list(type = c("cauchy","cauchy","cauchy","halfcauchy"), hypparams = c(2.5^2,2.5^2,2.5^2,2.5^2))) fit.lvm <- fit.lvm.PCX summary(fit.lvm)$ics # Estimated parameter values fit.lvm$hpdintervals # 95% credible intervals for model parameters # Examine coefficients: cbind(fit.lvm$lv.coefs.mean,fit.lvm$X.coefs.mean) ## In boral, the random intercepts are sampled from a non-zero normal distribution # Dunn-Smyth residual plots to check model assumption, outliers etc. par(mfrow=c(2,2)) plot(fit.lvm) # Residual ordination plot of sites lvsplot(fit.lvm) env.cors <- get.enviro.cor(fit.lvm) # the correlation matrix due to similarities in the response corrplot(env.cors$cor, diag = F, type = "lower", title = "Environmental correlations from LVM",mar=c(0,0,1,0),order = "FPC") res.cors <- get.residual.cor(fit.lvm) corrplot(res.cors$cor, diag = F, type = "lower", title = "Residual correlations from LVM", tl.srt = 45,order = "FPC") # Extract correlations upperTriangle <- upper.tri(env.cors$cor, diag=F) # turn into a upper triangle correlations.upperTriangle <- env.cors$cor # take a copy of the original cor-mat correlations.upperTriangle[!upperTriangle] <- NA # set everything not in upper triangle o NA correlations_melted <- na.omit(melt(correlations.upperTriangle, value.name ="correlationCoef")) # use melt to reshape the matrix into triplets, na.omit to get rid of the NA rows colnames(correlations_melted) <- c("X1", "X2", "correlation") # Extract significance upperTriangle <- upper.tri(env.cors$sig.cor, diag=F) correlations.upperTriangle <- env.cors$sig.cor # take a copy of the original cor-mat correlations.upperTriangle[!upperTriangle] <- NA # set everything not in upper triangle o NA correlations_melted_p <- na.omit(melt(correlations.upperTriangle, value.name ="correlationSign")) # use melt to reshape the matrix into triplets, na.omit to get rid of the NA rows colnames(correlations_melted_p) <- c("X1", "X2", "p.value"); correlations_melted_p <- correlations_melted_p[,3] correlations_melted <- cbind(correlations_melted,correlations_melted_p) # Merge datasets corr_PC1 <- correlations_melted df_PC1 <- mutate(corr_PC1,col1=as.character(X1)) df_PC1 <- mutate(df_PC1,col1=sapply(strsplit(df_PC1$col1, split='_', fixed=TRUE),function(x) (x[1]))) df_PC1 <- mutate(df_PC1,col2=as.character(X2)) df_PC1 <- mutate(df_PC1,col2=sapply(strsplit(df_PC1$col2, split='_', fixed=TRUE),function(x) (x[1]))) df_PC1$grp <- paste(df_PC1$col1,df_PC1$col2,sep="_") df_PC1_neg <- df_PC1[df_PC1$correlation < 0,] df_PC1_pos <- df_PC1[df_PC1$correlation >= 0,] names(df_PC1)[4] <- c("p.value") df_PC1 <- df_PC1[,c(1:4,7)] df_PC1$var <- rep("PC1",dim(df_PC1)[1]) df_PC1$coresp <- df_PC1$correlation df_PC1[df_PC1$coresp > 0,]$coresp <- 2 df_PC1[df_PC1$coresp < 0,]$coresp <- 1 df_PC1$coresp <- as.factor(df_PC1$coresp) df_PC1$var <- as.factor(df_PC1$var) # repeat above (lines 602 onwards) changing PC vatiable from PC1 to PC5 # and saving lines 667-684 as df_PC2, df_PC3, df_PC4, df_PC5, respectively #... #... # combine all 5 axes results df_PCs <- rbind(df_PC1,df_PC2,df_PC3,df_PC4,df_PC5) Fig_all <- ggplot(df_PCs,aes(var)) + scale_colour_brewer(palette = "Spectral") + scale_fill_brewer(palette = "Spectral") + ylab('Frequency of bivariate correlations \n among LVM fitted responses') + xlab('Distance') + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + geom_bar(aes(fill=coresp)) + facet_wrap(~grp,ncol=3) Fig_all # Examine significant coresponses only df_PCs.sign <- df_PCs[df_PCs$p.value != 0,] Fig_coresp <- ggplot(df_PCs.sign,aes(var)) + scale_colour_brewer(palette = "Spectral") + scale_fill_viridis(discrete=TRUE,option="magma") + ylab('Frequency of signif bivariate correlations \n among LVM fitted responses') + xlab('Distance') + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + geom_bar(aes(fill=coresp)) + facet_wrap(~grp,ncol=3) Fig_coresp Fig_freq_stack <- ggplot(df_PCs.sign,aes(var)) + scale_fill_viridis(discrete=TRUE,option="magma") + scale_colour_viridis(discrete=TRUE,option="magma") + ylab('Frequency of signif bivariate correlations \n among LVM fitted responses') + xlab('Distance') + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + geom_bar(aes(fill=grp)) Fig_freq_stack Fig_freq <- ggplot(df_PCs.sign,aes(grp)) + scale_fill_viridis(discrete=TRUE,option="magma") + scale_colour_viridis(discrete=TRUE,option="magma") + ylab('Frequency of signif bivariate correlations \n among LVM fitted responses') + xlab('Distance') + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + geom_bar(aes(fill=var),stat="count", position="dodge") Fig_freq # Make into proportion df_PCs.sign$counter <- c(rep(1,dim(df_PCs.sign)[1])) df_PCs$counter <- c(rep(1,dim(df_PCs)[1])) df.agg <- aggregate(counter ~ grp + var, data=df_PCs.sign, sum) df.agg.ns <- aggregate(counter ~ grp + var, data=df_PCs, sum) df.agg$prop <- df.agg$counter/df.agg.ns$counter Fig_prop <- ggplot(df.agg,aes(grp)) + scale_fill_viridis(discrete=TRUE,option="magma") + scale_colour_viridis(discrete=TRUE,option="magma") + ylab('Proportion of signif bivariate correlations \n among LVM fitted responses') + xlab('') + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + geom_bar(aes(y=prop,fill=var),stat="identity",position="dodge") Fig_prop # co-response plot df.agg <- aggregate(counter ~ grp + var + coresp, data=df_PCs.sign, sum) df.agg$ID <- paste(df.agg$grp,df.agg$var,df.agg$coresp,sep="_") df.agg.ns <- aggregate(counter ~ grp + var + coresp, data=df_PCs, sum) df.agg.ns$ID <- paste(df.agg.ns$grp,df.agg.ns$var,df.agg.ns$coresp,sep="_") df.merged <- merge(df.agg,df.agg.ns,by="ID") df.merged$prop <- df.merged$counter.x/df.merged$counter.y Fig_prop_coresp <- ggplot(df.merged,aes(var.x)) + scale_colour_brewer(palette = "Spectral") + scale_fill_grey() + ylab('Proportion of signif bivariate correlations \n among LVM fitted responses') + xlab('') + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + geom_bar(aes(y=prop,fill=coresp.x),stat="identity") + facet_wrap(~grp.x,ncol=3) Fig_prop_coresp # Saving data (*NOTE* files already saved, no need to run -> can simply load files in section 2.2.2) # phyto_res <- list(d1=PC1, d2=PC2, d3=PC3, d4=PC4, d5=PC5) # create a named list of your data frames # lapply(phyto_res, "[", 1) # names(phyto_res) <- c("PC1","PC2","PC3","PC4","PC5") # save(phyto_res, file = "allcorrs.phyto.corrected.RData") # # Repeat above section (starting on line 602) for zoop... then save # zoop_res <- list(d1=PC1, d2=PC2, d3=PC3, d4=PC4, d5=PC5) ## create a named list of your data frames # lapply(zoop_res, "[", 1) # names(zoop_res) <- c("PC1","PC2","PC3","PC4","PC5") # save(zoop_res, file = "allcorrs.zoop.corrected.RData") # # Repeat above section (starting on line 602) for fish then save # fish_res <- list(d1=PC1, d2=PC2, d3=PC3, d4=PC4, d5=PC5) ## create a named list of your data frames # lapply(fish_res, "[", 1) # names(fish_res) <- c("PC1","PC2","PC3","PC4","PC5") # save(fish_res, file = "allcorrs.fish.corrected.RData") # 2.2.1) ____ - Map of env PCs #### # plot limits xlim = c(-85,-52) ylim = c(42,63) worldmap = map_data("world","canada") setnames(worldmap, c("X","Y","PID","POS","region","subregion")) worldmap = clipPolys(worldmap, xlim=xlim,ylim=ylim, keepExtra=TRUE) colorRampPalette(brewer.pal(11,"Spectral"))(11) p_PC1 = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "gray96",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("Hardness-alkanity") + # geom_point(aes(fill=data_mNAr$PC1),colour="white",size=data_mNAr$PC1*3,shape=21) + scale_fill_gradient(low="white",high="#ed4040",name="PC1") # old col: #D53E4F geom_point(aes(fill=data_mNAr$PC1),colour="white",size=data_mNAr$PC1*4,shape=21) + scale_fill_gradient2(low="gray98",mid="gray80",high="gray25",name="PC1") # old col: #D53E4F p_PC1 <- p_PC1 + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_PC1 p_PC2 = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "gray96",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("Depth-volume-transparency") + # geom_point(aes(fill=data_mNAr$PC2),colour="white",size=data_mNAr$PC2*3,shape=21) + scale_fill_gradient(low="white",high="#f47a27",name="PC2") # old col: #FC8D59 geom_point(aes(fill=data_mNAr$PC2),colour="white",size=data_mNAr$PC2*3,shape=21) + scale_fill_gradient2(low="gray98",mid="gray80",high="gray25",name="PC2") # old col: #FC8D59 p_PC2 <- p_PC2 + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_PC2 p_PC3 = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "gray96",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("Dystrophy-colour") + # geom_point(aes(fill=data_mNAr$PC3),colour="white",size=data_mNAr$PC3*3,shape=21) + scale_fill_gradient(low="white",high="#ddd707",name="PC3") # old col: #FEE08B geom_point(aes(fill=data_mNAr$PC3),colour="white",size=data_mNAr$PC3*3,shape=21) + scale_fill_gradient2(low="gray98",mid="gray80",high="gray25",name="PC3") # old col: #FEE08B p_PC3 <- p_PC3 + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_PC3 p_PC4 = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "gray96",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("Lake size") + # geom_point(aes(fill=data_mNAr$PC4),colour="white",size=data_mNAr$PC4*3,shape=21) + scale_fill_gradient(low="white",high="#b5066a",name="PC4") # old col: #ABDDA4 geom_point(aes(fill=data_mNAr$PC4),colour="white",size=data_mNAr$PC4*3,shape=21) + scale_fill_gradient2(low="gray98",mid="gray80",high="gray25",name="PC4") p_PC4 <- p_PC4 + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_PC4 p_PC5 = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "gray96",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("Lake volume dev") + # geom_point(aes(fill=data_mNAr$PC5),colour="white",size=data_mNAr$PC5*3,shape=21) + scale_fill_gradient(low="white",high="#721f81",name="PC5") # old col: #66C2A5 geom_point(aes(fill=data_mNAr$PC5),colour="white",size=data_mNAr$PC5*3,shape=21) + scale_fill_gradient2(low="gray98",mid="gray80",high="gray25",name="PC5") p_PC5 <- p_PC5 + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_PC5 quartz() # if on Mac, or windows() if on Windows PC grid.arrange(p_PC1,p_PC2,p_PC3,p_PC4,p_PC5,ncol=3) # 2.2.2) ____ - network figures #### load(file <- "allcorrs.phyto.corrected.RData") allcorrs.names <- paste("phyto_res$", names(phyto_res), "$sig.cor", sep = "") # File for zoop # load(file <- "allcorrs.zoop.corrected.RData") # allcorrs.names <- paste("zoop_res$", names(zoop_res), "$sig.cor", sep = "") # File for fish # load(file <- "allcorrs.fish.corrected.RData") # allcorrs.names <- paste("fish_res$", names(fish_res), "$sig.cor", sep = "") # Circle plot: niche separation num.negs <- NULL # how many negative correlations for each factor? for (i in 1:length(allcorrs.names)){ num.negs[i] <- sum(c(eval(parse(text = allcorrs.names[i]))) < 0) } num.negs # order names by number of neg correlations for PCs allcorrs.shuff <- allcorrs.names[c(1,4,3,5,2)] # Phyto # allcorrs.shuff <- allcorrs.names[c(1,4,3,5,2)] # Zoop # allcorrs.shuff <- allcorrs.names[c(1,2,4,3,5)] # Fish unlist(strsplit(allcorrs.shuff[1], split = "$", fixed = TRUE))[2] # Names for PCs title_names <- c("Hardness-alkalinity (PC1)", "Lake SA (PC4)", "Dystrophy/ Lake order/ CA (PC3)", "lake devel/ altitude (PC5)","depth-vol-tr (PC2)") # Phyto # title_names <- c("Hardness-alkalinity (PC1)", "Lake SA (PC4)", "Dystrophy/ Lake order/ CA (PC3)", # "lake devel/ altitude (PC5)","depth-vol-tr (PC2)") # Zoop # title_names <- c("hardness-alkalinity", "depth-vol-tr (PC2)","Lake SA (PC4)", # "Dystrophy/ Lake order/ CA (PC3)","lake devel/ altitude (PC5)") # Fish dummy <- as.dist(phyto_res$PC4$sig.cor) # Phyto # dummy <- as.dist(zoop_res$PC1$sig.cor) # Zoop # dummy <- as.dist(fish_res$PC1$sig.cor) # Fish # Running the fig dummy.points <- (point.attr(dummy)) dummy.points$size <- 5 deep.circle <- my.circleplot(dummy, cluster = FALSE, plot.control=list(points = dummy.points, line.breaks = c(0:1, by = 0.01), line.curvature = 0)) par(mfrow = c(3,3), mar = c(0,0,0,0)) for (i in 1:length(allcorrs.shuff)){ # isolate negative connects only mydist <- eval(parse(text = allcorrs.shuff[i])) bluh <- mydist diag(bluh) <- 0 bluh[bluh > 0] <- 0 # set positive correlations to zero bluh[bluh < 0] <- 1 mydist <- mydist[rowSums(bluh)>0,] mydist <- mydist[,colSums(bluh)>0] mydist <- mydist[order(mydist[c(3,4,1,2,3),1]),] bluh <- bluh[rowSums(bluh)>0,] bluh <- bluh[,colSums(bluh)>0] deep.graph <- graph.adjacency(bluh, weighted = TRUE, mode = "undirected", diag = FALSE) # plot(deep.graph) circ.deep <- as.dist(mydist) circ.deep[circ.deep > 0] <- 0 my.points <- (point.attr(circ.deep)) my.points$size <- 5 mycol <- rep("grey",length(V(deep.graph))) mycol[largest.independent.vertex.sets(deep.graph)[[1]]] <- "white" my.points$colour <- mycol colfunc <- colorRampPalette(c("grey", "Black")) lwdscale <- sort(c(abs(circ.deep)))[sort(c(abs(circ.deep))) != 0] if(length(lwdscale) > 0){ my.circleplot(abs(circ.deep), cluster = FALSE, plot.control=list(points = my.points, line.breaks = c(seq(min(lwdscale)-0.01,max(lwdscale)+0.01, by = 0.01)), line.curvature = 0, line.cols = colfunc(length(seq(min(lwdscale)-0.01,max(lwdscale)+0.01, by = 0.01))-1), line.width = c(0.01,3))) title(sub = title_names[i], line = -1.2, cex.sub = 1.5)} else { plot(x = deep.circle$points$x, y = deep.circle$points$y, type = "p", ann = FALSE, axes = FALSE, asp = 1, pch = 21, bg = my.points$colour, cex = my.points$size, xlim = c(min(deep.circle$points$x)-0.5, max(deep.circle$points$x)+0.5)) text(deep.circle$points$x, deep.circle$points$y, label = deep.circle$points$label, col = "black", cex = 1, font = 4) title(sub = title_names[i], line = -1.2, cex.sub = 1.5)} } # Circle plot: co-response num.pos <- NULL # how many positive correlations for each factor? for (i in 1:length(allcorrs.names)){ num.pos[i] <- sum(c(eval(parse(text = allcorrs.names[i]))) > 0) } num.pos # order names by number of neg correlations allcorrs.shuff <- allcorrs.names[c(1,2,4,5,3)] # Phyto # allcorrs.shuff <- allcorrs.names[c(1,3,5,4,2)] # Zoop # allcorrs.shuff <- allcorrs.names[c(3,2,4,1,5)] # Fish unlist(strsplit(allcorrs.shuff[1], split = "$", fixed = TRUE))[2] # Names for PCs title_names <- c("hardness-alkalinity (PC1)", "depth-vol-tr (PC2)","lake SA (PC4)", "lake vol dev - altitude (PC5)","dystrophy-colour (PC3)") # Phyto # title_names <- c("hardness-alkalinity (PC1)", "dystrophy-colour (PC3)","lake vol dev - altitude (PC5)", # "lake SA (PC4)","depth-vol-tr (PC2)") # Zoop # title_names <- c("dystrophy-colour (PC3)","depth-vol-tr (PC2)","lake SA (PC4)", # "hardness-alkalinity (PC1)","lake vol dev - altitude (PC5)") # Fish dummy <- as.dist(phyto_res$PC1$sig.cor) # Phyto # dummy <- as.dist(zoop_res$PC1$sig.cor) # Zoop # dummy <- as.dist(fish_res$PC3$sig.cor) # Fish # Circle diagram dummy.points <- (point.attr(dummy)) dummy.points$size <- 5 deep.circle <- my.circleplot(dummy, cluster = FALSE, plot.control=list(points = dummy.points, line.breaks = c(0:1, by = 0.01), line.curvature = 0)) par(mfrow = c(3,3), mar = c(0,0,0,0)) for (i in 1:length(allcorrs.shuff)){ # isolate negative connects only mydist <- eval(parse(text = allcorrs.shuff[i])) bluh <- mydist diag(bluh) <- 0 bluh[bluh > 0] <- 1 bluh[bluh < 0] <- 0 # set negative correlations to zero mydist <- mydist[rowSums(bluh)>0,] mydist <- mydist[,colSums(bluh)>0] bluh <- bluh[rowSums(bluh)>0,] bluh <- bluh[,colSums(bluh)>0] deep.graph <- graph.adjacency(bluh, weighted = TRUE, mode = "undirected", diag = FALSE) # plot(deep.graph) circ.deep <- as.dist(mydist) circ.deep[circ.deep < 0] <- 0 # if set negative correlations to zero my.points <- (point.attr(circ.deep)) my.points$size <- 5 mycol <- rep("grey",length(V(deep.graph))) mycol[largest.independent.vertex.sets(deep.graph)[[1]]] <- "white" my.points$colour <- mycol colfunc <- colorRampPalette(c("grey", "Black")) lwdscale <- sort(c(abs(circ.deep)))[sort(c(abs(circ.deep))) != 0] if(length(lwdscale) > 0){ my.circleplot(abs(circ.deep), cluster = FALSE, plot.control=list(points = my.points, line.breaks = c(seq(min(lwdscale)-0.01,max(lwdscale)+0.01, by = 0.01)), line.curvature = 0, line.cols = colfunc(length(seq(min(lwdscale)-0.01,max(lwdscale)+0.01, by = 0.01))-1), line.width = c(0.01,3))) title(sub = title_names[i], line = -1.2, cex.sub = 1.5)} else { plot(x = deep.circle$points$x, y = deep.circle$points$y, type = "p", ann = FALSE, axes = FALSE, asp = 1, pch = 21, bg = my.points$colour, cex = my.points$size, xlim = c(min(deep.circle$points$x)-0.5, max(deep.circle$points$x)+0.5)) text(deep.circle$points$x, deep.circle$points$y, label = deep.circle$points$label, col = "black", cex = 0.5, font = 4) title(sub = "", line = -1.2, cex.sub = 1.5)} } # 2.2.3) ____ - nonlinear response figures #### library(ggthemes) library(scales) phyto.hell <- decostand(phyto,"hellinger") zoop.hell <- decostand(zoop,"hellinger") fish.hell <- decostand(fish,"hellinger") # Phyto plots clrs.hcl <- function(n) { hcl(h = seq(100,140, length.out = n), c = 40, l = seq(10, 90, length.out = n), fixup = TRUE) } abund_phyto_select <- phyto.hell[,c("ARIN","CHLI","CRTE","APPEL","APDE","OSLI", "SYNA","APNI","SCDE","DIEH","METE","MEMI", "MEDI","QUAD","MIIN","GLOE","NAVI","PERI", "CHRO","DESM","CHSK","DIBA","BICH","COEL", "MALL","DINO")] frame1 <- data.frame(abund_phyto_select, env[,"PC1"]) gg1 <- gather(frame1, species, value, -env....PC1..) p1 <- ggplot(gg1, aes(y = value, x = env....PC1..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(26)) + xlab("PC1 (Alkalinity/ Acidity)") + ylab("Phytoplankton (hellinger)") p1 <- p1 + guides(colour = guide_legend(nrow = 9)) p1 abund_phyto_select <- phyto.hell[,c("MICO","ANFC","TAFE","BICH","BOSU", "KEPH","OOPU")] frame5 <- data.frame(abund_phyto_select, env[,"PC5"]) gg5 <- gather(frame5, species, value, -env....PC5..) p5 <- ggplot(gg5, aes(y = value, x = env....PC5..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(7)) + xlab("PC5 (Altitude/ Dev Vol)") + ylab("Phytoplankton (hellinger)") p5 <- p5 + guides(colour = guide_legend(nrow = 9)) p5 abund_phyto_select <- phyto.hell[,c("METE","PEIN","SCEN","NITZ","APNI","SYNA","MEMI")] frame3 <- data.frame(abund_phyto_select, env[,"PC3"]) gg3 <- gather(frame3, species, value, -env....PC3..) p3 <- ggplot(gg3, aes(y = value, x = env....PC3..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(7)) + xlab("PC3 (Dystrophy/ Order/ CA)") + ylab("Phytoplankton (hellinger)") p3 <- p3 + guides(colour = guide_legend(nrow = 9)) p3 abund_phyto_select <- phyto.hell[,c("OOCY","COEL","OSLI","CYLI","SYNA","GLOE","MIIN")] frame4 <- data.frame(abund_phyto_select, env[,"PC4"]) gg4 <- gather(frame4, species, value, -env....PC4..) p4 <- ggplot(gg4, aes(y = value, x = env....PC4..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(7)) + xlab("PC4 (Size/ Sulfate)") + ylab("Phytoplankton (hellinger)") p4 <- p4 + guides(colour = guide_legend(nrow = 9)) p4 grid.arrange(p1, p5, p3, p4, ncol= 4) # Zoop plots clrs.hcl <- function(n) { hcl(h = seq(200,240, length.out = n), c = 40, l = seq(10, 90, length.out = n), fixup = TRUE) } abund_zoop_select <- zoop.hell[,c("KETA","KEQU","FILO","DIAH","MEED","ASPR", "COCO","DAGA","DIBT","HOGI","BOLO","SICR", "LESI","DARO","CYSC")] frame1 <- data.frame(abund_zoop_select, env[,"PC1"]) gg1 <- gather(frame1, species, value, -env....PC1..) p1 <- ggplot(gg1, aes(y = value, x = env....PC1..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(15)) + xlab("PC1 (Alkalinity/ Acidity)") + ylab("Zooplankton (hellinger)") p1 <- p1 + guides(colour = guide_legend(nrow = 4)) p1 abund_zoop_select <- zoop.hell[,c("MEED","TRMU","SYNC","ASPR","KETA","KELO")] frame3 <- data.frame(abund_zoop_select, env[,"PC3"]) gg3 <- gather(frame3, species, value, -env....PC3..) p3 <- ggplot(gg3, aes(y = value, x = env....PC3..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(6)) + xlab("PC3 (Dystrophy/ Order/ CA)") + ylab("Zooplankton (hellinger)") p3 <- p3 + guides(colour = guide_legend(nrow = 4)) p3 abund_zoop_select <- zoop.hell[,c("DASC","LEMI","DAGA","COCO","KECO","TRMU","FILO", "KETA","DIBT")] frame4 <- data.frame(abund_zoop_select, env[,"PC4"]) gg4 <- gather(frame4, species, value, -env....PC4..) p4 <- ggplot(gg4, aes(y = value, x = env....PC4..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(9)) + xlab("PC4 (Size/ Sulfate)") + ylab("Zooplankton (hellinger)") p4 <- p4 + guides(colour = guide_legend(nrow = 4)) p4 abund_zoop_select <- zoop.hell[,c("CERE","CYSC","KEHI")] frame2 <- data.frame(abund_zoop_select, env[,"PC2"]) gg2 <- gather(frame2, species, value, -env....PC2..) p2 <- ggplot(gg2, aes(y = value, x = env....PC2..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(3)) + xlab("PCA2 (Depth/ Vol/ MEI)") + ylab("Zooplankton (hellinger)") p2 <- p2 + guides(colour = guide_legend(nrow = 4)) p2 abund_zoop_select <- zoop.hell[,c("KECO","GAST","KETA")] frame5 <- data.frame(abund_zoop_select, env[,"PC5"]) gg5 <- gather(frame5, species, value, -env....PC5..) p5 <- ggplot(gg5, aes(y = value, x = env....PC5..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(3)) + xlab("PC5 (Altitude/ Dev Vol)") + ylab("Zooplankton (hellinger)") p5 <- p5 + guides(colour = guide_legend(nrow = 4)) p5 grid.arrange(p1, p3, p4, p2, p5, ncol= 5) # Fish plots clrs.hcl <- function(n) { hcl(h = seq(8,8, length.out = n), c = 70, l = seq(10, 90, length.out = n), fixup = TRUE) } abund_fish_select <- fish.hell[,c("STIVI","ANRO","SAFO","CYPR")] frame1 <- data.frame(abund_fish_select, env[,"PC1"]) gg1 <- gather(frame1, species, value, -env....PC1..) p1 <- ggplot(gg1, aes(y = value, x = env....PC1..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(4)) + xlab("PC1 (Alkalinity/ Acidity)") + ylab("Fish (hellinger)") p1 <- p1 + guides(colour = guide_legend(nrow = 4)) p1 abund_fish_select <- fish.hell[,c("STIVI","ESLU","SAFO")] frame2 <- data.frame(abund_fish_select, env[,"PC2"]) gg2 <- gather(frame2, species, value, -env....PC2..) p2 <- ggplot(gg2, aes(y = value, x = env....PC2..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(3)) + xlab("PC2 (Depth/ Vol/ MEI)") + ylab("Fish (hellinger)") p2 <- p2 + guides(colour = guide_legend(nrow = 4)) p2 abund_fish_select <- fish.hell[,c("STIVI","CACA","SAFO")] frame4 <- data.frame(abund_fish_select, env[,"PC4"]) gg4 <- gather(frame4, species, value, -env....PC4..) p4 <- ggplot(gg4, aes(y = value, x = env....PC4..)) + theme_bw() + theme(legend.position="bottom",legend.text=element_text(size=6)) + geom_smooth(method = "glm",formula = y ~ poly(x,2), aes(group = species,colour=species), size = 1.2, se = FALSE) + scale_colour_manual(values = clrs.hcl(3)) + xlab("PC4 (Size/ Sulfate)") + ylab("Fish (hellinger)") p4 <- p4 + guides(colour = guide_legend(nrow = 4)) p4 grid.arrange(p1, p2, p4, ncol= 3) # 3.1) BDtotal, SCBD and LCBD - species & interaction network #### # Build interaction matrix among sp pairs based on co-occurrence # where co-occurrence = possible competition btw sp construire.inter <- function(mat) { mat = ifelse(mat>0,1,0) n = nrow(mat) p = ncol(mat) out = matrix(NA,n,p*(p-1)/2) pairs = NA rownames(out) = rownames(mat) no.col = 0 for(i in 1:(p-1)) { for(j in (i+1):p) { no.col = no.col+1 # cat("no.col =",no.col,"\n") out[,no.col] = mat[,i] * mat[,j] pairs = c(pairs,paste("Sp",i,j,sep=".")) } } colnames(out) = pairs[-1] sp.sum = apply(out,2,sum) out[,(sp.sum>0)] } inter.phyto = construire.inter(phyto) inter.zoop = construire.inter(zoop) inter.fish = construire.inter(fish) # 1. Calculate BDtotal, SCBD et LCBD for phyto (beta.res1 = beta.div(phyto,nperm=9999)) # Hellinger par default # $beta # SStotal BDtotal ### DBtotal falls btw [0, 1] # 28.3437857 0.6161693 (res1.signif.LCBD = which(beta.res1$p.LCBD <= 0.05)) # [1] 17 28 31 39 # Maps of LCBD values and richness per quadrat QC.XY = geoXY(data_mNAr$Lat/1000,data_mNAr$Long/1000) plot(QC.XY, asp=1, type="n", xlab="Longitude (m)", ylab="Latitude (m)", main="LCBD indices, phytoplankton") points(QC.XY,pch=21, col="white", bg="slategray3", cex=18*sqrt(beta.res1$LCBD)) points(QC.XY[res1.signif.LCBD,],pch=21, col="white", bg="darkgreen", cex=18*sqrt(beta.res1$LCBD[res1.signif.LCBD])) text(QC.XY, labels=1:47, pos=4, cex=0.8, offset=0.7) data_mNAr$col <- beta.res1$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0 data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "green" data_mNAr[data_mNAr$col == 1,]$col <- "blue" p_phyto = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "white",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("LCBD indices, phytoplankton") + geom_point(aes(size=beta.res1$LCBD,fill=col),colour="white",shape=21) + scale_fill_manual(values=c("slategray3","darkgreen")) p_phyto <- p_phyto + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_phyto # 2. Calculate BDtotal, SCBD et SLCBD for phyto sp pairs (beta.res2 = beta.div(inter.phyto)) # Hellinger by default # $beta # SStotal BDtotal # 37.1095252 0.8067288 <= BDtotal greater than that of just community above (res2.signif.LCBD = which(beta.res2$p.LCBD <= 0.05)) # [1] 2 4 7 8 17 18 25 30 31 33 38 39 40 42 44 45 # Maps of LCBD values and richness per quadrat plot(QC.XY, asp=1, type="n", xlab="Longitude (m)", ylab="Latitude (m)", main="LCBD indices, pairs of phytoplankton") points(QC.XY,pch=21, col="white", bg="slategray3", cex=18*sqrt(beta.res2$LCBD)) points(QC.XY[res2.signif.LCBD,],pch=21, col="white", bg="darkgreen", cex=18*sqrt(beta.res2$LCBD[res2.signif.LCBD])) text(QC.XY, labels=1:47, pos=4, cex=0.8, offset=0.7) data_mNAr$col <- beta.res2$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0 data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "green" data_mNAr[data_mNAr$col == 1,]$col <- "blue" p_phyto_pairs = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "white",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("LCBD indices, pairs of phytoplankton") + geom_point(aes(size=beta.res2$LCBD,fill=col),colour="white",shape=21) + scale_fill_manual(values=c("slategray3","darkgreen")) p_phyto_pairs <- p_phyto_pairs + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_phyto_pairs # 3. BDtotal, SCBD et SLCBD for zoop (beta.res3 = beta.div(zoop)) # Hellinger # $SStotal_BDtotal # [1] 20.9507567 0.4554512 (res3.signif.LCBD = which(beta.res3$p.LCBD <= 0.05)) # [1] 4 9 10 13 23 42 # Maps of LCBD values and richness per quadrat plot(QC.XY, asp=1, type="n", xlab="Longitude (m)", ylab="Latitude (m)", main="LCBD indices, pairs of phytoplankton") points(QC.XY,pch=21, col="white", bg="slategray3", cex=18*sqrt(beta.res3$LCBD)) points(QC.XY[res3.signif.LCBD,],pch=21, col="white", bg="royal blue", cex=18*sqrt(beta.res3$LCBD[res3.signif.LCBD])) text(QC.XY, labels=1:47, pos=4, cex=0.8, offset=0.7) data_mNAr$col <- beta.res3$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0 data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "royalblue" data_mNAr[data_mNAr$col == 1,]$col <- "blue" p_zoop = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "white",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("LCBD indices, zooplankton") + geom_point(aes(size=beta.res3$LCBD,fill=col),colour="white",shape=21) + scale_fill_manual(values=c("slategray3","royalblue")) p_zoop <- p_zoop + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_zoop # 4. BDtotal, SCBD et SLCBD for zoop sp pairs (beta.res4 = beta.div(inter.zoop)) # Hellinger # $SStotal_BDtotal # [1] 31.9332797 0.6942017 <= BDtotal greater again than community (res4.signif.LCBD = which(beta.res4$p.LCBD <= 0.05)) # [1] 1 3 4 5 6 8 9 22 23 47 # Maps of LCBD values and richness per quadrat plot(QC.XY, asp=1, type="n", xlab="Longitude (m)", ylab="Latitude (m)", main="LCBD indices, zooplankton") points(QC.XY,pch=21, col="white", bg="slategray3", cex=18*sqrt(beta.res4$LCBD)) points(QC.XY[res4.signif.LCBD,],pch=21, col="white", bg="royal blue", cex=18*sqrt(beta.res4$LCBD[res4.signif.LCBD])) text(QC.XY, labels=1:47, pos=4, cex=0.8, offset=0.7) data_mNAr$col <- beta.res4$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0 data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "royalblue" data_mNAr[data_mNAr$col == 1,]$col <- "blue" p_zoop_pairs = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "white",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("LCBD indices, pairs of zooplankton") + geom_point(aes(size=beta.res4$LCBD,fill=col),colour="white",shape=21) + scale_fill_manual(values=c("slategray3","royal blue")) p_zoop_pairs <- p_zoop_pairs + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_zoop_pairs # 5. BDtotal, SCBD et SLCBD for fish (beta.res5 = beta.div(fish)) # $SStotal_BDtotal # [1] 32.1433251 0.6987679 (res5.signif.LCBD = which(beta.res5$p.LCBD <= 0.05)) # [1] 6 12 # Maps of LCBD values and richness per quadrat plot(QC.XY, asp=1, type="n", xlab="Longitude (m)", ylab="Latitude (m)", main="LCBD indices, fish") points(QC.XY,pch=21, col="white", bg="slategray3", cex=18*sqrt(beta.res5$LCBD)) points(QC.XY[res5.signif.LCBD,],pch=21, col="white", bg="firebrick", cex=18*sqrt(beta.res5$LCBD[res5.signif.LCBD])) text(QC.XY, labels=1:47, pos=4, cex=0.8, offset=0.7) data_mNAr$col <- beta.res5$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0 data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "red" data_mNAr[data_mNAr$col == 1,]$col <- "blue" p_fish = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "white",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("LCBD indices, fish") + geom_point(aes(size=beta.res5$LCBD,fill=col),colour="white",shape=21) + scale_fill_manual(values=c("slategray3","firebrick")) p_fish <- p_fish + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_fish # 6. BDtotal, SCBD et SLCBD for fish sp pairs (beta.res6 = beta.div(inter.fish)) # $SStotal_BDtotal # [1] 31.1310416 0.6767618 <= BDtotal comparable to BDtotal of the community (res6.signif.LCBD = which(beta.res6$p.LCBD <= 0.05)) # [1] 2 3 4 5 6 7 8 9 10 12 14 15 16 17 21 22 25 28 30 31 32 38 39 43 44 45 46 # [28] 47 # Maps of LCBD values and richness per quadrat plot(QC.XY, asp=1, type="n", xlab="Longitude (m)", ylab="Latitude (m)", main="LCBD indices, pairs of fish") points(QC.XY,pch=21, col="white", bg="slategray3", cex=18*sqrt(beta.res6$LCBD)) points(QC.XY[res6.signif.LCBD,],pch=21, col="white", bg="firebrick", cex=18*sqrt(beta.res6$LCBD[res6.signif.LCBD])) text(QC.XY, labels=1:47, pos=4, cex=0.8, offset=0.7) data_mNAr$col <- beta.res6$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0 data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "red" data_mNAr[data_mNAr$col == 1,]$col <- "blue" p_fish_pairs = ggplot(data=data_mNAr, aes(x=(-1)*Long, y=Lat)) + coord_map(xlim=xlim,ylim=ylim) + geom_polygon(data=worldmap,aes(X,Y,group=PID), fill = "white",color="slategray4") + labs(x="Longitude (DD)",y="Latitude (DD)") + theme_grey() + ggtitle("LCBD indices, pairs of fish") + geom_point(aes(size=beta.res6$LCBD,fill=col),colour="white",shape=21) + scale_fill_manual(values=c("slategray3","firebrick")) p_fish_pairs <- p_fish_pairs + theme(panel.background = element_rect(fill = "slategray3", colour = "slategray3", size = 0.5, linetype = "solid"), panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "slategray3"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "slategray3")) p_fish_pairs quartz() # or windows() grid.arrange(p_phyto,p_phyto_pairs,p_zoop,p_zoop_pairs,p_fish,p_fish_pairs,ncol=4) # Distribution of SCBDs par(mfrow=c(3,2)) hist(beta.res1$SCBD,xlab="SCBD for phytoplankton",main="", col="slategray3",border="slategray4") # 87 phyto sp hist(beta.res2$SCBD,xlab="SCBD for pairs of phytopankton", main="", col="slategray3",border="slategray4") # 3741 phyto sp pairs, of 2875 pairs (77%) are interacting # The other pairs were eliminated by construire.inter() hist(beta.res3$SCBD,xlab="SCBD for zooplankton",main="", col="slategray3",border="slategray4") # 34 zoop sp hist(beta.res4$SCBD,xlab="SCBD for pairs of zooplankton",main="", col="slategray3",border="slategray4") # 561 pairs, of which 497 pairs (89%) are interacting hist(beta.res5$SCBD,xlab="SCBD for fish",main="", col="slategray3",border="slategray4") # 18 fish sp hist(beta.res6$SCBD,xlab="SCBD for pairs of fish",main="", col="slategray3",border="slategray4") # 153 pairs, of which 54 pairs (35%) are interacting # SCBD plot for phytoplankton species keep.sp.1 = which(beta.res1$SCBD >= mean(beta.res1$SCBD)) # which species have SCBD greater than the mean SCBDall.1 <- as.data.frame(beta.res1$SCBD) SCBDall.1$Species <- rownames(SCBDall.1) SCBDall.1 <- arrange(SCBDall.1,desc(SCBDall.1[,1])) colnames(SCBDall.1) <- c("SCBD","Species") # SCBDall.1 <- as.data.frame(SCBDall.1[SCBDall.1$SCBD > 0.0044,]) SCBD.1 <- ggplot() + geom_point(aes(x=reorder(SCBDall.1[,2],desc(SCBDall.1[,1])),y=SCBDall.1[,1]+0.01,size=log10(SCBDall.1[,1]+0.01),col=SCBDall.1[,1] >= mean(SCBDall.1[,1]))) + theme(axis.text.x = element_text(colour="black", size=6,angle=90), legend.position = 'none') + ylab("log10(SCBD)") + xlab("Phytoplankton species") + scale_y_log10() SCBD.1 # SCBD plot for zooplankton species keep.sp.3 = which(beta.res3$SCBD >= mean(beta.res3$SCBD)) # which species have SCBD greater than the mean SCBDall.3 <- as.data.frame(beta.res3$SCBD) SCBDall.3$Species <- rownames(SCBDall.3) SCBDall.3 <- arrange(SCBDall.3,desc(SCBDall.3[,1])) colnames(SCBDall.3) <- c("SCBD","Species") # SCBDall.3 <- as.data.frame(SCBDall.3[SCBDall.3$SCBD > 0.0044,]) SCBD.3 <- ggplot() + geom_point(aes(x=reorder(SCBDall.3[,2],desc(SCBDall.3[,1])),y=SCBDall.3[,1]+0.01,size=log10(SCBDall.3[,1]+0.01),col=SCBDall.3[,1] >= mean(SCBDall.3[,1]))) + theme(axis.text.x = element_text(colour="black", size=6,angle=90), legend.position = 'none') + ylab("log10(SCBD)") + xlab("Zooplankton species") + scale_y_log10() SCBD.3 # SCBD plot for fish species keep.sp.5 = which(beta.res5$SCBD >= mean(beta.res5$SCBD)) # which species have SCBD greater than the mean SCBDall.5 <- as.data.frame(beta.res5$SCBD) SCBDall.5$Species <- rownames(SCBDall.5) SCBDall.5 <- arrange(SCBDall.5,desc(SCBDall.5[,1])) colnames(SCBDall.5) <- c("SCBD","Species") # SCBDall.5 <- as.data.frame(SCBDall.5[SCBDall.5$SCBD > 0.0044,]) SCBD.5 <- ggplot() + geom_point(aes(x=reorder(SCBDall.5[,2],desc(SCBDall.5[,1])),y=SCBDall.5[,1]+0.01,size=log10(SCBDall.5[,1]+0.01),col=SCBDall.5[,1] >= mean(SCBDall.5[,1]))) + theme(axis.text.x = element_text(colour="black", size=6,angle=90), legend.position = 'none') + ylab("log10(SCBD)") + xlab("Fish species") + scale_y_log10() SCBD.5 # Compare LCBD (community) and LCBD-interaction # Phytoplankton par(mfrow=c(1,3)) plot(beta.res1$LCBD, beta.res2$LCBD,xlab="LCBD phyto",ylab="LCBD pairs phyto", cex=2, pch=21, bg="slategray3", col="white") cor.test(beta.res1$LCBD, beta.res2$LCBD) # Zooplankton plot(beta.res3$LCBD, beta.res4$LCBD,xlab="LCBD zoop",ylab="LCBD pairs zoop", cex=2, pch=21, bg="slategray3", col="white") cor.test(beta.res3$LCBD, beta.res4$LCBD) # Fish plot(beta.res5$LCBD, beta.res6$LCBD,xlab="LCBD fish",ylab="LCBD pairs fish", cex=2, pch=21, bg="slategray3", col="white") cor.test(beta.res5$LCBD, beta.res6$LCBD) # Interpreting LCBD for the phytoplankton phyto.pa = ifelse(phyto>0,1,0) rich.phyto = apply(phyto.pa, 1, sum) # Richness = [8,39] cor(rich.phyto, beta.res1$LCBD) # [1] -0.51351 plot(rich.phyto, beta.res1$LCBD,cex=2, pch=21, bg="slategray3", col="white", xlab="Phytoplancton richness",ylab="LCBD phytoplancton") cbind(rich.phyto[res1.signif.LCBD], beta.res1$LCBD[res1.signif.LCBD]) # Sites with elevated LCBDs have a low sp richness reg.out1 <- lm(beta.res1$LCBD ~ ., env) reg.out1 <- lm(beta.res1$LCBD ~ PC1, env) summary(reg.out1) # => LCBD increases with PC1. R2 is weak (= 0.16). # Lakes with elevated LCBD have a low richness and tend to be more alkaline and acidic env2 <- data_mNAr[,c(123:154,156:163)] env2 <- as.data.frame(scale(env2)) m0 = lm(beta.res1$LCBD ~ 1, data=env2) mtot = lm(beta.res1$LCBD ~ ., data=env2) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) # => LCBD increases with CA, H, Long, Na and decreases with Tr et Lmax # Lakes with elevated LCBD have a low richness and tend to be more alkaline and acidic # but also more transparent and situated in the eastern part of the province # R2 = 0.60 # Interprating LCBD, phytoplankton sp pairs inter.phyto.pa = ifelse(inter.phyto>0,1,0) rich.inter.phyto = apply(inter.phyto.pa, 1, sum) # Richness = [28, 741] cor(rich.inter.phyto, beta.res2$LCBD) # [1] -0.6770544 plot(rich.inter.phyto, beta.res2$LCBD,cex=2, pch=21, bg="slategray3", col="white", xlab="Phytoplancton pairs richness",ylab="LCBD interaction phytoplancton") cbind(rich.inter.phyto[res2.signif.LCBD], beta.res2$LCBD[res2.signif.LCBD]) # Sites with elevated LCBDs have a low sp richness reg.out2 <- lm(beta.res2$LCBD ~ ., env) summary(reg.out2) # => LCBD.interaction not related to env PCs # R2 very weak (= 0.07). m0 = lm(beta.res2$LCBD ~ 1, data=env2) mtot = lm(beta.res2$LCBD ~ ., data=env2) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) # Interpreting LCBD, zooplankton zoop.pa = ifelse(zoop>0,1,0) rich.zoop = apply(zoop.pa, 1, sum) # Richness = [4,18] cor.test(rich.zoop, beta.res3$LCBD) # [1] -0.33734 plot(rich.zoop, beta.res3$LCBD,cex=2, pch=21, bg="slategray3", col="white", xlab="Zooplancton richness",ylab="LCBD zooplancton") cbind(rich.zoop[res3.signif.LCBD], beta.res3$LCBD[res3.signif.LCBD]) # Sites with elevated LCBDs tend to have a low sp richness reg.out3 = lm(beta.res3$LCBD ~ PC1, env) summary(reg.out3) # => LCBD decreases with PC1 (alkalinity) # R2 = 0.31 m0 = lm(beta.res3$LCBD ~ 1, data=env2) mtot = lm(beta.res3$LCBD ~ ., data=env2) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) # LCBD increases with MEI, Pb, Order, Alc-SO4 et Zn # R2 = 0.53 # Interpretating LCBD, zooplankton sp pairs inter.zoop.pa = ifelse(inter.zoop>0,1,0) rich.inter.zoop = apply(inter.zoop.pa, 1, sum) # Richness = [6,153] cor(rich.inter.zoop, beta.res4$LCBD) # [1] -0.336159 plot(rich.inter.zoop, beta.res4$LCBD,cex=2, pch=21, bg="slategray3", col="white", xlab="Zooplancton pairs richness",ylab="LCBD interaction zooplancton") cbind(rich.inter.zoop[res4.signif.LCBD], beta.res4$LCBD[res4.signif.LCBD]) # Sites with elevated LCBD-interactions tend to have a low sp pair richness reg.out4 = lm(beta.res4$LCBD ~ ., env[,c(1,5)]) summary(reg.out4) # => LCBD.interactions related to PC1 (alkalinity) and PC5 (lake dev vol) # R2 = 0.54 m0 = lm(beta.res4$LCBD ~ 1, data=env2) mtot = lm(beta.res4$LCBD ~ ., data=env2) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) # => LCBD.interaction increases with SO4, MEI # R2 = 0.68 # Interpretating LCBD, fish fish.pa = ifelse(fish>0,1,0) rich.fish = apply(fish.pa, 1, sum) # Richness = [1,6] cor.test(rich.fish, beta.res5$LCBD) # [1] -0.1770612 plot(rich.fish, beta.res5$LCBD,cex=2, pch=21, bg="slategray3", col="white", xlab="Fish richness",ylab="LCBD fish") cbind(rich.fish[res5.signif.LCBD], beta.res5$LCBD[res5.signif.LCBD]) # Sites with elevated LCBDs tend to have a low sp richness but large variance (wedge shaped relationship) reg.out5 = lm(beta.res5$LCBD ~ ., env) summary(reg.out5) # => LCBD not related to PCs, R2 = 0.17. m0 = lm(beta.res5$LCBD ~ 1, data=env2) mtot = lm(beta.res5$LCBD ~ ., data=env2) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) # LCBD increases with Lat, Long, Pb, NO3, Order, Zmax, Alt and decreases with Al # R2 = 0.80 # Interpretating LCBD, fish sp pairs inter.fish.pa = ifelse(inter.fish>0,1,0) rich.inter.fish = apply(inter.fish.pa, 1, sum) # Richness = [0,15] cor.test(rich.inter.fish, beta.res6$LCBD) # [1] 0.5122073 plot(rich.inter.fish, beta.res6$LCBD,cex=2, pch=21, bg="slategray3", col="white", xlab="Fish pairs richness",ylab="LCBD interactions fish") cbind(rich.inter.fish[res6.signif.LCBD], beta.res6$LCBD[res6.signif.LCBD]) reg.out6 = lm(beta.res6$LCBD ~ PC1 + PC5, env) summary(reg.out6) # => LCBD.interactions related to PC1 (alkalinity) and PC5 (lake dev vol) # R2 = 0.26 m0 = lm(beta.res6$LCBD ~ 1, data=env2) mtot = lm(beta.res6$LCBD ~ ., data=env2) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) # => LCBD.interactions increases with Wmax and H, but decreases with Tr and Lat # R2 = 0.64 # ________________ Additional analyses _______________ #### # MFA with PC1-PC5 as env vars #### tab1 <- cbind(phyto.esp[-c(7,12),],zoop.esp[-c(7,12),],fish.esp[-c(7,12),],scores(pca.env,choices=1:5)$sites[,1:5]) grn <- c(ncol(phyto.esp[-c(7,12),]),ncol(zoop.esp[-c(7,12),]),ncol(fish.esp[-c(7,12),]),ncol(scores(pca.env,choices=1:5)$sites[,1:5])) QC.mfa <- MFA(tab1, group=grn, type=c("c","c","c","s"), ncp=5, name.group=c("Phyto","Zoop","Fish","Env")) summary(QC.mfa) plot(QC.mfa,choix = "var",habillage = "group", shadowtext = FALSE) # Select the most characteristic variables aa <- dimdesc(QC.mfa, axes = 1:2, proba = 0.0001) # Plot varsig <- QC.mfa$quanti.var$cor[unique(c(rownames(aa$Dim.1$quanti), rownames(aa$Dim.2$quanti))), ] plot(varsig[, 1:2],asp = 1, type = "n", xlim = c(-1, 1), ylim = c(-1, 1)) abline(h = 0, lty = 3) abline(v = 0, lty = 3) symbols(0, 0, circles = 1, inches = FALSE, add = TRUE) arrows(0, 0, varsig[, 1], varsig[, 2], length = 0.08, angle = 20) for (v in 1 : nrow(varsig)) { if (abs(varsig[v, 1]) > abs(varsig[v, 2])) { if (varsig[v, 1] >= 0) pos <- 4 else pos <- 2 } else { if (varsig[v, 2] >= 0) pos <- 3 else pos <- 1 } text(varsig[v, 1], varsig[v, 2], labels = rownames(varsig)[v], pos = pos) } # Alternative plots # Partial axes plot(QC.mfa, choix="axes") fviz_mfa_axes(QC.mfa) # Numeric variables (correlation circle) plot(QC.mfa, choix="var") # Sites (global and partial viewpoints) plot(QC.mfa, choix="ind", habillage="group", partial="all", invisible="quali") # Proportion of variance explained eig.val <- get_eigenvalue(QC.mfa) fviz_screeplot(QC.mfa) ev <- QC.mfa$eig[, 1] names(ev) <- paste("MFA", 1 : length(ev)) screestick(ev, las = 2) # Groups of variables group <- get_mfa_var(QC.mfa, "group") fviz_mfa_var(QC.mfa, "group") # Contribution dimensions fviz_contrib(QC.mfa, "group", axes = 1) fviz_contrib(QC.mfa, "group", axes = 2) fviz_contrib(QC.mfa, "group", axes = 3) fviz_contrib(QC.mfa, "group", axes = 4) fviz_contrib(QC.mfa, "group", axes = 5) # Quantitative variables quanti.var <- get_mfa_var(QC.mfa, "quanti.var") quanti.var fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE) fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE,axes=c(2,3)) fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE,axes=c(3,4)) fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE,axes=c(4,5)) # Quantitative variables colored by groups plot(QC.mfa,choix="var",habillage="group",cex=0.8,select="contrib 14",axes=c(1,2)) # Plotting without fviz (so that arrows are equal) plot(QC.mfa,choix = "var",habillage = "group", shadowtext = FALSE,axes=c(2,3)) plot(QC.mfa,choix = "var",habillage = "group", shadowtext = FALSE,axes=c(3,4)) plot(QC.mfa,choix = "var",habillage = "group", shadowtext = FALSE,axes=c(4,5)) fviz_mfa_var(QC.mfa, "quanti.var", palette = "jco", col.var.sup = "violet", repel = TRUE, geom = c("point", "text"), legend = "bottom") # same plot but with points instead of arrows # Contribution of ind vars to dimensions fviz_contrib(QC.mfa, choice = "quanti.var", axes = 1, top = 20, palette = "jco") fviz_contrib(QC.mfa, choice = "quanti.var", axes = 2, top = 20, palette = "jco") fviz_contrib(QC.mfa, choice = "quanti.var", axes = 3, top = 20, palette = "jco") fviz_contrib(QC.mfa, choice = "quanti.var", axes = 4, top = 20, palette = "jco") fviz_contrib(QC.mfa, choice = "quanti.var", axes = 5, top = 20, palette = "jco") fviz_mfa_var(QC.mfa, "quanti.var", col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), col.var.sup = "violet", repel = TRUE, geom = c("point", "text")) fviz_mfa_var(QC.mfa, "quanti.var", col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), col.var.sup = "violet", repel = TRUE) # Graph of sites by dimensions fviz_mfa_ind(QC.mfa, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE) fviz_mfa_ind(QC.mfa, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE,axes=c(2,3)) fviz_mfa_ind(QC.mfa, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE,axes=c(3,4)) fviz_mfa_ind(QC.mfa, col.ind = "cos2",gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE,axes=c(4,5)) fviz_mfa_ind(QC.mfa, partial = c("231B", "20B", "230","228","229","215B","216","164")) # Variance partitioning #### # Construct dbMEM variables XY <- geoXY(latitude=data_m$Lat,longitude=data_m$Long,unit=1000) lacs.dbmem.tmp <- dbmem(XY, silent = FALSE) lacs.dbmem <- as.data.frame(lacs.dbmem.tmp) # Truncation distance used above: (thr <- give.thresh(dist(XY))) # Display and count the eigenvalues attributes(lacs.dbmem.tmp)$values length(attributes(lacs.dbmem.tmp)$values) # Test for broad spatial trend in spe data, if present, must detrend, otherwise the first few odd dbMEM # will be lost to simply reconstruct this spatial trend. If detrend first, then the dbMEMs will be # available to look for other, more fine scale spatial patterns (so a more interesting use of the dbMEMs) spe.all <- cbind(phyto.esp,zoop.esp,fish.esp) # spe.all <- phyto.esp # phyto dominant # spe.all <- zoop.esp # zoop dominant # spe.all <- fish.esp # fish dominant dbmem.quick <- quickMEM(spe.all,XY) XY.rda <- rda(spe.all, XY) anova(XY.rda) # Is there a linear trend in the spe data? # Result: yes, significant trend # Computation of linearly detrended data spe.all.det <- resid(XY.rda) # Selection of env vars env.rda <- rda(spe.all.det[-c(7,12),] ~., env.dis[-c(7,12),]) # rows 7 and 12 have NA for MEI and DevVol; delete across all datasets (env.R2a <- RsquareAdj(env.rda)$adj.r.squared) env.fwd <- forward.sel(spe.all.det[-c(7,12),], env.dis[-c(7,12),], adjR2thresh = env.R2a, nperm = 9999) env.backward <- ordistep(env.rda, permutations = how(nperm = 9999)) env.sign <- sort(env.backward$terminfo$ordered) env.red <- env.dis[ ,c(names(env.sign))] colnames(env.red) # Selection of dbMEM vars dbmem.rda <- rda(spe.all.det[-c(7,12),] ~ ., lacs.dbmem[-c(7,12),]) # to run on the *detrended* spe data # OR # dbmem.rda <- rda(spe.all[-c(7,12),] ~ . , lacs.dbmem[-c(7,12),]) # to run on the spe data anova(dbmem.rda) # Analysis is not significant... # and run a forward selection of the dbMEM variables (dbmem.R2a <- RsquareAdj(dbmem.rda)$adj.r.squared) (dbmem.fwd <- forward.sel(spe.all.det[-c(7,12),], as.matrix(lacs.dbmem[-c(7,12),]),adjR2thresh = dbmem.R2a)) dbmem.fwd (nb.sig.dbmem <- nrow(dbmem.fwd)) # Number of signif: 0 (dbmem.sign <- sort(dbmem.fwd[ ,2])) dbmem.red <- lacs.dbmem[ ,c(dbmem.sign)] # Run a backward selection dbmem.backward <- ordistep(dbmem.rda, permutations = how(nperm = 9999)) dbmem.sign <- sort(dbmem.backward$terminfo$ordered) dbmem.red <- lacs.dbmem[ ,c(names(dbmem.sign))] colnames(dbmem.red) # New dbMEM analysis with significant dbMEM variables (dbmem.rda2 <- rda(spe.all.det[-c(7,12),] ~ MEM3 + MEM4, data = lacs.dbmem[-c(7,12),])) (fwd.R2a <- RsquareAdj(dbmem.rda2)$adj.r.squared) anova(dbmem.rda2) (axes.test <- anova(dbmem.rda2, by = "axis")) # Number of significant axes (nb.ax <- length(which(axes.test[ , ncol(axes.test)] <= 0.05))) # Plot the significant canonical axes rda2.axes <- scores(dbmem.rda2, choices = c(1:nb.ax), display = "lc", scaling = 1) par(mfrow = c(1,nb.ax)) for(i in 1:nb.ax){ sr.value(XY[-c(7,12),], rda2.axes[ ,i], sub = paste("dbMEM",i), csub = 2) } # Other dbMEMs par(mfrow = c(1,ncol(lacs.dbmem))) for(i in 1:ncol(lacs.dbmem)){ sr.value(XY, lacs.dbmem[ ,i], sub = paste("MEM",i), csub = 2) } # Interpreting spatial variation # regression of the significant canonical axes on the environmental variables, rda2.axis1.env <- lm(rda2.axes[ ,1] ~ ., data = data_mNAr[,c("PC1","PC2","PC3","PC4","PC5")]) summary(rda2.axis1.env) # rda2.axis2.env <- lm(rda2.axes[ ,2] ~ ., data = data_mNAr[,c("PC1","PC2","PC3","PC4","PC5")]) # summary(rda2.axis2.env) # rda2.axis3.env <- lm(rda2.axes[ ,3] ~ ., data = data_mNAr[,c("PC1","PC2","PC3","PC4","PC5")]) # summary(rda2.axis3.env) # Scalogram # of the variance explained by all dbMEM eigenfunctions scalog(dbmem.rda) # Usually we would use the scalogram to separate the dbMEMs into broad, intermediate or small scale # (see p. 326 of NEwR), but here, with only 7 dbMEMs, and only 0 significant, # seems useless to bin into broad to fine. dbmem.broad <- lacs.dbmem[,1:3] dbmem.intermediate <- lacs.dbmem[,4:6] dbmem.fine <- lacs.dbmem[,7] # Combining dbMEM and variation partitioning (ponds.varpart <- varpart(spe.all.det[-c(7,12),], env.red[-c(7,12),], XY[-c(7,12),])) # Or using all env vars from the pars MFA (ponds.varpart <- varpart(spe.all.det[-c(7,12),], env.dis[-c(7,12),-c(9)], XY[-c(7,12),])) # Add in dbMEMs to see how things change (ponds.varpart <- varpart(spe.all.det[-c(7,12),], env.red[-c(7,12),], XY[-c(7,12),], lacs.dbmem[-c(7,12),3:4])) (ponds.varpart <- varpart(spe.all.det[-c(7,12),], env.red[-c(7,12),], XY[-c(7,12),], dbmem.broad[-c(7,12),], dbmem.intermediate[-c(7,12),])) # Show the symbols of the fractions and plot their values plot(ponds.varpart,digits = 2, bg = c("red", "yellow","blue")) # All species v <- venneuler(c(A=0.049, B=0.006, C=0.004, "A&B"=0,"B&C"=0, "A&C"=0,"A&B&C"=0.002)) v$colors <- c(.32, .6, .05, 0.15) v$labels <- c("PC3-4\n4.9%","dbMEM-med\n0.6%","dbMEM-fine\n0.4%\n\nAll\n0.2%") plot(v) # Fish v <- venneuler(c(A=0.108, B=0.002, C=0.007, "A&B"=0,"B&C"=0, "A&C"=0.003,"A&B&C"=0)) v$colors <- c(.32, .6, .05, 0.15) v$labels <- c("PC3-4\n10.8%","dbMEM-med\n0.2%","dbMEM-fine\n0.7%\n\ndbMEM-fine&PCs\n0.3%") plot(v) # Zoop v <- venneuler(c(A=0.016, B=0, C=0, "A&B"=0.005,"B&C"=0, "A&C"=0.012,"A&B&C"=0.002)) v$colors <- c(.32, .6, .05, 0.15) v$labels <- c("All\n0.2%\n\nPC4-5\n1.6%","dbMEM-broad&PCs\n0.5%","dbMEM-med&PCs\n1.2%") plot(v) # BDtotal, SCBD and LCBD - vs PC1 #### # LCBD phytoplankton vs PC1 data_mNAr$col <- beta.res1$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "darkolivegreen4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out1 <- lm(beta.res1$LCBD ~ PC1, env) summary(reg.out1)$adj.r.squared # [1] 0.0640529 plot(env$PC1, beta.res1$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="PC1",ylab="LCBD phytoplankton") abline(reg.out1) # LCBD pairs phytoplankton vs PC1 data_mNAr$col <- beta.res2$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "darkolivegreen4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out2 <- lm(beta.res2$LCBD ~ PC1, env) summary(reg.out2)$adj.r.squared # [1] -0.00752959 plot(env$PC1, beta.res2$LCBD,cex=2, pch=21, bg="slategray3", col="white", xlab="PC1",ylab="LCBD interaction phytoplancton") abline(reg.out2, lty=2, col="slategray3") # LCBD zooplankton vs PC1 data_mNAr$col <- beta.res3$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "deepskyblue4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out3 = lm(beta.res3$LCBD ~ PC1, env) summary(reg.out3)$adj.r.squared # [1] 0.1152872 plot(env$PC1, beta.res3$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="PC1",ylab="LCBD zooplankton") abline(reg.out3) # LCBD pairs zooplankton vs PC1 data_mNAr$col <- beta.res4$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "deepskyblue4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out4 = lm(beta.res4$LCBD ~ PC1, env) summary(reg.out4)$adj.r.squared # [1] 0.488109 plot(env$PC1, beta.res4$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="PC1",ylab="LCBD interaction zooplankton") abline(reg.out4) # LCBD fish vs PC1 data_mNAr$col <- beta.res5$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "firebrick"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out5 = lm(beta.res5$LCBD ~ PC1, env) summary(reg.out5)$adj.r.squared # [1] 0.05368598 plot(env$PC1, beta.res5$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="PC1",ylab="LCBD fish") abline(reg.out5,lty=2, col="slategray3") # LCBD pairs fish vs PC1 data_mNAr$col <- beta.res6$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "firebrick"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out6 = lm(beta.res6$LCBD[rich.inter.fish>0] ~ PC1, env[rich.inter.fish>0,]) summary(reg.out6)$adj.r.squared # [1] -0.0275324 plot(env$PC1[rich.inter.fish>0], beta.res6$LCBD[rich.inter.fish>0],cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="PC1",ylab="LCBD interactions fish") abline(reg.out6,lty=2, col="slategray3") # ______________________ - vs Spatial #### # LCBD phytoplankton vs spatial data_mNAr$col <- beta.res1$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "darkolivegreen4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" XY.rda <- rda(beta.res1$LCBD, XY[-c(7,12),]) anova(XY.rda) # Is there a linear trend in the spe data? # Result: No reg.out1 <- lm(beta.res1$LCBD ~ Lat, data_mNAr) summary(reg.out1) # NS reg.out1 <- lm(beta.res1$LCBD ~ Long, data_mNAr) summary(reg.out1) # Signif reg.out1 <- lm(beta.res1$LCBD ~ Altitude, data_mNAr) summary(reg.out1) # NS reg.out1 <- lm(beta.res1$LCBD ~ Long, data_mNAr) summary(reg.out1)$adj.r.squared # [1] 0.06772908 plot(data_mNAr$Long, beta.res1$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="Longitude",ylab="LCBD phytoplancton") abline(reg.out1) m0 = lm(beta.res1$LCBD ~ 1, data=lacs.dbmem[-c(7,12),]) mtot = lm(beta.res1$LCBD ~ ., data=lacs.dbmem[-c(7,12),]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) reg.out1 <- lm(beta.res1$LCBD ~ MEM3, lacs.dbmem[-c(7,12),]) summary(reg.out1)$adj.r.squared # [1] 0.1060835 plot(lacs.dbmem[-c(7,12),]$MEM3, beta.res1$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM3",ylab="LCBD phytoplankton") abline(reg.out1) # LCBD pairs phytoplankton vs spatial data_mNAr$col <- beta.res2$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "darkolivegreen4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" XY.rda <- rda(beta.res2$LCBD, XY[-c(7,12),]) anova(XY.rda) # Is there a linear trend in the spe data? # Result: No reg.out2 <- lm(beta.res2$LCBD ~ Lat, data_mNAr) summary(reg.out2) # NS reg.out2 <- lm(beta.res2$LCBD ~ Long, data_mNAr) summary(reg.out2) # NS reg.out2 <- lm(beta.res2$LCBD ~ Altitude, data_mNAr) summary(reg.out2) # NS m0 = lm(beta.res2$LCBD ~ 1, data=lacs.dbmem[-c(7,12),]) mtot = lm(beta.res2$LCBD ~ ., data=lacs.dbmem[-c(7,12),]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) reg.out2 <- lm(beta.res2$LCBD ~ MEM1, lacs.dbmem[-c(7,12),]) summary(reg.out2)$adj.r.squared # [1] 0.08269696 plot(lacs.dbmem[-c(7,12),]$MEM1, beta.res1$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM1",ylab="LCBD interaction phytoplankton") abline(reg.out2) # LCBD zooplancton vs spatial data_mNAr$col <- beta.res3$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "deepskyblue4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" XY.rda <- rda(beta.res3$LCBD, XY[-c(7,12),]) anova(XY.rda) # Is there a linear trend in the spe data? # Result: Yes # Computation of linearly detrended data beta.res3.det <- resid(XY.rda) reg.out3 <- lm(beta.res3$LCBD ~ Lat, data_mNAr) summary(reg.out3) # Signif reg.out3 <- lm(beta.res3$LCBD ~ Long, data_mNAr) summary(reg.out3) # Signif reg.out3 <- lm(beta.res3$LCBD ~ Altitude, data_mNAr) summary(reg.out3) # NS reg.out3 <- lm(beta.res3$LCBD ~ Long, data_mNAr) summary(reg.out3)$adj.r.squared # [1] 0.09602632 plot(data_mNAr$Long, beta.res3$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="Longitude",ylab="LCBD zooplankton") abline(reg.out3) reg.out3 <- lm(beta.res3$LCBD ~ Lat, data_mNAr) summary(reg.out3)$adj.r.squared # [1] 0.0823451 plot(data_mNAr$Lat, beta.res3$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="Latitude",ylab="LCBD zooplankton") abline(reg.out3) m0 = lm(beta.res3.det ~ 1, data=lacs.dbmem[-c(7,12),]) mtot = lm(beta.res3.det ~ ., data=lacs.dbmem[-c(7,12),]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) m0 = lm(beta.res3$LCBD ~ 1, data=lacs.dbmem[-c(7,12),]) mtot = lm(beta.res3$LCBD ~ ., data=lacs.dbmem[-c(7,12),]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) reg.out3 <- lm(beta.res3$LCBD ~ MEM4, lacs.dbmem[-c(7,12),]) summary(reg.out3)$adj.r.squared # [1] 0.139394 plot(lacs.dbmem[-c(7,12),]$MEM4, beta.res3$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM4",ylab="LCBD zooplankton") abline(reg.out3) # LCBD pairs zooplankton vs spatial data_mNAr$col <- beta.res4$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "deepskyblue4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" XY.rda <- rda(beta.res4$LCBD, XY[-c(7,12),]) anova(XY.rda) # Is there a linear trend in the spe data? # Result: Yes # Computation of linearly detrended data beta.res4.det <- resid(XY.rda) reg.out4 <- lm(beta.res4$LCBD ~ Lat, data_mNAr) summary(reg.out4) # Signif reg.out4 <- lm(beta.res4$LCBD ~ Long, data_mNAr) summary(reg.out4) # Signif reg.out4 <- lm(beta.res4$LCBD ~ Altitude, data_mNAr) summary(reg.out4) # NS reg.out4 <- lm(beta.res4$LCBD ~ Long, data_mNAr) summary(reg.out4)$adj.r.squared # [1] 0.09037523 plot(data_mNAr$Long, beta.res4$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="Longitude",ylab="LCBD interaction zooplankton") abline(reg.out4) reg.out4 <- lm(beta.res4$LCBD ~ Lat, data_mNAr) summary(reg.out4)$adj.r.squared # [1] 0.4612907 plot(data_mNAr$Lat, beta.res4$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="Latitude",ylab="LCBD interaction zooplankton") abline(reg.out4) m0 = lm(beta.res4.det ~ 1, data=lacs.dbmem[-c(7,12),]) mtot = lm(beta.res4.det ~ ., data=lacs.dbmem[-c(7,12),]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) m0 = lm(beta.res4$LCBD ~ 1, data=lacs.dbmem[-c(7,12),]) mtot = lm(beta.res4$LCBD ~ ., data=lacs.dbmem[-c(7,12),]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) reg.out4 <- lm(beta.res4$LCBD ~ MEM4, lacs.dbmem[-c(7,12),]) summary(reg.out4)$adj.r.squared # [1] 0.1745702 plot(lacs.dbmem[-c(7,12),]$MEM4, beta.res4$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM4",ylab="LCBD interaction zooplankton") abline(reg.out4) reg.out4 <- lm(beta.res4$LCBD ~ MEM2, lacs.dbmem[-c(7,12),]) summary(reg.out4)$adj.r.squared # [1] 0.2822928 plot(lacs.dbmem[-c(7,12),]$MEM2, beta.res4$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM2",ylab="LCBD interaction zooplankton") abline(reg.out4) # LCBD poisson vs spatial data_mNAr$col <- beta.res5$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "firebrick"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" XY.rda <- rda(beta.res5$LCBD, XY[-c(7,12),]) anova(XY.rda) # Is there a linear trend in the spe data? # Result: Yes # Computation of linearly detrended data beta.res5.det <- resid(XY.rda) reg.out5 <- lm(beta.res5$LCBD ~ Lat, data_mNAr) summary(reg.out5) # NS reg.out5 <- lm(beta.res5$LCBD ~ Long, data_mNAr) summary(reg.out5) # Signif reg.out5 <- lm(beta.res5$LCBD ~ Altitude, data_mNAr) summary(reg.out5) # NS reg.out5 <- lm(beta.res5$LCBD ~ Long, data_mNAr) summary(reg.out5)$adj.r.squared # [1] 0.1692258 plot(data_mNAr$Long, beta.res5$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="Longitude",ylab="LCBD fish") abline(reg.out5) m0 = lm(beta.res5.det ~ 1, data=lacs.dbmem[-c(7,12),]) mtot = lm(beta.res5.det ~ ., data=lacs.dbmem[-c(7,12),]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) m0 = lm(beta.res5$LCBD ~ 1, data=lacs.dbmem[-c(7,12),]) mtot = lm(beta.res5$LCBD ~ ., data=lacs.dbmem[-c(7,12),]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) reg.out5 <- lm(beta.res5$LCBD ~ MEM5, lacs.dbmem[-c(7,12),]) summary(reg.out5)$adj.r.squared # [1] 0.06589002 plot(lacs.dbmem[-c(7,12),]$MEM5, beta.res5$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM5",ylab="LCBD fish") abline(reg.out5) reg.out5 <- lm(beta.res5$LCBD ~ MEM1, lacs.dbmem[-c(7,12),]) summary(reg.out5)$adj.r.squared # [1] 0.140705 plot(lacs.dbmem[-c(7,12),]$MEM1, beta.res4$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM1",ylab="LCBD fish") abline(reg.out5) # LCBD pairs fish vs spatial data_mNAr$col <- beta.res6$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "firebrick"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" XY47 <- XY[-c(7,12),] XY.rda <- rda(beta.res6$LCBD[rich.inter.fish>0], XY47[rich.inter.fish>0,]) anova(XY.rda) # Is there a linear trend in the spe data? # Result: No reg.out6 <- lm(beta.res6$LCBD[rich.inter.fish>0] ~ Lat[rich.inter.fish>0], data_mNAr) summary(reg.out6) # NS reg.out6 <- lm(beta.res6$LCBD[rich.inter.fish>0] ~ Long[rich.inter.fish>0], data_mNAr) summary(reg.out6) # NS reg.out6 <- lm(beta.res6$LCBD[rich.inter.fish>0] ~ Altitude[rich.inter.fish>0], data_mNAr) summary(reg.out6) # NS lacs.dbmem47 <- lacs.dbmem[-c(7,12),] m0 = lm(beta.res6$LCBD[rich.inter.fish>0] ~ 1, data=lacs.dbmem47[rich.inter.fish>0,]) mtot = lm(beta.res6$LCBD[rich.inter.fish>0] ~ ., data=lacs.dbmem47[rich.inter.fish>0,]) res.step = step(m0, scope=formula(mtot), direction="both") summary(res.step) # ______________________ - vs Spatial Fig #### quartz() par(mfrow=c(3,2)) data_mNAr$col <- beta.res1$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "darkolivegreen4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out1 <- lm(beta.res1$LCBD ~ MEM3, lacs.dbmem[-c(7,12),]) summary(reg.out1)$adj.r.squared # [1] 0.1060835 plot(lacs.dbmem[-c(7,12),]$MEM3, beta.res1$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM3",ylab="LCBD phytoplancton") abline(reg.out1) data_mNAr$col <- beta.res2$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "darkolivegreen4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out2 <- lm(beta.res2$LCBD ~ MEM1, lacs.dbmem[-c(7,12),]) summary(reg.out2)$adj.r.squared # [1] 0.08269696 plot(lacs.dbmem[-c(7,12),]$MEM1, beta.res1$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM1",ylab="LCBD interaction phytoplancton") abline(reg.out2) data_mNAr$col <- beta.res3$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "deepskyblue4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out3 <- lm(beta.res3$LCBD ~ MEM4, lacs.dbmem[-c(7,12),]) summary(reg.out3)$adj.r.squared # [1] 0.139394 plot(lacs.dbmem[-c(7,12),]$MEM4, beta.res3$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM4",ylab="LCBD zooplankton") abline(reg.out3) data_mNAr$col <- beta.res4$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "deepskyblue4"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out4 <- lm(beta.res4$LCBD ~ MEM2, lacs.dbmem[-c(7,12),]) summary(reg.out4)$adj.r.squared # [1] 0.2822928 plot(lacs.dbmem[-c(7,12),]$MEM2, beta.res4$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="MEM2",ylab="LCBD interaction zooplankton") abline(reg.out4) data_mNAr$col <- beta.res5$p.LCBD data_mNAr[data_mNAr$col <= 0.05,]$col <- 0; data_mNAr[data_mNAr$col > 0.05,]$col <- 1 data_mNAr[data_mNAr$col == 0,]$col <- "firebrick"; data_mNAr[data_mNAr$col == 1,]$col <- "slategray3" reg.out5 <- lm(beta.res5$LCBD ~ Long, data_mNAr) summary(reg.out5)$adj.r.squared # [1] 0.1692258 plot(data_mNAr$Long, beta.res5$LCBD,cex=2, pch=21, bg=data_mNAr$col, col="white", xlab="Longitude",ylab="LCBD fish") abline(reg.out5) # ______________________ - Partial correlation #### pcor.test(beta.res1$LCBD,data_mNAr$PC1,data_mNAr$Long,method="kendall") pcor.test(beta.res1$LCBD,data_mNAr$PC1,lacs.dbmem[-c(7,12),]$MEM3,method="kendall") pcor.test(beta.res2$LCBD,data_mNAr$PC1,lacs.dbmem[-c(7,12),]$MEM1,method="kendall") pcor.test(beta.res3$LCBD,data_mNAr$PC1,data_mNAr$Long,method="kendall") pcor.test(beta.res3$LCBD,data_mNAr$PC1,data_mNAr$Lat,method="kendall") pcor.test(beta.res3$LCBD,data_mNAr$PC1,lacs.dbmem[-c(7,12),]$MEM4,method="kendall") pcor.test(beta.res4$LCBD,data_mNAr$PC1,data_mNAr$Long,method="kendall") pcor.test(beta.res4$LCBD,data_mNAr$PC1,data_mNAr$Lat,method="kendall") pcor.test(beta.res4$LCBD,data_mNAr$PC1,lacs.dbmem[-c(7,12),]$MEM2,method="kendall") pcor.test(beta.res4$LCBD,data_mNAr$PC5,lacs.dbmem[-c(7,12),]$MEM4,method="kendall")