--- title: "flexsdm: Overview of Modeling functions" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{flexsdm: Overview of Modeling functions} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ## Introduction Species distribution modeling (SDM) has become a standard tool in multiple research areas, including ecology, conservation biology, biogeography, paleobiogeography, and epidemiology. SDM is an area of active theoretical and methodological research. The *flexsdm* package provides users the ability to manipulate and parameterize models in a variety of ways that meet their unique research needs. This flexibility enables users to define their own complete or partial modeling procedure specific for their modeling situations (e.g., number of variables, number of records, different algorithms and ensemble methods, algorithms tuning, etc.). In this vignette, users will learn about the second set of functions in the *flexsdm* package that fall under the "modeling" umbrella. These functions were designed to construct and validate different types of models and can be grouped into fit_* , tune_* , and esm_* family functions. In addition there is a function to perform ensemble modeling. The fit_* functions construct and validate models with default hyper-parameter values. The tune_* functions construct and validate models by searching for the best combination of hyper-parameter values, and esm_ functions can be used for constructing and validating Ensemble of Small Models. Finally, the fit_ensemble() function is for fitting and validating ensemble models. These are the functions for model construction and validation: **fit_* functions family** + fit_gam() Fit and validate Generalized Additive Models + fit_gau() Fit and validate Gaussian Process models + fit_gbm() Fit and validate Generalized Boosted Regression models + fit_glm() Fit and validate Generalized Linear Models + fit_max() Fit and validate Maximum Entropy models + fit_net() Fit and validate Neural Networks models + fit_raf() Fit and validate Random Forest models + fit_svm() Fit and validate Support Vector Machine models **tune_* functions family** + tune_gbm() Fit and validate Generalized Boosted Regression models with exploration of hyper-parameters + tune_max() Fit and validate Maximum Entropy models with exploration of hyper-parameters + tune_net() Fit and validate Neural Networks models with exploration of hyper-parameters + tune_raf() Fit and validate Random Forest models with exploration of hyper-parameters + tune_svm() Fit and validate Support Vector Machine models with exploration of hyper-parameters **model ensemble** + fit_ensemble() Fit and validate ensemble models with different ensemble methods **esm_* functions family** + esm_gam() Fit and validate Generalized Additive Models with Ensemble of Small Model approach + esm_gau() Fit and validate Gaussian Process models Models with Ensemble of Small Model approach + esm_gbm() Fit and validate Generalized Boosted Regression models with Ensemble of Small Model approach + esm_glm() Fit and validate Generalized Linear Models with Ensemble of Small Model approach + esm_max() Fit and validate Maximum Entropy models with Ensemble of Small Model approach + esm_net() Fit and validate Neural Networks models with Ensemble of Small Model approach + esm_svm() Fit and validate Support Vector Machine models with Ensemble of Small Model approach ## Installation First, install the flexsdm package. You can install the released version of *flexsdm* from [github](https://github.com/sjevelazco/flexsdm) with: ``` r # devtools::install_github('sjevelazco/flexsdm') require(flexsdm) require(terra) require(dplyr) ``` ## Project directory setup Decide where on your computer you would like to store the inputs and outputs of your project (this will be your main directory). Use an existing one or use dir.create() to create your main directory. Then specify whether or not to include folders for projections, calibration areas, algorithms, ensembles, and thresholds. For more details see [Vignette 01_pre_modeling](https://sjevelazco.github.io/flexsdm/articles/v01_pre_modeling.html) ## Data, species occurrence and background data In this tutorial, we will be using species occurrences and environmental data that are available through the *flexsdm* package. The "abies" example dataset includes a pr_ab column (presence = 1, and absence = 0), location columns (x, y) and other environmental data. You can load the "abies" data into your local R environment by using the code below: (THIS EXAMPLE LOOKS A LITTLE STRANGE BECAUSE WE ARE ALSO USING BACKGROUND DATA, WHILE THE ABIES DATASET CLEARLY HAS ABSENCES...) ``` r data("abies") data("backg") # Use a subset of the data to keep this tutorial fast to run; # in a real analysis you may want to use the full dataset set.seed(10) abies <- abies %>% dplyr::group_by(pr_ab) %>% dplyr::slice_sample(prop = 0.3) %>% dplyr::ungroup() backg <- backg %>% dplyr::slice_sample(n = 1000) dplyr::glimpse(abies) #> Rows: 420 #> Columns: 13 #> $ id 12040, 10361, 9402, 9815, 10524, 8860, 6431, 11730, 808, 1105… #> $ pr_ab 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0… #> $ x -308908.77, -254286.44, -286978.67, -291848.83, -256658.45, 1… #> $ y 384247.811, 417157.885, 386206.009, 445594.587, 184437.725, -… #> $ aet 572.9367, 259.6567, 587.2900, 443.1700, 355.3867, 354.0600, 4… #> $ cwd 332.0133, 469.4567, 375.9467, 454.9833, 567.6433, 733.3933, 5… #> $ tmin 4.8400, 2.9333, 6.4533, 4.3933, 5.8667, 3.9733, 4.8733, 6.726… #> $ ppt_djf 521.4311, 151.2758, 332.6133, 331.5974, 303.1179, 181.9209, 1… #> $ ppt_jja 48.7567, 15.0839, 15.6589, 19.0647, 10.5549, 9.8277, 7.6569, … #> $ pH 5.631732, 6.202818, 5.500000, 6.000000, 5.200000, 0.000000, 5… #> $ awc 0.10841342, 0.09496477, 0.16000000, 0.07000000, 0.08000000, 0… #> $ depth 63.69832, 68.67968, 178.00000, 107.00000, 61.00000, 0.00000, … #> $ landform 6, 4, 7, 11, 6, 10, 14, 12, 7, 7, 14, 10, 7, 11, 14, 10, 10, … dplyr::glimpse(backg) #> Rows: 1,000 #> Columns: 13 #> $ pr_ab 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, … #> $ x 322239.164, 18219.164, 138909.164, -52520.836, -265550.83… #> $ y -554188.334, -319558.334, -407038.334, 71941.666, 234481.… #> $ aet 594.6933, 196.3567, 349.2067, 504.1167, 511.7433, 450.810… #> $ cwd 564.1600, 1172.8033, 1030.2833, 698.0667, 330.4567, 479.9… #> $ tmin 4.7500, 7.1067, 11.1400, 6.9300, 3.5200, 3.9767, 9.7767, … #> $ ppt_djf 129.3690, 36.1234, 104.8891, 187.1644, 312.0807, 86.2167,… #> $ ppt_jja 12.7164, 1.0201, 2.1754, 6.2898, 11.0941, 22.0211, 1.1523… #> $ pH 5.8000002, 7.8000002, 6.1700625, 5.6999998, 3.0788457, 5.… #> $ awc 0.13000000, 0.18000001, 0.17149687, 0.12000000, 0.0479597… #> $ depth 77.00000, 201.00000, 50.92476, 152.00000, 119.32758, 182.… #> $ percent_clay 17.100000, 25.400000, 22.295300, 21.000000, 13.046448, 41… #> $ landform 6, 13, 7, 14, 7, 7, 10, 13, 6, 11, 10, 10, 6, 10, 13, 11,… ``` If you want to replace the abies dataset with your own data, make sure that your dataset contains the environmental conditions related to presence-absence data. We use a pre-modeling family function for a k-fold partition method (to be used for cross-validation). In the partition method the number of folds or replications must be the same for presence-absence and for background points datasets. ``` r abies2 <- part_random( data = abies, pr_ab = "pr_ab", method = c(method = "kfold", folds = 3) ) dplyr::glimpse(abies2) #> Rows: 420 #> Columns: 14 #> $ id 12040, 10361, 9402, 9815, 10524, 8860, 6431, 11730, 808, 1105… #> $ pr_ab 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0… #> $ x -308908.77, -254286.44, -286978.67, -291848.83, -256658.45, 1… #> $ y 384247.811, 417157.885, 386206.009, 445594.587, 184437.725, -… #> $ aet 572.9367, 259.6567, 587.2900, 443.1700, 355.3867, 354.0600, 4… #> $ cwd 332.0133, 469.4567, 375.9467, 454.9833, 567.6433, 733.3933, 5… #> $ tmin 4.8400, 2.9333, 6.4533, 4.3933, 5.8667, 3.9733, 4.8733, 6.726… #> $ ppt_djf 521.4311, 151.2758, 332.6133, 331.5974, 303.1179, 181.9209, 1… #> $ ppt_jja 48.7567, 15.0839, 15.6589, 19.0647, 10.5549, 9.8277, 7.6569, … #> $ pH 5.631732, 6.202818, 5.500000, 6.000000, 5.200000, 0.000000, 5… #> $ awc 0.10841342, 0.09496477, 0.16000000, 0.07000000, 0.08000000, 0… #> $ depth 63.69832, 68.67968, 178.00000, 107.00000, 61.00000, 0.00000, … #> $ landform 6, 4, 7, 11, 6, 10, 14, 12, 7, 7, 14, 10, 7, 11, 14, 10, 10, … #> $ .part 3, 1, 3, 2, 2, 2, 1, 1, 2, 2, 1, 3, 3, 3, 3, 2, 2, 3, 2, 2, 2… ``` Now, in the abies2 object we have a new column called ".part" with the 3 k-folds (1, 2, 3), indicating which partition each record (row) is a member of. Next, we have to apply the same partition method and number of folds to the environmental conditions of the background points. ``` r backg2 <- part_random( data = backg, pr_ab = "pr_ab", method = c(method = "kfold", folds = 3) ) dplyr::glimpse(backg2) #> Rows: 1,000 #> Columns: 14 #> $ pr_ab 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, … #> $ x 322239.164, 18219.164, 138909.164, -52520.836, -265550.83… #> $ y -554188.334, -319558.334, -407038.334, 71941.666, 234481.… #> $ aet 594.6933, 196.3567, 349.2067, 504.1167, 511.7433, 450.810… #> $ cwd 564.1600, 1172.8033, 1030.2833, 698.0667, 330.4567, 479.9… #> $ tmin 4.7500, 7.1067, 11.1400, 6.9300, 3.5200, 3.9767, 9.7767, … #> $ ppt_djf 129.3690, 36.1234, 104.8891, 187.1644, 312.0807, 86.2167,… #> $ ppt_jja 12.7164, 1.0201, 2.1754, 6.2898, 11.0941, 22.0211, 1.1523… #> $ pH 5.8000002, 7.8000002, 6.1700625, 5.6999998, 3.0788457, 5.… #> $ awc 0.13000000, 0.18000001, 0.17149687, 0.12000000, 0.0479597… #> $ depth 77.00000, 201.00000, 50.92476, 152.00000, 119.32758, 182.… #> $ percent_clay 17.100000, 25.400000, 22.295300, 21.000000, 13.046448, 41… #> $ landform 6, 13, 7, 14, 7, 7, 10, 13, 6, 11, 10, 10, 6, 10, 13, 11,… #> $ .part 2, 2, 3, 1, 2, 3, 2, 2, 1, 3, 1, 1, 3, 1, 1, 3, 2, 2, 2, … ``` In backg2 object we have a new column called ".part" with the 3 k-folds (1, 2, 3). ### 1. Fit and validate models We fit and validate models: I. a maximum entropy model with default hyper-parameter values (flexsdm::fit_max) and II. a random forest model with exploration of hyper-parameters (flexsdm::tune_raf). I. Maximum Entropy models with default hyper-parameter values. ``` r max_t1 <- fit_max( data = abies2, response = "pr_ab", predictors = c("aet", "ppt_jja", "pH", "awc", "depth"), predictors_f = c("landform"), partition = ".part", background = backg2, thr = c("max_sens_spec", "equal_sens_spec", "max_sorensen"), clamp = TRUE, classes = "default", pred_type = "cloglog", regmult = 1 ) #> Formula used for model fitting: #> ~aet + ppt_jja + pH + awc + depth + I(aet^2) + I(ppt_jja^2) + I(pH^2) + I(awc^2) + I(depth^2) + hinge(aet) + hinge(ppt_jja) + hinge(pH) + hinge(awc) + hinge(depth) + ppt_jja:aet + pH:aet + awc:aet + depth:aet + pH:ppt_jja + awc:ppt_jja + depth:ppt_jja + awc:pH + depth:pH + depth:awc + categorical(landform) - 1 #> Replica number: 1/1 #> Partition number: 1/3 #> Partition number: 2/3 #> Partition number: 3/3 ``` This function returns a list object with the following elements: ``` r names(max_t1) #> [1] "model" "predictors" "performance" "performance_part" #> [5] "data_ens" ``` model: A "MaxEnt" class object. This object can be used for predicting. ``` r options(max.print = 20) max_t1$model #> #> Call: glmnet::glmnet(x = mm, y = as.factor(p), family = "binomial", weights = weights, lambda = 10^(seq(4, 0, length.out = 200)) * sum(reg)/length(reg) * sum(p)/sum(weights), standardize = F, penalty.factor = reg) #> #> Df %Dev Lambda #> 1 0 0.00 22.3600 #> 2 0 0.00 21.3500 #> 3 0 0.00 20.3800 #> 4 0 0.00 19.4600 #> 5 0 0.00 18.5800 #> 6 0 0.00 17.7400 #> [ reached 'max' / getOption("max.print") -- omitted 194 rows ] ``` predictors: A tibble with quantitative (c column names) and qualitative (f column names) variables use for modeling. ``` r max_t1$predictors #> # A tibble: 1 × 6 #> c1 c2 c3 c4 c5 f #> #> 1 aet ppt_jja pH awc depth landform ``` performance: The performance metric (see sdm_eval). Those metrics that are threshold dependent are calculated based on the threshold specified in the argument. We can see all the selected threshold values. ``` r max_t1$performance #> # A tibble: 3 × 33 #> model threshold thr_value n_presences n_absences TPR_mean TPR_sd TNR_mean #> #> 1 max equal_sens_s… 0.581 210 210 0.652 0.00825 0.652 #> 2 max max_sens_spec 0.491 210 210 0.767 0.0837 0.586 #> 3 max max_sorensen 0.468 210 210 0.943 0.0378 0.319 #> # ℹ 25 more variables: TNR_sd , W_TPR_TNR_mean , W_TPR_TNR_sd , #> # SORENSEN_mean , SORENSEN_sd , JACCARD_mean , #> # JACCARD_sd , FPB_mean , FPB_sd , OR_mean , OR_sd , #> # TSS_mean , TSS_sd , KAPPA_mean , KAPPA_sd , #> # MCC_mean , MCC_sd , AUC_mean , AUC_sd , #> # BOYCE_mean , BOYCE_sd , CRPS_mean , CRPS_sd , #> # IMAE_mean , IMAE_sd ``` Predicted suitability for each test partition (row) based on the best model. This database is used in fit_ensemble. ``` r max_t1$data_ens #> # A tibble: 420 × 5 #> rnames replicates part pr_ab pred #> #> 1 2 .part 1 0 0.498 #> 2 7 .part 1 0 0.297 #> 3 8 .part 1 0 0.413 #> 4 11 .part 1 0 0.552 #> 5 22 .part 1 0 0.269 #> 6 24 .part 1 0 0.149 #> 7 29 .part 1 0 0.0339 #> 8 35 .part 1 0 0.143 #> 9 36 .part 1 0 0.371 #> 10 37 .part 1 0 0.506 #> # ℹ 410 more rows ``` II- Random forest models with exploration of hyper-parameters. First, we create a data.frame that provides hyper-parameters values to be tested. It Is recommended to generate this data.frame. Hyper-parameter needed for tuning is 'mtry'. The maximum mtry must be equal to total number of predictors. ``` r tune_grid <- expand.grid( mtry = c(2, 4, 6), ntree = c(500, 800) ) ``` We use the same data object abies2, with the same k-fold partition method: ``` r rf_t <- tune_raf( data = abies2, response = "pr_ab", predictors = c( "aet", "cwd", "tmin", "ppt_djf", "ppt_jja", "pH", "awc", "depth" ), predictors_f = c("landform"), partition = ".part", grid = tune_grid, thr = "max_sens_spec", metric = "TSS", ) #> Formula used for model fitting: #> pr_ab ~ aet + cwd + tmin + ppt_djf + ppt_jja + pH + awc + depth + landform #> Tuning model... #> Replica number: 1/1 #> Formula used for model fitting: #> pr_ab ~ aet + cwd + tmin + ppt_djf + ppt_jja + pH + awc + depth + landform #> Replica number: 1/1 #> Partition number: 1/3 #> Partition number: 2/3 #> Partition number: 3/3 ``` Let's see what the output object contains. This function returns a list object with the following elements: ``` r names(rf_t) #> [1] "model" "predictors" "performance" #> [4] "performance_part" "hyper_performance" "data_ens" ``` model: A "randomForest" class object. This object can be used to see the formula details, a basic summary o fthe model, and for predicting. ``` r rf_t$model #> #> Call: #> randomForest(formula = formula1, data = data, mtry = mtry, ntree = ntree, importance = TRUE, ) #> Type of random forest: classification #> Number of trees: 500 #> No. of variables tried at each split: 2 #> #> OOB estimate of error rate: 14.05% #> Confusion matrix: #> 0 1 class.error #> 0 176 34 0.1619048 #> 1 25 185 0.1190476 ``` predictors: A tibble with quantitative (c column names) and qualitative (f column names) variables use for modeling. ``` r rf_t$predictors #> # A tibble: 1 × 9 #> c1 c2 c3 c4 c5 c6 c7 c8 f #> #> 1 aet cwd tmin ppt_djf ppt_jja pH awc depth landform ``` performance: The performance metric (see sdm_eval). Those metrics that are threshold dependent are calculated based on the threshold specified in the argument. We can see all the selected threshold values. ``` r rf_t$performance #> # A tibble: 1 × 35 #> mtry ntree model threshold thr_value n_presences n_absences TPR_mean TPR_sd #> #> 1 2 500 raf max_sens_s… 0.592 210 210 0.876 0.0595 #> # ℹ 26 more variables: TNR_mean , TNR_sd , W_TPR_TNR_mean , #> # W_TPR_TNR_sd , SORENSEN_mean , SORENSEN_sd , #> # JACCARD_mean , JACCARD_sd , FPB_mean , FPB_sd , #> # OR_mean , OR_sd , TSS_mean , TSS_sd , KAPPA_mean , #> # KAPPA_sd , MCC_mean , MCC_sd , AUC_mean , AUC_sd , #> # BOYCE_mean , BOYCE_sd , CRPS_mean , CRPS_sd , #> # IMAE_mean , IMAE_sd ``` Predicted suitability for each test partition (row) based on the best model. This database is used in fit_ensemble. ``` r rf_t$data_ens #> # A tibble: 420 × 5 #> rnames replicates part pr_ab pred #> #> 1 2 .part 1 0 0.734 #> 2 7 .part 1 0 0.172 #> 3 8 .part 1 0 0.1 #> 4 11 .part 1 0 0.642 #> 5 22 .part 1 0 0.028 #> 6 24 .part 1 0 0.516 #> 7 29 .part 1 0 0.03 #> 8 35 .part 1 0 0.516 #> 9 36 .part 1 0 0.1 #> 10 37 .part 1 0 0.352 #> # ℹ 410 more rows ``` These model objects can be used in flexsdm::fit_ensemble(). ### 2. Model Ensemble In this example we fit and validate and ensemble model using the two model objects that were just created. ``` r # Fit and validate ensemble model an_ensemble <- fit_ensemble( models = list(max_t1, rf_t), ens_method = "meansup", thr = NULL, thr_model = "max_sens_spec", metric = "TSS" ) #> | | | 0% | |======================================================================| 100% ``` ``` r # Outputs names(an_ensemble) #> [1] "models" "thr_metric" "predictors" "performance" #> [5] "performance_part" an_ensemble$thr_metric #> [1] "max_sens_spec" "TSS_mean" an_ensemble$predictors #> # A tibble: 2 × 9 #> c1 c2 c3 c4 c5 f c6 c7 c8 #> #> 1 aet ppt_jja pH awc depth landform #> 2 aet cwd tmin ppt_djf ppt_jja landform pH awc depth an_ensemble$performance #> # A tibble: 7 × 33 #> model threshold thr_value n_presences n_absences TPR_mean TPR_sd TNR_mean #> #> 1 meansup equal_sens… 0.538 210 210 0.852 0.00825 0.852 #> 2 meansup lpt 0.124 210 210 1 0 0.519 #> 3 meansup max_fpb 0.548 210 210 0.890 0.0644 0.838 #> 4 meansup max_jaccard 0.548 210 210 0.890 0.0644 0.838 #> 5 meansup max_sens_s… 0.56 210 210 0.876 0.0595 0.852 #> 6 meansup max_sorens… 0.548 210 210 0.890 0.0644 0.838 #> 7 meansup sensitivity 0.466 210 210 0.905 0.00825 0.771 #> # ℹ 25 more variables: TNR_sd , W_TPR_TNR_mean , W_TPR_TNR_sd , #> # SORENSEN_mean , SORENSEN_sd , JACCARD_mean , #> # JACCARD_sd , FPB_mean , FPB_sd , OR_mean , OR_sd , #> # TSS_mean , TSS_sd , KAPPA_mean , KAPPA_sd , #> # MCC_mean , MCC_sd , AUC_mean , AUC_sd , #> # BOYCE_mean , BOYCE_sd , CRPS_mean , CRPS_sd , #> # IMAE_mean , IMAE_sd ``` ### 3. Fit and validate models with Ensemble of Small Model approach This method consists of creating bivariate models with all pair-wise combinations of predictors and perform an ensemble based on the average of suitability weighted by Somers' D metric (D = 2 x (AUC -0.5)). ESM is recommended for modeling species with very few occurrences. This function does not allow categorical variables because the use of these types of variables could be problematic when applied to species with few occurrences. For more detail see Breiner et al. (2015, 2018) ``` r data("abies") library(dplyr) # Create a smaller subset of occurrences set.seed(10) abies2 <- abies %>% na.omit() %>% group_by(pr_ab) %>% dplyr::slice_sample(n = 10) %>% group_by() ``` We can use different methods in the flexsdm::part_random function according to our data. See [part_random](https://sjevelazco.github.io/flexsdm/reference/part_random.html) for more details. ``` r # Using k-fold partition method for model cross validation abies2 <- part_random( data = abies2, pr_ab = "pr_ab", method = c(method = "kfold", folds = 3) ) abies2 #> # A tibble: 20 × 14 #> id pr_ab x y aet cwd tmin ppt_djf ppt_jja pH awc #> #> 1 12787 0 -2.47e5 2.47e5 332. 569. 4.29 243. 10.7 6.72 0.0842 #> 2 7832 0 1.43e5 -1.71e5 138. 748. 1.68 98.6 5.92 0 0 #> 3 10645 0 -3.11e5 3.67e5 604. 332. 6.58 498. 41.3 5.30 0.130 #> 4 6835 0 1.40e5 -1.39e5 169. 385. -3.49 116. 12.3 0 0 #> 5 11397 0 -3.04e5 3.10e5 465. 406. 5.40 290. 15.0 5.80 0.0800 #> 6 9414 0 -2.94e5 3.84e5 487. 492. 5.65 277. 13.8 6.38 0.100 #> 7 1354 0 1.23e5 -9.96e4 169. 355. -3.61 121. 11.6 0 0 #> 8 10529 0 -2.66e5 1.71e5 331. 550. 5.28 354. 11.7 5.11 0.0926 #> 9 9921 0 -2.88e5 4.28e5 393. 534. 5.37 353. 15.8 5.99 0.0743 #> 10 420 0 -1.61e5 2.48e5 402. 757. 7.68 141. 8.51 6.10 0.0700 #> 11 3386 1 2.24e3 5.24e4 339. 401. -0.907 187. 13.4 4.80 0.0943 #> 12 2796 1 -1.58e4 1.33e5 302. 440. 0.143 214. 15.6 2.16 0.0275 #> 13 4392 1 1.26e5 -2.22e5 296. 650. 2.29 180. 9.64 6.10 0.0500 #> 14 2723 1 -1.21e5 2.23e5 494. 372. 2.02 416. 15.5 6 0.0900 #> 15 2896 1 -5.97e4 2.37e5 255. 535. 1.26 167. 9.04 6.22 0.0997 #> 16 479 1 -1.14e5 2.33e5 474. 412. 1.72 323. 15.4 5.70 0.110 #> 17 3570 1 -2.02e3 7.55e4 382. 332. -1.10 186. 15.9 0.0866 0.00117 #> 18 3505 1 -3.89e3 8.75e4 325. 345. -0.533 214. 23.1 6.04 0.0557 #> 19 3631 1 8.59e4 -3.83e4 254. 439. -3.04 134. 9.66 6.20 0.0500 #> 20 2833 1 -5.79e4 2.14e5 294. 523. 1.70 207. 12.5 6.02 0.0869 #> # ℹ 3 more variables: depth , landform , .part ``` This function constructs Generalized Additive Models using the Ensembles of Small Models (ESM) approach (Breiner et al., 2015, 2018). ``` r # We set the model without threshold specification and with the kfold created above esm_gam_t1 <- esm_gam( data = abies2, response = "pr_ab", predictors = c( "aet", "cwd", "tmin", "ppt_djf", "ppt_jja", "pH", "awc", "depth" ), partition = ".part", thr = NULL, k = 2 ) #> | | | 0% | |== | 4% | |===== | 7% | |======== | 11% | |========== | 14% | |============ | 18% | |=============== | 21% | |================== | 25% | |==================== | 29% | |====================== | 32% | |========================= | 36% | |============================ | 39% | |============================== | 43% | |================================ | 46% | |=================================== | 50% | |====================================== | 54% | |======================================== | 57% | |========================================== | 61% | |============================================= | 64% | |================================================ | 68% | |================================================== | 71% | |==================================================== | 75% | |======================================================= | 79% | |========================================================== | 82% | |============================================================ | 86% | |============================================================== | 89% | |================================================================= | 93% | |==================================================================== | 96% | |======================================================================| 100% ``` This function returns a list object with the following elements: ``` r names(esm_gam_t1) #> [1] "esm_model" "predictors" "performance" "performance_part" ``` esm_model: A list with "GAM" class object for each bivariate model. This object can be used for predicting using the ESM approachwith sdm_predict function. ``` r options(max.print = 10) # If you don't want to see printed all the output esm_gam_t1$esm_model #> $`0.37962962962963` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(aet, k = 2) + s(cwd, k = 2) #> #> Estimated degrees of freedom: #> 1.83 1.00 total = 3.83 #> #> UBRE score: 0.3478911 #> #> $`1` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(aet, k = 2) + s(tmin, k = 2) #> #> Estimated degrees of freedom: #> 1.00 1.86 total = 3.86 #> #> UBRE score: -0.6124228 #> #> $`0.726851851851852` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(aet, k = 2) + s(depth, k = 2) #> #> Estimated degrees of freedom: #> 1.90 1.94 total = 4.84 #> #> UBRE score: -0.5136501 #> #> $`0.736111111111111` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(cwd, k = 2) + s(tmin, k = 2) #> #> Estimated degrees of freedom: #> 1.72 1.15 total = 3.87 #> #> UBRE score: -0.6122497 #> #> $`0.263888888888889` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(cwd, k = 2) + s(ppt_djf, k = 2) #> #> Estimated degrees of freedom: #> 1.00 1.36 total = 3.36 #> #> UBRE score: 0.550581 #> #> $`0.0092592592592593` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(cwd, k = 2) + s(ppt_jja, k = 2) #> #> Estimated degrees of freedom: #> 1 1 total = 3 #> #> UBRE score: 0.5427465 #> #> $`0.115740740740741` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(cwd, k = 2) + s(pH, k = 2) #> #> Estimated degrees of freedom: #> 1 1 total = 3 #> #> UBRE score: 0.5553274 #> #> $`0.222222222222222` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(cwd, k = 2) + s(awc, k = 2) #> #> Estimated degrees of freedom: #> 1.00 1.72 total = 3.72 #> #> UBRE score: 0.5269209 #> #> $`0.189814814814815` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(cwd, k = 2) + s(depth, k = 2) #> #> Estimated degrees of freedom: #> 1.00 1.41 total = 3.41 #> #> UBRE score: 0.4707283 #> #> $`1` #> #> Family: binomial #> Link function: logit #> #> Formula: #> pr_ab ~ s(tmin, k = 2) + s(ppt_djf, k = 2) #> #> Estimated degrees of freedom: #> 1.79 1.08 total = 3.87 #> #> UBRE score: -0.6120194 #> #> [ reached 'max' / getOption("max.print") -- omitted 10 entries ] ``` predictors: A tibble with variables use for modeling. ``` r esm_gam_t1$predictors #> # A tibble: 1 × 8 #> c1 c2 c3 c4 c5 c6 c7 c8 #> #> 1 aet cwd tmin ppt_djf ppt_jja pH awc depth ``` performance: Performance metric (see sdm_eval). Those threshold dependent metrics are calculated based on the threshold specified in the argument. ``` r esm_gam_t1$performance #> # A tibble: 7 × 33 #> model threshold thr_value n_presences n_absences TPR_mean TPR_sd TNR_mean #> #> 1 esm_gam equal_sens_… 0.669 10 10 1 0 1 #> 2 esm_gam lpt 0.669 10 10 1 0 1 #> 3 esm_gam max_fpb 0.669 10 10 1 0 1 #> 4 esm_gam max_jaccard 0.669 10 10 1 0 1 #> 5 esm_gam max_sens_sp… 0.669 10 10 1 0 1 #> 6 esm_gam max_sorensen 0.669 10 10 1 0 1 #> 7 esm_gam sensitivity 0.712 10 10 1 0 1 #> # ℹ 25 more variables: TNR_sd , W_TPR_TNR_mean , W_TPR_TNR_sd , #> # SORENSEN_mean , SORENSEN_sd , JACCARD_mean , #> # JACCARD_sd , FPB_mean , FPB_sd , OR_mean , OR_sd , #> # TSS_mean , TSS_sd , KAPPA_mean , KAPPA_sd , #> # MCC_mean , MCC_sd , AUC_mean , AUC_sd , #> # BOYCE_mean , BOYCE_sd , CRPS_mean , CRPS_sd , #> # IMAE_mean , IMAE_sd ``` Now, we test the rep_kfold partition method. In this method 'folds' refers to the number of partitions for data partitioning and 'replicate' refers to the number of replicates. Both assume values >=1. ``` r # Remove the previous k-fold partition abies2 <- abies2 %>% select(-starts_with(".")) # Test with rep_kfold partition using 3 folds and 2 replicates set.seed(10) abies2 <- part_random( data = abies2, pr_ab = "pr_ab", method = c(method = "rep_kfold", folds = 3, replicates = 2) ) abies2 #> # A tibble: 20 × 15 #> id pr_ab x y aet cwd tmin ppt_djf ppt_jja pH awc #> #> 1 12787 0 -2.47e5 2.47e5 332. 569. 4.29 243. 10.7 6.72 0.0842 #> 2 7832 0 1.43e5 -1.71e5 138. 748. 1.68 98.6 5.92 0 0 #> 3 10645 0 -3.11e5 3.67e5 604. 332. 6.58 498. 41.3 5.30 0.130 #> 4 6835 0 1.40e5 -1.39e5 169. 385. -3.49 116. 12.3 0 0 #> 5 11397 0 -3.04e5 3.10e5 465. 406. 5.40 290. 15.0 5.80 0.0800 #> 6 9414 0 -2.94e5 3.84e5 487. 492. 5.65 277. 13.8 6.38 0.100 #> 7 1354 0 1.23e5 -9.96e4 169. 355. -3.61 121. 11.6 0 0 #> 8 10529 0 -2.66e5 1.71e5 331. 550. 5.28 354. 11.7 5.11 0.0926 #> 9 9921 0 -2.88e5 4.28e5 393. 534. 5.37 353. 15.8 5.99 0.0743 #> 10 420 0 -1.61e5 2.48e5 402. 757. 7.68 141. 8.51 6.10 0.0700 #> 11 3386 1 2.24e3 5.24e4 339. 401. -0.907 187. 13.4 4.80 0.0943 #> 12 2796 1 -1.58e4 1.33e5 302. 440. 0.143 214. 15.6 2.16 0.0275 #> 13 4392 1 1.26e5 -2.22e5 296. 650. 2.29 180. 9.64 6.10 0.0500 #> 14 2723 1 -1.21e5 2.23e5 494. 372. 2.02 416. 15.5 6 0.0900 #> 15 2896 1 -5.97e4 2.37e5 255. 535. 1.26 167. 9.04 6.22 0.0997 #> 16 479 1 -1.14e5 2.33e5 474. 412. 1.72 323. 15.4 5.70 0.110 #> 17 3570 1 -2.02e3 7.55e4 382. 332. -1.10 186. 15.9 0.0866 0.00117 #> 18 3505 1 -3.89e3 8.75e4 325. 345. -0.533 214. 23.1 6.04 0.0557 #> 19 3631 1 8.59e4 -3.83e4 254. 439. -3.04 134. 9.66 6.20 0.0500 #> 20 2833 1 -5.79e4 2.14e5 294. 523. 1.70 207. 12.5 6.02 0.0869 #> # ℹ 4 more variables: depth , landform , .part1 , .part2 ``` We use the new rep_kfold partition in the gam model ``` r esm_gam_t2 <- esm_gam( data = abies2, response = "pr_ab", predictors = c( "aet", "cwd", "tmin", "ppt_djf" ), partition = ".part", thr = NULL, k = 2 ) #> | | | 0% | |============ | 17% | |======================= | 33% | |=================================== | 50% | |=============================================== | 67% | |========================================================== | 83% | |======================================================================| 100% ``` Test with random bootstrap partitioning. In method 'replicate' refers to the number of replicates (assumes a value >=1), 'proportion' refers to the proportion of occurrences used for model fitting (assumes a value >0 and <=1). With this method we can configure the proportion of training and testing data according to the species occurrences. In this example, proportion='0.7' indicates that 70% of data will be used for model training, while 30% will be used for model testing. For this method, the function will return .partX columns with "train" or "test" words as the entries. ``` r # Remove the previous k-fold partition abies2 <- abies2 %>% select(-starts_with(".")) # Test with bootstrap partition using 3 replicates set.seed(10) abies2 <- part_random( data = abies2, pr_ab = "pr_ab", method = c(method = "boot", replicates = 3, proportion = 0.7) ) abies2 #> # A tibble: 20 × 16 #> id pr_ab x y aet cwd tmin ppt_djf ppt_jja pH awc #> #> 1 12787 0 -2.47e5 2.47e5 332. 569. 4.29 243. 10.7 6.72 0.0842 #> 2 7832 0 1.43e5 -1.71e5 138. 748. 1.68 98.6 5.92 0 0 #> 3 10645 0 -3.11e5 3.67e5 604. 332. 6.58 498. 41.3 5.30 0.130 #> 4 6835 0 1.40e5 -1.39e5 169. 385. -3.49 116. 12.3 0 0 #> 5 11397 0 -3.04e5 3.10e5 465. 406. 5.40 290. 15.0 5.80 0.0800 #> 6 9414 0 -2.94e5 3.84e5 487. 492. 5.65 277. 13.8 6.38 0.100 #> 7 1354 0 1.23e5 -9.96e4 169. 355. -3.61 121. 11.6 0 0 #> 8 10529 0 -2.66e5 1.71e5 331. 550. 5.28 354. 11.7 5.11 0.0926 #> 9 9921 0 -2.88e5 4.28e5 393. 534. 5.37 353. 15.8 5.99 0.0743 #> 10 420 0 -1.61e5 2.48e5 402. 757. 7.68 141. 8.51 6.10 0.0700 #> 11 3386 1 2.24e3 5.24e4 339. 401. -0.907 187. 13.4 4.80 0.0943 #> 12 2796 1 -1.58e4 1.33e5 302. 440. 0.143 214. 15.6 2.16 0.0275 #> 13 4392 1 1.26e5 -2.22e5 296. 650. 2.29 180. 9.64 6.10 0.0500 #> 14 2723 1 -1.21e5 2.23e5 494. 372. 2.02 416. 15.5 6 0.0900 #> 15 2896 1 -5.97e4 2.37e5 255. 535. 1.26 167. 9.04 6.22 0.0997 #> 16 479 1 -1.14e5 2.33e5 474. 412. 1.72 323. 15.4 5.70 0.110 #> 17 3570 1 -2.02e3 7.55e4 382. 332. -1.10 186. 15.9 0.0866 0.00117 #> 18 3505 1 -3.89e3 8.75e4 325. 345. -0.533 214. 23.1 6.04 0.0557 #> 19 3631 1 8.59e4 -3.83e4 254. 439. -3.04 134. 9.66 6.20 0.0500 #> 20 2833 1 -5.79e4 2.14e5 294. 523. 1.70 207. 12.5 6.02 0.0869 #> # ℹ 5 more variables: depth , landform , .part1 , .part2 , #> # .part3 ``` Use the new rep_kfold partition in the gam model ``` r esm_gam_t3 <- esm_gam( data = abies2, response = "pr_ab", predictors = c( "aet", "cwd", "tmin", "ppt_djf" ), partition = ".part", thr = NULL ) #> | | | 0% | |============ | 17% | |======================= | 33% | |=================================== | 50% | |=============================================== | 67% | |========================================================== | 83% | |======================================================================| 100% ```