ordreg_lasso_real_data.Rmd

---
title: 'Real Data Analysis: Ordinal Regression LASSO'
author: "Brian Chen"
output: html_document
date: "`r format(Sys.time(), '%m/%d/%Y')`"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(data.table)
library(devtools)
library(caret)
library(dplyr)
library(tictoc)
library(psych)
library(cvms)
load_all()

```

# Read in Data and Preprocess Data for Lasso

```{r read_scale}
#Read in Data
data <- fread('processed_data.csv')

#Create Key to make ordered factor and rename
#Reverse order to make interpretation easier
y_factor_key <- data.frame(original = sort(unique(data$`B1:SATISFACTION W/LIFE`), decreasing = F),
                           y = ordered(1:5))

#factor(1:length(unique(data$`B1:SATISFACTION W/LIFE`)), ordered = T )

#Merge and drop Case ID
data_ordered_y <- merge(y_factor_key, data, by.x = "original", by.y = "B1:SATISFACTION W/LIFE") %>% 
  arrange(`W1 Case Id`) %>%
  select(-c("original", "W1 Case Id"))

#Scale data for LASSO
data_ordered_y_scaled <- data_ordered_y
data_ordered_y_scaled[,2:199] <- as.data.frame(scale(data_ordered_y[,2:199]))

```

# Create Training and held out data / folds within training data for CV

```{r split_data}

set.seed(5)

#Create training and testing set
train_test_i <- createDataPartition(y = data_ordered_y_scaled$y, p = 0.6, list = F)

data_train <- data_ordered_y_scaled[train_test_i,]
data_heldout <- data_ordered_y_scaled[-train_test_i,]

#Within training set, create folds
k <- 5 #set number of folds 
fold_indicies <- createFolds(y = data_train$y, k=k)

```

# K fold cross validation to find best lambda

```{r cv} 

lambda_grid <- seq(0.005, 0.05, 0.005)

lambda_info <- list()

for (fold in 1:k){
  
  #Fold Training and testing set
  fold_train <- data_train[-fold_indicies[[fold]],]
  fold_test <- data_train[fold_indicies[[fold]],]
  
  #Fit Ordinal LASSO
  tic(paste0("Fold #", fold))
  res.ordreg <- suppressWarnings(ordreg.lasso(
    formula = y ~ .,
    data = fold_train,
    lambdas = lambda_grid
  ))
  toc()
  
  #Vector of Testing Kappas
  test_kappas <- rep(NA, length(res.ordreg$lambdas))
  
  #For each value of lambda of the fold, get kappa 
  for (i in 1:length(res.ordreg$lambdas)){
    
    #Predict
    predicted <- ordreg.predict(alpha = res.ordreg$alpha[i,], 
                                beta = res.ordreg$beta[i,], 
                                x = as.matrix(fold_test[,-1]))
    
    #Calculate Kappa
    test_kappas[i] <- cohen.kappa(x = matrix(c(fold_test[,1], predicted),
                                             ncol = 2, 
                                             byrow = F))$weighted.kappa
    
  }
  
  #Store lambda info
  lambda_info[[fold]] <- data.frame(fold = rep(fold, length(res.ordreg$lambdas)),
                                    lambda = res.ordreg$lambdas,
                                    n.nonzero = res.ordreg$n.nonzero,
                                    testing_kappa = test_kappas)
  
}


```
# Plots to Heruistically Select Lambda

```{r plot_lambda}

lambda_info_df <- do.call("rbind", lambda_info) %>%
  mutate(lambda = as.factor(lambda),
         fold = as.factor(fold))

lambda_avg <- lambda_info_df %>%
  group_by(lambda) %>%
  summarise(avg_testing_kappa = mean(testing_kappa),
            avg_n.nonzero = mean(n.nonzero))

p1 <- ggplot(lambda_info_df, aes(x = n.nonzero, y = testing_kappa, group = fold)) +
  geom_point(aes(color = fold)) +
  geom_line(aes(color = fold)) +
  labs(title = "No. of Nonzero Parameters vs. Testing Kappa by Fold",
       x = "No. of Nonzero Parameters",
       y = "Testing Kappa",
       color = "Fold")

p2 <- ggplot(data=lambda_avg, aes(x = avg_n.nonzero, y = avg_testing_kappa, label = lambda)) +
  geom_line() + 
  geom_point() +
  geom_label(alpha = 0, label.size = NA, nudge_y = 0.003, nudge_x = -5) +
  labs(title = "Average No. of Nonzero Parameters vs. Average Testing Kappa across Fold",
       x = "Average No. of Nonzero Parameters",
       y = "Average Testing Kappa")

print(p1)
print(p2)

```

# Fit LASSO with $\lambda = 0.02$ for variable selection 

```{r variable_selection}

#Ordinal regression w/ lambda = 0.02
tic("Variable selection model")
ordreg.vs <- ordreg.lasso(
    formula = y ~ .,
    data = data_train,
    lambdas = 0.02)
toc()

#Column names of nonzero Betas
nonzero_beta <- colnames(data_train[2:199])[ordreg.vs$beta != 0 ]
nonzero_beta

predicted_heldout <- ordreg.predict(alpha = ordreg.vs$alpha[1,],
                                    beta = ordreg.vs$beta[1,],
                                    x = as.matrix(data_heldout[,-1]))

kappa_heldout <- cohen.kappa(x = matrix(c(data_heldout[,1], predicted_heldout),
                                             ncol = 2, 
                                             byrow = F))$weighted.kappa

print(paste0("Testing Kappa of the model in held out data is ", round(kappa_heldout, digits = 3)))

# creating and plotting a confusion matrix
# save actual y and predicted y into a tibble
y = data_heldout[,1]
y.pred = predicted_heldout
d_multi<-tibble("target"=y,"prediction"=y.pred)
conf_mat<-confusion_matrix(targets = d_multi$target,
                             predictions = d_multi$prediction)
plot_confusion_matrix(conf_mat$`Confusion Matrix`[[1]],add_normalized=F)

```


# Fit oridnal regression on held out data with selected variables (no penalty)

```{r held_out_fit}

data_heldout_filter <- data_heldout[,c("y", nonzero_beta)] 

tic("Model on Held Out data")
ordreg.final <- ordreg.lasso(
    formula = y ~ .,
    data = data_heldout_filter,
    lambdas = 0,
    return.cov = TRUE)
toc()

```

# Testing for significance

```{r inference}

signif_tbl <- data.frame(Question = colnames(data_heldout_filter[-1]),
                         Estimate = unname(t(ordreg.final$beta)),
                         SE = unname(sqrt(diag(ordreg.final$cov[5:43,5:43]))))

signif_tbl <- signif_tbl %>%
  mutate(tvalue = Estimate / SE) %>%
  mutate(p = 2*pt(-abs(tvalue), 1434 - 43 - 1)) %>%
  arrange(p)

sum(signif_tbl$p < 0.05)

signif_tbl
```


```{r export, include =FALSE}
#Export workspace so we don't have to rerun analysis for report 
#save.image("C:/Users/brian/Documents/Grad School/BIOS 735/HOOVS/ordreg_lasso_real_data.RData")

```