Introduction

An auto insurer wants two things out of this book of business: (1) the probability that a policyholder crashes their car (TARGET_FLAG), and (2) if they do crash, how much that claim costs (TARGET_AMT). I build binary logistic regression models for the first question and multiple linear regression models for the second, using the ~8,000-record training panel and a held-out evaluation set of 1,633 policies.

I want to flag one thing up front: TARGET_AMT is zero for every policy that didn’t crash and strictly positive (and heavily right-skewed) for the ones that did. That single fact drives a lot of the modeling choices below, so I come back to it repeatedly rather than pretending it isn’t there.

data_dir   <- "C:/Users/dbely/OneDrive/Desktop/Econometrics/Assignment 2"
train_path <- file.path(data_dir, "insurance-training-data2-2.csv")
test_path  <- file.path(data_dir, "insurance-testing-data2-2.csv")

data_files <- c(train_path, test_path)
missing_files <- data_files[!file.exists(data_files)]
if (length(missing_files) > 0) {
  stop("Missing data file(s): ", paste(missing_files, collapse = ", "),
       ". Double-check these paths and filenames and re-knit.")
}

train_raw <- read.csv(train_path, stringsAsFactors = FALSE)
test_raw  <- read.csv(test_path,  stringsAsFactors = FALSE)

dim(train_raw)
## [1] 6528   26
dim(test_raw)
## [1] 1633   26

The training file has 6528 records and 26 columns (INDEX plus the two response variables plus 23 candidate predictors). The evaluation file has 1633 records in the same layout.

Data Exploration

Structure and a first look

str(train_raw)
## 'data.frame':    6528 obs. of  26 variables:
## $ INDEX : int 8023 4658 2005 4396 3706 3871 8670 833 5896 7432 ...
## $ TARGET_FLAG: int 0 0 0 0 0 0 0 1 0 0 ...
## $ TARGET_AMT : num 0 0 0 0 0 ...
## $ KIDSDRIV : int 0 0 0 0 0 2 0 0 0 0 ...
## $ AGE : int 40 48 51 56 51 45 65 56 53 49 ...
## $ HOMEKIDS : int 0 0 0 0 0 2 0 0 1 0 ...
## $ YOJ : int 9 0 14 10 9 11 15 13 10 12 ...
## $ INCOME : chr "$78,658" "$0" "$89,606" "$199,336" ...
## $ PARENT1 : chr "No" "No" "No" "No" ...
## $ HOME_VAL : chr "$236,005" "$75,688" "$288,266" "$501,076" ...
## $ MSTATUS : chr "z_No" "Yes" "z_No" "Yes" ...
## $ SEX : chr "z_F" "z_F" "z_F" "z_F" ...
## $ EDUCATION : chr "Masters" "Masters" "Bachelors" "PhD" ...
## $ JOB : chr "Lawyer" "Home Maker" "Clerical" "" ...
## $ TRAVTIME : int 69 39 32 23 46 14 39 48 51 26 ...
## $ CAR_USE : chr "Private" "Private" "Private" "Commercial" ...
## $ BLUEBOOK : chr "$17,040" "$1,500" "$10,460" "$41,280" ...
## $ TIF : int 4 6 1 1 6 3 4 1 10 1 ...
## $ CAR_TYPE : chr "Pickup" "Sports Car" "z_SUV" "Van" ...
## $ RED_CAR : chr "no" "no" "no" "no" ...
## $ OLDCLAIM : chr "$0" "$3,867" "$21,581" "$0" ...
## $ CLM_FREQ : int 0 3 2 0 0 2 0 1 0 0 ...
## $ REVOKED : chr "No" "No" "Yes" "No" ...
## $ MVR_PTS : int 1 5 0 1 0 3 1 7 0 0 ...
## $ CAR_AGE : int 15 NA 8 16 14 7 11 15 NA 10 ...
## $ URBANICITY : chr "z_Highly Rural/ Rural" "Highly Urban/ Urban"
##    "Highly Urban/ Urban" "Highly Urban/ Urban" ...

Several columns that are obviously numeric (INCOME, HOME_VAL, BLUEBOOK, OLDCLAIM) were read in as text because the raw values are formatted like "$78,658". I deal with that in Data Preparation. For now I look at the columns in their raw form to document what’s actually in the file.

Missing data

miss_tbl <- sapply(train_raw, function(x) sum(is.na(x) | x == ""))
miss_tbl <- miss_tbl[miss_tbl > 0]
kable(data.frame(Variable = names(miss_tbl), Missing = miss_tbl,
                  Pct = round(100 * miss_tbl / nrow(train_raw), 1)),
      row.names = FALSE, caption = "Missing values in the training data")
Missing values in the training data
Variable Missing Pct
AGE 3 0.0
YOJ 375 5.7
INCOME 354 5.4
HOME_VAL 368 5.6
JOB 419 6.4
CAR_AGE 399 6.1

Six variables have missing values: AGE (barely, 3 records), and then YOJ, INCOME, HOME_VAL, JOB, and CAR_AGE each missing somewhere between 5% and 6% of records. Nothing is missing more than about 6% of the time, which is a relief, since I don’t have to make any heroic imputation assumptions. CAR_AGE also has at least one clearly impossible value (a negative vehicle age), which I treat as a data-entry error and recode to missing before imputing it.

Summary statistics

parse_money <- function(x) as.numeric(gsub("[$,]", "", x))

num_preview <- train_raw %>%
  transmute(TARGET_FLAG, TARGET_AMT, AGE, YOJ,
            INCOME   = parse_money(INCOME),
            HOME_VAL = parse_money(HOME_VAL),
            BLUEBOOK = parse_money(BLUEBOOK),
            OLDCLAIM = parse_money(OLDCLAIM),
            TRAVTIME, TIF, CLM_FREQ, MVR_PTS, CAR_AGE, KIDSDRIV, HOMEKIDS)

desc <- psych::describe(num_preview)[, c("n","mean","sd","median","min","max","skew")]
kable(round(desc, 1), caption = "Summary statistics for key numeric variables")
Summary statistics for key numeric variables
n mean sd median min max skew
TARGET_FLAG 6528 0.3 0.4 0 0 1.0 1.1
TARGET_AMT 6528 1466.6 4545.7 0 0 107586.1 9.1
AGE 6525 44.9 8.7 45 16 81.0 0.0
YOJ 6153 10.5 4.1 11 0 23.0 -1.2
INCOME 6174 61649.3 47658.4 53276 0 367030.0 1.2
HOME_VAL 6160 154298.4 128874.6 160680 0 885282.0 0.5
BLUEBOOK 6528 15641.8 8381.5 14370 1500 69740.0 0.8
OLDCLAIM 6528 4119.3 8924.7 0 0 57037.0 3.1
TRAVTIME 6528 33.4 16.0 33 5 142.0 0.5
TIF 6528 5.4 4.1 4 1 25.0 0.9
CLM_FREQ 6528 0.8 1.2 0 0 5.0 1.2
MVR_PTS 6528 1.7 2.1 1 0 13.0 1.3
CAR_AGE 6129 8.3 5.7 8 -3 27.0 0.3
KIDSDRIV 6528 0.2 0.5 0 0 4.0 3.3
HOMEKIDS 6528 0.7 1.1 0 0 5.0 1.3

A few things jump out:

  • TARGET_AMT has a skew of roughly 9: almost all the mass sits at zero (no crash), with a long right tail of expensive claims. Any linear model on this variable in its raw form is going to fight the skew.
  • MVR_PTS, CLM_FREQ, and KIDSDRIV are also right-skewed count variables, which is expected: most drivers have zero moving violations, zero recent claims, and zero teenage drivers.
  • INCOME and HOME_VAL both have a meaningful mass of exact zeros (renters with HOME_VAL = 0, and some records with INCOME = 0), which is a legitimate value here, not obviously a missing-data placeholder.

Summary statistics: categorical variables

The numeric summary above doesn’t cover the 10 categorical predictors, so I build a matching table for them here: level counts, share of the sample, and, as a preview of what’s coming in the correlation and modeling sections, the raw crash rate within each level.

clean_disp <- function(x) str_remove(str_trim(x), "^z_")

cat_vars <- c("PARENT1","MSTATUS","SEX","EDUCATION","JOB","CAR_USE",
              "CAR_TYPE","RED_CAR","REVOKED","URBANICITY")

cat_summary <- lapply(cat_vars, function(v) {
  train_raw %>%
    mutate(level = clean_disp(.data[[v]]), level = ifelse(is.na(level) | level == "", "(Missing)", level)) %>%
    group_by(level) %>%
    summarise(N = n(), Pct = round(100 * n() / nrow(train_raw), 1),
              CrashRate = round(100 * mean(TARGET_FLAG), 1), .groups = "drop") %>%
    mutate(Variable = v, .before = 1)
}) %>% bind_rows()

kable(cat_summary, row.names = FALSE,
      caption = "Summary statistics for categorical variables: level, count, share of sample, and raw crash rate")
Summary statistics for categorical variables: level, count, share of sample, and raw crash rate
Variable level N Pct CrashRate
PARENT1 No 5672 86.9 23.7
PARENT1 Yes 856 13.1 44.3
MSTATUS No 2608 40.0 33.4
MSTATUS Yes 3920 60.0 21.7
SEX F 3493 53.5 27.2
SEX M 3035 46.5 25.4
EDUCATION <High School 973 14.9 32.0
EDUCATION Bachelors 1781 27.3 23.0
EDUCATION High School 1860 28.5 33.8
EDUCATION Masters 1331 20.4 20.4
EDUCATION PhD 583 8.9 17.2
JOB (Missing) 419 6.4 26.0
JOB Blue Collar 1452 22.2 35.2
JOB Clerical 1054 16.1 28.4
JOB Doctor 190 2.9 8.9
JOB Home Maker 512 7.8 28.1
JOB Lawyer 688 10.5 19.0
JOB Manager 758 11.6 14.2
JOB Professional 887 13.6 22.0
JOB Student 568 8.7 36.6
CAR_USE Commercial 2408 36.9 34.5
CAR_USE Private 4120 63.1 21.6
CAR_TYPE Minivan 1751 26.8 16.3
CAR_TYPE Panel Truck 526 8.1 25.9
CAR_TYPE Pickup 1106 16.9 32.3
CAR_TYPE SUV 1832 28.1 30.2
CAR_TYPE Sports Car 715 11.0 32.6
CAR_TYPE Van 598 9.2 26.1
RED_CAR no 4620 70.8 26.6
RED_CAR yes 1908 29.2 25.7
REVOKED No 5698 87.3 23.7
REVOKED Yes 830 12.7 44.5
URBANICITY Highly Rural/ Rural 1346 20.6 6.9
URBANICITY Highly Urban/ Urban 5182 79.4 31.4

A few takeaways that carry through to the models later: URBANICITY has by far the widest spread in crash rate across its two levels, single parents (PARENT1) and unmarried policyholders (MSTATUS) both show elevated crash rates, and JOB categories split fairly cleanly into “white collar” (Doctor, Lawyer, Manager, lower crash rate) versus everyone else, which lines up with the theoretical prior on the assignment sheet.

