Clarifying SAR and SpAR documentation, can suppress messages now

Author: kmartinet

R/create_SAR.R

@@ -3,15 +3,21 @@
 #' Use segmented regression to create a species-area relationship (SAR) plot.
 #' The X axis represents log(island area) and the Y axis represents log(number
 #' of species)
+#' 
+#' If the user would prefer to create their own plot of the 
+#' `ssarp::create_SAR()` output, the `aggDF` element of the returned list
+#' includes the raw points from the plot created here. They can be accessed
+#' as demonstrated in the Examples section.
 #' @param occurrences The dataframe output by `ssarp::find_areas()` (or if
 #' using  a custom dataframe, ensure that it has the following columns:
 #' specificEpithet, areas)
 #' @param npsi The maximum number of breakpoints to estimate for model
 #' selection.  Default: 1
 #' @param visualize (boolean) Whether the plot should be displayed when the 
 #' function is called. Default: FALSE
-#' @return A list of 3 including: the summary output, the segmented regression
-#' object, and the aggregated dataframe used to create the plot
+#' @return A list of 4 including: the summary output, the regression
+#' object, the aggregated dataframe used to create the plot, and the AIC scores
+#' used in model selection
 #' @examples
 #' # The GBIF key for the Anolis genus is 8782549
 #' # Read in example dataset filtered from:
@@ -28,7 +34,9 @@
 #'                   npsi = 1,
 #'                   visualize = FALSE)
 #' plot(seg)
