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.

train_path <- "C:/Users/dbely/OneDrive/Desktop/Econometrics/Assignment 2/insurance-training-data2-2.csv"
test_path  <- "C:/Users/dbely/OneDrive/Desktop/Econometrics/Assignment 2/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_histogram(bins = 40, fill = "steelblue", color = "white") +
  scale_x_continuous(labels = scales::comma) +
  labs(title = "Claim amount, conditional on a crash", x = "TARGET_AMT ($)", y = "Count") +
  theme_minimal()

Distribution of claim amount among crashes

ggplot(train_raw %>% filter(TARGET_AMT > 0), aes(x = log(TARGET_AMT))) +
  geom_histogram(bins = 40, fill = "darkorange", color = "white") +
  labs(title = "log(TARGET_AMT), conditional on a crash", x = "log(TARGET_AMT)", y = "Count") +
  theme_minimal()

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

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

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.7, tl.col = "black",
         addCoef.col = "black", number.cex = 0.55, diag = FALSE,
         title = "Correlation matrix: targets and numeric predictors", mar = c(0,0,2,0))

Correlation heatmap of numeric variables

corr_tbl <- data.frame(
  Variable = rownames(corr_all),
  With_TARGET_FLAG = round(corr_all[, "TARGET_FLAG"], 3),
  With_TARGET_AMT  = round(corr_all[, "TARGET_AMT"], 3)
) %>% filter(!Variable %in% c("TARGET_FLAG","TARGET_AMT")) %>%
  arrange(With_TARGET_FLAG)
rownames(corr_tbl) <- NULL
kable(corr_tbl, caption = "Correlation of numeric predictors with each target variable")

Correlation of numeric predictors with each target variable
Variable	With_TARGET_FLAG	With_TARGET_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. The same variables correlate with TARGET_AMT in the same direction but much more weakly, 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 below now shows the full 6528 observations (no more NAs pulling n down), and the log-transformed dollar variables have 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("n","mean","sd","median","min","max","skew")]
kable(round(post_desc, 1), caption = "Post-cleaning summary statistics (imputed + transformed)")

Post-cleaning summary statistics (imputed + transformed)
	n	mean	sd	median	min	max	skew
AGE	6528	44.9	8.6	45.0	16.0	81.0	0.0
YOJ	6528	10.5	4.0	11.0	0.0	23.0	-1.3
INCOME	6528	61195.3	46386.8	53276.0	0.0	367030.0	1.3
LOG_INCOME	6528	10.0	3.0	10.9	0.0	12.8	-2.8
HOME_VAL	6528	154658.2	125197.5	160680.0	0.0	885282.0	0.5
LOG_HOME_VAL	6528	8.8	5.5	12.0	0.0	13.7	-1.0
BLUEBOOK	6528	15641.8	8381.5	14370.0	1500.0	69740.0	0.8
LOG_BLUEBOOK	6528	9.5	0.6	9.6	7.3	11.2	-0.9
OLDCLAIM	6528	4119.3	8924.7	0.0	0.0	57037.0	3.1
LOG_OLDCLAIM	6528	3.4	4.3	0.0	0.0	11.0	0.5
CAR_AGE	6528	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.

Logistic Model 1: full “kitchen sink”

Every candidate predictor, untransformed, so I can see the whole field before narrowing anything down.

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.

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.

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)
  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,
      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	MVR_PTS	0.123	0.117	0.116
CLM_FREQ	CLM_FREQ	0.169	0.205	0.205
URBANICITYUrban	URBANICITY (Urban)	2.060	2.351	2.349
CAR_USECommercial	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).

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"],
    ErrorRate = 1 - 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"],
    ErrorRate = 1 - 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(Model1_Full = res1, Model2_Stepwise = res2, Model3_Parsimonious = res3,
              Bonus_Probit = res4, Bonus_LASSO = res5, Bonus_Ridge = res6)
kable(comp, digits = 4, row.names = TRUE,
      caption = "Logistic model comparison on training data, 0.5 threshold")

Logistic model comparison on training data, 0.5 threshold
	AIC	LogLik	Accuracy	ErrorRate	Sensitivity	Specificity	Precision	F1	AUC
Model1_Full	5905.352	-2910.676	0.7930	0.2070	0.4262	0.9245	0.6691	0.5208	0.8146
Model2_Stepwise	5892.661	-2914.331	0.7935	0.2065	0.4245	0.9257	0.6719	0.5203	0.8140
Model3_Parsimonious	6147.592	-3060.796	0.7774	0.2226	0.3374	0.9351	0.6506	0.4444	0.7873
Bonus_Probit	5897.858	-2916.929	0.7924	0.2076	0.4181	0.9266	0.6710	0.5152	0.8139
Bonus_LASSO	NA	NA	0.7851	0.2149	0.3188	0.9521	0.7047	0.4390	0.8069
Bonus_Ridge	NA	NA	0.7854	0.2146	0.3171	0.9532	0.7082	0.4380	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

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. LASSO and Ridge trade a little AUC (roughly 0.807-0.810 versus Model 2’s 0.814) for higher specificity, consistent with a penalty that shrinks marginal coefficients toward zero and makes the model slightly more conservative about flagging a crash. None of the three bonus models beats Model 2 outright, which is itself useful evidence that Model 2 was already a well-specified, non-overfit model rather than one a fancier estimator could visibly 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

Two examples, one positive and one negative coefficient, as requested:

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. Someone with 5 points on their record has odds of crashing roughly 1.79 times higher than an otherwise-identical driver with zero points, a large, intuitive effect.
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, consistent with the “homeowners drive more responsibly” prior in the assignment sheet.

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

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. Crash probability 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. Claim severity, once a crash has happened, is close to unexplainable by anything in this data set. 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.

ADEC 7320 - Homework #2: Modeling Auto Insurance Crash Risk and Cost

David N. Belyn

July 26, 2026