Compare different regression methods

1. The compare_methods()

2. Basic example

For the basic example, we will use the dataset with the Mean Vessel Area (MVA) chronology and the mean April temperature, the dataset is saved as dataset_MVA. All five regression methods will be compared with 10-fold cross-validation repeated 1 times. The optimize argument is set to TRUE, therefore all tuning parameters will be defined in a preliminary optimization phase. After the comparison, the output elements are retrieved with the “$” operator.

# Load the dendroTools R package
library(dendroTools)

# Load the data
data(dataset_MVA)

# Basic example
basic_example <- compare_methods(formula = T_Apr ~ MVA, dataset = dataset_MVA, k = 10, repeats = 1, optimize = TRUE, MT_committees_vector = c(1), RF_maxnodes_vector = c(5), RF_nodesize_vector = c(10))

# The data frame with mean and standard deviation of performance metrics for the calibration and the validation data
kable(basic_example$mean_std)

# The data frame with non-parametric estimation of different methods: average rank and share of rank one
kable(basic_example$rank)

# See the histogram of mean bias for the validation data
basic_example$bias_val

Histogram for the validation data for the basic_example

# See the histogram of mean bias for the calibration data
basic_example$bias_cal

Histogram for the calibration data for the basic_example

# See the transfer functions, separated by facets. This is a ggplot object and could be easily customized. 
library(ggplot2)
basic_example$transfer_functions +   
  xlab(expression(paste('MVA [',mm^2,']'))) +
  ylab("April Mean Temperature [°C]")

The transfer functions of different methods, facet is used to separate plots by method.

# See the transfer functions, plotted together. This is a ggplot object and could be easily customized. 
basic_example$transfer_functions_together +   
  xlab(expression(paste('MVA [',mm^2,']'))) +
  ylab("April Mean Temperature [°C]")

The transfer functions of different methods, all functions are on the same plot, therefore it is easy to see the differences among different methods.

# The data frame of optimized tuning parameters for different methods
kable(basic_example$parameter_values)

# For calibration data, there are residual diagnostic plots available. Similar plots are available also for holdout and edge data. 

basic_example$normal_QQ_cal

Residual diagnostic plots for calibration data: Normal Q-Q plot

# For calibration data, there are residual diagnostic plots available. Similar plots are available also for holdout and edge data. 

basic_example$residuals_vs_fitted_cal

## `geom_smooth()` using formula 'y ~ x'

Residual diagnostic plots for calibration data: residuals vs fitted plot

3. Principal component analysis in combination with the compare_methods()

Principal Component Analysis (PCA) is commonly used with tree-ring data to reduce the full set of original tree-ring chronologies to a more manageable set of transformed variables. These transformed variables, the set of principal component scores, are then used as predictors in climate reconstruction models. The PCA also acts to strengthen the common regional-scale climate response within a group of tree-ring chronologies by concentrating the common signal in the components with the largest eigenvalues.

To use PCA regression within the compare_methods(), set the argument PCA_transformation to TRUE. All independent variables in the dataset data frame will be transformed using the PCA transformation. If the parameter log_preprocess is set to TRUE, variables will be transformed with logarithmic transformation before used in PCA. With the argument components_selection, we specify how to select PC scores that will be used as predictors. There are three options: “automatic”, “manual” and “plot_selection”. If argument is set to “automatic”, all PC scores with eigenvalues above 1 will be selected. This threshold could be changed by changing the eigenvalues_threshold argument. If argument is set to “manual”, user should set the number of components with N_components argument. If components_selection is set to “plot_selection”, A scree plot will be shown, and a user must manually enter the number of components to be used as predictors. The latter seems to be the most reasonable choice.

For the example with PCA, we use dataset dataset_MVA_individual, which consist of 10 individual Mean Vessel Area (MVA) chronologies and mean April temperature for the Ljubljana region, Slovenia. The dataset has 56 observations. The selection of components is set to “manual”, N_components to be used in the later analysis is set to 2. In this example, the 5-fold cross-validation with 1 repeat will be used to compare MT, MLR and BRNN. The subset of methods could be set with the methods argument. The argument optimize is set to TRUE, therefore all tuning parameters will be set automatically.

# Load the dendroTools R package
library(dendroTools)

# Load data
data(dataset_MVA_individual)

# Example PCA
example_PCA <- compare_methods(formula = T_Apr ~ ., dataset = dataset_MVA_individual, k = 5, repeats = 1, optimize = TRUE, methods = c("MLR", "MT", "BRNN"), PCA_transformation = TRUE, components_selection = "manual", N_components = 2, seed_factor = 5, MT_committees_vector = c(1), RF_maxnodes_vector = c(5), RF_nodesize_vector = c(10))

# Get the summary statistics for the PCA
summary(example_PCA$PCA_output)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     1.8691406 1.5007125 1.1867201 0.87168877 0.79432916
## Proportion of Variance 0.3493686 0.2252138 0.1408305 0.07598413 0.06309588
## Cumulative Proportion  0.3493686 0.5745824 0.7154129 0.79139704 0.85449292
##                            Comp.6     Comp.7     Comp.8     Comp.9     Comp.10
## Standard deviation     0.68818504 0.61467322 0.57233447 0.49611738 0.173060035
## Proportion of Variance 0.04735987 0.03778232 0.03275667 0.02461325 0.002994978
## Cumulative Proportion  0.90185279 0.93963510 0.97239178 0.99700502 1.000000000

