#Step 5 - Conduct statistical analyses and make figures #A lot happens in this script. This is where all the analyses are done and figures are made. library(lme4) library(ggplot2) library(lmerTest) library(pbkrtest) library(raster) library(dismo) library(rgdal) library(rstanarm) library(tibble) library(shinystan) library(dplyr) library(grid) setwd("") #Used to evaluate f.normalityplots <- function(a) { par(mfrow = c(2, 2)) plot(a) boxplot(a) hist(a, main = "") qqnorm(a) qqline(a, lty = 2) shapiro.test(a) } #Used to plot multiple ggplot objects together f.multiplot <- function(l.plot, s.cols) { require(grid) s.numplots = length(l.plot) layout <- matrix( seq(1, s.cols * ceiling(s.numplots / s.cols)), ncol = s.cols, nrow = ceiling(s.numplots / s.cols) ) grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) for (i in 1:s.numplots) { matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print( l.plot[[i]], vp = viewport( layout.pos.row = matchidx$row, layout.pos.col = matchidx$col ) ) } } #The region-specific projections from Step 4 df.me.pred <- read.table("040-BigRun-NewFinalMaxentProjections2.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #MAXENT projections in climate space (not used) df.me.cpred <- read.table("040-BigRun-NewFinalMaxentClimProjections.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #Results of native-to-invasive ordination models from previous study df.ord <- read.table("025-RESULTS-Stats-NATUREEE-FINAL-NEE3.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #Native-range ordinations from previous study df.ord.n <- read.table("025-RESULTS-Stats-N-NATUREEE-FINAL-NEE3.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #MAXENT model evaluation statistics df.evals <- read.table("025-MAXENT-Evals.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #Center-of-gravity comparisons df.cogcomp <- read.table("040-BigRun-COGChangesPerComp.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #Species data (by sp.id) df.key <- read.table("025-PointData-KeyMaster-FINAL.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #The global projections from Step 4 df.me.pred.w <- read.table("040-BigRun-NewFinalMaxentProjections-W2.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #Global ordinations from previous study df.ord.w <- read.table("040-RESULTS-Stats-NATUREEE-W2-FINAL-NEE3.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #These lines make matrices that contain information about each region and how each comparison is named. #These help to merge the dataframes above and provided easily-retrieved data about what each region-code means. mx.use <- rbind( t(apply( matrix(c(2:6, 1, 3:6, 1:2, 4:6, 1:3, 5:6, 1:4, 6, 1:5, rep(1:6, each = 5)), ncol = 2), 1, function(i) c(paste("n", i[1], "-n", i[2], sep = ""), paste("n", i[2], "-n", i[2], sep = "")) )), t(apply( expand.grid(1:6, 1:6), 1, function(i) c(paste("n", i[1], "-i", i[2], sep = ""), paste("i", i[2], "-i", i[2], sep = "")) )) ) mx.use <- cbind( mx.use, c(2:6, 1, 3:6, 1:2, 4:6, 1:3, 5:6, 1:4, 6, 1:5, rep(1:6, times = 6)), c(rep(1:6, each = 5), rep(1:6, each = 6)) ) mx.combos <- as.matrix(expand.grid(1:6, 1:6)) mx.combos <- cbind(mx.combos, 1:36) mx.samecombos <- mx.combos[mx.combos[,1] != mx.combos[,2],] mx.use <- cbind(mx.use, c(rep("n", 30), rep("i", 36)), c(mx.samecombos[,3], mx.combos[,3]), 1:nrow(mx.use)) #The following block of commands clean up and merge the dataframes loaded above #Many data are transformed, scaled & centered #Certain functional groups are removed #As are records from region 4 (IP) due to problems with the ordinations, per. the previous study df.ord <- df.ord[sapply(1:nrow(df.ord), function(i) sum(is.na(df.ord[i,])) < (ncol(df.ord) - 2)),] df.ord.n <- df.ord.n[sapply(1:nrow(df.ord.n), function(i) sum(is.na(df.ord.n[i,])) < (ncol(df.ord.n) - 2)),] df.m.geog <- df.me.pred df.m.geog$id <- paste(df.m.geog$sp.id, "-", df.m.geog$key, sep = "") df.m.clim <- df.me.cpred df.m.clim$id <- paste(df.m.clim$sp.id, "-", df.m.clim$key, sep = "") df.m.clim <- subset(df.m.clim, select = -c(p1, p2, key, comptype, sp.id)) names(df.m.clim)[-which(names(df.m.clim) == "id")] <- paste("c.", names(df.m.clim)[-which(names(df.m.clim) == "id")], sep = "") df.o.clim.1 <- df.ord.n[, c("sp.id", "reg", "u.carea", "s.carea", "e.carea", "d", "a", "u", "s", "e", "p", "cent.x", "cent.y", "a.ratio", "u.ratio", "s.ratio", "e.ratio", "p.ratio", "dcc", "cor.dcw", "dcc.exp", "dcc.unf", "dc", "dc.e", "dc.u")] df.o.clim.1 <- df.o.clim.1[df.o.clim.1$reg %in% mx.use[1:30, 6],] df.o.clim.1$reguse <- sapply(df.o.clim.1$reg, function(i) which(mx.use[1:30, 6] == i)) df.o.clim.2 <- df.ord[, c("sp.id", "reg","d", "u.carea", "s.carea", "e.carea", "a", "u", "s", "e", "p", "cent.x", "cent.y", "a.ratio", "u.ratio", "s.ratio", "e.ratio", "p.ratio", "dcc", "cor.dcw", "dcc.exp", "dcc.unf", "dc", "dc.e", "dc.u")] df.o.clim.2$reguse <- df.o.clim.2$reg + 30 df.o.clim <- rbind(df.o.clim.1, df.o.clim.2) df.o.clim$id <- paste(df.o.clim$sp.id, "-", df.o.clim$reguse, sep = "") df.o.clim <- subset(df.o.clim, select = -c(sp.id, reg, reguse)) df.m.dcog <- df.cogcomp df.m.dcog$id <- paste(df.m.dcog$sp.id, "-", df.m.dcog$comp, sep = "") df.m.dcog <- subset(df.m.dcog, select = -c(sp.id, comp)) df.merge <- merge(df.m.geog, df.m.clim, by = "id", all = T) df.merge <- merge(df.merge, df.o.clim, by = "id", all = T) df.merge <- merge(df.merge, df.m.dcog, by = "id", all = T) df.merge$sp.id <- sapply(df.merge$id, function(i) substr(i, 1, gregexpr("-", i)[[1]][1] - 1)) df.merge$type <- as.numeric(sapply(df.merge$id, function(i) substr(i, gregexpr("-", i)[[1]][1] + 1, nchar(i)))) df.merge <- subset(df.merge, select = -c(key, p1, p2)) df.merge <- merge(df.merge, df.key[, c("sp.id", "cultivated", "growth", "life")], by = "sp.id", all.x = T) df.merge$cult <- df.merge$cultivated df.merge$cult[df.merge$cult == "widely"] <- "yes" df.merge$cult <- factor(df.merge$cult) df.merge$comptype <- sapply(df.merge$type, function(i) ifelse(i <= 30, "NN", "NI")) df.merge$reg <- as.numeric(sapply(df.merge$type, function(i) mx.use[i, 6])) df.w <- df.me.pred.w df.w$id <- paste(df.w$sp.id, "-", df.w$key, sep = "") df.w <- merge(df.w, df.ord.w, by = "sp.id", all = T) df.w <- merge(df.w, df.key[, c("sp.id", "cultivated", "growth", "life")], by = "sp.id", all.x = T) df.w$cult <- df.w$cultivated df.w$cult[df.w$cult == "widely"] <- "yes" df.w$cult <- factor(df.w$cult) df.merge.ni <- df.merge[df.merge$comptype == "NI",] df.merge.nn <- df.merge[df.merge$comptype == "NN",] df.use <- subset(df.merge, life != "Aquatic" & life != "NA" & life != "Parasitic" & growth != "Aquatic" & growth != "NA" & growth != "Parasitic" & reg %in% c(1:3, 5:9, 11:15, 17:18, 25:27, 29:33, 35:36)) df.use$life <- factor(df.use$life) df.use$growth <- factor(df.use$growth) df.use$reg <- factor(df.use$reg) df.use$type <- factor(df.use$type) df.use$sp.id <- factor(df.use$sp.id) df.use$wood <- factor(sapply(as.character(df.use$growth), function(i) if(i == "Shrub" | i == "Tree" | i == "Vine" | i == "Subshrub") "Woody" else i)) df.use$life2 <- df.use$life df.use[df.use$life2 == "Biennial", "life2"] <- "Perennial" df.use$life2 <- factor(df.use$life2) df.use.w <- subset(df.w, life != "Aquatic" & life != "NA" & life != "Parasitic" & growth != "Aquatic" & growth != "NA" & growth != "Parasitic") df.use.w$life <- factor(df.use.w$life) df.use.w$growth <- factor(df.use.w$growth) df.use.w$reg <- factor(df.use.w$reg) df.use.w$sp.id <- factor(df.use.w$sp.id) df.use.w$wood <- factor(sapply(as.character(df.use.w$growth), function(i) if(i == "Shrub" | i == "Tree" | i == "Vine" | i == "Subshrub") "Woody" else i)) df.use.w$life2 <- df.use.w$life df.use.w[df.use.w$life2 == "Biennial", "life2"] <- "Perennial" df.use.w$life2 <- factor(df.use.w$life2) df.use.ni <- df.use[df.use$comptype == "NI",] df.use.nn <- df.use[df.use$comptype == "NN",] df.use$z.bsd <- scale(df.use$bias.stat.d) df.use$z.bsi <- scale(df.use$bias.stat.i) df.use$z.rangediff <- scale(df.use$bias.rangediff) df.use$cz.bias.u <- scale(df.use$bias.u50^(1/3)) df.use$cz.bias.s <- scale(df.use$bias.s50^(1/3)) df.use$cz.bias.e <- scale(df.use$bias.e50^(1/3)) df.use$cz.bias.n <- scale(df.use$bias.n.50a^(1/3)) df.use$cz.bias.nn <- scale(df.use$bias.n.50b^(1/3)) df.use$z.d <- scale(df.use$d) df.use$z.cent.x <- scale(df.use$cent.x) df.use$cz.u <- scale(df.use$u^(1/3)) df.use$cz.s <- scale(df.use$s^(1/3)) df.use$cz.e <- scale(df.use$e^(1/3)) df.use$cz.uc <- scale(df.use$u.carea^(1/3)) df.use$cz.sc <- scale(df.use$s.carea^(1/3)) df.use$cz.ec <- scale(df.use$e.carea^(1/3)) df.use$bias.sr <- df.use$bias.s50 / (df.use$bias.s50 + df.use$bias.u50 + df.use$bias.e50) df.use$bias.ur <- df.use$bias.u50 / (df.use$bias.s50 + df.use$bias.u50 + df.use$bias.e50) df.use$bias.er <- df.use$bias.e50 / (df.use$bias.s50 + df.use$bias.u50 + df.use$bias.e50) df.use$cz.sr <- scale(asin(sqrt(df.use$bias.sr))) df.use$cz.ur <- scale(asin(sqrt(df.use$bias.ur))) df.use$cz.er <- scale(asin(sqrt(df.use$bias.er))) df.use.w$z.bsd <- scale(df.use.w$bias.stat.d) df.use.w$z.rangediff <- scale(df.use.w$bias.rangediff) df.use.w$cz.bias.u <- scale(df.use.w$bias.u50^(1/3)) df.use.w$cz.bias.s <- scale(df.use.w$bias.s50^(1/3)) df.use.w$cz.bias.e <- scale(df.use.w$bias.e50^(1/3)) df.use.w$cz.bias.n <- scale(df.use.w$bias.n.50a^(1/3)) df.use.w$z.d <- scale(df.use.w$d) df.use.w$z.cent.x <- scale(df.use.w$cent.x) df.use.w$cz.u <- scale(df.use.w$u^(1/3)) df.use.w$cz.s <- scale(df.use.w$s^(1/3)) df.use.w$cz.e <- scale(df.use.w$e^(1/3)) df.use.w$cz.uc <- scale(df.use.w$u.carea^(1/3)) df.use.w$cz.sc <- scale(df.use.w$s.carea^(1/3)) df.use.w$cz.ec <- scale(df.use.w$e.carea^(1/3)) #This is the working dataframe used for most analyses. #It is saved to it can be retrieved without re-running the steps above. write.table(df.use, "040-DFUseNew.csv", sep = ",", na = "", row.names = F) #Read df.use here to avoid running lines 107-207 df.use <- read.table("040-DFUseNew.csv", header = T, sep = ",", na.strings = "", dec = ".", strip.white = T) #This function runs Baysian models in which life-form is a predictor, along with comptype (NI vs NN). #Arguments include: # ch.y - the name of the y-variable # ch.x - the name of the continuous predictor to be used #The function will be iterated across multiple preditor variables. #Results are stored in a list: # [[1]] - ch.y (the name of the response variable) # [[2]] - cy.x (the name of the specified predictor) # [[3]] - m.do (the fitted model object) # [[4]] - b (the posterior predictions of 4000 models) # [[5]] - df.b (posterior marginal means) # [[6]] - df.seg (quantiles of the posterior marginal means, e.g. 95% confience interval) # [[7]] - df.bx (posterior marginal effects of the specified predictor) # [[8]] - df.segx (quantiles of the posterior marginal effects) f.m.life2 <- function(ch.y, ch.x, df.do) { m.do <- stan_lmer(data = df.do, get(ch.y) ~ get(ch.x) * comptype * life2 + (1|reg) + (1|sp.id), iter = 4000) b <- as_tibble(m.do) %>% rename(intercept = `(Intercept)`) %>% select(-sigma) df.b <- data.frame( ii.an = with(b, intercept), ii.pe = with(b, intercept + life2Perennial), ii.wo = with(b, intercept + life2Woody), ni.an = with(b, intercept + comptypeNN), ni.pe = with(b, intercept + comptypeNN + life2Perennial + `comptypeNN:life2Perennial`), ni.wo = with(b, intercept + comptypeNN + life2Woody + `comptypeNN:life2Woody`) ) df.bx <- data.frame( ii.an = with(b, `get(ch.x)`), ii.pe = with(b, `get(ch.x)` + `get(ch.x):life2Perennial`), ii.wo = with(b, `get(ch.x)` + `get(ch.x):life2Woody`), ni.an = with(b, `get(ch.x)` + `get(ch.x):comptypeNN`), ni.pe = with(b, `get(ch.x)` + `get(ch.x):comptypeNN` + `get(ch.x):life2Perennial` + `get(ch.x):comptypeNN:life2Perennial`), ni.wo = with(b, `get(ch.x)` + `get(ch.x):comptypeNN` + `get(ch.x):life2Woody` + `get(ch.x):comptypeNN:life2Woody`) ) df.seg <- do.call(rbind, lapply(1:ncol(df.b), function(i) data.frame( s.i = i, st.var = names(df.b)[i], st.type = substr(names(df.b)[i], 1, 2), s02.5 = quantile(df.b[, i], 0.025), s25.0 = quantile(df.b[, i], 0.250), s75.0 = quantile(df.b[, i], 0.750), s87.5 = quantile(df.b[, i], 0.875), s97.5 = quantile(df.b[, i], 0.975) ))) df.seg$pdif <- c(sum(df.b$ii.an > df.b$ni.an), sum(df.b$ii.pe > df.b$ni.pe), sum(df.b$ii.wo > df.b$ni.wo), rep(NA, 3)) / nrow(df.b) df.segx <- do.call(rbind, lapply(1:ncol(df.bx), function(i) data.frame( s.i = i, st.var = names(df.bx)[i], st.type = substr(names(df.bx)[i], 1, 2), s02.5 = quantile(df.bx[, i], 0.025), s25.0 = quantile(df.bx[, i], 0.250), s75.0 = quantile(df.bx[, i], 0.750), s87.5 = quantile(df.bx[, i], 0.875), s97.5 = quantile(df.bx[, i], 0.975) ))) df.segx$pdif <- c(sum(df.bx$ii.an > df.bx$ni.an), sum(df.bx$ii.pe > df.bx$ni.pe), sum(df.bx$ii.wo > df.bx$ni.wo), rep(NA, 3)) / nrow(df.bx) list(ch.y, ch.x, m.do, b, df.b, df.seg, df.bx, df.segx) } #This function runs Baysian models in which cultivation status is a predictor, along with comptype (NI vs NN). #Arguments include: # ch.y - the name of the y-variable # ch.x - the name of the continuous predictor to be used #The function will be iterated across multiple preditor variables. #Results are stored in a list: # [[1]] - ch.y (the name of the response variable) # [[2]] - cy.x (the name of the specified predictor) # [[3]] - m.do (the fitted model object) # [[4]] - b (the posterior predictions of 4000 models) # [[5]] - df.b (posterior marginal means) # [[6]] - df.seg (quantiles of the posterior marginal means, e.g. 95% confience interval) # [[7]] - df.bx (posterior marginal effects of the specified predictor) # [[8]] - df.segx (quantiles of the posterior marginal effects) f.m.cult <- function(ch.y, ch.x, df.do) { m.do <- stan_lmer(data = df.do, get(ch.y) ~ get(ch.x) * comptype * cultivated + (1|reg) + (1|sp.id), iter = 4000) b <- as_tibble(m.do) %>% rename(intercept = `(Intercept)`) %>% select(-sigma) df.b <- data.frame( ii.no = with(b, intercept), ii.ye = with(b, intercept + cultivatedyes), ii.wi = with(b, intercept + cultivatedwidely), ni.no = with(b, intercept + comptypeNN), ni.ye = with(b, intercept + comptypeNN + cultivatedyes + `comptypeNN:cultivatedyes`), ni.wi = with(b, intercept + comptypeNN + cultivatedwidely + `comptypeNN:cultivatedwidely`) ) df.bx <- data.frame( ii.no = with(b, `get(ch.x)`), ii.ye = with(b, `get(ch.x)` + `get(ch.x):cultivatedyes`), ii.wi = with(b, `get(ch.x)` + `get(ch.x):cultivatedwidely`), ni.no = with(b, `get(ch.x)` + `get(ch.x):comptypeNN`), ni.ye = with(b, `get(ch.x)` + `get(ch.x):comptypeNN` + `get(ch.x):cultivatedyes` + `get(ch.x):comptypeNN:cultivatedyes`), ni.wi = with(b, `get(ch.x)` + `get(ch.x):comptypeNN` + `get(ch.x):cultivatedwidely` + `get(ch.x):comptypeNN:cultivatedwidely`) ) df.seg <- do.call(rbind, lapply(1:ncol(df.b), function(i) data.frame( s.i = i, st.var = names(df.b)[i], st.type = substr(names(df.b)[i], 1, 2), s02.5 = quantile(df.b[, i], 0.025), s25.0 = quantile(df.b[, i], 0.250), s75.0 = quantile(df.b[, i], 0.750), s87.5 = quantile(df.b[, i], 0.875), s97.5 = quantile(df.b[, i], 0.975) ))) df.seg$pdif <- c(sum(df.b$ii.no > df.b$ni.no), sum(df.b$ii.ye > df.b$ni.ye), sum(df.b$ii.wi > df.b$ni.wi), rep(NA, 3)) / nrow(df.b) df.segx <- do.call(rbind, lapply(1:ncol(df.bx), function(i) data.frame( s.i = i, st.var = names(df.bx)[i], st.type = substr(names(df.bx)[i], 1, 2), s02.5 = quantile(df.bx[, i], 0.025), s25.0 = quantile(df.bx[, i], 0.250), s75.0 = quantile(df.bx[, i], 0.750), s87.5 = quantile(df.bx[, i], 0.875), s97.5 = quantile(df.bx[, i], 0.975) ))) df.segx$pdif <- c(sum(df.bx$ii.no > df.bx$ni.no), sum(df.bx$ii.ye > df.bx$ni.ye), sum(df.bx$ii.wi > df.bx$ni.wi), rep(NA, 3)) / nrow(df.bx) list(ch.y, ch.x, m.do, b, df.b, df.seg, df.bx, df.segx) } #Specify the response variables for Bayesian models #z.bsd (centered and scaled geographic overlap D) #z.rangediff (centered and scaled shift in habiatat suitability (deltaHS)) #cz.bias.s (centered and scaled range stability) #dz.bias.u (centered and scaled range unfilling) #dz.bias.e (centered and scaled range expansion) #cz.bias.n (centered and scaled range size (deltaRS) ch.m.use <- c("z.bsd", "z.rangediff", "cz.bias.s", "cz.bias.u", "cz.bias.e", "cz.bias.n") #Models take a long time to run, these are run and saved so they can be loaded directly into R in future sessions. #There is no reason to run the following, as they are directly loaded below #l.m.d <- lapply(ch.m.use, f.m.life2, ch.x = "z.d", df.do = df.use) #save(l.m.d, file = "040-ModList-OrdD-Effects.Rdata") #l.m.x <- lapply(ch.m.use, f.m.life2, ch.x = "z.cent.x", df.do = df.use) #save(l.m.x, file = "040-ModList-OrdCentX-Effects.Rdata") #l.m.u <- lapply(ch.m.use, f.m.life2, ch.x = "cz.u", df.do = df.use) #save(l.m.u, file = "040-ModList-OrdU-Effects.Rdata") #l.m.s <- lapply(ch.m.use, f.m.life2, ch.x = "cz.s", df.do = df.use) #save(l.m.s, file = "040-ModList-OrdS-Effects.Rdata") #l.m.e <- lapply(ch.m.use, f.m.life2, ch.x = "cz.e", df.do = df.use) #save(l.m.e, file = "040-ModList-OrdE-Effects.Rdata") #l.m.uc <- lapply(ch.m.use, f.m.life2, ch.x = "cz.uc", df.do = df.use) #save(l.m.uc, file = "040-ModList-OrdUc-Effects.Rdata") #l.m.sc <- lapply(ch.m.use, f.m.life2, ch.x = "cz.sc", df.do = df.use) #save(l.m.sc, file = "040-ModList-OrdSc-Effects.Rdata") #l.m.ec <- lapply(ch.m.use, f.m.life2, ch.x = "cz.ec", df.do = df.use) #save(l.m.ec, file = "040-ModList-OrdEc-Effects.Rdata") #l.m.nn <- lapply(c("z.d", "cz.ur", "cz.sr", "cz.er", "cz.bias.n"), f.m.life2, ch.x = "cz.bias.nn", df.do = df.use) #save(l.m.nn, file = "040-ModList-NatRangeSize-Effects.Rdata") #l.m.d.cult <- lapply(ch.m.use, f.m.cult, ch.x = "z.d", df.do = df.use) #save(l.m.d.cult, file = "040-ModList-OrdD-Effects-CULT.Rdata") #l.m.x.cult <- lapply(ch.m.use, f.m.cult, ch.x = "z.cent.x", df.do = df.use) #save(l.m.x.cult, file = "040-ModList-OrdCentX-Effects-CULT.Rdata") #l.m.u.cult <- lapply(ch.m.use, f.m.cult, ch.x = "cz.u", df.do = df.use) #save(l.m.u.cult, file = "040-ModList-OrdU-Effects-CULT.Rdata") #l.m.s.cult <- lapply(ch.m.use, f.m.cult, ch.x = "cz.s", df.do = df.use) #save(l.m.s.cult, file = "040-ModList-OrdS-Effects-CULT.Rdata") #l.m.e.cult <- lapply(ch.m.use, f.m.cult, ch.x = "cz.e", df.do = df.use) #save(l.m.e.cult, file = "040-ModList-OrdE-Effects-CULT.Rdata") #l.m.uc.cult <- lapply(ch.m.use, f.m.cult, ch.x = "cz.uc", df.do = df.use) #save(l.m.uc.cult, file = "040-ModList-OrdUc-Effects-CULT.Rdata") #l.m.sc.cult <- lapply(ch.m.use, f.m.cult, ch.x = "cz.sc", df.do = df.use) #save(l.m.sc.cult, file = "040-ModList-OrdSc-Effects-CULT.Rdata") #l.m.ec.cult <- lapply(ch.m.use, f.m.cult, ch.x = "cz.ec", df.do = df.use) #save(l.m.ec.cult, file = "040-ModList-OrdEc-Effects-CULT.Rdata") #l.m.d <- lapply(ch.m.use, f.m.life2, ch.x = "z.d", df.do = df.use) #save(l.m.d, file = "040-ModList-OrdD-Effects.Rdata") #Effects of niche overlap (climate D) load("040-ModList-OrdD-Effects.Rdata") load("040-ModList-OrdD-Effects-CULT.Rdata") #Effects of niche position (RPCA1 centroid) load("040-ModList-OrdCentX-Effects.Rdata") load("040-ModList-OrdCentX-Effects-CULT.Rdata") #Effects of niche unfilling (not used) load("040-ModList-OrdU-Effects.Rdata") load("040-ModList-OrdU-Effects-CULT.Rdata") #Effects of niche stability (not used) load("040-ModList-OrdS-Effects.Rdata") load("040-ModList-OrdS-Effects-CULT.Rdata") #Effects of niche expansion (not used) load("040-ModList-OrdE-Effects.Rdata") load("040-ModList-OrdE-Effects-CULT.Rdata") #Effects of niche unfilling (new way - used) load("040-ModList-OrdUc-Effects.Rdata") load("040-ModList-OrdUc-Effects-CULT.Rdata") #Effects of niche stability (new way - used) load("040-ModList-OrdSc-Effects.Rdata") load("040-ModList-OrdSc-Effects-CULT.Rdata") #Effects of niche expansion (new way - used) load("040-ModList-OrdEc-Effects.Rdata") load("040-ModList-OrdEc-Effects-CULT.Rdata") #Global models (less time needed to run these) #Effects of overlap (climate D) m.d.w <- stan_glm(data = df.use.w, z.bsd ~ z.d * life2) #Effects of stability (S) m.sc.w <- stan_glm(data = df.use.w, cz.bias.s ~ cz.sc * life2) #Effects of unfilling (U) m.uc.w <- stan_glm(data = df.use.w, cz.bias.u ~ cz.uc * life2) #Effects of expansion (E) m.ec.w <- stan_glm(data = df.use.w, cz.bias.e ~ cz.ec * life2) #List the global models together, and calculate posterior probabilities, #as above for the regional models #Results are stored in a list: # [[1]] - i (the fitted model object) # [[2]] - b (the posterior predictions of 4000 models) # [[3]] - df.b (posterior marginal means) # [[4]] - df.seg (quantiles of the posterior marginal means, e.g. 95% confience interval) # [[5]] - df.bx (posterior marginal effects of the specified predictor) # [[6]] - df.segx (quantiles of the posterior marginal effects) l.m.w <- lapply(list(m.d.w, m.sc.w, m.uc.w, m.ec.w), function(i) { b <- as_tibble(i) %>% rename(intercept = `(Intercept)`) %>% select(-sigma) names(b)[c(2, 5, 6)] <- c("v.x", "v.x:life2Perennial", "v.x:life2Woody") df.b <- data.frame( ii.an = with(b, intercept), ii.pe = with(b, intercept + life2Perennial), ii.wo = with(b, intercept + life2Woody) ) df.bx <- data.frame( ii.an = with(b, v.x), ii.pe = with(b, v.x + `v.x:life2Perennial`), ii.wo = with(b, v.x + `v.x:life2Woody`) ) df.seg <- do.call(rbind, lapply(1:ncol(df.b), function(i) data.frame( s.i = i, st.var = names(df.b)[i], st.type = substr(names(df.b)[i], 1, 2), s02.5 = quantile(df.b[, i], 0.025), s25.0 = quantile(df.b[, i], 0.250), s75.0 = quantile(df.b[, i], 0.750), s87.5 = quantile(df.b[, i], 0.875), s97.5 = quantile(df.b[, i], 0.975) ))) df.segx <- do.call(rbind, lapply(1:ncol(df.bx), function(i) data.frame( s.i = i, st.var = names(df.bx)[i], st.type = substr(names(df.bx)[i], 1, 2), s02.5 = quantile(df.bx[, i], 0.025), s25.0 = quantile(df.bx[, i], 0.250), s75.0 = quantile(df.bx[, i], 0.750), s87.5 = quantile(df.bx[, i], 0.875), s97.5 = quantile(df.bx[, i], 0.975) ))) list(i, b, df.b, df.seg, df.bx, df.segx) }) #Make a ggplot theme with positionable legend f.theme <- function(v.pos) theme( axis.text = element_text(size = 14, color = "black"), axis.title.y = element_text(size = 16), axis.title.x = element_text(size = 16), axis.line.x = element_line(color = "black"), axis.line.y = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.ticks.y = element_line(color = "black"), legend.position = v.pos, legend.key.width = unit(1, "cm"), legend.key = element_rect(fill = "white"), legend.title = element_blank(), legend.text = element_text(size = 14), panel.grid = element_blank(), panel.background = element_rect(fill = "white") ) ##### REMAINDER OF SCRIPT IS FIGURE PLOTTING ##### ##### No more annotations ##### f.cube <- function(i) i^(1/3) f.plot.slopes <- function(l.use) { df.x <- l.use[[1]][[7]] names(df.x) <- paste(names(df.x), ".x", sep = "") df.x <- sample_n(cbind(l.use[[1]][[5]], df.x), 1000) ggplot(data = l.use[[2]], aes(x = get(l.use[[3]]), y = get(l.use[[4]]), color = comptype)) + geom_point(size = l.use[[5]]) + scale_color_manual(values = l.use[[6]]) + geom_abline(data = df.x, aes(intercept = ni.an, slope = ni.an.x), color = l.use[[7]][4], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ni.pe, slope = ni.pe.x), color = l.use[[7]][5], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ni.wo, slope = ni.wo.x), color = l.use[[7]][6], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ii.an, slope = ii.an.x), color = l.use[[7]][1], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ii.pe, slope = ii.pe.x), color = l.use[[7]][2], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ii.wo, slope = ii.wo.x), color = l.use[[7]][3], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = median(ni.an), slope = median(ni.an.x)), color = l.use[[8]][4], size = 1, linetype = l.use[[17]][1]) + geom_abline(data = df.x, aes(intercept = median(ni.pe), slope = median(ni.pe.x)), color = l.use[[8]][5], size = 1, linetype = l.use[[17]][2]) + geom_abline(data = df.x, aes(intercept = median(ni.wo), slope = median(ni.wo.x)), color = l.use[[8]][6], size = 1, linetype = l.use[[17]][3]) + geom_abline(data = df.x, aes(intercept = median(ii.an), slope = median(ii.an.x)), color = l.use[[8]][1], size = 1, linetype = l.use[[17]][1]) + geom_abline(data = df.x, aes(intercept = median(ii.pe), slope = median(ii.pe.x)), color = l.use[[8]][2], size = 1, linetype = l.use[[17]][2]) + geom_abline(data = df.x, aes(intercept = median(ii.wo), slope = median(ii.wo.x)), color = l.use[[8]][3], size = 1, linetype = l.use[[17]][3]) + scale_x_continuous(breaks = (l.use[[18]](l.use[[10]]) - l.use[[20]]) / l.use[[21]], labels = l.use[[11]], expand = c(0, 0)) + scale_y_continuous(breaks = (l.use[[19]](l.use[[12]]) - l.use[[22]]) / l.use[[23]], labels = l.use[[13]], expand = c(0, 0)) + xlab(l.use[[14]]) + ylab(l.use[[15]]) + f.theme(l.use[[16]]) } f.plot.slopes.cult <- function(l.use) { df.x <- l.use[[1]][[7]] names(df.x) <- paste(names(df.x), ".x", sep = "") df.x <- sample_n(cbind(l.use[[1]][[5]], df.x), 1000) ggplot(data = l.use[[2]], aes(x = get(l.use[[3]]), y = get(l.use[[4]]), color = comptype)) + geom_point(size = l.use[[5]]) + scale_color_manual(values = l.use[[6]]) + geom_abline(data = df.x, aes(intercept = ni.no, slope = ni.no.x), color = l.use[[7]][4], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ni.ye, slope = ni.ye.x), color = l.use[[7]][5], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ni.wi, slope = ni.wi.x), color = l.use[[7]][6], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ii.no, slope = ii.no.x), color = l.use[[7]][1], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ii.ye, slope = ii.ye.x), color = l.use[[7]][2], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = ii.wi, slope = ii.wi.x), color = l.use[[7]][3], alpha = l.use[[9]]) + geom_abline(data = df.x, aes(intercept = median(ni.no), slope = median(ni.no.x)), color = l.use[[8]][4], size = 1, linetype = l.use[[17]][1]) + geom_abline(data = df.x, aes(intercept = median(ni.ye), slope = median(ni.ye.x)), color = l.use[[8]][5], size = 1, linetype = l.use[[17]][2]) + geom_abline(data = df.x, aes(intercept = median(ni.wi), slope = median(ni.wi.x)), color = l.use[[8]][6], size = 1, linetype = l.use[[17]][3]) + geom_abline(data = df.x, aes(intercept = median(ii.no), slope = median(ii.no.x)), color = l.use[[8]][1], size = 1, linetype = l.use[[17]][1]) + geom_abline(data = df.x, aes(intercept = median(ii.ye), slope = median(ii.ye.x)), color = l.use[[8]][2], size = 1, linetype = l.use[[17]][2]) + geom_abline(data = df.x, aes(intercept = median(ii.wi), slope = median(ii.wi.x)), color = l.use[[8]][3], size = 1, linetype = l.use[[17]][3]) + scale_x_continuous(breaks = (l.use[[18]](l.use[[10]]) - l.use[[20]]) / l.use[[21]], labels = l.use[[11]], expand = c(0, 0)) + scale_y_continuous(breaks = (l.use[[19]](l.use[[12]]) - l.use[[22]]) / l.use[[23]], labels = l.use[[13]], expand = c(0, 0)) + xlab(l.use[[14]]) + ylab(l.use[[15]]) + f.theme(l.use[[16]]) } v.slopebreaks <- seq(-3, 6, by = 0.1) v.slopelabels <- c(-3, rep("", 9), -2, rep("", 9), -1, rep("", 9), 0, rep("", 9), 1, rep("", 9), 2, rep("", 9), 3, rep("", 9), 4, rep("", 9), 5, rep("", 9), 6) ch.cols2 <- c("#d73027", "#4575b4") ch.cols6lite <- c(rep("#ea8f8a", 3), rep("#97b3d7", 3)) ch.cols6 <- c(rep("#d73027", 3), rep("#4575b4", 3)) ch.type <- c("solid", "longdash", "dotted") v.bigb <- c(seq(0, 1000, by = 100), seq(2000, 10000, by = 1000), seq(20000, 100000, by = 10000), seq(200000, 1000000, by = 100000), seq(2000000, 10000000, by = 1000000)) v.bigl <- c(0, rep("", 9), "1000", rep("", 8), "10000", rep("", 8), "100000", rep("", 8), "1000000", rep("", 8), "10000000") v.smb <- c(0, seq(0.1, 0.7, by = 0.2)) v.sml <- c(0, seq(0.1, 0.7, by = 0.2)) v.smd <- seq(0, 1, by = 0.1) v.cb <- c(seq(0, 10, by = 1), seq(20, 100, by = 10), seq(200, 1000, by = 100), seq(2000, 10000, by = 1000)) v.cl <- c(0, rep("", 9), 10, rep("", 8), 100, rep("", 8), 1000, rep("", 8), 10000) v.bigb2 <- c(seq(0, 10000, by = 1000), seq(20000, 100000, by = 10000), seq(200000, 1000000, by = 100000), seq(2000000, 10000000, by = 1000000)) v.bigl2 <- c(0, rep("", 9), "10000", rep("", 8), "100000", rep("", 8), "1000000", rep("", 8), "10000000") f.as <- function(i) asin(sqrt(i)) png("040-Fig-Slopes-RangeSize.png", width = 6, height = 25, units = "in", res = 320, pointsize = 15) f.multiplot(list( list(l.m.nn[[1]], df.use, "cz.bias.nn", "z.bsd", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.bigb2, v.bigl2, v.smd, v.smd, "Native Range Size", "Geographic Overlap (Dg)", c(10, 10), ch.type, f.cube, identity, mean(f.cube(df.use$bias.n.50b), na.rm = T), sd(f.cube(df.use$bias.n.50b), na.rm = T), mean(df.use$bias.stat.d, na.rm = T), sd(df.use$bias.stat.d, na.rm = T)), list(l.m.nn[[3]], df.use, "cz.bias.nn", "cz.sr", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.bigb2, v.bigl2, v.smd, v.smd, "Native Range Size", "Proportion Stability (%Sg)", c(10, 10), ch.type, f.cube, f.as, mean(f.cube(df.use$bias.n.50b), na.rm = T), sd(f.cube(df.use$bias.n.50b), na.rm = T), mean(df.use$bias.sr, na.rm = T), sd(df.use$bias.sr, na.rm = T)), list(l.m.nn[[2]], df.use, "cz.bias.nn", "cz.ur", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.bigb2, v.bigl2, v.smd, v.smd, "Native Range Size", "Proportion Unfilling (%Ug)", c(10, 10), ch.type, f.cube, f.as, mean(f.cube(df.use$bias.n.50b), na.rm = T), sd(f.cube(df.use$bias.n.50b), na.rm = T), mean(df.use$bias.ur, na.rm = T), sd(df.use$bias.ur, na.rm = T)), list(l.m.nn[[4]], df.use, "cz.bias.nn", "cz.er", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.bigb2, v.bigl2, v.smd, v.smd, "Native Range Size", "Proportion Expansion (%Eg)", c(10, 10), ch.type, f.cube, f.as, mean(f.cube(df.use$bias.n.50b), na.rm = T), sd(f.cube(df.use$bias.n.50b), na.rm = T), mean(df.use$bias.er, na.rm = T), sd(df.use$bias.er, na.rm = T)), list(l.m.nn[[5]], df.use, "cz.bias.nn", "cz.bias.n", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.bigb2, v.bigl2, v.bigb2, v.bigl2, "Native Range Size", "Introduced Range Size", c(10, 10), ch.type, f.cube, f.cube, mean(f.cube(df.use$bias.n.50b), na.rm = T), sd(f.cube(df.use$bias.n.50b), na.rm = T), mean(f.cube(df.use$bias.n.50a), na.rm = T), sd(f.cube(df.use$bias.n.50a), na.rm = T)) ), f.plot.slopes, 1 ) dev.off() png("040-Fig-Slopes.png", width = 12, height = 11, units = "in", res = 320, pointsize = 15) f.multiplot(list( list(l.m.d[[1]], df.use, "z.d", "z.bsd", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.smd, v.smd, v.smd, v.smd, "Climatic Overlap (Dc)", "Geographic Overlap (Dg)", c(10, 10), ch.type, identity, identity, mean(df.use$d, na.rm = T), sd(df.use$d, na.rm = T), mean(df.use$bias.stat.d, na.rm = T), sd(df.use$bias.stat.d, na.rm = T)), list(l.m.uc[[4]], df.use, "cz.uc", "cz.bias.u", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.cb, v.cl, v.bigb, v.bigl, "Climatic Unfilling (Uc)", "Geographic Unfilling (Ug)", c(10, 10), ch.type, f.cube, f.cube, mean(f.cube(df.use$u.carea), na.rm = T), sd(f.cube(df.use$u.carea), na.rm = T), mean(f.cube(df.use$bias.u50), na.rm = T), sd(f.cube(df.use$bias.u50), na.rm = T)), list(l.m.sc[[3]], df.use, "cz.sc", "cz.bias.s", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.cb, v.cl, v.bigb, v.bigl, "Climatic Stability (Sc)", "Geographic Stability (Sg)", c(10, 10), ch.type, f.cube, f.cube, mean(f.cube(df.use$s.carea), na.rm = T), sd(f.cube(df.use$s.carea), na.rm = T), mean(f.cube(df.use$bias.s50), na.rm = T), sd(f.cube(df.use$bias.s50), na.rm = T)), list(l.m.ec[[5]], df.use, "cz.ec", "cz.bias.e", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.cb, v.cl, v.bigb, v.bigl, "Climatic Expansion (Ec)", "Geographic Expansion (Eg)", c(10, 10), ch.type, f.cube, f.cube, mean(f.cube(df.use$e.carea), na.rm = T), sd(f.cube(df.use$e.carea), na.rm = T), mean(f.cube(df.use$bias.e50), na.rm = T), sd(f.cube(df.use$bias.e50), na.rm = T)) ), f.plot.slopes, 2 ) dev.off() png("040-Fig-Slopes-CULT.png", width = 12, height = 11, units = "in", res = 320, pointsize = 15) f.multiplot(list( list(l.m.d.cult[[1]], df.use, "z.d", "z.bsd", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.smd, v.smd, v.smd, v.smd, "Climatic Overlap (Dc)", "Geographic Overlap (Dg)", c(10, 10), ch.type, identity, identity, mean(df.use$d, na.rm = T), sd(df.use$d, na.rm = T), mean(df.use$bias.stat.d, na.rm = T), sd(df.use$bias.stat.d, na.rm = T)), list(l.m.uc.cult[[4]], df.use, "cz.uc", "cz.bias.u", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.cb, v.cl, v.bigb, v.bigl, "Climatic Unfilling (Uc)", "Geographic Unfilling (Ug)", c(10, 10), ch.type, f.cube, f.cube, mean(f.cube(df.use$u.carea), na.rm = T), sd(f.cube(df.use$u.carea), na.rm = T), mean(f.cube(df.use$bias.u50), na.rm = T), sd(f.cube(df.use$bias.u50), na.rm = T)), list(l.m.sc.cult[[3]], df.use, "cz.sc", "cz.bias.s", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.cb, v.cl, v.bigb, v.bigl, "Climatic Stability (Sc)", "Geographic Stability (Sg)", c(10, 10), ch.type, f.cube, f.cube, mean(f.cube(df.use$s.carea), na.rm = T), sd(f.cube(df.use$s.carea), na.rm = T), mean(f.cube(df.use$bias.s50), na.rm = T), sd(f.cube(df.use$bias.s50), na.rm = T)), list(l.m.ec.cult[[5]], df.use, "cz.ec", "cz.bias.e", 2, ch.cols2, ch.cols6lite, ch.cols6, 0.025, v.cb, v.cl, v.bigb, v.bigl, "Climatic Expansion (Ec)", "Geographic Expansion (Eg)", c(10, 10), ch.type, f.cube, f.cube, mean(f.cube(df.use$e.carea), na.rm = T), sd(f.cube(df.use$e.carea), na.rm = T), mean(f.cube(df.use$bias.e50), na.rm = T), sd(f.cube(df.use$bias.e50), na.rm = T)) ), f.plot.slopes.cult, 2 ) dev.off() t.plot <- theme( axis.text = element_text(size = 14, color = "black"), axis.title.y = element_blank(), axis.title.x = element_text(size = 16), axis.line.x = element_line(color = "black"), axis.line.y = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.ticks.y = element_blank(), panel.grid = element_blank(), panel.background = element_rect(fill = "white"), legend.position = "none" ) f.plot.seg <- function(l.use) ggplot(data = l.use[[1]]) + geom_segment(aes(x = s02.5, xend = s97.5, y = s.i, yend = s.i), size = l.use[[3]], color = l.use[[4]]) + geom_rect(aes(xmin = s25.0, xmax = s75.0, ymin = s.i - l.use[[6]], ymax = s.i + l.use[[6]], fill = st.type), size = l.use[[3]], color = l.use[[4]]) + scale_fill_manual(values = l.use[[5]]) + geom_vline(xintercept = 0, linetype = "dotted") + scale_y_continuous(breaks = 1:6, labels = l.use[[12]]) + scale_x_continuous(breaks = l.use[[8]], labels = l.use[[9]], limits = l.use[[10]]) + ylab("") + xlab(l.use[[11]]) + t.plot + geom_text(aes(x = s97.5 + s.bump, y = s.i, label = ifelse(is.na(pdif), "", format(round(pdif, digits = 3), nsmall = 3)))) v.breakuse <- seq(-2, 2, by = 0.5) v.limuse <- c(-0.75, 2) v.limuse2 <- c(-0.5, 2) s.size <- 0.5 s.bump <- 0.3 png("040-Fig-Coef-BEST.png", width = 15, height = 5, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.d[[1]][[8]], l.m.d[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse, "Overlap (D)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.d.cult[[1]][[8]], l.m.d.cult[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse, "Overlap (D)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.sc[[3]][[8]], l.m.sc[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse, "Stability (Sc)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.sc.cult[[3]][[8]], l.m.sc.cult[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse, "Stability (Sc)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.uc[[4]][[8]], l.m.uc[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Unfilling (Uc)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.uc.cult[[4]][[8]], l.m.uc.cult[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Unfilling (Uc)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.ec[[5]][[8]], l.m.ec[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Expansion (Ec)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.ec.cult[[5]][[8]], l.m.ec.cult[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Expansion (Ec)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")) ), f.plot.seg, 4 ) dev.off() png("040-Fig-Coef-Selected.png", width = 15, height = 5, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.d[[1]][[8]], l.m.d[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse, "Overlap (D)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.d[[1]][[8]], l.m.d[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse, "Overlap (D)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.sc[[3]][[8]], l.m.sc[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse, "Stability (Sc)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.s[[3]][[8]], l.m.s[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse, "Stability (Scg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.ec[[5]][[8]], l.m.ec[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Expansion (Ec)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.e[[5]][[8]], l.m.e[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Expansion (Ecg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.uc[[4]][[8]], l.m.uc[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Unfilling (Uc)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.u[[4]][[8]], l.m.u[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Unfilling (Ucg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")) ), f.plot.seg, 4 ) dev.off() png("040-Fig-Coef-Selected-CULT.png", width = 15, height = 5, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.d.cult[[1]][[8]], l.m.d.cult[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse2, "Overlap (D)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.d.cult[[1]][[8]], l.m.d.cult[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse2, "Overlap (D)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.sc.cult[[3]][[8]], l.m.sc.cult[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse2, "Stability (Sc)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.s.cult[[3]][[8]], l.m.s.cult[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse2, "Stability (Scg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.ec.cult[[5]][[8]], l.m.ec.cult[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Expansion (Ec)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.e.cult[[5]][[8]], l.m.e.cult[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Expansion (Ecg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.uc.cult[[4]][[8]], l.m.uc.cult[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Unfilling (Uc)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.u.cult[[4]][[8]], l.m.u.cult[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Unfilling (Ucg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")) ), f.plot.seg, 4 ) dev.off() png("040-Fig-Coef-D.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.d[[1]][[8]], l.m.d[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse, "Overlap (D)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.d[[2]][[8]], l.m.d[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.18, v.breakuse, v.breakuse, v.limuse, "Range Shift (R)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.d[[6]][[8]], l.m.d[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.11, v.breakuse, v.breakuse, v.limuse, "Inv. Range Size", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.d[[3]][[8]], l.m.d[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Stability (Sg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.d[[4]][[8]], l.m.d[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Unfilling (Ug)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.d[[5]][[8]], l.m.d[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Expansion (Eg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-CentX.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.x[[1]][[8]], l.m.x[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, c(-1.25, 0.75), "Overlap (D)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.x[[2]][[8]], l.m.x[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.18, v.breakuse, v.breakuse, c(-1.25, 0.75), "Range Shift (R)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.x[[6]][[8]], l.m.x[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.11, v.breakuse, v.breakuse, c(-1.25, 0.75), "Inv. Range Size", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.x[[3]][[8]], l.m.x[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, c(-1.25, 0.75), "Stability (Sg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.x[[4]][[8]], l.m.x[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, c(-1.25, 0.75), "Unfilling (Ug)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.x[[5]][[8]], l.m.x[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, c(-1.25, 0.75), "Expansion (Eg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-U.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.uc[[1]][[8]], l.m.uc[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse, "Overlap (D)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.uc[[2]][[8]], l.m.uc[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.18, v.breakuse, v.breakuse, v.limuse, "Range Shift (R)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.uc[[6]][[8]], l.m.uc[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.11, v.breakuse, v.breakuse, v.limuse, "Inv. Range Size", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.uc[[3]][[8]], l.m.uc[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Stability (Sg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.uc[[4]][[8]], l.m.uc[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Unfilling (Ug)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.uc[[5]][[8]], l.m.uc[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Expansion (Eg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-S.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.sc[[1]][[8]], l.m.sc[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.33, v.breakuse, v.breakuse, v.limuse, "Overlap (D)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.sc[[2]][[8]], l.m.sc[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.88, v.breakuse, v.breakuse, v.limuse, "Range Shift (R)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.sc[[6]][[8]], l.m.sc[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.38, v.breakuse, v.breakuse, v.limuse, "Inv. Range Size", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.sc[[3]][[8]], l.m.sc[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse, "Stability (Sg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.sc[[4]][[8]], l.m.sc[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse, "Unfilling (Ug)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.sc[[5]][[8]], l.m.sc[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse, "Expansion (Eg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-E.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.ec[[1]][[8]], l.m.ec[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Overlap (D)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.ec[[2]][[8]], l.m.ec[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.24, v.breakuse, v.breakuse, v.limuse, "Range Shift (R)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.ec[[6]][[8]], l.m.ec[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.17, v.breakuse, v.breakuse, v.limuse, "Inv. Range Size", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.ec[[3]][[8]], l.m.ec[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Stability (Sg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.ec[[4]][[8]], l.m.ec[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Unfilling (Ug)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")), list(l.m.ec[[5]][[8]], l.m.ec[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse, "Expansion (Eg)", c("NI-A", "NI-P", "NI-W", "NN-A", "NN-P", "NN-W")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-D-CULT.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.d.cult[[1]][[8]], l.m.d.cult[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse2, "Overlap (D)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.d.cult[[2]][[8]], l.m.d.cult[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.18, v.breakuse, v.breakuse, v.limuse2, "Range Shift (R)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.d.cult[[6]][[8]], l.m.d.cult[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.11, v.breakuse, v.breakuse, v.limuse2, "Inv. Range Size", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.d.cult[[3]][[8]], l.m.d.cult[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Stability (Sg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.d.cult[[5]][[8]], l.m.d.cult[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Expansion (Eg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.d.cult[[4]][[8]], l.m.d.cult[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Unfilling (Ug)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-CentX-CULT.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.x.cult[[1]][[8]], l.m.x.cult[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse2, "Overlap (D)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.x.cult[[2]][[8]], l.m.x.cult[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.18, v.breakuse, v.breakuse, v.limuse2, "Range Shift (R)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.x.cult[[6]][[8]], l.m.x.cult[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.11, v.breakuse, v.breakuse, v.limuse2, "Inv. Range Size", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.x.cult[[3]][[8]], l.m.x.cult[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Stability (Sg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.x.cult[[5]][[8]], l.m.x.cult[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Expansion (Eg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.x.cult[[4]][[8]], l.m.x.cult[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Unfilling (Ug)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-U-CULT.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.uc.cult[[1]][[8]], l.m.uc.cult[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.06, v.breakuse, v.breakuse, v.limuse2, "Overlap (D)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.uc.cult[[2]][[8]], l.m.uc.cult[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.18, v.breakuse, v.breakuse, v.limuse2, "Range Shift (R)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.uc.cult[[6]][[8]], l.m.uc.cult[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.11, v.breakuse, v.breakuse, v.limuse2, "Inv. Range Size", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.uc.cult[[3]][[8]], l.m.uc.cult[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Stability (Sg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.uc.cult[[5]][[8]], l.m.uc.cult[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Expansion (Eg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.uc.cult[[4]][[8]], l.m.uc.cult[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Unfilling (Ug)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-S-CULT.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.sc.cult[[1]][[8]], l.m.sc.cult[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.33, v.breakuse, v.breakuse, v.limuse2, "Overlap (D)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.sc.cult[[2]][[8]], l.m.sc.cult[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.88, v.breakuse, v.breakuse, v.limuse2, "Range Shift (R)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.sc.cult[[6]][[8]], l.m.sc.cult[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.38, v.breakuse, v.breakuse, v.limuse2, "Inv. Range Size", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.sc.cult[[3]][[8]], l.m.sc.cult[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse2, "Stability (Sg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.sc.cult[[5]][[8]], l.m.sc.cult[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse2, "Expansion (Eg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.sc.cult[[4]][[8]], l.m.sc.cult[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.61, v.breakuse, v.breakuse, v.limuse2, "Unfilling (Ug)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")) ), f.plot.seg, 2 ) dev.off() png("040-Fig-Coef-E-CULT.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.m.ec.cult[[1]][[8]], l.m.ec.cult[[1]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Overlap (D)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.ec.cult[[2]][[8]], l.m.ec.cult[[2]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.24, v.breakuse, v.breakuse, v.limuse2, "Range Shift (R)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.ec.cult[[6]][[8]], l.m.ec.cult[[6]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.17, v.breakuse, v.breakuse, v.limuse2, "Inv. Range Size", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.ec.cult[[3]][[8]], l.m.ec.cult[[3]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Stability (Sg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.ec.cult[[5]][[8]], l.m.ec.cult[[5]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Expansion (Eg)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")), list(l.m.ec.cult[[4]][[8]], l.m.ec.cult[[4]][[8]], s.size, "black", c("grey", "white"), 0.08, 0.16, v.breakuse, v.breakuse, v.limuse2, "Unfilling (Ug)", c("NI-No", "NI-Yes", "NI-Wide", "NN-No", "NN-Yes", "NN-Wide")) ), f.plot.seg, 2 ) dev.off() df.t1 <- df.use[, c("sp.id", "comptype", "cultivated", "growth", "life", "cult", "reg", "wood", "bias.stat.d", "bias.mean.e", "bias.mean.e.exp", "bias.mean.e.unf")] df.t2 <- df.use[, c("sp.id", "comptype", "cultivated", "growth", "life", "cult", "reg", "wood", "flat.stat.d", "flat.mean.e", "flat.mean.e.exp", "flat.mean.e.unf")] df.t3 <- df.use[, c("sp.id", "comptype", "cultivated", "growth", "life", "cult", "reg", "wood", "c.bias.stat.d", "c.bias.mean.e", "c.bias.mean.e.exp", "c.bias.mean.e.unf")] df.t4 <- df.use[, c("sp.id", "comptype", "cultivated", "growth", "life", "cult", "reg", "wood", "c.flat.stat.d", "c.flat.mean.e", "c.flat.mean.e.exp", "c.flat.mean.e.unf")] df.t5 <- df.use[, c("sp.id", "comptype", "cultivated", "growth", "life", "cult", "reg", "wood", "d", "dc", "dc.e", "dc.u")] names(df.t1) <- names(df.t2) <- names(df.t3) <- names(df.t4) <- names(df.t5) <- c( "sp.id", "comptype", "cultivated", "growth", "life", "cult", "reg", "wood", "d", "dc", "dc.e", "dc.u") df.t1$mod <- "max.bias" df.t2$mod <- "max.flat" df.t3$mod <- "max.c.bias" df.t4$mod <- "max.c.flat" df.t5$mod <- "ord.bias" df.test <- rbind(df.t1, df.t2, df.t3, df.t4, df.t5) df.test$comptype <- factor(df.test$comptype) df.test$mc <- factor(paste(df.test$mod, ".", df.test$comptype, sep = "")) df.test$mod <- factor(df.test$mod) df.test2 <- df.test[df.test$mod == "max.bias" | df.test$mod == "ord.bias",] df.test2$mod <- factor(df.test2$mod) df.test2$comptype <- factor(df.test2$comptype) m.test.d <- stan_lmer(data = df.test2, d ~ mod * comptype + (1|sp.id) + (1|reg)) pp_check(m.test.d) b.test.d <- as_tibble(m.test.d) %>% rename(intercept = `(Intercept)`) %>% select(-sigma) df.b.test.d <- data.frame( ni.max = with(b.test.d, intercept), ni.ord = with(b.test.d, intercept + modord.bias), nn.max = with(b.test.d, intercept + comptypeNN), nn.ord = with(b.test.d, intercept + modord.bias + comptypeNN + `modord.bias:comptypeNN`) ) df.seg.test.d <- do.call(rbind, lapply(1:ncol(df.b.test.d), function(i) data.frame( s.i = i, st.var = names(df.b.test.d)[i], st.type = substr(names(df.b.test.d)[i], 1, 2), s02.5 = quantile(df.b.test.d[, i], 0.025), s25.0 = quantile(df.b.test.d[, i], 0.250), s75.0 = quantile(df.b.test.d[, i], 0.750), s87.5 = quantile(df.b.test.d[, i], 0.875), s97.5 = quantile(df.b.test.d[, i], 0.975) ))) png("040-Fig-OrdVSMax.png", width = 9, height = 3, units = "in", res = 320, pointsize = 15) ggplot(data = df.seg.test.d) + geom_segment(aes(x = s02.5, xend = s97.5, y = s.i, yend = s.i), size = 0.5, color = "black") + geom_rect(aes(xmin = s25.0, xmax = s75.0, ymin = s.i - 0.08, ymax = s.i + 0.08, fill = st.type), size = 0.5, color = "black") + scale_fill_manual(values = c("grey", "white")) + geom_vline(xintercept = 0, linetype = "dotted") + scale_y_continuous(breaks = 1:4, labels = c("NI-Geo.", "NI-Clim.", "NN-Geo.", "NN-Clim.")) + scale_x_continuous(breaks = seq(-1, 1, by = 0.1), limits = c(-0.1, 0.7)) + ylab("") + xlab("Overlap (D)") + t.plot dev.off() mx.pdif.test.d <- matrix( c( 0, sum(df.b.test.d$ni.max > df.b.test.d$ni.ord), sum(df.b.test.d$ni.max > df.b.test.d$nn.max), sum(df.b.test.d$ni.max > df.b.test.d$nn.ord), sum(df.b.test.d$ni.ord > df.b.test.d$ni.max), 0, sum(df.b.test.d$ni.ord > df.b.test.d$nn.max), sum(df.b.test.d$ni.ord > df.b.test.d$nn.ord), sum(df.b.test.d$nn.max > df.b.test.d$ni.max), sum(df.b.test.d$nn.max > df.b.test.d$ni.ord), 0, sum(df.b.test.d$nn.max > df.b.test.d$nn.ord), sum(df.b.test.d$nn.ord > df.b.test.d$ni.max), sum(df.b.test.d$nn.ord > df.b.test.d$ni.ord), sum(df.b.test.d$nn.ord > df.b.test.d$nn.max), 0 ), ncol = 4 ) mx.pdif.test.d <- mx.pdif.test.d / nrow(df.b.test.d) mx.pdif.test.d[mx.pdif.test.d > 0.5] <- 1 - mx.pdif.test.d[mx.pdif.test.d > 0.5] load("040-MXDiffStats.RData") df.use$u.n <- df.use$bias.u50 / (df.use$bias.u50 + df.use$bias.s50 + df.use$bias.e50) df.use$s.n <- df.use$bias.s50 / (df.use$bias.u50 + df.use$bias.s50 + df.use$bias.e50) df.use$e.n <- df.use$bias.e50 / (df.use$bias.u50 + df.use$bias.s50 + df.use$bias.e50) f.as <- function(i) asin(sqrt(i)) l.f <- list(identity, f.cube, f.cube, f.cube, f.cube, identity) ch.get <- c("bias.stat.d", "bias.s50", "bias.u50", "bias.e50", "bias.n.50a", "bias.rangediff") df.qs <- do.call(rbind, lapply(1:length(l.f), function(i) data.frame( vname = rep(ch.get[i], 2), comptype = c("NN", "NI"), u.var = c(mean(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), na.rm = T), mean(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), na.rm = T)), u.sd = c(sd(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), na.rm = T), sd(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), na.rm = T)), q.025 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.025, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.025, na.rm = T)), q.050 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.050, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.050, na.rm = T)), q.125 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.125, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.125, na.rm = T)), q.250 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.250, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.250, na.rm = T)), q.375 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.375, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.375, na.rm = T)), q.500 = c(median(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), na.rm = T), median(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), na.rm = T)), q.625 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.625, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.625, na.rm = T)), q.750 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.750, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.750, na.rm = T)), q.875 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.875, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.875, na.rm = T)), q.950 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.950, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.950, na.rm = T)), q.975 = c(quantile(l.f[[i]](df.use[df.use$comptype == "NN", ch.get[i]]), 0.975, na.rm = T), quantile(l.f[[i]](df.use[df.use$comptype == "NI", ch.get[i]]), 0.975, na.rm = T)) ))) df.qs$s.i <- 12:1 png("040-Fig-RawData-D.png", width = 9, height = 2, units = "in", res = 320, pointsize = 15) ggplot(data = df.qs[df.qs$s.i > 10,]) + geom_segment(aes(x = q.025, xend = q.975, y = s.i, yend = s.i), size = 0.5, color = "black") + geom_rect(aes(xmin = q.250, xmax = q.750, ymin = s.i - 0.08, ymax = s.i + 0.08, fill = comptype), size = 0.5, color = "black") + geom_point(aes(x = q.500, y = s.i), size = 2, color = "black") + geom_text(aes(x = q.500, y = s.i + 0.25, label = format(round(q.500, digits = 2), nsmall = 2))) + scale_fill_manual(values = c("grey", "white")) + geom_vline(xintercept = 0, linetype = "dotted") + scale_y_continuous(breaks = 11.5, labels = "D", limits = c(10.5, 12.5)) + scale_x_continuous(breaks = seq(0, 1, by = 0.1), labels = seq(0, 1, by = 0.1)) + t.plot dev.off() png("040-Fig-RawData-USEN.png", width = 9, height = 6, units = "in", res = 320, pointsize = 15) ggplot(data = df.qs[df.qs$s.i > 2 & df.qs$s.i < 11,]) + geom_segment(aes(x = q.025, xend = q.975, y = s.i, yend = s.i), size = 0.5, color = "black") + geom_rect(aes(xmin = q.250, xmax = q.750, ymin = s.i - 0.08, ymax = s.i + 0.08, fill = comptype), size = 0.5, color = "black") + geom_point(aes(x = q.500, y = s.i), size = 2, color = "black") + geom_text(aes(x = q.500, y = s.i + 0.25, label = format(round(q.500^3, digits = 2), nsmall = 2))) + scale_fill_manual(values = c("grey", "white")) + geom_vline(xintercept = 0, linetype = "dotted") + scale_y_continuous(breaks = seq(3.5, 9.5, by = 2), labels = c("N", "E", "U", "S"), limits = c(2.75, 10.25)) + scale_x_continuous(breaks = v.bigb^(1/3), labels = v.bigl) + t.plot dev.off() png("040-Fig-RawData-R.png", width = 9, height = 2, units = "in", res = 320, pointsize = 15) ggplot(data = df.qs[df.qs$s.i < 3,]) + geom_segment(aes(x = q.025, xend = q.975, y = s.i, yend = s.i), size = 0.5, color = "black") + geom_rect(aes(xmin = q.250, xmax = q.750, ymin = s.i - 0.08, ymax = s.i + 0.08, fill = comptype), size = 0.5, color = "black") + geom_point(aes(x = q.500, y = s.i), size = 2, color = "black") + geom_text(aes(x = q.500, y = s.i + 0.25, label = format(round(q.500, digits = 2), nsmall = 2))) + scale_fill_manual(values = c("grey", "white")) + geom_vline(xintercept = 0, linetype = "dotted") + scale_y_continuous(breaks = 1.5, labels = "R", limits = c(0.5, 2.5)) + scale_x_continuous( breaks = seq(-1000000, 400000, by = 100000), labels = c("-1000000", "", "-800000", "", "-600000", "", "-400000", "", "-200000", "", "0", "", "200000", "", "400000") ) + t.plot dev.off() f.pairs <- function(l.use) { df.out <- data.frame( comp = paste(names(l.use[[7]])[c(1, 1, 2)], "-", names(l.use[[7]])[c(2, 3, 3)], sep = ""), st.y = rep(l.use[[1]], 3), st.x = rep(l.use[[2]], 3), v.difs = c( sum(l.use[[7]][,1] > l.use[[7]][,2]) / nrow(l.use[[7]]), sum(l.use[[7]][,1] > l.use[[7]][,3]) / nrow(l.use[[7]]), sum(l.use[[7]][,2] > l.use[[7]][,3]) / nrow(l.use[[7]]) ) ) df.out$v.difsp <- sapply(df.out$v.difs, function(j) ifelse(j < 0.5, j*2, (1 - j)*2)) df.out } df.lifepairs.c <- rbind(f.pairs(l.m.d[[1]]), f.pairs(l.m.sc[[3]]), f.pairs(l.m.uc[[4]]), f.pairs(l.m.ec[[5]])) df.cultpairs.c <- rbind(f.pairs(l.m.d.cult[[1]]), f.pairs(l.m.sc.cult[[3]]), f.pairs(l.m.uc.cult[[4]]), f.pairs(l.m.ec.cult[[5]])) write.table(df.lifepairs.c, "040-Results-LifePairsC.csv", sep = ",", na = "", row.names = F) write.table(df.cultpairs.c, "040-Results-CultPairsC.csv", sep = ",", na = "", row.names = F) f.nieffectcomp <- function(df.i, st.x, s.len, l.i, ch.y, st.xalt) { v.r <- seq(from = min(df.i[, st.x], na.rm = T), to = max(df.i[, st.x], na.rm = T), length.out = s.len) list( st.x, st.xalt, data.frame( v.r = rep(v.r, times = length(ch.y)), v.out = c( sapply(v.r, function(j) sum( (l.i[[5]][,paste("ii.", ch.y[1], sep = "")] + l.i[[7]][,paste("ii.", ch.y[1], sep = "")] * j) > (l.i[[5]][,paste("ni.", ch.y[1], sep = "")] + l.i[[7]][,paste("ni.", ch.y[1], sep = "")] * j) ) / nrow(l.i[[5]])), sapply(v.r, function(j) sum( (l.i[[5]][,paste("ii.", ch.y[2], sep = "")] + l.i[[7]][,paste("ii.", ch.y[2], sep = "")] * j) > (l.i[[5]][,paste("ni.", ch.y[2], sep = "")] + l.i[[7]][,paste("ni.", ch.y[2], sep = "")] * j) ) / nrow(l.i[[5]])), sapply(v.r, function(j) sum( (l.i[[5]][,paste("ii.", ch.y[3], sep = "")] + l.i[[7]][,paste("ii.", ch.y[3], sep = "")] * j) > (l.i[[5]][,paste("ni.", ch.y[3], sep = "")] + l.i[[7]][,paste("ni.", ch.y[3], sep = "")] * j) ) / nrow(l.i[[5]])) ), st.comp = rep(ch.y, each = s.len) ) ) } l.nief.d <- f.nieffectcomp(df.use, "z.d", 101, l.m.d[[1]], c("an", "pe", "wo"), "d") l.nief.s <- f.nieffectcomp(df.use, "cz.sc", 101, l.m.sc[[3]], c("an", "pe", "wo"), "s.carea") l.nief.u <- f.nieffectcomp(df.use, "cz.uc", 101, l.m.uc[[4]], c("an", "pe", "wo"), "u.carea") l.nief.e <- f.nieffectcomp(df.use, "cz.ec", 101, l.m.ec[[4]], c("an", "pe", "wo"), "e.carea") l.nief.d.cult <- f.nieffectcomp(df.use, "z.d", 101, l.m.d.cult[[1]], c("no", "ye", "wi"), "d") l.nief.s.cult <- f.nieffectcomp(df.use, "cz.sc", 101, l.m.sc.cult[[3]], c("no", "ye", "wi"), "s.carea") l.nief.u.cult <- f.nieffectcomp(df.use, "cz.uc", 101, l.m.uc.cult[[4]], c("no", "ye", "wi"), "u.carea") l.nief.e.cult <- f.nieffectcomp(df.use, "cz.ec", 101, l.m.ec.cult[[4]], c("no", "ye", "wi"), "e.carea") f.plotniefcomp <- function(l.i) { ggplot(data = l.i[[1]][[3]], aes(x = v.r, y = v.out, linetype = st.comp)) + geom_hline(yintercept = 0.025, size = 0.5, linetype = "dashed", color = "grey") + geom_hline(yintercept = 0.5, size = 0.5, linetype = "dashed", color = "grey") + geom_hline(yintercept = 0.975, size = 0.5, linetype = "dashed", color = "grey") + geom_line(size = 1) + scale_linetype_manual(values = c("solid", "longdash", "dotted")) + scale_y_continuous( breaks = seq(0, 1, by = 0.1), labels = c(0, "", 0.2, "", 0.4, "", 0.6, "", 0.8, "", 1), limits = c(0, 1), expand = c(0, 0) ) + scale_x_continuous( breaks = (l.i[[2]](l.i[[3]]) - l.i[[2]](mean(l.i[[5]][, l.i[[1]][[2]]], na.rm = T))) / sd(l.i[[2]](l.i[[5]][, l.i[[1]][[2]]]), na.rm = T), labels = l.i[[4]], expand = c(0, 0) ) + ylab(paste("Prop. N-I > N-N (", l.i[[7]], "g)", sep = "")) + xlab(l.i[[6]]) + theme( axis.text = element_text(size = 14, color = "black"), axis.title.y = element_text(size = 16), axis.title.x = element_text(size = 16), axis.line.x = element_line(color = "black"), axis.line.y = element_line(color = "black"), axis.ticks = element_line(color = "black"), panel.grid = element_blank(), panel.background = element_rect(fill = "white"), legend.position = "none" ) } png("040-Fig-PlotComp.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.nief.d, identity, seq(0, 1, by = 0.1), c(0, "", 0.2, "", 0.4, "", 0.6, "", 0.8, "", 1), df.use, "Overlap (Dc)", "D"), list(l.nief.u, f.cube, v.cb, v.cl, df.use, "Unfilling (Uc)", "U"), list(l.nief.s, f.cube, v.cb, v.cl, df.use, "Stability (Sc)", "S"), list(l.nief.e, f.cube, v.cb, v.cl, df.use, "Expansion (Ec)", "E") ), f.plotniefcomp, 2 ) dev.off() png("040-Fig-PlotComp-CULT.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.multiplot( l.plot = list( list(l.nief.d.cult, identity, seq(0, 1, by = 0.1), c(0, "", 0.2, "", 0.4, "", 0.6, "", 0.8, "", 1), df.use, "Overlap (Dc)", "D"), list(l.nief.u.cult, f.cube, v.cb, v.cl, df.use, "Unfilling (Uc)", "U"), list(l.nief.s.cult, f.cube, v.cb, v.cl, df.use, "Stability (Sc)", "S"), list(l.nief.e.cult, f.cube, v.cb, v.cl, df.use, "Expansion (Ec)", "E") ), f.plotniefcomp, 2 ) dev.off() df.evals.i <- df.evals[, c("sp.id", "mtype", "reg", "e.auc.i", "m.auc.i", "r1.i", "r2.i", "r3.i", "r4.i")] df.evals.n <- df.evals[, c("sp.id", "mtype", "reg", "e.auc.n", "m.auc.n", "r1.n", "r2.n", "r3.n", "r4.n")] names(df.evals.i) <- names(df.evals.n) <- c("sp.id", "mtype", "reg", "e.auc", "m.auc", "r1", "r2", "r3", "r4") df.evals.i$ni <- rep("i", nrow(df.evals.i)) df.evals.n$ni <- rep("n", nrow(df.evals.n)) df.evals.w <- rbind(df.evals.n, df.evals.i) df.evals.w$sp.fac <- factor(df.evals.w$sp.id) df.evals.w$reg.fac <- factor(df.evals.w$reg) #m.bias.m.auc <- stan_lmer(data = df.evals.w[df.evals.w$mtype == "bias",], m.auc ~ reg.fac * ni + (1|sp.id), iter = 4000) #save(m.bias.m.auc, file = "040-Mod_AUC.Rdata") load("040-Mod_AUC.Rdata") m.auc.b <- as_tibble(m.bias.m.auc) %>% rename(intercept = `(Intercept)`) %>% select(-sigma) df.auc.b <- data.frame( i.1 = with(m.auc.b, intercept), i.2 = with(m.auc.b, intercept + reg.fac2), i.3 = with(m.auc.b, intercept + reg.fac3), i.4 = with(m.auc.b, intercept + reg.fac4), i.5 = with(m.auc.b, intercept + reg.fac5), i.6 = with(m.auc.b, intercept + reg.fac6), n.1 = with(m.auc.b, intercept + nin), n.2 = with(m.auc.b, intercept + nin + reg.fac2 + `reg.fac2:nin`), n.3 = with(m.auc.b, intercept + nin + reg.fac3 + `reg.fac3:nin`), n.4 = with(m.auc.b, intercept + nin + reg.fac4 + `reg.fac4:nin`), n.5 = with(m.auc.b, intercept + nin + reg.fac5 + `reg.fac5:nin`), n.6 = with(m.auc.b, intercept + nin + reg.fac6 + `reg.fac6:nin`) ) df.auc.seg <- do.call(rbind, lapply(1:ncol(df.auc.b), function(i) data.frame( s.i = i, st.var = names(df.auc.b)[i], st.type = substr(names(df.auc.b)[i], 1, 2), s02.5 = quantile(df.auc.b[, i], 0.025), s25.0 = quantile(df.auc.b[, i], 0.250), s50.0 = quantile(df.auc.b[, i], 0.500), s75.0 = quantile(df.auc.b[, i], 0.750), s87.5 = quantile(df.auc.b[, i], 0.875), s97.5 = quantile(df.auc.b[, i], 0.975) ))) df.auc.seg$pdif <- c( sum(df.auc.b$i.1 > df.auc.b$n.1), sum(df.auc.b$i.2 > df.auc.b$n.2), sum(df.auc.b$i.3 > df.auc.b$n.3), sum(df.auc.b$i.4 > df.auc.b$n.4), sum(df.auc.b$i.5 > df.auc.b$n.5), sum(df.auc.b$i.6 > df.auc.b$n.6), rep(NA, 6) ) / nrow(df.auc.b) df.auc.seg$s.i2 <- c(seq(1, 11, by = 2), seq(2, 12, by = 2)) t.plot <- theme( axis.text = element_text(size = 14, color = "black"), axis.title.y = element_blank(), axis.title.x = element_text(size = 16), axis.line.x = element_line(color = "black"), axis.line.y = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.ticks.y = element_blank(), panel.grid = element_blank(), panel.background = element_rect(fill = "white"), legend.position = "none" ) s.size <- 0.5 s.bump <- 0.3 png("040-Fig-AUC.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) ggplot(data = df.auc.seg) + geom_segment(aes(x = s02.5, xend = s97.5, y = s.i2, yend = s.i2), size = s.size, color = "black") + geom_rect( aes(xmin = s25.0, xmax = s75.0, ymin = s.i2 - 0.08, ymax = s.i2 + 0.08, fill = st.type), size = s.size, color = "black" ) + geom_segment(aes(x = s50.0, xend = s50.0, y = s.i2 - 0.08, yend = s.i2 + 0.08), size = s.size, color = "black") + geom_text(aes(x = 0.94, y = s.i2 + 0.5, label = ifelse(is.na(pdif), "", format(round(pdif, digits = 3), nsmall = 3)))) + geom_hline(yintercept = 2.5, linetype = "longdash", size = 0.5, color = "black") + geom_hline(yintercept = 4.5, linetype = "longdash", size = 0.5, color = "black") + geom_hline(yintercept = 6.5, linetype = "longdash", size = 0.5, color = "black") + geom_hline(yintercept = 8.5, linetype = "longdash", size = 0.5, color = "black") + geom_hline(yintercept = 10.5, linetype = "longdash", size = 0.5, color = "black") + scale_fill_manual(values = c("grey", "white")) + scale_y_continuous( breaks = 1:12, labels = paste(rep(c("I-to-I", "N-to-I"), times = 6), rep(c("AF", "AU", "EA", "IP", "NA", "SA"), each = 2), sep = ": ") ) + scale_x_continuous( breaks = seq(0.8, 1, by = 0.01), labels = c(0.8, "", 0.82, "", 0.84, "", 0.86, "", 0.88, "", "0.90", "", 0.92, "", 0.94, "", 0.96, "", 0.98, "", 1.00) ) + ylab("") + xlab("AUC") + t.plot dev.off() df.wvsc <- cbind(c("an", "pe", "wo"), as.data.frame(matrix(c( sum(l.m.d[[1]][[7]]$ii.an > l.m.w[[1]][[5]]$ii.an, na.rm = T), sum(l.m.d[[1]][[7]]$ii.pe > l.m.w[[1]][[5]]$ii.pe, na.rm = T), sum(l.m.d[[1]][[7]]$ii.wo > l.m.w[[1]][[5]]$ii.wo, na.rm = T), sum(l.m.sc[[3]][[7]]$ii.an > l.m.w[[2]][[5]]$ii.an, na.rm = T), sum(l.m.sc[[3]][[7]]$ii.pe > l.m.w[[2]][[5]]$ii.pe, na.rm = T), sum(l.m.sc[[3]][[7]]$ii.wo > l.m.w[[2]][[5]]$ii.wo, na.rm = T), sum(l.m.uc[[4]][[7]]$ii.an > l.m.w[[3]][[5]]$ii.an, na.rm = T), sum(l.m.uc[[4]][[7]]$ii.pe > l.m.w[[3]][[5]]$ii.pe, na.rm = T), sum(l.m.uc[[4]][[7]]$ii.wo > l.m.w[[3]][[5]]$ii.wo, na.rm = T), sum(l.m.ec[[5]][[7]]$ii.an > l.m.w[[4]][[5]]$ii.an, na.rm = T), sum(l.m.ec[[5]][[7]]$ii.pe > l.m.w[[4]][[5]]$ii.pe, na.rm = T), sum(l.m.ec[[5]][[7]]$ii.wo > l.m.w[[4]][[5]]$ii.wo, na.rm = T) ), ncol = 4) / 8000)) names(df.wvsc) <- c("life2", "d", "s.c", "u.c", "e.c") write.table(df.wvsc, "040-WvsNoWComps.csv", sep = ",", na = "", row.names = F) df.wvsc.meds <- data.frame( d.g = c(median(df.use$bias.stat.d, na.rm = T), median(df.use.w$bias.stat.d, na.rm = T)), d.c = c(median(df.use$d, na.rm = T), median(df.use.w$d, na.rm = T)), s.g = c(median(df.use$bias.s50, na.rm = T), median(df.use.w$bias.s50, na.rm = T)), s.c = c(median(df.use$s.carea, na.rm = T), median(df.use.w$s.carea, na.rm = T)), u.g = c(median(df.use$bias.u50, na.rm = T), median(df.use.w$bias.u50, na.rm = T)), u.c = c(median(df.use$u.carea, na.rm = T), median(df.use.w$u.carea, na.rm = T)), e.g = c(median(df.use$bias.e50, na.rm = T), median(df.use.w$bias.e50, na.rm = T)), e.c = c(median(df.use$e.carea, na.rm = T), median(df.use.w$e.carea, na.rm = T)), rdiff = c(median(df.use$bias.rangediff, na.rm = T), median(df.use.w$bias.rangediff, na.rm = T)), n = c(median(df.use$bias.n.50a, na.rm = T), median(df.use.w$bias.n.50a, na.rm = T)), auc.i = c(median(df.evals[df.evals$reg < 7, "m.auc.i"], na.rm = T), median(df.evals[df.evals$reg == 7, "m.auc.i"], na.rm = T)), auc.n = c(median(df.evals[df.evals$reg < 7, "m.auc.n"], na.rm = T), median(df.evals[df.evals$reg == 7, "m.auc.n"], na.rm = T)) ) write.table(df.wvsc.meds, "040-WvsNoWCompsMedians.csv", sep = ",", na = "", row.names = F) mx.combos <- as.matrix(expand.grid(1:6, 1:6)) ch.reg <- c("AF", "AZ", "EA", "IP", "NA", "SA") load("040-MXDiffStats.RData") df.diffstat <- as.data.frame(mx.diffstat[, 1:7]) load("040-MXDiffStats2.RData") df.diffstat2 <- as.data.frame(mx.diffstat2[, 1:7]) names(df.diffstat) <- c("sp.id", "comp", "reg", "v.cor", "v.cor.null", "sd.null", "p.diff") names(df.diffstat2) <- c("sp.id", "comp", "reg", "v.cor2", "v.cor.null2", "sd.null2", "p.diff2") df.diffstat <- cbind(df.diffstat, df.diffstat2[, c("v.cor2", "v.cor.null2", "sd.null2", "p.diff2")]) df.diffstat <- merge(df.diffstat, df.key, by = "sp.id") df.diffstat$reg <- factor(df.diffstat$reg) df.diffstat$life <- factor(df.diffstat$life) df.diffstat$life2 <- df.diffstat$life df.diffstat[df.diffstat$life2 == "Biennial", "life2"] <- "Perennial" df.diffstat$life2 <- factor(df.diffstat$life2) df.diffstat$b.sig <- df.diffstat$b.sig2 <- rep("NO", nrow(df.diffstat)) df.diffstat[df.diffstat$p.diff >= 0.975 | df.diffstat$p.diff <= 0.025, "b.sig"] <- "YES" df.diffstat$b.sig <- factor(df.diffstat$b.sig) df.diffstat[df.diffstat$p.diff2 >= 0.975 | df.diffstat$p.diff2 <= 0.025, "b.sig2"] <- "YES" df.diffstat$b.sig2 <- factor(df.diffstat$b.sig2) df.diffstat$nuse <- sapply(1:nrow(df.diffstat), function(i) { df.diffstat[i, paste(ch.reg[mx.combos[df.diffstat[i, "comp"], 1]], ".nat", sep = "")] + df.diffstat[i, paste(ch.reg[mx.combos[df.diffstat[i, "comp"], 2]], ".inv", sep = "")] }) df.diffstat.use <- subset(df.diffstat, life != "Aquatic" & life != "NA" & life != "Parasitic" & growth != "Aquatic" & growth != "NA" & growth != "Parasitic" & reg %in% c(1:3, 5:9, 11:15, 17:18, 25:27, 29:33, 35:36)) df.use.i <- df.use[df.use$comptype == "NI", c("sp.id", "id", "reg", "type", "bias.stat.d", "bias.rangediff", "bias.s50", "d"),] df.diffstat.use$bias.stat.d <- sapply(1:nrow(df.diffstat.use), function(i) { s.out <- df.use.i[df.use.i$reg == df.diffstat.use[i, "comp"] & df.use.i$sp.id == df.diffstat.use[i, "sp.id"], "bias.stat.d"] if(length(s.out) == 0) s.out <- NA s.out }) df.diffstat.use$bias.rangediff <- sapply(1:nrow(df.diffstat.use), function(i) { s.out <- df.use.i[df.use.i$reg == df.diffstat.use[i, "comp"] & df.use.i$sp.id == df.diffstat.use[i, "sp.id"], "bias.rangediff"] if(length(s.out) == 0) s.out <- NA s.out }) df.diffstat.use$bias.s50 <- sapply(1:nrow(df.diffstat.use), function(i) { s.out <- df.use.i[df.use.i$reg == df.diffstat.use[i, "comp"] & df.use.i$sp.id == df.diffstat.use[i, "sp.id"], "bias.s50"] if(length(s.out) == 0) s.out <- NA s.out }) df.diffstat.use$d <- sapply(1:nrow(df.diffstat.use), function(i) { s.out <- df.use.i[df.use.i$reg == df.diffstat.use[i, "comp"] & df.use.i$sp.id == df.diffstat.use[i, "sp.id"], "d"] if(length(s.out) == 0) s.out <- NA s.out }) df.diffstat.use$diff.dg.dc <- df.diffstat.use$bias.stat.d - df.diffstat.use$d df.diffstat.use$prop.dg.dc <- df.diffstat.use$bias.stat.d / df.diffstat.use$d f.plot.cordif <- function(f.use, st.yuse, v.brk, ch.lab, st.xlab) { ggplot(data = df.diffstat.use, aes(x = f.use(get(st.yuse)), y = v.cor)) + geom_hline(yintercept = 0, linetype = "longdash", size = 0.5, color = "grey") + geom_point(data = df.diffstat.use[df.diffstat.use$b.sig == "NO",], size = 2, color = "grey") + geom_point(data = df.diffstat.use[df.diffstat.use$p.diff >= 0.975,], size = 2, color = "red") + geom_point(data = df.diffstat.use[df.diffstat.use$p.diff <= 0.025,], size = 2, color = "blue") + geom_smooth(color = "black", size = 1) + scale_x_continuous(breaks = f.use(v.brk), labels = ch.lab, expand = c(0, 0)) + xlab(st.xlab) + ylab("Correl. Niche Shift vs. Climate Avail.") + theme( axis.text = element_text(size = 14, color = "black"), axis.title.y = element_text(size = 16), axis.title.x = element_text(size = 16), axis.line.x = element_line(color = "black"), axis.line.y = element_line(color = "black"), axis.ticks = element_line(color = "black"), panel.grid = element_blank(), panel.background = element_rect(fill = "white"), legend.position = "none" ) } f.plot.cordif2 <- function(f.use, st.yuse, v.brk, ch.lab, st.xlab) { ggplot(data = df.diffstat.use, aes(x = f.use(get(st.yuse)), y = v.cor2)) + geom_hline(yintercept = 0, linetype = "longdash", size = 0.5, color = "grey") + geom_point(data = df.diffstat.use[df.diffstat.use$b.sig2 == "NO",], size = 2, color = "grey") + geom_point(data = df.diffstat.use[df.diffstat.use$p.diff2 >= 0.975,], size = 2, color = "red") + geom_point(data = df.diffstat.use[df.diffstat.use$p.diff2 <= 0.025,], size = 2, color = "blue") + geom_smooth(color = "black", size = 1) + scale_x_continuous(breaks = f.use(v.brk), labels = ch.lab, expand = c(0, 0)) + xlab(st.xlab) + ylab("Correl. Niche Shift Direction vs. Climate Avail.") + theme( axis.text = element_text(size = 14, color = "black"), axis.title.y = element_text(size = 16), axis.title.x = element_text(size = 16), axis.line.x = element_line(color = "black"), axis.line.y = element_line(color = "black"), axis.ticks = element_line(color = "black"), panel.grid = element_blank(), panel.background = element_rect(fill = "white"), legend.position = "none" ) } f.ihs <- function(i) log(i + sqrt(1 + i^2)) f.plot.cordif(identity, "diff.dg.dc", seq(-2, 2, by = 0.1), seq(-2, 2, by = 0.1), "DiffDcDg") f.plot.cordif2(identity, "diff.dg.dc", seq(-2, 2, by = 0.1), seq(-2, 2, by = 0.1), "DiffDcDg") f.plot.cordif(log, "prop.dg.dc", seq(-2, 2, by = 0.1), seq(-2, 2, by = 0.1), "PropDcDg") f.plot.cordif2(log, "prop.dg.dc", seq(-2, 2, by = 0.1), seq(-2, 2, by = 0.1), "PropDcDg") png("040-Fig-Cor.png", width = 9, height = 8, units = "in", res = 320, pointsize = 15) f.plot.cordif(log, "nuse", v.bigb, v.bigl, "Number of Occurrences") dev.off() f.cu <- function(i) sign(i) * abs(i)^(1/3) f.multiplot2 <- function(l.plot, s.cols) { require(grid) s.numplots = length(l.plot) layout <- matrix( seq(1, s.cols * ceiling(s.numplots / s.cols)), ncol = s.cols, nrow = ceiling(s.numplots / s.cols) ) grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) for (i in 1:s.numplots) { matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print( l.plot[[i]], vp = viewport( layout.pos.row = matchidx$row, layout.pos.col = matchidx$col ) ) } } png("040-Fig-Cor-All.png", width = 14, height = 12, units = "in", res = 320, pointsize = 15) f.multiplot2( list( f.plot.cordif(log, "nuse", v.bigb, v.bigl, "Number of Occurrences"), f.plot.cordif(identity, "bias.stat.d", seq(0, 1, by = 0.2), seq(0, 1, by = 0.2), "Range Overlap (Dg)") , f.plot.cordif( f.cu, "bias.rangediff", c(seq(-1000000, -200000, by = 100000), seq(-100000, 100000, by = 10000), seq(200000, 1000000, by = 100000)), c("-1000000", rep("", 8), "-100000", rep("", 9), "0", rep("", 9), "100000", rep("", 8), "1000000"), "Change in Habitat Suitability (dHS)" ), f.plot.cordif(f.cu, "bias.s50", v.bigb, v.bigl, "Range Stability (Sg)") ), 2 ) dev.off() f.multiplot2( list( f.plot.cordif2(log, "nuse", v.bigb, v.bigl, "Number of Occurrences"), f.plot.cordif2(identity, "bias.stat.d", seq(0, 1, by = 0.2), seq(0, 1, by = 0.2), "Range Overlap (Dg)") , f.plot.cordif2( f.cu, "bias.rangediff", c(seq(-1000000, -200000, by = 100000), seq(-100000, 100000, by = 10000), seq(200000, 1000000, by = 100000)), c("-1000000", rep("", 8), "-100000", rep("", 9), "0", rep("", 9), "100000", rep("", 8), "1000000"), "Change in Habitat Suitability (dHS)" ), f.plot.cordif2(f.cu, "bias.s50", v.bigb, v.bigl, "Range Stability (Sg)") ), 2 ) sum(df.diffstat.use$p.diff >= 0.975) / nrow(df.diffstat.use) sum(df.diffstat.use$p.diff <= 0.025) / nrow(df.diffstat.use) df.w <- read.table("025-PointData-AllValid-RPCA.csv", header = T, sep = ",", na.strings = "NA", dec = ".", strip.white = T, fill = T) df.301.5 <- df.w[df.w$sp.id == 301 & df.w$region == 5,] rm(df.w) r.301.3.5 <- raster("pred-bias-301-n3-i5.tif") r.301.5.5 <- raster("pred-bias-301-i5-i5.tif") load("025-Spline-301.Rdata") plot(f.make.rpk(r.301.3.5), col = colorRampPalette(v.rastercol)(100), axes = F, legend = F) points(df.301.5[, c("lon", "lat")], pch = 19, cex = 0.5) f.plot <- function(v.use, v.col) plot( raster(matrix(v.use, ncol = 100, byrow = T)[v.ax[4]:v.ax[3], v.ax[1]:v.ax[2]]), col = colorRampPalette(v.col)(100), axes = F, legend = F ) f.plot.er <- function(v.use, v.col) plot( raster(matrix(v.use, ncol = 100, byrow = T)[v.ax[4]:v.ax[3], v.ax[1]:v.ax[2]]), col = colorRampPalette(v.col)(100), breaks = seq(-1, 1, length.out = 101), axes = F, legend = F ) f.make.rpk <- function(i) i / sum(as.vector(i), na.rm = T) f.scale <- function(i) i / max(abs(as.vector(i)), na.rm = T) v.rastercol <- c("grey", "#0571b0", "#008837", "#ffff40", "#ca0020") v.rgb <- c("red", "grey", "blue") v.ax <- c(11, 60, 1, 40) png("040-Example301.png", width = 15, height = 13.1, units = "in", res = 320, pointsize = 15) par(mfrow = c(3, 2), mar = rep(1, 4)) plot(f.make.rpk(r.301.3.5), col = colorRampPalette(v.rastercol)(100), axes = F, legend = F) points(df.301.5[, c("lon", "lat")], pch = 3, cex = 0.5) f.plot(f.scale(l.katan[[3]][[3]][,121]) * (l.katan[[4]][[5]] > 0), v.rastercol) plot(f.make.rpk(r.301.5.5), col = colorRampPalette(v.rastercol)(100), axes = F, legend = F) points(df.301.5[, c("lon", "lat")], pch = 3, cex = 0.5) f.plot(f.scale(l.katan[[4]][[5]][,121]) * (l.katan[[3]][[3]] > 0), v.rastercol) plot( f.scale(f.make.rpk(r.301.3.5) - f.make.rpk(r.301.5.5)), col = colorRampPalette(v.rgb)(100), breaks = seq(-1, 1, length.out = 101), axes = F, legend = F ) f.plot.er(f.scale(l.katan[[3]][[3]][,121] - l.katan[[4]][[5]][,121]), v.rgb) dev.off()