Package 'spatialRF'

Description

Computes the area under the ROC curve in models with binary responses.

Usage

auc(o, p)

Arguments

Value

Numeric, AUC value.

Examples

if(interactive()){

 out <- auc(
   o = c(0, 0, 1, 1),
   p = c(0.1, 0.6, 0.4, 0.8)
   )

}

Description

Computes the correlation matrix among a set of predictors, orders the correlation matrix according to a user-defined preference order, and removes variables one by one, taking into account the preference order, until the remaining ones are below a given Pearson correlation threshold. Warning: variables in preference.order not in colnames(x), and non-numeric columns are removed silently from x and preference.order. The same happens with rows having NA values (na.omit() is applied). The function issues a warning if zero-variance columns are found.

Usage

auto_cor(
  x = NULL,
  preference.order = NULL,
  cor.threshold = 0.5,
  verbose = TRUE
)

Arguments

Details

Can be chained together with auto_vif() through pipes, see the examples below.

Value

List with three slots:

• cor: correlation matrix of the selected variables.

• selected.variables: character vector with the names of the selected variables.

• selected.variables.df: data frame with the selected variables.

See Also

auto_vif()

Examples

if(interactive()){

 #load data
 data(plant_richness_df)

 #on a data frame
 out <- auto_cor(x = plant_richness_df[, 5:21])

 #getting the correlation matrix
 out$cor

 #getting the names of the selected variables
 out$selected.variables

 #getting the data frame of selected variables
 out$selected.variables.df

 #on the result of auto_vif
 out <- auto_vif(x = plant_richness_df[, 5:21])
 out <- auto_cor(x = out)

 #with pipes
 out <- plant_richness_df[, 5:21] %>%
 auto_vif() %>%
 auto_cor()

}

Description

Selects predictors that are not linear combinations of other predictors by using computing their variance inflation factors (VIF). Allows the user to define an order of preference for the selection of predictors. Warning: variables in preference.order not in colnames(x), and non-numeric columns are removed silently from x and preference.order. The same happens with rows having NA values (na.omit() is applied). The function issues a warning if zero-variance columns are found.

Usage

auto_vif(
  x = NULL,
  preference.order = NULL,
  vif.threshold = 5,
  verbose = TRUE
)

Arguments

Details

This function has two modes of operation:

• 1. When the argument preference.order is NULL, the function removes on each iteration the variable with the highest VIF until all VIF values are lower than vif.threshold.

• 2. When preference.order is provided, the variables are selected by giving them priority according to their order in preference.order. If there are variables not in preference.order, these are selected as in option 1. Once both groups of variables have been processed, all variables are put together and selected by giving priority to the ones in preference.order. This method preserves the variables desired by the user as much as possible.

Can be chained together with auto_cor() through pipes, see the examples below.

Value

List with three slots:

• vif: data frame with the names of the selected variables and their respective VIF scores.

• selected.variables: character vector with the names of the selected variables.

• selected.variables.df: data frame with the selected variables.

See Also

auto_cor()

Examples

if(interactive()){

#loading data
data(plant_richness_df)

#on a data frame
out <- auto_vif(x = plant_richness_df[, 5:21])

#getting out the vif data frame
out$vif

#getting the names of the selected variables
out$selected.variables

#getting the data frame of selected variables
out$selected.variables.df

#on the result of auto_cor
out <- auto_cor(x = plant_richness_df[, 5:21])
out <- auto_vif(x = out)

#with pipes
out <- plant_richness_df[, 5:21] %>%
 auto_cor() %>%
 auto_vif()

}

Description

Defines a Beowulf cluster from the IPs of the machines in the cluster, the number of cores of each machine, and the user name. The returned cluster has to be registered with doParallel::registerDoParallel().

Usage

beowulf_cluster(
  cluster.ips = NULL,
  cluster.cores = NULL,
  cluster.user = Sys.info()[["user"]],
  cluster.port = "11000",
  outfile = NULL
)

Arguments

Value

A list ready to be used as input for the spec argument of the function makeCluster.

Examples

if(interactive()){

beowulf.cluster <- beowulf_cluster(
 cluster.ips = c(
   "10.42.0.1",
   "10.42.0.34",
   "10.42.0.104"
   ),
cluster.cores = c(7, 4, 4),
cluster.user = "blas",
cluster.port = "11000"
)


doParallel::registerDoParallel(cl = beowulf.cluster)

#PARALLELIZED foreach LOOP HERE

parallel::stopCluster(cl = beowulf.cluster)

}

Description

When the data is binary, setting the ranager argument case.weights helps to minimize the issues produced by class imbalance. This function takes a binary response variable and returns a vector of weights populated with the values ⁠1/#zeros⁠ and ⁠1/#ones⁠. It is used internally by the function rf().

Usage

case_weights(data = NULL, dependent.variable.name = NULL)

Arguments

Value

A vector with a length equal to nrow(data) with the respective weights of the cases.

Examples

if(interactive()){

 data <- data.frame(
   response = c(0, 0, 0, 1, 1)
 )

 case_weights(
   data = data,
   dependent.variable.name = "response"
 )

 }

Description

Generates four distance thresholds, from 0 to max(distance.matrix)/2.

Usage

default_distance_thresholds(distance.matrix = NULL)

Arguments

Value

A numeric vector with distance thresholds.

Examples

if(interactive()){

 #loading example distance matrix
 data(distance_matrix)

 #computing set of default distance thresholds
 default_distance_thresholds(distance_matrix)

 }

Description

Distance matrix (in km) among the edges of the American ecoregions described in the plant_richness_df dataset.

Usage

data(distance_matrix)

Format

A numeric matrix with 227 rows and columns.

See Also

plant_richness_df

Description

Generates a double-centered matrix (row and column means are zero) from the weights of a distance matrix (see weights_from_distance_matrix()) and a distance threshold. This is a required step before the computation of Moran's Eigenvector Maps.

Usage

double_center_distance_matrix (
  distance.matrix = NULL,
  distance.threshold = 0
)

Arguments

Value

A double-centered matrix of the same dimensions as x.

See Also

weights_from_distance_matrix(), mem(), mem_multithreshold()

Examples

if(interactive()){

 #loading the distance matrix
 data(distance_matrix)

 x <- double_center_distance_matrix(
   distance.matrix = distance_matrix
 )
 x

 }

Description

Removes spatial predictors that are pair-wise correlated with other spatial predictors (which happens when there are several close distance thresholds), and spatial predictors correlated with non-spatial predictors.

Usage

filter_spatial_predictors(
  data = NULL,
  predictor.variable.names = NULL,
  spatial.predictors.df = NULL,
  cor.threshold = 0.5
)

Arguments

Value

A data frame with non-redundant spatial predictors.

Examples

if(interactive()){

#loading data
data("distance_matrix")
data("plant_richness_df")

#computing Moran's Eigenvector Maps
spatial.predictors.df <- mem_multithreshold(
  distance_matrix = distance_matrix,
  distance.thresholds = c(0, 1000)
  )

#filtering spatial predictors
spatial.predictors.df <- filter_spatial_predictors(
  data = plant_richness_df,
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  spatial.predictors.df = spatial.predictors.df,
  cor.threshold = 0.50
 )


}

Description

Returns performance metrics produced by rf_evaluate().

Usage

get_evaluation(model)

Arguments

Value

A data frame with evaluation scores. The following columns are shown:

• model: Identifies the given model. The values are "Full", (original model introduced into rf_evaluate()), "Training" (model trained on an independent training spatial fold), and "Testing" (predictive performance of the training model on an independent testing spatial fold). The performance values of the "Testing" model represent the model performance on unseen data, and hence its ability to generalize.

• metric: Four values representing different evaluation metrics, "rmse", "nrmse", "r.squared", and "pseudo.r.squared".

• mean, sd, min, and max: Average, standard deviation, minimum, and maximum of each metric across the evaluation (cross-validation) iterations.

See Also

rf_evaluate(), plot_evaluation(), print_evaluation()

