```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = 'G:\\Il mio Drive\\Teaching\\Data Science Lab\\Lecture 3')
```

In some scenarios we may be interested in building a statistical model to predict an outcome. In this case, for instance, we may want to use different models to predict the location of a firm. The Compustat data entail US and Canadian enterprises. In the next chunks of code I will build five different models (logistic regression, CART, Conditional Inference Tree, Random Forest, Bayesian Additive Regression Trees) to predict the location of the firm.

Moreover, I will provide details on the most widely used performance measures in the case of a classification problem.

```{r, warning=FALSE}
library(readxl)
library(caret)
library(randomForest)
library(bartMachine)
library(PRROC)
library(rpart)
library(party)
```

First things first, let's upload the Compustat data and perform a naive trimming of the data, excluding all the missing observations.

```{r}
data <- read_excel("G:\\Il mio Drive\\Econometrics Lab\\Data\\Compustat Data.xlsx")
data <- data[, !names(data) %in% c("Interest Expense - Total (Financial Services)",
                            "Net Interest Income", "Nonperforming Assets - Total")]
data_clean <- na.omit(data)
```

Before running the analyses, I restict the set of predictors to the following variables and I create a dummy variable that assumes value 1 if the firm is located in the US and 0 if it is located in Canada.

```{r}
myvariables <- c("ISO Currency Code",
                 "Assets - Total", "Average Short-Term Borrowings",
                 "Current Assets - Total", "Long-Term Debt Due in One Year",
                 "Debt in Current Liabilities - Total", "Employees",
                 "Earnings Before Interest and Taxes", "Liabilities - Total",
                 "Net Income (Loss)", "In Process R&D Expense",
                 "GIC Sectors", "Standard Industry Classification Code")
data_prediction <- data_clean[myvariables]
data_prediction$iso_code <- ifelse(data_prediction$`ISO Currency Code`=="USD", 0, 1)
data_prediction <- data_prediction[, !names(data_prediction) %in% c("ISO Currency Code")]
```

In order to check how good are the five models, I randomly split the data into two disjoint sets: a training set that I will use to build the model and a test set that I will use to validate the quality of the model's prediction.

```{r}
set.seed(123)
index <- sample(seq_len(nrow(data_prediction)),
                    size = nrow(data_prediction)*0.5) 

train <- data_prediction[index,]
test <- data_prediction[-index,]
```

Moreover, I am renaming the variables in the dataset and constructing the "formula" that I will use for all the predictive models that I will run.

```{r}
colnames(train) <- c("assets", "short_term_borrow",
                     "current_assets", "debt",
                     "debt_liabilities", "employees", 
                     "EBIT", "liabilities",
                     "net_income", "r_d",
                     "gic", "SICC", "iso_code")
colnames(test) <- c("assets", "short_term_borrow",
                     "current_assets", "debt",
                     "debt_liabilities", "employees", 
                     "EBIT", "liabilities",
                     "net_income", "r_d",
                     "gic", "SICC", "iso_code")
predictors <- c("assets", "short_term_borrow",
                     "current_assets", "debt",
                     "debt_liabilities", "employees", 
                     "EBIT", "liabilities",
                     "net_income", "r_d",
                     "gic", "SICC")
formula <- as.formula(paste("as.factor(iso_code) ~",
                            paste(predictors, collapse="+")))
formula
```

## Logistic Regression

The first model that I run is a logistic regression with the inclusion of all the covariates.

```{r}
logit<-glm(formula, data= train, family=binomial(link='logit'))
summary(logit)
```

To get the accuracy of the model I first get the predicted probabilities, then impute the values for the outcome variable.

```{r}
# Accurancy from cv
fitted.results.logit <- predict(logit, newdata = test, type='response')
fitted.logit <- ifelse(fitted.results.logit >= 0.5, 1, 0)
head(fitted.logit)
```

Once I get the predicted (or fitted values) for this model, I can evaluate its performance using a number of different performance measures. Below the functions to compute the F-1 Score and the Balanced Accuracy.

