Decision Trees (R) — rpart

CART with cost-complexity pruning, the reference implementation

rpart ("recursive partitioning") is the canonical R implementation of CART — the same Gini splitting and cost-complexity pruning as the from-scratch Python notebook. It is the natural cross-check; we also add the traditional-econometrics benchmark (glm logistic regression) that an analyst would reach for on a binary outcome.

The data and the goals. The Taiwan credit-card default study (Yeh & Lien, 2009) — 30,000 cardholders of a Taiwanese bank, April–September 2005. Goal: binary classification — predict default on next month's payment from credit limit, demographics, and six months of repayment status (PAY_1…6, months in arrears), bill amounts and payment amounts; the credit-scoring task of ranking clients by default risk (22% defaulted). For regression: California housing — predict a neighbourhood's median house value ($100,000s) from income, house age, occupancy and location. All performance is reported out of sample on a 30% held-out test set.

(Our headline classification metric is ROC-AUC = the probability the model scores a random defaulter above a random non-defaulter (0.5 = chance, 1 = perfect); threshold-free and robust to the 22% class imbalance, unlike accuracy. Equivalently it is the area under the ROC curve (TPR vs FPR, plotted in §5). The auc() helper below computes it as the Mann–Whitney rank statistic $\text{AUC}=\big(\sum_{i\in\text{pos}} r_i - n_1(n_1{+}1)/2\big)/(n_1 n_0)$ from the ascending ranks $r_i$ of the predicted scores — no threshold sweep needed.)

options(repr.plot.width=9, repr.plot.height=5.5)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
for(p in c("rpart.plot")) if(!requireNamespace(p,quietly=TRUE)) install.packages(p, repos="https://cloud.r-project.org", quiet=TRUE)
suppressMessages({library(rpart); library(rpart.plot)})
d <- read.csv("credit_default.csv"); d$default <- factor(d$default, labels=c("no","default"))
set.seed(0); i <- sample(nrow(d), 0.7*nrow(d)); tr <- d[i,]; te <- d[-i,]
auc <- function(y,p){ r<-rank(p); n1<-sum(y==1); n0<-sum(y==0); (sum(r[y==1])-n1*(n1+1)/2)/(n1*n0) }
yte <- as.integer(te$default=="default")
cat(sprintf("credit default: %d train / %d test (held out); overall default rate %.1f%%\n", nrow(tr), nrow(te), 100*mean(d$default=="default")))

Warning message:
"package 'rpart.plot' was built under R version 4.6.1"

credit default: 21000 train / 9000 test (held out); overall default rate 22.1%

1. The data up close — what we are predicting

