## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 7, dpi = 100, out.width = "95%" ) ## ----get_ready, results='hide', message=FALSE, warning=FALSE------------------ # Load packages library(nicheR) #library(terra) # Current directory getwd() # Define new directory #setwd("YOUR/DIRECTORY") # modify if setting a new directory ## ----include=FALSE------------------------------------------------------------ # Saving original plotting parameters original_par <- par(no.readonly = TRUE) # Shared margin settings and palettes used throughout this vignette mars <- c(4, 4, 2, 1) marsr <- c(0.5, 0.5, 2, 4) vir_pal <- hcl.colors(100, palette = "Viridis") ## ----data--------------------------------------------------------------------- # Reference niche (nicheR_ellipsoid object) data("ref_ellipse", package = "nicheR") # Reference 3-dimentainal niche data("example_sp_4", package = "nicheR") # Background environmental data (data.frame) data("back_data", package = "nicheR") # Raster layers for geographic predictions (SpatRaster) ma_bios <- terra::rast(system.file("extdata", "ma_bios.tif", package = "nicheR")) ## ----data_check--------------------------------------------------------------- # Check reference niche ref_ellipse # Check background data head(back_data) # Check the raster layers ma_bios ## ----df_predict--------------------------------------------------------------- pred_df <- predict(ref_ellipse, newdata = back_data[, ref_ellipse$var_names]) class(pred_df) colnames(pred_df) head(pred_df) ## ----raster_predict_basic----------------------------------------------------- pred_rast <- predict(ref_ellipse, newdata = ma_bios[[ref_ellipse$var_names]], include_suitability = TRUE, suitability_truncated = TRUE, include_mahalanobis = TRUE, mahalanobis_truncated = TRUE, keep_data = TRUE) pred_rast names(pred_rast) ## ----concept_figure, echo=FALSE, fig.width=7, fig.height=5, out.width="95%"---- pred_ex <- predict(ref_ellipse, newdata = back_data[, ref_ellipse$var_names], include_suitability = TRUE, include_mahalanobis = TRUE, mahalanobis_truncated = TRUE, suitability_truncated = TRUE, verbose = FALSE) levels <- c(0.50, 0.90, 0.99) p <- ref_ellipse$dimensions d2_cuts <- qchisq(levels, df = p) suit_cuts <- exp(-0.5 * d2_cuts) cols <- c("#1b9e77", "#d95f02", "#7570b3") par(mfrow = c(2, 2), mar = c(4, 4, 2.5, 1), cex = 0.85) # Top-left: E-space scatter colored by Mahalanobis distance # Points outside the ellipsoid boundary are clipped intentionally, # showing only the region where D2 has been computed plot_ellipsoid(ref_ellipse, prediction = pred_ex, col_layer = "Mahalanobis_trunc", col_ell = "#e10000", lwd = 2, pch = 16, cex_bg = 0.4, xlab = "Bio1 (Mean Annual Temperature)", ylab = "Bio12 (Annual Precipitation)", main = "Mahalanobis distance", font.main = 1) add_data(as.data.frame(t(ref_ellipse$centroid)), x = "bio_1", y = "bio_12", pts_col = "#e10000", pch = 8, cex = 1.5, lwd = 2) # Top-right: histogram of D2 values with chi-square quantile lines hist(pred_ex$Mahalanobis[!is.na(pred_ex$Mahalanobis)], freq = FALSE, breaks = 60, col = "#a5d6a7", xlab = expression(paste("Mahalanobis Distance (", D^2, ")")), ylab = "Density", main = expression(paste("Distribution of ", D^2))) lines(density(pred_ex$Mahalanobis[!is.na(pred_ex$Mahalanobis)]), col = "#333333", lwd = 2) for (i in seq_along(d2_cuts)) { abline(v = d2_cuts[i], col = cols[i], lwd = 2, lty = 2) } legend("topright", legend = paste0(levels * 100, "% cl"), col = cols, lwd = 2, lty = 2, bty = "n", cex = 0.78) # Bottom-left: E-space scatter colored by suitability plot_ellipsoid(ref_ellipse, prediction = pred_ex, col_layer = "suitability_trunc", col_ell = "#e10000", lwd = 2, pch = 16, cex_bg = 0.4, xlab = "Bio1 (Mean Annual Temperature)", ylab = "Bio12 (Annual Precipitation)", main = "Suitability", font.main = 1) add_data(as.data.frame(t(ref_ellipse$centroid)), x = "bio_1", y = "bio_12", pts_col = "#e10000", pch = 8, cex = 1.5, lwd = 2) # Bottom-right: histogram of suitability values with matching threshold lines hist(pred_ex$suitability[!is.na(pred_ex$suitability)], freq = FALSE, breaks = 60, col = "#a5d6a7", xlab = "Suitability", ylab = "Density", main = "Distribution of suitability", font.main = 1) lines(density(pred_ex$suitability[!is.na(pred_ex$suitability)]), col = "#333333", lwd = 2) for (i in seq_along(suit_cuts)) { abline(v = suit_cuts[i], col = cols[i], lwd = 2, lty = 2) } legend("topright", legend = paste0(levels * 100, "% cl"), col = cols, lwd = 2, lty = 2, bty = "n", cex = 0.78) ## ----df_predict_all----------------------------------------------------------- pred_df_all <- predict(ref_ellipse, newdata = back_data[, ref_ellipse$var_names], include_mahalanobis = TRUE, include_suitability = TRUE, mahalanobis_truncated = TRUE, suitability_truncated = TRUE) colnames(pred_df_all) ## ----outside_example---------------------------------------------------------- # A point outside the ellipsoid outside_idx <- which(!is.na(pred_df_all$Mahalanobis) & is.na(pred_df_all$Mahalanobis_trunc))[1] pred_df_all[outside_idx, ] ## ----raster_predict_all------------------------------------------------------- pred_rast_all <- predict(ref_ellipse, newdata = ma_bios[[ref_ellipse$var_names]], include_mahalanobis = TRUE, include_suitability = TRUE, mahalanobis_truncated = TRUE, suitability_truncated = TRUE) pred_rast_all names(pred_rast_all) ## ----cl_example, eval=FALSE--------------------------------------------------- # # More conservative: only the innermost 90% of the MVN distribution # pred_conservative <- predict(ref_ellipse, # newdata = back_data[, ref_ellipse$var_names], # suitability_truncated = TRUE, # adjust_truncation_level = 0.90) # # # More permissive: the innermost 99% # pred_permissive <- predict(ref_ellipse, # newdata = back_data[, ref_ellipse$var_names], # suitability_truncated = TRUE, # adjust_truncation_level = 0.99) ## ----cl_comparison, echo=FALSE, fig.width=7, fig.height=2.5, out.width="95%"---- cls_compare <- c(0.90, 0.95, 0.99) par(mfrow = c(1, 3), cex = 0.75, mar = c(3.5, 3.5, 2, 0.5)) for (i in seq_along(cls_compare)) { pred_cl <- predict(ref_ellipse, newdata = back_data[, ref_ellipse$var_names], include_mahalanobis = FALSE, include_suitability = FALSE, suitability_truncated = TRUE, adjust_truncation_level = cls_compare[i], verbose = FALSE) plot_ellipsoid(ref_ellipse, xlab = "Bio1", ylab = "Bio12", main = paste0("cl = ", cls_compare[i])) add_data(pred_cl, x = "bio_1", y = "bio_12", pch = ".", pts_col = "grey", cex = 0.5) add_data(data = pred_cl, x = "bio_1", y = "bio_12", col_layer = "suitability_trunc", pch = 16) add_ellipsoid(ref_ellipse, lwd = 2, col_ell = "#e10000") } ## ----maha_espace-------------------------------------------------------------- maha_df <- predict(ref_ellipse, newdata = back_data[, ref_ellipse$var_names], include_mahalanobis = TRUE, include_suitability = FALSE, verbose = FALSE) par(mar = mars) plot_ellipsoid(ref_ellipse, prediction = maha_df, col_layer = "Mahalanobis", col_ell = "#e10000", lwd = 2, pch = 16, cex_bg = 0.5, xlab = "Bio1 (Mean Annual Temperature)", ylab = "Bio12 (Annual Precipitation)", main = "Mahalanobis Distance (D\u00b2) in E-space") add_data(as.data.frame(t(ref_ellipse$centroid)), x = "bio_1", y = "bio_12", pts_col = "#e10000", pch = 8, cex = 1.5, lwd = 2) legend("topright", legend = c("Ellipsoid boundary", "Centroid", "Low D\u00b2", "High D\u00b2"), col = c("#e10000", vir_pal[5], vir_pal[95]), pch = c(NA, 8, 16, 16), lty = c(1, NA, NA, NA), lwd = c(2, NA, NA, NA), cex = 0.75, bty = "n") ## ----suit_espace-------------------------------------------------------------- suit_df <- predict(ref_ellipse, newdata = back_data[, ref_ellipse$var_names], include_mahalanobis = FALSE, include_suitability = TRUE, verbose = FALSE) par(mar = mars) plot_ellipsoid(ref_ellipse, prediction = suit_df, col_layer = "suitability", pal = vir_pal, col_ell = "#e10000", lwd = 2, pch = 16, cex_bg = 0.5, xlab = "Bio1 (Mean Annual Temperature)", ylab = "Bio12 (Annual Precipitation)", main = "Suitability in E-space") add_data(as.data.frame(t(ref_ellipse$centroid)), x = "bio_1", y = "bio_12", pts_col = "#e10000", pch = 8, cex = 1.5, lwd = 2) legend("topright", legend = c("Ellipsoid boundary", "Centroid", "Low suitability", "High suitability"), col = c("#e10000", "#e10000", vir_pal[5], vir_pal[95]), pch = c(NA, 8, 16, 16), lty = c(1, NA, NA, NA), lwd = c(2, NA, NA, NA), cex = 0.75, bty = "n") ## ----trunc_espace, fig.width=7, fig.height=3.5, out.width="95%"--------------- trunc_df <- predict(ref_ellipse, newdata = back_data[, ref_ellipse$var_names], include_mahalanobis = FALSE, include_suitability = FALSE, mahalanobis_truncated = TRUE, suitability_truncated = TRUE, verbose = FALSE) par(mfrow = c(1, 2), cex = 0.7, mar = mars) plot_ellipsoid(ref_ellipse, prediction = trunc_df, col_layer = "Mahalanobis_trunc", col_bg = "#d4d4d4", col_ell = "#e10000", lwd = 2, pch = 16, cex_bg = 0.4, xlab = "Bio1 (Mean Annual Temperature)", ylab = "Bio12 (Annual Precipitation)", main = "Mahalanobis Trunc. (D\u00b2)") add_data(as.data.frame(t(ref_ellipse$centroid)), x = "bio_1", y = "bio_12", pts_col = "#e10000", pch = 8, cex = 1.2, lwd = 2) legend("topright", legend = c("Ellipsoid boundary", "Inside (low D\u00b2)", "Inside (high D\u00b2)", "Outside (NA)"), col = c("#e10000", vir_pal[5], vir_pal[95], "#d4d4d4"), pch = c(NA, 16, 16, 16), lty = c(1, NA, NA, NA), lwd = c(2, NA, NA, NA), cex = 0.65, bty = "n") plot_ellipsoid(ref_ellipse, prediction = trunc_df, col_layer = "suitability_trunc", col_bg = "#d4d4d4", col_ell = "#e10000", lwd = 2, pch = 16, cex_bg = 0.4, xlab = "Bio1 (Mean Annual Temperature)", ylab = "Bio12 (Annual Precipitation)", main = "Suitability Trunc.") add_data(as.data.frame(t(ref_ellipse$centroid)), x = "bio_1", y = "bio_12", pts_col = "#e10000", pch = 8, cex = 1.2, lwd = 2) legend("topright", legend = c("Ellipsoid boundary", "Low suitability", "High suitability", "Outside (zero)"), col = c("#e10000", vir_pal[5], vir_pal[95], "#d4d4d4"), pch = c(NA, 16, 16, 16), lty = c(1, NA, NA, NA), lwd = c(2, NA, NA, NA), cex = 0.65, bty = "n") ## ----binary_espace------------------------------------------------------------ par(mar = mars) plot_ellipsoid(ref_ellipse, prediction = trunc_df, col_layer = "suitability_trunc", pal = c("#d4d4d4", "#0004d5"), col_bg = "#d4d4d4", col_ell = "#e10000", lwd = 2, pch = 16, cex_bg = 0.5, xlab = "Bio1 (Mean Annual Temperature)", ylab = "Bio12 (Annual Precipitation)", main = "Binary Suitability in E-space") legend("topright", legend = c("Ellipsoid boundary", "Suitable (inside)", "Unsuitable (outside)"), col = c("#e10000", "#0004d5", "#d4d4d4"), pch = c(NA, 16, 16), lty = c(1, NA, NA), lwd = c(2, NA, NA), cex = 0.75, bty = "n") ## ----maha_gspace-------------------------------------------------------------- maha_rast <- predict(ref_ellipse, newdata = ma_bios[[ref_ellipse$var_names]], include_mahalanobis = TRUE, include_suitability = FALSE, verbose = FALSE) par(mar = marsr) terra::plot(maha_rast$Mahalanobis, axes = FALSE, box = TRUE, main = "Mahalanobis Distance (D\u00b2) in G-space") ## ----suit_gspace-------------------------------------------------------------- suit_rast <- predict(ref_ellipse, newdata = ma_bios[[ref_ellipse$var_names]], include_mahalanobis = FALSE, include_suitability = TRUE, verbose = FALSE) par(mar = marsr) terra::plot(suit_rast$suitability, axes = FALSE, box = TRUE, main = "Suitability in G-space") ## ----trunc_gspace, fig.width=7, fig.height=3.5, out.width="95%"--------------- trunc_rast <- predict(ref_ellipse, newdata = ma_bios[[ref_ellipse$var_names]], include_mahalanobis = FALSE, include_suitability = FALSE, mahalanobis_truncated = TRUE, suitability_truncated = TRUE, verbose = FALSE) par(mfrow = c(1, 2), cex = 0.7) terra::plot(trunc_rast$Mahalanobis_trunc, axes = FALSE, box = TRUE, mar = marsr, main = "Mahalanobis Trunc. (D\u00b2)") terra::plot(trunc_rast$suitability_trunc, axes = FALSE, box = TRUE, mar = marsr, main = "Suitability Trunc.") ## ----binary_gspace------------------------------------------------------------ binary_rast <- (trunc_rast$suitability_trunc > 0) * 1 par(mar = marsr) terra::plot(binary_rast, axes = FALSE, box = TRUE, col = c("#d4d4d4", "#0004d5"), main = "Potential Distribution (Binary)") ## ----ellipse_3d--------------------------------------------------------------- range_3d <- data.frame(bio_1 = c(27, 35), bio_12 = c(1000, 1500), bio_15 = c(60, 75)) ellipse_3d <- build_ellipsoid(range = range_3d) ellipse_3d ## ----pred_3d------------------------------------------------------------------ suit_3d <- predict(example_sp_4, newdata = ma_bios[[example_sp_4$var_names]], include_mahalanobis = TRUE, include_suitability = TRUE, suitability_truncated = TRUE, verbose = FALSE, keep_data = TRUE) colnames(suit_3d) head(suit_3d) ## ----pairs_3d, fig.width=7, fig.height=4, out.width="95%"--------------------- par(cex = 0.7) plot_ellipsoid_pairs(ellipse_3d, prediction = as.data.frame(suit_3d), col_layer = "suitability_trunc", col_bg = "#d4d4d4", col_ell = "#e10000", lwd = 2, pch = 16, cex_bg = 0.3) ## ----par_reset---------------------------------------------------------------- par(original_par) ## ----------------------------------------------------------------------------- # We draw 1,000 points from the mathematical distribution of the niche virt_bg <- virtual_data(ref_ellipse, n = 1000, truncate = FALSE, effect = "direct", seed = 1) # We score our background points so we can sample from them based on suitability pred_virt <- predict(ref_ellipse, newdata = virt_bg, include_suitability = TRUE, suitability_truncated = TRUE, keep_data = TRUE) head(pred_virt) ## ----------------------------------------------------------------------------- # We draw 1,000 points from the mathematical distribution of the 3D niche virt_3d_bg <- virtual_data(example_sp_4, n = 1000, truncate = FALSE, effect = "direct", seed = 1) # Predict suitability for the 3D virtual background pred_3d_virt <- predict(example_sp_4, newdata = virt_3d_bg, include_suitability = TRUE, suitability_truncated = TRUE, keep_data = TRUE) head(pred_3d_virt) ## ----save_import, eval=FALSE-------------------------------------------------- # temp_ellipse <- file.path(tempdir(), "ref_ellipse.rda") # save_nicheR(ref_ellipse, file = temp_ellipse) # # ref_ellipse_imported <- read_nicheR(temp_ellipse) ## ----save_preds, eval=FALSE--------------------------------------------------- # # 1. Define temporary file paths # temp_df <- file.path(tempdir(), "predictions_df.csv") # temp_rast <- file.path(tempdir(), "predictions_rast.tif") # # temp_3d_df <- file.path(tempdir(), "predictions_3d_df.csv") # temp_3d_rast <- file.path(tempdir(), "predictions_3d_rast.tif") # # temp_virt <- file.path(tempdir(), "predictions_virt.csv") # temp_virt_3d <- file.path(tempdir(), "predictions_virt_3d.csv") # # # # 2. Save predictions to disk # # Standard predictions # write.csv(pred_df, file = temp_df, row.names = FALSE) # terra::writeRaster(pred_rast, filename = temp_rast, overwrite = TRUE) # # # 3D predictions (saved as both a data frame and a raster) # write.csv(as.data.frame(suit_3d, xy = TRUE), file = temp_3d_df, row.names = FALSE) # terra::writeRaster(suit_3d, filename = temp_3d_rast, overwrite = TRUE) # # # Virtual data predictions (these are data frames) # write.csv(pred_virt, file = temp_virt, row.names = FALSE) # write.csv(pred_3d_virt, file = temp_virt_3d, row.names = FALSE) # # # # 3. Import predictions back into R # pred_df_imported <- read.csv(temp_df) # pred_rast_imported <- terra::rast(temp_rast) # # suit_3d_df_imported <- read.csv(temp_3d_df) # suit_3d_rast_imported <- terra::rast(temp_3d_rast) # # pred_virt_imported <- read.csv(temp_virt) # pred_3d_virt_imported <- read.csv(temp_virt_3d)