<jupyter_start><jupyter_text>Introduction to Predictive AnalyticsIn today's age of big data, insurance pricing models are becoming more dynamic. With better systems, companies are able to capture more comprehensive information about policyholders, leading to more rating factors becoming available. Regulators have also taken note of this and now require more data-driven rates. Given all this, traditional actuarial ratemaking methods may be due for an overhaul.In this workshop, we present the **Generalized Linear Model (GLM)** - a one-stop solution to generating highly interpretable and data-driven rating models. We will introduce how to build, diagnose and interpret your own GLMs in the R programming language by demonstrating a use case on sample motor insurance data.  The R statistical computing environment is one of the most comprehensive and widely used statistical modelling tools.   1. Understanding RLike most programming languages, R works through functions. These functions are stored within different libraries or **packages**. Hence, we need to load the required packages before we start to use R. In this tutorial we will use the following packages, all loaded using the `library` R command.<jupyter_code>library(MASS) #Used for Efficient Statistical Computing
library(ggplot2) #Used for Easy Data Visualization
library(readr) #Used for Faster Data Inputs and Outputs
library(dplyr) #Used for Easy Data Wrangling
library(statmod) #Used for Model Diagnostics<jupyter_output><empty_output><jupyter_text>Great, we have loaded all required packages. If you are a first-time R user, chances are that not all the above packages will be available to you at first. Hence, we may need to install certain packages into our environment. This is easy to do - just use the `install.packages("package_name")` R command. For example, to install the `MASS` package, use `install.packages("MASS")`.  2. Loading and Inspecting the DataFirst, we load the data into our environment. There are 2 datasets - *train* and *test*. The *train* dataset is used to fit the model, while the *test* dataset is used to evaluate its performance.<jupyter_code>train <- read_csv("./train.csv")
test <- read_csv("./test.csv")<jupyter_output><empty_output><jupyter_text>We will now look at a high-level summary of the *train* data.<jupyter_code>summary(train) #Summary Statistics of All Variables<jupyter_output><empty_output><jupyter_text>Next, we look at a structural summary of the *train* data. This tells us what type of variable each column in the data is.<jupyter_code>str(train) #Check Variable Types<jupyter_output><empty_output><jupyter_text>We also quickly look at the variable names in the *train* data. While the above 2 outputs already give us this, it might help sometimes to see them in a more readable manner.<jupyter_code>colnames(train) #Variable Names<jupyter_output><empty_output><jupyter_text>3. Exploratory Data Analysis (EDA)EDA is extremely important. It helps give a more clear picture of the data. Since we are fitting GLMs - models that are bound by a rigorous set of distributional assumptions - it would be useful to get an idea of how Claim Amounts and Claim Counts are distributed. In the EDA process, we therefore focus on this.We start with... 3.1 Histogram of Claim Amounts<jupyter_code>truehist(train$ClaimAmount)<jupyter_output><empty_output><jupyter_text>As expected, most policyholders have not made a claim. We do see some non-zero claims, however the above graph doesn't tell us much about them, so we'll need to look at these differently.  3.2 Histogram of Log of Claim Amounts<jupyter_code>truehist(log(train$ClaimAmount + 1)) #We add 1 to the ClaimAmount variable because log(0) doesn't exist.<jupyter_output><empty_output><jupyter_text>Pretty similar picture as the previous graph. Let's look at only the non-zero claims now. 3.3 Histogram of Non-Zero Claim Amounts<jupyter_code>truehist(train[train$ClaimAmount > 0, ]$ClaimAmount)<jupyter_output><empty_output><jupyter_text>This gives us a better idea of claims. Moving on to Claim Counts. 3.4 Bar Chart of Claim Counts<jupyter_code>truehist(train$ClaimNb)<jupyter_output><empty_output><jupyter_text>Again, as expected, most policyholders didn't file any claims, with some filing 1 and 2 claims.Great, now that we've got some idea of distributions, we can start creating some models! We look at the 2 components of aggregate claims - **Frequency** (number of claims per exposure) and **Severity** (average cost per claim). Both components are modelled separately.We start with **Frequency**. 4. Frequency ModellingSince Frequency deals with the number of claims, we need a distribution that is suited for counts. The **Poisson** distribution is an ideal candidate.Let's start with a simple model with just 2 rating factors - Driver Age Band and Vehicle Age Band. We will use the Poisson distribution with a **logarithmic link**.<jupyter_code>glm_freq1 <- glm(formula = ClaimNb ~ DrivAgeBand + VehAgeBand, data = train, family = poisson(link = "log"))
summary(glm_freq1)<jupyter_output><empty_output><jupyter_text>Let's see if this model conforms to the assumptions of a GLM. These are:1. The residuals must be normally distributed.2. The residuals must have constant variance and be centered around 03. The residuals must be random, i.e., there should be no discernable pattern to them.To check for these assumptions, we look at the **Residual Q-Q Plot** and a **Residuals vs. Fitted Values Plot**.<jupyter_code>#Generate Model Residuals
resid <- qresiduals(glm_freq1)
resid <- ifelse(is.na(resid) | is.infinite(resid), residuals(glm_freq1), resid)
#Residual Q-Q Plot
qqnorm(resid)
qqline(resid, col="blue",lwd=2)
#Plot of Residuals vs. Fitted Values
plot(predict(glm_freq1, train, type = "response"), resid)
abline(h = 0, col = "red")<jupyter_output><empty_output><jupyter_text>There is a clear right tail in the residuals, as seen in the Q-Q plot. This suggests that perhaps we have some outliers in the data. The plot of residuals vs. fitted values seems to show that the model conforms to Assumptions 2 and 3. Overall, this seems to be a decent model.However, Frequency isn't just the number of claims. It is the number of claims *per exposure*. Hence, we need to include the effect of exposure into this model. However, we do not necessarily want to model it. Since exposure is known, we do not need the GLM to estimate a parameter/coefficient for it. We need our model to estimate coefficients for all the rating factors, but then tack the effect of exposure on top of this. Luckily, GLMs allow for this through the **offset** parameter.<jupyter_code>glm_freq2 <- glm(formula = ClaimNb ~ DrivAgeBand + VehAgeBand, data = train, family = poisson(link = "log"), 
                 offset = log(Exposure))
