1 Executive summary

This report asks whether, for credit card customers with the same latest repayment status, rising bill balances and declining payments help predict default in the following month. Recent bill and payment trends are built from the six months of history in the dataset and added to a model that already contains repayment status and customer characteristics. The model is fitted on a training sample and judged on a separate test sample that was not used for fitting or selection.

All figures, tables and conclusions below are calculated from the data when the document is knitted; none are typed in by hand. The headline numbers (coefficients, R-squared, RMSE, AUC and the verdict on whether the trends improve prediction) are given in Section 5, and the plain-language answer to the business question is in Section 6.

2 Problem statement

2.1 Project overview and business question

Business question. Among credit card customers with the same latest repayment status, do rising bill balances and declining payments help predict default in the following month?

Purpose. We assess whether recent bill and payment trends provide additional predictive information beyond repayment status and other customer characteristics.

Nature of the study. The analysis is predictive and observational. Nothing in this document claims that a change in bill balances or payments causes default, or that an intervention would prevent it. The data are historical and come from Taiwan, so results may not generalise to present-day customers or to other markets.

Modelling framework.

  • The outcome is binary (default = 1, no default = 0).
  • The main model is a linear probability model (LPM) fitted by OLS with lm(). It is used because the course emphasises OLS and OLS diagnostics.

2.1.1 Package check

required_pkgs <- c("readxl", "ggplot2", "car", "lmtest", "sandwich", "knitr", "scales")
missing_pkgs <- required_pkgs[!vapply(required_pkgs, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_pkgs) > 0) {
  stop("Please install the missing packages once, then knit again:\n  install.packages(c(",
       paste0('"', missing_pkgs, '"', collapse = ", "), "))")
}
library(ggplot2)

# Consistent look for all charts
theme_set(theme_minimal(base_size = 12) +
            theme(plot.title = element_text(face = "bold"),
                  panel.grid.minor = element_blank()))
col_main <- "#2C7FB8"
col_def  <- c("No default" = "#2C7FB8", "Default" = "#D95F0E")
# ---- Formatting helpers (used for dynamic text and tables) ----
fmt_num <- function(x, d = 3) formatC(x, format = "f", digits = d, big.mark = ",")
fmt_pct <- function(x, d = 1) paste0(formatC(100 * x, format = "f", digits = d), "%")
fmt_p   <- function(p) ifelse(p < 0.001, "< 0.001", formatC(p, format = "f", digits = 3))
list_and <- function(x) {
  x <- x[!is.na(x) & nzchar(x)]
  if (length(x) == 0) return("no predictors listed")
  if (length(x) == 1) return(x)
  if (length(x) == 2) return(paste(x, collapse = " and "))
  paste0(paste(x[-length(x)], collapse = ", "), " and ", x[length(x)])
}

