Data Analysis - Predict Transitions

1 Introduction

Our model looks at a variety of ecological data, including organism abundance, human population, and temperature, to predict if there would be a regime shift. We got our information from a prior study found here: https://www.nature.com/articles/s41598-018-35057-4#data-availability

Import all libraries and data
- Some cleaning and manipulation for the analysis

library(dplyr)
library(reshape2)
library(ggplot2)
library(lubridate)
library(randomForest)
library(tidyr)
library(mlr3)
library(mlr3learners)
library(data.table)
library(mlr3tuning)
library(paradox)
library(mlr3viz)
library(mlr3mbo)
library(pROC)
library(mlr3pipelines)
library(mlr3torch)

analysisData <- read.csv("//Users/aidananderson/Documents/GitHub/ML-Workshop-Hawaii-Regime-Shifts/cleanedData1.csv")%>%
  select(-X)
scaledData <- read.csv("//Users/aidananderson/Documents/GitHub/ML-Workshop-Hawaii-Regime-Shifts/scaledData.csv")%>%
  select(-X)

analysisData$nextYearRegime<-as.factor(analysisData$nextYearRegime)
scaledData$nextYearRegime<-as.factor(scaledData$nextYearRegime)

analysisData <- analysisData%>%
    mutate(Transition = if_else(Regime!=nextYearRegime, 1, 0))
analysisData <- analysisData%>%
    select(-nextYearRegime)

scaledData <- scaledData%>%
    mutate(Transition = if_else(Regime!=nextYearRegime, 1, 0))
scaledData <- scaledData%>%
    select(-nextYearRegime)

analysisData$Transition<-as.factor(analysisData$Transition)
scaledData$Transition<-as.factor(scaledData$Transition)

deepScaledData <- scaledData

Create task and split data
- Set seed here for reproducible results

task = TaskClassif$new(id = "analysisData", backend = analysisData, target = "Transition")
NNtask = TaskClassif$new(id = "scaledData", backend = scaledData, target = "Transition")
split = partition(task, ratio = 0.8)

set.seed(42)

2 K Nearest Neighbors Modeling

Define the KNN learner

KNNlearner = lrn("classif.kknn", 
                 k = 5, 
                 predict_type = "prob" 
                 )

Note: Predict type is “response” when using general classification, and “prob” for Boolean predictions. We would also need to have the classification score as .auc. Usually use .acc for an area under the curve evaluation metric.

Train the KNN model

KNNlearner$train(task, row_ids = split$train)

Predict on testing data and generate results

KNNprediction = KNNlearner$predict(task, row_ids = split$test)
print(KNNprediction$score(msr("classif.auc")))

## classif.auc 
##   0.4333333

autoplot(KNNprediction)

KNNprediction$confusion

##         truth
## response  0  1
##        0 12  9
##        1  3  1

3 Random Forest Models

3.1 Random Forest Model using randomForest library:

Train the random forest

randomForestModel <- randomForest(Transition  ~ ., 
                   data = analysisData, 
                   ntree = 100,
                   mtry = 2,
                   max_features = 1,
                   importance=TRUE)
randomForestModel

## 
## Call:
##  randomForest(formula = Transition ~ ., data = analysisData, ntree = 100,      mtry = 2, max_features = 1, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 33.6%
## Confusion matrix:
##    0  1 class.error
## 0 79 13   0.1413043
## 1 29  4   0.8787879

Generate a variable importance plots to show which variable contributed most to the model’s predictions using randomForest function.

varImpPlot(randomForestModel)

3.2 Random Forest Model Using MLR3

Define the Ranger learner using the MLR3’s Random Forest model. The hyper parameters are set to “to_tune” since we later tune them using the instancing optimization method.

Rangerlearner = lrn("classif.ranger", 
                    num.trees = to_tune(10,500),
                    mtry = to_tune(1,10),
                    max.depth = to_tune(10,20),
                    importance = "permutation", 
                    predict_type = "prob"
                    )

Tune Hyper Parameters –> https://mlr3tuning.mlr-org.com/ We used an instance method to run through different combinations of hyper parameters and then assigned the most optimal to the learner.

Rangerinstance = ti(
  task = task$clone()$filter(split$train),
  learner = Rangerlearner,
  resampling = rsmp("cv", folds = 5),
  measures = msr("classif.auc"),
  terminator = trm("evals", n_evals = 25)
  )
# You can also use terminator = trm("run_time", secs = XX)
# However, this may lead to more variability due to the internal computer clock and background processes and could vary between computers. 

Rangertuner = tnr("random_search", batch_size=20) 
# Specify searching method, "random_search". Batch size is how many times to try a grouping of parameters before moving on to different groupings

# Run the optimization and set results to the learner:
Rangertuner$optimize(Rangerinstance)
Rangerlearner$param_set$values=Rangerinstance$result_learner_param_vals

Train the random forest model

Rangerlearner$train(task, row_ids = split$train)

Predict on test set

Rangerprediction = Rangerlearner$predict(task, row_ids = split$test)
print(Rangerprediction$score(msr("classif.auc")))

## classif.auc 
##   0.6533333

autoplot(Rangerprediction, type="roc")