Distributions and relationship to the target

ggplot(train_raw %>% filter(TARGET_AMT > 0), aes(x = TARGET_AMT)) +
  geom_density(fill = "steelblue", alpha = 0.5, color = "steelblue4") +
  scale_x_continuous(labels = scales::comma) +
  labs(title = "Claim amount, conditional on a crash", x = "TARGET_AMT ($)", y = "Density") +
  theme_minimal()
Distribution of claim amount among crashes

Distribution of claim amount among crashes 

ggplot(train_raw %>% filter(TARGET_AMT > 0), aes(x = log(TARGET_AMT))) +
  geom_density(fill = "darkorange", alpha = 0.5, color = "darkorange4") +
  labs(title = "log(TARGET_AMT), conditional on a crash", x = "log(TARGET_AMT)", y = "Density") +
  theme_minimal()
Log-claim amount among crashes

Log-claim amount among crashes

The log transform pulls the claim-amount distribution much closer to symmetric, which is the main justification for modeling log(TARGET_AMT) rather than the raw dollar figure later on.

ggplot(train_raw, aes(x = factor(TARGET_FLAG), y = MVR_PTS)) +
  geom_boxplot(fill = c("seagreen","firebrick"), alpha = 0.6) +
  labs(title = "MVR points vs. crash outcome", x = "TARGET_FLAG (1 = crash)", y = "MVR_PTS") +
  theme_minimal()
Motor vehicle record points by crash status

Motor vehicle record points by crash status 

train_raw %>%
  group_by(URBANICITY) %>%
  summarise(crash_rate = mean(TARGET_FLAG), n = n()) %>%
  ggplot(aes(x = URBANICITY, y = crash_rate)) +
  geom_col(fill = "slateblue") +
  labs(title = "Crash rate by home/work area", y = "Crash rate", x = NULL) +
  theme_minimal()
Crash rate by urbanicity

Crash rate by urbanicity

Urban policyholders crash far more often than rural ones, and this turns out to be one of the strongest single predictors in the whole data set, well beyond what I expected going in.

Correlation with the targets

corr_all <- cor(num_preview, use = "complete.obs")
corrplot(corr_all, method = "color", type = "upper", tl.cex = 0.75, tl.col = "black",
         diag = FALSE, title = "Correlation matrix: targets and numeric predictors",
         mar = c(0,0,2,0))
Correlation heatmap of numeric variables

Correlation heatmap of numeric variables

The exact values behind this heatmap, including significance markers, are in the table just below; the plot is here to make the overall pattern (strong blue cluster in the top-left among the targets, weaker and more scattered elsewhere) easy to see at a glance without needing to scan 15 rows of numbers.

sig_stars <- function(p) {
  case_when(p < 0.001 ~ "***", p < 0.01 ~ "**", p < 0.05 ~ "*", TRUE ~ "")
}
add_sig <- function(varname, target) {
  p <- cor.test(num_preview[[varname]], num_preview[[target]])$p.value
  sig_stars(p)
}

corr_vars <- setdiff(colnames(num_preview), c("TARGET_FLAG", "TARGET_AMT"))
corr_tbl <- data.frame(
  Variable = corr_vars,
  With_TARGET_FLAG = round(corr_all[corr_vars, "TARGET_FLAG"], 3),
  Sig_FLAG = sapply(corr_vars, add_sig, target = "TARGET_FLAG"),
  With_TARGET_AMT  = round(corr_all[corr_vars, "TARGET_AMT"], 3),
  Sig_AMT = sapply(corr_vars, add_sig, target = "TARGET_AMT")
) %>% arrange(With_TARGET_FLAG)
rownames(corr_tbl) <- NULL
kable(corr_tbl, row.names = FALSE,
      caption = "Correlation of numeric predictors with each target variable (*** p<.001, ** p<.01, * p<.05)")
Correlation of numeric predictors with each target variable (*** p<.001, ** p<.01, * p<.05)
Variable With_TARGET_FLAG Sig_FLAG With_TARGET_AMT Sig_AMT
HOME_VAL -0.168 *** -0.083 ***
INCOME -0.131 *** -0.052 ***
AGE -0.113 *** -0.052 **
BLUEBOOK -0.108 *** 0.001
CAR_AGE -0.103 *** -0.056 ***
TIF -0.078 *** -0.036 ***
YOJ -0.067 *** -0.028 *
TRAVTIME 0.044 *** 0.027 *
KIDSDRIV 0.110 *** 0.057 ***
HOMEKIDS 0.119 *** 0.066 ***
OLDCLAIM 0.144 *** 0.081 ***
CLM_FREQ 0.222 *** 0.130 ***
MVR_PTS 0.227 *** 0.146 ***

Reading the table: MVR_PTS and CLM_FREQ (past behavior) correlate positively with crashing again, which is exactly the Becker-style “past record predicts future record” logic the assignment sheet flags. HOME_VAL, INCOME, and vehicle AGE/BLUEBOOK all correlate negatively: wealthier, older-car households crash less, consistent with the theoretical priors in the assignment. With a training sample this large (n = 6,528), even fairly small correlations come back statistically significant; every variable with |r| > 0.04 against TARGET_FLAG clears p < .001, so significance alone isn’t the interesting part of this table so much as the pattern and relative ranking of the correlations is. The same variables correlate with TARGET_AMT in the same direction but much more weakly (BLUEBOOK doesn’t even clear p < .05 against TARGET_AMT, despite being a strong, highly significant predictor of TARGET_FLAG), which foreshadows a finding I come back to later: these variables predict whether a crash happens far better than they predict how expensive it is. Off the diagonal, the strongest predictor-predictor correlation is between HOME_VAL and INCOME (wealthier households own more expensive homes), which is worth keeping in mind for multicollinearity even though, as shown later, it doesn’t turn out to inflate any variance inflation factors past standard thresholds.

Data Preparation

I write a single prep_data() function and apply it identically to the training and evaluation sets, so whatever I do to one happens to the other. Imputation values (medians, modal job) are learned on the training set only and then reused on the evaluation set, to avoid leaking evaluation-set information into the model.

Cleaning currency fields and category labels

parse_money <- function(x) as.numeric(gsub("[$,]", "", x))

clean_base <- function(df) {
  df <- df %>%
    mutate(
      INCOME    = parse_money(INCOME),
      HOME_VAL  = parse_money(HOME_VAL),
      BLUEBOOK  = parse_money(BLUEBOOK),
      OLDCLAIM  = parse_money(OLDCLAIM),
      PARENT1    = str_trim(PARENT1),
      MSTATUS    = str_remove(str_trim(MSTATUS), "^z_"),
      SEX        = str_remove(str_trim(SEX), "^z_"),
      EDUCATION  = str_remove(str_trim(EDUCATION), "^z_"),
      EDUCATION  = ifelse(EDUCATION == "<High School", "Under High School", EDUCATION),
      JOB        = str_remove(str_trim(JOB), "^z_"),
      JOB        = na_if(JOB, ""),
      CAR_USE    = str_trim(CAR_USE),
      CAR_TYPE   = str_remove(str_trim(CAR_TYPE), "^z_"),
      RED_CAR    = str_trim(RED_CAR),
      REVOKED    = str_trim(REVOKED),
      URBANICITY = str_remove(str_trim(URBANICITY), "^z_"),
      URBANICITY = case_when(
        URBANICITY == "Highly Rural/ Rural" ~ "Rural",
        URBANICITY == "Highly Urban/ Urban" ~ "Urban",
        TRUE ~ URBANICITY)
    )
  # a negative vehicle age is a data-entry error, not a real value
  df$CAR_AGE[df$CAR_AGE < 0 & !is.na(df$CAR_AGE)] <- NA
  df
}

train_c <- clean_base(train_raw)
test_c  <- clean_base(test_raw)

Missing-value flags, then imputation

Per the assignment suggestions, I create a “was this missing” flag before imputing, so the model can pick up any signal in missingness itself, separate from whatever value I plug in.

miss_vars <- c("AGE","YOJ","INCOME","HOME_VAL","CAR_AGE","JOB")

add_miss_flags <- function(df) {
  for (v in miss_vars) {
    df[[paste0("M_", v)]] <- as.integer(is.na(df[[v]]))
  }
  df
}
train_c <- add_miss_flags(train_c)
test_c  <- add_miss_flags(test_c)

# learn imputation values on TRAIN only
med_age      <- median(train_c$AGE, na.rm = TRUE)
med_yoj      <- median(train_c$YOJ, na.rm = TRUE)
med_income   <- median(train_c$INCOME, na.rm = TRUE)
med_homeval  <- median(train_c$HOME_VAL, na.rm = TRUE)
med_carage   <- median(train_c$CAR_AGE, na.rm = TRUE)
mode_job     <- names(sort(table(train_c$JOB), decreasing = TRUE))[1]

impute_all <- function(df) {
  df$AGE[is.na(df$AGE)]           <- med_age
  df$YOJ[is.na(df$YOJ)]           <- med_yoj
  df$INCOME[is.na(df$INCOME)]     <- med_income
  df$HOME_VAL[is.na(df$HOME_VAL)] <- med_homeval
  df$CAR_AGE[is.na(df$CAR_AGE)]   <- med_carage
  df$JOB[is.na(df$JOB)]           <- mode_job
  df
}
train_c <- impute_all(train_c)
test_c  <- impute_all(test_c)

cat("Imputed with -> AGE:", med_age, " YOJ:", med_yoj, "\n",
    "  INCOME:", med_income, " HOME_VAL:", med_homeval, "\n",
    "  CAR_AGE:", med_carage, " JOB:", mode_job, "\n")
## Imputed with -> AGE: 45  YOJ: 11 
##    INCOME: 53276  HOME_VAL: 160680 
##    CAR_AGE: 8  JOB: Blue Collar

Note on M_AGE: only 3 training records are missing AGE. That’s too few for a stable coefficient on the flag (it produces a huge standard error / quasi-separation if included), so I impute AGE but drop M_AGE from the candidate variable list in modeling below.

Engineered variables and transformations

engineer <- function(df) {
  df %>%
    mutate(
      HOME_OWNER      = as.integer(HOME_VAL > 0),
      HAS_TEEN_DRIVER = as.integer(KIDSDRIV > 0),
      LOG_INCOME      = log1p(INCOME),
      LOG_HOME_VAL    = log1p(HOME_VAL),
      LOG_BLUEBOOK    = log(BLUEBOOK),
      LOG_OLDCLAIM    = log1p(OLDCLAIM),
      LOG_TARGET_AMT  = log1p(TARGET_AMT)
    )
}
train_c <- engineer(train_c)
test_c  <- engineer(test_c)
  • HOME_OWNER collapses the very right-skewed HOME_VAL into a simple owns-a-home flag, on top of keeping the continuous log-transformed value.
  • HAS_TEEN_DRIVER mirrors KIDSDRIV, letting a model use whichever form (continuous count vs. any-teen-driver-at-all) fits better.
  • LOG_INCOME, LOG_HOME_VAL, LOG_BLUEBOOK, and LOG_OLDCLAIM address the skew documented above; log1p (i.e. log(1+x)) handles the variables that can legitimately be zero, while BLUEBOOK (never zero) gets a plain log.