```{r}
## F1- Score
# predicted: vector of predicted values
# expected: vector of observed value
# positive.class: class of binary predictions we are mostly interested in (e.g., "1", "0")

f1_score <- function(predicted, expected, positive.class) {
  
  # Generate Confusion Matrix
  c.matrix = as.matrix(table(expected, predicted))
  
  # Compute Precision
  precision <- diag(c.matrix) / colSums(c.matrix)
  
  # Compute Recall
  recall <- diag(c.matrix) / rowSums(c.matrix)
  
  # Compute F-1 Score
  f1 <-  ifelse(precision + recall == 0, 0, 2*precision*recall/(precision + recall))
  
  # Extract F1-score for the pre-defined "positive class"
  f1 <- f1[positive.class]
  
  # Assuming that F1 is zero when it's not possible compute it
  f1[is.na(f1)] <- 0
  
  # Return F1-score
  return(f1)
}

## Balanced Accuracy (BACC)
# predicted: vector of predicted values
# expected: vector of observed value

balanced_accuracy <- function(predicted, expected) {
  
  # Generate Confusion Matrix
  c.matrix = as.matrix(table(predicted, expected))
  
  # First Row Generation
  first.row <- c.matrix[1,1] / (c.matrix[1,1] + c.matrix[1,2])  
  
  # Second Row Generation
  second.row <- c.matrix[2,2] / (c.matrix[2,1] + c.matrix[2,2])  
  
  # # "Balanced" proportion correct (you can use different weighting if needed)
  acc <- (first.row + second.row)/2 
  
  # Return Balanced Accuracy
  return(acc)
}
```

```{r}
# RMSE
caret::postResample(fitted.logit, test$iso_code)
# For good predictive model the MAE and RMSE values should be low 


# Confusion Matrix
confusionMatrix(data = as.factor(fitted.logit),
                reference = as.factor(test$iso_code))

# Balanced Accuracy
balanced_accuracy_logit<-balanced_accuracy(fitted.logit, test$iso_code)
balanced_accuracy_logit

# F1-Score
f1_logit_1 <- f1_score(fitted.logit,
                     test$iso_code,
                     positive.class="1")
f1_logit_1
f1_logit_0 <- f1_score(fitted.logit,
                     test$iso_code,
                     positive.class="0")
f1_logit_0

# ROC Curve and PR- Curve
fg.logit <- fitted.logit[test$iso_code==1]
bg.logit <- fitted.logit[test$iso_code==0]
roc_logit <- roc.curve(scores.class0 = fg.logit,
                       scores.class1 = bg.logit,
                       curve = T)
plot(roc_logit)
pr_logit <- pr.curve(scores.class0 = fg.logit,
                     scores.class1 = bg.logit,
                     curve = T)
plot(pr_logit) 
```

## Classification and Regression Tree

The second model that I run is a classification and regression tree from the "rpart" package in R.

```{r}
rpart <- rpart(formula, data=train, method="class") 
printcp(rpart) # display the results 
plotcp(rpart) # visualize cross-validation results 
# summary(rpart) # detailed summary of splits
```

You can depict the classification tree by using the "plot()" function.

```{r}
#Plot tree 
rpart.plot::rpart.plot(rpart)
```

To get the accuracy of the model I first get the predicted probabilities, then impute the values for the outcome variable.

```{r}
fitted.results.rpart <- predict(rpart, newdata=test,type='prob')
fitted.rpart <- ifelse(fitted.results.rpart[,2] >= 0.5, 1, 0)
```

Below, I depict the predictive performance of the model.

```{r}
# RMSE
caret::postResample(fitted.rpart, test$iso_code)
# For good predictive model the MAE and RMSE values should be low 


# Confusion Matrix
confusionMatrix(data = as.factor(fitted.rpart),
                reference = as.factor(test$iso_code))

# Balanced Accuracy
balanced_accuracy_rpart<-balanced_accuracy(fitted.rpart, test$iso_code)
balanced_accuracy_rpart

# F1-Score
f1_rpart_1 <- f1_score(fitted.rpart,
                       test$iso_code,
                       positive.class="1")
f1_rpart_1
f1_rpart_0 <- f1_score(fitted.rpart,
                       test$iso_code,
                       positive.class="0")
f1_rpart_0

# ROC Curve and PR- Curve
fg.rpart <- fitted.rpart[test$iso_code==1]
bg.rpart <- fitted.rpart[test$iso_code==0]
roc_rpart <- roc.curve(scores.class0 = fg.rpart,
                       scores.class1 = bg.rpart,
                       curve = T)
plot(roc_rpart)
pr_rpart <- pr.curve(scores.class0 = fg.rpart,
                     scores.class1 = bg.rpart,
                     curve = T)
plot(pr_rpart) 
```

