Bayesian Nonparametric Prediction (R) — mgcv · kernlab

Gaussian processes and additive splines on the two datasets

The R cross-check for the Python BNP benchmark: kernlab::gausspr (Gaussian process regression and classification) and mgcv (Simon Wood's package — the reference implementation of penalized additive splines / GAMs, and the same package behind Bayesian Penalised Splines & Additive Models). Same Taiwan credit-default (classification) and California-housing (regression) data as the whole ML arc. The question is the Python notebook's: how do Bayesian nonparametric predictors fare, and what do they add? GPs are $O(n^3)$ so they use a subsample (as BART — Bayesian Additive Regression Trees did); GAMs use the full data. ROC-AUC is defined in the CART notebook.

options(repr.plot.width=13, repr.plot.height=4.6)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(kernlab); library(mgcv)})
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) }
rmse <- function(a,b) sqrt(mean((a-b)^2))
cat("kernlab", as.character(packageVersion("kernlab")), "| mgcv", as.character(packageVersion("mgcv")), "\n")

kernlab 0.9.33 | mgcv 1.9.4 

1. Gaussian process regression — California

kernlab::gausspr fits a GP with an RBF kernel, choosing the kernel width automatically. On a 1,200-point subsample it returns a posterior mean for every block group; we score out-of-sample RMSE and plot predicted vs actual with decile means. The GP is the Bayesian kernel machine — its posterior mean is the kernel ridge regression of Support Vector Machines & Kernel Methods, where the two are shown to agree to 5.5×10⁻¹⁴ using kernlab::gausspr and a hand-built kernel-ridge solve.

h<-read.csv("cali_housing.csv"); hf<-setdiff(names(h),"MedHouseVal")
set.seed(0); j<-sample(nrow(h),0.7*nrow(h)); htr<-h[j,]; hte<-h[-j,]
ctr<-colMeans(htr[,hf]); scl<-apply(htr[,hf],2,sd)
Xtr<-scale(htr[,hf],ctr,scl); Xte<-scale(hte[,hf],ctr,scl)
set.seed(1); s<-sample(nrow(Xtr),1200)
gpr<-gausspr(Xtr[s,], htr$MedHouseVal[s], kernel="rbfdot", kpar="automatic")
mu<-as.numeric(predict(gpr, Xte)); gp_rmse<-rmse(hte$MedHouseVal, mu)
cat(sprintf("GP regression (subsample 1200): test RMSE %.3f\n", gp_rmse))
options(repr.plot.width=6.4, repr.plot.height=5)
plot(mu, hte$MedHouseVal, pch=".", col="#a0aec0", xlab="GP predicted value ($100k)", ylab="actual value",
     main=sprintf("GP regression -- predicted vs actual (RMSE %.3f)", gp_rmse)); abline(0,1,lty=2)
