GAM

Hoy he dado una charla en la Carlos III. En la comida me han preguntado, algo extrañados, por un ejemplo que había enseñado en el que ajustaba un modelo usando GAMs.

El motivo era que quienes preguntaban —que trabajan con ese tipo de modelos— encuentran muy difícil, se ve, convencer a otros usuarios de los métodos estadísticos (economistas, etc.) de adoptarlos. Yo he contestado que hace unos pocos días a unos primíparos que acababan de ajustar sus tres primeros lms con R les invité a probar GAMs con sus datos. ¿Por qué no?

La situación que teníamos planteada era parecida aunque un poco más complicada que esta:

irradiation <- read.table("http://rredc.nrel.gov/solar/old_data/nsrdb/1961-1990/hourly/1990/13994_90.txt", 
                          header = F, skip = 1)
 
irradiation <- irradiation[,4:5]
names(irradiation) <- c("hour", "irradiation")
 
plot(irradiation$hour, irradiation$irradiation, xlab = "hora", ylab = "irradiación")

que produce

irradiacion00

Obviamente, la irradiación solar (en San Luis, Misuri, durante el año 1990) depende grandemente de la hora del día.

Si uno prueba lo de siempre, lo que enseñan al principio, lo que parece que solo le está permitido a uno si no consigue cinturón negro en las artes del chamanismo estadístico, obtiene

mi.lm <- lm(irradiation ~ hour, data = irradiation)
 
summary(mi.lm)
 
# Call:
#   lm(formula = irradiation ~ hour, data = irradiation)
# 
# Residuals:
#   Min     1Q Median     3Q    Max 
# -341.2 -336.9 -260.4  336.0  930.4 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 335.4736     9.1814  36.538   <2e-16 ***
#   hour          0.2408     0.6426   0.375    0.708    
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 416.3 on 8758 degrees of freedom
# Multiple R-squared:  1.604e-05,	Adjusted R-squared:  -9.814e-05 
# F-statistic: 0.1405 on 1 and 8758 DF,  p-value: 0.7078
 
plot(irradiation$hour, irradiation$irradiation, xlab = "hora", ylab = "irradiación")
abline(mi.lm, col = "red")

El gráfico obtenido, que muestra la predicción del modelo en rojo es

irradiacion01

En resumen, la hora del día, aunque manifiestamente relacionada con la predicción, tiene un coeficiente nulo según los manuales al uso y quedaría fuera del estudio según los más de ellos. ¡Fijaos en ese p-valor tan tirrioso!

Contraintiutivo, ¿verdad?

¿Por qué no probar entonces con un modelo que recoja el efecto no lineal de la hora? En el menú tampoco hay tantas opciones y con suerte uno desemboca muy naturalmente en

library(mgcv)
 
mi.gam <- gam(irradiation ~ s(hour, bs = "cr"), data = irradiation)
 
summary(mi.gam)
 
# Family: gaussian 
# Link function: identity 
# 
# Formula:
#   irradiation ~ s(hour, bs = "cr")
# 
# Parametric coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  338.484      1.647   205.5   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Approximate significance of smooth terms:
#   edf Ref.df    F p-value    
# s(hour) 8.919  8.998 6122  <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# R-sq.(adj) =  0.863   Deviance explained = 86.3%
# GCV =  23799  Scale est. = 23772     n = 8760
 
plot(mi.gam, xlab = "hora", ylab = "irradiación")

que, ahora sí, produce

irradiacion02

El resto de la discusión es poco relevante para los que tenemos prisa.

Y las notas:

  • Sí, podemos en este caso usar la hora como variable categórica (¿pero qué si tenemos minutos, etc.?)
  • Efectivamente, hay otras opciones en el menú. ¡Pero nos gusta la A de GAM!
  • Quienes no descabalgan de los (g)lms puede que igual se dediquen a analizar datos demasiado tontos; o puede que igual no sepan lo que hacen… ¡pobrecicos!
  • Quienes se quejan de que los del punto anterior no descabalgan tal vez no hayan sabido dar con el conjunto de datos —les regalo el de hoy— para convencerlos.

