Grados de libertad efectivos extraños (suavidad) seleccionados para el componente suave en el modelo GAM con mgcv

Question

Considere el siguiente ejemplo muy sencillo:

set.seed( 1 )
SimData <- data.frame( X = runif( 1000, 0, 1 ) )
SimData$Y <- rnbinom( nrow( SimData ), mu = 100*sin( SimData$X*2*pi )+100, size = 10 )
plot( gam( Y ~ s( X ), data = SimData, family = nb( link = log ) ) )

Aquí está el resultado:

Los grados de libertad efectivos son muy-muy cercanos a 9, lo que hace sospechar que la dimensión de la base por defecto no es lo suficientemente grande. Aumentémosla:

plot( gam( Y ~ s( X, k = 20 ), data = SimData, family = nb( link = log ) ) )

Hm. La imagen general es esencialmente la misma, pero la FED sugiere que la dimensión de la base todavía no es lo suficientemente grande.

Incluso si aumentamos k a 30, la FED sigue siendo mayor (24,8), por lo que, esencialmente, la FED parece seguir simplemente la k límite, lo cual es bastante extraño... (especialmente para una forma funcional tan simple, y sobre todo que ya estaba bien capturado por el modelo por defecto).

EDIT (07 Sep, 2017): Según una respuesta a la pregunta original, la aplicación de splines adaptativos ( bs="ad" ) puede ser la solución a este problema. Sin embargo...

Pongamos otro ejemplo sencillo:

set.seed( 1 )
SimData <- data.frame( X = runif( 1000, 0, 1 ) )
SimData$Y <- rnbinom( nrow( SimData ), mu = 100*sin( SimData$X*2*pi*3 )+1000, size = 10 )
plot( gam( Y ~ s( X ), data = SimData, family = nb( link = log ), method = "REML" ) )

¡Parece perfecto! Ahora vamos a "estropearlo":

SimData$Y[ SimData$X<=0.05|SimData$X>=0.95 ] <- 0

Esto da lugar al problema original: el FED es de 8,92 con el k 18,91 si k=20 , 46,74 si k=50 etc. A modo de ilustración, aquí está el k=50 caso:

Así que, efectivamente, tenemos el problema original. Intentemos por tanto bs="ad" :

plot( gam( Y ~ s( X, bs = "ad" ), data = SimData, family = nb( link = log ), method = "REML" ) )

Así que, por desgracia, el problema seguía siendo, incluso con bs="ad" ¡! (En realidad es la misma situación: aumentar k a 50 da un FED de 44,72, k=100 da 78,03. Por eso he decidido editar esta pregunta en lugar de empezar una nueva: parece que se trata en cierto modo de la misma historia...)

Answer 1

El ajuste en la escala logarítmica está causando el problema; en esa escala, como se muestra en su gráfico, el grado de suavidad no es realmente el mismo en toda la función; en valores bajos de x hay poca curvatura en la función verdadera, pero hay mucha curvatura alrededor del mínimo en ~ x = 7.5 . Como, por defecto, la ondulación de la spline está controlada por un único parámetro de suavidad, esto tiene el efecto de asumir el mismo grado de suavidad o curvatura en todas partes. Al aumentar k El modelo es capaz de adaptarse a los datos en torno a la subida y la bajada rápidas, pero al hacerlo también se adapta, porque también ha dado el parámetro de suavidad/penalización, al ruido de muestreo en otras partes de los datos.

Con estos datos he podido conseguir un ajuste aceptable cambiando la función de enlace por el enlace sqrt mediante family = nb(link = "sqrt") . En esta escala, las subidas y bajadas no son tan bruscas ni se limitan a un estrecho rango de x ya que está en la escala logarítmica.

m <-  gam(Y ~ s(X, k = 20), data = SimData, family = nb(link = "sqrt"), method = "REML")

La spline ajustada es:

Todavía hay algunas pruebas de adaptación al ruido de muestreo en la región donde x es baja, pero esto es bastante trivial y por cualquier métrica razonable sería indistinguible de la verdadera función, supongo.

La prueba de la dimensión de la base aplicada en gam.check() también parece razonablemente feliz con esto - nota que aumenté k a 20 pero esto parece suficiente para este modelo, mientras que k = 10 no lo era:

> gam.check(m)

Method: REML   Optimizer: outer newton
full convergence after 4 iterations.
Gradient range [2.520842e-06,0.0003287577]
(score 4439.901 & scale 1).
Hessian positive definite, eigenvalue range [5.53589,345.0859].
Model rank =  20 / 20 

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) 19.0 14.9    1.01    0.76

Los gráficos de diagnóstico también son bastante razonables:

Si se quiere ajustar en la escala logarítmica, se puede utilizar un spline adaptativo, que también parece dar buenos resultados para estos datos. Un spline adaptativo utiliza múltiples penalizaciones de suavidad (parámetros) para cada spline y, en efecto, es como ajustar varios splines separados a partes de los datos, cada uno con su propio parámetro de suavidad. Sin embargo, estos splines requieren muchos más datos para obtener buenos resultados.

mad <-  gam(Y ~ s(X, k = 20, bs = "ad"), data = SimData, family = nb, method = "REML")

La spline ajustada es ahora mucho más suave para valores bajos de x pero quizás haya un poco de sobreajuste al ruido de muestreo en grandes x . Sin embargo, en la escala de la respuesta esto va a ser trivial.

y gam.check() también se alegra de que k = 20 es suficiente para esta base/modelo/datos:

> gam.check(mad)

Method: REML   Optimizer: outer newton
full convergence after 11 iterations.
Gradient range [-0.0002868044,0.001512152]
(score 4461.377 & scale 1).
Hessian positive definite, eigenvalue range [0.0005701379,344.5224].
Model rank =  20 / 20 

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) 19.0 12.9    0.97    0.24

y de nuevo, los diagnósticos tampoco parecen tan malos

aunque hay un indicio de varianza no constante en los residuos en el gráfico de la parte superior derecha.

También añadiré que la prueba en gam.check() para la dimensión de la base no es infalible. Como mucho, es una guía. Si el modelo parece bueno en otros aspectos, yo pondría menos énfasis en los resultados de la prueba de dimensión. Sin embargo, para el modelo k = 10 que has ajustado, el modelo no es lo suficientemente flexible como para captar los bajos recuentos en torno a la caída y está bastante sesgado allí.

Yo diría que mientras sólo se aumenta k no ayudaba aquí, el modelo inicial tampoco era adecuado para estos datos.