GP Classification & Log-Gaussian Cox Processes (R)¶

mgcv — smooth logistic and Poisson surfaces, with glm as the parametric baseline¶

The R counterpart to gplink_python.ipynb. Where the Python notebook builds the latent-GP classifier and Cox process from scratch (Laplace approximation) and cross-checks in PyMC, this notebook uses the standard smoothing package mgcv: a smooth binomial GAM is the frequentist cousin of the GP classifier, and a smooth Poisson GAM on binned event counts is the cousin of the log-Gaussian Cox process. Plain glm gives the parametric baseline. Data: a non-monotone binary curve, the Pima diabetes surface, and the coal-mining disasters.

In [1]:
options(repr.plot.width=12, repr.plot.height=4.4)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages(library(mgcv))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; PURP<-"#6b46c1"; GREY<-"#718096"
cat("R engine: mgcv smooth GAMs -- binomial (GP classifier cousin), Poisson (Cox-process cousin)\n")

R engine: mgcv smooth GAMs -- binomial (GP classifier cousin), Poisson (Cox-process cousin)

1. A non-monotone probability — smooth GAM vs logistic¶

A binary outcome with $p(x)=\sigma(1.5\sin(1.2x))$: the probability rises and falls, which glm logistic (linear log-odds) cannot bend to but a smooth mgcv binomial GAM can — the frequentist analogue of the GP classifier.

In [2]:
set.seed(1); x<-sort(runif(220,-3.2,3.2)); ptrue<-function(t) plogis(1.5*sin(1.2*t)); y<-as.numeric(runif(220)<ptrue(x))
xg<-seq(-3.2,3.2,length=200)
g<-gam(y~s(x), family=binomial); pg<-predict(g,data.frame(x=xg),type="response",se.fit=TRUE)
glm1<-glm(y~x, family=binomial); pl<-predict(glm1,data.frame(x=xg),type="response")
eta<-predict(g,data.frame(x=xg),se.fit=TRUE); lo<-plogis(eta$fit-1.96*eta$se.fit); hi<-plogis(eta$fit+1.96*eta$se.fit)
cat(sprintf("GAM RMSE(p) %.3f  vs logistic RMSE %.3f\n", sqrt(mean((pg$fit-ptrue(xg))^2)), sqrt(mean((pl-ptrue(xg))^2))))
par(mar=c(4,4,3,1)); plot(xg, ptrue(xg), type="l", lty=2, lwd=1.6, ylim=c(0,1), xlab="x", ylab="P(y=1)", main="Smooth binomial GAM recovers the wave; logistic cannot")
polygon(c(xg,rev(xg)), c(lo,rev(hi)), col=adjustcolor(BLUE,.18), border=NA)
lines(xg, pg$fit, col=BLUE, lwd=2.4); lines(xg, pl, col=RED, lwd=1.8, lty=3); points(x, y*1.03-.015, pch=19, col=adjustcolor(GREY,.4), cex=.5)
legend("topright", c("true p(x)","mgcv binomial GAM","glm logistic"), col=c("black",BLUE,RED), lwd=2, lty=c(2,1,3), bty="n", cex=.8)
cat("The smooth GAM tracks the rise-and-fall (like the GP classifier in Python); logistic is forced to one monotone S.\n")

GAM RMSE(p) 0.084  vs logistic RMSE 0.205

The smooth GAM tracks the rise-and-fall (like the GP classifier in Python); logistic is forced to one monotone S.
No description has been provided for this image

2. The Pima data — a 2-D probability surface¶

Diabetes against glucose and BMI. mgcv's tensor smooth s(glu, bmi) draws a smooth 2-D probability surface; glm gives the linear-boundary baseline. On these features the classes are near-linearly separable, so accuracies nearly match — the value is the smooth surface and its uncertainty.

In [3]:
d<-read.csv("pima.csv")
gm<-gam(diabetes~s(glu,bmi), data=d, family=binomial); glm2<-glm(diabetes~glu+bmi, data=d, family=binomial)
ari_acc<-function(p,y) mean((p>.5)==y)
cat(sprintf("accuracy  mgcv-GAM %.3f   glm logistic %.3f\n", ari_acc(fitted(gm),d$diabetes), ari_acc(fitted(glm2),d$diabetes)))
g1<-seq(min(d$glu),max(d$glu),length=60); g2<-seq(min(d$bmi),max(d$bmi),length=60); gr<-expand.grid(glu=g1,bmi=g2)
P<-matrix(predict(gm,gr,type="response"),60,60)
par(mar=c(4,4,3,1))
image(g1,g2,P,col=hcl.colors(24,"Blue-Red 3"),xlab="glucose",ylab="BMI",main="P(diabetes) surface (mgcv GAM) + 0.5 boundary")
contour(g1,g2,P,levels=0.5,add=TRUE,lwd=2)
points(d$glu[d$diabetes==1],d$bmi[d$diabetes==1],pch=19,col=adjustcolor("darkred",.5),cex=.5)
points(d$glu[d$diabetes==0],d$bmi[d$diabetes==0],pch=19,col=adjustcolor("navy",.4),cex=.5)
cat("The GAM and logistic accuracies nearly match -- Pima is close to linearly separable in these two features; the\n")
cat("nonparametric payoff is the smooth surface and the uncertainty away from the data, as the Python GP shows.\n")

accuracy  mgcv-GAM 0.773   glm logistic 0.773

The GAM and logistic accuracies nearly match -- Pima is close to linearly separable in these two features; the

