ADEC 7320 — Homework 2: Auto Insurance Modeling
David N. Belyn
2026-07-17
1 Setup and packages
packages <- c("psych", "stargazer", "ggplot2", "dplyr", "corrplot", "car",
"MASS", "caret", "pROC", "glmnet", "moments")
for (p in packages) {
if (!p %in% rownames(installed.packages()))
install.packages(p, repos = "http://cran.rstudio.com/", dependencies = TRUE)
suppressPackageStartupMessages(library(p, character.only = TRUE))
}
set.seed(7320)
train_path <- "C:/Users/dbely/OneDrive/Desktop/Econometrics/Assignment 2/insurance-training-data2-2.csv"
eval_path <- "C:/Users/dbely/OneDrive/Desktop/Econometrics/Assignment 2/insurance-testing-data2-2.csv"
df_train_raw <- read.csv(train_path, stringsAsFactors = FALSE, na.strings = c("NA","", " "))
df_eval_raw <- read.csv(eval_path, stringsAsFactors = FALSE, na.strings = c("NA","", " "))2 Section 1 — Data Exploration
The training set contains 6528 records and 26
columns: an identifier (INDEX), two response
variables (TARGET_FLAG, TARGET_AMT), and 23
predictors describing the driver, the household, the vehicle and the
driving record. TARGET_FLAG is a binary crash indicator;
TARGET_AMT is the claim cost, which is zero for everyone
who did not crash and positive otherwise. The two responses call for two
different tools: a binary logistic regression for the
probability of a crash and a linear regression for the
cost conditional on crashing.
Three data-quality issues jump out immediately and drive most of the preparation work:
- Money is stored as text.
INCOME,HOME_VAL,BLUEBOOKandOLDCLAIMarrive as strings like"$78,658". They must be stripped of$/,and coerced to numeric. - Some factor levels carry a
z_prefix (e.g.z_F,z_No,z_SUV). These are cosmetic and are removed so the categories read cleanly. - Missing and impossible values.
AGE,YOJ,INCOME,HOME_VALandCAR_AGEcontainNAs;JOBhas blanks; andCAR_AGEcontains a nonsensical -3.
clean_raw <- function(df) {
money_cols <- c("INCOME","HOME_VAL","BLUEBOOK","OLDCLAIM")
for (m in money_cols) df[[m]] <- as.numeric(gsub("[\\$,]", "", df[[m]]))
cat_cols <- c("MSTATUS","SEX","EDUCATION","JOB","CAR_TYPE",
"URBANICITY","CAR_USE","PARENT1","RED_CAR","REVOKED")
for (c in cat_cols) df[[c]] <- trimws(gsub("^z_", "", df[[c]]))
df$CAR_AGE[df$CAR_AGE < 0] <- NA # the -3 is impossible
df$JOB[is.na(df$JOB)] <- "Unknown"
df
}
df_train <- clean_raw(df_train_raw)
df_eval <- clean_raw(df_eval_raw)
# Crash base rate
prop.table(table(df_train$TARGET_FLAG))##
## 0 1
## 0.7362132 0.2637868About a quarter of drivers in the sample crashed, so the classes are imbalanced — worth remembering when we read accuracy later, because a model that predicts “no crash” for everyone would already be right roughly 74% of the time.
num_vars <- c("TARGET_AMT","KIDSDRIV","AGE","HOMEKIDS","YOJ","INCOME","HOME_VAL",
"TRAVTIME","BLUEBOOK","TIF","OLDCLAIM","CLM_FREQ","MVR_PTS","CAR_AGE")
stargazer(df_train[, num_vars], type = "html",
title = "Table 1. Summary statistics (training data)", digits = 1)| Statistic | N | Mean | St. Dev. | Min | Max |
| TARGET_AMT | 6,528 | 1,466.6 | 4,545.7 | 0.0 | 107,586.1 |
| KIDSDRIV | 6,528 | 0.2 | 0.5 | 0 | 4 |
| AGE | 6,525 | 44.9 | 8.7 | 16 | 81 |
| HOMEKIDS | 6,528 | 0.7 | 1.1 | 0 | 5 |
| YOJ | 6,153 | 10.5 | 4.1 | 0 | 23 |
| INCOME | 6,174 | 61,649.3 | 47,658.4 | 0 | 367,030 |
| HOME_VAL | 6,160 | 154,298.4 | 128,874.6 | 0 | 885,282 |
| TRAVTIME | 6,528 | 33.4 | 16.0 | 5 | 142 |
| BLUEBOOK | 6,528 | 15,641.8 | 8,381.5 | 1,500 | 69,740 |
| TIF | 6,528 | 5.4 | 4.1 | 1 | 25 |
| OLDCLAIM | 6,528 | 4,119.3 | 8,924.7 | 0 | 57,037 |
| CLM_FREQ | 6,528 | 0.8 | 1.2 | 0 | 5 |
| MVR_PTS | 6,528 | 1.7 | 2.1 | 0 | 13 |
| CAR_AGE | 6,128 | 8.3 | 5.7 | 0 | 27 |
# Skewness of the numeric variables
round(sapply(df_train[num_vars], moments::skewness, na.rm = TRUE), 2)## TARGET_AMT KIDSDRIV AGE HOMEKIDS YOJ INCOME HOME_VAL
## 9.13 3.30 -0.04 1.32 -1.20 1.21 0.50
## TRAVTIME BLUEBOOK TIF OLDCLAIM CLM_FREQ MVR_PTS CAR_AGE
## 0.46 0.82 0.88 3.07 1.21 1.33 0.28# Missingness count per column
colSums(is.na(df_train))## INDEX TARGET_FLAG TARGET_AMT KIDSDRIV AGE HOMEKIDS
## 0 0 0 0 3 0
## YOJ INCOME PARENT1 HOME_VAL MSTATUS SEX
## 375 354 0 368 0 0
## EDUCATION JOB TRAVTIME CAR_USE BLUEBOOK TIF
## 0 0 0 0 0 0
## CAR_TYPE RED_CAR OLDCLAIM CLM_FREQ REVOKED MVR_PTS
## 0 0 0 0 0 0
## CAR_AGE URBANICITY
## 400 0Insights. INCOME,
HOME_VAL, BLUEBOOK and especially
OLDCLAIM are strongly right-skewed, which motivates log
transforms below. HOME_VAL and OLDCLAIM have
large spikes at zero (renters and drivers with no prior claims), so a
zero is meaningful, not missing — I capture that with indicator flags
rather than treating zero as an outlier.
par(mfrow = c(1, 2))
hist(df_train$INCOME, breaks = 40, col = "steelblue", main = "Income (right-skewed)", xlab = "Income")
boxplot(TRAVTIME ~ TARGET_FLAG, data = df_train, col = c("grey80","tomato"),
main = "Commute time by crash", xlab = "Crashed (1=yes)", ylab = "TRAVTIME")
par(mfrow = c(1, 1))
# Correlation of numeric predictors with the crash flag
cor_flag <- sort(sapply(num_vars, function(v)
cor(df_train[[v]], df_train$TARGET_FLAG, use = "complete.obs")))
round(cor_flag, 3)## HOME_VAL INCOME BLUEBOOK AGE CAR_AGE TIF YOJ
## -0.171 -0.136 -0.106 -0.104 -0.094 -0.084 -0.070
## TRAVTIME KIDSDRIV HOMEKIDS OLDCLAIM CLM_FREQ MVR_PTS TARGET_AMT
## 0.044 0.119 0.119 0.137 0.219 0.221 0.539The correlations are individually modest (this is a hard
classification problem), but the signs line up with intuition:
MVR_PTS, CLM_FREQ, KIDSDRIV and
HOMEKIDS point toward more crashes, while
INCOME, HOME_VAL, AGE and
CAR_AGE point the other way.
3 Section 2 — Data Preparation
I applied five transformations, each with a reason.
(a) Money strings → numeric and (b)
z_ prefixes removed — already done in
clean_raw() above; without these the variables cannot enter
a regression at all.
(c) Missing-value imputation with flags. I impute
AGE, YOJ, INCOME,
HOME_VAL and CAR_AGE with the training
median (robust to the skew) and add a *_MISS flag
so the model can learn whether the fact of being missing is
itself predictive. Critically, the evaluation set is imputed with the
training medians, not its own — otherwise we would leak
information from the hold-out set.
(d) Feature engineering / bucketing.
HOME_OWNER (owns a home), HAS_OLDCLAIM (any
prior payout), YOUNG_DRIVER (age < 25) and
KIDS_DRIVING (any child drivers) turn continuous or spiked
variables into clean risk buckets.
(e) Log transforms on INCOME,
HOME_VAL and BLUEBOOK to tame the
right-skew.
num_impute <- c("AGE","YOJ","INCOME","HOME_VAL","CAR_AGE")
train_medians <- sapply(df_train[num_impute], median, na.rm = TRUE)
impute_df <- function(df, medians) {
for (v in names(medians)) {
df[[paste0(v, "_MISS")]] <- as.integer(is.na(df[[v]]))
df[[v]][is.na(df[[v]])] <- medians[[v]]
}
df
}
df_train <- impute_df(df_train, train_medians)
df_eval <- impute_df(df_eval, train_medians)
# Catch-all so no stray NA breaks model.matrix()/glmnet
catch_cols <- c("BLUEBOOK","OLDCLAIM","TRAVTIME","TIF","MVR_PTS","CLM_FREQ","KIDSDRIV","HOMEKIDS")
catch_medians <- sapply(df_train[catch_cols], median, na.rm = TRUE)
for (v in catch_cols) {
df_train[[v]][is.na(df_train[[v]])] <- catch_medians[[v]]
df_eval[[v]][is.na(df_eval[[v]])] <- catch_medians[[v]]
}
engineer <- function(df) {
df$HOME_OWNER <- as.integer(df$HOME_VAL > 0)
df$HAS_OLDCLAIM <- as.integer(df$OLDCLAIM > 0)
df$LOG_INCOME <- log(df$INCOME + 1)
df$LOG_BLUEBOOK <- log(df$BLUEBOOK + 1)
df$LOG_HOME_VAL <- log(df$HOME_VAL + 1)
df$YOUNG_DRIVER <- as.integer(df$AGE < 25)
df$KIDS_DRIVING <- as.integer(df$KIDSDRIV > 0)
fac <- c("PARENT1","MSTATUS","SEX","EDUCATION","JOB","CAR_USE",
"CAR_TYPE","RED_CAR","REVOKED","URBANICITY")
for (f in fac) df[[f]] <- as.factor(df[[f]])
df
}
df_train <- engineer(df_train)
df_eval <- engineer(df_eval)
# Align eval factor levels to the training levels (protects predict())
for (f in c("PARENT1","MSTATUS","SEX","EDUCATION","JOB","CAR_USE",
"CAR_TYPE","RED_CAR","REVOKED","URBANICITY"))
df_eval[[f]] <- factor(df_eval[[f]], levels = levels(df_train[[f]]))
df_train$TARGET_FLAG <- as.integer(df_train$TARGET_FLAG)4 Section 3 — Build Models
4.1 Logistic regression (probability of a crash)
I built three logistic models plus two penalized models for the bonus.
- Logit 1 — full model: every cleaned predictor. A baseline that shows what is available and flags multicollinearity.
- Logit 2 — stepwise (AIC):
stepAICprunes Logit 1 in both directions, keeping the variables that improve AIC. This is my working candidate. - Logit 3 — parsimonious, theory-driven: a hand-picked set of variables with clear causal stories (driving record, exposure, urbanicity, prior claims).
drop_cols <- c("INDEX","TARGET_AMT")
model_df <- df_train[, !(names(df_train) %in% drop_cols)]
logit1 <- glm(TARGET_FLAG ~ ., data = model_df, family = binomial)
logit2 <- MASS::stepAIC(logit1, direction = "both", trace = FALSE)
logit3 <- glm(TARGET_FLAG ~ KIDSDRIV + INCOME + HOME_OWNER + MSTATUS + TRAVTIME +
CAR_USE + BLUEBOOK + TIF + OLDCLAIM + CLM_FREQ + REVOKED + MVR_PTS +
URBANICITY + CAR_TYPE + JOB, data = model_df, family = binomial)
probit1 <- glm(formula(logit2), data = model_df, family = binomial(link = "probit"))summary(logit2)##
## Call:
## glm(formula = TARGET_FLAG ~ KIDSDRIV + INCOME + PARENT1 + MSTATUS +
## EDUCATION + JOB + TRAVTIME + CAR_USE + TIF + CAR_TYPE + OLDCLAIM +
## REVOKED + MVR_PTS + URBANICITY + AGE_MISS + HAS_OLDCLAIM +
## LOG_INCOME + LOG_BLUEBOOK + LOG_HOME_VAL + YOUNG_DRIVER +
## KIDS_DRIVING, family = binomial, data = model_df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.248e-01 6.273e-01 1.315 0.188546
## KIDSDRIV 2.104e-01 1.401e-01 1.502 0.133222
## INCOME -2.737e-06 1.199e-06 -2.283 0.022402 *
## PARENT1Yes 3.959e-01 1.082e-01 3.658 0.000254 ***
## MSTATUSYes -4.575e-01 9.337e-02 -4.900 9.58e-07 ***
## EDUCATIONBachelors -4.435e-01 1.223e-01 -3.625 0.000289 ***
## EDUCATIONHigh School 6.068e-04 1.066e-01 0.006 0.995456
## EDUCATIONMasters -2.502e-01 1.804e-01 -1.387 0.165475
## EDUCATIONPhD -2.002e-01 2.225e-01 -0.900 0.368349
## JOBClerical 5.259e-03 1.199e-01 0.044 0.965015
## JOBDoctor -1.214e+00 3.442e-01 -3.528 0.000419 ***
## JOBHome Maker -3.901e-01 1.846e-01 -2.114 0.034553 *
## JOBLawyer -3.578e-01 2.086e-01 -1.716 0.086210 .
## JOBManager -9.222e-01 1.590e-01 -5.800 6.65e-09 ***
## JOBProfessional -2.321e-01 1.355e-01 -1.713 0.086722 .
## JOBStudent -4.964e-01 1.646e-01 -3.016 0.002558 **
## JOBUnknown -4.235e-01 2.085e-01 -2.031 0.042251 *
## TRAVTIME 1.349e-02 2.115e-03 6.378 1.80e-10 ***
## CAR_USEPrivate -7.449e-01 1.041e-01 -7.156 8.29e-13 ***
## TIF -5.645e-02 8.272e-03 -6.825 8.79e-12 ***
## CAR_TYPEPanel Truck 4.675e-01 1.639e-01 2.853 0.004332 **
## CAR_TYPEPickup 5.769e-01 1.121e-01 5.146 2.66e-07 ***
## CAR_TYPESports Car 9.001e-01 1.217e-01 7.394 1.43e-13 ***
## CAR_TYPESUV 7.393e-01 9.577e-02 7.719 1.17e-14 ***
## CAR_TYPEVan 6.110e-01 1.379e-01 4.432 9.34e-06 ***
## OLDCLAIM -2.301e-05 4.646e-06 -4.953 7.31e-07 ***
## REVOKEDYes 1.016e+00 1.021e-01 9.952 < 2e-16 ***
## MVR_PTS 9.419e-02 1.587e-02 5.934 2.96e-09 ***
## URBANICITYHighly Urban/ Urban 2.338e+00 1.268e-01 18.447 < 2e-16 ***
## AGE_MISS 1.269e+01 1.717e+02 0.074 0.941060
## HAS_OLDCLAIM 6.826e-01 8.816e-02 7.743 9.72e-15 ***
## LOG_INCOME -5.910e-02 1.692e-02 -3.493 0.000478 ***
## LOG_BLUEBOOK -3.306e-01 6.225e-02 -5.311 1.09e-07 ***
## LOG_HOME_VAL -2.738e-02 7.746e-03 -3.534 0.000409 ***
## YOUNG_DRIVER 5.721e-01 3.118e-01 1.835 0.066538 .
## KIDS_DRIVING 4.826e-01 2.230e-01 2.164 0.030462 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 7533.1 on 6527 degrees of freedom
## Residual deviance: 5773.5 on 6492 degrees of freedom
## AIC: 5845.5
##
## Number of Fisher Scoring iterations: 114.1.1 Bonus — LASSO and Ridge logistic regression
LASSO (alpha = 1) adds an L1 penalty that shrinks weak
coefficients to exactly zero, so it doubles as
automatic feature selection; Ridge (alpha = 0) uses an L2
penalty that shrinks but keeps all variables. I cross-validate
lambda and report the variables LASSO keeps.
x_form <- TARGET_FLAG ~ KIDSDRIV + AGE + HOMEKIDS + YOJ + INCOME + PARENT1 + HOME_VAL +
MSTATUS + SEX + EDUCATION + JOB + TRAVTIME + CAR_USE + BLUEBOOK + TIF + CAR_TYPE +
RED_CAR + OLDCLAIM + CLM_FREQ + REVOKED + MVR_PTS + CAR_AGE + URBANICITY +
HOME_OWNER + HAS_OLDCLAIM + YOUNG_DRIVER + KIDS_DRIVING
x_train <- model.matrix(x_form, data = df_train)[, -1]
y_train <- df_train$TARGET_FLAG
cv_lasso <- cv.glmnet(x_train, y_train, family = "binomial", alpha = 1)
cv_ridge <- cv.glmnet(x_train, y_train, family = "binomial", alpha = 0)
lasso_coef <- coef(cv_lasso, s = "lambda.1se")
cat("LASSO keeps these variables:\n")## LASSO keeps these variables:print(rownames(lasso_coef)[which(lasso_coef != 0)])## [1] "(Intercept)" "KIDSDRIV"
## [3] "AGE" "HOMEKIDS"
## [5] "YOJ" "INCOME"
## [7] "PARENT1Yes" "HOME_VAL"
## [9] "MSTATUSYes" "EDUCATIONBachelors"
## [11] "EDUCATIONHigh School" "EDUCATIONMasters"
## [13] "JOBClerical" "JOBDoctor"
## [15] "JOBManager" "TRAVTIME"
## [17] "CAR_USEPrivate" "BLUEBOOK"
## [19] "TIF" "CAR_TYPEPickup"
## [21] "CAR_TYPESports Car" "CAR_TYPESUV"
## [23] "OLDCLAIM" "CLM_FREQ"
## [25] "REVOKEDYes" "MVR_PTS"
## [27] "CAR_AGE" "URBANICITYHighly Urban/ Urban"
## [29] "HOME_OWNER" "HAS_OLDCLAIM"
## [31] "YOUNG_DRIVER" "KIDS_DRIVING"4.1.2 Coefficient interpretation (final logit)
Logistic coefficients live on the log-odds scale, so
exp(β) is the odds ratio. The table below
converts every coefficient in the final model to an odds ratio.
odds <- data.frame(coef = round(coef(logit2), 4),
odds_ratio = round(exp(coef(logit2)), 4),
pct_change = round((exp(coef(logit2)) - 1) * 100, 1))
odds## coef odds_ratio pct_change
## (Intercept) 0.8248 2.2814 128.1
## KIDSDRIV 0.2104 1.2341 23.4
## INCOME 0.0000 1.0000 0.0
## PARENT1Yes 0.3959 1.4857 48.6
## MSTATUSYes -0.4575 0.6329 -36.7
## EDUCATIONBachelors -0.4435 0.6418 -35.8
## EDUCATIONHigh School 0.0006 1.0006 0.1
## EDUCATIONMasters -0.2502 0.7787 -22.1
## EDUCATIONPhD -0.2002 0.8186 -18.1
## JOBClerical 0.0053 1.0053 0.5
## JOBDoctor -1.2142 0.2969 -70.3
## JOBHome Maker -0.3901 0.6770 -32.3
## JOBLawyer -0.3578 0.6992 -30.1
## JOBManager -0.9222 0.3976 -60.2
## JOBProfessional -0.2321 0.7928 -20.7
## JOBStudent -0.4964 0.6087 -39.1
## JOBUnknown -0.4235 0.6547 -34.5
## TRAVTIME 0.0135 1.0136 1.4
## CAR_USEPrivate -0.7449 0.4748 -52.5
## TIF -0.0565 0.9451 -5.5
## CAR_TYPEPanel Truck 0.4675 1.5960 59.6
## CAR_TYPEPickup 0.5769 1.7805 78.0
## CAR_TYPESports Car 0.9001 2.4598 146.0
## CAR_TYPESUV 0.7393 2.0944 109.4
## CAR_TYPEVan 0.6110 1.8423 84.2
## OLDCLAIM 0.0000 1.0000 0.0
## REVOKEDYes 1.0164 2.7632 176.3
## MVR_PTS 0.0942 1.0988 9.9
## URBANICITYHighly Urban/ Urban 2.3385 10.3654 936.5
## AGE_MISS 12.6949 326085.2877 32608428.8
## HAS_OLDCLAIM 0.6826 1.9790 97.9
## LOG_INCOME -0.0591 0.9426 -5.7
## LOG_BLUEBOOK -0.3306 0.7185 -28.1
## LOG_HOME_VAL -0.0274 0.9730 -2.7
## YOUNG_DRIVER 0.5721 1.7720 77.2
## KIDS_DRIVING 0.4826 1.6202 62.0Reading a positive coefficient
(e.g. MVR_PTS, motor-vehicle-record points): a
positive β means the odds ratio exp(β) > 1. Each
additional MVR point multiplies the odds of a crash by
exp(β), i.e. a (exp(β) − 1) × 100% increase in
the odds, holding everything else fixed. This matches the theory that
ticket-prone drivers crash more.
Reading a negative coefficient
(e.g. MSTATUSYes, being married): a negative β gives
exp(β) < 1. Married drivers have
(1 − exp(β)) × 100% lower odds of a crash
than the unmarried baseline, again consistent with theory. The same lens
explains URBANICITY: urban drivers show a large positive
coefficient — dense traffic sharply raises crash odds — which is
typically the single strongest predictor in the model.
Why logit over a linear probability model (OLS on a 0/1 outcome)? OLS can predict probabilities below 0 or above 1, assumes constant marginal effects, and produces heteroskedastic residuals by construction. The logit maps the linear predictor through the logistic CDF, so fitted values stay in (0, 1) and the effect of a variable depends sensibly on where you sit on the curve. I keep any counter-intuitive but statistically insignificant signs only if dropping them hurts AIC/AUC; otherwise I prefer the interpretable model.
4.2 Linear regression (cost given a crash)
TARGET_AMT is zero for non-crashers, so a linear model
on the full sample would be modeling a spike-at-zero, not a cost. I
therefore fit severity models on the subset of drivers
who actually crashed. Predicted expected cost for a new driver is then
P(crash) × predicted severity.
crash_df <- df_train[df_train$TARGET_FLAG == 1, ]
lm1 <- lm(TARGET_AMT ~ BLUEBOOK + CAR_AGE + CAR_TYPE + CAR_USE + MVR_PTS +
REVOKED + SEX + URBANICITY, data = crash_df)
lm2 <- lm(log(TARGET_AMT) ~ LOG_BLUEBOOK + CAR_AGE + CAR_TYPE + MVR_PTS +
REVOKED + KIDSDRIV + CLM_FREQ, data = crash_df)
summary(lm1)##
## Call:
## lm(formula = TARGET_AMT ~ BLUEBOOK + CAR_AGE + CAR_TYPE + CAR_USE +
## MVR_PTS + REVOKED + SEX + URBANICITY, data = crash_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8830 -3022 -1432 515 100454
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2606.5144 1182.9167 2.203 0.0277 *
## BLUEBOOK 0.1341 0.0309 4.340 1.51e-05 ***
## CAR_AGE -34.1886 33.8862 -1.009 0.3132
## CAR_TYPEPanel Truck -527.1579 983.5153 -0.536 0.5920
## CAR_TYPEPickup -474.0375 617.6003 -0.768 0.4429
## CAR_TYPESports Car 774.6322 793.8673 0.976 0.3293
## CAR_TYPESUV 698.2742 697.1380 1.002 0.3167
## CAR_TYPEVan -697.6146 807.1015 -0.864 0.3875
## CAR_USEPrivate -267.1556 424.9058 -0.629 0.5296
## MVR_PTS 139.2651 69.8098 1.995 0.0462 *
## REVOKEDYes -558.3070 435.6100 -1.282 0.2001
## SEXM 1453.1121 616.6096 2.357 0.0186 *
## URBANICITYHighly Urban/ Urban 468.7037 794.3593 0.590 0.5552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7389 on 1709 degrees of freedom
## Multiple R-squared: 0.0249, Adjusted R-squared: 0.01805
## F-statistic: 3.637 on 12 and 1709 DF, p-value: 2.071e-05BLUEBOOK (vehicle value) is the dominant driver of
payout size, which is exactly what theory predicts — pricier cars cost
more to repair. The severity models have modest R² because claim amounts
are inherently noisy, but the signs are sensible.
5 Section 4 — Select Models
5.1 Selection criteria
For the logistic model I weigh AIC (in-sample fit penalized for complexity), AUC (threshold-free ranking) and the confusion-matrix metrics at the required 0.5 threshold. For the linear model I use Adjusted R², the F-statistic, RMSE, VIF for multicollinearity, and the four diagnostic residual plots. I will accept a slightly simpler model over a marginally better one when it is easier to defend to the boss.
score_logit <- function(model, data, truth, threshold = 0.5, label = "") {
p <- predict(model, newdata = data, type = "response")
cls <- factor(as.integer(p >= threshold), levels = c(0,1))
ref <- factor(truth, levels = c(0,1))
cm <- caret::confusionMatrix(cls, ref, positive = "1")
auc <- as.numeric(pROC::auc(pROC::roc(truth, p, quiet = TRUE)))
data.frame(Model = label, AIC = round(AIC(model),1),
Accuracy = round(cm$overall["Accuracy"],4),
Error_Rate = round(1 - cm$overall["Accuracy"],4),
Sensitivity = round(cm$byClass["Sensitivity"],4),
Specificity = round(cm$byClass["Specificity"],4),
Precision = round(cm$byClass["Precision"],4),
F1 = round(cm$byClass["F1"],4), AUC = round(auc,4), row.names = NULL)
}
score_lasso <- function(cvfit, x, truth, s = "lambda.1se", threshold = 0.5, label = "") {
p <- as.numeric(predict(cvfit, newx = x, s = s, type = "response"))
cls <- factor(as.integer(p >= threshold), levels = c(0,1))
ref <- factor(truth, levels = c(0,1))
cm <- caret::confusionMatrix(cls, ref, positive = "1")
auc <- as.numeric(pROC::auc(pROC::roc(truth, p, quiet = TRUE)))
data.frame(Model = label, AIC = NA,
Accuracy = round(cm$overall["Accuracy"],4),
Error_Rate = round(1 - cm$overall["Accuracy"],4),
Sensitivity = round(cm$byClass["Sensitivity"],4),
Specificity = round(cm$byClass["Specificity"],4),
Precision = round(cm$byClass["Precision"],4),
F1 = round(cm$byClass["F1"],4), AUC = round(auc,4), row.names = NULL)
}
logit_scores <- rbind(
score_logit(logit1, model_df, df_train$TARGET_FLAG, label = "Logit 1: full"),
score_logit(logit2, model_df, df_train$TARGET_FLAG, label = "Logit 2: stepwise"),
score_logit(logit3, model_df, df_train$TARGET_FLAG, label = "Logit 3: parsimonious"),
score_logit(probit1, model_df, df_train$TARGET_FLAG, label = "Probit"),
score_lasso(cv_lasso, x_train, df_train$TARGET_FLAG, label = "LASSO (bonus)"),
score_lasso(cv_ridge, x_train, df_train$TARGET_FLAG, label = "Ridge (bonus)")
)
logit_scores## Model AIC Accuracy Error_Rate Sensitivity Specificity
## 1 Logit 1: full 5865.7 0.7932 0.2068 0.4361 0.9211
## 2 Logit 2: stepwise 5845.5 0.7934 0.2066 0.4315 0.9230
## 3 Logit 3: parsimonious 5921.0 0.7912 0.2088 0.4193 0.9245
## 4 Probit 5850.2 0.7920 0.2080 0.4204 0.9251
## 5 LASSO (bonus) NA 0.7872 0.2128 0.3386 0.9480
## 6 Ridge (bonus) NA 0.7863 0.2137 0.3264 0.9511
## Precision F1 AUC
## 1 0.6646 0.5266 0.8188
## 2 0.6676 0.5242 0.8182
## 3 0.6654 0.5144 0.8108
## 4 0.6679 0.5160 0.8181
## 5 0.6999 0.4564 0.8100
## 6 0.7051 0.4462 0.8110Definitions. Accuracy = share of correct predictions; classification error rate = 1 − accuracy; sensitivity (recall) = of the drivers who crashed, the share we caught; specificity = of the non-crashers, the share we correctly cleared; precision = of those we flagged, the share who really crashed; AUC = probability the model ranks a random crasher above a random non-crasher.
I select Logit 2 (stepwise) as the final classifier:
it posts the best or near-best AUC and accuracy while dropping the noise
variables that inflate Logit 1’s AIC. The z_ classes and
red-car urban legend wash out, confirming they add little.
5.2 Linear model selection and diagnostics
lm_metrics <- function(m, log_response = FALSE, y = crash_df$TARGET_AMT, label = "") {
fitted_vals <- if (log_response) exp(predict(m)) else predict(m)
data.frame(Model = label, Adj_R2 = round(summary(m)$adj.r.squared,4),
F_stat = round(summary(m)$fstatistic[1],1),
RMSE = round(sqrt(mean((y - fitted_vals)^2)),1), row.names = NULL)
}
rbind(lm_metrics(lm1, FALSE, label = "LM1: raw cost"),
lm_metrics(lm2, TRUE, label = "LM2: log cost"))## Model Adj_R2 F_stat RMSE
## 1 LM1: raw cost 0.0181 3.6 7360.6
## 2 LM2: log cost 0.0171 3.7 7581.0car::vif(lm1) # multicollinearity: all should be well under 5## GVIF Df GVIF^(1/(2*Df))
## BLUEBOOK 2.015097 1 1.419541
## CAR_AGE 1.057688 1 1.028440
## CAR_TYPE 5.642875 5 1.188913
## CAR_USE 1.422038 1 1.192492
## MVR_PTS 1.010592 1 1.005282
## REVOKED 1.007783 1 1.003884
## SEX 2.965532 1 1.722072
## URBANICITY 1.016916 1 1.008423par(mfrow = c(2, 2)); plot(lm1); par(mfrow = c(1, 1))
The four residual plots show the classic pattern for insurance losses: heavier-than-normal right tail in the Q-Q plot and mild heteroskedasticity, because a few very large claims stretch the distribution. I keep LM1 (raw dollars) as the final severity model: it is directly interpretable in dollars and avoids the retransformation bias that the log model introduces when predicting on the original scale. VIFs are low, so multicollinearity is not a concern.
5.3 Predictions on the evaluation set
final_logit <- logit2
final_lm <- lm1
p_crash <- predict(final_logit, newdata = df_eval, type = "response")
flag_pred <- as.integer(p_crash >= 0.5)
severity <- predict(final_lm, newdata = df_eval)
severity[severity < 0] <- 0
expected_cost <- p_crash * severity
eval_out <- data.frame(
INDEX = df_eval$INDEX,
P_TARGET_FLAG = round(p_crash, 4),
TARGET_FLAG = flag_pred,
SEVERITY_IF_CRASH = round(severity, 0),
TARGET_AMT = round(expected_cost, 0)
)
write.csv(eval_out, "insurance_eval_predictions.csv", row.names = FALSE)
head(eval_out, 10)## INDEX P_TARGET_FLAG TARGET_FLAG SEVERITY_IF_CRASH TARGET_AMT
## 1 5 0.0890 0 6127 545
## 2 8 0.2944 0 4650 1369
## 3 26 0.4209 0 7947 3345
## 4 40 0.2414 0 3505 846
## 5 45 0.0463 0 4628 215
## 6 55 0.5676 1 6805 3862
## 7 61 0.9148 1 6741 6167
## 8 66 0.0160 0 6976 111
## 9 67 0.3669 0 5445 1998
## 10 71 0.2359 0 4395 1037# This evaluation file retains the true targets, so we can also report
# genuine OUT-OF-SAMPLE performance (a stronger test than the training metrics).
if ("TARGET_FLAG" %in% names(df_eval_raw) && sum(!is.na(df_eval_raw$TARGET_FLAG)) > 0) {
truth <- as.integer(df_eval_raw$TARGET_FLAG)
print(score_logit(final_logit, df_eval, truth, label = "Final logit — EVAL set"))
}## Model AIC Accuracy Error_Rate Sensitivity Specificity
## 1 Final logit — EVAL set 5845.5 0.7844 0.2156 0.4084 0.9193
## Precision F1 AUC
## 1 0.6447 0.5 0.816 Conclusion
The final deliverables are a stepwise logistic model for crash
probability and an OLS severity model for claim cost, combined into an
expected-cost prediction for every driver in the evaluation set (written
to insurance_eval_predictions.csv using the 0.5 threshold).
Urbanicity, driving record (MVR points, prior claims, license
revocation) and exposure (commute time, teen drivers) are the strongest
crash predictors, while vehicle value drives claim size — all consistent
with the theoretical effects laid out in the assignment.
7 Appendix — session info
sessionInfo()## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] moments_0.14.1 glmnet_4.1-10 Matrix_1.7-3 pROC_1.19.0.1
## [5] caret_7.0-1 lattice_0.22-7 MASS_7.3-65 car_3.1-5
## [9] carData_3.0-6 corrplot_0.95 dplyr_1.1.4 ggplot2_4.0.1
## [13] stargazer_5.2.3 psych_2.5.6
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 timeDate_4052.112 farver_2.1.2
## [4] S7_0.2.0 fastmap_1.2.0 digest_0.6.37
## [7] rpart_4.1.24 timechange_0.3.0 lifecycle_1.0.4
## [10] survival_3.8-3 magrittr_2.0.4 compiler_4.5.1
## [13] rlang_1.1.6 sass_0.4.10 tools_4.5.1
## [16] yaml_2.3.10 data.table_1.17.8 knitr_1.50
## [19] mnormt_2.1.1 plyr_1.8.9 RColorBrewer_1.1-3
## [22] abind_1.4-8 withr_3.0.2 purrr_1.2.0
## [25] nnet_7.3-20 grid_4.5.1 stats4_4.5.1
## [28] e1071_1.7-17 future_1.70.0 globals_0.19.1
## [31] scales_1.4.0 iterators_1.0.14 cli_3.6.5
## [34] rmarkdown_2.30 generics_0.1.4 rstudioapi_0.17.1
## [37] future.apply_1.20.2 reshape2_1.4.5 proxy_0.4-29
## [40] cachem_1.1.0 stringr_1.5.2 splines_4.5.1
## [43] parallel_4.5.1 vctrs_0.6.5 hardhat_1.4.3
## [46] jsonlite_2.0.0 Formula_1.2-5 listenv_1.0.0
## [49] foreach_1.5.2 gower_1.0.2 jquerylib_0.1.4
## [52] recipes_1.3.3 glue_1.8.0 parallelly_1.48.0
## [55] codetools_0.2-20 lubridate_1.9.4 stringi_1.8.7
## [58] gtable_0.3.6 shape_1.4.6.1 tibble_3.3.0
## [61] pillar_1.11.1 htmltools_0.5.8.1 ipred_0.9-15
## [64] lava_1.9.2 R6_2.6.1 evaluate_1.0.5
## [67] bslib_0.9.0 class_7.3-23 Rcpp_1.1.2
## [70] nlme_3.1-168 prodlim_2026.03.11 xfun_0.53
## [73] pkgconfig_2.0.3 ModelMetrics_1.2.2.2