# ---- Evaluation helpers ----
auc_rank <- function(y, p) {            # AUC via the rank (Mann-Whitney) identity
  r <- rank(p); n1 <- sum(y == 1); n0 <- sum(y == 0)
  (sum(r[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}
log_loss_vec <- function(y, p, eps = 1e-15) {
  p <- pmin(pmax(p, eps), 1 - eps)
  -(y * log(p) + (1 - y) * log(1 - p))
}
predict_prob <- function(model, newdata) {
  if (inherits(model, "glm")) as.numeric(predict(model, newdata = newdata, type = "response"))
  else as.numeric(predict(model, newdata = newdata))
}

3 Data description

3.1 Dataset import and initial validation

The importer reads every sheet as text, finds the row that contains the variable name LIMIT_BAL, and treats that row as the header. If the header sits inside the data (the X1 ... X24 layout), the rows above it are discarded and every row below it is kept. If the file has proper headers, the header is simply row 1.

find_data_file <- function(user_path) {
  if (!is.null(user_path) && nzchar(user_path)) {
    if (file.exists(user_path)) return(normalizePath(user_path))
    warning("The file given in params$data_file was not found (", user_path,
            "). Falling back to automatic detection of an Excel file.")
  }
  dirs <- c(".", "data", "../data")
  cand <- as.character(unlist(lapply(dirs, function(d) {
    if (dir.exists(d)) list.files(d, pattern = "\\.xlsx?$", full.names = TRUE, ignore.case = TRUE)
  })))
  cand <- cand[!grepl("^~\\$", basename(cand))]            # ignore Excel lock files
  cand <- unique(normalizePath(cand))
  if (length(cand) == 0) {
    stop("No Excel file found. Fix params$data_file in the YAML header, or put the dataset ",
         "next to this .Rmd or in a 'data' folder.")
  }
  pref <- cand[grepl("default|credit", basename(cand), ignore.case = TRUE)]
  if (length(pref) >= 1) cand <- pref
  if (length(cand) > 1) {
    stop("Several candidate Excel files were found:\n  ", paste(cand, collapse = "\n  "),
         "\nSet params$data_file to the correct one.")
  }
  cand
}

norm_names <- function(x) {                      # "default payment next month" -> "DEFAULT_PAYMENT_NEXT_MONTH"
  x <- toupper(trimws(as.character(x)))
  x <- gsub("[^A-Z0-9]+", "_", x)
  gsub("^_+|_+$", "", x)
}

locate_header <- function(path) {
  for (sh in readxl::excel_sheets(path)) {
    top <- suppressMessages(readxl::read_excel(path, sheet = sh, col_names = FALSE,
                                               col_types = "text", n_max = 10))
    top <- as.data.frame(top, stringsAsFactors = FALSE)
    if (nrow(top) == 0) next
    for (i in seq_len(nrow(top))) {
      if ("LIMIT_BAL" %in% norm_names(unlist(top[i, ], use.names = FALSE))) {
        return(list(sheet = sh, header_row = i))
      }
    }
  }
  stop("Could not find a header row containing LIMIT_BAL in the first 10 rows of any sheet.")
}

data_path <- find_data_file(params$data_file)
hdr       <- locate_header(data_path)

raw_all <- as.data.frame(
  suppressMessages(readxl::read_excel(data_path, sheet = hdr$sheet, col_names = FALSE, col_types = "text")),
  stringsAsFactors = FALSE
)
col_names <- norm_names(unlist(raw_all[hdr$header_row, ], use.names = FALSE))
dat_text  <- raw_all[-seq_len(hdr$header_row), , drop = FALSE]

keep_cols <- !(is.na(col_names) | col_names == "")
dat_text  <- dat_text[, keep_cols, drop = FALSE]
names(dat_text) <- col_names[keep_cols]
if (anyDuplicated(names(dat_text))) stop("Duplicate column names after normalisation: ",
                                         paste(names(dat_text)[duplicated(names(dat_text))], collapse = ", "))

# Convert to numeric and count any values that fail to convert
dat <- dat_text
dat[] <- lapply(dat_text, function(x) suppressWarnings(as.numeric(x)))
n_new_na <- sum(is.na(dat)) - sum(is.na(dat_text))

# Drop completely empty rows (if any)
blank_rows <- rowSums(!is.na(dat)) == 0
n_blank <- sum(blank_rows)
dat <- dat[!blank_rows, , drop = FALSE]

# Harmonise column names across the common versions of this dataset
if (!"PAY_0" %in% names(dat) && "PAY_1" %in% names(dat)) names(dat)[names(dat) == "PAY_1"] <- "PAY_0"
target_col <- grep("^DEFAULT", names(dat), value = TRUE)
if (length(target_col) != 1) stop("Could not identify exactly one target column starting with DEFAULT.")
names(dat)[names(dat) == target_col] <- "default_next_month"

needed <- c("ID", "LIMIT_BAL", "SEX", "EDUCATION", "MARRIAGE", "AGE",
            paste0("PAY_", c(0, 2:6)), paste0("BILL_AMT", 1:6), paste0("PAY_AMT", 1:6),
            "default_next_month")
missing_cols <- setdiff(needed, names(dat))
if (length(missing_cols) > 0) stop("Required columns are missing: ", paste(missing_cols, collapse = ", "))
if (!all(dat$default_next_month %in% c(0, 1, NA))) stop("The target contains values other than 0 and 1.")

import_summary <- data.frame(
  Item = c("File used", "Sheet used", "Header found in row", "Header row inside the data",
           "Data rows imported", "Columns imported", "Values that failed numeric conversion",
           "Completely empty rows dropped", "First customer ID retained"),
  Value = c(basename(data_path), hdr$sheet, hdr$header_row,
            ifelse(hdr$header_row > 1, "Yes (rows above it were discarded)", "No (proper headers in row 1)"),
            nrow(dat), ncol(dat), n_new_na, n_blank, dat$ID[1])
)
knitr::kable(import_summary, caption = "Import summary")
Import summary
Item Value
File used default of credit card clients.xls
Sheet used Data
Header found in row 2
Header row inside the data Yes (rows above it were discarded)
Data rows imported 30000
Columns imported 25
Values that failed numeric conversion 0
Completely empty rows dropped 0
First customer ID retained 1 
dim_tab <- data.frame(Rows = nrow(dat), Columns = ncol(dat))
knitr::kable(dim_tab, caption = "Dimensions of the imported data")
Dimensions of the imported data
Rows Columns
30000 25 
type_tab <- data.frame(Variable = names(dat),
                       Type = vapply(dat, function(x) class(x)[1], character(1)),
                       Missing = vapply(dat, function(x) sum(is.na(x)), numeric(1)),
                       Distinct = vapply(dat, function(x) length(unique(x)), numeric(1)),
                       Min = vapply(dat, function(x) min(x, na.rm = TRUE), numeric(1)),
                       Max = vapply(dat, function(x) max(x, na.rm = TRUE), numeric(1)),
                       row.names = NULL)
knitr::kable(type_tab, digits = 0, format.args = list(big.mark = ","),
             caption = "Data types, missing values and ranges")
Data types, missing values and ranges
Variable Type Missing Distinct Min Max
ID numeric 0 30,000 1 30,000
LIMIT_BAL numeric 0 81 10,000 1,000,000
SEX numeric 0 2 1 2
EDUCATION numeric 0 7 0 6
MARRIAGE numeric 0 4 0 3
AGE numeric 0 56 21 79
PAY_0 numeric 0 11 -2 8
PAY_2 numeric 0 11 -2 8
PAY_3 numeric 0 11 -2 8
PAY_4 numeric 0 11 -2 8
PAY_5 numeric 0 10 -2 8
PAY_6 numeric 0 10 -2 8
BILL_AMT1 numeric 0 22,723 -165,580 964,511
BILL_AMT2 numeric 0 22,346 -69,777 983,931
BILL_AMT3 numeric 0 22,026 -157,264 1,664,089
BILL_AMT4 numeric 0 21,548 -170,000 891,586
BILL_AMT5 numeric 0 21,010 -81,334 927,171
BILL_AMT6 numeric 0 20,604 -339,603 961,664
PAY_AMT1 numeric 0 7,943 0 873,552
PAY_AMT2 numeric 0 7,899 0 1,684,259
PAY_AMT3 numeric 0 7,518 0 896,040
PAY_AMT4 numeric 0 6,937 0 621,000
PAY_AMT5 numeric 0 6,897 0 426,529
PAY_AMT6 numeric 0 6,939 0 528,666
default_next_month numeric 0 2 0 1 
dup_id      <- sum(duplicated(dat$ID))
dup_content <- sum(duplicated(dat[, setdiff(names(dat), "ID")]))
code_table <- function(x, nm) {
  t <- table(x, useNA = "ifany")
  data.frame(Variable = nm, Code = names(t), n = as.integer(t),
             Share = fmt_pct(as.integer(t) / length(x)))
}
knitr::kable(do.call(rbind, list(code_table(dat$SEX, "SEX"),
                                 code_table(dat$EDUCATION, "EDUCATION"),
                                 code_table(dat$MARRIAGE, "MARRIAGE"),
                                 code_table(dat$PAY_0, "PAY_0"))),
             caption = "Category codes as found in the file (before any recoding)")
Category codes as found in the file (before any recoding)
Variable Code n Share
SEX 1 11888 39.6%
SEX 2 18112 60.4%
EDUCATION 0 14 0.0%
EDUCATION 1 10585 35.3%
EDUCATION 2 14030 46.8%
EDUCATION 3 4917 16.4%
EDUCATION 4 123 0.4%
EDUCATION 5 280 0.9%
EDUCATION 6 51 0.2%
MARRIAGE 0 54 0.2%
MARRIAGE 1 13659 45.5%
MARRIAGE 2 15964 53.2%
MARRIAGE 3 323 1.1%
PAY_0 -2 2759 9.2%
PAY_0 -1 5686 19.0%
PAY_0 0 14737 49.1%
PAY_0 1 3688 12.3%
PAY_0 2 2667 8.9%
PAY_0 3 322 1.1%
PAY_0 4 76 0.3%
PAY_0 5 26 0.1%
PAY_0 6 11 0.0%
PAY_0 7 9 0.0%
PAY_0 8 19 0.1% 
unusual <- data.frame(
  Check = c("Duplicated customer IDs", "Duplicate rows ignoring ID",
            "Customers with a negative bill in at least one month",
            "Credit limits that are zero or negative", "Ages below 18 or above 100",
            "Repayment codes below -2 (PAY_0 to PAY_6)"),
  Count = c(dup_id, dup_content,
            sum(rowSums(dat[, paste0("BILL_AMT", 1:6)] < 0) > 0),
            sum(dat$LIMIT_BAL <= 0), sum(dat$AGE < 18 | dat$AGE > 100),
            sum(rowSums(dat[, paste0("PAY_", c(0, 2:6))] < -2) > 0))
)
knitr::kable(unusual, caption = "Checks for duplicates and unusual values")
Checks for duplicates and unusual values
Check Count
Duplicated customer IDs 0
Duplicate rows ignoring ID 35
Customers with a negative bill in at least one month 1930
Credit limits that are zero or negative 0
Ages below 18 or above 100 0
Repayment codes below -2 (PAY_0 to PAY_6) 0

Reading the checks.

  • Duplicate IDs. 0 customer ID(s) are duplicated. Rows are only removed if the customer ID itself is repeated, because the same customer must not appear twice.
  • Duplicate content. 35 rows are identical to another row once the ID is ignored. They are kept, because different customers can legitimately share the same values (for example, inactive accounts with all-zero bills). Removing them would be an unjustified deletion.
  • Negative bills. Negative bill amounts are kept as valid. They typically indicate a credit balance, that is, an overpayment.

3.2 Data cleaning and variable definitions

3.2.1 Cleaning and recoding decisions

cc <- dat
if (dup_id > 0) cc <- cc[!duplicated(cc$ID), ]      # only remove repeated IDs

ref_most_common <- function(f) {                    # reference level = largest category
  f <- factor(f)
  relevel(f, ref = names(which.max(table(f))))
}

# EDUCATION: 1 graduate school, 2 university, 3 high school, 4 others; 0, 5, 6 undocumented
edu_lab <- ifelse(cc$EDUCATION == 1, "Graduate school",
           ifelse(cc$EDUCATION == 2, "University",
           ifelse(cc$EDUCATION == 3, "High school", "Other or unknown")))
cc$EDUCATION <- ref_most_common(edu_lab)

# MARRIAGE: 1 married, 2 single, 3 others; 0 undocumented
mar_lab <- ifelse(cc$MARRIAGE == 1, "Married",
           ifelse(cc$MARRIAGE == 2, "Single", "Other or unknown"))
cc$MARRIAGE <- ref_most_common(mar_lab)

# Latest repayment status PAY_0:
#   -2 no consumption, -1 paid in full, 0 revolving credit (minimum paid),
#   1, 2, ... = months of payment delay. Negative codes are valid and are NOT recoded to 0.
#   Delays of 3 or more months are pooled because the individual categories are sparse
#   (the pooling rule uses category counts only, never the outcome).
pay_levels <- c("No consumption", "Paid in full", "Revolving credit, minimum paid",
                "1 month delay", "2 months delay", "3 or more months delay")
pay0_lab <- ifelse(cc$PAY_0 <= -2, pay_levels[1],
            ifelse(cc$PAY_0 == -1, pay_levels[2],
            ifelse(cc$PAY_0 == 0,  pay_levels[3],
            ifelse(cc$PAY_0 == 1,  pay_levels[4],
            ifelse(cc$PAY_0 == 2,  pay_levels[5], pay_levels[6])))))
cc$pay0_status <- ref_most_common(factor(pay0_lab, levels = pay_levels))

# Binary label for plots
cc$default_label <- factor(ifelse(cc$default_next_month == 1, "Default", "No default"),
                           levels = c("No default", "Default"))

recode_tab <- rbind(
  data.frame(Variable = "EDUCATION", Rule = "1 = Graduate school, 2 = University, 3 = High school; 4 and any undocumented code (0, 5, 6) pooled into 'Other or unknown'",
             Pooled_observations = sum(!dat$EDUCATION %in% 1:3)),
  data.frame(Variable = "MARRIAGE", Rule = "1 = Married, 2 = Single; 3 and undocumented 0 pooled into 'Other or unknown'",
             Pooled_observations = sum(!dat$MARRIAGE %in% 1:2)),
  data.frame(Variable = "PAY_0", Rule = "-2, -1, 0 kept as separate valid categories; 1, 2 kept; 3 or more months pooled",
             Pooled_observations = sum(dat$PAY_0 >= 3)),
  data.frame(Variable = "ID", Rule = "Kept only for tracing observations; never used as a predictor", Pooled_observations = 0)
)
knitr::kable(recode_tab, caption = "Recoding decisions and the number of observations affected")
Recoding decisions and the number of observations affected
Variable Rule Pooled_observations
EDUCATION 1 = Graduate school, 2 = University, 3 = High school; 4 and any undocumented code (0, 5, 6) pooled into ‘Other or unknown’ 468
MARRIAGE 1 = Married, 2 = Single; 3 and undocumented 0 pooled into ‘Other or unknown’ 377
PAY_0 -2, -1, 0 kept as separate valid categories; 1, 2 kept; 3 or more months pooled 463
ID Kept only for tracing observations; never used as a predictor 0

Why these decisions.

  • The undocumented education and marital codes are small groups with no reliable meaning. Pooling them with the documented “other” category keeps every observation without inventing a meaning for them.
  • The repayment status is not treated as a plain continuous number. The codes mix three ideas (no use, paid in full or minimum, and months late), and the step from -1 to 0 is not the same kind of step as from 1 to 2. As a factor, each status gets its own effect.
  • The reference level of each factor is its most common category, so coefficients compare with the largest group.
  • SEX is not used as a predictor. Using sex in a credit decision model is ethically and legally problematic in many jurisdictions, and it is not needed to answer the business question.

3.2.2 Variable definitions

var_defs <- data.frame(
  Variable = c("default_next_month", "pay0_status", "LIMIT_BAL", "AGE", "EDUCATION", "MARRIAGE",
               "earlier_delays_n", "BILL_AMT1 to BILL_AMT6", "PAY_AMT1 to PAY_AMT6",
               "bill_trend", "payment_trend"),
  Meaning = c("Target: 1 = default next month, 0 = no default",
              "Latest repayment status (factor built from PAY_0)",
              "Credit limit in NT dollars",
              "Age in years",
              "Education (factor)", "Marital status (factor)",
              "Number of the five earlier months (PAY_2 to PAY_6) with a payment delay of one month or more (0 to 5)",
              "Monthly bill statement amounts in NT dollars (1 = most recent)",
              "Monthly payment amounts in NT dollars (1 = most recent)",
              "Average of the 3 most recent bills minus the average of the 3 earlier bills (NT dollars)",
              "Average of the 3 most recent payments minus the average of the 3 earlier payments (NT dollars)")
)
knitr::kable(var_defs, caption = "Variables used in the analysis")
Variables used in the analysis
Variable Meaning
default_next_month Target: 1 = default next month, 0 = no default
pay0_status Latest repayment status (factor built from PAY_0)
LIMIT_BAL Credit limit in NT dollars
AGE Age in years
EDUCATION Education (factor)
MARRIAGE Marital status (factor)
earlier_delays_n Number of the five earlier months (PAY_2 to PAY_6) with a payment delay of one month or more (0 to 5)
BILL_AMT1 to BILL_AMT6 Monthly bill statement amounts in NT dollars (1 = most recent)
PAY_AMT1 to PAY_AMT6 Monthly payment amounts in NT dollars (1 = most recent)
bill_trend Average of the 3 most recent bills minus the average of the 3 earlier bills (NT dollars)
payment_trend Average of the 3 most recent payments minus the average of the 3 earlier payments (NT dollars)

3.3 Feature engineering

3.3.1 Confirming the monthly ordering

The UCI documentation states that BILL_AMT1 and PAY_AMT1 refer to September (the most recent month), and that higher numbers refer to earlier months back to April (BILL_AMT6, PAY_AMT6). The data cannot prove the calendar, but two simple patterns should hold if the numbering runs from recent to older:

  1. Bills in adjacent months are more strongly correlated than bills six months apart.
  2. The most recent repayment status should be more strongly related to next month’s default than the oldest one.
bill_cols <- paste0("BILL_AMT", 1:6)
pay_cols  <- paste0("PAY_AMT", 1:6)
month_lab <- c("Sep (latest)", "Aug", "Jul", "Jun", "May", "Apr")

adj_cor <- sapply(1:5, function(k) cor(cc[[bill_cols[k]]], cc[[bill_cols[k + 1]]]))
far_cor <- cor(cc$BILL_AMT1, cc$BILL_AMT6)
pay_status_cols <- paste0("PAY_", c(0, 2:6))
cor_default <- sapply(pay_status_cols, function(v) cor(cc[[v]], cc$default_next_month))

order_ok_1 <- far_cor < min(adj_cor)
order_ok_2 <- abs(cor_default["PAY_0"]) > abs(cor_default["PAY_6"])

ord_tab <- data.frame(
  Check = c(paste0("Correlation of ", bill_cols[1:5], " with ", bill_cols[2:6]),
            "Correlation of BILL_AMT1 with BILL_AMT6",
            paste0("Correlation of default with ", pay_status_cols)),
  Correlation = c(adj_cor, far_cor, cor_default)
)
knitr::kable(ord_tab, digits = 3, caption = "Evidence on the temporal ordering of the monthly columns")
Evidence on the temporal ordering of the monthly columns
Check Correlation
Correlation of BILL_AMT1 with BILL_AMT2 0.951
Correlation of BILL_AMT2 with BILL_AMT3 0.928
Correlation of BILL_AMT3 with BILL_AMT4 0.924
Correlation of BILL_AMT4 with BILL_AMT5 0.940
Correlation of BILL_AMT5 with BILL_AMT6 0.946
Correlation of BILL_AMT1 with BILL_AMT6 0.803
Correlation of default with PAY_0 0.325
Correlation of default with PAY_2 0.264
Correlation of default with PAY_3 0.235
Correlation of default with PAY_4 0.217
Correlation of default with PAY_5 0.204
Correlation of default with PAY_6 0.187

Adjacent months are more strongly correlated than the most distant pair, and the latest repayment status is more strongly related to default than the oldest. Both patterns are consistent with the documented ordering (1 = most recent), so the documented ordering is used.

3.3.2 Creating the trend features

cc$bill_recent    <- rowMeans(cc[, bill_cols[1:3]])
cc$bill_earlier   <- rowMeans(cc[, bill_cols[4:6]])
cc$bill_trend     <- cc$bill_recent - cc$bill_earlier          # NT dollars

cc$payment_recent  <- rowMeans(cc[, pay_cols[1:3]])
cc$payment_earlier <- rowMeans(cc[, pay_cols[4:6]])
cc$payment_trend   <- cc$payment_recent - cc$payment_earlier    # NT dollars

# Earlier repayment behaviour: delayed months among PAY_2 to PAY_6
for (k in 2:6) cc[[paste0("delay_", k)]] <- as.integer(cc[[paste0("PAY_", k)]] >= 1)
cc$earlier_delays_n <- rowSums(cc[, paste0("delay_", 2:6)])

# Model-scale versions in units of NT$10,000 so coefficients are readable
cc$limit_10k          <- cc$LIMIT_BAL / 1e4
cc$bill_latest_10k    <- cc$BILL_AMT1 / 1e4
cc$pay_latest_10k     <- cc$PAY_AMT1 / 1e4
cc$bill_trend_10k     <- cc$bill_trend / 1e4
cc$payment_trend_10k  <- cc$payment_trend / 1e4
  • Bill trend = mean(BILL_AMT1, BILL_AMT2, BILL_AMT3) minus mean(BILL_AMT4, BILL_AMT5, BILL_AMT6). The unit is NT dollars. A positive value means recent bills are higher than earlier bills.
  • Payment trend = mean(PAY_AMT1, PAY_AMT2, PAY_AMT3) minus mean(PAY_AMT4, PAY_AMT5, PAY_AMT6). The unit is NT dollars. A negative value means recent payments are lower than earlier payments.
  • In the regression models both trends are expressed in units of NT$10,000 (bill_trend_10k, payment_trend_10k), so a coefficient describes the effect of a NT$10,000 change in the trend.
knitr::kable(data.frame(
  Feature = c("bill_trend (NT$)", "payment_trend (NT$)"),
  `Share positive` = c(fmt_pct(mean(cc$bill_trend > 0)), fmt_pct(mean(cc$payment_trend > 0))),
  `Share negative` = c(fmt_pct(mean(cc$bill_trend < 0)), fmt_pct(mean(cc$payment_trend < 0))),
  `Share exactly zero` = c(fmt_pct(mean(cc$bill_trend == 0)), fmt_pct(mean(cc$payment_trend == 0))),
  check.names = FALSE),
  caption = "Direction of the engineered trend features")
Direction of the engineered trend features
Feature Share positive Share negative Share exactly zero
bill_trend (NT\() |56.5% |39.1% |4.3% | |payment_trend (NT\)) 62.3% 31.4% 6.3%

3.4 Exploratory data analysis

The exploratory work below is descriptive and uses all rows. No modelling choice (variable selection, thresholds, transformations of predictors) is based on it. All of those choices are made later on the training data only.

3.4.1 Summary statistics and the target

sumstat <- function(x) c(Mean = mean(x), SD = sd(x), Min = min(x),
                         Q1 = unname(quantile(x, .25)), Median = median(x),
                         Q3 = unname(quantile(x, .75)), Max = max(x))
key_df <- data.frame(`Credit limit (NT$)` = cc$LIMIT_BAL, `Age (years)` = cc$AGE,
                     `Latest bill (NT$)` = cc$BILL_AMT1, `Latest payment (NT$)` = cc$PAY_AMT1,
                     `Bill trend (NT$)` = cc$bill_trend, `Payment trend (NT$)` = cc$payment_trend,
                     `Earlier delayed months (0 to 5)` = cc$earlier_delays_n, check.names = FALSE)
knitr::kable(t(sapply(key_df, sumstat)), digits = 1, format.args = list(big.mark = ","),
             caption = "Summary statistics for key numerical variables")
Summary statistics for key numerical variables
Mean SD Min Q1 Median Q3 Max
Credit limit (NT\() | 167,484.3| 129,747.7| 10,000| 50,000.0| 140,000.0| 240,000| 1,000,000| |Age (years) | 35.5| 9.2| 21| 28.0| 34.0| 41| 79| |Latest bill (NT\)) 51,223.3 73,635.9 -165,580 3,558.8 22,381.5 67,091 964,511
Latest payment (NT\() | 5,663.6| 16,563.3| 0| 1,000.0| 2,100.0| 5,006| 873,552| |Bill trend (NT\)) 8,323.2 29,831.2 -291,275 -1,899.8 658.7 12,647 432,428
Payment trend (NT$) 656.5 13,892.6 -177,931 -316.7 403.3 1,513 730,689
Earlier delayed months (0 to 5) 0.6 1.3 0 0.0 0.0 1 5 
tgt <- table(cc$default_label)
tgt_tab <- data.frame(Class = names(tgt), Count = as.integer(tgt), Share = fmt_pct(as.integer(tgt) / sum(tgt)))
knitr::kable(tgt_tab, caption = "Target class counts and observed default rate")
Target class counts and observed default rate
Class Count Share
No default 23364 77.9%
Default 6636 22.1% 
default_rate_all <- mean(cc$default_next_month)

ggplot(tgt_tab, aes(Class, Count, fill = Class)) +
  geom_col(width = 0.6, show.legend = FALSE) +
  geom_text(aes(label = paste0(scales::comma(Count), " (", Share, ")")), vjust = -0.4) +
  scale_fill_manual(values = col_def) + scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.12))) +
  labs(title = "Next-month default: class balance", x = NULL, y = "Customers")

