METHODS AND FEATURES
The ‘sars’ (species–area relationships) package has been programmed using standard S3 methods and is available on CRAN (version 1.2.1) and the development version on GitHub (txm676/sars), meaning researchers can easily add in their own models and functions and integrate these into the multi-model inference framework. Thus, the package represents a resource for future SAR work that can be built on and expanded by workers in the field.
Fitting individual SAR models and a set of SAR models The package provides functions to fit each of the 20 SAR models (see also Triantis, Guilhaumon, and Whittaker (2012)). The ‘sar_models’ function can be used to bring up a list of the 20 model names and a Table can be generated in the package using the ‘display_sars_models’ function. With the exception of the linear model (which is fitted using standard linear regression), all models are fitted using non-linear regression and the model parameters are estimated by minimizing the residual sum of squares with an unconstrained Nelder-Mead optimization algorithm and the ‘optim’ R function.
Choosing starting parameter values for non-linear regression optimisation algorithms is not always straight forward, depending on the data at hand. The starting values for the parameter estimates used as defaults in the package are carefully chosen to avoid numerical problems and to speed up the convergence process. Custom starting values can be provided for any of the 20 models using the ‘start’ argument in the model fit functions. In addition, we also use a grid search / brute force process (see Archontoulis & Miguez, 2014). ‘grid_start’ creates 100 starting values for each parameter within sensible bounds and then creates a large matrix storing all possible combinations of starting values across the n parameters in the model. For example, for a three parameter model this matrix contains one million rows. A random number of these rows are then selected and used as starting parameter values. There are three options for the grid_start argument to control this process. The default (grid_start = “partial”) randomly samples 500 different sets of starting parameter values for each model, adds these to the model’s default starting values and tests all of these. A more comprehensive set of starting parameter estimates can be used (grid_start = “exhaustive”) - this option allows the user to choose the number of starting parameter sets to be tested (using the grid_n argument) and includes a range of additional starting parameter estimates, e.g. very small values and particular values we have found to be useful for individual models. A large ‘grid_n’ will ensure a more extensive search of starting parameter value space, but will increase the time taken to fit the model. More generally, using grid_start = “exhaustive” in combination with a large grid_n can be very time consuming; however, we would recommend it as it makes it more likely that the optimal model fit will be found, particularly for the more complex models. This is particularly true if any of the model fits does not converge, returns a singular gradient at parameter estimates, or the plot of the model fit does not look optimum. The grid start procedure can also be turned off (grid_start = “none”), meaning just the default starting parameter estimates are used. It should be noted that for certain models (e.g. Weibull 3 par.) ‘grid_start’ has been found to occasionally push model fits into strange parameter space (e.g. step curves) and so for these models ‘grid_start’ just reverts back to our provided parameter starting values. However, for these models the ‘start’ argument can still be used. For certain models, there are also checks to ensure that the returned parameters are logical. For example, the optimisation procedure can sometimes return a negative z-value for the asymptotic model, meaning in the term z^A, the base is negative. Thus, if a fitting process returns a negative z-value we return that the model could not be fitted. Remember that, as grid_start has a random component, when grid_start != “none”, you can get slightly different results each time you fit a model or run sar_average().
Each individual model fit returns an object of class ‘sars’, which is a list of 22 elements containing relevant model fit information, such as the model parameter estimates, the fitted values, the residuals, model fit statistics (e.g. AIC, R2), the observed model shape (linear, convex or sigmoidal), whether or not the fit is asymptotic, and convergence information. The returned object can easily be plotted using the ‘plot.sars’ generic function; as this function is based on the base R plotting framework, the plot aesthetics can be edited using standard plotting arguments (see the ‘plot.sars’ documentation in the package). Summary and print generic functions are also provided for class ‘sars’; these functions follow the output of the standard ‘lm’ function in the ‘stats’ R package. Multiple SAR models can be fitted to the same dataset using the ‘sar_multi’ function and the resultant n model fit objects stored together as a ‘fit_collection’ object. This object is a list of class ‘sars’, where each of the n elements contains an individual SAR model fit. Using the ‘plot.sars’ generic function on a ‘fit_collection’ object generates a grid of the n individual model fit plots (Fig. 1).
#load an example dataset (Preston, 1962), fit the logarithmic SAR model using
#the grid_start method of selecting starting parameter values, return a model
#fit summary and plot the model fit.
data(galap)
fit <- sar_loga(data = galap, grid_start = "partial")
summary(fit) ##
## Model:
## Logarithmic
##
## Call:
## S == c + z * log(A)
##
## Did the model converge: TRUE
##
## Residuals:
## 0% 25% 50% 75% 100%
## -115.68753 -38.74990 -11.36383 25.67070 163.95994
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|) 2.5% 97.5%
## c 0.2915341 34.5985767 0.0084262 0.9933958 -68.9056193 69.489
## z 30.2802630 8.3203931 3.6392827 0.0026812 13.6394768 46.921
##
## R-squared: 0.49, Adjusted R-squared: 0.41
## AIC: 187.19, AICc: 189.19, BIC: 189.51
## Observed shape: convex up, Asymptote: FALSE
##
##
## Warning: The fitted values of the model contain negative values (i.e. negative species richness values)
#Create a fit_collection object containing multiple SAR model fits, and
#plot all fits.
fitC <- sar_multi(data = galap, obj = c("power", "loga", "monod"))##
## Now attempting to fit the 3 SAR models:
##
## -- multi_sars ---------------------------------------------- multi-model SAR --
## > power : v
## > loga : v
## > monod : v
Model fit validation
Model fits can be evaluated through tests of the normality and homoscedasticity of the residuals. In certain cases this is advised, given that the maximum likelihood parameter estimates and least-squares estimates are equivalent only under the assumption of normal errors with constant variance (see the Potential issues section, below). Any of three tests can be selected to test the normality of the residuals: 1) the Lilliefors extension of the Kolmogorov normality test (normaTest = “lillie”), 2) the Shapiro-Wilk test of normality (to be preferred when sample size is small; “shapiro”), and 3) the Kolmogorov-Smirnov test (“kolmo”). As a default, the option to omit a residuals normality test (“none”) is used - thus it is up to the user to select a test if appropriate. Three options are provided to check for the homogeneity of the residuals: 1) a correlation of the squared residuals with the model fitted values (“cor.fitted”), 2) a correlation of the squared residuals with the area values (“cor.area”), or 3) no homogeneity test (“none”; the default). If a test is selected and is significant at the 5% level a warning is provided in the model summary; alternatively, the full results of the three tests can be accessed in the model fit output. A third model validation check for negative predicted richness values (i.e. when at least one of the fitted values is negative) can be undertaken and a warning is provided in the model summary if negative values are predicted. It should be noted that in earlier versions of the package (pre Version 1.3.1) there was a bug in the test checking for homogeneity of the residual variance - the correlation used the raw residuals rather than the squared residuals. This has now been corrected. Finally, the convergence code from the optim algorithm can be checked manually. Occasionally the optim algorithm may not fully converge (e.g. due to degeneracy of the Nelder–Mead simplex) but still return a fitted model. As such, for each model fit we provide the optim model convergence code (anything other than 0 can indicate an issue) and a logical value (verge) specifying whether or not the code was 0.
#load an example dataset, fit the linear SAR model whilst running residual
#normality and homogeneity tests, and return the results of the residual
#normality test
data(galap)
fit <- sar_linear(data = galap, normaTest ="lillie", homoTest = "cor.fitted")
summary(fit) #a warning is provided indicating the normality test failed ##
## Model:
## Linear model
##
## Call:
## S == c + m*A
##
## Did the model converge: TRUE
##
## Residuals:
## 0% 25% 50% 75% 100%
## -107.94802 -52.11499 -24.16226 31.72842 217.95711
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|) 2.5 % 97.5 %
## c 70.314003 25.743130 2.731370 0.016228 15.100481 125.5275
## m 0.185382 0.072884 2.543504 0.023410 0.029060 0.3417
##
## R-squared: 0.32, Adjusted R-squared: 0.21
## AIC: 191.77, AICc: 193.77, BIC: 194.08
## Observed shape: linear, Asymptote: FALSE
##
##
## Warning: The normality test selected indicated the model residuals are not normally distributed (i.e. P < 0.05)
##
## Warning: The homogeneity test selected indicated a signifiant correlation between the residuals and the fitted values (i.e. P < 0.05)## $test
## [1] "lillie"
##
## [[2]]
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: res
## D = 0.2146, p-value = 0.04725Observed model shape and identifying an asymptote
Whilst each of the 20 models has a general shape (Triantis, Guilhaumon, and Whittaker (2012)), the actual observed shape of the model fit can be different, for some models, depending on the parameter estimates. This is important as the shape of the curve has significant implications for conservation applications and the testing of macroecological theory (Rosenzweig (1995)). In ‘sars’, the observed shape of a model fit is determined using the sequential algorithm outlined in Triantis, Guilhaumon, and Whittaker (2012). The shape is calculated using the model fit within the observed range of area values. Briefly, the algorithm works by first determining whether the fit is a straight line. Then, if the fit is classified as not being linear, the observed shape is classified as either convex (upwards or downwards) or sigmoidal by analysis of the second derivative (with respect to area) of the model fit (the full algorithm is detailed in Triantis, Guilhaumon, and Whittaker (2012), p. 220). It should be noted that a small number of models (e.g. Extended Power 2) are listed in previous papers (e.g. Tjørve (2009)) as being sigmoidal without an upper asymptote (and we have followed these classifications). However, a recent study (Godeau et al. (2020)) has argued that sigmoidal models must have an upper asymptote by definition; thus, the classification of these models is currently uncertain.
There has also been considerable debate in the SAR literature as to whether or not the SAR is asymptotic. In the ‘sars’ package, to determine whether a fit is asymptotic, for the relevant models the fitted model parameters are analysed to check whether the estimated asymptote is within the range of the sample data (Triantis, Guilhaumon, and Whittaker (2012)).
Multimodel SAR curve
As well as fitting individual models, the package provides a function (‘sar_average’) to fit up to 20 models, compare the resultant fits using information criteria, and construct a multimodel-averaged SAR curve based on information criteria weight (see Guilhaumon, Mouillot, and Gimenez (2010)). The multimodel average curve is constructed as a linear combination of individual model fits by multiplying the predicted richness values of each of the successfully fitted models by the model’s information criterion weight, and then summing the resultant values across all models (Burnham and Anderson (2002)). Three information criteria are available in the package: AIC, AICc and BIC. Confidence intervals around the multimodel averaged curve can be calculated using a non-parametric bootstrap algorithm described in Guilhaumon, Mouillot, and Gimenez (2010). Briefly, each of the SAR models used in the ‘sar_average’ function is fitted to the data, and the fitted values and residuals stored. The residuals are then transformed using the approach in Davison and Hinkley (1997) (p.259). For each bootstrap sample, an individual model fit is selected with the probability of selection being equal to that model’s information criterion weight. The transformed residuals from this fit are then sampled with replacement and added to the model’s fitted richness values. The ‘sar_average’ function is then used to fit all candidate SAR models to this bootstrapped set of response values, and the multimodel averaged fitted values stored. Percentile confidence intervals are then calculated using all bootstrapped fitted values.
The ‘sar_average’ function can be used without specifying any models, in which case the package attempts to fit each of the 20 models in Table 1; alternatively, a vector of model names or a ‘fit_collection’ object (generated using the ‘sar_multi’ function) can be provided using the ‘obj’ argument. The three model validation tests listed above (normality and homogeneity of residuals, and negative predicted values) can be selected (by default none are selected, so it is up to the user to select any that are appropriate); if any model fails one or more of the tests during the fitting process it is removed from the resultant multimodel SAR curve. The output of the ‘sar_average’ function is a list of class ‘multi’ and class ‘sars’, with two elements. The first element (‘mmi’) contains the multi model inference (fitted values of the multimodel SAR curve), and the second element (‘details’) contains a range of information regarding the fitting process, including the successfully fitted models, the models removed due to failing any of the validation tests, and the information criterion values, delta values and weights for the successfully fitted models. The returned object can easily be plotted using the ‘plot.multi’ generic function, and multiple plot options are available (using the ‘type’ argument; see Fig. 2). The fits of all the successfully fitted models and the multimodel SAR curve (with or without confidence intervals) can be plotted together (‘type’ = “multi”), and a barplot of the information criterion weights of each model can also be produced (‘type’ = “bar”).
#load an example dataset (Niering, 1963), run the ‘sar_average’ function
#using a vector of model names and with no model validation tests, and
#produce the plots in Figure 2 of the paper
data(niering)
#run the ‘sar_average’ function using a vector of model names, and with simply
#using the default model starting parameter estimates (grid_start = "none")
fit <- sar_average(data= niering, obj =c("power","loga","koba","logistic","monod",
"negexpo","chapman","weibull3","asymp"),
grid_start = "none", normaTest = "none", homoTest = "none", neg_check = FALSE,
confInt = TRUE, ciN = 50) #a message is provided indicating that one model ##
## Now attempting to fit the 9 SAR models:
##
## -- multi_sars ---------------------------------------------- multi-model SAR --
## > power : v
## > loga : v
## > koba : v
## > logistic : v
## > monod : v
## > negexpo : v
## > chapman : v
## > weibull3 : v
## > asymp : v
##
## No model validation checks selected
##
## 9 remaining models used to construct the multi SAR:
## Power, Logarithmic, Kobayashi, Logistic(Standard), Monod, Negative exponential, Chapman Richards, Cumulative Weibull 3 par., Asymptotic regression
## --------------------------------------------------------------------------------
##
## Calculating sar_multi confidence intervals - this may take some time:## Warning in sar_conf_int(fit = res, n = ciN, crit = IC, obj_all = obj_all, : The following model(s) has been removed from the confidence interval calculations as NAs/Infs are
## present in the transformed residuals: asymp##
|
| | 0%
|
|= | 2%
|
|=== | 4%
|
|==== | 6%
|
|====== | 8%
|
|======= | 10%
|
|======== | 12%
|
|========== | 14%
|
|=========== | 16%
|
|============= | 18%
|
|============== | 20%
|
|=============== | 22%
|
|================= | 24%
|
|================== | 26%
|
|==================== | 28%
|
|===================== | 30%
|
|====================== | 32%
|
|======================== | 34%
|
|========================= | 36%
|
|=========================== | 38%
|
|============================ | 40%
|
|============================= | 42%
|
|=============================== | 44%
|
|================================ | 46%
|
|================================== | 48%
|
|=================================== | 50%
|
|==================================== | 52%
|
|====================================== | 54%
|
|======================================= | 56%
|
|========================================= | 58%
|
|========================================== | 60%
|
|=========================================== | 62%
|
|============================================= | 64%
|
|============================================== | 66%
|
|================================================ | 68%
|
|================================================= | 70%
|
|================================================== | 72%
|
|==================================================== | 74%
|
|===================================================== | 76%
|
|======================================================= | 78%
|
|======================================================== | 80%
|
|========================================================= | 82%
|
|=========================================================== | 84%
|
|============================================================ | 86%
|
|============================================================== | 88%
|
|=============================================================== | 90%
|
|================================================================ | 92%
|
|================================================================== | 94%
|
|=================================================================== | 96%
|
|===================================================================== | 98%
|
|======================================================================| 100%#(asymp) could not be used in the confidence interval calculation
par(mfrow = c(3,1)) #plot all model fits and the multimodel SAR curve as a separate curve on top
plot(fit, ModTitle = "a) Multimodel SAR", mmSep = TRUE)
#plot the multimodel SAR curve (with confidence intervals; see explanation
#in the main text, above) on its own
plot(fit, allCurves = FALSE, ModTitle =
"c) Multimodel SAR with confidence intervals", confInt = TRUE)
#Barplot of the information criterion weights of each model
plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)
One thing to note is how the sar_average() function deals with model non-convergence. There are different types of non-convergence and these are dealt with differently. If the optimisation algorithm fails to return any solution, the model fit is defined as NA and is then removed, and so does not appear in the model summary table or multi-model curve etc. However, the optimisation algorithm (e.g. Nelder-Mead) can also return non-NA model fits but where the solution is potentially non-optimal (e.g. degeneracy of the Nelder–Mead simplex) - these cases are identified by any optim convergence code that is not zero. We have decided not to remove these fits (i.e. they are kept in the model summary table and multimodel curve), but any instances can be checked using the returned ‘details$converged’ vector and then the model fitting re-run without these models, if preferred. Increasing the starting parameters grid search (see above) may also help avoid this issue.