Categorical variables as factors with a sensible reference level

Rather than hand-building dummy columns, I let R’s factor/glm/lm machinery do it, but I explicitly set the reference level for each factor so the coefficients have an intuitive baseline (e.g., “Blue Collar” worker, “Minivan” driver, non-urban).

set_factors <- function(df) {
  df$JOB        <- relevel(factor(df$JOB), ref = "Blue Collar")
  df$EDUCATION  <- relevel(factor(df$EDUCATION), ref = "Under High School")
  df$CAR_TYPE   <- relevel(factor(df$CAR_TYPE), ref = "Minivan")
  df$PARENT1    <- relevel(factor(df$PARENT1), ref = "No")
  df$MSTATUS    <- relevel(factor(df$MSTATUS), ref = "No")
  df$SEX        <- relevel(factor(df$SEX), ref = "F")
  df$RED_CAR    <- relevel(factor(df$RED_CAR), ref = "no")
  df$REVOKED    <- relevel(factor(df$REVOKED), ref = "No")
  df$CAR_USE    <- relevel(factor(df$CAR_USE), ref = "Private")
  df$URBANICITY <- relevel(factor(df$URBANICITY), ref = "Rural")
  df
}
train_p <- set_factors(train_c)
test_p  <- set_factors(test_c)

sapply(train_p[c("JOB","EDUCATION","CAR_TYPE","URBANICITY")], levels)
## $JOB
## [1] "Blue Collar"  "Clerical"     "Doctor"       "Home Maker"  
## [5] "Lawyer"       "Manager"      "Professional" "Student"     
## 
## $EDUCATION
## [1] "Under High School" "Bachelors"         "High School"      
## [4] "Masters"           "PhD"              
## 
## $CAR_TYPE
## [1] "Minivan"     "Panel Truck" "Pickup"      "Sports Car" 
## [5] "SUV"         "Van"        
## 
## $URBANICITY
## [1] "Rural" "Urban"

Final prepared training set: 6528 rows, 39 columns.

Updated summary statistics, post-cleaning

Re-running the same summary statistics from Data Exploration confirms the missing-value and skewness fixes actually landed. Every variable in this table now has the full n = 6528 (no more NAs pulling the count down), so I drop the n column entirely rather than repeat the same number in every row, and the log-transformed dollar variables show visibly lower skew than their raw versions.

post_vars <- train_p %>%
  transmute(AGE, YOJ, INCOME, LOG_INCOME, HOME_VAL, LOG_HOME_VAL,
            BLUEBOOK, LOG_BLUEBOOK, OLDCLAIM, LOG_OLDCLAIM, CAR_AGE)
post_desc <- psych::describe(post_vars)[, c("mean","sd","median","min","max","skew")]
kable(round(post_desc, 1),
      caption = paste0("Post-cleaning summary statistics, imputed + transformed (n = ", nrow(train_p), " for every variable)"))
Post-cleaning summary statistics, imputed + transformed (n = 6528 for every variable)
mean sd median min max skew
AGE 44.9 8.6 45.0 16.0 81.0 0.0
YOJ 10.5 4.0 11.0 0.0 23.0 -1.3
INCOME 61195.3 46386.8 53276.0 0.0 367030.0 1.3
LOG_INCOME 10.0 3.0 10.9 0.0 12.8 -2.8
HOME_VAL 154658.2 125197.5 160680.0 0.0 885282.0 0.5
LOG_HOME_VAL 8.8 5.5 12.0 0.0 13.7 -1.0
BLUEBOOK 15641.8 8381.5 14370.0 1500.0 69740.0 0.8
LOG_BLUEBOOK 9.5 0.6 9.6 7.3 11.2 -0.9
OLDCLAIM 4119.3 8924.7 0.0 0.0 57037.0 3.1
LOG_OLDCLAIM 3.4 4.3 0.0 0.0 11.0 0.5
CAR_AGE 8.3 5.5 8.0 0.0 27.0 0.3

INCOME skew drops from roughly 1.2 to well under 1 on the log scale, and HOME_VAL and OLDCLAIM show the same pattern, so the log transforms are doing what they’re supposed to.

Build Models

I build three binary logistic regression models for TARGET_FLAG and two multiple linear regression models (plus one supplementary severity model) for TARGET_AMT, then add a bonus LASSO/Ridge/probit check on top of the logistic side. All three logistic specifications share the same estimating form, just with a different right-hand side; letting \(\pi_i = P(\text{TARGET\_FLAG}_i = 1)\), every logit model below is a version of:

\[\text{logit}(\pi_i) = \ln\!\left(\frac{\pi_i}{1-\pi_i}\right) = \beta_0 + \sum_{j=1}^{k}\beta_j X_{ji}\]

where the \(X_j\) are whatever that model’s specific right-hand-side variables are.

Logistic Model 1: full “kitchen sink”

Every candidate predictor, untransformed, so I can see the whole field before narrowing anything down. Written out explicitly (grouping the ten categorical predictors and five missingness flags into blocks of dummy variables rather than one beta per factor level, since several are multi-level factors):

\[ \begin{aligned} \text{logit}(\pi_i) = \;& \beta_0 + \beta_1 \text{KIDSDRIV}_i + \beta_2 \text{AGE}_i + \beta_3 \text{HOMEKIDS}_i + \beta_4 \text{YOJ}_i \\ & + \beta_5 \text{INCOME}_i + \beta_6 \text{HOME\_VAL}_i + \beta_7 \text{TRAVTIME}_i + \beta_8 \text{BLUEBOOK}_i \\ & + \beta_9 \text{TIF}_i + \beta_{10} \text{OLDCLAIM}_i + \beta_{11} \text{CLM\_FREQ}_i + \beta_{12} \text{MVR\_PTS}_i + \beta_{13}\text{CAR\_AGE}_i \\ & + \underbrace{\boldsymbol\beta_{14}^\top \mathbf{D}^{cat}_i}_{\text{10 categorical blocks}} + \underbrace{\boldsymbol\beta_{15}^\top \mathbf{D}^{miss}_i}_{\text{5 missingness flags}} + \varepsilon_i \end{aligned} \]

where \(\mathbf{D}^{cat}_i\) collects the dummy variables for PARENT1, MSTATUS, SEX, RED_CAR, REVOKED, URBANICITY, CAR_USE, EDUCATION, JOB, and CAR_TYPE, and \(\mathbf{D}^{miss}_i\) collects the five missingness flags (M_YOJ, M_INCOME, M_HOME_VAL, M_CAR_AGE, M_JOB). That block notation is standard practice once a model has this many dummy-coded factors; writing out all 28+ individual coefficient terms by hand would be less readable, not more informative, than the R formula in the code below, which is the literal, unambiguous version of the same equation.

logit1 <- glm(TARGET_FLAG ~ KIDSDRIV + AGE + HOMEKIDS + YOJ + INCOME + HOME_VAL +
                TRAVTIME + BLUEBOOK + TIF + OLDCLAIM + CLM_FREQ + MVR_PTS + CAR_AGE +
                PARENT1 + MSTATUS + SEX + RED_CAR + REVOKED + URBANICITY + CAR_USE +
                EDUCATION + JOB + CAR_TYPE + M_YOJ + M_INCOME + M_HOME_VAL + M_CAR_AGE + M_JOB,
              data = train_p, family = binomial)
summary(logit1)
## 
## Call:
## glm(formula = TARGET_FLAG ~ KIDSDRIV + AGE + HOMEKIDS + YOJ + 
##     INCOME + HOME_VAL + TRAVTIME + BLUEBOOK + TIF + OLDCLAIM + 
##     CLM_FREQ + MVR_PTS + CAR_AGE + PARENT1 + MSTATUS + SEX + 
##     RED_CAR + REVOKED + URBANICITY + CAR_USE + EDUCATION + JOB + 
##     CAR_TYPE + M_YOJ + M_INCOME + M_HOME_VAL + M_CAR_AGE + M_JOB, 
##     family = binomial, data = train_p)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -3.297e+00  3.141e-01 -10.495  < 2e-16 ***
## KIDSDRIV              4.452e-01  6.939e-02   6.416 1.40e-10 ***
## AGE                  -1.717e-03  4.516e-03  -0.380 0.703737    
## HOMEKIDS              4.371e-02  4.227e-02   1.034 0.301144    
## YOJ                  -1.424e-02  9.588e-03  -1.485 0.137522    
## INCOME               -3.146e-06  1.207e-06  -2.606 0.009149 ** 
## HOME_VAL             -1.092e-06  3.835e-07  -2.848 0.004397 ** 
## TRAVTIME              1.345e-02  2.103e-03   6.396 1.60e-10 ***
## BLUEBOOK             -1.999e-05  5.905e-06  -3.386 0.000710 ***
## TIF                  -5.648e-02  8.234e-03  -6.859 6.91e-12 ***
## OLDCLAIM             -1.535e-05  4.287e-06  -3.580 0.000343 ***
## CLM_FREQ              2.046e-01  3.180e-02   6.434 1.24e-10 ***
## MVR_PTS               1.158e-01  1.525e-02   7.592 3.16e-14 ***
## CAR_AGE               3.995e-03  8.433e-03   0.474 0.635691    
## PARENT1Yes            3.981e-01  1.230e-01   3.236 0.001213 ** 
## MSTATUSYes           -4.988e-01  9.393e-02  -5.311 1.09e-07 ***
## SEXM                  1.571e-01  1.253e-01   1.254 0.209832    
## RED_CARyes           -1.374e-02  9.648e-02  -0.142 0.886721    
## REVOKEDYes            9.378e-01  1.003e-01   9.351  < 2e-16 ***
## URBANICITYUrban       2.349e+00  1.251e-01  18.775  < 2e-16 ***
## CAR_USECommercial     7.540e-01  1.036e-01   7.277 3.42e-13 ***
## EDUCATIONBachelors   -4.585e-01  1.290e-01  -3.554 0.000380 ***
## EDUCATIONHigh School -1.250e-02  1.062e-01  -0.118 0.906263    
## EDUCATIONMasters     -3.022e-01  1.990e-01  -1.518 0.128937    
## EDUCATIONPhD         -1.706e-01  2.363e-01  -0.722 0.470313    
## JOBClerical           3.291e-02  1.196e-01   0.275 0.783261    
## JOBDoctor            -1.198e+00  3.448e-01  -3.473 0.000515 ***
## JOBHome Maker        -1.265e-01  1.718e-01  -0.737 0.461415    
## JOBLawyer            -3.029e-01  2.081e-01  -1.455 0.145614    
## JOBManager           -8.963e-01  1.588e-01  -5.645 1.65e-08 ***
## JOBProfessional      -1.971e-01  1.353e-01  -1.456 0.145281    
## JOBStudent           -1.259e-01  1.444e-01  -0.872 0.383313    
## CAR_TYPEPanel Truck   4.441e-01  1.822e-01   2.437 0.014801 *  
## CAR_TYPEPickup        5.561e-01  1.120e-01   4.963 6.94e-07 ***
## CAR_TYPESports Car    1.029e+00  1.460e-01   7.048 1.82e-12 ***
## CAR_TYPESUV           8.337e-01  1.246e-01   6.693 2.18e-11 ***
## CAR_TYPEVan           5.528e-01  1.417e-01   3.902 9.53e-05 ***
## M_YOJ                 1.325e-01  1.381e-01   0.959 0.337574    
## M_INCOME             -8.301e-02  1.465e-01  -0.567 0.570928    
## M_HOME_VAL            1.367e-02  1.390e-01   0.098 0.921654    
## M_CAR_AGE             1.250e-01  1.344e-01   0.930 0.352591    
## M_JOB                -3.594e-01  2.082e-01  -1.726 0.084306 .  
## ---
## 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: 5821.4  on 6486  degrees of freedom
## AIC: 5905.4
## 
## Number of Fisher Scoring iterations: 5
cat("AIC:", AIC(logit1), " LogLik:", as.numeric(logLik(logit1)), "\n")
## AIC: 5905.352  LogLik: -2910.676