The observed default rate is 22.1% across 30,000 customers. The classes are imbalanced, which is why the train/test split is stratified and why accuracy alone would be a poor metric.

3.4.2 Credit limits, bills and payments

Monetary variables are heavily right-skewed. For visualisation only:

  • Credit limits are always positive, so a plain log10 axis is used.
  • Bills, payments and trends can be zero or negative, so a signed log transform is used: sign(x) * log10(1 + abs(x)). It handles zero and negative values, and it is symmetric around zero.

The models use the original NT dollar scale (rescaled to NT$10,000), not these transformed values.

ggplot(cc, aes(LIMIT_BAL)) +
  geom_histogram(bins = 40, fill = col_main, colour = "white") +
  geom_vline(xintercept = median(cc$LIMIT_BAL), linetype = "dashed") +
  scale_x_log10(labels = scales::comma) +
  labs(title = "Distribution of credit limits (log10 axis)",
       subtitle = paste0("Dashed line: median = NT$", scales::comma(median(cc$LIMIT_BAL))),
       x = "Credit limit (NT$, log scale)", y = "Customers")

slog10 <- function(x) sign(x) * log10(1 + abs(x))
long_df <- function(prefix, labels) {
  out <- do.call(rbind, lapply(1:6, function(i)
    data.frame(month = labels[i], value = cc[[paste0(prefix, i)]])))
  out$month <- factor(out$month, levels = labels)
  out
}
bills_long <- long_df("BILL_AMT", month_lab)
pays_long  <- long_df("PAY_AMT", month_lab)

ggplot(bills_long, aes(slog10(value))) +
  geom_histogram(bins = 40, fill = col_main, colour = "white") +
  facet_wrap(~ month, ncol = 3) +
  labs(title = "Distribution of monthly bill amounts (signed log10 scale)",
       x = "Signed log10 of bill amount (NT$)", y = "Customers")

ggplot(pays_long, aes(slog10(value))) +
  geom_histogram(bins = 40, fill = "#41AB5D", colour = "white") +
  facet_wrap(~ month, ncol = 3) +
  labs(title = "Distribution of monthly payment amounts (signed log10 scale)",
       x = "Signed log10 of payment amount (NT$)", y = "Customers")

3.4.4 Correlation of numerical predictors

num_df <- data.frame(`Credit limit` = cc$LIMIT_BAL, Age = cc$AGE,
                     `Latest bill` = cc$BILL_AMT1, `Latest payment` = cc$PAY_AMT1,
                     `Bill trend` = cc$bill_trend, `Payment trend` = cc$payment_trend,
                     `Earlier delays` = cc$earlier_delays_n, `PAY_0 (code)` = cc$PAY_0,
                     Default = cc$default_next_month, check.names = FALSE)
cm <- cor(num_df, use = "pairwise.complete.obs")
cm_long <- as.data.frame(as.table(cm)); names(cm_long) <- c("Var1", "Var2", "r")

ggplot(cm_long, aes(Var1, Var2, fill = r)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = sprintf("%.2f", r)), size = 3) +
  scale_fill_gradient2(low = "#B2182B", mid = "white", high = "#2166AC", limits = c(-1, 1), name = "r") +
  labs(title = "Correlation of numerical variables", x = NULL, y = NULL,
       subtitle = "PAY_0 is shown as its raw code only for this correlation view") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

3.4.5 Default rate by latest repayment status

pay_tab <- do.call(rbind, lapply(levels(cc$pay0_status), function(l) {
  s <- cc$default_next_month[cc$pay0_status == l]
  p <- mean(s); n <- length(s)
  data.frame(Status = l, Customers = n, Defaults = sum(s), Rate = p,
             Low = max(0, p - 1.96 * sqrt(p * (1 - p) / n)), High = min(1, p + 1.96 * sqrt(p * (1 - p) / n)))
}))
pay_tab$Status <- factor(pay_tab$Status, levels = pay_levels)
pay_tab <- pay_tab[order(pay_tab$Status), ]
show_pay <- pay_tab; show_pay$Rate <- fmt_pct(show_pay$Rate); show_pay$Low <- NULL; show_pay$High <- NULL
knitr::kable(show_pay, row.names = FALSE, format.args = list(big.mark = ","),
             caption = "Observed default rate by latest repayment status")
Observed default rate by latest repayment status
Status Customers Defaults Rate
No consumption 2,759 365 13.2%
Paid in full 5,686 954 16.8%
Revolving credit, minimum paid 14,737 1,888 12.8%
1 month delay 3,688 1,252 33.9%
2 months delay 2,667 1,844 69.1%
3 or more months delay 463 333 71.9% 
ggplot(pay_tab, aes(Status, Rate)) +
  geom_col(fill = col_main, width = 0.7) +
  geom_errorbar(aes(ymin = Low, ymax = High), width = 0.2) +
  geom_hline(yintercept = default_rate_all, linetype = "dashed") +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Default rate by latest repayment status",
       subtitle = "Error bars: approximate 95% intervals; dashed line: overall default rate",
       x = NULL, y = "Observed default rate") +
  theme(axis.text.x = element_text(angle = 25, hjust = 1))

4 Modeling process

4.1 Training and testing split

The data are split once, before any model selection, using a stratified 80/20 random split so that both sets have the same default rate (up to rounding). The same split is used for every model.

model_vars <- c("ID", "default_next_month", "pay0_status", "limit_10k", "AGE", "EDUCATION", "MARRIAGE",
                "earlier_delays_n", paste0("delay_", 2:6), "bill_latest_10k", "pay_latest_10k",
                "bill_trend_10k", "payment_trend_10k")
model_df_all <- cc[, model_vars]
n_dropped_na <- sum(!complete.cases(model_df_all))
model_df <- model_df_all[complete.cases(model_df_all), ]