Examples

if(interactive()){

#loading data
data(plant_richness_df)
data(distance_matrix)

#fitting a random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1,
  verbose = FALSE
)

#evaluating the model with spatial cross-validation
rf.model <- rf_evaluate(
  model = rf.model,
  xy = plant_richness_df[, c("x", "y")],
  n.cores = 1,
  verbose = FALSE
)

#getting evaluation results from the model
x <- get_evaluation(rf.model)
x

}

Description

Gets variable importance scores from rf(), rf_repeat(), and rf_spatial() models.

Usage

get_importance(model)

Arguments

Value

A data frame with variable names and importance scores.

See Also

rf(), rf_repeat(), rf_spatial(), plot_importance(), print_importance().

Examples

if(interactive()){

data(plant_richness_df)
data(distance_matrix)

rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1,
  verbose = FALSE
)

x <- get_importance(rf.model)
x

}

Description

Gets local importance scores from rf(), rf_repeat(), and rf_spatial() models.

Usage

get_importance_local(model)

Arguments

Value

A data frame with variable names and local importance scores.

See Also

rf(), rf_repeat(), rf_spatial(), plot_importance(), print_importance().

Examples

if(interactive()){

#loading example data
data(plant_richness_df)
data(distance_matrix)

#fittinga random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1,
  verbose = FALSE
)

#getting importance scores
x <- get_importance_local(rf.model)
x

}

Description

Returns the Moran's I test on the residuals of a model produced by rf(), rf_repeat(), or rf_spatial().

Usage

get_moran(model)

Arguments

Value

A data frame with Moran's I test results produced by moran_multithreshold().

See Also

moran(), moran_multithreshold(), plot_moran(), print_moran().

Examples

if(interactive()){

 #loading example data
 data(plant_richness_df)
 data(distance_matrix)

 #fitting a random forest model
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = c(0, 1000, 2000),
   n.cores = 1,
   verbose = FALSE
 )

 #getting Moran's I of the residuals
 x <- get_moran(rf.model)

}

Description

Returns the performance slot of an rf(), rf_repeat(), or rf_spatial() model computed on the out-of-bag data.

Usage

get_performance(model)

Arguments

Value

A data frame with four columns:

• metric Name of the performance metric.

• median Value of the performance metric. Truly a median only if the model is fitted with rf_repeat().

• median_absolute_deviation median absolute deviation (MAD), only if the model is fitted with rf_repeat(), and NA otherwise.

See Also

print_performance()

Examples

if(interactive()){

 #loading example data
 data(plant_richness_df)
 data(distance.matrix)

 #fitting random forest model
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   n.cores = 1,
   verbose = FALSE
 )

 #getting model performance
 x <- get_performance(rf.model)
 x

}

Description

Returns model predictions from a model fitted with rf(), rf_repeat(), or rf_spatial().

Usage

get_predictions(model)

Arguments

Value

A vector with predictions, or median of the predictions across repetitions if the model was fitted with rf_repeat().

Examples

if(interactive()){

#loading example data
data(plant_richness_df)

#fitting a random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  n.cores = 1,
  verbose = FALSE
)

#get vector of predictions
x <- get_predictions(rf.model)
x

}

Description

Returns the residuals of models fitted with rf(), rf_repeat(), or rf_spatial().

Usage

get_residuals(model)

Arguments

Value

A vector with model residuals, or the median of model residuals across repetitions if the model was fitted with rf_repeat().

Examples

if(interactive()){

#load example data
data(plant_richness_df)

#fit random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  n.cores = 1,
  verbose = FALSE
)

#getting vector with residuals
x <- get_residuals(rf.model)
x

}

Description

Generates and returns the data required to plot the response curves of a model fitted with rf(), rf_repeat(), or rf_spatial().

Usage

get_response_curves(
  model = NULL,
  variables = NULL,
  quantiles = c(0.1, 0.5, 0.9),
  grid.resolution = 200,
  verbose = TRUE
)

Arguments

Details

All variables that are not plotted in a particular response curve are set to the values of their respective quantiles, and the response curve for each one of these quantiles is shown in the plot.

Value

A data frame with the following columns:

• response: Predicted values of the response, obtained with stats::predict().

• predictor: Values of the given predictor.

• quantile: Grouping column, values of the quantiles at which the other predictors are set to generate the response curve.

• model: Model number, only relevant if the model was produced with rf_repeat().

• predictor.name: Grouping variable, name of the predictor.

• response.name: Grouping variable, name of the response variable.

See Also

plot_response_curves()

Examples

if(interactive()){

#loading example data
data(plant_richness_df)

#fitting random forest model
out <- rf(
 data = plant_richness_df,
 dependent.variable.name = "richness_species_vascular",
 predictor.variable.names = colnames(plant_richness_df)[5:21],
 n.cores = 1,
 verbose = FALSE
)

#getting data frame with response curves
p <- get_response_curves(out)
head(p)

}

Description

Returns spatial predictors from a model fitted with rf_spatial() in order of importance.

Usage

get_spatial_predictors(model)

Arguments

Value

A data frame with the spatial predictors included in the model.

Examples

if(interactive()){

 #loading example data
 data(distance_matrix)
 data(plant_richness_df)

 #fittind spatial model
 model <- rf_spatial(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = c(0, 1000),
   n.cores = 1,
   method = "mem.moran.sequential"
 )

#getting data frame with the selected spatial predictors
spatial.predictors <- get_spatial_predictors(model)
head(spatial.predictors)

}

Description

Returns TRUE if dependent.variable.name is a binary variable with the values 1 and 0.

Usage

is_binary(data = NULL, dependent.variable.name = NULL)

Arguments

Value

Logical.

Examples

if(interactive()){

 #dummy data frame
 data <- data.frame(
   response = c(0, 0, 0, 1, 1)
 )

 #checking if response is binary
 is_binary(
   data = data,
   dependent.variable.name = "response"
 )

}

Description

Used internally by make_spatial_folds() and rf_evaluate(). Uses the coordinates of a point xy.i to generate two spatially independent data folds from the data frame xy. It does so by growing a rectangular buffer from xy.i until a number of records defined by training.fraction is inside the buffer. The indices of these records are then stored as "training" in the output list. The indices of the remaining records outside of the buffer are stored as "testing". These training and testing records can be then used to evaluate a model on independent data via cross-validation.

Usage

make_spatial_fold(
  data = NULL,
  dependent.variable.name = NULL,
  xy.i = NULL,
  xy = NULL,
  distance.step.x = NULL,
  distance.step.y = NULL,
  training.fraction = 0.8
)

Arguments

Value

A list with two slots named training and testing with the former having the indices of the training records selected from xy, and the latter having the indices of the testing records.

See Also

make_spatial_folds(), rf_evaluate()

Examples

if(interactive()){

 #loading example data
 data(plant_richness_df)

 #getting case coordinates
 xy <- plant_richness_df[, 1:3]
 colnames(xy) <- c("id", "x", "y")

 #building a spatial fold centered in the first pair of coordinates
 out <- make_spatial_fold(
   xy.i = xy[1, ],
   xy = xy,
   training.fraction = 0.6
 )

 #indices of the training and testing folds
 out$training
 out$testing

 #plotting the data
 plot(xy[ c("x", "y")], type = "n", xlab = "", ylab = "")
 #plots training points
 points(xy[out$training, c("x", "y")], col = "red4", pch = 15)
 #plots testing points
 points(xy[out$testing, c("x", "y")], col = "blue4", pch = 15)
 #plots xy.i
 points(xy[1, c("x", "y")], col = "black", pch = 15, cex = 2)

}

Description

Applies make_spatial_fold() to every record in a data frame xy.selected to generate as many spatially independent folds over the dataset xy as rows are in xy.selected.

Usage