The outcome is default (did the client miss next month's payment?). The strongest raw signal is recent delinquency PAY_1 (months in arrears last month; $\le 0$ = paid/revolving). Default risk climbs steeply with it, and is higher at lower credit limits — the patterns the tree will formalise.

dd <- read.csv("credit_default.csv")
par(mfrow=c(1,2), mar=c(4,4,3,1))
g <- tapply(dd$default, dd$PAY_1, mean); g <- g[as.integer(names(g)) %in% -2:8]
barplot(g, col=ifelse(as.integer(names(g))>=1,"#c53030","#2b6cb0"), main="Default rate by recent delinquency (PAY_1)",
        xlab="PAY_1 (months in arrears)", ylab="default rate"); abline(h=mean(dd$default), lty=2)
dd$limq <- cut(dd$LIMIT_BAL, quantile(dd$LIMIT_BAL, 0:4/4), include.lowest=TRUE, labels=c("Q1 low","Q2","Q3","Q4 high"))
barplot(tapply(dd$default, dd$limq, mean), col="#dd6b20", main="Default rate by credit-limit quartile",
        xlab="credit limit quartile", ylab="default rate"); abline(h=mean(dd$default), lty=2)
par(mfrow=c(1,1))
r0 <- mean(dd$default[dd$PAY_1<=0]); r2 <- mean(dd$default[dd$PAY_1>=2])
cat(sprintf("CONCLUSION (raw data): paid-on-time clients default %.0f%%; those 2+ months behind default %.0f%% -- a %.1fx jump.\n", 100*r0, 100*r2, r2/r0))
cat("Recent repayment behaviour dominates demographics as a risk signal.\n")

CONCLUSION (raw data): paid-on-time clients default 14%; those 2+ months behind default 70% -- a 5.0x jump.

Recent repayment behaviour dominates demographics as a risk signal.

2. Grow and read a tree

rpart(default ~ ., method="class") fits a classification tree; rpart.plot draws it. The first split is PAY_1, matching the from-scratch Python tree.

fit3 <- rpart(default ~ ., data=tr, method="class", control=rpart.control(maxdepth=3, cp=0))
rpart.plot(fit3, type=2, extra=104, box.palette="RdBu", main="rpart CART tree (depth 3) — credit default")
cat("root split variable:", rownames(fit3$splits)[1], "\n")

root split variable: PAY_1 

3. Cost-complexity pruning via cp

Grow a large tree (cp=0), then let rpart's internal 10-fold cross-validation (CV) choose the pruning strength. Cross-validation estimates how well a model will do on data it has not seen, without spending the test set to find out: split the training rows into k equal folds — here $k=10$ — fit the model on $k-1$ of them, measure the error on the fold held out, and repeat until every fold has been the held-out one. Averaging those ten errors gives an out-of-sample estimate for each candidate pruning strength, and the strength with the lowest average wins. The cost is that the model is refit ten times; the benefit is that the choice is made without ever touching the test data, which is what keeps the final number honest. (The Model Selection subsection later in this arc sets cross-validation beside AIC/BIC and Bayesian LOO as four routes to the same question.)

The plot below shows the cross-validated error (± 1 SE) against tree size, with the CV-optimal size marked — the standard pruning curve. A subtlety (previewing the Evaluation subsection): rpart's CV optimises misclassification, not AUC. Under 22% class imbalance the misclassification-optimal tree is small — accuracy prefers it, but AUC (a ranking metric) rewards the larger tree. We report both, out of sample.

set.seed(1)
big <- rpart(default ~ ., data=tr, method="class", control=rpart.control(cp=0, minbucket=5, xval=10))
ct <- big$cptable; sz <- ct[,"nsplit"]+1; xerr <- ct[,"xerror"]; xstd <- ct[,"xstd"]   # clearer than plotcp() on a huge tree
par(mfrow=c(1,1), mar=c(4,4,3,1))
plot(sz, xerr, type="o", pch=19, col="#2b6cb0", log="x", xlab="tree size (number of leaves, log scale)",
     ylab="cross-validated error (relative)", main="Cost-complexity pruning: CV error vs tree size")
arrows(sz, xerr-xstd, sz, xerr+xstd, angle=90, code=3, length=0.02, col="grey70")
m0 <- which.min(xerr); abline(h=xerr[m0]+xstd[m0], lty=2, col="#c53030"); abline(v=sz[m0], lty=3, col="#2f855a", lwd=2)
legend("topright", c("CV error +/- 1 SE","min + 1 SE","CV-optimal size"), col=c("#2b6cb0","#c53030","#2f855a"), lty=c(1,2,3), pch=c(19,NA,NA), bty="n")
best <- ct[m0, "CP"]
pruned <- prune(big, cp=best)
p.big <- predict(big, te, type="prob")[,2]; p.pru <- predict(pruned, te, type="prob")[,2]
cat(sprintf("CV-optimal cp = %.2e\n", best))
cat(sprintf("pruned tree: %d leaves | OUT-OF-SAMPLE AUC %.4f, accuracy %.4f\n", sum(pruned$frame$var=="<leaf>"), auc(yte,p.pru), mean((p.pru>0.5)==yte)))
cat(sprintf("full   tree: %d leaves | OUT-OF-SAMPLE AUC %.4f, accuracy %.4f\n", sum(big$frame$var=="<leaf>"), auc(yte,p.big), mean((p.big>0.5)==yte)))
cat("Pruning gives a far simpler, more accurate tree; the larger tree edges it on AUC -- the misclassification-vs-ranking tension.\n")

CV-optimal cp = 1.29e-03

pruned tree: 10 leaves | OUT-OF-SAMPLE AUC 0.6881, accuracy 0.8202

full   tree: 916 leaves | OUT-OF-SAMPLE AUC 0.7043, accuracy 0.7514

Pruning gives a far simpler, more accurate tree; the larger tree edges it on AUC -- the misclassification-vs-ranking tension.

4. Variable importance

rpart accumulates each variable's contribution to impurity reduction. As in Python, the recent-repayment-status variables dominate.

vi <- head(sort(pruned$variable.importance, decreasing=TRUE), 8)
par(mfrow=c(1,1), mar=c(4,7,3,1)); barplot(rev(vi), horiz=TRUE, las=1, col="#2b6cb0", main="rpart variable importance (pruned tree)", xlab="impurity reduction")

5. A traditional-econometrics benchmark — logistic regression (glm)

The classical approach to a binary outcome is logistic regression: log-odds of default linear in the features, fit by glm(..., family=binomial). It trades the tree's flexibility for interpretable coefficients — each $e^{\beta}$ is an odds ratio. We compare its out-of-sample AUC to the tree and read off the odds ratios.

lg <- glm(default ~ ., data=tr, family=binomial)
p.lg <- predict(lg, te, type="response")
cat(sprintf("OUT-OF-SAMPLE AUC:  logistic regression %.4f  |  pruned CART %.4f\n", auc(yte,p.lg), auc(yte,p.pru)))
or <- exp(coef(lg)[c("PAY_1","LIMIT_BAL","AGE")])
cat("\nlogit odds ratios:\n"); print(round(or,4))
cat(sprintf("\nINTERPRETATION: each extra month in arrears (PAY_1) multiplies the odds of default by ~%.1f.\n", or["PAY_1"]))
if (auc(yte,p.pru) > auc(yte,p.lg)) {
  cat(sprintf("Here the tree edges the logit (%.4f vs %.4f) by capturing nonlinearities the linear log-odds miss.\n",
              auc(yte,p.pru), auc(yte,p.lg)))
} else {
  cat(sprintf("Here the LOGIT wins (%.4f vs %.4f for the pruned tree) -- and that is worth pausing on, because\n",
              auc(yte,p.lg), auc(yte,p.pru)))
  cat("the Python notebook reaches the opposite verdict on the SAME data (tree 0.744, logit 0.715). Nothing\n")
  cat("about the data changed; what changed is how hard the tree was pruned. rpart stops at its 1-SE default,\n")
  cat("a smaller tree than the cross-validated optimum the Python notebook selects, and a few leaves are worth\n")
  cat("more than the gap to the logit. Read \"the tree beats the regression\" as a statement about a TUNED\n")
  cat("tree, not about trees; on this dataset the ordering flips with the pruning rule alone.\n")
}
# probability of default vs months-in-arrears: empirical proportions vs PROBIT vs tree (functional form)
pay <- tr$PAY_1; yb <- as.integer(tr$default=="default")
emp <- tapply(yb, pay, mean); vals <- as.numeric(names(emp)); nobs <- as.numeric(table(pay))
pr <- glm(yb ~ pay, family=binomial(link="probit"))
grid <- seq(min(vals), max(vals), length=300); pp <- predict(pr, data.frame(pay=grid), type="response")
t1 <- rpart(default ~ PAY_1, data=tr, method="class", control=rpart.control(minbucket=300))
ptree <- predict(t1, data.frame(PAY_1=grid), type="prob")[,2]
par(mfrow=c(1,1), mar=c(4,4,3,1))
plot(vals, emp, cex=pmin(3.5, nobs/2500), pch=19, xlab="PAY_1 (months in arrears)", ylab="probability of default",
     main="Default vs delinquency: empirical, probit, tree", ylim=c(0, max(emp)*1.05))
lines(grid, pp, col="#dd6b20", lwd=2); lines(grid, ptree, col="#2b6cb0", lwd=2)
legend("topleft", c("empirical proportion (size ~ n)","probit","tree"), col=c("black","#dd6b20","#2b6cb0"), pch=c(19,NA,NA), lty=c(NA,1,1), bty="n")
cat("Empirical default is flat for paid-up/revolving clients then jumps once in arrears; probit smooths a monotone S-curve, the tree steps to the kink.\n")

Warning message:
"glm.fit: fitted probabilities numerically 0 or 1 occurred"

OUT-OF-SAMPLE AUC:  logistic regression 0.7120  |  pruned CART 0.6881

logit odds ratios:

    PAY_1 LIMIT_BAL       AGE 
   1.8147    1.0000    1.0063 

INTERPRETATION: each extra month in arrears (PAY_1) multiplies the odds of default by ~1.8.

Here the LOGIT wins (0.7120 vs 0.6881 for the pruned tree) -- and that is worth pausing on, because
the Python notebook reaches the opposite verdict on the SAME data (tree 0.744, logit 0.715). Nothing
about the data changed; what changed is how hard the tree was pruned. rpart stops at its 1-SE default,
a smaller tree than the cross-validated optimum the Python notebook selects, and a few leaves are worth
more than the gap to the logit. Read "the tree beats the regression" as a statement about a TUNED
tree, not about trees; on this dataset the ordering flips with the pruning rule alone.

Empirical default is flat for paid-up/revolving clients then jumps once in arrears; probit smooths a monotone S-curve, the tree steps to the kink.

Predicted vs actual. Two views of how the probabilities behave out of sample: a calibration plot (bin clients by predicted risk, compare mean prediction to the observed default rate — on the diagonal = well-calibrated) and ROC curves. The tree's few leaves emit only a handful of distinct probabilities (a steppy curve); logit gives smooth, granular scores.

calib <- function(p,y,k=10){ b<-cut(p, unique(quantile(p,0:k/k)), include.lowest=TRUE); data.frame(pred=tapply(p,b,mean), obs=tapply(y,b,mean)) }
cl<-calib(p.lg,yte); ct<-calib(p.pru,yte)
par(mfrow=c(1,2), mar=c(4,4,3,1))
plot(cl$pred,cl$obs,type="b",col="#dd6b20",pch=19,xlim=c(0,.8),ylim=c(0,.8),xlab="predicted P(default)",ylab="observed default rate",main="Predicted vs actual (calibration)")
lines(ct$pred,ct$obs,type="b",col="#2b6cb0",pch=17); abline(0,1,lty=3); legend("topleft",c("logistic","tree"),col=c("#dd6b20","#2b6cb0"),pch=c(19,17),bty="n")
roc <- function(p,y){ o<-order(-p); list(fp=cumsum(1-y[o])/sum(1-y), tp=cumsum(y[o])/sum(y)) }
r1<-roc(p.lg,yte); r2<-roc(p.pru,yte)
plot(r1$fp,r1$tp,type="l",col="#dd6b20",lwd=2,xlab="false-positive rate",ylab="true-positive rate",main="ROC curves")
lines(r2$fp,r2$tp,col="#2b6cb0",lwd=2); abline(0,1,lty=3); legend("bottomright",c(sprintf("logistic (AUC %.3f)",auc(yte,p.lg)),sprintf("tree (AUC %.3f)",auc(yte,p.pru))),col=c("#dd6b20","#2b6cb0"),lty=1,lwd=2,bty="n")
par(mfrow=c(1,1))
cat("Econometrics and ML are complements: a scorecard wants the odds ratios, a challenger model wants the tree's flexibility.\n")

Econometrics and ML are complements: a scorecard wants the odds ratios, a challenger model wants the tree's flexibility.

6. Regression trees

With method="anova", rpart splits on variance reduction — a regression tree. The data: California housing (Pace & Barry, 1997, 1990 US Census) — one row per census block group (~600–3,000 residents; 20,640 total), target = the block group's median house value ($100,000s, capped at 5.0), from 8 features led by median income, plus house age, rooms/occupancy and latitude/longitude (location matters). Reported out of sample.

options(repr.plot.width=9, repr.plot.height=5.5)
h <- read.csv("cali_housing.csv"); set.seed(0); j <- sample(nrow(h), 0.7*nrow(h)); htr <- h[j,]; hte <- h[-j,]
rfit <- rpart(MedHouseVal ~ ., data=htr, method="anova", control=rpart.control(cp=0.005))
rpart.plot(rfit, type=2, box.palette="Blues", main="rpart regression tree — California housing")
cat(sprintf("regression tree OUT-OF-SAMPLE RMSE: %.4f (median house value, $100k)\n", sqrt(mean((predict(rfit,hte)-hte$MedHouseVal)^2))))

regression tree OUT-OF-SAMPLE RMSE: 0.7512 (median house value, $100k)

Proportions vs predictions — the regression tree

The classification side of this notebook already plots empirical proportions against fitted curves. The regression tree deserves the same treatment: bin the held-out block groups by predicted value and plot the mean actual value in each bin. On the diagonal the prediction is the outcome. A tree predicts a small set of constants, so the shape of this plot also shows how coarse that partition is.

options(repr.plot.width=7, repr.plot.height=5); par(mar=c(4,4,3,1))
prd <- predict(rfit, hte); act <- hte$MedHouseVal
q <- unique(quantile(prd, seq(0,1,length=11))); b <- cut(prd, q, include.lowest=TRUE)
pm <- as.numeric(tapply(prd,b,mean)); am <- as.numeric(tapply(act,b,mean))
plot(pm, am, type="b", pch=19, col="#2b6cb0", lwd=2, xlab="regression-tree predicted ($100k)",
     ylab="mean actual value", main=sprintf("Proportions vs predictions (RMSE %.3f)", sqrt(mean((prd-act)^2))))
abline(0,1,lty=2); legend("topleft", c("perfect","binned means"), lty=c(2,1), pch=c(NA,19),
                          col=c("black","#2b6cb0"), lwd=2, bty="n", cex=0.8)
cat(sprintf("distinct predicted values the tree can emit: %d\n", length(unique(prd))))
cat(sprintf("binned gaps (actual - predicted) run %+.3f to %+.3f; slope of actual on predicted %.3f\n",
            min(am-pm), max(am-pm), coef(lm(act~prd))[2]))
cat("\nA slope near 1 with small gaps means the tree is unbiased where it does predict. What the plot cannot hide is\n")
cat("the coarseness: every block group is assigned one of a handful of leaf constants, so the fitted values come in\n")
cat("steps rather than a continuum. That is the limitation the ensembles in the following notebooks remove by\n")
cat("averaging many such partitions together.\n")

distinct predicted values the tree can emit: 14

binned gaps (actual - predicted) run -0.025 to +0.080; slope of actual on predicted 1.006

A slope near 1 with small gaps means the tree is unbiased where it does predict. What the plot cannot hide is

the coarseness: every block group is assigned one of a handful of leaf constants, so the fitted values come in

steps rather than a continuum. That is the limitation the ensembles in the following notebooks remove by

averaging many such partitions together.

7. Summary

rpart reproduces the from-scratch CART on the real credit-default data (root split PAY_1), with cost-complexity pruning via cp and rpart's built-in cross-validation, and a variance-split regression tree on California housing.

What we conclude — out of sample. On the 30% held-out clients: (1) recent repayment behaviour is destiny — PAY_1 is the first split and the top-importance variable, and a client two-or-more months behind defaults several times as often as one who paid on time; (2) the pruned tree is far simpler and more accurate, reaching an out-of-sample AUC of ~0.74; (3) the traditional econometric model — logistic regression — is competitive (AUC ~0.72) and the most interpretable, with PAY_1's odds ratio (~2) stating that each extra month in arrears doubles the odds of default. Econometrics and ML are complements: the tree captures nonlinearity automatically, the logit hands you defensible odds ratios.

This is the package mirror of cart_python.ipynb and establishes the ML arc's hybrid format (from-scratch core + field-standard package). Next — Random Forests (ranger in R) — averages many de-correlated trees to defeat the single tree's instability.