set.seed(params$seed)
idx1 <- which(model_df$default_next_month == 1)
idx0 <- which(model_df$default_next_month == 0)
train_idx <- sort(c(sample(idx1, round(params$train_share * length(idx1))),
                    sample(idx0, round(params$train_share * length(idx0)))))
train <- model_df[train_idx, ]
test  <- model_df[-train_idx, ]

split_tab <- data.frame(
  Set = c("Training", "Testing", "Total"),
  Observations = c(nrow(train), nrow(test), nrow(model_df)),
  Defaults = c(sum(train$default_next_month), sum(test$default_next_month), sum(model_df$default_next_month)),
  `Default rate` = fmt_pct(c(mean(train$default_next_month), mean(test$default_next_month), mean(model_df$default_next_month)), 2),
  check.names = FALSE)
knitr::kable(split_tab, format.args = list(big.mark = ","), caption = "Stratified training and testing split")
Stratified training and testing split
Set Observations Defaults Default rate
Training 24,000 5,309 22.12%
Testing 6,000 1,327 22.12%
Total 30,000 6,636 22.12% 
y_train_mean <- mean(train$default_next_month)

0 observation(s) were removed for missing model variables. The seed is 2026.

trend_terms <- c("bill_trend_10k", "payment_trend_10k")

coef_table_lm <- function(model, robust = TRUE) {
  vc <- if (robust) sandwich::vcovHC(model, type = "HC3") else stats::vcov(model)
  ct <- lmtest::coeftest(model, vcov. = vc)
  ci <- lmtest::coefci(model, vcov. = vc)
  out <- data.frame(Term = rownames(ct), Estimate = ct[, 1], SE = ct[, 2], `t value` = ct[, 3],
                    `p value` = ct[, 4], `CI low` = ci[, 1], `CI high` = ci[, 2],
                    check.names = FALSE, row.names = NULL)
  out$`p value` <- fmt_p(out$`p value`)
  out
}
get_coef <- function(model, term) {
  vc <- sandwich::vcovHC(model, type = "HC3")
  ct <- lmtest::coeftest(model, vcov. = vc)
  ci <- lmtest::coefci(model, vcov. = vc)
  data.frame(estimate = ct[term, 1], se = ct[term, 2], p = ct[term, 4], lo = ci[term, 1], hi = ci[term, 2])
}

4.2 Baseline model

4.2.1 Specification and the question of earlier repayment history

Baseline predictors (pre-specified): latest repayment status (pay0_status), credit limit, age, education, marital status, and the count of earlier delayed months (earlier_delays_n).

The five earlier repayment statuses (PAY_2 to PAY_6) are highly correlated with each other and with PAY_0. Adding all of them mechanically would add redundancy and variance. The compromise is one compact summary: how many of the five earlier months were delayed. To check that this choice is defensible, three specifications are compared on the training data only:

base_core  <- c("pay0_status", "limit_10k", "AGE", "EDUCATION", "MARRIAGE")
spec_a <- lm(reformulate(base_core, "default_next_month"), data = train)
spec_b <- lm(reformulate(c(base_core, "earlier_delays_n"), "default_next_month"), data = train)
spec_c <- lm(reformulate(c(base_core, paste0("delay_", 2:6)), "default_next_month"), data = train)

spec_tab <- data.frame(
  Specification = c("(a) No earlier history", "(b) Count of earlier delayed months", "(c) Five separate delay indicators"),
  Parameters = sapply(list(spec_a, spec_b, spec_c), function(m) length(coef(m))),
  AIC = sapply(list(spec_a, spec_b, spec_c), AIC),
  BIC = sapply(list(spec_a, spec_b, spec_c), BIC),
  `Adjusted R2 (training)` = sapply(list(spec_a, spec_b, spec_c), function(m) summary(m)$adj.r.squared),
  check.names = FALSE)
knitr::kable(spec_tab, digits = 4, caption = "Training-data comparison of ways to include earlier repayment history")
Training-data comparison of ways to include earlier repayment history
Specification Parameters AIC BIC Adjusted R2 (training)
(a) No earlier history 13 20861 20974 0.1900
(b) Count of earlier delayed months 14 20353 20474 0.2070
(c) Five separate delay indicators 18 20358 20512 0.2069 
best_bic_spec <- spec_tab$Specification[which.min(spec_tab$BIC)]

The lowest BIC belongs to (b) Count of earlier delayed months. Specification (b) is used as the baseline because it captures earlier delinquency with a single, interpretable variable and avoids the redundancy of five correlated indicators. The BIC table above shows how it compares with the alternatives.

4.2.2 Baseline linear probability model

base_terms <- c(base_core, "earlier_delays_n")
lpm_base <- lm(reformulate(base_terms, "default_next_month"), data = train)
knitr::kable(coef_table_lm(lpm_base), digits = 4,
             caption = "Baseline LPM, training data (HC3 robust standard errors; coefficients are changes in probability, so 0.01 = 1 percentage point)")
Baseline LPM, training data (HC3 robust standard errors; coefficients are changes in probability, so 0.01 = 1 percentage point)
Term Estimate SE t value p value CI low CI high
(Intercept) 0.1209 0.0107 11.2793 < 0.001 0.0999 0.1419
pay0_statusNo consumption 0.0370 0.0080 4.6517 < 0.001 0.0214 0.0527
pay0_statusPaid in full 0.0450 0.0063 7.1128 < 0.001 0.0326 0.0574
pay0_status1 month delay 0.1468 0.0095 15.4245 < 0.001 0.1281 0.1654
pay0_status2 months delay 0.4476 0.0125 35.7549 < 0.001 0.4231 0.4721
pay0_status3 or more months delay 0.3868 0.0252 15.3513 < 0.001 0.3374 0.4362
limit_10k -0.0024 0.0002 -12.6027 < 0.001 -0.0027 -0.0020
AGE 0.0008 0.0003 2.4642 0.014 0.0002 0.0014
EDUCATIONGraduate school -0.0020 0.0055 -0.3588 0.720 -0.0127 0.0088
EDUCATIONHigh school -0.0143 0.0074 -1.9320 0.053 -0.0288 0.0002
EDUCATIONOther or unknown -0.0924 0.0145 -6.3882 < 0.001 -0.1208 -0.0641
MARRIAGEMarried 0.0212 0.0056 3.7951 < 0.001 0.0102 0.0321
MARRIAGEOther or unknown 0.0114 0.0239 0.4763 0.634 -0.0354 0.0582
earlier_delays_n 0.0536 0.0030 18.1402 < 0.001 0.0479 0.0594

Training fit: R-squared = 0.2074, adjusted R-squared = 0.2070. Training R-squared is not used as evidence of predictive success. Out-of-sample metrics are in the validation section.

4.3 Expanded model

The expanded model adds the two engineered trend features to the baseline.

exp_terms <- c(base_terms, trend_terms)
lpm_exp <- lm(reformulate(exp_terms, "default_next_month"), data = train)
knitr::kable(coef_table_lm(lpm_exp), digits = 4,
             caption = "Expanded LPM, training data (HC3 robust standard errors)")
Expanded LPM, training data (HC3 robust standard errors)
Term Estimate SE t value p value CI low CI high
(Intercept) 0.1188 0.0107 11.0649 < 0.001 0.0977 0.1398
pay0_statusNo consumption 0.0414 0.0081 5.1274 < 0.001 0.0256 0.0573
pay0_statusPaid in full 0.0495 0.0065 7.6244 < 0.001 0.0367 0.0622
pay0_status1 month delay 0.1492 0.0096 15.6143 < 0.001 0.1304 0.1679
pay0_status2 months delay 0.4480 0.0125 35.7806 < 0.001 0.4235 0.4726
pay0_status3 or more months delay 0.3860 0.0252 15.3247 < 0.001 0.3366 0.4354
limit_10k -0.0024 0.0002 -12.8923 < 0.001 -0.0028 -0.0021
AGE 0.0008 0.0003 2.3965 0.017 0.0001 0.0014
EDUCATIONGraduate school -0.0016 0.0055 -0.2891 0.773 -0.0123 0.0091
EDUCATIONHigh school -0.0143 0.0074 -1.9297 0.054 -0.0288 0.0002
EDUCATIONOther or unknown -0.0943 0.0145 -6.5249 < 0.001 -0.1227 -0.0660
MARRIAGEMarried 0.0213 0.0056 3.8178 < 0.001 0.0104 0.0322
MARRIAGEOther or unknown 0.0124 0.0239 0.5191 0.604 -0.0344 0.0592
earlier_delays_n 0.0542 0.0030 18.2769 < 0.001 0.0483 0.0600
bill_trend_10k 0.0024 0.0008 3.0614 0.002 0.0009 0.0039
payment_trend_10k -0.0027 0.0012 -2.2451 0.025 -0.0051 -0.0003
Model Predictors
Baseline latest repayment status, credit limit, age, education, marital status, earlier delayed months
Expanded baseline predictors plus bill_trend and payment_trend (NT$10,000 units)

4.4 Variable selection

4.4.1 Method

  • Procedure: bidirectional stepwise regression (step(), direction "both") on the training data only, applied to the linear probability model.
  • Criterion: BIC (k = log(n)). BIC penalises complexity more heavily than AIC. With roughly 24,000 training rows, AIC would keep almost every candidate even when its gain is trivial.
  • Starting model (and lower bound): latest repayment status plus the two trend features. These are the variables in the business question, so they are forced to stay. The procedure therefore selects the other controls, and the research question about the trends is never settled by the selection algorithm itself.
  • Upper scope (candidates): the baseline predictors plus the latest bill and the latest payment (in NT$10,000). Individual monthly bill and payment columns are not offered, because they are strongly correlated with each other and with the trend features, which would create redundant predictors.
  • SEX and the customer ID are not candidates.
lower_terms <- c("pay0_status", trend_terms)
upper_terms <- c(base_terms, "bill_latest_10k", "pay_latest_10k", trend_terms)
f_lower <- reformulate(lower_terms, "default_next_month")
f_upper <- reformulate(upper_terms, "default_next_month")

start_lpm <- lm(f_lower, data = train)
step_lpm  <- step(start_lpm, scope = list(lower = f_lower, upper = f_upper),
                  direction = "both", k = log(nrow(train)), trace = 0)

path <- step_lpm$anova
names(path)[names(path) == "AIC"] <- "BIC"
knitr::kable(path, digits = 2, caption = "Stepwise path (criterion shown is BIC)")
Stepwise path (criterion shown is BIC)
Step Df Deviance Resid. Df Resid. Dev BIC
23992 3388 -46909
+ earlier_delays_n -1 85.99 23991 3302 -47516
+ limit_10k -1 17.84 23990 3284 -47636
+ MARRIAGE -2 4.07 23988 3280 -47645
+ pay_latest_10k -1 2.06 23987 3278 -47650 
sel_terms <- attr(terms(step_lpm), "term.labels")

# Sensitivity: does AIC give a different answer?
step_aic <- step(start_lpm, scope = list(lower = f_lower, upper = f_upper), direction = "both", k = 2, trace = 0)
aic_terms <- attr(terms(step_aic), "term.labels")

Result. Starting from pay0_status, bill_trend_10k and payment_trend_10k, the BIC procedure ended with: pay0_status, bill_trend_10k, payment_trend_10k, earlier_delays_n, limit_10k, MARRIAGE and pay_latest_10k. Candidates not retained: AGE, EDUCATION and bill_latest_10k. As a sensitivity check, AIC-based selection also retains: EDUCATION, bill_latest_10k and AGE.

Stepwise selection is a data-driven procedure. The p-values from the selected model are therefore somewhat optimistic about the controls, which is another reason the test set is used for the predictive claims.