Logistic Model 2: backward-stepwise (AIC)

Starting from Model 1, I let MASS::stepAIC remove variables one at a time as long as doing so improves AIC. This is the “let the algorithm narrow it down” model, contrasted with Model 3’s manual, theory-first approach below. Its estimating equation has the identical functional form shown above for Model 1, just with whichever terms stepAIC keeps on the right-hand side; since that subset is only known after fitting, I show the actual kept-variable list as part of the output below rather than writing it out twice.

logit2 <- stepAIC(logit1, direction = "backward", trace = FALSE)
summary(logit2)
## 
## Call:
## glm(formula = TARGET_FLAG ~ KIDSDRIV + INCOME + HOME_VAL + TRAVTIME + 
##     BLUEBOOK + TIF + OLDCLAIM + CLM_FREQ + MVR_PTS + PARENT1 + 
##     MSTATUS + REVOKED + URBANICITY + CAR_USE + EDUCATION + JOB + 
##     CAR_TYPE + M_JOB, family = binomial, data = train_p)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -3.319e+00  2.184e-01 -15.197  < 2e-16 ***
## KIDSDRIV              4.699e-01  6.218e-02   7.556 4.16e-14 ***
## INCOME               -3.257e-06  1.199e-06  -2.716 0.006607 ** 
## HOME_VAL             -1.133e-06  3.808e-07  -2.976 0.002924 ** 
## TRAVTIME              1.334e-02  2.099e-03   6.358 2.05e-10 ***
## BLUEBOOK             -2.398e-05  5.317e-06  -4.510 6.48e-06 ***
## TIF                  -5.662e-02  8.224e-03  -6.884 5.82e-12 ***
## OLDCLAIM             -1.552e-05  4.283e-06  -3.625 0.000289 ***
## CLM_FREQ              2.053e-01  3.176e-02   6.463 1.03e-10 ***
## MVR_PTS               1.169e-01  1.520e-02   7.688 1.49e-14 ***
## PARENT1Yes            4.679e-01  1.060e-01   4.415 1.01e-05 ***
## MSTATUSYes           -4.813e-01  8.926e-02  -5.393 6.94e-08 ***
## REVOKEDYes            9.446e-01  1.001e-01   9.439  < 2e-16 ***
## URBANICITYUrban       2.351e+00  1.251e-01  18.794  < 2e-16 ***
## CAR_USECommercial     7.504e-01  1.032e-01   7.275 3.48e-13 ***
## EDUCATIONBachelors   -4.367e-01  1.213e-01  -3.602 0.000316 ***
## EDUCATIONHigh School -7.829e-03  1.057e-01  -0.074 0.940938    
## EDUCATIONMasters     -2.629e-01  1.793e-01  -1.466 0.142611    
## EDUCATIONPhD         -1.356e-01  2.207e-01  -0.614 0.538901    
## JOBClerical           3.041e-02  1.193e-01   0.255 0.798795    
## JOBDoctor            -1.187e+00  3.438e-01  -3.454 0.000552 ***
## JOBHome Maker        -8.201e-02  1.623e-01  -0.505 0.613341    
## JOBLawyer            -3.099e-01  2.074e-01  -1.494 0.135183    
## JOBManager           -9.029e-01  1.583e-01  -5.704 1.17e-08 ***
## JOBProfessional      -2.049e-01  1.349e-01  -1.520 0.128623    
## JOBStudent           -5.754e-02  1.375e-01  -0.418 0.675601    
## CAR_TYPEPanel Truck   5.389e-01  1.706e-01   3.158 0.001589 ** 
## CAR_TYPEPickup        5.494e-01  1.117e-01   4.918 8.76e-07 ***
## CAR_TYPESports Car    9.220e-01  1.206e-01   7.643 2.12e-14 ***
## CAR_TYPESUV           7.289e-01  9.559e-02   7.625 2.44e-14 ***
## CAR_TYPEVan           6.032e-01  1.371e-01   4.399 1.09e-05 ***
## M_JOB                -3.618e-01  2.079e-01  -1.740 0.081822 .  
## ---
## 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: 5828.7  on 6496  degrees of freedom
## AIC: 5892.7
## 
## Number of Fisher Scoring iterations: 5
cat("AIC:", AIC(logit2), " LogLik:", as.numeric(logLik(logit2)), "\n")
## AIC: 5892.661  LogLik: -2914.331
cat("Variables retained:", length(coef(logit2)) - 1, "\n")
## Variables retained: 31

Logistic Model 3: parsimonious, theory-driven, log-transformed

Instead of algorithmic selection, I hand-pick a small set of variables that the assignment’s theoretical-effect column flags as the strongest candidates, and I use the log-transformed versions of the skewed dollar variables instead of their raw dollar form. Small enough this time to write out every single term:

\[ \begin{aligned} \text{logit}(\pi_i) = \;& \beta_0 + \beta_1 \text{KIDSDRIV}_i + \beta_2 \text{MVR\_PTS}_i + \beta_3 \text{CLM\_FREQ}_i + \beta_4 \text{REVOKED}_i \\ & + \beta_5 \text{URBANICITY}_i + \beta_6 \text{CAR\_USE}_i + \beta_7 \text{LOG\_INCOME}_i + \beta_8 \text{HOME\_OWNER}_i \\ & + \beta_9 \text{LOG\_BLUEBOOK}_i + \beta_{10} \text{TIF}_i + \beta_{11} \text{PARENT1}_i + \beta_{12} \text{MSTATUS}_i + \varepsilon_i \end{aligned} \]

logit3 <- glm(TARGET_FLAG ~ KIDSDRIV + MVR_PTS + CLM_FREQ + REVOKED + URBANICITY +
                CAR_USE + LOG_INCOME + HOME_OWNER + LOG_BLUEBOOK + TIF +
                PARENT1 + MSTATUS,
              data = train_p, family = binomial)
summary(logit3)
## 
## Call:
## glm(formula = TARGET_FLAG ~ KIDSDRIV + MVR_PTS + CLM_FREQ + REVOKED + 
##     URBANICITY + CAR_USE + LOG_INCOME + HOME_OWNER + LOG_BLUEBOOK + 
##     TIF + PARENT1 + MSTATUS, family = binomial, data = train_p)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        2.58585    0.48653   5.315 1.07e-07 ***
## KIDSDRIV           0.46597    0.06029   7.729 1.09e-14 ***
## MVR_PTS            0.12338    0.01482   8.327  < 2e-16 ***
## CLM_FREQ           0.16906    0.02754   6.138 8.37e-10 ***
## REVOKEDYes         0.77203    0.08552   9.028  < 2e-16 ***
## URBANICITYUrban    2.06050    0.12428  16.579  < 2e-16 ***
## CAR_USECommercial  0.88093    0.06639  13.269  < 2e-16 ***
## LOG_INCOME        -0.08220    0.01092  -7.527 5.19e-14 ***
## HOME_OWNER        -0.27372    0.08410  -3.255  0.00114 ** 
## LOG_BLUEBOOK      -0.52521    0.05257  -9.991  < 2e-16 ***
## TIF               -0.05039    0.00798  -6.315 2.70e-10 ***
## PARENT1Yes         0.51817    0.10215   5.072 3.93e-07 ***
## MSTATUSYes        -0.34923    0.08671  -4.028 5.63e-05 ***
## ---
## 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: 6121.6  on 6515  degrees of freedom
## AIC: 6147.6
## 
## Number of Fisher Scoring iterations: 5
cat("AIC:", AIC(logit3), " LogLik:", as.numeric(logLik(logit3)), "\n")
## AIC: 6147.591  LogLik: -3060.796

Checking multicollinearity (Model 3)

vif(logit3)
##     KIDSDRIV      MVR_PTS     CLM_FREQ      REVOKED   URBANICITY 
##     1.083902     1.144096     1.159780     1.002933     1.097885 
##      CAR_USE   LOG_INCOME   HOME_OWNER LOG_BLUEBOOK          TIF 
##     1.078922     1.163527     1.541232     1.147466     1.005685 
##      PARENT1      MSTATUS 
##     1.401222     1.863930

All VIFs are comfortably below 5, so multicollinearity isn’t a concern for the parsimonious model.

Bonus: penalized logistic regression (LASSO, Ridge) and a probit alternative

The assignment allows substituting a ridge or LASSO fit for one of the required models as a bonus. Rather than replace anything, I fit both alongside the three required logistic models, using the same full variable set as Model 1, plus a probit link on Model 2’s selected variables as a second, smaller bonus check.

x_form <- TARGET_FLAG ~ KIDSDRIV + AGE + HOMEKIDS + YOJ + INCOME + HOME_VAL +
  TRAVTIME + BLUEBOOK + TIF + OLDCLAIM + CLM_FREQ + MVR_PTS + CAR_AGE +
  PARENT1 + MSTATUS + SEX + RED_CAR + REVOKED + URBANICITY + CAR_USE +
  EDUCATION + JOB + CAR_TYPE + M_YOJ + M_INCOME + M_HOME_VAL + M_CAR_AGE + M_JOB

x_train <- model.matrix(x_form, data = train_p)[, -1]
y_train <- train_p$TARGET_FLAG

set.seed(7320)
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")
kept <- setdiff(rownames(lasso_coef)[which(lasso_coef[, 1] != 0)], "(Intercept)")
cat("LASSO (lambda.1se) keeps", length(kept), "of", ncol(x_train), "candidate predictors.\n")
## LASSO (lambda.1se) keeps 23 of 41 candidate predictors.
probit1 <- glm(formula(logit2), data = train_p, family = binomial(link = "probit"))