nonparametric payoff is the smooth surface and the uncertainty away from the data, as the Python GP shows.
No description has been provided for this image

3. Coal-mining disasters as a Cox process¶

Bin the 191 events and fit a smooth Poisson GAM to the counts — mgcv's version of the log-Gaussian Cox process: it recovers the intensity and the famous post-1890 fall with a confidence band. (The same series is modelled as a change-point process — abrupt piecewise-constant jumps rather than a smooth rate — in Reversible-Jump MCMC — Inference over Model Dimension; both find the fall near 1890.)

In [4]:
ev<-read.csv("coal.csv",header=FALSE)[,1]
edges<-seq(1851,1963,length=57); w<-diff(edges)[1]; cx<-head(edges,-1)+w/2; cy<-as.numeric(table(cut(ev,edges)))
gp<-gam(cy~s(cx,k=20), family=poisson, offset=rep(log(w),length(cy)))
pr<-predict(gp,se.fit=TRUE)
# The exposure offset is supplied through gam()'s `offset` ARGUMENT, which predict.gam() does
# not carry into its returned linear predictor. exp(pr$fit) is therefore ALREADY the per-year
# rate -- dividing by w again would halve it. (Check: fitted(gp) = exp(pr$fit)*w, the counts.)
lam<-exp(pr$fit); lo<-exp(pr$fit-1.96*pr$se.fit); hi<-exp(pr$fit+1.96*pr$se.fit)
par(mar=c(4,4,3,1))
plot(cx, cy/w, type="h", col=GREY, lwd=3, xlab="year", ylab="disasters per year",
     ylim=c(0, max(c(cy/w, hi))),      # include the bars, not just the band, or the tall ones clip
     main="Log-Gaussian Cox process (mgcv Poisson smooth): coal-mining intensity")
polygon(c(cx,rev(cx)), c(lo,rev(hi)), col=adjustcolor(BLUE,.2), border=NA); lines(cx, lam, col=BLUE, lwd=2.4)
legend("topright", c("binned rate","Poisson-smooth intensity","95% band"), col=c(GREY,BLUE,adjustcolor(BLUE,.4)), lwd=c(3,2,6), bty="n")
cat(sprintf("intensity ~%.1f/yr around 1860 falling to ~%.2f/yr by 1950 -- the sharp drop near 1890 (mine-safety laws),\n",
    lam[which.min(abs(cx-1860))], lam[which.min(abs(cx-1950))]))
cat(sprintf("For reference the raw event rate is %.2f/yr over 1851-1890 and %.2f/yr over 1900-1962, so the\n",
            sum(ev<1890)/(1890-1851), sum(ev>=1900)/(1962-1900)))
cat("smooth is on the right scale. This matches the from-scratch LGCP in the Python notebook, which is\n")
cat("the point of running two engines: the exposure offset is the easiest thing to get wrong in a Cox\n")
cat("process, and a factor-of-two error survives every internal check while looking perfectly plausible.\n")
cat("The smooth Poisson GAM is the Cox process's package cousin.\n")

# The factor-of-two, demonstrated rather than asserted: predict.gam() does not carry the
# offset ARGUMENT into its linear predictor, but fitted() does apply it.
i60 <- which.min(abs(cx-1860))
cat(sprintf("\nOffset check at %.0f -- exp(eta) = %.3f per year;  fitted() = %.3f events per bin;\n",
            cx[i60], lam[i60], fitted(gp)[i60]))
cat(sprintf("ratio = %.3f, which is exactly the bin width w = %.0f years. So exp(eta) is ALREADY a\n",
            fitted(gp)[i60]/lam[i60], w))
cat("per-year rate and dividing it by w again halves it. That is the factor-of-two error, and the\n")
cat(sprintf("raw rate of %.2f/yr over 1851-1890 confirms which of the two figures was the right one.\n",
            sum(ev<1890)/(1890-1851)))

intensity ~3.1/yr around 1860 falling to ~0.71/yr by 1950 -- the sharp drop near 1890 (mine-safety laws),

For reference the raw event rate is 3.15/yr over 1851-1890 and 0.90/yr over 1900-1962, so the

smooth is on the right scale. This matches the from-scratch LGCP in the Python notebook, which is

the point of running two engines: the exposure offset is the easiest thing to get wrong in a Cox

process, and a factor-of-two error survives every internal check while looking perfectly plausible.

The smooth Poisson GAM is the Cox process's package cousin.

Offset check at 1860 -- exp(eta) = 3.055 per year;  fitted() = 6.109 events per bin;

ratio = 2.000, which is exactly the bin width w = 2 years. So exp(eta) is ALREADY a

per-year rate and dividing it by w again halves it. That is the factor-of-two error, and the

raw rate of 3.15/yr over 1851-1890 confirms which of the two figures was the right one.
No description has been provided for this image

4. Summary¶

mgcv's smooth GAMs are the standard-package view of the linked GP: a smooth binomial GAM recovered a non-monotone probability that glm logistic cannot bend to, drew the Pima diabetes surface, and a smooth Poisson GAM recovered the coal-mining intensity and its post-1890 fall — the frequentist cousins of GP classification and the log-Gaussian Cox process. They cross-check the from-scratch Laplace-approximation GP and PyMC gp.Latent in gplink_python.ipynb.

GP classification is logistic regression with a nonparametric predictor; the log-Gaussian Cox process is Poisson regression with a GP log-rate — and mgcv fits both by penalised likelihood. The coal series is also modelled as a change-point process in the variable-selection arc. Next the arc turns to nonparametric survival and hazards.