Rangerprediction$confusion

##         truth
## response  0  1
##        0 13  9
##        1  2  1

Generate a variable importance plots to show which variable contributed most to the model’s predictions using MLR3 functions.

# Get variable importance as named numeric vector
imp = Rangerlearner$importance()

# Convert to data frame
imp_df = data.frame(
  variable = names(imp),
  importance = as.numeric(imp)
)

# Sort by importance
imp_df = imp_df[order(imp_df$importance, decreasing = TRUE), ]

# Plot with ggplot2
ggplot(imp_df, aes(x = reorder(variable, importance), y = importance)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Variable Importance (Permutation)", x = "Variable", y = "Importance") +
  theme_minimal()

4 Gradient Boosting Model

Model Setup

GBlearner = lrn("classif.xgboost",
                nrounds = 500,
                eta = to_tune(0.01, 0.3),
                max_depth = to_tune(3, 10),
                min_child_weight = to_tune(1, 10),
                lambda = to_tune(0, 1),
                subsample = to_tune(0.5, 1),
                colsample_bytree = to_tune(0.5, 1),
                predict_type = "prob",
                booster = "gbtree"
                )

GBinstance = ti(
  task = task$clone()$filter(split$train),
  learner = GBlearner,
  resampling = rsmp("cv", folds = 5),
  measures = msr("classif.auc"),
  terminator = trm("run_time", secs=10)
)

GBtuner = tnr("grid_search")
GBtuner$optimize(GBinstance)
GBlearner$param_set$values = GBinstance$result_learner_param_vals

GBlearner$train(task, row_ids = split$train)

GB Model Results

GBprediction = GBlearner$predict(task, row_ids = split$test)
print(GBprediction$score(msr("classif.auc")))

## classif.auc 
##   0.8733333

autoplot(GBprediction)

GBprediction$confusion

##         truth
## response  0  1
##        0 15 10
##        1  0  0

# Get variable importance as named numeric vector
imp = GBlearner$importance()

# Convert to data frame
imp_df = data.frame(
  variable = names(imp),
  importance = as.numeric(imp)
)

# Sort by importance
imp_df = imp_df[order(imp_df$importance, decreasing = TRUE), ]

# Plot with ggplot2
ggplot(imp_df, aes(x = reorder(variable, importance), y = importance)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Variable Importance (Permutation)", x = "Variable", y = "Importance") +
  theme_minimal()

5 Neural Network Model

Note: This is a mathematical model, so scaled data would be more beneficial here as opposed to decision tree models such as randomForest or gradientBoosted. This is why this model calls a unique task, NNtask, which pulls from the scaled data.

NNlearner = lrn("classif.nnet",
                predict_type = "prob",
                decay = to_tune(10,200),
                size = to_tune(15,30)
                )
# Size is the neuron hyper parameter and must be greater than the amount of inputs going into the data

NNinstance = TuningInstanceSingleCrit$new(
  task = NNtask$clone()$filter(split$train),
  learner = NNlearner,
  resampling = rsmp("cv", folds = 5),
  measure = msr("classif.auc"),
  terminator = trm("evals", n_evals = 100) 
  )
NNtuner = tnr("random_search", batch_size=5) 
NNtuner$optimize(NNinstance)
NNlearner$param_set$values=NNinstance$result_learner_param_vals

NNlearner$train(NNtask, row_ids = split$train)

## # weights:  331
## initial  value 599.358230 
## iter  10 value 74.660774
## iter  20 value 55.775589
## iter  30 value 55.766230
## final  value 55.766202 
## converged

Neural Network Results

NNprediction= NNlearner$predict(NNtask, row_ids = split$test)
print(NNprediction$score(msr("classif.auc")))

## classif.auc 
##        0.58

autoplot(NNprediction, type="roc")

NNprediction$confusion

##         truth
## response  0  1
##        0 15 10
##        1  0  0

autoplot(NNprediction)

6 Deep Neural Network Model

Preparing the Data using one-hot encoding

# Create a classification task using scaledData
deepNNTask = TaskClassif$new("deepScaledData", backend = deepScaledData, target = "Transition")

# Define one-hot encoder pipeline operator
poe = po("encode", method = "one-hot")

# Train the pipeline (apply one-hot encoding)
encoded_task = poe$train(list(deepNNTask))[[1]]

Define the deep NN learner using “LearnerTorchMLP” and specify use for classification

deepNNlearner = LearnerTorchMLP$new("classif")
deepNNlearner$param_set$set_values(
  epochs = 100,
  device = "cpu",
  neurons = c(10,10),
  batch_size = 20
  )
deepNNlearner$predict_type = "prob"

set.seed(42)
split = partition(task, ratio = 0.8)

deepNNlearner$train(encoded_task, row_ids = split$train)

Deep NN Model Predictions

deepNNprediction = deepNNlearner$predict(encoded_task, row_ids = split$test)
print(deepNNprediction$score(msr("classif.auc")))

## classif.auc 
##   0.6597222

autoplot(deepNNprediction)

7 Results

After analyzing our data through multiple models, we came to inconclusive results. No model stood out from the rest and we believe more data would be needed for training to create more accurate models in the future.