LASSO’s own feature selection is a useful cross-check on the stepwise selection in Model 2: if the two methods agree on which variables matter, that’s good evidence Model 2 isn’t just an artifact of AIC’s particular penalty.

ridge_coef <- coef(cv_ridge, s = "lambda.1se")
compare_tbl <- data.frame(
  Variable   = kept,
  LASSO_coef = round(as.numeric(lasso_coef[kept, 1]), 4),
  Ridge_coef = round(as.numeric(ridge_coef[kept, 1]), 4)
)
kable(compare_tbl, row.names = FALSE,
      caption = "Variables retained by LASSO (lambda.1se) and their Ridge counterparts")
Variables retained by LASSO (lambda.1se) and their Ridge counterparts
Variable LASSO_coef Ridge_coef
KIDSDRIV 0.3140 0.3062
HOMEKIDS 0.0231 0.0501
INCOME 0.0000 0.0000
HOME_VAL 0.0000 0.0000
TRAVTIME 0.0075 0.0082
BLUEBOOK 0.0000 0.0000
TIF -0.0341 -0.0386
CLM_FREQ 0.1366 0.1710
MVR_PTS 0.0985 0.1004
CAR_AGE -0.0062 -0.0052
PARENT1Yes 0.3804 0.3449
MSTATUSYes -0.2846 -0.3330
REVOKEDYes 0.6011 0.6565
URBANICITYUrban 1.7497 1.3795
CAR_USECommercial 0.6691 0.5351
EDUCATIONBachelors -0.0796 -0.2123
EDUCATIONHigh School 0.1210 0.1100
EDUCATIONMasters -0.0162 -0.1487
JOBDoctor -0.2916 -0.6183
JOBManager -0.5051 -0.5467
CAR_TYPEPickup 0.0405 0.2696
CAR_TYPESports Car 0.3054 0.4599
CAR_TYPESUV 0.2292 0.3540

Every variable LASSO keeps also survived the backward-AIC selection in Model 2 (motor-vehicle points, claim history, revoked license, urbanicity, commercial use, and several of the same job/car-type dummies), which is reassuring: two quite different selection procedures, one penalty-based and one information-criterion-based, land on essentially the same set of predictors.

Side-by-side comparison of all three logistic models

Rather than four separate walls of summary() output, stargazer puts all four specifications (including the bonus probit) in one compact table so coefficient sign, magnitude, and significance can be compared directly across models.

stargazer(logit1, logit2, logit3, probit1, type = "latex", header = FALSE,
          title = "Logistic (and probit) regression models: crash probability (TARGET\\_FLAG)",
          column.labels = c("Model 1: Full","Model 2: Stepwise","Model 3: Parsimonious","Bonus: Probit"),
          dep.var.labels = "P(crash)", single.row = TRUE, no.space = TRUE,
          font.size = "tiny", omit.stat = c("bic"),
          star.cutoffs = c(0.05, 0.01, 0.001))

The probit column uses the exact same variables as Model 2, just with a different link function, and the coefficients tell the same story on a rescaled axis (probit coefficients are always smaller in magnitude than the corresponding logit ones, roughly by a factor of 1.6 to 1.8, since the two links are fit on different underlying scales), with the same signs and essentially the same significance pattern throughout.

Coefficient evolution across the three specifications

To illustrate omitted-variable bias, the table below tracks a few variables that appear in all three models as the specification grows from the parsimonious model to the full kitchen-sink model.

# pull coefficients by name, matching each model's actual dummy-coded names
get_coef <- function(model, name) {
  cn <- names(coef(model))
  hit <- cn[grepl(name, cn, fixed = TRUE)]
  if (length(hit) == 0) return(NA)
  unname(coef(model)[hit[1]])
}
coef_table <- data.frame(
  Variable = c("MVR_PTS", "CLM_FREQ", "URBANICITY (Urban)", "CAR_USE (Commercial)"),
  Model3_parsimonious = c(get_coef(logit3,"MVR_PTS"), get_coef(logit3,"CLM_FREQ"),
                          get_coef(logit3,"URBANICITY"), get_coef(logit3,"CAR_USE")),
  Model2_stepwise     = c(get_coef(logit2,"MVR_PTS"), get_coef(logit2,"CLM_FREQ"),
                          get_coef(logit2,"URBANICITY"), get_coef(logit2,"CAR_USE")),
  Model1_full         = c(get_coef(logit1,"MVR_PTS"), get_coef(logit1,"CLM_FREQ"),
                          get_coef(logit1,"URBANICITY"), get_coef(logit1,"CAR_USE"))
)
kable(coef_table, digits = 3, row.names = FALSE,
      caption = "Coefficient stability across specifications (logit, log-odds scale)")
Coefficient stability across specifications (logit, log-odds scale)
Variable Model3_parsimonious Model2_stepwise Model1_full
MVR_PTS 0.123 0.117 0.116
CLM_FREQ 0.169 0.205 0.205
URBANICITY (Urban) 2.060 2.351 2.349
CAR_USE (Commercial) 0.881 0.750 0.754

The coefficients on MVR_PTS, CLM_FREQ, URBANICITY, and CAR_USE are remarkably stable across all three specifications: they don’t collapse or flip sign as I add controls, which is reassuring: it suggests these effects aren’t artifacts of omitted-variable bias from the handful of variables each specification leaves out.

Linear Model 1: raw TARGET_AMT, full variable set

The naive approach: regress the dollar claim amount directly on everything, using every record (crash or not). Same variable set as logistic Model 1 (minus the missingness flags, which weren’t significant predictors of cost in early exploration and are dropped here to keep the two model families easier to compare side by side), estimated by OLS instead of maximum likelihood:

\[ \text{TARGET\_AMT}_i = \beta_0 + \beta_1 \text{KIDSDRIV}_i + \beta_2 \text{AGE}_i + \dots + \boldsymbol\beta^\top \mathbf{D}^{cat}_i + \varepsilon_i, \quad \varepsilon_i \sim (0, \sigma^2) \]

using the identical right-hand-side variable block notation as Model 1’s logit equation above. Model 2 below is the same equation with \(\ln(1+\text{TARGET\_AMT}_i)\) on the left-hand side instead of the raw dollar figure.

mlr_vars_formula <- TARGET_AMT ~ KIDSDRIV + AGE + HOMEKIDS + YOJ + INCOME + HOME_VAL +
  TRAVTIME + BLUEBOOK + TIF + OLDCLAIM + CLM_FREQ + MVR_PTS + CAR_AGE +
  PARENT1 + MSTATUS + SEX + RED_CAR + REVOKED + URBANICITY + CAR_USE +
  EDUCATION + JOB + CAR_TYPE

mlr1 <- lm(mlr_vars_formula, data = train_p)
summ1 <- summary(mlr1)
cat("R2:", summ1$r.squared, " Adj R2:", summ1$adj.r.squared,
    " F-stat:", summ1$fstatistic[1], "\n")
## R2: 0.07325004  Adj R2: 0.06811015  F-stat: 14.2513
cat("RMSE:", sqrt(mean(residuals(mlr1)^2)), "\n")
## RMSE: 4375.669

Linear Model 2: log1p(TARGET_AMT), full variable set

The same variable set, but the response is log(1 + TARGET_AMT) to address the severe right skew documented in Data Exploration.

mlr2 <- lm(update(mlr_vars_formula, log1p(TARGET_AMT) ~ .), data = train_p)
summ2 <- summary(mlr2)
cat("R2:", summ2$r.squared, " Adj R2:", summ2$adj.r.squared,
    " F-stat:", summ2$fstatistic[1], "\n")
## R2: 0.2270173  Adj R2: 0.2227303  F-stat: 52.95396
pred_log  <- predict(mlr2)
pred_back <- expm1(pred_log)
cat("RMSE (log scale):", sqrt(mean(residuals(mlr2)^2)), "\n")
## RMSE (log scale): 3.223216
cat("RMSE (back-transformed to $):", sqrt(mean((train_p$TARGET_AMT - pred_back)^2)), "\n")
## RMSE (back-transformed to $): 4748.125
summary(mlr2)
## 
## Call:
## lm(formula = update(mlr_vars_formula, log1p(TARGET_AMT) ~ .), 
##     data = train_p)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.6982 -2.3347 -0.8817  2.0829 10.2750 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           4.145e-02  3.672e-01   0.113 0.910113    
## KIDSDRIV              5.876e-01  9.061e-02   6.485 9.50e-11 ***
## AGE                  -1.732e-03  5.620e-03  -0.308 0.757893    
## HOMEKIDS              5.291e-02  5.241e-02   1.010 0.312768    
## YOJ                  -2.062e-02  1.192e-02  -1.730 0.083666 .  
## INCOME               -3.474e-06  1.419e-06  -2.448 0.014377 *  
## HOME_VAL             -1.145e-06  4.683e-07  -2.444 0.014554 *  
## TRAVTIME              1.535e-02  2.556e-03   6.006 2.01e-09 ***
## BLUEBOOK             -1.777e-05  6.835e-06  -2.600 0.009345 ** 
## TIF                  -6.528e-02  9.690e-03  -6.736 1.76e-11 ***
## OLDCLAIM             -2.153e-05  5.789e-06  -3.719 0.000201 ***
## CLM_FREQ              2.704e-01  4.370e-02   6.187 6.49e-10 ***
## MVR_PTS               1.855e-01  2.066e-02   8.977  < 2e-16 ***
## CAR_AGE               2.906e-03  1.016e-02   0.286 0.774789    
## PARENT1Yes            6.535e-01  1.607e-01   4.066 4.83e-05 ***
## MSTATUSYes           -6.045e-01  1.156e-01  -5.231 1.74e-07 ***
## SEXM                  2.072e-01  1.458e-01   1.421 0.155439    
## RED_CARyes           -2.846e-02  1.184e-01  -0.240 0.810077    
## REVOKEDYes            1.309e+00  1.351e-01   9.685  < 2e-16 ***
## URBANICITYUrban       2.397e+00  1.103e-01  21.726  < 2e-16 ***
## CAR_USECommercial     1.027e+00  1.305e-01   7.872 4.07e-15 ***
## EDUCATIONBachelors   -6.207e-01  1.603e-01  -3.871 0.000109 ***
## EDUCATIONHigh School  5.686e-03  1.359e-01   0.042 0.966620    
## EDUCATIONMasters     -5.832e-01  2.119e-01  -2.753 0.005928 ** 
## EDUCATIONPhD         -5.071e-01  2.545e-01  -1.993 0.046348 *  
## JOBClerical           8.688e-02  1.500e-01   0.579 0.562458    
## JOBDoctor            -9.021e-01  3.118e-01  -2.893 0.003831 ** 
## JOBHome Maker         1.549e-02  2.073e-01   0.075 0.940440    
## JOBLawyer            -1.538e-01  2.036e-01  -0.756 0.449927    
## JOBManager           -9.565e-01  1.637e-01  -5.842 5.39e-09 ***
## JOBProfessional      -9.801e-02  1.570e-01  -0.624 0.532521    
## JOBStudent           -5.784e-02  1.842e-01  -0.314 0.753487    
## CAR_TYPEPanel Truck   2.265e-01  2.162e-01   1.047 0.294957    
## CAR_TYPEPickup        5.804e-01  1.346e-01   4.312 1.64e-05 ***
## CAR_TYPESports Car    1.189e+00  1.734e-01   6.862 7.44e-12 ***
## CAR_TYPESUV           9.691e-01  1.423e-01   6.810 1.07e-11 ***
## CAR_TYPEVan           4.900e-01  1.678e-01   2.920 0.003508 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.232 on 6491 degrees of freedom
## Multiple R-squared:  0.227,  Adjusted R-squared:  0.2227 
## F-statistic: 52.95 on 36 and 6491 DF,  p-value: < 2.2e-16