summary(glm_freq2)<jupyter_output><empty_output><jupyter_text>We look at the residual diagnostics of this model.<jupyter_code>#Generate Model Residuals
resid <- qresiduals(glm_freq2)
resid <- ifelse(is.na(resid) | is.infinite(resid), residuals(glm_freq1), resid)
#Residual Q-Q Plot
qqnorm(resid)
qqline(resid, col="blue",lwd=2)
#Plot of Residuals vs. Fitted Values
plot(predict(glm_freq2, train, type = "response"), resid)
abline(h = 0, col = "red")<jupyter_output><empty_output><jupyter_text>Now we get the true picture. We have seen that there is a heavy right tail in the data. The plot of residuals vs. fitted values now shows that the true model actually does not satisfy the constant variance assumption. Hence, we need to improve it. In a bid to do so, we add a third rating factor to our model - Vehicle Brand.<jupyter_code>glm_freq3 <- glm(formula = ClaimNb ~ DrivAgeBand + VehAgeBand + VehBrand, data = train, family = poisson(link = "log"), 
                 offset = log(Exposure))
summary(glm_freq3)<jupyter_output><empty_output><jupyter_text>Vehicle Brand doesn't look as promising as Vehicle Age Band and Driver Age Band. This is because not all levels of Vehicle Brand are *statistically significant*. This means that, for some levels (eg. VehBrandB5 and VehBrand B13), there is a high likelihood our model came up with the coefficient it estimated purely by chance. However since this appears to be the case for only a few levels, we can still accept this rating factor.We look at the residual diagnostics of this model.<jupyter_code>#Generate Model Residuals
resid <- qresiduals(glm_freq3)
resid <- ifelse(is.na(resid) | is.infinite(resid), residuals(glm_freq1), resid)
#Residual Q-Q Plot
qqnorm(resid)
qqline(resid, col="blue",lwd=2)
#Plot of Residuals vs. Fitted Values
plot(predict(glm_freq3, train, type = "response"), resid)
abline(h = 0, col = "red")<jupyter_output><empty_output><jupyter_text>From the residual diagnostics, it doesn't look like this model was much of an improvement. Perhaps adding some more complexity might help.We now look at a model where an interaction effect between Driver Age Band and Vehicle Brand is added.<jupyter_code>glm_freq4 <- glm(formula = ClaimNb ~ DrivAgeBand + VehAgeBand + VehBrand + DrivAgeBand:VehBrand, data = train, family = poisson(link = "log"), 
                 offset = log(Exposure))