make_spatial_folds(
  data = NULL,
  dependent.variable.name = NULL,
  xy.selected = NULL,
  xy = NULL,
  distance.step.x = NULL,
  distance.step.y = NULL,
  training.fraction = 0.75,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Value

A list with as many slots as rows are in xy.selected. Each slot has two slots named training and testing, with the former having the indices of the training records selected from xy, and the latter having the indices of the testing records.

See Also

make_spatial_fold(), rf_evaluate()

Examples

if(interactive()){

 #loading example data
 data(plant_richness_df)

 #getting case coordinates
 xy <- plant_richness_df[, 1:3]
 colnames(xy) <- c("id", "x", "y")

 #thining til 20 cases
 xy.selected <- thinning_til_n(
   xy = xy,
   n = 20
   )

 #making spatial folds centered on these 20 cases
 out <- make_spatial_folds(
   xy.selected = xy.selected,
   xy = xy,
   distance.step = 0.05, #degrees
   training.fraction = 0.6,
   n.cores = 1
 )

 #plotting training and testing folds
 plot(xy[ c("x", "y")], type = "n", xlab = "", ylab = "")
 #plots training points
 points(xy[out[[10]]$training, c("x", "y")], col = "red4", pch = 15)
 #plots testing points
 points(xy[out[[10]]$testing, c("x", "y")], col = "blue4", pch = 15)
 #plots xy.i
 points(xy[10, c("x", "y")], col = "black", pch = 15, cex = 2)

}

Description

Computes the positive Moran's Eigenvector Maps of a distance matrix.

Usage

mem(
  distance.matrix = NULL,
  distance.threshold = 0,
  colnames.prefix = "mem"
)

Arguments

Details

Takes the distance matrix x, double-centers it with double_center_distance_matrix(), applies eigen, and returns eigenvectors with positive normalized eigenvalues (a.k.a Moran's Eigenvector Maps, or MEMs). These MEMs are later used as spatial predictors by rf_spatial().

Value

A data frame with positive Moran's Eigenvector Maps.

See Also

mem_multithreshold(), rf_spatial()

Examples

if(interactive()){

 #loading example distance matrix
 data(distance_matrix)

 #Moran's Eigenvector Maps of the distance matrix
 mem <- mem(x = distance_matrix)

}

Description

Computes Moran's Eigenvector Maps of a distance matrix (using mem()) over different distance thresholds.

Usage

mem_multithreshold(
  distance.matrix = NULL,
  distance.thresholds = NULL,
  max.spatial.predictors = NULL
)

Arguments

Details

The function takes the distance matrix x, computes its weights at difference distance thresholds, double-centers the resulting weight matrices with double_center_distance_matrix(), applies eigen to each double-centered matrix, and returns eigenvectors with positive normalized eigenvalues for different distance thresholds.

Value

A data frame with as many rows as the distance matrix x containing positive Moran's Eigenvector Maps. The data frame columns are named "spatial_predictor_DISTANCE_COLUMN", where DISTANCE is the given distance threshold, and COLUMN is the column index of the given spatial predictor.

Examples

if(interactive()){

 #loading example data
 data(distance_matrix)

 #computing Moran's eigenvector maps for 0, 1000, and 2000 km
 mem.df <- mem_multithreshold(
   distance.matrix = distance_matrix,
   distance.thresholds = c(0, 1000, 2000)
   )
 head(mem.df)

}

Description

Computes the spatial correlation coefficient (Moran's I) of a vector given a distance matrix, and a distance threshold used to define "neighborhood".

Usage

moran(
  x = NULL,
  distance.matrix = NULL,
  distance.threshold = NULL,
  verbose = TRUE
)

Arguments

Details

Inspired in the Moran.I() function of the ape package.

Value

A list with three named slots:

• test: Data frame with observed and expected Moran's I values, p-value, and interpretation.

• plot: Moran's plot of the vector x against the spatial lags of x.

• plot.df: Data used in the Moran's plot.

See Also

moran_multithreshold()

Examples

if(interactive()){

 #loading example data
 data(distance_matrix)
 data(plant_richness)

 #Moran's I of the response variable
 out <- moran(
   x = plant_richness$richness_species_vascular,
   distance.matrix = distance_matrix
   )
 out

}

Description

Applies moran() to different distance thresholds at the same time.

Usage

moran_multithreshold(
  x = NULL,
  distance.matrix = NULL,
  distance.thresholds = NULL,
  verbose = TRUE
)

Arguments

Details

Using different distance thresholds helps to take into account the uncertainty about what "neighborhood" means in ecological systems (1000km in geological time means little, but 100m might be quite a long distance for a tree to disperse seeds over), and allows to explore spatial autocorrelation of model residuals for several minimum-distance criteria at once.

Value

A named list with the slots:

• df: Data frame with the results of moran per distance threshold.

• plot: A plot of Moran's I across distance thresholds.

• max.moran: Maximum value of Moran's I across thresholds.

• max.moran.distance.threshold: Distance threshold with the maximum Moran's I value.

See Also

moran()

Examples

if(interactive()){

 #loading example data
 data(distance_matrix)
 data(plant_richness)

 #computing Moran's I for the response variable at several reference distances
 out <- moran_multithreshold(
   x = plant_richness$richness_species_vascular,
   distance.matrix = distance_matrix,
   distance.thresholds = c(0, 100, 1000, 10000),
   plot = TRUE
   )
 out

}

Description

Shows the size of the objects currently in the R environment. Helps to locate large objects cluttering the R environment and/or causing memory problems during the execution of large workflows.

Usage

objects_size(n = 10)

Arguments

Value

A data frame with the row names indicating the object name, the field 'Type' indicating the object type, 'Size' indicating the object size, and the columns 'Length/Rows' and 'Columns' indicating the object dimensions if applicable.

Examples

if(interactive()){

 #creating dummy objects
 x <- matrix(runif(100), 10, 10)
 y <- matrix(runif(10000), 100, 100)

 #reading their in-memory size
 objects_size()

}

Description

Optimizes the selection of spatial predictors using two different methods: "moran.i", and "p.value".

Usage

optimization_function(
  x = NULL,
  weight.r.squared = NULL,
  weight.penalization.n.predictors = NULL,
  optimization.method = "moran.i"
)

Arguments

Details

The method "moran.i" tries to maximize ⁠1 - Moran's⁠ I while taking into account the R-squared of the model and a penalization on the number of introduced spatial predictors through the expression

(1 - Moran's I) + w1 * r.squared - w2 * penalization

The method "p.value" uses a binary version of the p-values of Moran's I (1 if >= 0.05, 0 otherwise), and uses the expression

max(1 - Moran's I, binary p-value) + w1 * r.squared - w2 * penalization

The "moran.i" method generally selects more spatial predictors than the "p.value" method.

Value

A numeric vector with the optimization criteria.

See Also

select_spatial_predictors_recursive(), select_spatial_predictors_sequential()

Description

Extracts all factors of a principal component analysis of a matrix or data frame. Just a convenient wrapper for prcomp.

Usage

pca(
  x = NULL,
  colnames.prefix = "pca_factor"
)

Arguments

Details

Columns in x with zero variance are removed before computing the PCA.

Value

A data frame with the PCA factors of x.

See Also

pca_multithreshold()

Examples

if(interactive()){

 #load example distance matrix
 data(distance_matrix)

 #PCA of the distance matrix
 out <- pca(x = distance_matrix)
 out

}

Description

Computes PCA factors of a distance matrix over different distance thresholds to generate spatial predictors for a model fitted with rf_spatial().

Usage

pca_multithreshold(
  distance.matrix = NULL,
  distance.thresholds = NULL,
  max.spatial.predictors = NULL
)

Arguments

Details

The distance matrix is converted into weights with weights_from_distance_matrix() before computing the PCA. This produces more meaningful spatial predictors than using the distance matrix as is.

Value

A data frame with the PCA factors of the thresholded matrix. The data frame columns are named "spatial_predictor_DISTANCE_COLUMN", where DISTANCE is the given distance threshold, and COLUMN is the column index of the given predictor.

See Also

pca()

Examples

if(interactive()){

 #loading example distance matrix
 load(distance_matrix)

 #PCA factors of the distance matrix for two reference distances
 x <- pca_multithreshold(
   distance.matrix = distance_matrix,
   distance.thresholds = c(0, 1000)
   )
 head(x)

}

Description

Richness of vascular plants of the American ecoregions as defined in Ecoregions 2017.

Usage

data(plant_richness_df)

Format

A data frame with 227 rows and 22 columns:

• ecoregion_id: Id of the ecoregion).

• x: Longitude in degrees (WGS84).

• y: Latitude in degrees (WGS84).

• richness_species_vascular: Number of vascular species found in the ecoregion. Response variable.

• bias_area_km2: Area of the ecoregion in squared kilometers.

• bias_species_per_record: Number of species divided by the number of spatial GBIF records available in the ecoregion as a measure of sampling bias.

• climate_aridity_index_average: Average of the ecoregion.

• climate_hypervolume: Volume of the climatic envelope of the ecoregion, computed with the hypervolume package.

• climate_velocity_lgm_average: Average climate velocity of the ecoregion since the Last Glacial Maximum.

• neighbors_count: Number of immediate neighbors of the ecoregion as a measure of connectivity/isolation.

• neighbors_percent_shared_edge: Percentage of shared edge with the neighbors as a measure of connectivity/isolation.

• human_population_density: Population density of the ecoregion.

• topography_elevation_average: Average elevation of the ecoregion.

• landcover_herbs_percent_average: Average cover percentage of herbs extracted from MODIS Vegetation Continuous Fields.

• fragmentation_cohesion: Geographic fragmentation index of the ecoregion as computed with the R package landscapemetrics.

• fragmentation_division: Another fragmentation index.

• neighbors_area: Total area of the ecoregions's immediate neighbors.

• human_population: Human population in the ecoregion.

• human_footprint_average: Average human footprint in the ecoregion.

• climate_bio1_average: Average mean annual temperature.

• climate_bio15_minimum: Average precipitation seasonality.

See Also

distance_matrix

Description

Plots the results of an spatial cross-validation performed with rf_evaluate().

Usage

plot_evaluation(
  model,
  fill.color = viridis::viridis(
    3,
    option = "F",
    alpha = 0.8,
    direction = -1
    ),
  line.color = "gray30",
  verbose = TRUE,
  notch = TRUE
)

Arguments

Value

A ggplot.

See Also

rf_evaluate(), get_evaluation(), print_evaluation().

Examples

if(interactive()){

#loading example data
data(plant_richness_df)
data(distance_matrix)

#fitting a random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1,
  verbose = FALSE
)

#evaluating the model with spatial cross-validation
rf.model <- rf_evaluate(
  model = rf.model,
  xy = plant_richness_df[, c("x", "y")],
  n.cores = 1
)

#plotting the evaluation results
plot_evaluation(rf.model)

}

Description

Plots variable importance scores of rf(), rf_repeat(), and rf_spatial() models. Distributions of importance scores produced with rf_repeat() are plotted using ggplot2::geom_violin, which shows the median of the density estimate rather than the actual median of the data. However, the violin plots are ordered from top to bottom by the real median of the data to make small differences in median importance easier to spot. Ths function does not plot the result of rf_importance() yet, but you can find it under model$importance$cv.per.variable.plot.

Usage

plot_importance(
  model,
  fill.color = viridis::viridis(
    100,
    option = "F",
    direction = -1,
    alpha = 1,
    end = 0.9
   ),
  line.color = "white",
  verbose = TRUE
)

Arguments

Value

A ggplot.

See Also

print_importance(), get_importance()

Examples

if(interactive()){

#loading example data
data(plant_richness_df)
data(distance_matrix)

#fitting a random forest model
rf.model <- rf(
 data = plant_richness_df,
 dependent.variable.name = "richness_species_vascular",
 predictor.variable.names = colnames(plant_richness_df)[5:21],
 distance.matrix = distance_matrix,
 distance.thresholds = 0,
 n.cores = 1,
 verbose = FALSE
)

#plotting variable importance scores
plot_importance(model = rf.model)

}

Description

Plots the results of spatial autocorrelation tests for a variety of functions within the package. The x axis represents the Moran's I estimate, the y axis contains the values of the distance thresholds, the dot sizes represent the p-values of the Moran's I estimate, and the red dashed line represents the theoretical null value of the Moran's I estimate.

Usage

plot_moran(
  model,
  point.color = viridis::viridis(
    100,
    option = "F",
    direction = -1
   ),
  line.color = "gray30",
  option = 1,
  ncol = 1,
  verbose = TRUE
)

Arguments

Value

A ggplot.

See Also

moran(), moran_multithreshold()

Examples

if(interactive()){

 #loading example data
 data(plant_richness_df)
 data(distance.matrix)

 #fitting a random forest model
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = c(0, 1000, 2000),
   n.cores = 1,
   verbose = FALSE
 )

 #Incremental/multiscale Moran's I
 plot_moran(rf.model)

 #Moran's scatterplot
 plot_moran(rf.model, option = 2)

}

Description

Plots optimization data frames produced by select_spatial_predictors_sequential() and select_spatial_predictors_recursive().

Usage

plot_optimization(
  model,
  point.color = viridis::viridis(
    100,
    option = "F",
    direction = -1
  ),
  verbose = TRUE
)

Arguments

Details

If the method used to fit a model with rf_spatial() is "hengl", the function returns nothing, as this method does not require optimization.

Value

A ggplot.

Examples

if(interactive()){

 #loading example data
 data(distance_matrix)
 data(plant_richness_df)

 #names of the response and predictors
 dependent.variable.name <- "richness_species_vascular"
 predictor.variable.names <- colnames(plant_richness_df)[5:21]

 #spatial model
 model <- rf_spatial(
   data = plant_richness_df,
   dependent.variable.name = dependent.variable.name,
   predictor.variable.names = predictor.variable.names,
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   method = "mem.moran.sequential",
   n.cores = 1,
   seed = 1
 )

 #plotting selection of spatial predictors
 plot_optimization(model = model)


}

Description

Plots normality and autocorrelation tests of model residuals.

Usage

plot_residuals_diagnostics(
  model,
  point.color = viridis::viridis(100, option = "F"),
  line.color = "gray10",
  fill.color = viridis::viridis(4, option = "F", alpha = 0.95)[2],
  option = 1,
  ncol = 1,
  verbose = TRUE
)

Arguments

Value

A patchwork object.

Examples

if(interactive()){

 #load example data
 data(plant_richness_df)
 data(distance_matrix)

 #fit a random forest model
 x <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   n.cores = 1
 )

 #residuals diagnostics
 plot_residuals_diagnostics(x)

}

Description

Plots the response curves of models fitted with rf(), rf_repeat(), or rf_spatial().

Usage

plot_response_curves(
  model = NULL,
  variables = NULL,
  quantiles = c(0.1, 0.5, 0.9),
  grid.resolution = 200,
  line.color = viridis::viridis(length(quantiles), option = "F", end = 0.9),
  ncol = 2,
  show.data = FALSE,
  verbose = TRUE
)

Arguments

Details

All variables that are not plotted in a particular response curve are set to the values of their respective quantiles, and the response curve for each one of these quantiles is shown in the plot. When the input model was fitted with rf_repeat() with keep.models = TRUE, then the plot shows the median of all model runs, and each model run separately as a thinner line. The output list can be plotted all at once with patchwork::wrap_plots(p) or cowplot::plot_grid(plotlist = p), or one by one by extracting each plot from the list.

Value

A list with slots named after the selected variables, with one ggplot each.

See Also

plot_response_surface()

Examples

if(interactive()){

#loading example data
data(plant_richness_df)

#fitting a random forest model
m <- rf(
 data = plant_richness_df,
 dependent.variable.name = "richness_species_vascular",
 predictor.variable.names = colnames(plant_richness_df)[5:21],
 n.cores = 1,
 verbose = FALSE
)

#response curves of most important predictors
plot_response_curves(model = m)

}

Description

Plots response surfaces for any given pair of predictors in a rf(), rf_repeat(), or rf_spatial() model.

Usage

plot_response_surface(
  model = NULL,
  a = NULL,
  b = NULL,
  quantiles = 0.5,
  grid.resolution = 100,
  point.size.range = c(0.5, 2.5),
  point.alpha = 1,
  fill.color = viridis::viridis(100, option = "F", direction = -1, alpha = 0.9),
  point.color = "gray30",
  verbose = TRUE
)

Arguments

Details

All variables that are not a or b in a response curve are set to the values of their respective quantiles to plot the response surfaces. The output list can be plotted all at once with patchwork::wrap_plots(p) or cowplot::plot_grid(plotlist = p), or one by one by extracting each plot from the list.

Value

A list with slots named after the selected quantiles, each one with a ggplot.

See Also

plot_response_curves()

Examples

if(interactive()){

#load example data
data(plant_richness_df)

#fit random forest model
out <- rf(
 data = plant_richness_df,
 dependent.variable.name = "richness_species_vascular",
 predictor.variable.names = colnames(plant_richness_df)[5:21],
 n.cores = 1,
 verbose = FALSE
)

#plot interactions between most important predictors
plot_response_surfaces(x = out)


}

Description

Plots the dependent variable against each predictor.

Usage

plot_training_df(
  data = NULL,
  dependent.variable.name = NULL,
  predictor.variable.names = NULL,
  ncol = 4,
  method = "loess",
  point.color = viridis::viridis(100, option = "F"),
  line.color = "gray30"
)

Arguments

Value

A wrap_plots object.

Examples

if(interactive()){

   #load example data
   data(plant_richness_df)

   #scatterplot of the training data
   plot_training_data(
     data = plant_richness_df,
     dependent.variable.name = "richness_species_vascular",
     predictor.variable.names = colnames(plant_richness_df)[5:21]
     )
 }

Description

Plots the the Moran's I test of the response and the predictors in a training data frame.

Usage

plot_training_df_moran(
  data = NULL,
  dependent.variable.name = NULL,
  predictor.variable.names = NULL,
  distance.matrix = NULL,
  distance.thresholds = NULL,
  fill.color = viridis::viridis(100, option = "F", direction = -1),
  point.color = "gray30"
)

Arguments

Value

A ggplot2 object.

Examples

if(interactive()){

   #load example data
   data(plant_richness_df)
   data(distance_matrix)

   #plot Moran's I of training data
   plot_moran_training_data(
     data = plant_richness_df,
     dependent.variable.name = "richness_species_vascular",
     predictor.variable.names = colnames(plant_richness_df)[5:21],
     distance.matrix = distance_matrix,
     distance.thresholds = c(
       0,
       2000,
       4000,
       6000,
       8000
       )
     )
}

Plots a tuning object produced by rf_tuning()

Description

Plots the tuning of the hyperparameters num.trees, mtry, and min.node.size performed by rf_tuning().

Usage

plot_tuning(
  model,
  point.color = viridis::viridis(
    100,
    option = "F"
  ),
  verbose = TRUE
)

Arguments

Value

A ggplot.

See Also

rf_tuning()

Examples

if(interactive()){

#load example data
data(plant_richness_df)

#fit random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1,
  verbose = FALSE
)

#tune random forest model
rf.model <- rf_tuning(
 model = rf.model,
 xy = plant_richness_df[, c("x", "y")],
 n.cores = 1,
 verbose = FALSE
)

#generate tuning plot
plot_tuning(model = rf.model)

}

