Gradient Boosting (R) — gbm¶

Friedman's gradient boosting machine, the reference implementation¶

gbm (Ridgeway) is the canonical R implementation of Friedman's gradient boosting — the package cross-check for the from-scratch gboost.py. It fits the same additive, gradient-driven sequence of shallow trees, chooses the number of trees by built-in cross-validation (gbm.perf), and reports relative influence (its variable importance). Same Taiwan credit-default and California-housing data; out of sample on a 30% held-out test set. ROC-AUC = P(model ranks a random defaulter above a random non-defaulter); 0.5 = chance, 1 = perfect (defined in the CART notebook).

In [1]:
options(repr.plot.width=8, repr.plot.height=5)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
for(p in c("gbm")) if(!requireNamespace(p,quietly=TRUE)) install.packages(p, repos="https://cloud.r-project.org", quiet=TRUE)
suppressMessages(library(gbm))
d <- read.csv("credit_default.csv")                                   # default stays numeric 0/1 for bernoulli
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) }
cat(sprintf("credit default: %d train / %d test; default rate %.1f%%\n", nrow(tr), nrow(te), 100*mean(d$default)))

Warning message:
"package 'gbm' was built under R version 4.6.1"

credit default: 21000 train / 9000 test; default rate 22.1%

1. Fit gbm and choose the number of trees¶

gbm(default ~ ., distribution="bernoulli", shrinkage=0.05, interaction.depth=3, cv.folds=3) grows the boosted sequence; gbm.perf plots training vs cross-validated error against iteration and returns the CV-optimal tree count — boosting's early-stopping point, exactly the peak the Python notebook showed.

In [2]:
set.seed(1)
invisible(capture.output(m <- gbm(default ~ ., data=tr, distribution="bernoulli", n.trees=600, interaction.depth=3,
         shrinkage=0.05, bag.fraction=0.8, cv.folds=3, n.cores=1, verbose=FALSE)))
best <- gbm.perf(m, method="cv", plot.it=TRUE)
p <- predict(m, te, n.trees=best, type="response")
cat(sprintf("CV-optimal trees: %d;  out-of-sample AUC %.4f\n", best, auc(te$default, p)))
cat("Training error keeps falling (black); the CV error (green) bottoms out and turns up -- the point past which more trees overfit.\n")

CV-optimal trees: 356;  out-of-sample AUC 0.7780

Training error keeps falling (black); the CV error (green) bottoms out and turns up -- the point past which more trees overfit.
No description has been provided for this image

2. Relative influence¶

gbm measures each variable's contribution to loss reduction across all trees ("relative influence"). As in every tree model here, recent repayment status (PAY_*) dominates.

In [3]:
ri <- summary(m, n.trees=best, plotit=FALSE)
par(mar=c(4,7,3,1)); barplot(rev(head(ri$rel.inf,8)), names.arg=rev(head(ri$var,8)), horiz=TRUE, las=1,
        col="#2b6cb0", main="gbm relative influence", xlab="relative influence (%)")
No description has been provided for this image

3. How boosting fits — graphically¶

One-feature views (as in the forest/BART notebooks): $P(\text{default})$ vs credit limit and California value vs median income, comparing the linear/logistic baseline, a single rpart tree, and the gbm booster. Boosting traces a smooth, bent curve where the parametric fit is straight and the single tree steps.

In [4]:
suppressMessages(library(rpart)); par(mfrow=c(1,2), mar=c(4,4,3,1))
# classification: P(default) vs credit limit
lim <- tr$LIMIT_BAL/1000; qb <- unique(quantile(lim,seq(0,1,length=16))); b <- findInterval(lim,qb,all.inside=TRUE)
ctr <- tapply(lim,b,mean); emp <- tapply(tr$default,b,mean); g <- seq(min(lim),quantile(lim,0.99),length=300)
lo <- predict(glm(default~LIMIT_BAL, tr, family=binomial), data.frame(LIMIT_BAL=g*1000), type="response")
t1 <- predict(rpart(factor(default)~LIMIT_BAL, tr, control=rpart.control(maxdepth=4)), data.frame(LIMIT_BAL=g*1000), type="prob")[,2]
set.seed(1); mg <- gbm(default~LIMIT_BAL, data=tr, distribution="bernoulli", n.trees=400, interaction.depth=3, shrinkage=0.05, verbose=FALSE)
gbx <- predict(mg, data.frame(LIMIT_BAL=g*1000), n.trees=400, type="response")
plot(ctr,emp,pch=19,xlab="credit limit (NT$ thousands)",ylab="P(default)",main="Classification: P(default) vs credit limit",ylim=c(0,max(emp)*1.1))
lines(g,lo,col="#dd6b20",lwd=2); lines(g,t1,col="#a0aec0",lwd=2); lines(g,gbx,col="#2b6cb0",lwd=2)
legend("topright",c("empirical","logistic","single tree","gradient boosting"),col=c("black","#dd6b20","#a0aec0","#2b6cb0"),pch=c(19,NA,NA,NA),lty=c(NA,1,1,1),lwd=2,bty="n",cex=0.8)
# regression: California value vs income
h <- read.csv("cali_housing.csv"); set.seed(0); jj<-sample(nrow(h),0.7*nrow(h)); htr<-h[jj,]
xi <- htr$MedInc; gh <- seq(min(xi),quantile(xi,0.99),length=300)
ols <- predict(lm(MedHouseVal~MedInc, htr), data.frame(MedInc=gh))
trg <- predict(rpart(MedHouseVal~MedInc, htr, control=rpart.control(maxdepth=4)), data.frame(MedInc=gh))
set.seed(1); mgr <- gbm(MedHouseVal~MedInc, data=htr, distribution="gaussian", n.trees=400, interaction.depth=3, shrinkage=0.05, verbose=FALSE)
grx <- predict(mgr, data.frame(MedInc=gh), n.trees=400)
plot(xi,htr$MedHouseVal,pch=".",col="grey",xlab="median income",ylab="median house value ($100k)",main="Regression: value vs income",xlim=c(min(xi),quantile(xi,0.99)))
lines(gh,ols,col="#dd6b20",lwd=2); lines(gh,trg,col="#a0aec0",lwd=2); lines(gh,grx,col="#2b6cb0",lwd=2)
legend("topleft",c("linear","single tree","gradient boosting"),col=c("#dd6b20","#a0aec0","#2b6cb0"),lty=1,lwd=2,bty="n",cex=0.8)
par(mfrow=c(1,1))
No description has been provided for this image

