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.