01-countrywide-analysis.Rmd

---
title: "Countrywide Analysis"
output: html_document
date: "2023-12-29"
---

This markdown follows monograph section 7.3 Countrywide Model Results

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

Package setup

```{r packages, echo=FALSE, include=FALSE}
# Function to install and load packages
install_and_load <- function(package) {
  if (!require(package, character.only = TRUE)) {
    install.packages(package)
    library(package, character.only = TRUE)
  }
}

# List of packages
packages <- c("data.table", "statmod", "ggplot2", "tweedie", 
              "glmnet", "HDtweedie", "plotly", "arrow", "knitr")

# Install and load all packages
sapply(packages, install_and_load)
```

Load helper functions to create various graphs

```{r echo=FALSE, include=FALSE}
source( "helper_functions.R")
```

Load the modeling data. Set your working directlry as needed. 

```{r pressure, echo=FALSE, include=FALSE}
getwd()
# setwd()
modeling_data <- read_parquet("output_file.parquet")

```

Here, we will build a lasso model on the countrywide modeling data

```{r prep-data, echo=FALSE, include = FALSE}

# filter to training rows only
training_data <- modeling_data[stratification %in% c(1, 2, 3, 4), ]

# select modeling variables for Lasso
# The HDtweedie packages uses a different format for predictors, but the included predictors are identical to our GLM
lasso_modeling_variables <- c(
  "driver_age_18_38_hinge", "driver_age_38_76_hinge", "driver_age_76_99_hinge",
  "vehicle_age_0_10_hinge", "vehicle_age_10_99_hinge",
  "multi_yes", "xTreme_yes",
  "weight_heavy", "weight_light", "weight_extra_light",
  "ind_construction", "ind_farming", "ind_finance_and_insurance",
  "ind_fine_arts", "ind_fireworks", "ind_food_services",
  "ind_health_care", "ind_real_estate", "ind_retail",
  "ind_large_state", "ind_medium_state",
  "ind_small_state", "ind_small_state_cw_dist"
)

# HDTweedie requires a matrix format for predictor variables. First, we will subset to all of the columns that we need.
train_cols_for_matrix <- modeling_data[stratification %in% c(1, 2, 3, 4), lasso_modeling_variables, with = FALSE]
validate_cols_for_matrix <- modeling_data[stratification %in% c(5), lasso_modeling_variables, with = FALSE]

# Reformatting to a matrix for HDtweedie so that Lasso Tweedie can work.
training_data_matrix <- as.matrix(train_cols_for_matrix, rownames = FALSE)
validate_data_matrix <- as.matrix(validate_cols_for_matrix, rownames = FALSE)

# This is the data format needed for creating our loss vector and CV vector below
training_data <- modeling_data[stratification %in% c(1, 2, 3, 4), ]
validate_data <- modeling_data[stratification %in% c(5), ]

# Losses and CV need to be a vector
loss_vector <- as.vector(training_data$incurred_loss)
cv_vector <- as.vector(training_data$stratification)

# As exposures are all 1, we will create this vector.
weights <- rep(1, length(drop(loss_vector)))

```


We will fit a GLM to our modeing data. This was completed in the earlier markdown script, but we are fitting it again as only the data was saved from the prior file. This will not take a long time. 

```{r glm-fit}

# run our initial GLM model
glm_model <- glm(
  formula = incurred_loss ~
    driver_age_18_38_hinge + driver_age_38_76_hinge + driver_age_76_99_hinge +
    vehicle_age_0_10_hinge + vehicle_age_10_99_hinge +
    C(MultiPolicy) + C(xTreme_TurnSignal) + C(car_weight) + C(industry_code) + C(subset),
  family = tweedie(1.6, link.power = 0),
  offset = log(exposure),
  data = training_data
)

```


Fit a lasso model to the CW data. This may take some time, so we will save this model to avoid needing to re-create it. 


```{r}

# after the initial run, set recompute to false. 
recompute <- TRUE
# recompute <- FALSE

if (recompute) {
  # fit the lasso model
  lasso_model <- cv.HDtweedie(
    x = training_data_matrix,
    y = loss_vector, group = NULL, p = 1.6, weights, lambda = NULL,
    pred.loss = "deviance",
    nfolds = 4, foldid = cv_vector,
    standardize = TRUE
  )

  saveRDS(lasso_model, file = "lasso_model_on_training_data_matrix.rds")
} else {
  lasso_model <- readRDS("lasso_model_on_training_data_matrix.rds")
}
```

Let's see CV plot

```{r analyse-error}

# Preparing the data for plotting
data <- data.frame(
  lambda = lasso_model$lambda,
  cvm = lasso_model$cvm,
  cvsd = lasso_model$cvsd
)
data$cvlower = data$cvm - data$cvsd
data$cvupper = data$cvm + data$cvsd

# Plotting
ggplot(data, aes(x = lambda, y = cvm)) +
  geom_line() +
  geom_ribbon(aes(ymin = cvlower, ymax = cvupper), alpha = 0.2) +
  labs(x = "Lambda", y = "Mean Cross-Validated Error") +
  theme_minimal()
```

Add the GLM and Lasso predictions to the dataset as new columns. 

```{r}
coef(lasso_model, s = c("lambda.min")) # NB if you don't tell it to pick the best lambda, it won't.
coef(lasso_model) # Please note that this does NOT give you the best model.

# get the coefficient table to compare GLM and Lasso predictions
cw_coeff_table <- create_coeff_table(glm_model, lasso_model)

# Base Model Coefficient Table
print(cw_coeff_table)
write.csv(cw_coeff_table, "../Graphs/FullData/CSV/cw_coeff_table.csv")
summary(glm_model)

# add the predictions to our training data
training_data[, glm_prediction := predict.glm(glm_model, training_data, type = "response")]
training_data[, lasso_prediction := predict(lasso_model, training_data_matrix, s = "lambda.min")]
```

