¿Por qué este GAM logístico se ajusta tan mal?

Question

Estoy intentando crear un modelo de regresión logística con mgcv::gam con lo que creo que es un límite de decisión simple, pero el modelo que construyo tiene un rendimiento muy pobre. Un modelo de regresión local construido usando locfit::locfit en los mismos datos encuentra el límite muy fácilmente. Quiero añadir regresores paramétricos adicionales a mi modelo real, así que No quiero cambiar a una regresión puramente local.

Quiero entender por qué el GAM tiene problemas para ajustar los datos, y si hay formas de especificar los suavizados que puede funcionar mejor.

He aquí un ejemplo simplificado y reproducible:

La verdad sobre el terreno es 1 = el punto se encuentra dentro del círculo unitario, 0 si está fuera

por ejemplo, z = 1 si sqrt(x^2 + y^2) <= 1, 0 en caso contrario

Los datos observados son ruidosos, con falsos positivos y falsos negativos

Construye una regresión logística para predecir si un punto está dentro del círculo o no, basándose en las coordenadas cartesianas del punto coordenadas cartesianas.

La regresión local puede encontrar bien el límite (el contorno del 50% de probabilidad está muy cerca del círculo unitario), pero un GAM logístico sobreestima sistemáticamente sobreestima el tamaño del círculo para la misma banda de probabilidad.

library(ggplot2)
library(locfit)
library(mgcv)
library(plotrix)

set.seed(0)
radius <- 1 # actual boundary
n <- 10000 # data points
jit <- 0.5 # noise factor

# Simulate random data, add polar coordinates
df <- data.frame(x=runif(n,-3,3), y=runif(n,-3,3))
df$r <- with(df, sqrt(x^2+y^2))
df$theta <- with(df, atan(y/x))

# Noisy indicator for inside the boundary
df$inside <- with(df, ifelse(r < radius + runif(nrow(df),-jit,jit),1,0))

# Plot data, shows ragged edge
(ggplot(df, aes(x=x, y=y, color=inside)) + geom_point() + coord_fixed() + xlim(-4,4) + ylim(-4,4))

### Model boundary condition using x,y coordinates

### local regression finds the boundary pretty accurately
m.locfit <- locfit(inside ~ lp(x,y, nn=0.3), data=df, family="binomial")
plot(m.locfit, asp=1, xlim=c(-2,-2,2,2))
draw.circle(0,0,1, border="red")

### But GAM fits very poorly, also tried with fx=TRUE but didn't help
m.gam <- gam(inside ~ s(x,y), data=df, family=binomial)
plot(m.gam, trans=plogis, se=FALSE, rug=FALSE)
draw.circle(0,0,1, border="red")

### gam.check doesn't indicate a problem with the model itself
gam.check(m.gam)

Method: UBRE   Optimizer: outer newton
full convergence after 8 iterations.
Gradient range [5.41668e-10,5.41668e-10]
(score -0.815746 & scale 1).
Hessian positive definite, eigenvalue range [0.0002169789,0.0002169789].

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

           k'    edf k-index p-value
s(x,y) 29.000 13.795   0.973    0.08

#### Try using polar coordinates

### Again, locfit works well
m.locfit2 <- locfit(inside ~ lp(r, nn=0.3), data=df, family="binomial")
plot(m.locfit2)
abline(v=1, col="red")

### But GAM misses again
m.gam2 <- gam(inside ~ s(r, k=50), data=df, family=binomial)
plot(m.gam2, se=FALSE, rug=FALSE, trans=plogis)
abline(v=1, col="red")

### Can also plot gam on link scale for alternate view
plot(m.gam2, se=FALSE, rug=FALSE)
abline(v=1, col="red")

gam.check(m.gam2)

Method: UBRE   Optimizer: outer newton
full convergence after 4 iterations.
Gradient range [-3.29203e-08,-3.29203e-08]
(score -0.8240065 & scale 1).
Hessian positive definite, eigenvalue range [7.290233e-05,7.290233e-05].

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

         k'    edf k-index p-value
s(r) 49.000 10.537   0.979    0.06

Answer 1

Estás ignorando el intercepto del modelo al evaluar el ajuste del mismo. El plot muestra el spline ajustado, pero el modelo incluye un término constante paramétrico, al igual que el intercepto en un modelo de regresión logística estándar.

En su lugar, prediga a partir del modelo ajustado utilizando el predict() para ubicaciones en una cuadrícula de ubicaciones sobre el intervalo. Por ejemplo:

m.gam <- gam(inside ~ te(x, y), data=df, family=binomial, method = "REML")
locs <- with(df,
             data.frame(x = seq(min(x), max(x), length = 100),
                        y = seq(min(y), max(y), length = 100)))
pred <- expand.grid(locs)
pred <- transform(pred,
                  fitted = predict(m.gam, newdata = pred, type = "response"))
contour(locs$x, locs$y, matrix(pred$fitted, ncol = 100))
draw.circle(0, 0, 1, border="red")

que da

Utilizando un te() más suave parece hacer un poco mejor que s() y utilicé method = "REML" ya que esto puede ayudar en situaciones en las que la función objetivo en la selección basada en GCV/UBRE puede llegar a ser plana (y por lo tanto estos métodos pueden ser poco suaves), en caso de que ese fuera el problema aquí.