XGBoost (R) — xgboost

Production gradient boosting, the R side

xgboost is the R interface to the same second-order (Newton) gradient-boosting engine as the Python notebook — the package cross-check for XGBoost/LightGBM/CatBoost (LightGBM and CatBoost also ship R packages; XGBoost is the most widely used, so we feature it here). We fit it on the Taiwan credit-default data with a DMatrix, use a watch-list for early stopping, read off gain importance, and compare its fitted shape to a single tree and the logistic baseline. Out of sample on a 30% test set. ROC-AUC = P(rank a random defaulter above a random non-defaulter); 0.5 = chance, 1 = perfect (defined in the CART notebook).

options(repr.plot.width=8, repr.plot.height=5)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages(library(xgboost))
d <- read.csv("credit_default.csv"); feat <- setdiff(names(d),"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) }
dtr <- xgb.DMatrix(as.matrix(tr[,feat]), label=tr$default)
dte <- xgb.DMatrix(as.matrix(te[,feat]), label=te$default)
cat(sprintf("xgboost %s;  credit default: %d train / %d test\n", as.character(packageVersion("xgboost")), nrow(tr), nrow(te)))

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

xgboost 3.2.1.1;  credit default: 21000 train / 9000 test

1. Fit with early stopping

xgb.cv runs cross-validation while boosting and, with early_stopping_rounds, records the train-vs-CV-test AUC each round so we can read off the tree count where the CV score peaks — the production default that removes manual tree-count tuning. We then refit at that optimal number of rounds and score out of sample.

params <- list(objective="binary:logistic", eval_metric="auc", max_depth=4, eta=0.05, nthread=2)
set.seed(1)
cv <- xgb.cv(params, dtr, nrounds=600, nfold=4, early_stopping_rounds=30, verbose=0)
el <- cv$evaluation_log; best <- which.max(el$test_auc_mean)
plot(el$iter, el$train_auc_mean, type="l", col="#2b6cb0", lwd=2, ylim=range(c(el$train_auc_mean, el$test_auc_mean)),
     xlab="boosting round", ylab="AUC (4-fold CV)", main="XGBoost: early stopping at the CV-AUC peak")
lines(el$iter, el$test_auc_mean, col="#c53030", lwd=2); abline(v=best, col="#2f855a", lty=2)
legend("bottomright", c("train AUC","CV-test AUC",sprintf("optimal @ %d rounds",best)),
       col=c("#2b6cb0","#c53030","#2f855a"), lty=c(1,1,2), lwd=2, bty="n")
m <- xgb.train(params, dtr, nrounds=best, verbose=0)
cat(sprintf("CV-optimal rounds %d;  out-of-sample AUC %.4f\n", best, auc(te$default, predict(m, dte))))

CV-optimal rounds 97;  out-of-sample AUC 0.7756

2. Gain importance

xgb.importance ranks features by gain — the total loss reduction their splits deliver. As in every tree model in this subsection, recent repayment status (PAY_*) dominates.

imp <- xgb.importance(model=m)
par(mar=c(4,7,3,1)); xgb.plot.importance(imp[1:8,], main="XGBoost gain importance")

3. Calibration — proportions vs predictions

AUC only measures ranking; it says nothing about whether a predicted 0.30 really defaults 30% of the time. Left — credit: sort the test clients into ten equal-count bins by predicted probability and plot the observed default proportion against the mean predicted probability; points on the 45° line are perfectly calibrated (the Brier score — mean squared probability error — summarises it). Right — California: the regression analogue, predicted vs actual value with decile-binned means. This is the "proportions vs predictions" view of both data sets, matching the Python notebook.

par(mfrow=c(1,2), mar=c(4,4,3,1))
# credit reliability: xgboost vs logistic
pxg <- predict(m, dte)
plg <- predict(glm(default~., data=tr, family=binomial), te, type="response")
brier <- function(y,p) mean((p-y)^2)
rel <- function(y,p){ q<-quantile(p,seq(0,1,length=11)); b<-cut(p,unique(q),include.lowest=TRUE)
                      data.frame(pp=tapply(p,b,mean), ot=tapply(y,b,mean)) }
rx<-rel(te$default,pxg); rl<-rel(te$default,plg); mx0<-max(pxg,plg)
plot(c(0,mx0),c(0,mx0),type="l",lty=2,xlab="predicted P(default)  [decile bins]",ylab="observed default proportion",
     main="Credit default -- reliability curve")
points(rx$pp,rx$ot,type="b",pch=19,col="#2b6cb0",lwd=2); points(rl$pp,rl$ot,type="b",pch=19,col="#dd6b20",lwd=2)
legend("topleft",c(sprintf("xgboost (Brier %.3f)",brier(te$default,pxg)),sprintf("logistic (Brier %.3f)",brier(te$default,plg))),
       col=c("#2b6cb0","#dd6b20"),pch=19,lwd=2,bty="n",cex=0.85)
# california predicted vs actual
hh<-read.csv("cali_housing.csv"); set.seed(0); j<-sample(nrow(hh),0.7*nrow(hh)); htr2<-hh[j,]; hte2<-hh[-j,]
dhtr<-xgb.DMatrix(as.matrix(htr2[,setdiff(names(hh),"MedHouseVal")]),label=htr2$MedHouseVal)
dhte<-xgb.DMatrix(as.matrix(hte2[,setdiff(names(hh),"MedHouseVal")]))
mh<-xgb.train(list(objective="reg:squarederror",max_depth=4,eta=0.05,nthread=2),dhtr,nrounds=300,verbose=0)
phr<-predict(mh,dhte); rmse<-sqrt(mean((phr-hte2$MedHouseVal)^2))
plot(phr,hte2$MedHouseVal,pch=".",col="#a0aec0",xlab="predicted value ($100k)",ylab="actual value ($100k)",
     main=sprintf("California -- predicted vs actual (RMSE %.3f)",rmse))
abline(0,1,lty=2)
qb<-quantile(phr,seq(0,1,length=11)); bb<-cut(phr,unique(qb),include.lowest=TRUE)
points(tapply(phr,bb,mean),tapply(hte2$MedHouseVal,bb,mean),type="b",pch=19,col="#c53030",lwd=2)
legend("topleft",c("perfect","decile means"),col=c("black","#c53030"),lty=c(2,1),pch=c(NA,19),lwd=2,bty="n",cex=0.85)
par(mfrow=c(1,1))
cat("Credit: xgboost hugs the 45-degree line (low Brier); the logistic is slightly less calibrated in the risky deciles.\n")
cat("California: decile means track the diagonal, flattening at the top from the data set's $500k price CAP.\n")

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

Credit: xgboost hugs the 45-degree line (low Brier); the logistic is slightly less calibrated in the risky deciles.

California: decile means track the diagonal, flattening at the top from the data set's $500k price CAP.

4. Fitted shape — graphically

One-feature views (as throughout the subsection): $P(\text{default})$ vs credit limit and California value vs median income, comparing the logistic/linear baseline, a single rpart tree, and xgboost. Boosting bends to the nonlinear trend where the parametric fit is straight and the single tree steps.

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); mgc <- xgb.train(list(objective="binary:logistic",max_depth=4,eta=0.05,min_child_weight=50,nthread=2), xgb.DMatrix(matrix(tr$LIMIT_BAL,ncol=1),label=tr$default), nrounds=200, verbose=0)
gbx <- predict(mgc, xgb.DMatrix(matrix(g*1000,ncol=1)))
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","xgboost"),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 <- xgb.train(list(objective="reg:squarederror",max_depth=4,eta=0.05,min_child_weight=20,nthread=2), xgb.DMatrix(matrix(xi,ncol=1),label=htr$MedHouseVal), nrounds=300, verbose=0)
grx <- predict(mgr, xgb.DMatrix(matrix(gh,ncol=1)))
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","xgboost"),col=c("#dd6b20","#a0aec0","#2b6cb0"),lty=1,lwd=2,bty="n",cex=0.8)
par(mfrow=c(1,1))

5. Summary

xgboost in R reproduces the production booster: second-order gradient boosting fit through a DMatrix, early stopping on a validation watch-list (no manual tree count), gain importance dominated by recent repayment status, and the familiar smooth-and-bent fit that beats the single tree and the linear/logistic baseline. It is the package mirror of xgb_python.ipynb (LightGBM and CatBoost also have R packages; XGBoost is the workhorse).

This closes the Tree Ensembles subsection across both languages: CART → Random Forests → BART → Gradient Boosting → XGBoost/LightGBM/CatBoost. On clean tabular data their accuracy converges; the real choice is speed, tooling, categoricals, and — with BART — calibrated uncertainty.