Produce relativity plots and prediction plots for all risk characteristics. 

```{r}
if(TRUE){
# Full Model Prediction Plots and Relativity Plots
# Figure 7.8
relativity_plot(get_cat_var_table(training_data, cw_coeff_table, "industry_code"))
# Figure 7.9
prediction_plot(get_cat_var_table(training_data, cw_coeff_table, "industry_code"))
relativity_plot(get_cat_var_table(training_data, cw_coeff_table, "MultiPolicy"))
prediction_plot(get_cat_var_table(training_data, cw_coeff_table, "MultiPolicy"))
relativity_plot(get_cat_var_table(training_data, cw_coeff_table, "xTreme_TurnSignal"))
prediction_plot(get_cat_var_table(training_data, cw_coeff_table, "xTreme_TurnSignal"))
relativity_plot(get_cat_var_table(training_data, cw_coeff_table, "car_weight"))
prediction_plot(get_cat_var_table(training_data, cw_coeff_table, "car_weight"))

# Get tables for continuous variables.
driver_age_table <- get_driver_var_table(training_data, cw_coeff_table, variable_name = "driver_age")
vehicle_age_table <- get_vehicle_var_table(training_data, cw_coeff_table, variable_name = "vehicle_age")

# Full Model Continuous Prediction Plots and Relativity Plots
relativity_plot(driver_age_table)
prediction_plot(driver_age_table)
relativity_plot(vehicle_age_table)
prediction_plot(vehicle_age_table)

# Save the tables as csv
write.csv(cw_coeff_table, file = "../Graphs/FullData/CSV/cw_coeff_table.csv", row.names = FALSE)
write.csv(get_cat_var_table(training_data, cw_coeff_table, "industry_code"), file = "../Graphs/FullData/CSV/industry_code_table.csv", row.names = FALSE)
write.csv(get_cat_var_table(training_data, cw_coeff_table, "MultiPolicy"), file = "../Graphs/FullData/CSV/multipolicy_table.csv", row.names = FALSE)
write.csv(get_cat_var_table(training_data, cw_coeff_table, "xTreme_TurnSignal"), file = "../Graphs/FullData/CSV/xtreme_turnsignal_table.csv", row.names = FALSE)
write.csv(get_cat_var_table(training_data, cw_coeff_table, "car_weight"), file = "../Graphs/FullData/CSV/car_weight_table.csv", row.names = FALSE)
write.csv(get_driver_var_table(training_data, cw_coeff_table, variable_name = "driver_age"), file = "../Graphs/FullData/CSV/driver_age_table.csv", row.names = FALSE)
write.csv(get_vehicle_var_table(training_data, cw_coeff_table, variable_name = "vehicle_age"), file = "../Graphs/FullData/CSV/vehicle_age_table.csv", row.names = FALSE)
}
```

Create prediction plots and lift charts on validate data. Note that we compare both to experienced relativities as well as "true" relativities. 

```{r}

# Now, we'll see how well we do on validate data. Here, we will only look at predictions as the
# relativities will give us no new information

validate_data[, glm_prediction := predict.glm(glm_model, validate_data, type = "response")]
validate_data[, lasso_prediction := predict(lasso_model, validate_data_matrix, s = "lambda.min")]

prediction_plot(get_cat_var_table(validate_data, cw_coeff_table, "industry_code"))
prediction_plot(get_cat_var_table(validate_data, cw_coeff_table, "MultiPolicy"))
prediction_plot(get_cat_var_table(validate_data, cw_coeff_table, "xTreme_TurnSignal"))
prediction_plot(get_cat_var_table(validate_data, cw_coeff_table, "car_weight"))

driver_age_table <- get_driver_var_table(validate_data, cw_coeff_table, variable_name = "driver_age")
vehicle_age_table <- get_vehicle_var_table(validate_data, cw_coeff_table, variable_name = "vehicle_age")

prediction_plot(driver_age_table)
prediction_plot(vehicle_age_table)

# Lift Charts for Full Dataset
# Figure 7.11
double_lift_chart(
  dataset = training_data,
  lasso_credibility = "lasso_prediction",
  glm = "glm_prediction",
  actual = "true_risk",
  normalize = TRUE
)

double_lift_chart(
  dataset = training_data,
  lasso_credibility = "lasso_prediction",
  glm = "glm_prediction",
  actual = "incurred_loss",
  normalize = TRUE
)

selected_columns <- training_data[, c("lasso_prediction", "glm_prediction", "true_risk", "incurred_loss")]
write.csv(selected_columns, file = "../Graphs/FullData/CSV/training_data_predictions.csv", row.names = FALSE)

# Figure 7.10
double_lift_chart(
  dataset = validate_data,
  lasso_credibility = "lasso_prediction",
  glm = "glm_prediction",
  actual = "incurred_loss",
  normalize = TRUE
)

double_lift_chart(
  dataset = validate_data,
  lasso_credibility = "lasso_prediction",
  glm = "glm_prediction",
  actual = "true_risk",
  normalize = TRUE
)

selected_columns <- validate_data[, c("lasso_prediction", "glm_prediction", "true_risk", "incurred_loss")]
write.csv(selected_columns, file = "../Graphs/FullData/CSV/validation_data_predictions.csv", row.names = FALSE)

```

Do not clear your working environment after this script to ensure that all data is available for the next script.