Description

Prepares variable importance data frames and plots for models fitted with rf_spatial().

Usage

prepare_importance_spatial(model)

Arguments

Value

A list with importance data frames in different formats depending on whether the model was fitted with rf() or rf_repeat().

Examples

if(interactive()){

 #loading example data
 data(distance_matrix)
 data(plant_richness_df)

 #fittind spatial model
 model <- rf_spatial(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds =  0,
   n.cores = 1
 )

 #preparing the importance data frame
 importance <- prepare_importance_spatial(model)
 names(importance)

}

Description

Prints the results of an spatial cross-validation performed with rf_evaluate().

Usage

print_evaluation(model)

Arguments

Value

A table printed to the standard output.

See Also

plot_evaluation(), get_evaluation()

Examples

if(interactive()){

#loading example data
data(plant_richness_df)
data(distance_matrix)

#fitting random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1,
  verbose = FALSE
)

#evaluation with spatial cross-validation
rf.model <- rf_evaluate(
  model = rf.model,
  xy = plant_richness_df[, c("x", "y")],
  n.cores = 1
)

#checking evaluation results
print_evaluation(rf.model)

}

Description

Prints variable importance scores from rf, rf_repeat, and rf_spatial models.

Usage

print_importance(
  model,
  verbose = TRUE
)

