mgcv: la base p-spline no tiene soporte local

Question

Entiendo que una de las propiedades atractivas de la base B-spline es que las funciones base individuales tienen soporte local, es decir, será >0 en el intervalo entre $d+2$ nudos adyacentes (dejando $d$ sea el grado de la base del spline), pero cero en el resto. Sin embargo, los p-spines utilizados en mgcv no parecen tener esta propiedad.

library(MASS)
library(mgcv)
gamfit <- gam(accel ~ s(times, bs="ps"), data=mcycle)
mc_ord <- mcycle[order(mcycle$times), ]
X <- predict(gamfit, mc_ord, type='lpmatrix')
matplot(y=X, x=mc_ord$times, type='l')

Después de extraer la proyección de la base del spline utilizando predict(gamfit, mc_ord, type='lpmatrix') y trazando las columnas contra la variable original no proyectada, estoy viendo que los splines usados en mgcv (usando bs="ps" ) no parecen ser no negativos.

Simon Wood es un tipo súper inteligente, así que asumo que estoy mal informado/confundido/ignorante de alguna manera. Puede alguien ayudarme a estar menos desinformado/confuso/ignorante?

Mis principales preguntas:

1. Después de buscar en Google, me parece que los diferentes tipos de estrías son confusos. ¿Existe una nomenclatura estándar? ¿Cómo llamamos a este tipo de spline / variante de B-spline?
2. ¿Cuál es la ventaja de utilizar este tipo de función de base spline frente a una que sea distinta de cero?
3. ¿Hay alguna manera de exigir mgcv para utilizar splines con soporte local?

Answer 1

Es debido a las restricciones de la suma a cero. Si sólo necesita suavizar una sola vez en su modelo hay un par de opciones...

library(mgcv)
## set up p-spline smoother without constraint
sm <- smoothCon(s(times,bs="ps",k=20),data=mcycle)[[1]]
X <- sm$X ## no constraint => compact support
S <- sm$S[[1]]
accel <- mcycle$accel
m <- gam(accel~X-1,paraPen=list(X=list(S)))
plot(mcycle)
lines(mcycle$times,fitted(m))

## or fool gam into not using a constraint... 
n <- nrow(mcycle)
## dum is basically a constant, but just variable enough for gam to 
## not impose a constraint... 
dum <- rep(1,n) + (runif(n)-.5) * 1e-8
m <- gam(accel~s(times,bs="ps",k=20,by=dum)-1,data=mcycle)
plot(m,scheme=1,residuals=T,pch=19,cex=.5,by.resids=T)
model.matrix(m)

Si necesita varios alisados, puede clonar los constructores de alisados p-spline y establecer restricciones que mantengan un soporte compacto - el problema es que tienden a conducir a CIs muy amplios ya que se confunde a través de la intercepción...

Simon