## Conditional Inference Tree

One potential drawback of the classification and regression trees 

```{r}
c.tree <- ctree(formula, data=train,
               control = ctree_control(testtype = "MonteCarlo",
               mincriterion = 0.99, nresample = 1000))
```

You can plot the tree by running the following chunk of code.

```{r}
plot(c.tree, gp = gpar(fontsize = 6))
```

To get the accuracy of the model I first get the predicted probabilities, then impute the values for the outcome variable.

```{r}
fitted.results.tree <- as.matrix(unlist(predict(c.tree,
                                 newdata = test, type='prob')))
fitted.prob.tree <- fitted.results.tree[seq_along(fitted.results.tree) %%2 == 0]
fitted.tree <- ifelse(fitted.prob.tree >= 0.5, 1, 0)
```

Below, I depict the predictive performance of the model.

```{r}
# RMSE
caret::postResample(fitted.tree, test$iso_code)
# For good predictive model the MAE and RMSE values should be low 


# Confusion Matrix
confusionMatrix(data = as.factor(fitted.tree),
                reference = as.factor(test$iso_code))

# Balanced Accuracy
balanced_accuracy_tree<-balanced_accuracy(fitted.tree, test$iso_code)
balanced_accuracy_tree

# F1-Score
f1_tree_1 <- f1_score(fitted.tree,
                       test$iso_code,
                       positive.class="1")
f1_tree_1
f1_tree_0 <- f1_score(fitted.tree,
                       test$iso_code,
                       positive.class="0")
f1_tree_0

# ROC Curve and PR- Curve
fg.tree <- fitted.tree[test$iso_code==1]
bg.tree <- fitted.tree[test$iso_code==0]
roc_tree <- roc.curve(scores.class0 = fg.tree,
                       scores.class1 = bg.tree,
                       curve = T)
plot(roc_tree)
pr_tree <- pr.curve(scores.class0 = fg.tree,
                     scores.class1 = bg.tree,
                     curve = T)
plot(pr_tree) 
```

## Random Forest

The last model that I build is a random forest from the "randomForest" package.

```{r}
set.seed(133234)
rf <- randomForest(formula, data=train, importance=TRUE, ntree=200)
```

```{r}
print(rf)
plot(rf)
varImpPlot(rf)
```

To get the accuracy of the model I first get the predicted probabilities, then impute the values for the outcome variable.

```{r}
fitted.rf <- predict(rf, test)
fitted.rf <- as.numeric(matrix(fitted.rf))
```

Below, I depict the predictive performance of the model.

```{r}
# RMSE
caret::postResample(fitted.rf, test$iso_code)
# For good predictive model the MAE and RMSE values should be low 


# Confusion Matrix
confusionMatrix(data = as.factor(fitted.rf),
                reference = as.factor(test$iso_code))

# Balanced Accuracy
balanced_accuracy_rf<-balanced_accuracy(fitted.rf, test$iso_code)
balanced_accuracy_rf

# F1-Score
f1_rf_1 <- f1_score(fitted.rf,
                       test$iso_code,
                       positive.class="1")
f1_rf_1
f1_rf_0 <- f1_score(fitted.rf,
                       test$iso_code,
                       positive.class="0")
f1_rf_0

# ROC Curve and PR- Curve
fg.rf <- fitted.rf[test$iso_code==1]
bg.rf <- fitted.rf[test$iso_code==0]
roc_rf <- roc.curve(scores.class0 = fg.rf,
                       scores.class1 = bg.rf,
                       curve = T)
plot(roc_rf)
pr_rf <- pr.curve(scores.class0 = fg.rf,
                     scores.class1 = bg.rf,
                     curve = T)
plot(pr_rf) 
```