Arguments

Value

A table printed to the standard output.

See Also

plot_importance(), get_importance()

Examples

if(interactive()){

#loading example data
data(plant_richness_df)
data(distance.matrix)

#fitting a random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1,
  verbose = FALSE
)

#printing variable importance scores
print_importance(model = rf.model)

}

Description

Prints the results of a Moran's I test on the residuals of a model.

Usage

print_moran(
  model,
  caption = NULL,
  verbose = TRUE
)

Arguments

Value

Prints a table in the console using the huxtable package.

See Also

moran(), moran_multithreshold(), get_moran(), plot_moran()

Examples

if(interactive()){

 #loading example data
 data(plant_richness_df)
 data(distance.matrix)

 #fitting random forest model
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = c(0, 1000, 2000),
   n.cores = 1,
   verbose = FALSE
 )

 #printing Moran's I of model's residuals
 print_moran(rf.model)

}

Description

Prints the performance slot of a model fitted with rf(), rf_repeat(), or rf_spatial(). For models fitted with rf_repeat() it shows the median and the median absolute deviation of each performance measure.

Usage

print_performance(model)

Arguments

Value

Prints model performance scores to the console.

See Also

print_performance(), get_performance()

Examples

if(interactive()){

 #loading example data
 data(plant_richness_df)
 data(distance.matrix)

 #fitting a random forest model
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   n.cores = 1,
   verbose = FALSE
 )

 #printing performance scores
 print_performance(rf.model)

}

Description

Custom print method for models fitted with rf(), rf_repeat(), and rf_spatial().

Usage

## S3 method for class 'rf'
print(x, ...)

Arguments

Value

Prints model details to the console.

See Also

print_evaluation(), print_importance(), print_moran(), print_performance()

Examples

if(interactive()){

 #loading example data
 data("plant_richness_df")
 data("distance_matrix")

 #fitting random forest model
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   n.cores = 1
 )

 #printing model summary
 print(rf.model)
}

Description

Ranks spatial predictors generated by mem_multithreshold() or pca_multithreshold() by their effect in reducing the Moran's I of the model residuals (ranking.method = "effect"), or by their own Moran's I (ranking.method = "moran").

In the former case, one model of the type y ~ predictors + spatial_predictor_X is fitted per spatial predictor, and the Moran's I of this model's residuals is compared with the one of the model without spatial predictors (y ~ predictors), to finally rank the spatial predictor from maximum to minimum difference in Moran's I.

In the latter case, the spatial predictors are ordered by their Moran's I alone (this is the faster option).

In both cases, spatial predictors that are redundant with others at a Pearson correlation > 0.5 and spatial predictors with no effect (no reduction of Moran's I or Moran's I of the spatial predictor equal or lower than 0) are removed.

This function has been designed to be used internally by rf_spatial() rather than directly by a user.

Usage

