7.3 Modelos aditivos

Se supone que: \[Y= \beta_{0} + f_1(X_1) + f_2(X_2) + \ldots + f_p(X_p) + \varepsilon\] con $f_{i},$ $i=1,...,p,$ funciones cualesquiera. De esta forma se consigue mucha mayor flexibilidad que con los modelos lineales pero manteniendo la interpretabilidad de los efectos de los predictores. Adicionalmente se puede considerar una función de enlace, obteniéndose los denominados modelos aditivos generalizados (GAM). Para más detalles sobre este tipo modelos ver por ejemplo T. Hastie y Tibshirani (1990) o Wood (2017).

Los modelos lineales (generalizados) serían un caso particular considerando $f_{i}(x) = \beta_{i}x$. Además, se podrían considerar cualquiera de los métodos de suavizado descritos anteriormente para construir las componentes no paramétricas (por ejemplo si se emplean splines naturales de regresión el ajuste se reduciría al de un modelo lineal). Se podrían considerar distintas aproximaciones para el modelado de cada componente (modelos semiparamétricos) y realizar el ajuste mediante backfitting (se ajusta cada componente de forma iterativa, empleando los residuos obtenidos al mantener las demás fijas). Si en las componentes no paramétricas se emplea únicamente splines de regresión (con o sin penalización), se puede reformular el modelo como un GLM (regularizado si hay penalización) y ajustarlo fácilmente adaptando herramientas disponibles (penalized re-weighted iterative least squares, PIRLS).

De entre todos los paquetes de R que implementan estos modelos destacan:

gam: Admite splines de suavizado (univariantes, s()) y regresión polinómica local (multivariante, lo()), pero no dispone de un método para la selección automática de los parámetros de suavizado (se podría emplear un criterio por pasos para la selección de componentes). Sigue la referencia T. Hastie y Tibshirani (1990).
mgcv: Admite una gran variedad de splines de regresión y splines penalizados (s(); por defecto emplea thin plate regression splines penalizados multivariantes), con la opción de selección automática de los parámetros de suavizado mediante distintos criterios. Además de que se podría emplear un método por pasos, permite la selección de componentes mediante regularización. Al ser más completo que el anterior sería el recomendado en la mayoría de los casos (ver ?mgcv::mgcv.package para una introducción al paquete). Sigue la referencia Wood (2017).

La función gam() del paquete mgcv permite ajustar modelos aditivos generalizados empleando suavizado mediante splines:

ajuste <- gam(formula, family = gaussian, data, method = "GCV.Cp", select = FALSE, ...)

(también dispone de la función bam() para el ajuste de estos modelos a grandes conjuntos de datos y de la función gamm() para el ajuste de modelos aditivos generalizados mixtos, incluyendo dependencia en los errores). El modelo se establece a partir de la formula empleando s() para especificar las componentes “suaves” (ver help(s) y Sección 7.3.3).

Algunas posibilidades de uso son las que siguen:

Modelo lineal:
```
ajuste <- gam(y ~ x1 + x2 + x3)
```
Modelo (semiparamétrico) aditivo con efectos no paramétricos para x1 y x2, y un efecto lineal para x3:
```
ajuste <- gam(y ~ s(x1) + s(x2) + x3)
```
Modelo no aditivo (con interacción):
```
ajuste <- gam(y ~ s(x1, x2))
```
Modelo (semiparamétrico) con distintas combinaciones :
```
ajuste <- gam(y ~ s(x1, x2) + s(x3) + x4)
```

En esta sección utilizaremos como ejemplo el conjunto de datos Prestige de la librería carData. Se tratará de explicar prestige (puntuación de ocupaciones obtenidas a partir de una encuesta) a partir de income (media de ingresos en la ocupación) y education (media de los años de educación).

library(mgcv)
data(Prestige, package = "carData")
modelo <- gam(prestige ~ s(income) + s(education), data = Prestige)
summary(modelo)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## prestige ~ s(income) + s(education)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  46.8333     0.6889   67.98   <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(income)    3.118  3.877 14.61  <2e-16 ***
## s(education) 3.177  3.952 38.78  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.836   Deviance explained = 84.7%
## GCV = 52.143  Scale est. = 48.414    n = 102

