En realidad, no hay funciones de base en el mismo sentido que las bases de spline de regresión de placa fina o de spline de regresión cúbica que mencionas. Las funciones base son funciones indicadoras binarias que devuelven un 1
si la observación pertenece a un grupo específico o 0
de lo contrario.
Para estos modelos, todo lo que tenemos que hacer es evaluar las funciones base en los valores observados de la(s) covariable(s) que forman la función suave: estos valores crean la matriz del modelo para la que se estimarán los coeficientes. También necesitamos una matriz de penalización $\mathbf{S}$ que cuando se escribe en forma cuadrática, $\boldsymbol{\beta}^{\mathsf{T}}\mathbf{S}\boldsymbol{\beta}$ da la penalización por ondulación deseada; por defecto, la matriz de penalización se construye de forma que penaliza la segunda derivada cuadrada integrada del spline. Esta forma es conveniente ya que toda la penalización puede escribirse como una función de la matriz de penalización y los coeficientes del modelo solamente. Los valores estimados para $\boldsymbol{\beta}$ que minimicen el criterio de logaritmo-verosimilitud penalizado.
Para un suavizador de efecto aleatorio, las columnas de la matriz del modelo que pertenecen al suavizador de efecto aleatorio son sólo un conjunto de vectores indicadores binarios, que indexan si una observación pertenece al $j$ nivel del factor/categoría. La matriz de penalización $\mathbf{S}$ es una matriz de identidad de tamaño $J$ (donde $J$ es el número de niveles del factor que forma el efecto aleatorio). Esto tiene el efecto de penalizar cada coeficiente a nivel de grupo para el efecto aleatorio en proporción a su desviación al cuadrado de 0. En otras palabras, los coeficientes a nivel de grupo para el efecto aleatorio se reducen hacia cero, tal y como lo harían si se ajustara un efecto aleatorio en un modelo de efectos mixtos.
Para ver esto, esto es lo que mgcv establece para el modelo simple con 1 efecto aleatorio suave y 1 regular suave:
library('mgcv')
set.seed(1)
dat <- gamSim(1, n = 400, scale = 2, verbose = FALSE) ## 4 term additive truth
## add a factor for the random effect and include it's effect in the simulated data
fac <- as.factor(sample(1:10, 400, replace = TRUE))
dat$X <- model.matrix(~ fac - 1)
b <- rnorm(10) * 0.5
dat <- transform(dat, y = y + X %*% b)
## fit the model
m <- gam(y ~ s(fac, bs = "re") + s(x0), data = dat, method = "ML")
En este ejemplo tenemos un factor de 10 niveles por lo que tenemos 10 coeficientes a nivel de grupo (para s(fac)
) más 9 funciones de base para el término TPRS regular (9 porque una de las funciones de base está confundida con el término de intercepción del modelo):
> coef(m)
(Intercept) s(fac).1 s(fac).2 s(fac).3 s(fac).4 s(fac).5
7.71917475 0.24762873 0.48122511 -0.17534482 0.03670067 -0.03800310
s(fac).6 s(fac).7 s(fac).8 s(fac).9 s(fac).10 s(x0).1
-0.52197990 0.53661911 -0.23833765 -0.26088376 -0.06762439 0.02740419
s(x0).2 s(x0).3 s(x0).4 s(x0).5 s(x0).6 s(x0).7
0.15183166 -0.06347520 -0.15041864 -0.04918018 0.15609539 -0.03564953
s(x0).8 s(x0).9
0.78256472 -0.26510925
cada una de las cuales corresponde a una función base que es una variable indicadora binaria en la matriz del modelo
> head(model.matrix(m)[, 2:11])
s(fac).1 s(fac).2 s(fac).3 s(fac).4 s(fac).5 s(fac).6 s(fac).7 s(fac).8
1 1 0 0 0 0 0 0 0
2 1 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0 1
5 1 0 0 0 0 0 0 0
6 0 0 0 0 0 1 0 0
s(fac).9 s(fac).10
1 0 0
2 0 0
3 0 1
4 0 0
5 0 0
6 0 0
La matriz de penalización para el s(fac, bs = 're')
efecto aleatorio suave es una matriz de identidad
> m[['smooth']][[1L]][['S']]
[[1]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0
[9,] 0 0 0 0 0 0 0 0 1 0
[10,] 0 0 0 0 0 0 0 0 0 1