9 comentarios sobre “GAM

  1. ffernandez 13 noviembre, 2015 11:38

    Si por lo menos aceptasen los glm en toda su extensión… Un día discutí con un cliente durante horas porque metí una interacción en una regresión porque «metía multicolinealidad en el modelo»… reducir el error de predicción a la mitad no parecía ser argumento suficientemente convincente…

  2. daniel 13 noviembre, 2015 20:16

    No se si la pregunta es un poco tonta, pero el gráfico parece una parábola y por tanto está claro que nos e puede aproximar bien mediante una recta, y claramente es preciso utilizar un modelo lineal que tenga en cuenta la hora al cuadrado. En fin, que no me parece justo descalificar a los modelos lineales usando una recta cuando se ve gráficamente que se necesita tener en cuenta el cuadrado de la variable explicativa, sería algo así como descalificar a todos los polinomios diciendo que sólo existen las rectas.

  3. Carlos J. Gil Bellosta 13 noviembre, 2015 20:23

    Hummm… tienes los datos: ¡prueba! Y a ver qué parábola ajustas que tenga una «base» plana (la correspondiente a las horas de la noche). De todos modos, los modelos GAM son lineales: la A significa «additive». Solo que viene precedida de una G de «generalized».

    Los polinomios no están mal en general. Pero tienen un problema: se van al infitino por los bordes. Si fuesen suficientes para todo fin, nadie se habría preocupado por desarrollar teoría que los dejase atrás para ciertas aplicaciones.

  4. David 13 noviembre, 2015 21:34

    En este ejemplo tan sencillo propongo una alternativa que permite una mayor interpretabilidad del modelo basada en la segmentación del modelo lineal original (aunque el ajuste a los datos es más discutible que en el caso del GAM):

    library(segmented)
    mi.seglm <- segmented(mi.lm, ~hour, psi=c(5, 12, 20)) #psi son la aproximación incial a los puntos de cambio de tendencia
    summary(mi.seglm) #Los coeficientes son mucho más interpretables que el GAM
    plot(mi.seglm) #El ajuste es bueno, aunque los cambios tan bruscos de pendiente en los puntos de cambio de tendencia obviamente son mucho menos creíbles que los cambios suaves del GAM

  5. daniel 14 noviembre, 2015 0:19

    Seleccionando el intervalo de horas entre las 5 y las 20 horas para eliminar la zona plana obtenemos un coeficiente de regresión cercano al 80 por ciento que es peor que el otro modelo pero no muy malo. Se ha ajustado un polinomio de cuarto grado, i = irradiation para abreviar, h = hora,

    j =5 & i$h summary(lm(i ~ h + I(h^2) + I(h^3) + I(h^4) + I(h^5), j))

    Call:
    lm(formula = i ~ h + I(h^2) + I(h^3) + I(h^4) + I(h^5), data = j)

    Residuals:
    Min 1Q Median 3Q Max
    -367.53 -134.42 -8.56 164.39 296.43

    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 3.767e+03 4.258e+02 8.847 < 2e-16 ***
    h -1.832e+03 2.072e+02 -8.842 < 2e-16 ***
    I(h^2) 2.967e+02 3.804e+01 7.799 7.36e-15 ***
    I(h^3) -1.796e+01 3.315e+00 -5.417 6.29e-08 ***
    I(h^4) 3.666e-01 1.379e-01 2.659 0.00785 **
    I(h^5) -2.475e-04 2.202e-03 -0.112 0.91051

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Residual standard error: 188.8 on 5834 degrees of freedom
    Multiple R-squared: 0.7953, Adjusted R-squared: 0.7951
    F-statistic: 4532 on 5 and 5834 DF, p-value: < 2.2e-16

  6. Alanítico 15 noviembre, 2015 14:41

    Sin ánimo de crear una discusión metodológica, parece que algunos perdéis perspectiva, esto no va de aplicar un modelo para un conjunto de datos, esto va de explicar los datos. La modelización de los datos no es el fin del estudio, ese es un medio para interpretarlos.

    Es poco inteligente correr todos los modelos habidos y por haber (como pollo sin cabeza) para encontrar el que tenga el p-valor o el R^2 o el estadístico que sea de la manera que más nos guste, antes que nada bien vale la pena saber para qué sirve cada modelo y en qué situaciones resulta más útil su uso, no se busca adaptar o pulir el resultado a un algoritmo, una metodología o un modelo por capricho onanista, se busca explicar de manera sencilla los DATOS y si no se hace por lo menos que parezca que no fue deliberado o en todo caso como aquí nos presenta el Sr. Gil Bellosta sin más ambición que explicar el funcionamiento de GAM. Yo no usaría GAM para los datos de irradiación solar, pero aún menos haría OLS, eso es comulgar con ruedas de molino. ¿Por qué no GAM? Porque la estructura de los datos no es TAN propicia como con otras alternativas y porque es un capricho (aunque es un capricho más honesto que el modelo lineal). Qué es sensato en este caso, la solución más sencilla renunciar a estimaciones, sabemos que lo que tenemos entre manos es una https://es.wikipedia.org/wiki/Proyecci%C3%B3n_cil%C3%ADndrica Lo más sensato es usar trigonometría ¿es más costoso? solo si nos hemos vuelto tan gandules como para no ver más allá del modelo lineal ¿Qué vamos a hacer calcular eclipses al 95% de confianza? ¿No nos gusta, estamos atados a los estadísticos o es preciso hacerlo de otra manera? Entonces series temporales (tampoco es la panacea), qué duda cabe que la rotación se trata de un movimiento armónico.

    Como colofón a este chorizo infumable que os he soltado, toca el sermón, de verdad preguntaos qué sentido tiene lo que estáis haciendo, ¿no os estáis saltando los principios más elementales del método científico por empeñaros en hacer algo de determinada manera? ¿que buscáis segmentando los datos para explicar la irradiación, me vais a decir que se trata de una variable cualitativa, una dummy? por favor… por supuesto que hay irradiación cuando hay sol, solo faltaría. Usar el modelo lineal es renunciar de antemano a explicar las estaciones, la declinación solar… y dejarlo todo a merced de la media. Si tienes un reloj roto estate seguro de que va a dar la hora correcta dos veces al día, cuando usas el modelo lineal sucede algo parecido, es el argumento que minimiza los errores y punto, solo eso (mirad el primer gráfico del artículo y pensad que solo habéis trazado una línea en mitad de esa maraña), el modelo lineal sirve cuando hay pocos datos, cuando hay mucho ruido y cuando no sabemos absolutamente nada de los datos, cuando nos enfrentamos a ellos por primera vez es un buen punto de partida (probablemente el mejor), pero quedarse ahí también es una renuncia expresa a indagar y a querer saber más. Sin nada más que objetar y sin el más mínimo ánimo de ofender: un saludo y a estrujarse más la cabeza.

  7. Carlos J. Gil Bellosta 15 noviembre, 2015 15:52

    Gracias por el comentario, que es mayormente acertado. Solo que creo que responde en parte a cuestiones distintas de las que aquí se quieren tratar.

    Una nota previa: hay dos motivos por los que se construyen modelos: uno, para interpretar datos (como apuntas); otra, para predecir. Hay casos totalmente legítimos de unos, otros y, cómo no, intermedios. A veces la interpretación es más importante que la predicción; otras, ocurre lo contrario.

    Ahora, al grano. Es cierto que la irradiación teórica en un ahora determinada se puede calcular usando cálculos… digamos que astronómicos. En un modelo serio, es lo que habría que hacer. De hecho, existen tablas que la proporcionan a futuro. Otra cosa es la medida en superficie, que suele ser más interesante: existen las nubes, etc. En ese caso hay que ir más allá. Un ejemplo/problema más próximo a la realidad sería el de estimar la irradiación real en superficie a una hora concreta dados los valores de irradiación teóricos y las predicciones meteorológicas. Existen muchos modelos estadísticos que pueden emplearse para tal fin. Los mencionados en la entrada (y muchos más) valdrían.

    Quise ilustrar el uso de los GAM con unos datos en los que la relación entre una variable importante y la objetivo fuese manifiestamente no lineal. Elegí esos por motivos diversos (entre los que hay que contar la falta de imaginación). El objetivo era exclusivamente mostrar que hay vida más allá de los modelos más tradicionales de los que son, de hecho, generalizaciones. Generalizaciones que han aparecido, precisamente, para afrontar los problemas que crean variables como las que aparecen en el conjunto de datos en cuestión.

    Y, finalmente, no es que aquí esté tratando de saltarme los principios del método científico para hacer algo de una determinada manera. Estoy simplemente ilustrando otra manera, otra técnica, tal vez poco conocida, tanto para aplicar a rajatabla los principios del método científico como para lo contrario. Queda eso a la entera discreción del usuario, como con todo tipo de técnica.

  8. Alanítico 15 noviembre, 2015 19:07

    No tengo queja alguna contra su artículo ni pongo en duda su buena intención, entendí desde el primer momento que el propósito era ilustrar el funcionamiento de GAM (creo que lo digo en el tostón de arriba) como un recurso más con el que contar, pero cuando vi los comentarios empeñándose en utilizar el modelo lineal se me cayó el alma al suelo.

    Dicho lo cual, fue deliberadamente por eso por lo que quise dar un paso más, saltarme la temática del blog y traer a colación que todo se puede hacer de otra forma que a veces ni nos lo planteamos, a lo mejor me pasé de frenada.

    Sobre el método lo que reprocho es que en lugar de modelar se moldee, que se pretenda reproducir la forma de una función sin inferir lo más mínimo su significado como si esculpiésemos el David de Miguel Ángel, cuenta la leyenda que el artista le pegó un martillado implorándole que hablase, bueno pues aquí igual, he visto Frankensteins sin pies ni cabeza, y la gente tan orgullosa. Que un polinomio de grado 5 te explique el 80 % no es casualidad, pero ¿significa algo?¿vale en todas las regiones del planeta? Y aunque la respuesta fuese sí, ¿de verdad el proceso ha sido honesto o solo estás buscando un R^2 alto?

  9. David 15 noviembre, 2015 20:06

    No sé si parte de los comentarios de Alanítico iban referidos hacia mi aproximación de segmentar el modelo lineal: «¿que buscáis segmentando los datos para explicar la irradiación, me vais a decir que se trata de una variable cualitativa, una dummy? »

    En cualquier caso, aclarar que la segmentación que he propuesto en mi comentario no consiste en categorizar la variable predictora sino dividir la recta de regresión en varios tramos en los que se le permite tener pendientes distintas, de manera muy similar a lo que hace GAM. La ventaja de este método es que en cada tramo seguimos teniendo una regresión lineal y, por lo tanto, los resultados son fácilmente interpretables:

    Estimated Break-Point(s):
    Est. St.Err
    psi1.hour 5.66 0.05876
    psi2.hour 12.51 0.02620
    psi3.hour 19.41 0.05912

    t value for the gap-variable(s) V: 0 0 0

    Meaningful coefficients of the linear terms:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) -2.671 8.666 -0.308 0.758
    hour 1.335 2.613 0.511 0.609
    U1.hour 169.273 3.044 55.608 NA
    U2.hour -339.736 2.208 -153.840 NA
    U3.hour 167.127 3.044 54.903 NA

    Residual standard error: 157.9 on 8752 degrees of freedom
    Multiple R-Squared: 0.8563, Adjusted R-squared: 0.8562

Los comentarios están desabilitados.