-#' summary <- seg[1]
+#' summary <- seg$summary
+#' model_object <- seg$segObj
+#' points <- seg$aggDF
 #' @export
 
 create_SAR <- function(occurrences, npsi = 1, visualize = FALSE) {

R/create_SpAR.R

@@ -3,6 +3,11 @@
 #' Use segmented regression to create a speciation-area relationship plot. The
 #' X axis represents log(island area) and the Y axis represents
 #' log(speciation rate)
+#' 
+#' If the user would prefer to create their own plot of the 
+#' `ssarp::create_SpAR()` output, the `aggDF` element of the returned list
+#' includes the raw points from the plot created here. They can be accessed
+#' as demonstrated in the Examples section.
 #' @param occurrences The dataframe output by one of ssarp's speciation
 #' methods (`ssarp::estimate_BAMM()`, `ssarp::estimate_DR()`,
 #' `ssarp::estimate_MS()`), or if using a custom dataframe, ensure that it
@@ -11,8 +16,9 @@
 #' selection.  Default: 1
 #' @param visualize (boolean) Whether the plot should be displayed when the 
 #' function is called. Default: FALSE
-#' @return A list of 3 including: the summary output, the segmented regression
-#' object, and the aggregated dataframe used to create the plot
+#' @return A list of 4 including: the summary output, the regression
+#' object, the aggregated dataframe used to create the plot, and the AIC scores
+#' used in model selection
 #' @examples
 #' # The GBIF key for the Anolis genus is 8782549
 #' # Read in example dataset filtered from:
@@ -38,7 +44,9 @@
 #'                    npsi = 1,
 #'                    visualize = FALSE)
 #' plot(seg)
-#' summary <- seg[1]
+#' summary <- seg$summary
+#' model_object <- seg$segObj
+#' points <- seg$aggDF
 #' @export
 
 create_SpAR <- function(occurrences, npsi = 1, visualize = FALSE) {

R/find_areas.R

@@ -93,8 +93,10 @@ find_areas <- function(occs, area_custom = NULL,
     minus <- rep(NA, nrow(occs))
     # Loop through dataframe
     for (i in seq_len(nrow(occs))) {
-      if (nrow(occs) == 0) {
-        cli::cli_alert_warning("No data in occurrence record dataframe")
+      if (nrow(occs) == 0) { 
+        if(!getOption("ssarp.silent", FALSE)){
+          cli::cli_alert_warning("No data in occurrence record dataframe")
+        }
         break
       }
       if (
@@ -128,10 +130,14 @@ find_areas <- function(occs, area_custom = NULL,
     # Next, go through the occs dataframe and see if the Third column has a name.
     # If yes, add to the island list. If NA, go to the Second column.
     # If Second column is NA, go to the First column.
-    cli::cli_alert_info("Recording island names...")
+    if(!getOption("ssarp.silent", FALSE)){
+      cli::cli_alert_info("Recording island names...")
+    }
     for (i in seq_len(nrow(occs))) {
       if (nrow(occs) == 0) {
-        cli::cli_alert_warning("No data in occurrence record dataframe")
+        if(!getOption("ssarp.silent", FALSE)){
+          cli::cli_alert_warning("No data in occurrence record dataframe")
+        }
         break
       }
       if (!is.na(occs[i, "Third"])) {
@@ -156,7 +162,9 @@ find_areas <- function(occs, area_custom = NULL,
     }
 
     # Look through the island area file and find the names in uniq_islands list
-    cli::cli_alert_info("Assembling island dictionary...")
+    if(!getOption("ssarp.silent", FALSE)){
+      cli::cli_alert_info("Assembling island dictionary...")
+    }
     # Initialize vector of island names from island area dataset with
     #  "Island" appended
     area_file_append <- paste0(area_file$Name, " Island")
@@ -197,7 +205,9 @@ find_areas <- function(occs, area_custom = NULL,
     }
 
     # Use the dictionary to add the areas to the final dataframe
-    cli::cli_alert_info("Adding areas to final dataframe...")
+    if(!getOption("ssarp.silent", FALSE)){
+      cli::cli_alert_info("Adding areas to final dataframe...")
+    }
     areas <- rep(0, times = nrow(occs))
 
     for (i in seq_len(nrow(occs))) {
@@ -227,8 +237,10 @@ find_areas <- function(occs, area_custom = NULL,
       if(!is.null(names)){
         polygons <- terra::subset(shapefile, shapefile$name %in% names)
       } else {
-        cli::cli_alert_info(
+        if(!getOption("ssarp.silent", FALSE)){
+          cli::cli_alert_info(
           "Using all names in the shapefile, this might extend processing time")
+        }
         # If the user did not input a "names" vector, use 
         #   the full list of polygon names
         # If there are any NAs in shapefile$name, remove them

R/find_land.R

@@ -64,7 +64,9 @@ find_land <- function(occurrences, fillgaps = FALSE) {
   # Check if there is anything in the given occurrences
   # If not, return NULL to stop the function call
   if (nrow(occurrences) == 0) {
-    cli::cli_alert_warning("Occurrence record dataframe has no entries")
+    if(!getOption("ssarp.silent", FALSE)){
+      cli::cli_alert_warning("Occurrence record dataframe has no entries")
+    }
     return(NULL)
   }
 
@@ -119,10 +121,14 @@ find_land <- function(occurrences, fillgaps = FALSE) {
   if (fillgaps == TRUE) {
     # There might still be a lot of NA entries, so use Photon to try to
     #  fill in gaps
-    cli::cli_alert_info("Filling gaps using Photon...")
+    if(!getOption("ssarp.silent", FALSE)){
+      cli::cli_alert_info("Filling gaps using Photon...")
+    }
     for (i in seq_len(nrow(occs))) {
       if (nrow(occs) == 0) {
-        cli::cli_alert_warning("Occurrence record dataframe has no entries")
+        if(!getOption("ssarp.silent", FALSE)){
+          cli::cli_alert_warning("Occurrence record dataframe has no entries")
+        }
         break
       }
       if (is.na(occs[i, "First"])) {

R/get_sources.R

@@ -43,7 +43,9 @@ get_sources <- function(occs) {
 
     return(count_df)
   } else {
-    cli::cli_alert_warning("datasetKey column not found.")
+    if(!getOption("ssarp.silent", FALSE)){
+      cli::cli_alert_warning("datasetKey column not found.")
+    }
     return(NULL)
   }
 }

README.Rmd

@@ -39,6 +39,16 @@ library(devtools)
 install_github("kmartinet/ssarp")
 ```
 
+## Suppressing Messages
+
+To suppress all messages output by *ssarp*, run
+
+``` r
+options(ssarp.silent = TRUE)
+```
+
+before using any of *ssarp*'s functions.
+
 ## Example: Creating a Species-Area Relationship
 
 A species-area relationship (SAR) visualizes the relationship between species richness (the number of species) and the area of the land mass on which the species live. This brief example covers the *ssarp* workflow for creating a SAR, and more detailed explanations of the code and methods can be found [in the Articles on the ssarp pkgdown website](https://kmartinet.github.io/ssarp/index.html).

README.md

@@ -261,23 +261,81 @@ statistical information about the model.
 
 ### Workflow Summary for using data from GBIF to create a species-area relationship plot
 
-1. Use `rgbif` to gather occurrence records, or input your own dataframe of occurrence records.
-2. Use `find_land(occ, fillgaps)` with the dataframe obtained in Step 1 to figure out the names of landmasses using the occurrence record GPS points and the [*maps* R package](https://cran.r-project.org/web/packages/maps/index.html).  Setting the "fillgaps" parameter to `TRUE` will enable the use of [Photon API](https://photon.komoot.io/) to fill in any missing landmass names left by the *maps* R package.
-3. Use `find_areas(occ, area_custom)` with the dataframe obtained in Step 2 to match the landmass names to a dataset that includes names of most islands on the planet and their areas. If the user would like to use a custom island area dataset instead of the built-in one, the "area_custom" parameter can be set to the name of the custom island area dataframe.
-3a. If you'd like to only include occurrence records from islands, you can remove continental records by using `remove_continents(occ)` with the dataframe returned by `find_areas()`
-5. Use `create_SAR(occ, npsi)` with the dataframe obtained in Step 3 to create a species-area relationship plot that reports information important to the associated regression. The "npsi" parameter indicates the maximum number of breakpoints the user would like to compare for model selection. The returned model and plot correspond with the best-fit model.
+1. Use `rgbif` to gather occurrence records, or input your own dataframe of 
+occurrence records.
+2. Use `find_land(occ, fillgaps)` with the dataframe obtained in Step 1 to 
+figure out the names of landmasses using the occurrence record GPS points and 
+the [*maps* R package](https://cran.r-project.org/web/packages/maps/index.html).
+Setting the "fillgaps" parameter to `TRUE` will enable the use of 
+[Photon API](https://photon.komoot.io/) to fill in any missing landmass names 
+left by the *maps* R package.
+3. Use `find_areas(occ, area_custom)` with the dataframe obtained in Step 2 to 
+match the landmass names to a dataset that includes names of most islands on 
+the planet and their areas. If the user would like to use a custom island area 
+dataset instead of the built-in one, the "area_custom" parameter can be set to 
+the name of the custom island area dataframe.
+3a. If you'd like to only include occurrence records from islands, you can 
+remove continental records by using `remove_continents(occ)` with the dataframe
+returned by `find_areas()`
+4. Use `create_SAR(occ, npsi)` with the dataframe obtained in Step 3 to create 
+a species-area relationship plot that reports information important to the 
+associated regression. The "npsi" parameter indicates the maximum number of 
+breakpoints the user would like to compare for model selection. The returned 
+model and plot correspond with the best-fit model.
 
 ### Workflow summary for using data from GBIF and a user-provided phylogenetic tree to create a speciation-area relationship plot
 
-1. Use `rgbif` to gather occurrence records, or input your own dataframe of occurrence records.
-2. Use `find_land(occ, fillgaps)` with the dataframe obtained in Step 1 to figure out the names of landmasses using the occurrence record GPS points and the [*maps* R package](https://cran.r-project.org/web/packages/maps/index.html).  Setting the "fillgaps" parameter to `TRUE` will enable the use of [Photon API](https://photon.komoot.io/) to fill in any missing landmass names left by the *maps* R package.
-3. Use `find_areas(occ, area_custom)` with the dataframe obtained in Step 2 to match the landmass names to a dataset that includes names of most islands on the planet and their areas. If the user would like to use a custom island area dataset instead of the built-in one, the "area_custom" parameter can be set to the name of the custom island area dataframe.
-3a. If you'd like to only include occurrence records from islands, you can remove continental records by using `remove_continents(occ)` with the dataframe returned by `find_areas()`
-4. Use either `estimate_DR(tree, label_type, occ)` or `estimate_MS(tree, label_type, occ)` with your own phylogenetic tree that corresponds with the taxa signified in previous steps, a classifier that describes your tip labels (whether the tip labels are simply species epithets or full scientific names), and the dataframe obtained in Step 3 to add tip speciation rates using the DR statistic (Jetz et al. 2012) or the lambda calculation for crown groups from Magallόn and Sanderson (2001) respectively to the occurrence dataframe. The user may also choose to estimate tip speciation rates from a BAMM analysis (Rabosky 2014) by using `estimate_BAMM(label_type, occ, edata)` with a classifier that describes your tip labels (whether the tip labels are simply species epithets or full scientific names), the occurrence record dataframe obtained in Step 3, and a bammdata object generated by reading the event data file from a BAMM analysis with the *BAMMtools* package (Rabosky et al. 2014).
-5. Use `create_SpAR(occ, npsi)` with the dataframe obtained in Step 4 to create a speciation-area relationship plot that reports information important to the associated regression. The "npsi" parameter indicates the maximum number of breakpoints the user would like to compare for model selection. The returned model and plot correspond with the best-fit model.
+1. Use `rgbif` to gather occurrence records, or input your own dataframe of 
+occurrence records.
+2. Use `find_land(occ, fillgaps)` with the dataframe obtained in Step 1 to 
+figure out the names of landmasses using the occurrence record GPS points and 
+the [*maps* R package](https://cran.r-project.org/web/packages/maps/index.html).
+Setting the "fillgaps" parameter to `TRUE` will enable the use of 
+[Photon API](https://photon.komoot.io/) to fill in any missing landmass names 
+left by the *maps* R package.
+3. Use `find_areas(occ, area_custom)` with the dataframe obtained in Step 2 to 
+match the landmass names to a dataset that includes names of most islands on 
+the planet and their areas. If the user would like to use a custom island area 
+dataset instead of the built-in one, the "area_custom" parameter can be set to 
+the name of the custom island area dataframe.
+3a. If you'd like to only include occurrence records from islands, you can 
+remove continental records by using `remove_continents(occ)` with the dataframe 
+returned by `find_areas()`
+4. Use either `estimate_DR(tree, label_type, occ)` or 
+`estimate_MS(tree, label_type, occ)` with your own phylogenetic tree that 
+corresponds with the taxa signified in previous steps, a classifier that 
+describes your tip labels (whether the tip labels are simply species epithets 
+or full scientific names), and the dataframe obtained in Step 3 to add tip 
+speciation rates using the DR statistic (Jetz et al. 2012) or the lambda 
+calculation for crown groups from Magallόn and Sanderson (2001) respectively 
+to the occurrence dataframe. The user may also choose to estimate tip 
+speciation rates from a BAMM analysis (Rabosky 2014) by using 
+`estimate_BAMM(label_type, occ, edata)` with a classifier that describes your 
+tip labels (whether the tip labels are simply species epithets or full 
+scientific names), the occurrence record dataframe obtained in Step 3, and a 
+bammdata object generated by reading the event data file from a BAMM analysis 
+with the *BAMMtools* package (Rabosky et al. 2014).
+5. Use `create_SpAR(occ, npsi)` with the dataframe obtained in Step 4 to create 
+a speciation-area relationship plot that reports information important to the 
+associated regression. The "npsi" parameter indicates the maximum number of 
+breakpoints the user would like to compare for model selection. The returned 
+model and plot correspond with the best-fit model.
 
 ### Some helpful notes about well-known text (WKT) representation of geometry
-When running `getData()`, the user can specify a well-known text (WKT) representation of geometry to restrict the geographic location of the returned occurrence records. The rgbif::occ_search function that `getData()` calls requires a counter-clockwise winding order for WKT. I find it helpful to think about WKT polygons in this way: imagine a square around your geographic area of interest and pick one of the corners as a starting point. The order of points in WKT format should follow counter-clockwise from the corner you picked first, and the final entry in the WKT string needs to be the same as the first entry. Additionally, while GPS points are typically represented in "latitude, longitude" format, WKT expects them in "longitude latitude" format with commas separating the points rather than individual longitude and latitude values. WKT polygons can have more specified points than included in this simple square example, and even include polygons nested within others or polygons with holes in the middle. 
+When running `getData()`, the user can specify a well-known text (WKT) 
+representation of geometry to restrict the geographic location of the returned
+occurrence records. The rgbif::occ_search function that `getData()` calls 
+requires a counter-clockwise winding order for WKT. I find it helpful to think 
+about WKT polygons in this way: imagine a square around your geographic area of 
+interest and pick one of the corners as a starting point. The order of points 
+in WKT format should follow counter-clockwise from the corner you picked first,
+and the final entry in the WKT string needs to be the same as the first entry. 
+Additionally, while GPS points are typically represented in "latitude, 
+longitude" format, WKT expects them in "longitude latitude" format with commas 
+separating the points rather than individual longitude and latitude values. 
+WKT polygons can have more specified points than included in this simple square
+example, and even include polygons nested within others or polygons with holes
+in the middle. 
 
 #### Literature Cited