4.5 OLS diagnostics and remedies

The diagnostics are for the selected linear probability model. The selected specification is first checked for multicollinearity, and the final model is the result of that check.

4.5.1 Multicollinearity

  • VIF is computed with car::vif().
  • Factors have several coefficients, so the generalised VIF (GVIF) is used and converted to a VIF-equivalent scale, GVIF^(1/Df), which is directly comparable with the usual threshold.
  • Threshold: a VIF-equivalent above 5 is flagged. Flagged variables that are not part of the research question are removed one at a time (largest first), the model is refitted, and VIF is recalculated. The repayment status and the two trend features are protected and are never dropped automatically.
vif_table <- function(model, threshold = 5) {
  v <- car::vif(model)
  if (is.matrix(v)) {
    out <- data.frame(term = rownames(v), GVIF = v[, "GVIF"], Df = v[, "Df"], row.names = NULL)
  } else {
    out <- data.frame(term = names(v), GVIF = as.numeric(v), Df = 1, row.names = NULL)
  }
  out$vif_equiv <- out$GVIF^(1 / out$Df)
  out$flag <- out$vif_equiv > threshold
  out
}

prune_vif <- function(model, protected, threshold = 5) {
  dropped <- NULL
  repeat {
    vt <- vif_table(model, threshold)
    bad <- vt[vt$flag & !(vt$term %in% protected), ]
    if (nrow(bad) == 0) break
    worst <- bad$term[which.max(bad$vif_equiv)]
    dropped <- rbind(dropped, data.frame(Dropped = worst, `VIF equivalent` = max(bad$vif_equiv), check.names = FALSE))
    model <- update(model, as.formula(paste(". ~ . -", worst)))
  }
  list(model = model, dropped = dropped)
}

vif_before <- vif_table(step_lpm)
vif_show <- function(vt) data.frame(Term = vt$term, GVIF = vt$GVIF, Df = vt$Df,
                                    `VIF equivalent` = vt$vif_equiv,
                                    Flag = ifelse(vt$flag, "Above 5", "OK"), check.names = FALSE)
knitr::kable(vif_show(vif_before), digits = 3, caption = "VIF for the stepwise-selected model (before any remedy)")
VIF for the stepwise-selected model (before any remedy)
Term GVIF Df VIF equivalent Flag
pay0_status 1.786 5 1.123 OK
bill_trend_10k 1.095 1 1.095 OK
payment_trend_10k 1.305 1 1.305 OK
earlier_delays_n 1.610 1 1.610 OK
limit_10k 1.197 1 1.197 OK
MARRIAGE 1.017 2 1.009 OK
pay_latest_10k 1.365 1 1.365 OK 
pr <- prune_vif(step_lpm, protected = c("pay0_status", trend_terms))
final_lpm <- pr$model
final_terms <- attr(terms(final_lpm), "term.labels")

vif_after <- vif_table(final_lpm)
knitr::kable(vif_show(vif_after), digits = 3, caption = "VIF for the final linear probability model (after remedy, if any)")
VIF for the final linear probability model (after remedy, if any)
Term GVIF Df VIF equivalent Flag
pay0_status 1.786 5 1.123 OK
bill_trend_10k 1.095 1 1.095 OK
payment_trend_10k 1.305 1 1.305 OK
earlier_delays_n 1.610 1 1.610 OK
limit_10k 1.197 1 1.197 OK
MARRIAGE 1.017 2 1.009 OK
pay_latest_10k 1.365 1 1.365 OK

No predictor exceeded the threshold of 5 (largest VIF equivalent: 1.61), so no variable was removed and the stepwise-selected model is the final model. After this step no term exceeds the threshold.

knitr::kable(coef_table_lm(final_lpm), digits = 4,
             caption = "Final LPM, training data, HC3 robust standard errors")
Final LPM, training data, HC3 robust standard errors
Term Estimate SE t value p value CI low CI high
(Intercept) 0.1378 0.0049 28.1156 < 0.001 0.1282 0.1474
pay0_statusNo consumption 0.0404 0.0080 5.0245 < 0.001 0.0246 0.0561
pay0_statusPaid in full 0.0506 0.0065 7.8423 < 0.001 0.0380 0.0632
pay0_status1 month delay 0.1482 0.0096 15.5132 < 0.001 0.1295 0.1669
pay0_status2 months delay 0.4472 0.0125 35.6735 < 0.001 0.4226 0.4717
pay0_status3 or more months delay 0.3851 0.0252 15.2756 < 0.001 0.3357 0.4345
bill_trend_10k 0.0025 0.0008 3.2604 0.001 0.0010 0.0041
payment_trend_10k 0.0008 0.0016 0.4998 0.617 -0.0024 0.0040
earlier_delays_n 0.0546 0.0030 18.4346 < 0.001 0.0488 0.0604
limit_10k -0.0022 0.0002 -11.7391 < 0.001 -0.0025 -0.0018
MARRIAGEMarried 0.0262 0.0049 5.3854 < 0.001 0.0167 0.0357
MARRIAGEOther or unknown 0.0190 0.0237 0.8014 0.423 -0.0275 0.0656
pay_latest_10k -0.0063 0.0014 -4.6070 < 0.001 -0.0090 -0.0036

4.5.2 Heteroscedasticity

Why it is expected. For a binary outcome, the variance of the error depends on the predicted probability, p(1 - p). It is largest near 0.5 and smaller near 0 or 1, so constant error variance cannot hold for an LPM unless every predicted probability is the same.

diag_df <- data.frame(fitted = fitted(final_lpm), resid = resid(final_lpm))
diag_df$abs_resid <- abs(diag_df$resid)
diag_df$bin <- cut(rank(diag_df$fitted, ties.method = "first"),
                   breaks = seq(0, nrow(diag_df), length.out = 21), labels = FALSE, include.lowest = TRUE)
bin_summary <- aggregate(cbind(fitted, resid, abs_resid) ~ bin, data = diag_df, FUN = mean)

set.seed(params$seed)
samp <- diag_df[sample(nrow(diag_df), min(5000, nrow(diag_df))), ]

ggplot(samp, aes(fitted, resid)) +
  geom_point(alpha = 0.15, colour = "grey40") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(data = bin_summary, aes(fitted, resid), colour = "#D95F0E", linewidth = 1) +
  geom_point(data = bin_summary, aes(fitted, resid), colour = "#D95F0E", size = 2) +
  labs(title = "Residuals versus fitted values (final LPM, training data)",
       subtitle = "Grey: random sample of 5,000 points; orange: mean residual within 20 bins of fitted values",
       x = "Fitted probability of default", y = "Residual")

The two bands of points are expected: for each fitted value the residual can only take two values (1 - fitted if the customer defaulted, -fitted if not). The scatter is therefore much less informative than it would be for a continuous outcome, and the formal test and the theory carry more weight.

bp <- lmtest::bptest(final_lpm)
bp_tab <- data.frame(Test = "Breusch-Pagan (studentised)", Statistic = unname(bp$statistic),
                     Df = unname(bp$parameter), `p value` = fmt_p(bp$p.value), check.names = FALSE)
knitr::kable(bp_tab, digits = 3, caption = "Breusch-Pagan test on the final LPM")
Breusch-Pagan test on the final LPM
Test Statistic Df p value
BP Breusch-Pagan (studentised) 1094 12 < 0.001 
ggplot(diag_df, aes(sample = resid)) +
  stat_qq(alpha = 0.2, colour = "grey30") + stat_qq_line(colour = "#D95F0E") +
  labs(title = "Normal Q-Q plot of LPM residuals (training data)",
       x = "Theoretical quantiles", y = "Sample quantiles")

Result. The Breusch-Pagan test gives a statistic of 1,093.59 with 12 degrees of freedom and a p-value of < 0.001. The null hypothesis of constant error variance is rejected, which agrees with the theoretical expectation for a binary outcome.

Remedy. Heteroscedasticity-robust (HC3) standard errors are used for all inference, as in the coefficient tables above. The comparison of classical and robust standard errors for the key terms:

cls <- coef_table_lm(final_lpm, robust = FALSE); rob <- coef_table_lm(final_lpm, robust = TRUE)
se_cmp <- data.frame(Term = cls$Term, `Classical SE` = cls$SE, `HC3 robust SE` = rob$SE,
                     `Ratio (robust / classical)` = rob$SE / cls$SE, check.names = FALSE)
knitr::kable(se_cmp, digits = 4, caption = "Classical versus robust standard errors")
Classical versus robust standard errors
Term Classical SE HC3 robust SE Ratio (robust / classical)
(Intercept) 0.0051 0.0049 0.9618
pay0_statusNo consumption 0.0089 0.0080 0.9017
pay0_statusPaid in full 0.0067 0.0065 0.9606
pay0_status1 month delay 0.0081 0.0096 1.1773
pay0_status2 months delay 0.0100 0.0125 1.2525
pay0_status3 or more months delay 0.0206 0.0252 1.2211
bill_trend_10k 0.0008 0.0008 0.9321
payment_trend_10k 0.0019 0.0016 0.8342
earlier_delays_n 0.0024 0.0030 1.2507
limit_10k 0.0002 0.0002 0.9258
MARRIAGEMarried 0.0048 0.0049 1.0038
MARRIAGEOther or unknown 0.0221 0.0237 1.0721
pay_latest_10k 0.0016 0.0014 0.8437

What this does and does not fix. Heteroscedasticity does not bias the OLS coefficients and does not by itself change the fitted values, so point predictions are unaffected. It makes the classical standard errors, t-tests and confidence intervals unreliable, which is what the robust standard errors address. Robust standard errors do not repair the other LPM limitations discussed below, notably predicted values outside the 0 to 1 range.

4.5.3 Autocorrelation

The Durbin-Watson test is designed for time-ordered observations. Here each row is a different customer and the row order is arbitrary (by customer ID), so there is no meaningful chronology. The statistic is reported only for completeness:

dw <- lmtest::dwtest(final_lpm)
knitr::kable(data.frame(Statistic = unname(dw$statistic), `p value` = fmt_p(dw$p.value), check.names = FALSE),
             digits = 4, caption = "Durbin-Watson test (rows in file order; not a time series)")
Durbin-Watson test (rows in file order; not a time series)
Statistic p value
2.014 0.858

The statistic is 2.014 (p = 0.858). It must not be read as proof that customers are independent. It only looks at correlation between neighbouring rows in the file. Dependence between customers (for example customers of the same branch, product or region, or shared economic conditions) could exist without leaving any trace in file order, and no information on such groupings is available here.

4.5.4 Influential observations and other assumptions

cd  <- cooks.distance(final_lpm)
hv  <- hatvalues(final_lpm)
p_k <- length(coef(final_lpm)); n_tr <- nrow(train)
cd_thr <- 4 / n_tr; hv_thr <- 2 * p_k / n_tr

infl_summary <- data.frame(
  Measure = c("Cook's distance above 4/n", "Leverage above 2p/n", "Largest Cook's distance"),
  Value = c(sum(cd > cd_thr), sum(hv > hv_thr), max(cd)))
knitr::kable(infl_summary, digits = 4, caption = "Influence summary (flagged counts and the maximum)")
Influence summary (flagged counts and the maximum)
Measure Value
Cook’s distance above 4/n 1244.0000
Leverage above 2p/n 1486.0000
Largest Cook’s distance 0.0188 
ggplot(data.frame(i = seq_along(cd), cd = cd), aes(i, cd)) +
  geom_point(alpha = 0.3, size = 0.6, colour = "grey30") +
  geom_hline(yintercept = cd_thr, linetype = "dashed", colour = "#D95F0E") +
  labs(title = "Cook's distance (final LPM, training data)",
       subtitle = "Dashed line: 4/n rule of thumb", x = "Training observation index", y = "Cook's distance")

