Modelos aditivos generalizados (GAM)
Marly Valencia
2024-05-15
#suppressMessages(install.packages("gamair"))
#install.packages("mgcv")
#install.packages("splines")
#install.packages("knitr")
#install.packages("broom")
#install.packages("kableExtra")
#install.packages("ISLR")
#install.packages("GGally")
#install.packages("glmtoolbox")
#install.packages("plotly")
suppressMessages(suppressWarnings(library(plotly)))
suppressMessages(suppressWarnings(library(glmtoolbox)))
suppressMessages(suppressWarnings(library(GGally)))
suppressMessages(suppressWarnings(library(readr)))
suppressMessages(suppressWarnings(library(ggplot2)))
suppressMessages(suppressWarnings(library(mgcv)))
suppressMessages(suppressWarnings(library(gamair)))
suppressMessages(suppressWarnings(library(splines)))
suppressMessages(suppressWarnings(library(ISLR)))
suppressMessages(suppressWarnings(library(MASS)))
suppressMessages(suppressWarnings(library(dplyr)))MODELOS ADITIVOS GENERALIZADOS (GAM)
Comunmente las relaciones entre diferentes variables no presentan una relación lineal, en este sentido, los modelos aditivos generalizados (GAM, por sus siglas en inglés) son modelos que permiten ajustar estas relaciones cuando no se tienen no son lineales. Tal como se muestra en la figura 1, en la cual una regresión lineal simple no se ajusta a la relación de los datos. Al respecto, la diferencia con el modelo lineal general y los modelos lineales generalizados se tiene en la especificación del predictor lineal:
\[ \begin{gathered} y \sim \text { LEF }(\mu, \phi/\omega_i) \\ E(y)=\mu \\ g(\mu)=b_0+f_1\left(x_1\right)+f_2\left(x_2\right) \ldots+f_i\left(x_p\right) \end{gathered} \]
Donde \(g(\mu)\) es la función de enlace y en el predictor lineal se emplean funciones de suavizado que se presentan como \(f_1 \ldots f_i\) para cada una de las \(x_1 \ldots x_p\) covariables.
set.seed(124)
x <- seq(0,10, length.out = 500)
y <- sin(x) + cos(4 * x) + rnorm(500, sd = 0.2)
df <- data.frame(x = x, y = y)
# Crear el gráfico
ggplot(df, aes(x = x, y = y)) +
#geom_line(color = "purple") +
geom_point(color = "red", size = 1) +
geom_smooth(method = "lm", se = FALSE, color = "black", linetype = "dashed") + #
labs(title = "Modelo lineal general (Fig.1)",
x = "X",
y = "Y") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
Modelo con una sola covariable
El caso más sencillo es el modelo en el cual solo se tiene una covariable. Siguiendo a Wood(2006) la especificación es la siguiente:
\[y_i=f(x)+\epsilon_i\] Para estimar \(f(x)\) se puede definir de una forma esta de forma lineal donde \(b_i\) es una base elegida, lo cual es quivalente a elegir funciones de base conocidas.
\[ f(x)=\sum_{i=1}^q b_i(x) \beta_i \] Una de las formas más simples que puede tomar \(f(x)\) es la de un polinomio, por ejemplo, si este es de grado cinco la base para este espacio sería la siguiente: \(b_1(x)=1\), \(b_2(x)=x\),\(b_3(x)=x^2\),\(b_4(x)=x^3\),\(b_5(x)=x^4\), \(b_6(x)=x^5\) así:
\[f(x)=\beta_1+x\beta_2+x^2\beta_3+x^3\beta_4+x^4\beta_5+x^5\beta_6\]
Splines
Estos consisten en una curva comformada por trozos unidos por puntos (también conocidos en inglés como knots), los cuales pueden tomar alguna forma funcional con el objetivo de acercase lo máximo posible a la forma que toma la relación de los datos. Estos son utiles cuando el uso de un único polinomio no genera una buena aproximación.
A continuación, se puede ver la diferencia entre cada uno de los splines para el ejemplo anterior, lo que se observa es que entre mayor sea el grado, mayor es el ondulamiento.
plot(x, y, t="p", ylim=range(y)*c(1,1.2), main="Splines", col="gray")
fit<-smooth.spline(x, y, df=2)
fit5<-smooth.spline(x, y, df=8)
fit10<-smooth.spline(x, y, df=20)
lines(fit, col="red")
lines(fit5, col="green")
lines(fit10, col="orange")
Ejemplo: Spline cúbico
Como se puede observar en la figura 2, se tiene una curva conformada por partes polinomios cúbicos. Al respecto, cada curva es continua hasta su segunda derivada:
¿Cómo encontrar el grado adecuado?
Para ello se fija la dimensión de una base en un tamaño mayor que un valor que se creería razonablemente necesario y se emplea un término de penalización, de esta forma, el problema de minimización se reduce a:
\[ \|\mathbf{y}-\mathbf{X} \boldsymbol{\beta}\|^2+\lambda \int_0^1\left[f^{\prime \prime}(x)\right]^2 d x \]
Donde,
\[ \int_0^1\left[f^{\prime \prime}(x)\right]^2 d x=\boldsymbol{\beta}^{\top} \mathbf{S} \boldsymbol{\beta}, \]
Este segundo término penaliza los modelos que son muy “ondulantes” y el equilibrio entre la suavidad y ajuste del modelo se controla mediante el parámetro de suavizado $$, de esta forma, la estimación del grado de suavidad del modelo consiste en estimar este parámetro. En este sentido, el estimador de minimos cuadrados penalizados es:
\[ \hat\beta=(X^{\top}X+\lambda S)^{-1}X^{\top}y \]
Estimación penalizada
Mientras que en un modelo lineal generalizado se maximiza la logverosimilitud, aquí se realiza el mismo procedimiento, añadiendo el término de penalización, de tal forma que:
\[l_p(\beta)=l(\beta)- penalización \] \[l_p(\beta)=l(\beta)- penalización \]
Estimación con una y más variables
Para una variable tenemos una sola función y una base para la covariable, por lo cual, la estimación estará dada por un GLM con la siguiente estimación:
\[l_p(\beta)=l(\beta)-\frac{1}{2} \lambda \boldsymbol{\beta}^{\top} \mathbf{S} \boldsymbol{\beta} \]
Para cuando se tiene más de una variable (caso multivariado) dado que se tienen covariables y por ende más funciones:
\[l_p(\beta)=l(\beta)- \frac{1}{2} \sum_j \lambda_j \boldsymbol{\beta}^{\top} \mathbf{S}_j \boldsymbol{\beta} \]
APLICACIÓN
Ajustando un GAM para los datos simulados tenemos lo siguiente con la función gam del paquete mgcv:
gamaj<-gam(y~s(x), data=df)
summary(gamaj)##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.18866 0.03303 5.712 1.93e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(x) 7.128 8.174 45.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.431 Deviance explained = 43.9%
## GCV = 0.55439 Scale est. = 0.54537 n = 500plot(gamaj)
A diferencia de una regresión por medio de un modelo lineal generalizado, aquí se observa una estimación del intercepto y por otro lado, se tiene información sobre los términos de suavizado. Al respecto, el primer término hace referencia a los grados de libertad los cuales indican el ondulamiento, es decir, un grado de libertad igual a 1 indica que se tiene una relación lineal, y un edf de 2 equivale a una curva cuadrática.
En las siguientes dos columnas se presentan dos estadisticos (grado de libertad de referencia y F) para elaborar el p-valor que se encuentra al final, por medio del cual se realiza una prueba de significancia del coeficiente de suavizado.
Ej: Retinopatía en pacientes diabéticos
Empleando los datos disponibles en la libreria gamair sobre pacientes diabéticos que desarrollaron o no retinopatía, la cual es una enfermedad que afecta en gran parte a estos pacientes generando daños en los vasos sanguíneos de la retina ocasionando así ceguera o pérdida de la visión. Entre las variables de la base de datos wesdr se tiene el índice de masa corporal (bmi), porcentaje de hemoglobina unida a la glucosa de la sangre (gly) y duración de la enfermedad (dur).
A continuación se presenta el correlograma de los datos, en el cual se observa que no hay una relación lineal entre las variables.
# Crear el correlograma
data("wesdr")
ggpairs(wesdr,
columns = 1:3, # Especifica las columnas a usar
labs(title = "Correlograma de las variables del conjunto de datos wesdr"))
Se implementan dos modelos, un GLM y un GAM empleando una distribución binomial dado que la variable de interés tiene respuesta binaria (presencia o no de retinopatía).
data(wesdr)
modelo2 <- gam(ret ~ s(dur)+ s(gly)+s(bmi),
family = binomial(link = logit), data =wesdr, method = "ML")
glm_model2=glm(ret ~ dur+ gly+bmi,
family = binomial(link = logit), data = wesdr)Se calcula el AIC y se realiza una prueba de análisis de varianza (ANOVA) para comparar los dos modelos. Lo que se observa es que el modelo GAM presenta una menor devianza y a su vez un menor AIC.
#Comparar los dos modelos
anova(modelo2, glm_model2, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: ret ~ s(dur) + s(gly) + s(bmi)
## Model 2: ret ~ dur + gly + bmi
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 657.28 747.27
## 2 669.00 780.98 -11.724 -33.713 0.0006337 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(modelo2, glm_model2)## df AIC
## modelo2 10.07012 767.4112
## glm_model2 4.00000 788.9838summary(modelo2)##
## Family: binomial
## Link function: logit
##
## Formula:
## ret ~ s(dur) + s(gly) + s(bmi)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.41768 0.08909 -4.688 2.76e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(dur) 3.665 4.573 15.78 0.00589 **
## s(gly) 1.000 1.000 90.17 < 2e-16 ***
## s(bmi) 2.752 3.497 13.94 0.00505 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.216 Deviance explained = 17.7%
## -ML = 384 Scale est. = 1 n = 669plot(modelo2)
Empleando un Step criterion en el modelo de glm se llega a que se debe ajustar uno con las variables gly y bmi.
stepCriterion(glm_model2, criterion ="aic")##
## Family: binomial
## Link function: logit
## Criterion: AIC
##
## Initial model:
## ~ 1
##
##
## Step 0 :
## df AIC BIC adj.R-squared P(Chisq>)(*)
## + gly 1 795.34 804.36 0.1274 < 2e-16
## + bmi 1 909.45 918.46 0.0016 0.09556
## <none> 910.25 914.76 0.0000
## + dur 1 912.10 921.11 -0.0013 0.69633
##
## Step 1 : + gly
##
## df AIC BIC adj.R-squared P(Chisq>)(*)
## + bmi 1 787.55 801.06 0.1369 0.002234
## <none> 795.34 804.36 0.1274
## + dur 1 797.31 810.82 0.1261 0.845635
##
## Step 2 : + bmi
##
## df AIC BIC adj.R-squared P(Chisq>)(*)
## <none> 787.55 801.06 0.1369
## + dur 1 788.98 807.01 0.1362 0.454
## - gly 1 909.45 918.46 0.0016 <2e-16
##
## Final model:
## ~ gly + bmi
##
## ****************************************************************************
## (*) p-values of the Wald testSe ajustan nuevamente los dos modelos para compararlos y se calcula nuevamente el AIC, en el cual el modelo GAM sigue siendo mejor.
modelo3 <- gam(ret ~ s(gly)+s(bmi),
family = binomial(link = logit), data =wesdr, method = "ML")
glm_model3=glm(ret ~ gly+bmi,
family = binomial(link = logit), data = wesdr)
AIC(modelo3, glm_model3)## df AIC
## modelo3 5.891956 777.1440
## glm_model3 3.000000 787.5469vis.gam(modelo3,type = 'response',plot.type = 'persp',phi = 30,theta= 30,
n.grid= 500,border= NA)
par(mfrow=c(2,2))
gam.check(modelo3, type = "deviance", old.style = FALSE)
##
## Method: ML Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-9.696999e-05,4.155108e-07]
## (score 387.5268 & scale 1).
## Hessian positive definite, eigenvalue range [9.695475e-05,0.856485].
## 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(gly) 9.00 1.00 0.96 0.17
## s(bmi) 9.00 3.11 1.00 0.55Concurvidad
Esta se puede entender o ver como una generalización o una extensión de lo que en los modelos lineales se entiende por multicolinealidad, por lo tanto, existe concurvidad si la función \(f_i\) asociada a la variable \(i\) es una combinación de la funciones de las otras variables del modelo. Al respecto, se emplea la función concurvity, la cual presenta tres medidas diferentes de concurvidad para cada una de las variables, donde cero indica no concurvidad y uno que la función es una combinación de una o más funciones de otras variables.
En el ejemplo, se observa que la concurvidad es baja para cada uno de los términos en el modelo.
concurvity(modelo3)## para s(gly) s(bmi)
## worst 6.530498e-23 0.04329072 0.04329072
## observed 6.530498e-23 0.02103373 0.01933658
## estimate 6.530498e-23 0.01709244 0.01544374CONCLUSIONES
Los modelos aditivos generalizados presentan un mejor ajuste para aquellos escenarios en los cuales las variables no tienen relación lineal, sin embargo, la interpretación de los resultados no es tan clara como lo es en un modelo lineal general o lineal generalizado.
data(wesdr)
modelop <- gam(ret ~ s(dur)+ s(gly)+s(bmi),
family = binomial(link = logit), data =wesdr, method = "ML")
glm_modelp=glm(ret ~ dur+ gly+bmi,
family = binomial(link = logit), data = wesdr)
AIC(modelo2, glm_model2)## df AIC
## modelo2 10.07012 767.4112
## glm_model2 4.00000 788.9838