Modelos Aditivos Generalizados

Where the wild things are

(slides en verde claro pueden obviarse)

Modelos Aditivos Generalizados

\[g(\mu_i) = \eta(X_{i1}, ..., X_{ik})\] GLM \[g(\mu_i) = \beta_{0} + \beta_{1}X_{1i} + ... +\beta_{k}X_{ki}\]

GAM \[g(\mu_i) = \beta_{0} + f_{1}X_{1i} + ... +f_{k}X_{ki}\]

Primero: polinomios por piezas

Basis functions

¿Qué colocamos en f? la idea es tener una familia de funciones o trasformaciones que puedan ser aplicadas a las variables X: \(f_{1}X_{1i} + f_{2}X_{2} ... +f_{k}X_{ki}\). En vez de ajustar un modelo lineal en X, ajustamos el modelo incluyendo estas funciones.
Ejemplo: picewise cubic spline polynomial

\[y_i = \left\{\begin{matrix} \beta_{01} + \beta_{11}x_{i} + \beta_{21}x_{i}^2 + \beta_{31}x_{i}^3 + \epsilon_i \; si \: x_i < \xi \\ \beta_{02} + \beta_{12}x_{i} + \beta_{22}x_{i}^2 + \beta_{32}x_{i}^3 + \epsilon_i \; si \: x_i \geq \xi \end{matrix}\right.\]

\(\xi\) es un nudo (knot)

claramente es horrible

Existen requisitos a cumplir por las funciones básicas

  • Debe ser continuo.
  • Debe ser suave.

En la práctica esto se logra imponiendo contricciones sobre la continuidad de la primera y la segunda derivada de la curva. El resultado de estas contricciones es un cubic spline.
Existen muchas funciones básicas, veremos los splines porque están entre las más "populares" y para entender su lógica.

basic truncated functions in cubic splines

\(y_i = \beta_0 + \beta_1b_1(x_i) + \beta_2b_2(x_i) + ... + \beta_{k+3}b_{k+3}(x_i) + \epsilon_i\)
\(b_1 = x_i\)
\(b_2 = x_i^2\)
\(b_3 = x_i^3\)
\(b_{k+3} = (x_i - \xi)_+^3 \;\; k = 1, 2... K\)
\((x_i - \xi)_+^3 = \left\{\begin{matrix} (x_i - \xi)^3 \quad if \quad x > \xi \\ 0 \quad otherwise \end{matrix}\right.\)

basic truncated functions in cubic splines

con un knot, \(x_i < \xi\)
\(y_i = \beta_0 + \beta_1x_i+\beta_2x_i^2 + \beta_3x_i^3\)
\(\frac{\partial y}{\partial x} = \beta_1 + 2\beta_2x + 3\beta_3x^2\)

cuando \(x_i = \xi\)
\(y_\xi = \beta_0 + \beta_1\xi+\beta_2\xi^2 + \beta_3\xi^3 + \beta_4(\xi-\xi)^3\)
\(\frac{\partial y}{\partial x} = \beta_1 + 2\beta_2x + 3\beta_3x^2 + 3\beta_4(\xi - \xi)^2\)

cuando \(x_i > \xi\)
\(y_i = \beta_0 + \beta_1x_i+\beta_2x_i^2 + \beta_3x_i^3+ \beta_4(x_i-\xi)^3\)
\(\frac{\partial y}{\partial x} = \beta_1 + 2\beta_2x + 3\beta_3x^2 + 3\beta_4(x_i - \xi)^2\)

las bases truncadas aseguran que las funciones sean continuas (que empiecen donde terminó la anterior) y que esta condición se cumpla también en la primera (y segunda) derivada, originando una curva suave.

Segundo: Smoothing splines

¿Cómo determinar la complejidad de la curva que ajustamos? Si usamos más nudos (y por lo tanto, más grados de libertad), la curva se va a ajustar cada vez mejor a los datos, pero va a resultar un modelo demasiado complejo.

Bias-variance trade-off

En general, los modelos se ajustan minimizando algún tipo de medida del error (por ejemplo, los mínimos cuadrados.

\[1/n \sum(y_i - f(x_i))^2 \]

Pero esta medida decrece a medida que se incrementa la flexibilidad del modelo (medida de su complejidad, como los grados de libertad).
¿Qué pasa con un modelo complejo si luego intentamos utilizarlo para predecir?


negro, d.f. = 2; rojo d.f. = 3; azul d.f. = 10.

Esto sucede porque estamos utilizando los residuos del set de DATOS DE ENTRENAMIENTO (training data).

Si usamos los residuso sobre DATOS DE PRUEBA (test data)

El valor esperado de test MSE

\[E(MSE_{test}) = Var(\hat f(x_0)) + [Bias(\hat f(x_0))]^2 + Var(\epsilon)\] La varianza surge de los cambios en \(\hat f\) con diferentes datos de entrenamiento. Métodos más flexibles tienen más varianza.
El sesgo es el error introducido por tratar de aproximar el comportamiento de la vida real con un modelo. Modelos más sencillos tienen más sesgo.
Este trade-off es la base de la PENALIZACIÓN.


* de James et al. 2014

Volviendo a los smoothing splines…

minimizar \[\sum(y_i - f(x_i))^2 + \lambda\int f''(t)^2 d(t)\] donde el segúndo término es una penalización contra la variabilidad en \(f\) (\(f''\) es la segunda derivada, el cambio en la pendiente). \(\lambda\) es un parámetro de suavizado.

  • Si \(\lambda = 0\) la penalidad se anula y el modelo puede ser tan complejo como sea necesario para que los residuos tengan el mínimo valor posible (en la práctica, puede volver los residuos iguales a cero).
  • Si \(\lambda \rightarrow \infty\) la penalidad se vuelve muy grande y el modelo debe ser tan simple como sea posible.
  • Si un spline es minimizado de esta manera se vuelve un smoothing spline.
  • \(\lambda\) controla el trade-off entre sesgo y varianza.
  • \(\lambda\) controla los efective degreees of freedom.

Selección de \(\lambda\). Cross Validation.

Ordinary Cross Validation (OCV) - Leave One Out Cross Validation (LOOCV) Para cada valor de \(\lambda\) en una serie se calculan los errores de la predicción como… \[RSS_{\lambda} = \sum (y_i - f_{\lambda}^{[-i]}(x_i))^2\]

Donde el -i indica que se excluye la observación i para calcular el residuo respecto a \(y-i\).

A tener en cuenta…
- Existe una simplificación que evita hacer n ajustes.
- AIC es asintóticamente LOOCV. BIC es asintóticamente (para cierto k) leave-k-cross validation.
- GCV (generalized cross validation) es una versión invariante de OCV.

Selección de \(\lambda\). ML y REML.

Métodos basados en verosimilitud
- Los métodos basados en los errores de la predicción tienen cierta tendencia a encontrar más de un mínimo local.
- Los métodos basados en ML o REML tratan las funciones de suavizado como efectos aleatorios, de tal forma que \(\lambda\) puede ser tomado como un parámatro de varianza.
- Esto surge de una interpretación Bayesiana de los GAMs, ya que estamos imponiendo un prior de que las funciones son suaves.
- (sin fórmulas aquí)

Selección de \(\lambda\)

Tercero: GAMs en R

(veremos gran parte de los detalles en el ejercicio de ejemplo en el práctico)

Paquete mgcv de Simon Wood

  • Amplia variedad de splines.
  • Selección de variables por shrinkage.
  • Soporta métodos de verosimilitud y de validación cruzada para obtener \(\lambda\).
  • Tiene una función especial (bam) para grandes sets de datos.
  • Omite las observaciones sin datos.
  • Soporta suavizados multidimensionales utilizando tensores y thin plate splines.
  • Posee herramientas de diagnóstico generales y de diagnóstico de la concurvidad.

library(mgcv)
fit <- gam(y ~ s(x), method = "REML")
summary(fit)
Family: gaussian 
Link function: identity 

Formula:
y ~ s(x)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.01746    0.04851    0.36     0.72

Approximate significance of smooth terms:
      edf Ref.df     F p-value    
s(x) 5.45  6.588 32.88  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.777   Deviance explained = 79.6%
-REML = 38.474  Scale est. = 0.14827   n = 63

  • Los valores P en el summary (y en anova) son aproximados y han cambiado según la versión de mgcv.

gam.check(fit)
Method: REML   Optimizer: outer newton
full convergence after 5 iterations.
Gradient range [-2.152141e-05,1.58881e-05]
(score 38.47406 & scale 0.1482711).
Hessian positive definite, eigenvalue range [2.207732,30.67498].
Model rank =  10 / 10 

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) 9.00 5.45    1.05    0.65

END