User Guide
Marlon E. Cobos, Hannah L. Owens, and A. Townsend Peterson
Some of the main functions of nichevol use data that need to be loaded from a local directory and others produce results that need to be written in a local directory. Loading data from a local directory and writing the results outside the R environment helps to avoid problems related to RAM limitations. The figure below shows how the directory needs to be structured when working on a project considering multiple variables and species.
Figure 1. Suggested structure of data in a working directory for starting analyses, when more than one environmental variable needs to be considered. An example of how each file should look like is presented.
The code below helps to set a working directory before starting analyses:
directory <- "DRIVE:/YOUR/DIRECTORY" # change the characters accordingly
setwd(directory) The stable version of nichevol is in CRAN and it can be installed using the code below:
install.packages("nichevol")The most recent version of nichevol is available at a GitHub repository and can be installed using the code below. Please, have in mind that updates will be done on this version continuously, which may introduce problems and instabilities.
Note: Try the code below first… If you have any problem during the installation, restart your R session, close other sessions you may have open, and try again. If during the installation you are asked to update packages, please do it. If any of the packages gives an error, please install it alone using install.packages(), then try re-installing nichevol again. Also, it may be a good idea to update R and RStudio (if you are using it).
# Installing and loading packages
if (!require(devtools)) {
install.packages("devtools")
}
if (!require(nichevol)) {
devtools::install_github("marlonecobos/nichevol")
}Once nichevol is installed you can load the package with the following line.
library(nichevol)The following code helps in getting the variable Annual Mean Temperature from the WorldClim database. The example data for species occurrences and spatial polygons representing their accessible areas will be loaded here. If you have your own data in your working directory you may want to load them using raster::raster() for a variable, rgdal::readOGR() for spatial polygons, ape::read.nexus() for phylogenetic trees, and read.csv() for occurrences. Accessible areas and occurrence records for multiple species need to be put in two distinct lists (see examples below).
# example data
## list of species records
data("occ_list", package = "nichevol")
## list of species accessible areas
m_files <- list.files(system.file("extdata", package = "nichevol"),
pattern = "m\\d.gpkg", full.names = TRUE)
m_list <- lapply(m_files, terra::vect)
## raster variable
temp <- rast(system.file("extdata", "temp.tif", package = "nichevol"))
## a simple phylogenetic tree for demonstrations
data("tree", package = "nichevol")Knowing what the data looks like is an important first step before running further analyses.
First, let's see some descriptive statistics of environmental values in the accessible areas and in the species occurrences.
# learn how to use the function
help(stats_eval)
# running stats
stat <- stats_eval(stats = c("mean", "sd", "median", "range"),
Ms = m_list, occurrences = occ_list, species = "species",
longitude = "x", latitude = "y", variable = temp)
# see stats
stat
#> $M_stats
#> Species mean sd median range1 range2
#> 1 RD 9830 218.4736 55.10421 240 14 263
#> 2 RD 3351 244.6283 14.77594 249 175 272
#> 3 RD 6933 152.9771 76.81137 168 -1 269
#> 4 RD 761 222.2735 52.64984 238 -27 263
#> 5 RD 6773 245.2930 10.68543 248 160 269
#> 6 RD 7516 211.0188 78.19016 256 -31 288
#>
#> $Occurrence_stats
#> Species mean sd median range1 range2
#> 1 RD 9830 255.5667 4.279523 255 245 263
#> 2 RD 3351 254.9320 4.120879 254 244 266
#> 3 RD 6933 259.9392 5.271185 260 243 269
#> 4 RD 761 255.4200 4.279293 255 246 263
#> 5 RD 6773 250.1370 5.173714 250 240 268
#> 6 RD 7516 259.0946 6.296892 260 241 278Now, let's produce a graphical representation of environmental values in the accessible areas and in the species occurrences. The plots are produced by species.
# learn how to use the function
help(hist_evalues)
# creating the plot for one species
hist_evalues(M = m_list[[1]], occurrences = occ_list[[1]], species = "species",
longitude = "x", latitude = "y", variable = temp,
CL_lines = c(95, 99), col = c("blue", "red"))
The code below helps in creating tables of niche characters considering one environmental variable. Potential values of characters are: 0 = not present, 1 = present, and ? = uncertain.
# learn how to use the function
help(bin_table)
# preparing binned characters of niche for all species
bin_tab <- bin_table(Ms = m_list, occurrences = occ_list, species = "species",
longitude = "x", latitude = "y", variable = temp,
percentage_out = 5, bin_size = 10)
# the table (only 8 columns are shown here)
bin_tab
#> 21 to 30 31 to 40 41 to 50 51 to 60 61 to 70 71 to 80 81 to 90 91 to 100
#> RD 9830 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 3351 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 6933 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 761 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 6773 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 7516 "0" "0" "0" "0" "0" "0" "0" "0" Tow options for ancestral reconstruction of niche characters are currently available, maximum parsimony and maximum likelihood. After reconstructions, a smoothing process is necessary to avoid suitable environments to be discontinuous because of absences or uncertainties present in between the ranges of suitable values. These functions work with one variable at the time, but users can perform several analyses in a loop if needed.
Data for reconstructions
# renaming tree tips
tree$tip.label <- rownames(bin_tab)
# joining tree with data
tree_data <- geiger::treedata(tree, bin_tab)Maximum parsimony reconstructions
# learn how to use the functions
help(bin_par_rec)
help(smooth_rec)
# reconstruction
par_rec_table <- bin_par_rec(tree_data)
# smoothing
s_par_rec_table <- smooth_rec(par_rec_table)
# ancestral reconstruction results (only 8 columns are shown here)
s_par_rec_table
#> 21 to 30 31 to 40 41 to 50 51 to 60 61 to 70 71 to 80 81 to 90 91 to 100
#> RD 9830 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 3351 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 6933 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 761 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 6773 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 7516 "0" "0" "0" "0" "0" "0" "0" "0"
#> 7 "0" "0" "0" "0" "0" "0" "0" "0"
#> 8 "0" "0" "0" "0" "0" "0" "0" "0"
#> 9 "0" "0" "0" "0" "0" "0" "0" "0"
#> 10 "0" "0" "0" "0" "0" "0" "0" "0"
#> 11 "0" "0" "0" "0" "0" "0" "0" "0" Maximum likelihood reconstructions
# learn how to use the functions
help(bin_ml_rec)
# reconstruction
ml_rec_table <- bin_ml_rec(tree_data)
# smoothing
s_ml_rec_table <- smooth_rec(ml_rec_table)
# ancestral reconstruction results (only 8 columns are shown here)
s_ml_rec_table
#> 21 to 30 31 to 40 41 to 50 51 to 60 61 to 70 71 to 80 81 to 90 91 to 100
#> RD 9830 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 3351 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 6933 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 761 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 6773 "0" "0" "0" "0" "0" "0" "0" "0"
#> RD 7516 "0" "0" "0" "0" "0" "0" "0" "0"
#> 7 "0" "0" "0" "0" "0" "0" "0" "0"
#> 8 "0" "0" "0" "0" "0" "0" "0" "0"
#> 9 "0" "0" "0" "0" "0" "0" "0" "0"
#> 10 "0" "0" "0" "0" "0" "0" "0" "0"
#> 11 "0" "0" "0" "0" "0" "0" "0" "0"
#> LogLik NA NA NA NA NA NA NA NA
#> Rates NA NA NA NA NA NA NA NA
#> SE NA NA NA NA NA NA NA NAThe figure to be produced with the following lines of code has two panels: (1) the left panel shows (as bar-plots) the niches of the species and the reconstructed niches of the ancestors; and (2) the right panel shows how these niches have evolved after comparing the ancestral niches aginst the niches of their descendants.
# learn how to use the functions
help(plot.phylo)
help(niche_labels)
help(niche_legend)
# plotting results from maximum parsimony
par(mfrow = c(1, 2))
plot.phylo(tree, label.offset = 0.04)
niche_labels(tree, s_par_rec_table, label_type = "tip_node", height = 0.6, width = 1.3)
niche_legend(position = "topleft", cex = 0.6)
plot.phylo(tree, label.offset = 0.04)
niche_labels(tree, s_par_rec_table, label_type = "tip", height = 0.6, width = 1.3)
nichevol_labels(tree, s_par_rec_table, height = 0.6, width = 1.3)
nichevol_legend(position = "topleft", cex = 0.6)