top_infl <- order(cd, decreasing = TRUE)[1:10]
top_tab <- data.frame(ID = train$ID[top_infl], `Cook's distance` = cd[top_infl], Leverage = hv[top_infl],
                      Default = train$default_next_month[top_infl],
                      `Latest repayment status` = as.character(train$pay0_status[top_infl]),
                      `Bill trend (NT$10,000)` = train$bill_trend_10k[top_infl],
                      `Payment trend (NT$10,000)` = train$payment_trend_10k[top_infl], check.names = FALSE)
knitr::kable(top_tab, digits = 4, caption = "Ten most influential training observations")
Ten most influential training observations
ID Cook’s distance Leverage Default Latest repayment status Bill trend (NT$10,000) Payment trend (NT$10,000)
28719 28717 0.0188 0.1541 0 Paid in full -13.0215 73.0689
6776 6774 0.0129 0.0182 1 Revolving credit, minimum paid 21.4055 -1.6184
21010 21008 0.0045 0.0083 1 Paid in full 0.9334 6.2864
27104 27102 0.0044 0.0093 1 1 month delay -3.7650 7.8540
29924 29922 0.0037 0.0090 1 Revolving credit, minimum paid 43.2428 1.1968
11571 11569 0.0037 0.0067 1 Revolving credit, minimum paid -22.8800 11.1389
13480 13478 0.0026 0.0049 1 No consumption -17.8052 -9.2691
8175 8173 0.0025 0.0054 1 Revolving credit, minimum paid 33.5331 0.6375
27543 27541 0.0024 0.0044 1 Revolving credit, minimum paid 1.9943 5.7834
2910 2908 0.0023 0.0053 1 Revolving credit, minimum paid 33.2040 0.4399

No observation is removed just because it is influential. These are real customers, and no data error was identified. Instead, a sensitivity analysis refits the model without the observations above the 4/n Cook’s distance rule and compares the trend coefficients:

lpm_noinf <- update(final_lpm, data = train[cd <= cd_thr, ])
sens <- rbind(
  data.frame(Model = "Final LPM (all training rows)", Term = trend_terms,
             Estimate = sapply(trend_terms, function(t) get_coef(final_lpm, t)$estimate)),
  data.frame(Model = "Refit without flagged rows", Term = trend_terms,
             Estimate = sapply(trend_terms, function(t) get_coef(lpm_noinf, t)$estimate)))
rownames(sens) <- NULL
knitr::kable(sens, digits = 5, caption = "Sensitivity of the trend coefficients to influential observations")
Sensitivity of the trend coefficients to influential observations
Model Term Estimate
Final LPM (all training rows) bill_trend_10k 0.00254
Final LPM (all training rows) payment_trend_10k 0.00081
Refit without flagged rows bill_trend_10k 0.00081
Refit without flagged rows payment_trend_10k 0.00052

Normality and other assumptions. For a binary outcome the errors cannot be normally distributed, because each customer’s residual takes only two possible values, as the Q-Q plot shows. Normality is not needed for unbiased coefficients. With this many observations, inference on coefficients is reasonable using robust standard errors, but prediction intervals based on normal errors are not meaningful and are not reported.

fit_tr <- fitted(final_lpm)
share_out_train <- mean(fit_tr < 0 | fit_tr > 1)

Another LPM limitation: fitted values outside 0 to 1. 0.30% of fitted values in the training data fall outside the 0 to 1 range (minimum -0.394, maximum 0.882). A probability cannot be below 0 or above 1. This is one reason the logistic model is included as a benchmark.

4.5.5 Joint test of the trend features

jt <- car::linearHypothesis(final_lpm, c("bill_trend_10k = 0", "payment_trend_10k = 0"),
                            vcov. = sandwich::vcovHC(final_lpm, type = "HC3"))
jt_F <- jt[2, "F"]; jt_p <- jt[2, "Pr(>F)"]
knitr::kable(data.frame(Test = "Joint test that both trend coefficients are zero (HC3)",
                        `F statistic` = jt_F, `p value` = fmt_p(jt_p), check.names = FALSE),
             digits = 3, caption = "Robust joint Wald test")
Robust joint Wald test
Test F statistic p value
Joint test that both trend coefficients are zero (HC3) 5.464 0.004

The joint test of both trend coefficients gives F = 5.46 and p = 0.004. This is statistical evidence about association in the training sample. It does not by itself show better prediction. That is assessed on the held-out test data below.

4.6 Logistic regression supplementary benchmark

The logistic model uses the same predictor set as the final LPM (pay0_status, bill_trend_10k, payment_trend_10k, earlier_delays_n, limit_10k, MARRIAGE and pay_latest_10k), so any difference in performance reflects the model type and not the choice of variables. A model without the two trend features is fitted for the same comparison. As a robustness check, BIC-based stepwise selection is also run for the logistic model.

controls_terms <- setdiff(final_terms, trend_terms)

logit_final <- glm(reformulate(final_terms, "default_next_month"), data = train, family = binomial)
logit_base  <- glm(reformulate(controls_terms, "default_next_month"), data = train, family = binomial)

f_lower_g <- reformulate(lower_terms, "default_next_month")
step_logit <- step(glm(f_lower_g, data = train, family = binomial),
                   scope = list(lower = f_lower_g, upper = f_upper),
                   direction = "both", k = log(nrow(train)), trace = 0)
logit_sel_terms <- attr(terms(step_logit), "term.labels")

or_tab <- function(m) {
  co <- summary(m)$coefficients; ci <- confint.default(m)
  data.frame(Term = rownames(co), `Odds ratio` = exp(co[, 1]), `CI low` = exp(ci[, 1]),
             `CI high` = exp(ci[, 2]), `p value` = fmt_p(co[, 4]), check.names = FALSE, row.names = NULL)
}
knitr::kable(or_tab(logit_final), digits = 3,
             caption = "Logistic regression with trend features: odds ratios with Wald 95% confidence intervals (training data)")
Logistic regression with trend features: odds ratios with Wald 95% confidence intervals (training data)
Term Odds ratio CI low CI high p value
(Intercept) 0.164 0.152 0.177 < 0.001
pay0_statusNo consumption 1.356 1.179 1.559 < 0.001
pay0_statusPaid in full 1.508 1.365 1.667 < 0.001
pay0_status1 month delay 2.474 2.231 2.744 < 0.001
pay0_status2 months delay 9.075 8.057 10.222 < 0.001
pay0_status3 or more months delay 6.201 4.843 7.941 < 0.001
bill_trend_10k 1.021 1.008 1.035 0.002
payment_trend_10k 1.005 0.971 1.039 0.785
earlier_delays_n 1.316 1.279 1.355 < 0.001
limit_10k 0.983 0.980 0.986 < 0.001
MARRIAGEMarried 1.215 1.134 1.303 < 0.001
MARRIAGEOther or unknown 1.136 0.835 1.547 0.416
pay_latest_10k 0.898 0.860 0.938 < 0.001 
knitr::kable(vif_show(vif_table(logit_final)), digits = 3, caption = "VIF for the logistic model")
VIF for the logistic model
Term GVIF Df VIF equivalent Flag
pay0_status 1.605 5 1.099 OK
bill_trend_10k 1.107 1 1.107 OK
payment_trend_10k 1.071 1 1.071 OK
earlier_delays_n 1.435 1 1.435 OK
limit_10k 1.208 1 1.208 OK
MARRIAGE 1.019 2 1.009 OK
pay_latest_10k 1.155 1 1.155 OK

Convergence: the logistic model converged. The BIC-based logistic selection chose a different set: pay0_status, bill_trend_10k, payment_trend_10k, earlier_delays_n, limit_10k, pay_latest_10k, MARRIAGE and EDUCATION. The shared predictor set is kept for a like-for-like comparison.

An odds ratio above 1 means higher odds of default for a one-unit increase in that variable, holding the others fixed. Odds are not probabilities, so an odds ratio should not be read as a percentage-point change in default probability.

5 Results and interpretation

5.1 Model validation and comparison

All models are evaluated on the same held-out test set (6,000 customers). Six models are compared:

  1. LPM 1: pre-specified baseline.
  2. LPM 2: pre-specified expanded (baseline plus trends).
  3. LPM 3: stepwise-selected controls without the trends.
  4. LPM 4: final selected LPM with trends (main model).
  5. Logit 1: selected controls without the trends.
  6. Logit 2: selected controls plus the trends.
lpm_sel_base <- lm(reformulate(controls_terms, "default_next_month"), data = train)

models <- list(
  "LPM 1: pre-specified baseline"       = lpm_base,
  "LPM 2: pre-specified expanded"       = lpm_exp,
  "LPM 3: selected, no trends"          = lpm_sel_base,
  "LPM 4: selected + trends (main)"     = final_lpm,
  "Logit 1: selected, no trends"        = logit_base,
  "Logit 2: selected + trends"          = logit_final)

y_te <- test$default_next_month
eval_model <- function(model, label) {
  p <- as.numeric(predict_prob(model, test))
  is_glm <- inherits(model, "glm")
  data.frame(
    Model = label, Family = ifelse(is_glm, "Logistic", "Linear probability"),
    RMSE = sqrt(mean((y_te - p)^2)),
    MAE = mean(abs(y_te - p)),
    `Test R2` = 1 - sum((y_te - p)^2) / sum((y_te - y_train_mean)^2),
    Brier = mean((y_te - p)^2),
    `Log loss` = if (is_glm) mean(log_loss_vec(y_te, p)) else NA_real_,
    AUC = auc_rank(y_te, p),
    `RMSE (clipped to 0 to 1)` = sqrt(mean((y_te - pmin(pmax(p, 0), 1))^2)),
    `Share outside 0 to 1` = mean(p < 0 | p > 1),
    check.names = FALSE)
}
metrics <- do.call(rbind, Map(eval_model, models, names(models)))
rownames(metrics) <- NULL
knitr::kable(metrics, digits = 4, caption = "Test-set performance of all models")
Test-set performance of all models
Model Family RMSE MAE Test R2 Brier Log loss AUC RMSE (clipped to 0 to 1) Share outside 0 to 1
LPM 1: pre-specified baseline Linear probability 0.3721 0.2765 0.1963 0.1384 0.7653 0.3720 0.0032
LPM 2: pre-specified expanded Linear probability 0.3718 0.2764 0.1976 0.1382 0.7680 0.3717 0.0037
LPM 3: selected, no trends Linear probability 0.3721 0.2768 0.1964 0.1384 0.7648 0.3720 0.0025
LPM 4: selected + trends (main) Linear probability 0.3718 0.2766 0.1975 0.1382 0.7675 0.3718 0.0030
Logit 1: selected, no trends Logistic 0.3724 0.2764 0.1950 0.1387 0.4415 0.7661 0.3724 0.0000
Logit 2: selected + trends Logistic 0.3722 0.2762 0.1959 0.1385 0.4406 0.7690 0.3722 0.0000

