library(RColorBrewer) library(mgcv) library(vegan) library(MuMIn) #Load diversity results alphaPlants <- readRDS("alphaResultsPlants.rds") alphaAnts <- readRDS("alphaResultsAnts.rds") #Load environmental variables (plants): envData <- read.table("Water_content.txt", header = T) #Load environmental data (Ants) abun_ants.raw <- read.csv("ant.freq.csv") # load abundance dataset abun_ants <- abun_ants.raw[,2:30] # identify the abundances invasionAnts <- abun_ants.raw$treatment row.names(abun_ants) <- abun_ants.raw$sites # identify the sites envir <- read.csv("env.data.csv", row.names = 1) # load abundance dataset envir <- envir[rownames(abun_ants),] envir$invasionAnts <- invasionAnts #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### PLANTS #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### Colores<-brewer.pal(8, "Set1") Colores <- Colores[c(-6, -8)] alphaPoints<-0.25 ColTransp<-c(rgb(t(col2rgb(Colores)/255), alpha=alphaPoints, maxColorValue=1)) varPlot <- c("FRic" #,"FEve" #, "FDiv" , "Rao") nombresVar <- c("FRic" # , "Functional evenness" # , "Functional divergence" , "Rao") namesMethods <- c("M1. TPD; individuals\nwithin populations", "M2. TPD;\npopulation mean", "M3. Classic;\npopulation mean", "M4. TPD; individuals\nacross populations", "M5a. TPD;\nspecies mean", "M6. Classic;\nspecies mean") lwdLines <- 5 #### Patterns along the gradient #### waterGrad <- envData$WaterCont waterGrad2 <- waterGrad**2 waterGrad3 <- waterGrad**3 pdf("Manuscript/Figures/FiguresR1/Figure_plants_gradient.pdf", height = 0.8*14*(length(varPlot)), width = 0.8*14*length(namesMethods), pointsize=50) par(mfrow=c(length(varPlot) , length(namesMethods)), mgp=c(1, 0.1, 0), mar = c(2, 3, 2, 0), oma=c(2,2,3,2)) #### #### plotNum <- 1 for(i in 1:(length(varPlot))){ #For each component of diversity one row for(j in 1:length(alphaPlants)){ # For each method, one column labX <- "Soil water content (%)" labY <- ifelse(j == 1, nombresVar[i], "") # alphaUse <- scale(alphaPlants[[i]][, j]) if(i != 3){ # FRichness and Rao alphaUse <- (alphaPlants[[j]][, i]) allPreds <- allPredsPlot <- c("waterGrad", "waterGrad2", "waterGrad3") modelo<-lm(alphaUse ~ waterGrad + waterGrad2 + waterGrad3, na.action="na.fail") modeloDredge<-dredge(modelo, rank = "AICc", subset = dc(waterGrad, waterGrad2, waterGrad3)) selectedModel <- get.models(modeloDredge, 1)[[1]] if(j == 1){ selectedVarsBestModel <- rownames(anova(selectedModel)) selectedVarsBestModel <- selectedVarsBestModel[ - length(selectedVarsBestModel)] if(length(selectedVarsBestModel) > 0){ nonSelectedVarsBestModel <- allPreds[- which(allPreds %in% selectedVarsBestModel)] }else{ nonSelectedVarsBestModel <- allPreds } selectedVarsModel <- rownames(anova(selectedModel)) selectedVarsModel <- selectedVarsModel[ - length(selectedVarsModel)] positionBestModel <- as.numeric(which(rowSums(!is.na(modeloDredge[, selectedVarsBestModel, drop=F])) == length(selectedVarsBestModel) & rowSums(!is.na(modeloDredge[, nonSelectedVarsBestModel, drop = F])) == 0 )) deltaBestModel <- modeloDredge[positionBestModel, "delta"] } else{ selectedVarsModel <- rownames(anova(selectedModel)) selectedVarsModel <- selectedVarsModel[ - length(selectedVarsModel)] if(length(selectedVarsBestModel) > 0){ positionBestModel <- as.numeric(which(rowSums(!is.na(modeloDredge[, selectedVarsBestModel, drop = F])) == length(selectedVarsBestModel) & rowSums(!is.na(modeloDredge[, nonSelectedVarsBestModel, drop = F])) == 0)) } else{ positionBestModel <- which(rowSums(!is.na(modeloDredge[, nonSelectedVarsBestModel])) == 0) } deltaBestModel <- modeloDredge[positionBestModel, "delta"] } if(length(selectedVarsModel) == 0){ leg <- vector("expression", 2) # One for R2, one for no predictors, one for deltaAIC leg[1]<- substitute(expression(R^2== MYVALUE), list( MYVALUE = max(0, round(summary(selectedModel)$r.sq,3))))[2] leg[2]<- substitute(expression(Delta[AICc]== MYVALUE), list( MYVALUE = round(deltaBestModel, 3)))[2] } else{ leg <- vector("expression", 2 + length(selectedVarsModel)) leg[1]<- substitute(expression(R^2== MYVALUE), list( MYVALUE = max(0, round(summary(selectedModel)$r.sq,3))))[2] leg[2]<- substitute(expression(Delta[AICc]== MYVALUE), list( MYVALUE = round(deltaBestModel, 3)))[2] } limY <- range(alphaUse); limY[2] <- limY[2] + 0.5 * (limY[2] - limY[1]) plot(alphaUse ~ waterGrad, ylab="", axes=F, xlab="", pch = 16, col=Colores[j], bg = "grey85" , ylim= limY, cex=1.5) mtext(side=1, adj=0.5, labX, outer=F, line=1.5) mtext(side=2, adj=0.5, labY, outer=F, line=1.8) axis(1, tcl=0.25, lwd=0.8) axis(2, las=1, tcl=0.25,lwd=0.8) box( bty = "l", lwd =0.8 ) posLeg <- "topright" legend(posLeg, leg, adj=0, bty="n", cex=1.4) if(i==1){ mtext(side=3, adj=0.5, paste0(namesMethods[j]), outer=F, line=0.5) } mtext(side=3, adj=0, paste0(letters[plotNum], ") "), outer=F, line=-0.5) plotNum <- plotNum + 1 waterGradSeq <- seq(min(waterGrad), max(waterGrad), length.out = 100) waterGradSeq2 <- waterGradSeq**2 waterGradSeq3 <- waterGradSeq**3 predMat <- data.frame(waterGrad = waterGradSeq, waterGrad2 = waterGradSeq2, waterGrad3 = waterGradSeq3) lines <- predict(selectedModel, newdata = predMat) lines(waterGradSeq, lines, lwd=lwdLines, col = Colores[j], lty= 1) } } } #### #### dev.off() ################## ################## ################## ################## ################## CORRELATIONS BETWEEN METHODS PLANTS ################## ################## ################## ################## cexLeg <- 1.3 cexMethod <- 1.4 varPlot <- c("FRic" #,"FEve" #, "FDiv" , "Rao") nombresVar <- c("FRic" # , "Functional evenness" # , "Functional divergence" , "Rao") namesMethods <- c("M1. TPD; individuals\nwithin populations", "M2. TPD;\npopulation mean", "M3. Classic;\npopulation mean", "M4. TPD; individuals\nacross populations", "M5a. TPD;\nspecies mean", "M6. Classic;\nspecies mean") pdf("Manuscript/Figures/FiguresR1/Figure_plants_correlationsAlphaFRic.pdf", width = 20, height = 10, pointsize=23) #### #### matFRicPlot <- matrix(1:36, nrow=6, ncol=6, byrow=T) matRaoPlot <- matrix((1:36)+36, nrow=6, ncol=6, byrow=T) matBlank <- matrix(0, nrow=6, ncol=1) layoutMat <- cbind(matFRicPlot, matBlank, matRaoPlot) widthsCols <- c(rep(1, 6), 0.5, rep(1, 6)) widthsCols <- widthsCols/sum(widthsCols) layout(layoutMat, widths = widthsCols) par(mgp=c(1, 0.1, 0), mar = c(0.1, 0.1, 0.1, 0.1), oma=c(1, 1, 3 , 1)) ##### FUNCTIONAL RICHNESS ---- rangeFRic <- rangeRao <- numeric() for(i in 1:length(alphaPlants)){ rangeFRic <- c(rangeFRic, alphaPlants[[i]]$FRic) rangeRao <- c(rangeRao, alphaPlants[[i]]$Rao) } rangeFRic <- c(0, 1.1 * max(rangeFRic)) rangeRao <- c(0, 1.1 * max(rangeRao)) for(i in 1:length(alphaPlants)){#Rows for(j in 1:length(alphaPlants)){#Columns if(i == j){#Diag xAxis <- alphaPlants[[i]]$FRic dens <- density(xAxis) plot(dens$x, dens$y, type="l", ylab="", axes=F, xlab="", xlim=rangeFRic, ylim = 1.5 * range(dens$y), col=Colores[i], lwd=lwdLines) points(xAxis, rep(0, length(xAxis)), bg=ColTransp[i], col=Colores[i], pch=21, cex=0.8) box(lwd =0.8) leg <- ifelse(i != 5, paste0(names(alphaPlants)[i]), "M5a") legend("topright", leg = leg, bty="n", cex =cexMethod) if(i==6){axis(1,tcl=0.25,lwd=0.8)} if(i == 1) mtext(side=3, adj=0, "a) FRic", outer=F, line=0.5, cex=1.5) } else{ if(j < i){#Lower trian xAxis <- (alphaPlants[[j]]$FRic) yAxis <- (alphaPlants[[i]]$FRic) densX <- density(xAxis) densY <- density(yAxis) plot(xAxis, yAxis, ylab="", axes=F, xlab="", pch=21, bg="grey80", col="grey30", xlim=rangeFRic, ylim=rangeFRic) box(lwd =0.8) # abline(lm(yAxis~xAxis), lty=2) abline(0,1, lty =2) } if(j > i){#Upper trian xAxis <- (alphaPlants[[j]]$FRic) yAxis <- (alphaPlants[[i]]$FRic) plot(0, ylab="", axes=F, xlab="", type ="n", pch=21, bg="grey80", col="grey30", xlim=c(0, 1), ylim=c(0, 1)) box(lwd =0.8) corr <-cor.test(xAxis, yAxis) leg<-vector("expression",2) leg[1]<- substitute(expression(rho== MYVALUE), list( MYVALUE = max(0, round(corr$estimate,2))))[2] pValor<-round(corr$p.value,3) leg[2] <- ifelse(pValor < 0.001, "p < 0.001", paste0("p = ", pValor)) legend(0.4, 0.5, leg=leg, bty="n", cex = cexLeg, xjust = 0.5, yjust = 0.5) } if(j==1){axis(2, las=1, tcl=0.25,lwd=0.8, at =c (0, 20, 40), labels = c(0, 20, 40))} if(i==6){axis(1,tcl=0.25,lwd=0.8, at =c (0, 20, 40), labels = c(0, 20, 40))} } } } ##### RAO ---- for(i in 1:length(alphaPlants)){#Rows for(j in 1:length(alphaPlants)){#Columns if(i == j){#Diag xAxis <- alphaPlants[[i]]$Rao dens <- density(xAxis) plot(dens$x, dens$y, type="l", ylab="", axes=F, xlab="", xlim=rangeRao, ylim = 1.5 * range(dens$y), col=Colores[i], lwd=lwdLines) points(xAxis, rep(0, length(xAxis)), bg=ColTransp[i], col=Colores[i], pch=21, cex=0.8) box(lwd =0.8) leg <- ifelse(i != 5, paste0(names(alphaPlants)[i]), "M5a") legend("topright", leg=leg, bty="n", cex = cexMethod) if(i==6){axis(1,tcl=0.25,lwd=0.8)} if(i == 1) mtext(side=3, adj=0, "b) Rao", outer=F, line=0.5, cex=1.5) } else{ if(j < i){#Lower trian xAxis <- (alphaPlants[[j]]$Rao) yAxis <- (alphaPlants[[i]]$Rao) densX <- density(xAxis) densY <- density(yAxis) plot(xAxis, yAxis, ylab="", axes=F, xlab="", pch=21, bg="grey80", col="grey30", xlim=rangeRao, ylim=rangeRao) box(lwd =0.8) # abline(lm(yAxis~xAxis), lty=2) abline(0,1, lty =2) } if(j > i){#Upper trian xAxis <- (alphaPlants[[j]]$Rao) yAxis <- (alphaPlants[[i]]$Rao) plot(0, ylab="", axes=F, xlab="", type ="n", pch=21, bg="grey80", col="grey30", xlim=c(0, 1), ylim=c(0, 1)) box(lwd =0.8) corr <-cor.test(xAxis, yAxis) leg<-vector("expression",2) leg[1]<- substitute(expression(rho== MYVALUE), list( MYVALUE = max(0, round(corr$estimate,2))))[2] pValor<-round(corr$p.value,3) leg[2] <- ifelse(pValor < 0.001, "p < 0.001", paste0("p = ", pValor)) legend(0.4, 0.5, leg=leg, bty="n", cex = cexLeg, xjust = 0.5, yjust = 0.5) } if(j==1){axis(2, las=1, tcl=0.25,lwd=0.8, at =c(0, 0.3, 0.6), labels = c(0, 0.3, 0.6))} if(i==6){axis(1,tcl=0.25,lwd=0.8, at =c (0, 0.3, 0.6), labels = c(0, 0.3, 0.6))} } } } #### #### dev.off() #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### ANTS #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### #### invasion <- envir$invasionAnts invasion <- ifelse(invasion == "intact", 0, 1) gcover <- envir$gcover gcoverCenter <- gcover gcoverCenter2 <- gcoverCenter**2 gcoverCenter3 <- gcoverCenter**3 pchUse <- c(16, 21)[invasion + 1] library(RColorBrewer) library(mgcv) library(vegan) Colores<-brewer.pal(8, "Set1") Colores <- Colores[c(4, 5, 8, 7)] alphaPoints<-0.25 ColTransp<-c(rgb(t(col2rgb(Colores)/255), alpha=alphaPoints, maxColorValue=1)) varPlot <- c("FRic" #,"FEve" #, "FDiv" , "Rao") nombresVar <- c("FRic" # , "Functional evenness" # , "Functional divergence" , "Rao") namesMethods <- c("M4. TPD; individuals\nacross populations", "M5a. TPD; species mean\n(sd for each species)", "M5b. TPD; species mean\n(sd common across species)", "M6. Classic;\nspecies mean") lwdLines <- 8 ###### Patterns along the gradient ###### pdf("Manuscript/Figures/FiguresR1/Figure_ants_invasion.pdf", height = 14*(length(varPlot)), width = 14*length(namesMethods), pointsize=55) par(mfrow=c(length(varPlot), length(namesMethods)), mgp=c(1, 0.1, 0), mar = c(3, 3, 2, 1), oma=c(2,2,3,2), mgp = c(2.8, 0.1, 0)) plotNum <- 1 ###### ###### for(i in 1:(length(varPlot))){ #For each component of diversity, one row for(j in 1:length(alphaAnts)){ # For each method, one column labX <- ifelse(i == length(varPlot), "Ground cover dissimilarity", "Ground cover (%)") labY <- ifelse(j == 1, nombresVar[i], "") if(i != 3){ # FRichness and Rao alphaUse <- (alphaAnts[[j]][, i]) allPreds <- c("gcoverCenter", "gcoverCenter2", "invasion", "gcoverCenter:invasion", "gcoverCenter2:invasion") allPredsPlot <- c("ground cover", "ground cover2", "invasion", "ground cover:invasion", "ground cover2:invasion") modelo<-lm(alphaUse ~ (gcoverCenter + gcoverCenter2 ) * invasion, na.action="na.fail") modeloDredge<-dredge(modelo, rank = "AICc", subset = dc(gcoverCenter, gcoverCenter2, gcoverCenter3)) selectedModel <- get.models(modeloDredge, 1)[[1]] if(j == 1){ selectedVarsBestModel <- rownames(anova(selectedModel)) selectedVarsBestModel <- selectedVarsBestModel[ - length(selectedVarsBestModel)] if(length(selectedVarsBestModel) > 0){ nonSelectedVarsBestModel <- allPreds[- which(allPreds %in% selectedVarsBestModel)] }else{ nonSelectedVarsBestModel <- allPreds } selectedVarsModel <- rownames(anova(selectedModel)) selectedVarsModel <- selectedVarsModel[ - length(selectedVarsModel)] positionBestModel <- which(rowSums(!is.na(modeloDredge[, selectedVarsBestModel])) == length(selectedVarsBestModel) & rowSums(!is.na(modeloDredge[, nonSelectedVarsBestModel])) == 0) deltaBestModel <- modeloDredge[positionBestModel, "delta"] } else{ selectedVarsModel <- rownames(anova(selectedModel)) selectedVarsModel <- selectedVarsModel[ - length(selectedVarsModel)] if(length(selectedVarsBestModel) > 0){ positionBestModel <- which(rowSums(!is.na(modeloDredge[, selectedVarsBestModel])) == length(selectedVarsBestModel) & rowSums(!is.na(modeloDredge[, nonSelectedVarsBestModel])) == 0) } else{ positionBestModel <- which(rowSums(!is.na(modeloDredge[, nonSelectedVarsBestModel])) == 0) } deltaBestModel <- modeloDredge[positionBestModel, "delta"] } if(length(selectedVarsModel) == 0){ leg <- vector("expression", 2 + 1) # One for R2, one for no predictors, one for deltaAIC leg[1]<- substitute(expression(R^2== MYVALUE), list( MYVALUE = max(0, round(summary(selectedModel)$r.sq,3))))[2] leg[2]<- "no predictors retained" leg[3]<- substitute(expression(Delta[AICc]== MYVALUE), list( MYVALUE = round(deltaBestModel, 3)))[2] } else{ leg <- vector("expression", 2 + length(selectedVarsModel)) leg[1]<- substitute(expression(R^2== MYVALUE), list( MYVALUE = max(0, round(summary(selectedModel)$r.sq,3))))[2] leg[2]<- substitute(expression(Delta[AICc]== MYVALUE), list( MYVALUE = round(deltaBestModel, 3)))[2] } ### NOW, PLOT PREDICTIONS FOR EACH INVASION TYPE limY <- range(alphaUse); limY[2] <- limY[2] + 0.33 * (limY[2] - limY[1]) plot(alphaUse ~ gcover, ylab="", axes=F, xlab="", pch = pchUse, col=Colores[j], bg = "grey85" , ylim= limY, cex=1.5) mtext(side=1, adj=0.5, labX, outer=F, line=1.5) mtext(side=2, adj=0.5, labY, outer=F, line=1.8) axis(1, tcl=0.25, lwd=0.8) axis(2, las=1, tcl=0.25,lwd=0.8) box( bty = "l", lwd =0.8 ) posLeg <- "topleft" legend(posLeg, leg, adj=0, bty="n", cex=1.25) if(i==1){ mtext(side=3, adj=0.5, paste0(namesMethods[j]), outer=F, line=1.5) } mtext(side=3, adj=0, paste0(letters[plotNum], ") "), outer=F, line=0) plotNum <- plotNum + 1 selectedModel gcoverGrad <- seq(min(gcoverCenter), max(gcoverCenter), length.out = 100) gcoverGrad2 <- gcoverGrad**2 gcoverGrad3 <- gcoverGrad**3 predMatInv0 <- data.frame(gcoverCenter = gcoverGrad, gcoverCenter2 = gcoverGrad2, gcoverCenter3 = gcoverGrad3, invasion = rep(0, length(gcoverGrad))) linesInv0 <- predict(selectedModel, newdata = predMatInv0) predMatInv1 <- data.frame(gcoverCenter = gcoverGrad, gcoverCenter2 = gcoverGrad2, gcoverCenter3 = gcoverGrad3, invasion = rep(1, length(gcoverGrad))) linesInv1 <- predict(selectedModel, newdata = predMatInv1) lines(gcoverGrad, linesInv0, lwd=lwdLines, col = Colores[j], lty= 1) lines(gcoverGrad, linesInv1, lwd=lwdLines, col = Colores[j], lty= 2) } } } ###### ###### dev.off() ################## CORRELATIONS BETWEEN METHODS ANTS cexLeg <- 1.3 cexMethod <- 1.4 varPlot <- c("FRic" #,"FEve" #, "FDiv" , "Rao") nombresVar <- c("Functional richness" # , "Functional evenness" # , "Functional divergence" , "Rao") namesMethods <- c("M4. TPD; individuals\nacross populations", "M5. TPD; species mean\n(sd for each species)", "M5b. TPD; species mean\n(sd common across species)", "M6. Classic;\nspecies mean") lwdLines <- 8 pdf("Manuscript/Figures/FiguresR1/Figure_ants_correlationsAlphaFRic.pdf", width = 20, height = 10, pointsize=23) #### #### matFRicPlot <- matrix(1:16, nrow=4, ncol=4, byrow=T) matRaoPlot <- matrix((1:16)+16, nrow=4, ncol=4, byrow=T) matBlank <- matrix(0, nrow=4, ncol=1) layoutMat <- cbind(matFRicPlot, matBlank, matRaoPlot) widthsCols <- c(rep(1, 4), 0.5, rep(1, 4)) widthsCols <- widthsCols/sum(widthsCols) layout(layoutMat, widths = widthsCols) par(mgp=c(1, 0.1, 0), mar = c(0.1, 0.1, 0.1, 0.1), oma=c(1, 1, 3 , 1)) ##### FUNCTIONAL RICHNESS ---- rangeFRic <- rangeRao <- numeric() for(i in 1:length(alphaAnts)){ rangeFRic <- c(rangeFRic, alphaAnts[[i]]$FRic) rangeRao <- c(rangeRao, alphaAnts[[i]]$Rao) } rangeFRic <- c(0, 1.1 * max(rangeFRic)) rangeRao <- c(0, 1.1 * max(rangeRao)) for(i in 1:length(alphaAnts)){#Rows for(j in 1:length(alphaAnts)){#Columns if(i == j){#Diag xAxis <- alphaAnts[[i]]$FRic dens <- density(xAxis) plot(dens$x, dens$y, type="l", ylab="", axes=F, xlab="", xlim=rangeFRic, ylim = 1.5 * range(dens$y), col=Colores[i], lwd=lwdLines) points(xAxis, rep(0, length(xAxis)), bg=ColTransp[i], col=Colores[i], pch=21, cex=0.8) box(lwd =0.8) legend("topright", leg=paste0(names(alphaAnts)[i]), bty="n", cex = cexMethod) if(i==4){axis(1,tcl=0.25,lwd=0.8)} if(i == 1) mtext(side=3, adj=0, "a) FRic", outer=F, line=0.5, cex=1.5) } else{ if(j < i){#Lower trian xAxis <- (alphaAnts[[j]]$FRic) yAxis <- (alphaAnts[[i]]$FRic) densX <- density(xAxis) densY <- density(yAxis) plot(xAxis, yAxis, ylab="", axes=F, xlab="", pch=21, bg="grey80", col="grey30", xlim=rangeFRic, ylim=rangeFRic) box(lwd =0.8) # abline(lm(yAxis~xAxis), lty=2) abline(0,1, lty =2) } if(j > i){#Upper trian xAxis <- (alphaAnts[[j]]$FRic) yAxis <- (alphaAnts[[i]]$FRic) plot(0, ylab="", axes=F, xlab="", type ="n", pch=21, bg="grey80", col="grey30", xlim=c(0,1), ylim=c(0,1)) box(lwd =0.8) corr <-cor.test(xAxis, yAxis) leg<-vector("expression",2) leg[1]<- substitute(expression(rho== MYVALUE), list( MYVALUE = max(0, round(corr$estimate,2))))[2] pValor<-round(corr$p.value,3) leg[2] <- ifelse(pValor < 0.001, "p < 0.001", paste0("p = ", pValor)) legend(0.4, 0.5, leg=leg, bty="n", cex = cexLeg, xjust = 0.5, yjust = 0.5) } if(j==1){axis(2, las=1, tcl=0.25,lwd=0.8, at =c (0, 40, 80), labels = c(0, 40, 80))} if(i==4){axis(1,tcl=0.25,lwd=0.8, at =c (0, 40, 80), labels = c(0, 40, 80))} } } } ##### RAO ---- for(i in 1:length(alphaAnts)){#Rows for(j in 1:length(alphaAnts)){#Columns if(i == j){#Diag xAxis <- alphaAnts[[i]]$Rao dens <- density(xAxis) plot(dens$x, dens$y, type="l", ylab="", axes=F, xlab="", xlim=rangeRao, ylim = 1.5 * range(dens$y), col=Colores[i], lwd=lwdLines) points(xAxis, rep(0, length(xAxis)), bg=ColTransp[i], col=Colores[i], pch=21, cex=0.8) box(lwd =0.8) legend("topright", leg=paste0(names(alphaAnts)[i]), bty="n", cex = cexMethod) if(i==4){axis(1,tcl=0.25,lwd=0.8)} if(i == 1) mtext(side=3, adj=0, "b) Rao", outer=F, line=0.5, cex=1.5) } else{ if(j < i){#Lower trian xAxis <- (alphaAnts[[j]]$Rao) yAxis <- (alphaAnts[[i]]$Rao) densX <- density(xAxis) densY <- density(yAxis) plot(xAxis, yAxis, ylab="", axes=F, xlab="", pch=21, bg="grey80", col="grey30", xlim=rangeRao, ylim=rangeRao) box(lwd =0.8) # abline(lm(yAxis~xAxis), lty=2) abline(0,1, lty =2) } if(j > i){#Upper trian xAxis <- (alphaAnts[[j]]$Rao) yAxis <- (alphaAnts[[i]]$Rao) plot(0, ylab="", axes=F, xlab="", type ="n", pch=21, bg="grey80", col="grey30", xlim=c(0,1), ylim=c(0,1)) box(lwd =0.8) corr <-cor.test(xAxis, yAxis) leg<-vector("expression",2) leg[1]<- substitute(expression(rho== MYVALUE), list( MYVALUE = max(0, round(corr$estimate,2))))[2] pValor<-round(corr$p.value,3) leg[2] <- ifelse(pValor < 0.001, "p < 0.001", paste0("p = ", pValor)) legend(0.4, 0.5, leg=leg, bty="n", cex = cexLeg, xjust = 0.5, yjust = 0.5) } if(j==1){axis(2, las=1, tcl=0.25,lwd=0.8, at =c(0, 0.4, 0.8), labels = c(0, 0.4, 0.8))} if(i==4){axis(1,tcl=0.25,lwd=0.8, at =c (0, 0.4, 0.8), labels = c(0, 0.4, 0.8))} } } } #### #### dev.off()