## Bayesian Forest (Bayesian Additive Regression Trees)

```{r}
set.seed(133234)
bart_machine <- bartMachine(X = as.data.frame(train[predictors]),
                          y = as.factor(train$iso_code),
                          use_missing_data=FALSE)
```


To get the accuracy of the model I first get the predicted probabilities, then impute the values for the outcome variable.

```{r}
fitted.prob.bart <- 1- round(predict(bart_machine,
                                        as.data.frame(test[predictors]),
                                        type='prob'), 6)
fitted.bart <- ifelse(fitted.prob.bart> 0.5, 1, 0)
```

```{r}
# RMSE
caret::postResample(fitted.bart, test$iso_code)
# For good predictive model the MAE and RMSE values should be low 


# Confusion Matrix
confusionMatrix(data = as.factor(fitted.bart),
                reference = as.factor(test$iso_code))

# Balanced Accuracy
balanced_accuracy_bart<-balanced_accuracy(fitted.bart, test$iso_code)
balanced_accuracy_bart

# F1-Score
f1_bart_1 <- f1_score(fitted.bart,
                       test$iso_code,
                       positive.class="1")
f1_bart_1
f1_bart_0 <- f1_score(fitted.bart,
                       test$iso_code,
                       positive.class="0")
f1_bart_0

# ROC Curve and PR- Curve
fg.bart <- fitted.bart[test$iso_code==1]
bg.bart <- fitted.bart[test$iso_code==0]
roc_bart <- roc.curve(scores.class0 = fg.bart,
                       scores.class1 = bg.bart,
                       curve = T)
plot(roc_bart)
pr_bart <- pr.curve(scores.class0 = fg.bart,
                     scores.class1 = bg.bart,
                     curve = T)
plot(pr_bart) 
```

The package bartMachine provides tools for investigation of variables' importance and variables' selection

```{r}
investigate_var_importance(bart_machine, num_replicates_for_avg = 20)
```

```{r, eval=FALSE}
vs <- var_selection_by_permute(bart_machine, 
                         num_permute_samples = 10)
```

As well as for partial dependency plots.

```{r}
pd_plot(bart_machine, j = "gic")
pd_plot(bart_machine, j = "current_assets")
```

We can incorporate information on the best predictors in the model by using a different set of priors.

```{r}
predictors
prior <- c(rep(1, times = 10), rep(2, times = 2))
```

```{r}
set.seed(133234)
bart_prior <- bartMachine(X = as.data.frame(train[predictors]),
                          cov_prior_vec = prior,
                          y = as.factor(train$iso_code),
                          use_missing_data=FALSE)
fitted.prob.bart <- 1- round(predict(bart_prior,
                             as.data.frame(test[predictors]),
                             type='prob'), 6)
fitted.prior <- ifelse(fitted.prob.bart> 0.5, 1, 0)
```

Evaluate the performance of this new model.

```{r}
# RMSE
caret::postResample(fitted.prior, test$iso_code)
# For good predictive model the MAE and RMSE values should be low 


# Confusion Matrix
confusionMatrix(data = as.factor(fitted.prior),
                reference = as.factor(test$iso_code))

# Balanced Accuracy
balanced_accuracy_prior<-balanced_accuracy(fitted.prior, test$iso_code)
balanced_accuracy_prior

# F1-Score
f1_prior_1 <- f1_score(fitted.prior,
                       test$iso_code,
                       positive.class="1")
f1_prior_1
f1_prior_0 <- f1_score(fitted.prior,
                       test$iso_code,
                       positive.class="0")
f1_prior_0

# ROC Curve and PR- Curve
fg.prior <- fitted.prior[test$iso_code==1]
bg.prior <- fitted.prior[test$iso_code==0]
roc_prior <- roc.curve(scores.class0 = fg.prior,
                       scores.class1 = bg.prior,
                       curve = T)
plot(roc_prior)
pr_prior <- pr.curve(scores.class0 = fg.prior,
                     scores.class1 = bg.prior,
                     curve = T)
plot(pr_prior) 
```