rank_spatial_predictors(
  data = NULL,
  dependent.variable.name = NULL,
  predictor.variable.names = NULL,
  distance.matrix = NULL,
  distance.thresholds = NULL,
  ranger.arguments = NULL,
  spatial.predictors.df = NULL,
  ranking.method = c("moran", "effect"),
  reference.moran.i = 1,
  verbose = FALSE,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Value

A list with four slots:

• method: Character, name of the method used to rank the spatial predictors.

• criteria: Data frame with two different configurations depending on the ranking method. If ranking.method = "effect", the columns contain the names of the spatial predictors, the r-squared of the model, the Moran's I of the model residuals, the difference between the Moran's I of the model including the given spatial predictor, and the Moran's I of the model fitted without spatial predictors, and the interpretation of the Moran's I value. If ranking.method = "moran", only the name of the spatial predictor and it's Moran's I are in the output data frame.

• ranking: Ordered character vector with the names of the spatial predictors selected.

• spatial.predictors.df: data frame with the selected spatial predictors in the order of the ranking.

Examples

if(interactive()){

 #loading distance matrix
 data(distance_matrix)

 #computing Moran's Eigenvector Maps
 mem.df <- mem(
  distance.matrix = distance_matrix[1:50, 1:50],
  distance.threshold = 0
 )

 #ranking by the Moran's I of the spatial predictor
 rank <- rank_spatial_predictors(
  distance.matrix = distance_matrix[1:50, 1:50],
  distance.thresholds = 0,
  spatial.predictors.df = mem.df,
  ranking.method = "moran",
  n.cores = 1
 )

 #checking Moran's I of MEMs
 rank$criteria

 #checking rank of MEMs
 rank$ranking
}

Description

Rescales a numeric vector to a new range.

Usage

rescale_vector(
  x = NULL,
  new.min = 0,
  new.max = 1,
  integer = FALSE
)

Arguments

Value

A numeric vector of the same length as x, but with its values rescaled between new.min and new.max.

Examples

if(interactive()){

 out <- rescale_vector(
   x = rnorm(100),
   new.min = 0,
   new.max = 100,
   integer = TRUE
   )
   out

}

Description

Applies a Shapiro-Wilks test to a numeric vector, and plots the qq plot and the histogram.

Usage

residuals_diagnostics(residuals, predictions)

Arguments

Details

The function shapiro.test() has a hard limit of 5000 cases. If the model residuals have more than 5000 cases, then sample(x = residuals, size = 5000) is applied to the model residuals before the test.

Value

A list with four slots:

/item w W statistic returned by shapiro.test(). /item p.value p-value of the Shapiro test. /item interpretation Character vector, one of "x is normal", "x is not normal". /item plot A patchwork plot with the qq plot and the histogram of x.

See Also

ggplot,aes,geom_qq_line,ggtheme,labs,geom_freqpoly,geom_abline plot_annotation

Examples

if(interactive()){

 residuals_diagnostics(
   residuals = runif(100),
   predictions = runif(100)
 )

 }

Description

Applies a Shapiro-Wilks test to a numeric vector, and returns a list with the statistic W, its p-value, and a character string with the interpretation.

Usage

residuals_test(residuals)

Arguments

Value

A list with four slots:

/item w W statistic returned by shapiro.test(). /item p.value p-value of the Shapiro test. /item interpretation Character vector, one of "x is normal", "x is not normal". /item plot A patchwork plot with the qq plot and the histogram of x.

See Also

ggplot,aes,geom_qq_line,ggtheme,labs,geom_freqpoly,geom_abline plot_annotation

Examples

if(interactive()){

 residuals_test(residuals = runif(100))

}

Description

A convenient wrapper for ranger that completes its output by providing the Moran's I of the residuals for different distance thresholds, the rmse and nrmse (as computed by root_mean_squared_error()), and variable importance scores based on a scaled version of the data generated by scale.

Usage

rf(
  data = NULL,
  dependent.variable.name = NULL,
  predictor.variable.names = NULL,
  distance.matrix = NULL,
  distance.thresholds = NULL,
  xy = NULL,
  ranger.arguments = NULL,
  scaled.importance = FALSE,
  seed = 1,
  verbose = TRUE,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Details

Please read the help file of ranger for further details. Notice that the formula interface of ranger is supported through ranger.arguments, but variable interactions are not allowed (but check the_feature_engineer()).

Value

A ranger model with several extra slots:

• ranger.arguments: Stores the values of the arguments used to fit the ranger model.

• importance: A list containing a data frame with the predictors ordered by their importance, a ggplot showing the importance values, and local importance scores (difference in accuracy between permuted and non permuted variables for every case, computed on the out-of-bag data).

• performance: performance scores: R squared on out-of-bag data, R squared (cor(observed, predicted) ^ 2), pseudo R squared (cor(observed, predicted)), RMSE, and normalized RMSE (NRMSE).

• residuals: residuals, normality test of the residuals computed with residuals_test(), and spatial autocorrelation of the residuals computed with moran_multithreshold().

Examples

if(interactive()){

 #loading example data
 data("plant_richness_df")
 data("distance_matrix")

 #fittind random forest model
 out <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   n.cores = 1
 )

 class(out)

 #data frame with ordered variable importance
 out$importance$per.variable

 #variable importance plot
 out$importance$per.variable.plot

 #performance
 out$performance

 #spatial correlation of the residuals
 out$spatial.correlation.residuals$per.distance

 #plot of the Moran's I of the residuals for different distance thresholds
 out$spatial.correlation.residuals$plot

 #predictions for new data as done with ranger models:
 predicted <- stats::predict(
   object = out,
   data = plant_richness_df,
   type = "response"
 )$predictions

 #alternative data input methods
 ###############################

 #ranger.arguments can contain ranger arguments and any other rf argument
 my.ranger.arguments <- list(
 data = plant_richness_df,
 dependent.variable.name = "richness_species_vascular",
 predictor.variable.names = colnames(plant_richness_df)[8:21],
 distance.matrix = distance_matrix,
 distance.thresholds = c(0, 1000)
 )

 #fitting model with these ranger arguments
 out <- rf(
   ranger.arguments = my.ranger.arguments,
   n.cores = 1
   )

}

Description

Uses rf_evaluate() to compare the performance of several models on independent spatial folds via spatial cross-validation.

Usage

rf_compare(
  models = NULL,
  xy = NULL,
  repetitions = 30,
  training.fraction = 0.75,
  metrics = c("r.squared", "pseudo.r.squared", "rmse", "nrmse", "auc"),
  distance.step = NULL,
  distance.step.x = NULL,
  distance.step.y = NULL,
  fill.color = viridis::viridis(100, option = "F", direction = -1, alpha = 0.8),
  line.color = "gray30",
  seed = 1,
  verbose = TRUE,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Value

A list with three slots:

• comparison.df: Data frame with one performance value per spatial fold, metric, and model.

• spatial.folds: List with the indices of the training and testing records for each evaluation repetition.

• plot: Violin-plot of comparison.df.

See Also

rf_evaluate()

Examples

if(interactive()){

 #loading example data
 data(distance_matrix)
 data(plant_richness_df)

 #fitting random forest model
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   n.cores = 1
 )

 #fitting a spatial model with Moran's Eigenvector Maps
 rf.spatial <- rf_spatial(
 model = rf.model,
 n.cores = 1
 )

 #comparing the spatial and non spatial models
 comparison <- rf_compare(
 models = list(
   `Non spatial` = rf.model,
   Spatial = rf.spatial
 ),
 xy = plant_richness_df[, c("x", "y")],
 metrics = c("r.squared", "rmse"),
 n.cores = 1
 )

}

Description

Evaluates the performance of random forest on unseen data over independent spatial folds.

Usage

rf_evaluate(
  model = NULL,
  xy = NULL,
  repetitions = 30,
  training.fraction = 0.75,
  metrics = c("r.squared", "pseudo.r.squared", "rmse", "nrmse", "auc"),
  distance.step = NULL,
  distance.step.x = NULL,
  distance.step.y = NULL,
  grow.testing.folds = FALSE,
  seed = 1,
  verbose = TRUE,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Details

The evaluation algorithm works as follows: the number of repetitions and the input dataset (stored in model$ranger.arguments$data) are used as inputs for the function thinning_til_n(), that applies thinning() to the input data until as many cases as repetitions are left, and as separated as possible. Each of these remaining records will be used as a "fold center". From that point, the fold grows, until a number of points equal (or close) to training.fraction is reached. The indices of the records within the grown spatial fold are stored as "training" in the output list, and the remaining ones as "testing". Then, for each spatial fold, a "training model" is fitted using the cases corresponding with the training indices, and predicted over the cases corresponding with the testing indices. The model predictions on the "unseen" data are compared with the observations, and the performance measures (R squared, pseudo R squared, RMSE and NRMSE) computed.

Value

A model of the class "rf_evaluate" with a new slot named "evaluation", that is a list with the following slots:

• training.fraction: Value of the argument training.fraction.

• spatial.folds: Result of applying make_spatial_folds() on the data coordinates. It is a list with as many slots as repetitions are indicated by the user. Each slot has two slots named "training" and "testing", each one having the indices of the cases used on the training and testing models.

• per.fold: Data frame with the evaluation results per spatial fold (or repetition). It contains the ID of each fold, it's central coordinates, the number of training and testing cases, and the training and testing performance measures: R squared, pseudo R squared (cor(observed, predicted)), rmse, and normalized rmse.

• per.model: Same data as above, but organized per fold and model ("Training", "Testing", and "Full").

• aggregated: Same data, but aggregated by model and performance measure.

Examples

if(interactive()){

#loading example data
data(plant_richness_df)
data(distance_matrix)

#fitting random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1,
  verbose = FALSE
)

#evaluation with spatial cross-validation
rf.model <- rf_evaluate(
  model = rf.model,
  xy = plant_richness_df[, c("x", "y")],
  n.cores = 1
)

#checking evaluation results
plot_evaluation(rf.model)
print_evaluation(rf.model)
x <- get_evaluation(rf.model)

}

Description

Evaluates the contribution of the predictors to model transferability via spatial cross-validation. The function returns the median increase or decrease in a given evaluation metric (R2, pseudo R2, RMSE, nRMSE, or AUC) when a variable is introduced in a model, by comparing and evaluating via spatial cross-validation models with and without the given variable. This function was devised to provide importance scores that would be less sensitive to spatial autocorrelation than those computed internally by random forest on the out-of-bag data. This function is experimental.

Usage

rf_importance(
  model = NULL,
  xy = NULL,
  repetitions = 30,
  training.fraction = 0.75,
  metric = c("r.squared", "pseudo.r.squared", "rmse", "nrmse", "auc"),
  distance.step = NULL,
  distance.step.x = NULL,
  distance.step.y = NULL,
  fill.color = viridis::viridis(100, option = "F", direction = -1, alpha = 1, end = 0.9),
  line.color = "white",
  seed = 1,
  verbose = TRUE,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Value

The input model with new data in its "importance" slot. The new importance scores are included in the data frame model$importance$per.variable, under the column names "importance.cv" (median contribution to transferability over spatial cross-validation repetitions), "importance.cv.mad" (median absolute deviation of the performance scores over spatial cross-validation repetitions), "importance.cv.percent" ("importance.cv" expressed as a percent, taking the full model's performance as baseline), and "importance.cv.mad" (median absolute deviation of "importance.cv"). The plot is stored as "cv.per.variable.plot".

Examples

if(interactive()){

#loading example data
data(plant_richness_df)
data(distance_matrix)
xy <- plant_richness_df[, c("x", "y")]

#fitting random forest model
rf.model <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  xy = xy,
  n.cores = 1,
  verbose = FALSE
)

#computing predictor contribution to model transferability
rf.model <- rf_importance(rf.model)

}

Description

Fits several random forest models on the same data in order to capture the effect of the algorithm's stochasticity on the variable importance scores, predictions, residuals, and performance measures. The function relies on the median to aggregate performance and importance values across repetitions. It is recommended to use it after a model is fitted (rf() or rf_spatial()), tuned (rf_tuning()), and/or evaluated (rf_evaluate()). This function is designed to be used after fitting a model with rf() or rf_spatial(), tuning it with rf_tuning() and evaluating it with rf_evaluate().

Usage

rf_repeat(
  model = NULL,
  data = NULL,
  dependent.variable.name = NULL,
  predictor.variable.names = NULL,
  distance.matrix = NULL,
  distance.thresholds = NULL,
  xy = NULL,
  ranger.arguments = NULL,
  scaled.importance = FALSE,
  repetitions = 10,
  keep.models = TRUE,
  seed = 1,
  verbose = TRUE,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Value

A ranger model with several new slots:

• ranger.arguments: Stores the values of the arguments used to fit the ranger model.

• importance: A list containing a data frame with the predictors ordered by their importance, a ggplot showing the importance values, and local importance scores.

• performance: out-of-bag performance scores: R squared, pseudo R squared, RMSE, and normalized RMSE (NRMSE).

• pseudo.r.squared: computed as the correlation between the observations and the predictions.

• residuals: residuals, normality test of the residuals computed with residuals_test(), and spatial autocorrelation of the residuals computed with moran_multithreshold().

Examples

if(interactive()){

 #loading example data
 data(plant_richness_df)
 data(distance_matrix)

 #fitting 5 random forest models
 out <- rf_repeat(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   repetitions = 5,
   n.cores = 1
 )

 #data frame with ordered variable importance
 out$importance$per.variable

 #per repetition
 out$importance$per.repetition

 #variable importance plot
 out$importance$per.repetition.plot

 #performance
 out$performance


 #spatial correlation of the residuals for different distance thresholds
 out$spatial.correlation.residuals$per.distance

 #plot of the Moran's I of the residuals for different distance thresholds
 out$spatial.correlation.residuals$plot

 #using a model as an input for rf_repeat()
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[8:21],
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   n.cores = 1
   )

 #repeating the model 5 times
 rf.repeat <- rf_repeat(
   model = rf.model,
   n.cores = 1
   )


 rf.repeat$performance
 rf.repeat$importance$per.repetition.plot

}

Description

Fits spatial random forest models using different methods to generate, rank, and select spatial predictors acting as proxies of spatial processes not considered by the non-spatial predictors. The end goal is providing the model with information about the spatial structure of the data to minimize the spatial correlation (Moran's I) of the model residuals and generate honest variable importance scores.

Usage

rf_spatial(
  model = NULL,
  data = NULL,
  dependent.variable.name = NULL,
  predictor.variable.names = NULL,
  distance.matrix = NULL,
  distance.thresholds = NULL,
  xy = NULL,
  ranger.arguments = NULL,
  scaled.importance = TRUE,
  method = c("mem.moran.sequential", "mem.effect.sequential", "mem.effect.recursive",
    "hengl", "hengl.moran.sequential", "hengl.effect.sequential",
    "hengl.effect.recursive", "pca.moran.sequential", "pca.effect.sequential",
    "pca.effect.recursive"),
  max.spatial.predictors = NULL,
  weight.r.squared = NULL,
  weight.penalization.n.predictors = NULL,
  seed = 1,
  verbose = TRUE,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Details

The function uses three different methods to generate spatial predictors ("hengl", "pca", and "mem"), two methods to rank them in order to define in what order they are introduced in the model ("effect" and "moran), and two methods to select the spatial predictors that minimize the spatial correlation of the model residuals ("sequential" and "recursive"). All method names but "hengl" (that uses the complete distance matrix as predictors in the spatial model) are named by combining a method to generate the spatial predictors, a method to rank them, and a method to select them, separated by a point. Examples are "mem.moran.sequential" or "mem.effect.recursive". All combinations are not possible, since the ranking method "moran" cannot be used with the selection method "recursive" (because the logics behind them are very different, see below). Methods to generate spatial predictors:

• "hengl": named after the method RFsp presented in the paper "Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables", by Hengl et al. (2018), where the authors propose to use the distance matrix among records as predictors in spatial random forest models (RFsp method). In this function, all methods starting with "hengl" use either the complete distance matrix, or select columns of the distance matrix as spatial predictors.

• "mem": Generates Moran's Eigenvector Maps, that is, the eigenvectors of the double-centered weights of the distance matrix. The method is described in "Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM)", by Dray et al. (2006), and "Statistical methods for temporal and space–time analysis of community composition data", by Legendre and Gauthier (2014).

• "pca": Computes spatial predictors from the principal component analysis of a weighted distance matrix (see weights_from_distance_matrix()). This is an experimental method, use with caution.

Methods to rank spatial predictors (see rank_spatial_predictors()):

• "moran": Computes the Moran's I of each spatial predictor, selects the ones with positive values, and ranks them from higher to lower Moran's I.

• "effect": If a given non-spatial random forest model is defined as y = p1 + ... + pn, being p1 + ... + pn the set of predictors, for every spatial predictor generated (spX) a spatial model y = p1 + ... + pn + spX is fitted, and the Moran's I of its residuals is computed. The spatial predictors are then ranked by how much they help to reduce spatial autocorrelation between the non-spatial and the spatial model.

Methods to select spatial predictors:

• "sequential" (see select_spatial_predictors_sequential()): The spatial predictors are added one by one in the order they were ranked, and once all spatial predictors are introduced, the best first n predictors are selected. This method is similar to the one employed in the MEM methodology (Moran's Eigenvector Maps) described in the paper "Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM)", by Dray et al. (2006), and "Statistical methods for temporal and space–time analysis of community composition data", by Legendre and Gauthier (2014). This method generally introduces tens of predictors into the model, but usually offers good results.

• "recursive" (see select_spatial_predictors_recursive()): This method tries to find the smallest combination of spatial predictors that reduce the spatial correlation of the model's residuals the most. The algorithm goes as follows: 1. The first ranked spatial predictor is introduced into the model; 2. the remaining predictors are ranked again using the "effect" method, using the model in 1. as reference. The first spatial predictor in the resulting ranking is then introduced into the model, and the steps 1. and 2. are repeated until spatial predictors stop having an effect in reducing the Moran's I of the model residuals. This method takes longer to compute, but generates smaller sets of spatial predictors. This is an experimental method, use with caution.

Once ranking procedure is completed, an algorithm is used to select the minimal subset of spatial predictors that reduce the most the Moran's I of the residuals: for each new spatial predictor introduced in the model, the Moran's I of the residuals, it's p-value, a binary version of the p-value (0 if < 0.05 and 1 if >= 0.05), the R-squared of the model, and a penalization linear with the number of spatial predictors introduced (computed as ⁠(1 / total spatial predictors) * introduced spatial predictors⁠) are rescaled between 0 and 1. Then, the optimization criteria is computed as ⁠max(1 - Moran's I, p-value binary) + (weight.r.squared * R-squared) - (weight.penalization.n.predictors * penalization)⁠. The predictors from the first one to the one with the highest optimization criteria are then selected as the best ones in reducing the spatial correlation of the model residuals, and used along with data to fit the final spatial model.

Value

A ranger model with several new slots:

• ranger.arguments: Values of the arguments used to fit the ranger model.

• importance: A list containing the vector of variable importance as originally returned by ranger (scaled or not depending on the value of 'scaled.importance'), a data frame with the predictors ordered by their importance, and a ggplot showing the importance values.

• performance: With the out-of-bag R squared, pseudo R squared, RMSE and NRMSE of the model.

• residuals: residuals, normality test of the residuals computed with residuals_test(), and spatial autocorrelation of the residuals computed with moran_multithreshold().

• spatial: A list with four slots:

  • method: Character, method used to generate, rank, and select spatial predictors.

  • names: Character vector with the names of the selected spatial predictors. Not returned if the method is "hengl".

  • optimization: Criteria used to select the spatial predictors. Not returned if the method is "hengl".

  • plot: Plot of the criteria used to select the spatial predictors. Not returned if the method is "hengl".

Examples

if(interactive()){

 #loading example data
 data(distance_matrix)
 data(plant_richness_df)

 #names of the response and predictors
 dependent.variable.name <- "richness_species_vascular"
 predictor.variable.names <- colnames(plant_richness_df)[5:21]

 #hengl
 model <- rf_spatial(
   data = plant_richness_df,
   dependent.variable.name = dependent.variable.name,
   predictor.variable.names = predictor.variable.names,
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   method = "hengl",
   n.cores = 1
 )

 #mem.moran.sequential
 model <- rf_spatial(
   data = plant_richness_df,
   dependent.variable.name = dependent.variable.name,
   predictor.variable.names = predictor.variable.names,
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   method = "mem.moran.sequential",
   n.cores = 1
 )

 #fitting an rf_spatial model from an rf model
 rf.model <- rf(
   data = plant_richness_df,
   dependent.variable.name = "richness_species_vascular",
   predictor.variable.names = colnames(plant_richness_df)[5:21],
   distance.matrix = distance_matrix,
   distance.thresholds = 0,
   n.cores = 1,
   verbose = FALSE
 )
 rf.model$spatial.correlation.residuals$plot

 #spatial version of the rf model
 rf.spatial <- rf_spatial(model = rf.model)
 rf.spatial$spatial.correlation.residuals$plot


}

Description

Finds the optimal set of random forest hyperparameters num.trees, mtry, and min.node.size via grid search by maximizing the model's R squared, or AUC, if the response variable is binomial, via spatial cross-validation performed with rf_evaluate().

Usage

rf_tuning(
  model = NULL,
  num.trees = NULL,
  mtry = NULL,
  min.node.size = NULL,
  xy = NULL,
  repetitions = 30,
  training.fraction = 0.75,
  seed = 1,
  verbose = TRUE,
  n.cores = parallel::detectCores() - 1,
  cluster = NULL
)

Arguments

Value

A model with a new slot named tuning, with a data frame with the results of the tuning analysis.

See Also

rf_evaluate()

Examples

if(interactive()){

#loading example data
data(plant_richness_df)
data(distance_matrix)

#fitting model to tune
out <- rf(
  data = plant_richness_df,
  dependent.variable.name = "richness_species_vascular",
  predictor.variable.names = colnames(plant_richness_df)[5:21],
  distance.matrix = distance_matrix,
  distance.thresholds = 0,
  n.cores = 1
)

#model tuning
tuning <- rf_tuning(
  model = out,
  num.trees = c(100, 500),
  mtry = c(2, 8),
  min.node.size = c(5, 10),
  xy = plant_richness_df[, c("x", "y")],
  n.cores = 1
)

}

Description

Computes the rmse or normalized rmse (nrmse) between two numeric vectors of the same length representing observations and model predictions.

Usage

root_mean_squared_error(
  o,
  p,
  normalization = c("rmse", "all", "mean", "sd", "maxmin", "iq")
)

Arguments

Details

The normalization methods go as follows:

• "rmse": RMSE with no normalization.

• "mean": RMSE dividied by the mean of the observations (rmse/mean(o)).

• "sd": RMSE dividied by the standard deviation of the observations (rmse/sd(o)).

• "maxmin": RMSE divided by the range of the observations (rmse/(max(o) - min(o))).

• "⁠iq"⁠: RMSE divided by the interquartile range of the observations (rmse/(quantile(o, 0.75) - quantile(o, 0.25)))

Value

Named numeric vector with either one or 5 values, as selected by the user.

Examples

if(interactive()){

 root_mean_squared_error(
   o = runif(10),
   p = runif(10)
   )

}

Description

Selects spatial predictors following these steps:

1. Gets the spatial predictors ranked by rank_spatial_predictors() and fits a model of the form y ~ predictors + best_spatial_predictor_1. The Moran's I of the residuals of this model is used as reference value for the next step.

2. The remaining spatial predictors are introduced again into rank_spatial_predictors(), and the spatial predictor with the highest ranking is introduced in a new model of the form y ~ predictors + best_spatial_predictor_1 + best_spatial_predictor_2.

3. Steps 1 and 2 are repeated until the Moran's I doesn't improve for a number of repetitions equal to the 20 percent of the total number of spatial predictors introduced in the function.