qb<-quantile(mu,seq(0,1,length=11)); bb<-cut(mu,unique(qb),include.lowest=TRUE)
points(tapply(mu,bb,mean), tapply(hte$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")

Using automatic sigma estimation (sigest) for RBF or laplace kernel 

GP regression (subsample 1200): test RMSE 0.611

2. Gaussian process classification — credit default

gausspr with a factor response fits GP classification (a latent GP squashed through a link, Laplace-approximated). On a 1,200-client subsample we score out-of-sample AUC and the reliability curve. It should beat the linear logistic baseline (0.715) despite the tiny training sample, the kernel supplying nonlinearity.

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,]
cc<-colMeans(tr[,feat]); cs<-apply(tr[,feat],2,sd)
Ctr<-scale(tr[,feat],cc,cs); Cte<-scale(te[,feat],cc,cs)
set.seed(1); s2<-sample(nrow(Ctr),1200)
gpc<-gausspr(Ctr[s2,], factor(tr$default[s2]), kernel="rbfdot", kpar="automatic")
pg<-predict(gpc, Cte, type="probabilities")[,"1"]; gp_auc<-auc(te$default, pg)
cat(sprintf("GP classification (subsample 1200): test AUC %.3f\n", gp_auc))
options(repr.plot.width=13, repr.plot.height=4.6); par(mfrow=c(1,2), mar=c(4,4,3,1))
# ROC
th<-sort(unique(pg)); tpr<-sapply(th,function(t)mean(pg[te$default==1]>=t)); fpr<-sapply(th,function(t)mean(pg[te$default==0]>=t))
plot(fpr,tpr,type="l",col="#2b6cb0",lwd=2,xlab="false positive rate",ylab="true positive rate",main="Credit -- ROC"); abline(0,1,lty=3)
legend("bottomright", sprintf("GP classifier (AUC %.3f)",gp_auc), col="#2b6cb0", lwd=2, bty="n")
# reliability
q<-quantile(pg,seq(0,1,length=11)); b<-cut(pg,unique(q),include.lowest=TRUE)
plot(tapply(pg,b,mean), tapply(te$default,b,mean), type="b", pch=19, col="#c53030", lwd=2,
     xlab="predicted P(default) [decile bins]", ylab="observed default proportion", main="Credit -- reliability")
abline(0,1,lty=2); par(mfrow=c(1,1))
cat(sprintf("GP classification AUC %.3f beats the logistic baseline (0.715) from only 1200 of %d training rows.\n", gp_auc, nrow(tr)))

Using automatic sigma estimation (sigest) for RBF or laplace kernel 

GP classification (subsample 1200): test AUC 0.739

GP classification AUC 0.739 beats the logistic baseline (0.715) from only 1200 of 21000 training rows.

3. Additive splines — mgcv GAM

mgcv fits a GAM as a sum of penalized smooths $g(E[y])=\beta_0+\sum_j f_j(x_j)$, choosing the smoothness by (restricted) marginal likelihood — the penalty is a Bayesian random-walk prior on the spline coefficients, so mgcv reports genuine credible bands. We fit gam on the full California data (Gaussian) and a bam binomial GAM on credit, then read each feature's effect straight off its smooth, drawn over the central 99% of each covariate so that a handful of extreme block groups cannot compress the informative range. This is the reference version of the Python notebook's pyGAM.

# California: additive smooths on all 8 features
fH<-as.formula(paste("MedHouseVal ~", paste(sprintf("s(%s)",hf),collapse="+")))
gamH<-bam(fH, data=htr); gam_rmse<-rmse(hte$MedHouseVal, predict(gamH, hte))
# credit: smooths on continuous features, linear terms for discrete
cont<-c("LIMIT_BAL","AGE",paste0("BILL_AMT",1:6),paste0("PAY_AMT",1:6))
disc<-setdiff(feat,cont)
fC<-as.formula(paste("default ~", paste(c(sprintf("s(%s)",cont),disc),collapse="+")))
gamC<-bam(fC, data=tr, family=binomial); gam_auc<-auc(te$default, predict(gamC, te, type="response"))
cat(sprintf("mgcv GAM  -- California test RMSE %.3f | credit test AUC %.3f\n", gam_rmse, gam_auc))
# Drawn from predict(type="terms") rather than plot.gam. Two reasons: plot.gam's scale=-1 default puts
# every smooth on ONE shared y-axis, which here spans [-38, 43] because the AveRooms and AveBedrms bands
# explode out in the sparse tail -- on that axis MedInc occupies 4% of the height and looks flat. And its
# `shade` argument is not a formal of plot.gam in mgcv 1.9.4, so shade=TRUE is silently ignored.
smooth_band <- function(model, var, data, col, main, ylab, hi=0.99, n=200) {
  nd <- data[rep(1, n), , drop=FALSE]
  for (v in names(nd)) if (is.numeric(nd[[v]])) nd[[v]] <- median(data[[v]])
  nd[[var]] <- seq(min(data[[var]]), quantile(data[[var]], hi), length.out=n)
  pr <- predict(model, nd, type="terms", se.fit=TRUE); k <- paste0("s(", var, ")")
  fit <- pr$fit[,k]; se <- pr$se.fit[,k]; lo <- fit-1.96*se; up <- fit+1.96*se
  plot(nd[[var]], fit, type="n", ylim=range(lo,up), xlab=var, ylab=ylab, main=main)
  polygon(c(nd[[var]], rev(nd[[var]])), c(up, rev(lo)), col=paste0(col,"33"), border=NA)
  lines(nd[[var]], fit, col=col, lwd=2); abline(h=0, lty=3, col="grey50")
  rug(quantile(data[[var]], seq(0, hi, length.out=40)), col="grey60")
  invisible(data.frame(x=nd[[var]], fit=fit))
}
options(repr.plot.width=14, repr.plot.height=4.2); par(mfrow=c(1,3), mar=c(4,4,3,1))
s1<-smooth_band(gamH,"MedInc",   htr,"#2b6cb0","GAM smooth: MedInc (California)",  "partial effect on value")
s2<-smooth_band(gamH,"AveRooms", htr,"#2b6cb0","GAM smooth: AveRooms (California)","partial effect on value")
s3<-smooth_band(gamC,"LIMIT_BAL",tr, "#c53030","GAM smooth: credit limit (credit)","partial effect on log-odds")
par(mfrow=c(1,1))
rmin<-s2$x[which.min(s2$fit)]; occ_ratio<-max(htr$AveOccup)/quantile(htr$AveOccup,0.99)
cat(sprintf("Each smooth is a readable curve with a credible band. MedInc climbs steeply and flattens above ~8; AveRooms\n"))
cat(sprintf("is U-shaped with its minimum near %.1f rooms, so both the smallest and the largest units carry a premium over\n", rmin))
cat(sprintf("the middle once income is held fixed; and the credit-limit effect is near-linear (edf %.1f), higher limits\n", summary(gamC)$edf[which(cont=="LIMIT_BAL")]))
cat(sprintf("going with lower default odds. The x-axes stop at the 99th percentile deliberately: AveOccup's maximum is\n"))
cat(sprintf("%.0fx its p99, and letting a handful of such blocks set the range compresses every curve into a flat line.\n", occ_ratio))

mgcv GAM  -- California test RMSE 0.607 | credit test AUC 0.743

Each smooth is a readable curve with a credible band. MedInc climbs steeply and flattens above ~8; AveRooms

is U-shaped with its minimum near 5.0 rooms, so both the smallest and the largest units carry a premium over

the middle once income is held fixed; and the credit-limit effect is near-linear (edf 1.4), higher limits

going with lower default odds. The x-axes stop at the 99th percentile deliberately: AveOccup's maximum is

234x its p99, and letting a handful of such blocks set the range compresses every curve into a flat line.

4. Scoreboard and summary

The Bayesian nonparametric predictors, computed here, beside the frequentist field. The logistic and linear baselines are refitted in this notebook; the forest, XGBoost and BART figures are carried from the R tree notebooks, which share this notebook's set.seed(0) 70/30 split. They are deliberately not the Python notebook's numbers — that split is a different one, and mixing the two would make the bars incomparable.

# baselines computed HERE on this notebook's split; ensemble numbers carried from the R tree notebooks,
# which use the same set.seed(0) 70/30 split -- Python's figures are NOT interchangeable with these.
glm_auc<-auc(te$default, predict(glm(default~., tr, family=binomial), te, type="response"))
lin_rmse<-rmse(hte$MedHouseVal, predict(lm(MedHouseVal~., htr), hte))
credit<-c("logistic"=glm_auc,"random forest"=0.774,"XGBoost"=0.776,"BART (dbarts)"=0.783,
          "GP (kernlab)"=gp_auc,"GAM (mgcv)"=gam_auc)
cali  <-c("linear"=lin_rmse,"random forest"=0.491,"XGBoost"=0.484,"BART (dbarts)"=0.526,
          "GP (kernlab)"=gp_rmse,"GAM (mgcv)"=gam_rmse)
bnp<-c("GP (kernlab)","GAM (mgcv)","BART (dbarts)")
options(repr.plot.width=14, repr.plot.height=4.6); par(mfrow=c(1,2), mar=c(4,7,3,1))
barplot(rev(credit),horiz=TRUE,las=1,xlim=c(0.5,0.8),col=rev(ifelse(names(credit)%in%bnp,"#6b46c1","#a0aec0")),main="Credit -- test AUC (higher better)")
barplot(rev(cali),horiz=TRUE,las=1,col=rev(ifelse(names(cali)%in%bnp,"#6b46c1","#a0aec0")),main="California -- test RMSE (lower better)")
par(mfrow=c(1,1))
cat(sprintf("On California all three Bayesian nonparametric methods (purple) clear the linear model decisively: GAM %.3f and\n", gam_rmse))
cat(sprintf("GP %.3f against %.3f, with dbarts BART at 0.526 -- close enough to the forest's 0.491 to be a fair fight. Credit is\n", gp_rmse, lin_rmse))
cat(sprintf("the harder verdict: the mgcv GAM reaches %.3f and the GP %.3f against ensembles at 0.774-0.776, so they trail by\n", gam_auc, gp_auc))
cat("0.03 rather than matching. Part of that is a modelling choice -- the smooths here are fitted only to the continuous\n")
cat("features, with the discrete PAY_* and demographic terms entered linearly, where pyGAM splines every column and gets\n")
cat("0.765. What all three carry and none of the frequentist winners do is calibrated uncertainty: credible bands on every\n")
cat("smooth, a posterior variance at every input. Boosting still leads raw accuracy at full scale.\n")

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

On California all three Bayesian nonparametric methods (purple) clear the linear model decisively: GAM 0.607 and

GP 0.611 against 0.734, with dbarts BART at 0.526 -- close enough to the forest's 0.491 to be a fair fight. Credit is

the harder verdict: the mgcv GAM reaches 0.743 and the GP 0.739 against ensembles at 0.774-0.776, so they trail by

0.03 rather than matching. Part of that is a modelling choice -- the smooths here are fitted only to the continuous

features, with the discrete PAY_* and demographic terms entered linearly, where pyGAM splines every column and gets

0.765. What all three carry and none of the frequentist winners do is calibrated uncertainty: credible bands on every

smooth, a posterior variance at every input. Boosting still leads raw accuracy at full scale.

Summary

kernlab and mgcv reproduce the Python benchmark: a Gaussian process (the Bayesian kernel machine, posterior mean = kernel ridge) and an additive-spline GAM (penalized smooths = a Bayesian random-walk prior, hence credible bands) are competitive with the tree ensembles — clearing the linear model decisively on California, and trailing the ensembles by about 0.03 AUC on credit — while adding the calibrated uncertainty the frequentist models lack. Their cost is the GP's $O(n^3)$ scaling (hence subsampling), the reason boosting remains the default for large tabular data. This is the package mirror of bnp_benchmark_python.ipynb, cross-linking the ML arc to the Bayesian Nonparametrics section — Gaussian-Process Regression, GP Classification & Log-Gaussian Cox Processes, Bayesian Penalised Splines & Additive Models — and to BART — Bayesian Additive Regression Trees.