summary(glm_freq4)<jupyter_output><empty_output><jupyter_text>This is interesting. Notice that for the interaction between DrivAgeBand(65,Inf] and VehBrand B14 the model returned an `NA` coefficient. This means that this combination is perfectly correlated in our data. Hence, the effect of this combination is already contained in the individual variables itself and having this particular interaction does not add any information to the model, thereby making it redundant. This is not true for other combinations of Driver Age Band and Vehicle Brand. Let's look at the residual diagnostics of this model.<jupyter_code>#Generate Model Residuals
resid <- qresiduals(glm_freq4)
resid <- ifelse(is.na(resid) | is.infinite(resid), residuals(glm_freq1), resid)
#Residual Q-Q Plot
qqnorm(resid)
qqline(resid, col="blue",lwd=2)
#Plot of Residuals vs. Fitted Values
plot(predict(glm_freq4, train, type = "response"), resid)
abline(h = 0, col = "red")<jupyter_output><empty_output><jupyter_text>The above plots tell a similar story for this model as the previous 2 models. Further, we observed a redundancy in our interaction term. Hence, adding this layer of complexity does not seem to help our model all that much.All the diagnostics we have looked at does not tell us one important thing - **how will do the models fare against data that they have not seen before?** The true test of a predictive model lies in how well it can generalize the patterns it has learnt on the data. The only way to evaluate this is to predict the model on completely unseen data and compare it to the actuals. Notice that in the beginning, we imported 2 datasets - *train* and *test*. From the PoV of our models, the *train* dataset is data that the model has fitted on. The *test* dataset is data the model has never seen before. Hence now we predict each of our models on the *test* data and compare these predictions to the actual claim counts.<jupyter_code>test$ClaimNb_Predict1 <- predict(glm_freq1, test, type = "response") #Model 1
test$ClaimNb_Predict2 <- predict(glm_freq2, test, type = "response") #Model 2
test$ClaimNb_Predict3 <- predict(glm_freq3, test, type = "response") #Model 3
test$ClaimNb_Predict4 <- predict(glm_freq4, test, type = "response") #Model 4
#Comparing Predicted Claim Counts to Actual Claim Counts
summary(test$ClaimNb)
summary(test$ClaimNb_Predict1)
summary(test$ClaimNb_Predict2)
summary(test$ClaimNb_Predict3)
summary(test$ClaimNb_Predict4)

sum(test$ClaimNb)
sum(test$ClaimNb_Predict1)
sum(test$ClaimNb_Predict2)
sum(test$ClaimNb_Predict3)
sum(test$ClaimNb_Predict4)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07469797515739322, 'fraction_numerical': 0.007316658158924621, 'mean_word_length': 5.476584022038567, 'pattern_counts': {'":': 0, '<': 91, '<?xml version=': 0, '>': 74, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '292056', 'license': 'MIT', 'n_tokens_mistral': 3484, 'n_tokens_neox': 3266, 'n_words': 1505, 'path': 'Introduction to Predictive Analytics.ipynb', 'repo_name': 'Lux-Actuaries/rvrbinder', 'size': '11754'}