How to read the metrics.

  • Test R-squared is 1 - SSE_test / SST, where the benchmark is a constant prediction equal to the training default rate. It measures improvement over that naive forecast and can be negative.
  • RMSE and MAE measure the typical size of the error between predicted probability and the 0/1 outcome. RMSE punishes large errors more than MAE.
  • Brier score is the mean squared error of the predicted probability. It equals RMSE squared, so it carries the same information as RMSE for probabilistic predictions.
  • Log loss is defined only for predictions strictly between 0 and 1, so it is reported for the logistic models only. It penalises confident wrong predictions much more heavily than RMSE does, so RMSE/Brier (a quadratic score) and log loss (a logarithmic score) emphasise different aspects of probability quality and should not be treated as interchangeable.
  • AUC measures how well a model ranks defaulters above non-defaulters, regardless of calibration.
  • The linear probability and logistic models are compared on the metrics they share (RMSE, MAE, Brier, AUC). The logistic model is a benchmark and is not claimed to satisfy the course’s continuous-outcome requirement.

5.1.1 Does adding the trend features improve prediction?

Differences are evaluated customer by customer on the test set: for each customer the squared error (or log loss) of the model without trends minus that of the model with trends, with a 95% confidence interval on the mean difference. A positive difference means the trend features reduced the error. A paired bootstrap interval is given for AUC.

pair_compare <- function(m_base, m_exp, loss = c("sq", "log")) {
  loss <- match.arg(loss)
  pb <- as.numeric(predict_prob(m_base, test)); pe <- as.numeric(predict_prob(m_exp, test))
  if (loss == "sq") { lb <- (y_te - pb)^2; le <- (y_te - pe)^2
  } else { lb <- log_loss_vec(y_te, pb); le <- log_loss_vec(y_te, pe) }
  d <- lb - le; tt <- t.test(d)
  rel <- if (loss == "sq") 1 - sqrt(mean(le)) / sqrt(mean(lb)) else 1 - mean(le) / mean(lb)
  list(mean_diff = mean(d), lo = tt$conf.int[1], hi = tt$conf.int[2], p = tt$p.value,
       rel_change = rel, auc_base = auc_rank(y_te, pb), auc_exp = auc_rank(y_te, pe), pb = pb, pe = pe)
}
boot_auc_diff <- function(pb, pe, B = 500) {
  set.seed(params$seed + 1)
  n <- length(y_te)
  d <- replicate(B, { i <- sample.int(n, replace = TRUE); auc_rank(y_te[i], pe[i]) - auc_rank(y_te[i], pb[i]) })
  quantile(d, c(0.025, 0.975))
}

pc_pre  <- pair_compare(lpm_base, lpm_exp, "sq")
pc_sel  <- pair_compare(lpm_sel_base, final_lpm, "sq")
pc_logi <- pair_compare(logit_base, logit_final, "log")
auc_ci_pre  <- boot_auc_diff(pc_pre$pb,  pc_pre$pe)
auc_ci_sel  <- boot_auc_diff(pc_sel$pb,  pc_sel$pe)
auc_ci_logi <- boot_auc_diff(pc_logi$pb, pc_logi$pe)

practical <- 0.01   # analyst-chosen benchmark: 1% relative reduction in error counts as practically material

verdict <- function(pc, what) {
  if (pc$lo > 0) {
    if (pc$rel_change >= practical) paste0("a statistically detectable and practically material improvement (", what, " fell by ", fmt_pct(pc$rel_change, 2), ")")
    else paste0("a statistically detectable but small improvement (", what, " fell by only ", fmt_pct(pc$rel_change, 2), ", below the ", fmt_pct(practical, 0), " benchmark)")
  } else if (pc$hi < 0) {
    paste0("a worse out-of-sample result (", what, " rose by ", fmt_pct(-pc$rel_change, 2), ")")
  } else {
    paste0("no reliable difference (the interval for the mean paired difference includes zero; ", what, " changed by ", fmt_pct(pc$rel_change, 2), ")")
  }
}

pair_tab <- data.frame(
  Comparison = c("LPM 1 to LPM 2 (pre-specified): squared error", "LPM 3 to LPM 4 (selected): squared error",
                 "Logit 1 to Logit 2 (selected): log loss"),
  `Mean reduction in loss` = c(pc_pre$mean_diff, pc_sel$mean_diff, pc_logi$mean_diff),
  `CI low` = c(pc_pre$lo, pc_sel$lo, pc_logi$lo), `CI high` = c(pc_pre$hi, pc_sel$hi, pc_logi$hi),
  `p value` = fmt_p(c(pc_pre$p, pc_sel$p, pc_logi$p)),
  `Relative reduction (RMSE or log loss)` = fmt_pct(c(pc_pre$rel_change, pc_sel$rel_change, pc_logi$rel_change), 2),
  `AUC without trends` = c(pc_pre$auc_base, pc_sel$auc_base, pc_logi$auc_base),
  `AUC with trends` = c(pc_pre$auc_exp, pc_sel$auc_exp, pc_logi$auc_exp),
  `AUC difference, bootstrap 95% CI` = c(sprintf("[%.4f, %.4f]", auc_ci_pre[1], auc_ci_pre[2]),
                                         sprintf("[%.4f, %.4f]", auc_ci_sel[1], auc_ci_sel[2]),
                                         sprintf("[%.4f, %.4f]", auc_ci_logi[1], auc_ci_logi[2])),
  check.names = FALSE)
knitr::kable(pair_tab, digits = 5, caption = "Paired test-set comparison: with versus without the trend features")
Paired test-set comparison: with versus without the trend features
Comparison Mean reduction in loss CI low CI high p value Relative reduction (RMSE or log loss) AUC without trends AUC with trends AUC difference, bootstrap 95% CI
LPM 1 to LPM 2 (pre-specified): squared error 0.00022 0.00008 0.00037 0.003 0.08% 0.7653 0.7680 [0.0010, 0.0042]
LPM 3 to LPM 4 (selected): squared error 0.00020 0.00006 0.00034 0.005 0.07% 0.7648 0.7674 [0.0009, 0.0044]
Logit 1 to Logit 2 (selected): log loss 0.00083 0.00026 0.00140 0.004 0.19% 0.7661 0.7690 [0.0012, 0.0047]

The 1% relative-error benchmark for “practically material” is an analyst-chosen illustration, not a universal standard. The reader can apply a different benchmark to the table above.

5.1.2 Calibration: predicted probabilities versus observed outcomes

Test customers are ranked by predicted probability and split into ten equal groups. For a well-calibrated model, the average predicted probability in each group is close to the observed default rate.

calib <- function(p, name, bins = 10) {
  b <- cut(rank(p, ties.method = "first"), breaks = seq(0, length(p), length.out = bins + 1),
           labels = FALSE, include.lowest = TRUE)
  d <- data.frame(bin = b, observed = y_te, predicted = p)
  out <- aggregate(cbind(observed, predicted) ~ bin, data = d, FUN = mean)
  out$n <- as.vector(table(b)); out$Model <- name; out
}
p_lpm   <- as.numeric(predict_prob(final_lpm, test))
p_logit <- as.numeric(predict_prob(logit_final, test))
cal_df  <- rbind(calib(p_lpm, "LPM 4 (main)"), calib(p_logit, "Logit 2"))

cal_show <- reshape(cal_df[, c("Model", "bin", "predicted", "observed")], idvar = "bin", timevar = "Model", direction = "wide")
knitr::kable(cal_show, digits = 3, row.names = FALSE,
             caption = "Mean predicted probability and observed default rate by test-set decile of predicted risk")