# coef(modelo) 
# El resultado es un modelo lineal en transformaciones de los predictores

En este caso el método plot() representa los efectos (parciales) estimados de cada predictor (ver Figura 7.8):

plot(modelo, shade = TRUE, pages = 1) # residuals = FALSE por defecto

Figura 7.8: Estimaciones de los efectos parciales de income (izquierda) y education (derecha).

En general se representa cada componente no paramétrica (salvo que se especifique all.terms = TRUE), incluyendo gráficos de contorno para el caso de componentes bivariantes (correspondientes a interacciones entre predictores).

Se dispone también de un método predict() para calcular las predicciones de la forma habitual (por defecto devuelve las correspondientes a las observaciones modelo$fitted.values y para nuevos datos hay que emplear el argumento newdata).

7.3.1 Superficies de predicción

En el caso bivariante, para representar las estimaciones (la superficie de predicción) obtenidas con el modelo se pueden utilizar las funciones persp() o versiones mejoradas como plot3D::persp3D(). Estas funciones requieren que los valores de entrada estén dispuestos en una rejilla bidimensional. Para generar esta rejilla se puede emplear la función expand.grid(x,y) que crea todas las combinaciones de los puntos dados en x e y (ver Figura 7.9):

inc <- with(Prestige, seq(min(income), max(income), len = 25))
ed <- with(Prestige, seq(min(education), max(education), len = 25))
newdata <- expand.grid(income = inc, education = ed)
# Representamos la rejilla
plot(income ~ education, Prestige, pch = 16)
abline(h = inc, v = ed, col = "grey")

Figura 7.9: Observaciones y rejilla de predicción (para los predictores education e income).

y usaríamos estos valores para obtener la superficie de predicción, que en este caso⁴⁹ representamos con la función plot3D::persp3D() (ver Figura 7.10):

# Se calculan las predicciones
pred <- predict(modelo, newdata)
# Se representan
pred <- matrix(pred, nrow = 25)
# persp(inc, ed, pred, theta = -40, phi = 30)
# contour(inc, ed, pred, xlab = "Income", ylab = "Education")
# filled.contour(inc, ed, pred, xlab = "Income", ylab = "Education", 
#                key.title = title("Prestige"))
plot3D::persp3D(inc, ed, pred, theta = -40, phi = 30, ticktype = "detailed",
                xlab = "Income", ylab = "Education", zlab = "Prestige")

Figura 7.10: Superficie de predicción obtenida con el modelo GAM.

Puede ser más cómodo emplear el paquete modelr (emplea gráficos ggplot2) para trabajar con modelos y predicciones.

7.3.2 Comparación y selección de modelos

Además de las medidas de bondad de ajuste como el coeficiente de determinación ajustado, también se puede emplear la función anova() para la comparación de modelos (y seleccionar las componentes por pasos de forma interactiva). Por ejemplo, viendo el gráfico de los efectos se podría pensar que el efecto de education podría ser lineal:

# plot(modelo)
modelo0 <- gam(prestige ~ s(income) + education, data = Prestige)
summary(modelo0)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## prestige ~ s(income) + education
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.2240     3.7323   1.132    0.261    
## education     3.9681     0.3412  11.630   <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(income) 3.58  4.441 13.6  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.825   Deviance explained = 83.3%
## GCV = 54.798  Scale est. = 51.8      n = 102

anova(modelo0, modelo, test="F")

## Analysis of Deviance Table
## 
## Model 1: prestige ~ s(income) + education
## Model 2: prestige ~ s(income) + s(education)
##   Resid. Df Resid. Dev     Df Deviance      F Pr(>F)  
## 1    95.559     4994.6                                
## 2    93.171     4585.0 2.3886   409.58 3.5418 0.0257 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso aceptaríamos que el modelo original es significativamente mejor.