# The mean and standard deviation data frame 
kable(example_PCA$mean_std)

4. Example of multiproxy analysis

The compare_methods() enables the comparison of methods for regression problems with two or more independent variables. However, users should select multiple proxies reasonably and with caution, since there is nothing to prevent from including colinear variables. To perform the comparison of methods with multiproxy variables, simply include dataset with more than one independent variable and specify the relationship with the formula argument. If metrics on validation data are much lower than on calibration data, there is a problem of overfitting and users should exclude some independent variables and repeat the analysis.

For the multiproxy_example, we will use example_dataset_1, which consist of the Mean Vessel Area (MVA) chronology and two temperature variables, the mean April temperature (T_APR) and the mean temperature from August to September (T_aug_sep) from the previous growing season. To compare methods with multiproxy approach, specify formula with two independent variables, such as formula = MVA ~ T_APR + T_aug_sep. Here, we will compare MT, BRNN and RF with 10-fold cross-validation and 1 repeat.

# Load the dendroTools R package
library(dendroTools)

# Load data
data(example_dataset_1)

# Example multiproxy
example_multiproxy <- compare_methods(formula = MVA ~ T_APR + T_aug_sep, dataset = example_dataset_1, k = 10, repeats = 1, optimize = FALSE, methods = c("MT", "BRNN", "RF"))

# The mean and standard deviation data frame 
kable(example_multiproxy$mean_std)

5. Example of climate reconstruction

Reconstructions of past climate conditions include reconstructions of past temperature, precipitation, vegetation, stream flow, sea surface temperature, and other climatic or climate-dependent conditions. With the compare_methods() it is possible to directly reconstruct the dependent variable specified with the formula argument. To do so, supply additional complete dataset with tree-ring chronology that goes beyond the observations of instrumental records.

For the example_reconstruction, we use data_TRW dataset, which includes a tree-ring width (TRW) chronology of Pinus nigra from Albania and mean June-July temperature from Albania. The complete TRW chronology is supplied with the dataset_TRW_complete. In this example, we will compare RF and MLR models with 3-fold cross-validation and 1 repeat.

# Load the dendroTools R package
library(dendroTools)

# Load the data
data(dataset_TRW)
data(dataset_TRW_complete)

# Example reconstruction
example_reconstruction <- compare_methods(formula = T_Jun_Jul ~ TRW, dataset = dataset_TRW, k = 3, optimize = FALSE, methods = c("MLR", "BRNN", "MT", "RF"), dataset_complete = dataset_TRW_complete)

example_reconstruction$reconstructions

The reconstructed June-July temperatures based on the dataset_complete argument, facet is used to split plots by methods.

example_reconstruction$reconstructions_together

The reconstructed June-July temperatures based on the dataset_complete argument, all reconstructions are on the same plot. The RF model reconstructed temperatures with much lower variance than the MLR model.

The comparison among different methods usually shows, that different regression techniques perform relatively similar for the central-calibration data, while for the extreme values, predictions differ. Therefore, we implemented additional central-edge test, where the central part of the dataset represents the calibration data, while the edge data, i.e. extreme values, represent the edge data. Different regression models are calibrated using the central data and validated using the edge (extreme) data. This test is particularly important to assess the performance of different methods for the prediction of the extreme and out of calibration data. The share of the edge (extreme) data is defined with the edge_share argument. To get the results for the central-edge test, use example_reconstruction$edge_results.

# The central-edge test
kable(example_reconstruction$edge_results)

In this example, MLR and MT had worse validation results than RF and BRNN. This indicates that linear interpolation for extreme values is not the best choice.

6. Tuning the machine learning parameters

Machine learning methods have several tuning parameters that need to be set, e.g. the number of neurons for the BRNN (parameter BRNN_neurons = 2). By default the optimize argument is set to TRUE, therefore all parameters will be automatically optimized in a preliminary cross-validation phase, where different combinations of parameters are tested and the best combination for each method is later used for the final model comparison. Each parameter has a pre-defined vector of possible values, however, this vector of possible values could be extended and therefore a wider space of tuned values could be explored. To change the vector of possible values for the BRNN_neurons parameter, use e.g. BRNN_neurons_vector = c(1, 2, 3, 4, 5). Bellow, see the table of all tuning parameters together with the vectors of possible values.

To set the tuning parameters manually, set the parameter optimize to FALSE and supply the selected value of each tuning parameter with the corresponding argument (if not, the default value will be used). Here is a simple example, where tuning parameters are set manually.

# Load the dendroTools R package
library(dendroTools)

# Load the data
data(example_dataset_1)

example_optimize <- compare_methods(formula = MVA ~  T_APR, dataset = example_dataset_1, k = 5, repeats = 2, optimize = FALSE, BRNN_neurons = 1, MT_committees = 1, MT_neighbors = 0, MT_rules = 100, MT_unbiased = FALSE, MT_extrapolation = 100, MT_sample = 0, RF_mtry = 1, RF_ntree = 100, RF_maxnodes = 20, seed_factor = 5)

# The data frame of tuning parameters, as defined by the user
kable(example_optimize$parameter_values)