The log transform roughly triples R^2 relative to Model 1 and substantially improves the F-statistic, which is the expected payoff from addressing the skew.

Supplementary model: claim severity conditional on a crash

Models 1 and 2 both regress TARGET_AMT on the whole training set, where roughly three-quarters of records are structural zeros (no crash, no cost). That mixes two very different questions together: “does a crash happen?” and “if it happens, how expensive is it?” As a diagnostic, I also fit the severity question on its own, restricting to the 1722 records that actually crashed.

crash_only <- train_p %>% filter(TARGET_FLAG == 1)
mlr3 <- lm(update(mlr_vars_formula, log(TARGET_AMT) ~ .), data = crash_only)
summ3 <- summary(mlr3)
cat("n =", nrow(crash_only), " R2:", summ3$r.squared, " Adj R2:", summ3$adj.r.squared,
    " F-stat:", summ3$fstatistic[1], " p-value:",
    pf(summ3$fstatistic[1], summ3$fstatistic[2], summ3$fstatistic[3], lower.tail = FALSE), "\n")
## n = 1722  R2: 0.02694921  Adj R2: 0.006159998  F-stat: 1.296307  p-value: 0.1130199
summary(mlr3)
## 
## Call:
## lm(formula = update(mlr_vars_formula, log(TARGET_AMT) ~ .), data = crash_only)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9263 -0.3991  0.0424  0.4066  3.2587 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           7.955e+00  1.856e-01  42.870  < 2e-16 ***
## KIDSDRIV             -4.299e-02  3.627e-02  -1.185  0.23614    
## AGE                   1.373e-03  2.474e-03   0.555  0.57898    
## HOMEKIDS              1.075e-02  2.454e-02   0.438  0.66153    
## YOJ                  -3.729e-03  5.739e-03  -0.650  0.51596    
## INCOME               -8.007e-07  7.858e-07  -1.019  0.30838    
## HOME_VAL              1.812e-07  2.328e-07   0.778  0.43653    
## TRAVTIME             -3.331e-04  1.276e-03  -0.261  0.79409    
## BLUEBOOK              1.213e-05  3.510e-06   3.457  0.00056 ***
## TIF                   7.578e-04  4.871e-03   0.156  0.87639    
## OLDCLAIM              3.349e-06  2.536e-06   1.320  0.18687    
## CLM_FREQ             -2.453e-02  1.804e-02  -1.360  0.17401    
## MVR_PTS               1.565e-02  7.896e-03   1.982  0.04760 *  
## CAR_AGE              -4.423e-03  5.053e-03  -0.875  0.38154    
## PARENT1Yes            2.637e-02  6.796e-02   0.388  0.69800    
## MSTATUSYes           -8.124e-02  5.686e-02  -1.429  0.15326    
## SEXM                  1.087e-01  7.531e-02   1.444  0.14896    
## RED_CARyes            1.689e-02  5.724e-02   0.295  0.76798    
## REVOKEDYes           -8.433e-02  5.778e-02  -1.459  0.14464    
## URBANICITYUrban       6.170e-02  8.690e-02   0.710  0.47781    
## CAR_USECommercial    -5.328e-03  6.047e-02  -0.088  0.92980    
## EDUCATIONBachelors    4.567e-03  7.403e-02   0.062  0.95082    
## EDUCATIONHigh School  4.298e-02  5.935e-02   0.724  0.46909    
## EDUCATIONMasters      1.435e-01  1.030e-01   1.394  0.16360    
## EDUCATIONPhD          1.656e-01  1.302e-01   1.271  0.20381    
## JOBClerical          -4.507e-03  6.725e-02  -0.067  0.94657    
## JOBDoctor             6.890e-02  2.208e-01   0.312  0.75503    
## JOBHome Maker        -5.480e-02  9.949e-02  -0.551  0.58186    
## JOBLawyer             4.400e-03  1.056e-01   0.042  0.96675    
## JOBManager            2.554e-03  9.276e-02   0.028  0.97804    
## JOBProfessional       6.552e-02  7.496e-02   0.874  0.38220    
## JOBStudent           -1.047e-02  8.126e-02  -0.129  0.89747    
## CAR_TYPEPanel Truck   1.627e-02  1.100e-01   0.148  0.88247    
## CAR_TYPEPickup        3.966e-02  6.779e-02   0.585  0.55862    
## CAR_TYPESports Car    3.331e-02  8.687e-02   0.383  0.70140    
## CAR_TYPESUV           1.018e-01  7.627e-02   1.335  0.18207    
## CAR_TYPEVan          -4.338e-02  8.874e-02  -0.489  0.62499    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7909 on 1685 degrees of freedom
## Multiple R-squared:  0.02695,    Adjusted R-squared:  0.00616 
## F-statistic: 1.296 on 36 and 1685 DF,  p-value: 0.113

This is the most important substantive finding in the cost side of the assignment: conditional on a crash occurring, almost none of these variables explain how expensive the claim is. Adjusted R^2 is close to zero and the overall F-test is not significant. The variables that predict whether someone crashes (MVR_PTS, URBANICITY, CLM_FREQ, etc.) tell you almost nothing about how much the claim costs once a crash has happened; those appear to be two separate data-generating processes. I keep this model in the write-up specifically because it’s counter-intuitive and, to the boss’s likely question of “why doesn’t reckless driving predict expensive claims,” the honest answer is that the size of a claim is driven mostly by things this data set doesn’t observe (impact severity, what was hit, injury involvement), not by the driver characteristics that predict crash frequency.

Side-by-side comparison of all three linear models

stargazer(mlr1, mlr2, mlr3, type = "latex", header = FALSE,
          title = "Linear regression models: claim cost",
          column.labels = c("Model 1: raw \\$ (all)","Model 2: log1p (all)","Model 3: log severity (crash-only)"),
          dep.var.labels.include = FALSE,
          single.row = TRUE, no.space = TRUE, font.size = "tiny",
          omit.stat = c("bic"), star.cutoffs = c(0.05, 0.01, 0.001))

The stargazer table makes the earlier point visually obvious: almost every significance star in Model 1 and Model 2 (fit on the full, zero-inflated population) disappears in Model 3 (fit on crashes only). The same variables that light up when predicting whether a claim happens go dark when predicting how big it is.

Select Models

Logistic regression: comparing all five (plus the bonus penalized pair)

I evaluate each model on the training data at the assignment’s specified 0.5 classification threshold, using AIC/log-likelihood (in-sample fit, penalized for complexity, where applicable) plus the full confusion-matrix battery. LASSO and Ridge don’t have a comparable AIC/log-likelihood (they’re fit by cross-validated penalty, not maximum likelihood), so those cells are left blank rather than filled with a misleading number.

score_logit <- function(model, data, y) {
  p <- predict(model, newdata = data, type = "response")
  pred <- factor(ifelse(p >= 0.5, 1, 0), levels = c(0,1))
  cm <- confusionMatrix(pred, factor(y, levels = c(0,1)), positive = "1")
  roc_obj <- roc(y, p, quiet = TRUE)
  data.frame(
    AIC = AIC(model),
    LogLik = as.numeric(logLik(model)),
    Accuracy = cm$overall["Accuracy"],
    Sensitivity = cm$byClass["Sensitivity"],
    Specificity = cm$byClass["Specificity"],
    Precision = cm$byClass["Precision"],
    F1 = cm$byClass["F1"],
    AUC = as.numeric(auc(roc_obj))
  )
}

score_penalized <- function(cvfit, x, y) {
  p <- as.numeric(predict(cvfit, newx = x, s = "lambda.1se", type = "response"))
  pred <- factor(ifelse(p >= 0.5, 1, 0), levels = c(0,1))
  cm <- confusionMatrix(pred, factor(y, levels = c(0,1)), positive = "1")
  roc_obj <- roc(y, p, quiet = TRUE)
  data.frame(
    AIC = NA_real_, LogLik = NA_real_,
    Accuracy = cm$overall["Accuracy"],
    Sensitivity = cm$byClass["Sensitivity"],
    Specificity = cm$byClass["Specificity"],
    Precision = cm$byClass["Precision"],
    F1 = cm$byClass["F1"],
    AUC = as.numeric(auc(roc_obj))
  )
}

y_train <- train_p$TARGET_FLAG
res1 <- score_logit(logit1, train_p, y_train)
res2 <- score_logit(logit2, train_p, y_train)
res3 <- score_logit(logit3, train_p, y_train)
res4 <- score_logit(probit1, train_p, y_train)
res5 <- score_penalized(cv_lasso, x_train, y_train)
res6 <- score_penalized(cv_ridge, x_train, y_train)

comp <- rbind(res1, res2, res3, res4, res5, res6)
rownames(comp) <- c("Model 1 (Full)", "Model 2 (Stepwise)", "Model 3 (Parsimonious)",
                     "Bonus (Probit)", "Bonus (LASSO)", "Bonus (Ridge)")
# transpose: metrics as rows, models as columns, so no row label has to
# share a line with a number the way "Model3_Parsimonious" and its AIC did
comp_t <- as.data.frame(round(t(comp), 4))
kable(comp_t, row.names = TRUE,
      caption = "Logistic model comparison on training data, 0.5 threshold (rows = metrics, columns = models). Error rate = 1 - Accuracy, omitted as redundant.")