Alternativamente, podríamos pensar que hay interacción:

modelo2 <- gam(prestige ~ s(income, education), data = Prestige)
summary(modelo2)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## prestige ~ s(income, education)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  46.8333     0.7138   65.61   <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(income,education) 4.94  6.303 75.41  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.824   Deviance explained = 83.3%
## GCV = 55.188  Scale est. = 51.974    n = 102

En este caso el coeficiente de determinación ajustado es menor y ya no tendría sentido realizar el contraste.

Ademas se pueden seleccionar componentes del modelo (mediante regularización) empleando el parámetro select = TRUE. Para más detalles consultar la ayuda help(gam.selection) o ejecutar example(gam.selection).

7.3.3 Diagnosis del modelo

La función gam.check() realiza una diagnosis descriptiva y gráfica del modelo ajustado (ver Figura 7.11):

gam.check(modelo)

Figura 7.11: Gráficas de diagnóstico del modelo aditivo ajustado.

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 4 iterations.
## The RMS GCV score gradient at convergence was 9.783945e-05 .
## The Hessian was positive definite.
## Model rank =  19 / 19 
## 
## 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(income)    9.00 3.12    0.98    0.41
## s(education) 9.00 3.18    1.03    0.56

Lo ideal sería observar normalidad en los dos gráficos de la izquierda, falta de patrón en el superior derecho, y ajuste a una recta en el inferior derecho. En este caso parece que el modelo se comporta adecuadamente. Como se deduce del resultado anterior, podría ser recomendable modificar la dimensión k de la base utilizada construir la componente no paramétrica, este valor se puede interpretar como el grado máximo de libertad permitido en ese componente, aunque normalmente no influye demasiado en el resultado (puede influir en el tiempo de computación).

También se podría chequear concurvidad (generalización de la colinealidad; función concurvity()) entre las componentes del modelo:

concurvity(modelo)

##                  para s(income) s(education)
## worst    3.107241e-23 0.5931528    0.5931528
## observed 3.107241e-23 0.4065402    0.4398639
## estimate 3.107241e-23 0.3613674    0.4052251

Esta función devuelve tres medidas por componente, que tratan de medir la proporción de variación de esa componente que está contenida en el resto (similares al complementario de la tolerancia), un valor próximo a 1 indicaría que puede haber problemas de concurvidad.

También se puede ajustar modelos GAM empleando caret. Por ejemplo con los métodos "gam" y "gamLoess":

library(caret)
# names(getModelInfo("gam")) # 4 métodos
modelLookup("gam")

##   model parameter             label forReg forClass probModel
## 1   gam    select Feature Selection   TRUE     TRUE      TRUE
## 2   gam    method            Method   TRUE     TRUE      TRUE

modelLookup("gamLoess")

##      model parameter  label forReg forClass probModel
## 1 gamLoess      span   Span   TRUE     TRUE      TRUE
## 2 gamLoess    degree Degree   TRUE     TRUE      TRUE

Ejercicio 7.1 Continuando con los datos de MASS:mcycle, emplear mgcv::gam() para ajustar un spline penalizado para predecir accel a partir de times con las opciones por defecto y representar el ajuste obtenido. Comparar el ajuste con el obtenido empleando un spline penalizado adaptativo (bs="ad"; ver ?adaptive.smooth).

Ejercicio 7.2 Empleando el conjunto de datos airquality, crear una muestra de entrenamiento y otra de test, buscar un modelo aditivo que resulte adecuado para explicar sqrt(Ozone) a partir de Temp, Wind y Solar.R. ¿Es preferible suponer que hay una interacción entre Temp y Wind?

References

Hastie, T., y Tibshirani, R. (1990). Generalized Additive Models. Chapman; Hall. https://doi.org/10.1201/9780203753781

Wood, S. N. (2017). Generalized Additive Models: An Introduction with R, Second Edition. CRC Press.

Alternativamente se podrían emplear las funciones contour(), filled.contour(), plot3D::image2D() o similares.↩︎