Mean predicted probability and observed default rate by test-set decile of predicted risk
bin predicted.LPM 4 (main) observed.LPM 4 (main) predicted.Logit 2 observed.Logit 2
1 0.065 0.057 0.078 0.055
2 0.104 0.080 0.107 0.073
3 0.121 0.095 0.121 0.097
4 0.131 0.130 0.131 0.118
5 0.144 0.132 0.141 0.150
6 0.158 0.152 0.156 0.157
7 0.180 0.180 0.174 0.185
8 0.241 0.278 0.224 0.265
9 0.405 0.437 0.406 0.445
10 0.724 0.672 0.732 0.667 
ggplot(cal_df, aes(predicted, observed, colour = Model)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_line() + geom_point(size = 2) +
  scale_x_continuous(labels = scales::percent) + scale_y_continuous(labels = scales::percent) +
  labs(title = "Calibration on the test set",
       subtitle = "Dashed line: perfect calibration",
       x = "Mean predicted probability of default", y = "Observed default rate")

5.1.3 Classification view of the logistic benchmark (optional)

A probability model can be turned into a yes/no flag by choosing a threshold. The threshold here is chosen on the training data only, by maximising sensitivity plus specificity minus one (Youden’s J). A 0.5 threshold would flag very few customers because the default rate is low. The threshold is a business choice. This one only illustrates the trade-off.

best_thr <- function(y, p) {
  grid <- seq(0.05, 0.95, by = 0.01)
  j <- sapply(grid, function(t) {
    pred <- p >= t
    sum(pred & y == 1) / sum(y == 1) + sum(!pred & y == 0) / sum(y == 0) - 1
  })
  grid[which.max(j)]
}
class_metrics <- function(y, p, thr) {
  pred <- p >= thr
  TP <- sum(pred & y == 1); FP <- sum(pred & y == 0); TN <- sum(!pred & y == 0); FN <- sum(!pred & y == 1)
  c(Threshold = thr, Sensitivity = TP / (TP + FN), Specificity = TN / (TN + FP),
    Precision = ifelse(TP + FP > 0, TP / (TP + FP), NA), Accuracy = (TP + TN) / length(y),
    `Balanced accuracy` = (TP / (TP + FN) + TN / (TN + FP)) / 2)
}
p_logit_base <- as.numeric(predict_prob(logit_base, test))
thr_base  <- best_thr(train$default_next_month, fitted(logit_base))
thr_final <- best_thr(train$default_next_month, fitted(logit_final))
cls_tab <- rbind(
  data.frame(Model = "Logit 1: selected, no trends", t(class_metrics(y_te, p_logit_base, thr_base)), check.names = FALSE),
  data.frame(Model = "Logit 2: selected + trends",   t(class_metrics(y_te, p_logit, thr_final)), check.names = FALSE))
knitr::kable(cls_tab, digits = 3, caption = "Test-set classification metrics at the training-selected threshold")
Test-set classification metrics at the training-selected threshold
Model Threshold Sensitivity Specificity Precision Accuracy Balanced accuracy
Logit 1: selected, no trends 0.19 0.626 0.788 0.456 0.752 0.707
Logit 2: selected + trends 0.19 0.626 0.788 0.456 0.752 0.707 
pred_class <- factor(ifelse(p_logit >= thr_final, "Predicted default", "Predicted no default"),
                     levels = c("Predicted no default", "Predicted default"))
obs_class  <- factor(ifelse(y_te == 1, "Observed default", "Observed no default"),
                     levels = c("Observed no default", "Observed default"))
knitr::kable(as.data.frame.matrix(table(pred_class, obs_class)), caption = "Confusion matrix for Logit 2 on the test set")
Confusion matrix for Logit 2 on the test set
Observed no default Observed default
Predicted no default 3681 496
Predicted default 992 831

5.2 Dynamic results interpretation

5.2.1 Sample, selected model and fit

The analysis uses 30,000 customers, split into 24,000 for training and 6,000 for testing. The observed default rate is 22.12% overall (22.12% in training and 22.12% in testing).

The final linear probability model contains: latest repayment status (PAY_0), bill trend (NT$10,000), payment trend (NT$10,000), number of earlier delayed months (PAY_2 to PAY_6), credit limit (NT$10,000), marital status and latest payment (NT$10,000). On the training data its R-squared is 0.2072 and adjusted R-squared 0.2068. On the test data, the main LPM has RMSE 0.3718, MAE 0.2766, test R-squared 0.1975 and AUC 0.7675. The logistic benchmark with the same predictors has Brier score 0.1385, log loss 0.4406 and AUC 0.7690.

5.2.2 The trend features in the final LPM (HC3 robust inference)

Bill trend. The estimated coefficient is 0.00254 (robust SE 0.00078, 95% CI 0.00101 to 0.00407, p = 0.001). In plain terms: for customers with the same values of the other model variables, an increase of NT$10,000 in the bill trend (recent bills rising relative to earlier bills) is associated with a default probability that is 0.254 percentage points higher. A one standard deviation difference in the bill trend (NT$29,866) corresponds to about 0.759 percentage points. The coefficient is statistically significant at the 5% level.

Payment trend. The estimated coefficient is 0.00081 (robust SE 0.00161, 95% CI -0.00235 to 0.00397, p = 0.617). An increase of NT$10,000 in the payment trend is associated with a default probability that is 0.081 percentage points higher. Because declining payments are negative trend values, a NT$10,000 decline in recent payments relative to earlier payments is associated with a probability that is 0.081 percentage points lower. One standard deviation of the payment trend (NT$14,107) corresponds to about 0.114 percentage points. The coefficient is not statistically significant at the 5% level. The data did not provide sufficient evidence of an association under this model. This is not proof that the effect is exactly zero.

Joint test. The robust joint test of both trend coefficients gives F = 5.46, p = 0.004, so the two trends together are statistically significant at the 5% level.

Logistic benchmark. The odds ratio for a NT$10,000 increase in the bill trend is 1.021 (95% CI 1.008 to 1.035), and for the payment trend it is 1.005 (95% CI 0.971 to 1.039). An odds ratio of 1 means no association.

5.2.3 Statistical significance versus predictive value

A small p-value shows that an association is unlikely to be exactly zero in the training sample. It does not show that the feature improves prediction for new customers. The answer to the business question therefore rests on the test-set comparison:

  • Linear probability models (selected specification). Adding the trends gives a statistically detectable but small improvement (RMSE fell by only 0.07%, below the 1% benchmark). Test RMSE moved from 0.37205 to 0.37179 and AUC from 0.7648 to 0.7675 (bootstrap 95% interval for the AUC change: 0.0009 to 0.0044).
  • Linear probability models (pre-specified baseline and expanded). Adding the trends gives a statistically detectable but small improvement (RMSE fell by only 0.08%, below the 1% benchmark).
  • Logistic models. Adding the trends gives a statistically detectable but small improvement (log loss fell by only 0.19%, below the 1% benchmark). Test log loss moved from 0.44147 to 0.44064 and AUC from 0.7661 to 0.7690.

Across all six models, the lowest test RMSE is LPM 2: pre-specified expanded, the highest AUC is Logit 2: selected + trends, and the lowest log loss among the logistic models is Logit 2: selected + trends.

5.2.4 Charts of the main results

# Chart 1: size of the two trend effects (percentage points of default probability per NT$10,000), robust 95% CIs
eff <- data.frame(
  Feature = c("Bill trend (+NT$10,000)", "Payment trend (+NT$10,000)"),
  est = c(bt$estimate, pt$estimate) * 100,
  lo  = c(bt$lo, pt$lo) * 100,
  hi  = c(bt$hi, pt$hi) * 100)
ggplot(eff, aes(est, Feature)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_segment(aes(x = lo, xend = hi, yend = Feature), linewidth = 1) +
  geom_point(size = 3, colour = "#D95F0E") +
  labs(title = "Estimated association of each trend with next-month default",
       subtitle = "Final linear probability model; points: estimates; bars: HC3 robust 95% confidence intervals",
       x = "Change in default probability (percentage points)", y = NULL)

# Chart 2: test-set error of every model (lower is better)
ggplot(metrics, aes(reorder(Model, -RMSE), RMSE)) +
  geom_col(fill = col_main, width = 0.7) +
  geom_text(aes(label = fmt_num(RMSE, 4)), hjust = -0.1, size = 3.3) +
  coord_flip(ylim = c(min(metrics$RMSE) * 0.9, max(metrics$RMSE) * 1.03)) +
  labs(title = "Test-set RMSE by model",
       subtitle = "RMSE = typical error of the predicted default probability (probability units); lower is better. Axis does not start at zero.",
       x = NULL, y = "Test RMSE (probability units)")

# Chart 3: how well each model ranks defaulters above non-defaulters
ggplot(metrics, aes(reorder(Model, AUC), AUC)) +
  geom_col(fill = "#41AB5D", width = 0.7) +
  geom_text(aes(label = fmt_num(AUC, 3)), hjust = -0.1, size = 3.3) +
  coord_flip(ylim = c(0.5, max(metrics$AUC) * 1.03)) +
  labs(title = "Test-set AUC by model",
       subtitle = "AUC: chance that a random defaulter gets a higher predicted risk than a random non-defaulter; 0.5 = no skill. Axis starts at 0.5.",
       x = NULL, y = "Test AUC")

The first chart shows the direction and size of each trend effect, and whether the interval covers zero (no detectable association). The other two charts compare all six models on the same test customers: RMSE measures error size and AUC measures ranking ability. Small gaps between bars should be judged with the paired comparison table above, not by eye.

5.2.5 Diagnostic summary

Diagnostic Result for the final LPM
Multicollinearity Largest VIF equivalent 1.61 (threshold 5); no variable removed
Heteroscedasticity Breusch-Pagan p < 0.001; HC3 robust standard errors used for all inference
Autocorrelation Durbin-Watson 2.014 (p = 0.858); not meaningful for unordered customers
Influence 1244 training rows above the 4/n Cook’s distance rule; sensitivity refit reported
Fitted values outside 0 to 1 0.30% of training fitted values; 0.30% of test predictions
Normality of errors Not possible for a binary outcome; inference relies on robust standard errors and large samples

6 Conclusion and recommendations

6.1 Business implications and recommendations

Main finding. Recent bill and payment trends added a detectable but small amount of predictive information beyond the latest repayment status and the other customer characteristics. The improvement on the held-out customers was below the 1% benchmark used here, so its practical value is limited. The latest repayment status remains the dominant signal: observed default rates range from 12.8% to 71.9% across its categories.

Recommendations.

  1. Prioritise repayment status and earlier delinquency history in early-warning monitoring. The trend features are not a reliable substitute for them in this data.
  2. If trend indicators are still of operational interest, treat them as supplementary and monitor whether their value changes on more recent data.
  3. Because the model describes association and not cause, do not assume that changing a customer’s bill or payment pattern would change their default risk. Any intervention should be piloted and evaluated separately, for example with a controlled trial.
  4. Choose an operating threshold from the business costs of missed defaults versus unnecessary reviews. The threshold shown above is only an illustration.
  5. Review the use of age and marital status with compliance before any operational use, since regulations on these characteristics vary by jurisdiction.

6.2 Limitations

  • Observational and predictive. The results describe associations. They do not show that bill or payment changes cause default.
  • Historical, single-market data. The data come from Taiwan and describe an earlier period. Customer behaviour, products and the economy may differ now or elsewhere, so the findings may not generalise.
  • Random split, not a time split. The test set is drawn from the same period and population as the training set. Performance on future customers may be lower, and no temporal validation is possible with this single snapshot.
  • Linear probability model. It has non-constant error variance, non-normal errors and 0.30% of test predictions outside the 0 to 1 range. Robust standard errors address inference only, not the range problem or non-linearity. The logistic benchmark is a better fit for a binary outcome in principle.
  • Dependence between customers. The Durbin-Watson test cannot rule out clustering of customers. Standard errors assume independent customers.
  • Constructed features. Trends compare two three-month averages of six months of data. They are crude, can be unstable for customers with little activity, and are expressed in a single currency scale without adjusting for credit limit or inflation.
  • Overlap with repayment status. The repayment status variables are themselves partly determined by earlier bills and payments, so the trends can add little once repayment status is known.
  • Data-driven selection. Stepwise selection can overfit, which is why the held-out test set is used for all predictive conclusions.
  • Pooled categories. Undocumented education and marital codes were pooled, and delays of 3 or more months were grouped. Other reasonable choices might change coefficients slightly.
  • Excluded characteristic. SEX was deliberately excluded. Fairness across groups was not assessed.

6.3 Final conclusion

Using 30,000 customers (24,000 for training and 6,000 for testing; observed default rate 22.12%), a linear probability model with the predictors latest repayment status (PAY_0), bill trend (NT$10,000), payment trend (NT$10,000), number of earlier delayed months (PAY_2 to PAY_6), credit limit (NT$10,000), marital status and latest payment (NT$10,000) was selected on the training data and checked with OLS diagnostics. Adding the bill trend and payment trend gave a statistically detectable but small improvement (RMSE fell by only 0.07%, below the 1% benchmark) in the linear probability model, and a statistically detectable but small improvement (log loss fell by only 0.19%, below the 1% benchmark) in the logistic benchmark. The bill trend coefficient was 0.254 percentage points per NT$10,000 (robust 95% CI 0.101 to 0.407), and the payment trend coefficient was 0.081 percentage points per NT$10,000 (robust 95% CI -0.235 to 0.397). These are associations in historical Taiwanese data and should not be read as causal effects. Any use for present-day customers or other markets would require fresh validation.

7 Appendix

7.1 Appendix A: Pre-specified models versus the selected model

knitr::kable(metrics[, c("Model", "RMSE", "MAE", "Test R2", "AUC", "RMSE (clipped to 0 to 1)", "Share outside 0 to 1")],
             digits = 4, caption = "Test-set metrics for the linear probability models and the logistic benchmark")
Test-set metrics for the linear probability models and the logistic benchmark
Model RMSE MAE Test R2 AUC RMSE (clipped to 0 to 1) Share outside 0 to 1
LPM 1: pre-specified baseline 0.3721 0.2765 0.1963 0.7653 0.3720 0.0032
LPM 2: pre-specified expanded 0.3718 0.2764 0.1976 0.7680 0.3717 0.0037
LPM 3: selected, no trends 0.3721 0.2768 0.1964 0.7648 0.3720 0.0025
LPM 4: selected + trends (main) 0.3718 0.2766 0.1975 0.7675 0.3718 0.0030
Logit 1: selected, no trends 0.3724 0.2764 0.1950 0.7661 0.3724 0.0000
Logit 2: selected + trends 0.3722 0.2762 0.1959 0.7690 0.3722 0.0000

7.2 Appendix B: Final model formulas

formulas <- data.frame(
  Model = names(models),
  Formula = vapply(models, function(m) paste(deparse(formula(m), width.cutoff = 500), collapse = " "), character(1)),
  row.names = NULL)
knitr::kable(formulas, caption = "Formulas of all models compared")
Formulas of all models compared
Model Formula
LPM 1: pre-specified baseline default_next_month ~ pay0_status + limit_10k + AGE + EDUCATION + MARRIAGE + earlier_delays_n
LPM 2: pre-specified expanded default_next_month ~ pay0_status + limit_10k + AGE + EDUCATION + MARRIAGE + earlier_delays_n + bill_trend_10k + payment_trend_10k
LPM 3: selected, no trends default_next_month ~ pay0_status + earlier_delays_n + limit_10k + MARRIAGE + pay_latest_10k
LPM 4: selected + trends (main) default_next_month ~ pay0_status + bill_trend_10k + payment_trend_10k + earlier_delays_n + limit_10k + MARRIAGE + pay_latest_10k
Logit 1: selected, no trends default_next_month ~ pay0_status + earlier_delays_n + limit_10k + MARRIAGE + pay_latest_10k
Logit 2: selected + trends default_next_month ~ pay0_status + bill_trend_10k + payment_trend_10k + earlier_delays_n + limit_10k + MARRIAGE + pay_latest_10k