Logistic model comparison on training data, 0.5 threshold (rows = metrics, columns = models). Error rate = 1 - Accuracy, omitted as redundant.
Model 1 (Full) Model 2 (Stepwise) Model 3 (Parsimonious) Bonus (Probit) Bonus (LASSO) Bonus (Ridge)
AIC 5905.3518 5892.6612 6147.5915 5897.8581 NA NA
LogLik -2910.6759 -2914.3306 -3060.7957 -2916.9290 NA NA
Accuracy 0.7930 0.7935 0.7774 0.7924 0.7851 0.7854
Sensitivity 0.4262 0.4245 0.3374 0.4181 0.3188 0.3171
Specificity 0.9245 0.9257 0.9351 0.9266 0.9521 0.9532
Precision 0.6691 0.6719 0.6506 0.6710 0.7047 0.7082
F1 0.5208 0.5203 0.4444 0.5152 0.4390 0.4380
AUC 0.8146 0.8140 0.7873 0.8139 0.8069 0.8102 
roc1 <- roc(y_train, predict(logit1, type = "response"), quiet = TRUE)
roc2 <- roc(y_train, predict(logit2, type = "response"), quiet = TRUE)
roc3 <- roc(y_train, predict(logit3, type = "response"), quiet = TRUE)
roc4 <- roc(y_train, predict(probit1, type = "response"), quiet = TRUE)
roc5 <- roc(y_train, as.numeric(predict(cv_lasso, newx = x_train, s = "lambda.1se", type = "response")), quiet = TRUE)
roc6 <- roc(y_train, as.numeric(predict(cv_ridge, newx = x_train, s = "lambda.1se", type = "response")), quiet = TRUE)

plot(roc1, col = "steelblue", main = "ROC: all six logistic specifications")
lines(roc2, col = "darkorange")
lines(roc3, col = "seagreen")
lines(roc4, col = "purple")
lines(roc5, col = "black")
lines(roc6, col = "grey50")
legend("bottomright", cex = 0.7,
       legend = c("Model 1 (full)","Model 2 (stepwise)","Model 3 (parsimonious)",
                  "Bonus: probit","Bonus: LASSO","Bonus: Ridge"),
       col = c("steelblue","darkorange","seagreen","purple","black","grey50"), lwd = 2)
ROC curves for all six models

ROC curves for all six models

The probit model lands within a rounding error of Model 2’s AUC and accuracy, which makes sense since it uses the identical variable set, just a different link function. LASSO and Ridge are the more interesting case, and worth addressing directly rather than glossing over: penalized regression is usually framed as theoretically superior to plain stepwise selection, precisely because it’s supposed to guard against overfitting. So why does it come in a shade behind here (AUC 0.807-0.810 versus Model 2’s 0.814)?

Penalization earns its keep when there’s overfitting risk for it to correct: many predictors relative to sample size, meaningful multicollinearity, or an unconstrained variable set prone to fitting noise. None of those conditions hold much here. The training sample is large relative to the variable count (n = 6,528 against roughly 28-40 terms), the VIF checks on both Model 1 and Model 2 earlier came back comfortably under 5, and Model 2’s variables were already whittled down by a principled AIC-based procedure rather than kept in an unconstrained kitchen-sink form. With little instability left to fix, the LASSO/Ridge penalty still does what it’s designed to do (shrink coefficients toward zero to trade a little bias for lower variance), but there’s not much variance left to trade away, so the shrinkage mostly just adds bias, showing up as marginally lower in-sample fit. There’s also a deliberate design choice behind that gap: I tuned both penalized models to lambda.1se rather than lambda.min, which is the standard convention for picking the most parsimonious model within one cross-validated standard error of the best one, explicitly trading a sliver of fit for a simpler, more stable model, rather than the tuning that would most favor LASSO/Ridge on raw performance. None of the three bonus models beats Model 2 outright, and given the reasoning above that’s the expected result on this data, not a red flag: it’s good supporting evidence that Model 2 was already well-specified and not overfit, rather than one a fancier estimator could meaningfully improve on.

Model selection: Model 2 (backward-stepwise). Model 2 actually has the lowest AIC of the three despite dropping several variables from Model 1, a cleaner example of parsimony paying off than I expected going in. Its accuracy, sensitivity, specificity, and AUC are statistically indistinguishable from the full kitchen-sink model, so there’s no accuracy cost to dropping the variables stepAIC removed. Model 3 is the most interpretable and by far the smallest, but it gives up close to 3 points of AUC relative to Models 1 and 2, and the trade for parsimony there is real, not free. Given the assignment explicitly asks whether I’d trade a bit of performance for parsimony, my answer is: not this time, because Model 2’s extra variables aren’t just noise. The job category and car type variables it retains are individually significant and carry sensible signs (see interpretation below).

Coefficient interpretation, Model 2 (selected)

or_table <- exp(cbind(OddsRatio = coef(logit2), confint.default(logit2)))
kable(round(or_table, 3), caption = "Odds ratios, selected logistic model (Model 2)")
Odds ratios, selected logistic model (Model 2)
OddsRatio 2.5 % 97.5 %
(Intercept) 0.036 0.024 0.056
KIDSDRIV 1.600 1.416 1.807
INCOME 1.000 1.000 1.000
HOME_VAL 1.000 1.000 1.000
TRAVTIME 1.013 1.009 1.018
BLUEBOOK 1.000 1.000 1.000
TIF 0.945 0.930 0.960
OLDCLAIM 1.000 1.000 1.000
CLM_FREQ 1.228 1.154 1.307
MVR_PTS 1.124 1.091 1.158
PARENT1Yes 1.597 1.297 1.965
MSTATUSYes 0.618 0.519 0.736
REVOKEDYes 2.572 2.114 3.129
URBANICITYUrban 10.496 8.214 13.412
CAR_USECommercial 2.118 1.730 2.593
EDUCATIONBachelors 0.646 0.509 0.820
EDUCATIONHigh School 0.992 0.807 1.221
EDUCATIONMasters 0.769 0.541 1.093
EDUCATIONPhD 0.873 0.566 1.346
JOBClerical 1.031 0.816 1.302
JOBDoctor 0.305 0.155 0.598
JOBHome Maker 0.921 0.670 1.266
JOBLawyer 0.734 0.489 1.101
JOBManager 0.405 0.297 0.553
JOBProfessional 0.815 0.625 1.061
JOBStudent 0.944 0.721 1.236
CAR_TYPEPanel Truck 1.714 1.227 2.395
CAR_TYPEPickup 1.732 1.392 2.156
CAR_TYPESports Car 2.514 1.985 3.185
CAR_TYPESUV 2.073 1.719 2.500
CAR_TYPEVan 1.828 1.397 2.392
M_JOB 0.696 0.463 1.047

A note on reading this table: INCOME and HOME_VAL both show an odds ratio of 1.000. That is not a rounding error and it does not mean these variables have no effect; it means their odds ratio per single dollar is indistinguishable from 1 at 3 decimal places, because a single dollar is a tiny move relative to the scale these variables live on (HOME_VAL ranges up into the hundreds of thousands). The table below rescales the same two coefficients to a more meaningful $10,000 increment.

rescale_tbl <- data.frame(
  Variable = c("INCOME", "HOME_VAL"),
  OddsRatio_per_dollar   = round(exp(coef(logit2)[c("INCOME","HOME_VAL")]), 6),
  OddsRatio_per_10k      = round(exp(coef(logit2)[c("INCOME","HOME_VAL")] * 10000), 4),
  Pct_change_per_10k     = round((1 - exp(coef(logit2)[c("INCOME","HOME_VAL")] * 10000)) * 100, 2)
)
kable(rescale_tbl, row.names = FALSE,
      caption = "INCOME and HOME_VAL odds ratios rescaled to a $10,000 increment")
INCOME and HOME_VAL odds ratios rescaled to a $10,000 increment
Variable OddsRatio_per_dollar OddsRatio_per_10k Pct_change_per_10k
INCOME 0.999997 0.9680 3.20
HOME_VAL 0.999999 0.9887 1.13

Two examples, one positive and one negative coefficient, as requested, each touching on the same four things: what the raw coefficient means, its practical/economic magnitude, its statistical significance, and whether the sign matches the assignment’s theoretical prior:

  • MVR_PTS (positive, log-odds coefficient around 0.117): each additional motor-vehicle-record point multiplies the odds of a crash by about 1.124, i.e. roughly a 12.4% increase in the odds of a crash per point, economically a large, intuitive effect: someone with 5 points on their record has odds of crashing roughly 1.79 times higher than an otherwise-identical driver with zero points. It’s also highly statistically significant, p = <0.001, and the positive sign matches the assignment sheet’s prior exactly (“the more tickets you get, the more crashes you tend to have”).
  • HOME_VAL (negative, coefficient near zero on the dollar scale because home value is measured in single dollars): a $10,000 increase in home value is associated with roughly a 1.1% decrease in the odds of a crash, a modest but real economic effect that compounds over larger swings in home value. It’s statistically significant at p = 0.0029, and the negative sign is exactly what the “homeowners drive more responsibly” prior in the assignment sheet predicts.

Odds ratios are the standard way to report a logistic coefficient, but they don’t always convey practical size on their own, so here’s the same model translated into predicted probabilities for two illustrative policyholders who differ only on the handful of variables below:

mk_profile <- function(...) {
  base <- data.frame(
    KIDSDRIV = 0, AGE = 45, HOMEKIDS = 0, YOJ = 11, INCOME = 55000,
    HOME_VAL = 150000, TRAVTIME = 33, BLUEBOOK = 15000, TIF = 6,
    OLDCLAIM = 0, CLM_FREQ = 0, MVR_PTS = 0, CAR_AGE = 8,
    PARENT1 = factor("No", levels = levels(train_p$PARENT1)),
    MSTATUS = factor("Yes", levels = levels(train_p$MSTATUS)),
    SEX = factor("F", levels = levels(train_p$SEX)),
    RED_CAR = factor("no", levels = levels(train_p$RED_CAR)),
    REVOKED = factor("No", levels = levels(train_p$REVOKED)),
    URBANICITY = factor("Urban", levels = levels(train_p$URBANICITY)),
    CAR_USE = factor("Private", levels = levels(train_p$CAR_USE)),
    EDUCATION = factor("Bachelors", levels = levels(train_p$EDUCATION)),
    JOB = factor("Clerical", levels = levels(train_p$JOB)),
    CAR_TYPE = factor("Minivan", levels = levels(train_p$CAR_TYPE)),
    M_YOJ = 0, M_INCOME = 0, M_HOME_VAL = 0, M_CAR_AGE = 0, M_JOB = 0
  )
  modifyList(base, list(...)) %>% as.data.frame()
}

low_risk  <- mk_profile()
high_risk <- mk_profile(MVR_PTS = 6, CLM_FREQ = 3, OLDCLAIM = 8000,
                         REVOKED = factor("Yes", levels = levels(train_p$REVOKED)),
                         CAR_USE = factor("Commercial", levels = levels(train_p$CAR_USE)),
                         MSTATUS = factor("No", levels = levels(train_p$MSTATUS)),
                         PARENT1 = factor("Yes", levels = levels(train_p$PARENT1)),
                         KIDSDRIV = 1, TIF = 1)

profiles <- rbind(low_risk, high_risk)
profiles$P_crash <- round(predict(logit2, newdata = profiles, type = "response"), 3)
kable(data.frame(Profile = c("Married homeowner, clean record, private use",
                              "Single parent, 6 MVR points, revoked license, commercial use"),
                  P_crash = profiles$P_crash),
      caption = "Predicted crash probability: two contrasting policyholder profiles")
Predicted crash probability: two contrasting policyholder profiles
Profile P_crash
Married homeowner, clean record, private use 0.078
Single parent, 6 MVR points, revoked license, commercial use 0.893

