trainByCrossValid: Calibrate a distribution/niche model using cross-validation

In adamlilith/enmSdm: Tools for Modeling Niches and Distributions of Species

Calibrate a distribution/niche model using cross-validation

Description

This function is an extension of any of the "trainXYZ" functions for calibrating species distribution and ecological niche models. This function uses the "trainXYZ" function to calibrate and evaluate a suite of models using cross-validation. The models are evaluated against withheld data to determine the optimal settings for a "final" model using all available data.

Usage

trainByCrossValid(
  data,
  resp = names(data)[1],
  preds = names(data)[2:ncol(data)],
  folds = dismo::kfold(data),
  trainFx = enmSdm::trainGlm,
  ...,
  metrics = c("logLoss", "cbi", "auc", "fpb", "tss", "msss", "mdss", "minTrainPres",
    "trainSe95", "trainSe90"),
  weightEvalTrain = TRUE,
  weightEvalTest = TRUE,
  na.rm = FALSE,
  out = c("models", "tuning"),
  verbose = 1
)

Arguments

data	Data frame or matrix. Environmental predictors (and no other fields) for presences and background sites.
resp	Character or integer. Name or column index of response variable. Default is to use the first column in data.
preds	Character list or integer list. Names of columns or column indices of predictors. Default is to use the second and subsequent columns in data.
folds	Either a numeric vector, or matrix or data frame: If a vector, there must be one value per row in data. If there are K unique values in the vector, then K unique models will be trained. Each model will use all of the data except for rows that match a particular value in the folds vector. For example, if folds = c(1, 1, 1, 2, 2, 2, 3, 3, 3), then three models will be trained, one with all rows that match the 2s and 3s, one with all rows matching 1s and 2s, and one will all rows matching 1s and 3s. The models will be evaluated against the withheld data and against the training data. Use NA to exclude rows from all testing/training. The default is to construct 5 folds of roughly equal size. If a matrix or data frame, there must be one row per row in data. Each column corresponds to a different model to be trained. For a given column there should be only two unique values, plus possibly NAs. Of the two values, the lesser value will be used to identify the calibration data and the greater value the evaluation data. Rows with NAs will be ignored. For example, a particular column could contain 1s, 2, and NAs. Data rows corresponding to 1s will be used as training data, data rows corresponding to 2s as test data, and rows with NA are dropped. The NA flag is useful for creating spatially-structured cross-validation folds where training and test sites are separated (spatially) by censored (ignored) data.
trainFx	Function, name of the "trainXYZ" function to use. Currently the functions/algorithms supported are trainGlm, trainBrt, trainMaxEnt, trainMaxNet, and trainNs.
...	Arguments to pass to the "trainXYZ" function.
metrics	Character vector, names of evaluation metrics to calculate. If weightEvalTrain and/or weightEvalTest is TRUE, then the "train" and "test" version of each metric will be weighted versions of each. The default is to calculate all of: 'logLoss': Log loss. Higher (less negative) values imply better fit. 'cbi': Continuous Boyce Index (CBI). Calculated with contBoyce. 'auc': Area under the receiver-operator characteristic curve (AUC). Calculated with aucWeighted. 'tss': Maximum value of the True Skill Statistic. Calculated with tssWeighted. 'msss': Sensitivity and specificity calculated at the threshold that maximizes sensitivity (true presence prediction rate) plus specificity (true absence prediction rate). 'mdss': Sensitivity (se) and specificity (sp) calculated at the threshold that minimizes the difference between sensitivity and specificity. 'minTrainPres': Sensitivity and specificity at the greatest threshold at which all training presences are classified as "present". 'trainSe95' and/or 'trainSe90': Sensitivity at the threshold that ensures either 95 or 90 percent of all training presences are classified as "present" (training sensitivity = 0.95 or 0.9).
weightEvalTrain	Logical, if TRUE (default) and an argument named w is specified in ..., then evaluation statistics that support weighting will use the weights specified by w for the "train" version of evaluation statistics. If FALSE, there will be no weighting of test sites. Note that this applies only to the calculation of evaluation statistics. If w is supplied weights they will be used for model calibration.
weightEvalTest	Logical, if TRUE (default) and an argument named w is specified in ..., then evaluation statistics that support weighting will use the weights specified by w for the "test" version of evaluation statistics. If FALSE, there will be no weighting of test sites. Note that this applies only to the calculation of evaluation statistics. If w is supplied then weights will be used for model calibration.
na.rm	Logical, if TRUE then remove NA predictions before calculating evaluation statistics. If FALSE (default), propagate NAs (meaning if predictions contain NAs, then the evaluation statistic will most likely also be NA.)
out	Character. Indicates type of value returned. If 'models' then returns a list of a list of candidate models (one sublist per fold. If 'tuning' then just return the evaluation table for candidate models of each fold. If both then return a 2-item list with all candidate models and tuning tables. WARNING: Depending on the type of model, using 'models' may produce objects that are very large in memory.
verbose	Numeric. If 0 show no progress updates. If > 0 then show minimal progress updates for this function only. If > 1 show detailed progress for this function. If > 2 show detailed progress plus detailed progress for the "trainXYZ" function.

Details

In some cases models do not converge (e.g., boosted regression trees and generalized additive models sometimes suffer from this issue). In this case the model will be skipped, but a data frame with the k-fold and model number in the fold will be returned in the $meta element in the output. If all models converged, then this data frame will be empty.

Value

A list object with several named elements:

• meta: Meta-data on the model call.

• folds: The folds object.

• models (if 'models' is in argument out): A list of model objects, one per data fold

• tuning (if 'tuning' is in argument out): One data frame per k-fold, each containing evaluation statistics for all candidate models in the fold.

References

Fielding, A.H. and J.F. Bell. 1997. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation 24:38-49.

La Rest, K., Pinaud, D., Monestiez, P., Chadoeuf, J., and Bretagnolle, V. 2014. Spatial leave-one-out cross-validation for variable selection in the presence of spatial autocorrelation. Global Ecology and Biogeography 23:811-820.

Wunderlich, R.F., Lin, P-Y., Anthony, J., and Petway, J.R. 2019. Two alternative evaluation metrics to replace the true skill statistic in the assessment of species distribution models. Nature Conservation 35:97-116.

See Also

trainBrt, trainCrf, trainGam, trainGlm, trainMaxEnt, trainMaxNet, trainLars, trainMaxNet, trainRf, trainNs

Examples

## Not run: 
set.seed(123)
### contrived example
# generate training/testing data
n <- 10000
x1 <- seq(-1, 1, length.out=n) + rnorm(n)
x2 <- seq(10, 0, length.out=n) + rnorm(n)
x3 <- rnorm(n)
y <- 2 * x1 + x1^2 - 10 * x2 - x1 * x2
y <- statisfactory::invLogitAdj(y, 0.001)
presAbs <- as.integer(runif(10000) > (1 - y))
data <- data.frame(presAbs=presAbs, x1=x1, x2=x2, x3=x3)

model <- trainGlm(data, verbose=TRUE)
summary(model) # most parsimonious model

folds <- dismo::kfold(data, 3)
out <- trainByCrossValid(data, folds=folds, verbose=1)

str(out, 1)
summaryByCrossValid(out)

str(out, 1)
head(out$tuning[[1]])
head(out$tuning[[2]])
head(out$tuning[[3]])

# can do following for each fold (3 of them)
lapply(out$models[[1]], coefficients)
sapply(out$models[[1]], logLik)
sapply(out$models[[1]], AIC)

# select model for k = 1 with greatest CBI
top <- which.max(out$tuning[[1]]$cbiTest)
summary(out$models[[1]][[top]])

# in fold k = 1, which models perform well but are not overfit?
plot(out$tuning[[1]]$cbiTrain, out$tuning[[1]]$cbiTest, col='white',
		main='Model Numbers for k = 1')
abline(0, 1, col='red')
numModels <- nrow(out$tuning[[1]])
text(out$tuning[[1]]$cbiTrain, out$tuning[[1]]$cbiTest, labels=1:numModels)
usr <- par('usr')
x <- usr[1] + 0.9 * (usr[4] - usr[3])
y <- usr[3] + 0.1 * (usr[4] - usr[3])
text(x, y, labels='overfit', col='red', xpd=NA)
x <- usr[1] + 0.1 * (usr[4] - usr[3])
y <- usr[3] + 0.9 * (usr[4] - usr[3])
text(x, y, labels='suspicious', col='red', xpd=NA)

# other algorithms
# boosted regression trees (with "fast" set of parameters... not recommended
# for normal use)
brt <- trainByCrossValid(data, folds=folds, verbose=2, trainFx=trainBrt,
	maxTrees=2000, treeComplexity=2, learningRate=c(0.01, 0.001))

# MaxEnt with "fast" set of settings (not recommended for normal use)
mx <- trainByCrossValid(data, folds=folds, verbose=2, trainFx=trainMaxEnt,
	regMult=c(1, 2), classes='lp')


## End(Not run)

adamlilith/enmSdm documentation built on Jan. 6, 2023, 11 a.m.