Are the boosted probabilities honest?¶

gbm with type="response" returns a probability, and AUC never checks whether it is one. Binning the test set by predicted probability and plotting the observed default rate answers that directly — the same diagnostic the Python notebook applies, where the ordering turned out to contradict the usual expectation about boosting.

In [5]:
rel <- function(p, y, lab, col) {
  q <- unique(quantile(p, seq(0,1,length=11))); b <- cut(p, q, include.lowest=TRUE)
  pp <- tapply(p,b,mean); oo <- tapply(y,b,mean); w <- as.numeric(table(b))/length(y)
  ece <- sum(w*abs(oo-pp), na.rm=TRUE)
  list(pp=as.numeric(pp), oo=as.numeric(oo), ece=ece, lab=lab, col=col)
}
draw <- function(L, main) {
  xm <- max(sapply(L, function(z) max(z$pp))); ym <- max(sapply(L, function(z) max(z$oo)))
  plot(NA, xlim=c(0,xm), ylim=c(0,max(xm,ym)), xlab="predicted P(default)",
       ylab="observed default rate", main=main); abline(0,1,lty=2)
  for (z in L) lines(z$pp, z$oo, type="b", pch=19, col=z$col, lwd=2)
  legend("topleft", sapply(L, function(z) sprintf("%s (ECE %.3f)", z$lab, z$ece)),
         col=sapply(L, function(z) z$col), lwd=2, pch=19, bty="n", cex=0.8)
}
options(repr.plot.width=7, repr.plot.height=5); par(mar=c(4,4,3,1))
pg  <- predict(m, te, n.trees=best, type="response")
pgl <- predict(glm(default~., tr, family=binomial), te, type="response")
L <- list(rel(pg,  te$default, "gbm",      "#2b6cb0"),
          rel(pgl, te$default, "logistic", "#a0aec0"))
draw(L, "Reliability: do the probabilities mean what they say?")
cat(sprintf("gbm      AUC %.3f  ECE %.3f  mean predicted %.3f\n", auc(te$default,pg),  L[[1]]$ece, mean(pg)))
cat(sprintf("logistic AUC %.3f  ECE %.3f  mean predicted %.3f\n", auc(te$default,pgl), L[[2]]$ece, mean(pgl)))
cat(sprintf("base default rate in the test set: %.3f\n", mean(te$default)))
cat("\nThe folklore says boosting is the badly calibrated one, because stage-wise fitting of the log-odds pushes\n")
cat("confident cases outward. Read the two ECE figures before accepting that: the Python notebook finds the\n")
cat("opposite on this data, with a regularised boosting fit better calibrated than the logistic, which is underfit\n")
cat("against the nonlinearity in PAY_1. Miscalibration is a property of a fit, not of a model class.\n")

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

gbm      AUC 0.778  ECE 0.008  mean predicted 0.222

logistic AUC 0.712  ECE 0.052  mean predicted 0.221

base default rate in the test set: 0.221

The folklore says boosting is the badly calibrated one, because stage-wise fitting of the log-odds pushes

confident cases outward. Read the two ECE figures before accepting that: the Python notebook finds the

opposite on this data, with a regularised boosting fit better calibrated than the logistic, which is underfit

against the nonlinearity in PAY_1. Miscalibration is a property of a fit, not of a model class.
No description has been provided for this image

4. Summary¶

gbm reproduces the from-scratch gradient booster: an additive sequence of shallow trees whose CV-optimal count is chosen by gbm.perf (the early-stopping peak the Python notebook drew), relative influence dominated by recent repayment status, and a smooth, bent fit that beats the single tree and the linear/logistic baseline out of sample. It is the package mirror of gboost_python.ipynb.

Boosting cuts bias by adding dependent shallow trees, the counterpart to the random forest's variance-cutting independent deep trees. The next notebook takes this to its industrial form — XGBoost / LightGBM / CatBoost — with second-order boosting, histogram splitting and native categorical handling. The Bayesian cousin is BART.