Holding everything else at the same baseline values, moving from the low-risk profile to the high-risk one (a handful of realistic changes, not an extreme hypothetical) moves the predicted crash probability from 8% to 89%, roughly an 11.4x increase in probability. That’s the practical, real-world size of what the odds ratios above represent: this isn’t a model where the “significant” variables move risk by a percentage point or two. A realistic combination of a few of them is the difference between a policyholder who’s very unlikely to file a claim and one who’s more likely than not to.

Diagnostics on the selected model

vif2 <- vif(logit2)
# vif() returns a matrix (with GVIF^(1/(2*Df))) for models containing factors
if (is.matrix(vif2)) {
  vif2_tbl <- data.frame(Variable = rownames(vif2), GVIF_adj = round(vif2[, 3], 2))
} else {
  vif2_tbl <- data.frame(Variable = names(vif2), GVIF_adj = round(vif2, 2))
}
kable(vif2_tbl, row.names = FALSE,
      caption = "Variance inflation (Model 2, selected logistic model)")
Variance inflation (Model 2, selected logistic model)
Variable GVIF_adj
KIDSDRIV 1.04
INCOME 1.57
HOME_VAL 1.36
TRAVTIME 1.02
BLUEBOOK 1.32
TIF 1.01
OLDCLAIM 1.28
CLM_FREQ 1.21
MVR_PTS 1.08
PARENT1 1.20
MSTATUS 1.37
REVOKED 1.14
URBANICITY 1.07
CAR_USE 1.57
EDUCATION 1.28
JOB 1.20
CAR_TYPE 1.10
M_JOB 1.66 
cat("Max adjusted GVIF:", round(max(vif2_tbl$GVIF_adj), 2), "\n")
## Max adjusted GVIF: 1.66

Even in the larger, 26-plus-variable stepwise model, the worst adjusted VIF stays well under the usual rule-of-thumb cutoff of 5, so multicollinearity isn’t distorting these coefficient estimates despite the number of dummy variables in play.

Linear regression: comparing the two required models plus the severity check

mlr_metrics <- function(model, label) {
  s <- summary(model)
  data.frame(Model = label, R2 = s$r.squared, AdjR2 = s$adj.r.squared,
             F = unname(s$fstatistic[1]), RMSE = sqrt(mean(residuals(model)^2)))
}
mlr_comp <- rbind(
  mlr_metrics(mlr1, "Model 1: raw $ (full data)"),
  mlr_metrics(mlr2, "Model 2: log1p (full data)"),
  mlr_metrics(mlr3, "Model 3: log severity (crash-only)")
)
kable(mlr_comp, digits = 4, row.names = FALSE, caption = "Linear model comparison")
Linear model comparison
Model R2 AdjR2 F RMSE
Model 1: raw $ (full data) 0.0733 0.0681 14.2513 4375.6688
Model 2: log1p (full data) 0.2270 0.2227 52.9540 3.2232
Model 3: log severity (crash-only) 0.0269 0.0062 1.2963 0.7824 
vif_mlr2 <- vif(mlr2)
if (is.matrix(vif_mlr2)) {
  vif_mlr2_tbl <- data.frame(Variable = rownames(vif_mlr2), GVIF_adj = round(vif_mlr2[, 3], 2))
} else {
  vif_mlr2_tbl <- data.frame(Variable = names(vif_mlr2), GVIF_adj = round(vif_mlr2, 2))
}
kable(vif_mlr2_tbl, row.names = FALSE,
      caption = "Variance inflation (Model 2, selected linear model)")
Variance inflation (Model 2, selected linear model)
Variable GVIF_adj
KIDSDRIV 1.15
AGE 1.21
HOMEKIDS 1.45
YOJ 1.19
INCOME 1.65
HOME_VAL 1.47
TRAVTIME 1.02
BLUEBOOK 1.43
TIF 1.01
OLDCLAIM 1.29
CLM_FREQ 1.27
MVR_PTS 1.11
CAR_AGE 1.41
PARENT1 1.36
MSTATUS 1.41
SEX 1.82
RED_CAR 1.35
REVOKED 1.13
URBANICITY 1.12
CAR_USE 1.57
EDUCATION 1.28
JOB 1.19
CAR_TYPE 1.18 
cat("Max adjusted GVIF:", round(max(vif_mlr2_tbl$GVIF_adj), 2), "\n")
## Max adjusted GVIF: 1.82

All VIFs for the selected linear model stay comfortably under 5 as well.

par(mfrow = c(2,2))
plot(mlr2)
Diagnostic plots for the selected linear model

Diagnostic plots for the selected linear model 

par(mfrow = c(1,1))

Model selection: Model 2 (log1p, full data). R^2 nearly triples relative to the raw-dollar Model 1, and the F-statistic is far more significant, so on the assignment’s own suggested criteria (R^2, RMSE, F-statistic) Model 2 wins outright. The residual diagnostics aren’t clean, though, and I want to be upfront about that rather than bury it: the Scale-Location plot still shows a funnel shape and the Q-Q plot has a heavy tail, both echoes of the fact that roughly three-quarters of the outcome variable is exactly zero. A log transform narrows that problem, it doesn’t eliminate it, because no monotonic transform can make a zero-inflated variable look normal.

That’s exactly why I keep Model 3 (the crash-only severity regression) as a companion piece rather than a competing candidate for “the” cost model: it isolates the “how much, given a crash” question from the “zero vs. positive” question, and its near-zero R^2 is itself the finding worth reporting: claim severity is close to unpredictable from these particular policyholder characteristics, even though claim frequency (Model 2 of the logistic section) is quite predictable from the very same variables. For generating point predictions on the evaluation set below, I use the two-part combination that’s standard for this kind of zero-inflated outcome: predicted cost = P(crash) x E[cost | crash], using the selected logistic model for the first term and Model 2 for the second, rather than either the raw full-population OLS or the crash-only severity model alone.

Predictions on the evaluation set

Using the selected models (logistic Model 2 for probability/classification, and linear Model 2 for the dollar amount, applied only to records classified as a crash), I score all 1633 evaluation records at the assignment’s required 0.5 threshold.

p_eval    <- predict(logit2, newdata = test_p, type = "response")
flag_eval <- ifelse(p_eval >= 0.5, 1, 0)

amt_if_crash <- expm1(predict(mlr2, newdata = test_p))
amt_if_crash <- pmax(amt_if_crash, 0)
amt_eval     <- ifelse(flag_eval == 1, amt_if_crash, 0)

predictions <- data.frame(
  INDEX            = test_p$INDEX,
  P_TARGET_FLAG    = round(p_eval, 4),
  TARGET_FLAG_PRED = flag_eval,
  TARGET_AMT_PRED  = round(amt_eval, 2)
)
write.csv(predictions, "evaluation_predictions.csv", row.names = FALSE)
kable(head(predictions, 10), caption = "First 10 evaluation-set predictions")
First 10 evaluation-set predictions
INDEX P_TARGET_FLAG TARGET_FLAG_PRED TARGET_AMT_PRED
5 0.0904 0 0.00
8 0.3121 0 0.00
26 0.4720 0 0.00
40 0.1913 0 0.00
45 0.0566 0 0.00
55 0.6258 1 108.36
61 0.9079 1 1265.51
66 0.0184 0 0.00
67 0.2879 0 0.00
71 0.3191 0 0.00 
cat("Predicted crash rate on evaluation set:", mean(flag_eval), "\n")
## Predicted crash rate on evaluation set: 0.1702388
cat("Actual crash rate on training set:     ", mean(train_p$TARGET_FLAG), "\n")
## Actual crash rate on training set:      0.2637868

The predicted crash rate on the evaluation set lines up closely with the training rate, which is a basic sanity check that the model isn’t systematically over- or under-calling crashes on new data.

A bonus honesty check: this evaluation file happens to include the true outcomes

Most versions of this assignment ship an evaluation file with blank targets. This copy of insurance-testing-data2.csv happens to still carry the true TARGET_FLAG and TARGET_AMT values, so, purely as an extra validation step and not part of the required deliverable, I use them to check out-of-sample performance:

true_flag <- test_raw$TARGET_FLAG
true_amt  <- test_raw$TARGET_AMT

holdout_auc <- as.numeric(auc(roc(true_flag, p_eval, quiet = TRUE)))
holdout_acc <- mean(flag_eval == true_flag)
cat("Holdout AUC:", round(holdout_auc, 4), " Holdout accuracy:", round(holdout_acc, 4), "\n")
## Holdout AUC: 0.8082  Holdout accuracy: 0.7936
holdout_rmse <- sqrt(mean((amt_eval - true_amt)^2))
cat("Holdout $ RMSE:", round(holdout_rmse, 2), "\n")
## Holdout $ RMSE: 5508.76

Out-of-sample AUC and accuracy come in essentially identical to the training-set numbers reported above, which is reassuring evidence that Model 2 isn’t overfit to the training data.

Conclusion

The two questions this assignment asks, will someone crash and how much will it cost, turn out to have very different answers in terms of how predictable they are from the variables provided.

On the data side, the exploration and preparation work mattered as much as the modeling: fixing the currency-formatted fields, imputing the roughly 5-6% of records missing AGE, YOJ, INCOME, HOME_VAL, JOB, and CAR_AGE while flagging that missingness explicitly, and log-transforming the heavily skewed dollar variables were all necessary before any model coefficient could be trusted. The categorical breakdown was just as informative as the numeric one: urbanicity, single-parent status, and job category all showed real spread in raw crash rate well before any model touched them.

On the modeling side, crash probability (TARGET_FLAG) is well-explained by past behavior (MVR_PTS, CLM_FREQ, prior REVOKED status), vehicle use, and, unexpectedly to me, urban/rural location, which alone carries close to a 10x odds ratio. I selected the backward-stepwise logistic model (Model 2, AUC 0.814) over the full kitchen-sink and parsimonious alternatives because it matched the full model’s performance with fewer variables, and the bonus LASSO/Ridge/probit check confirmed it wasn’t leaving performance on the table: none of the three beat it, for reasons the penalized-regression discussion above lays out explicitly rather than treating as a mystery. The practical size of that model is large in absolute terms too: a realistic shift from a low-risk to a high-risk profile moves predicted crash probability from single digits to roughly ninety percent.

Claim severity, once a crash has happened, is a different story: it’s close to unexplainable by anything in this data set (the crash-only severity model’s R^2 is close to zero and its F-test isn’t significant), even though the same variables predict whether a crash happens quite well. The honest conclusion, and one I’d flag to a manager rather than paper over, is that better cost prediction would need information this file doesn’t contain (impact type, injury involvement, vehicle repair complexity) rather than more clever modeling of the variables already here. For the evaluation-set predictions, I used the two-part combination this implies: the selected logistic model for P(crash), multiplied through by the log1p linear model’s severity estimate, which the holdout validation check confirmed generalizes about as well as the training fit suggested it would.