Parte práctica
MLG
2024-09-22
Modelos Aditivos Generalizados
Descripción de los datos
Librerías a utilizar
suppressMessages(suppressWarnings(library(glmtoolbox)))
suppressMessages(suppressWarnings(library(GLMsData)))
suppressMessages(suppressWarnings(library(aplore3)))
suppressMessages(suppressWarnings(library(FSAdata)))
suppressMessages(suppressWarnings(library(ISLR)))
suppressMessages(suppressWarnings(library(plotly)))
suppressMessages(suppressWarnings(library(leaps)))
suppressMessages(suppressWarnings(library(gamair)))
suppressMessages(suppressWarnings(library(mgcv)))
suppressMessages(suppressWarnings(library(psych)))
data(wesdr)Lectura de los datos
Los datos, originarios del paquete gss, registran si los pacientes diabéticos desarrollaron o no retinopatía junto con tres posibles predictores. Cuenta con 669 registros.
El marco de datos wesdr tiene las siguientes columnas:
ret Variable binaria: 1 = retinopatía, 0 = no (variable respuesta).
bmi Índice de masa corporal (peso en kg dividido por el cuadrado de la altura en metros).
gly Hemoglobina glicosilada: porcentaje de hemoglobina unida al glucosa en la sangre. Esto refleja los niveles medios de glucosa en sangre a largo plazo: menos del 6% es típico de los no diabéticos, pero rara vez lo alcanzan los pacientes diabéticos.
dur Duración de la enfermedad en años.
Así entonces nuestra variable de interes: \[Y_i \sim Bernoulli(\pi). \]
data(wesdr)
w <- wesdr
head(w)## dur gly bmi ret
## 1 10.3 13.7 23.8 0
## 2 9.9 13.5 23.5 0
## 3 15.6 13.8 24.8 0
## 4 26.0 13.0 21.6 1
## 5 13.8 11.1 24.6 1
## 6 31.1 11.3 24.6 1summary(w)## dur gly bmi ret
## Min. : 1.20 Min. : 6.00 Min. :14.40 Min. :0.0000
## 1st Qu.: 5.40 1st Qu.:10.70 1st Qu.:20.60 1st Qu.:0.0000
## Median : 9.20 Median :12.20 Median :22.90 Median :0.0000
## Mean :11.46 Mean :12.49 Mean :23.16 Mean :0.4155
## 3rd Qu.:15.30 3rd Qu.:14.10 3rd Qu.:25.20 3rd Qu.:1.0000
## Max. :55.20 Max. :22.60 Max. :50.80 Max. :1.0000cor.plot(cor(w[,1:3]),
main="Mapa de Calor",
diag=T,
show.legend = T,numbers=F,upper=F) 
library(ggcorrplot)
library(corrr)
model.matrix(~0+., data=w[,1:3]) %>%
cor(use="pairwise.complete.obs") %>%
ggcorrplot(show.diag = F, type="lower", lab=TRUE, lab_size=2)
library(GGally)
# Crear la matriz de gráficos de dispersión
w[,4]<-as.factor(w[,4])
ggpairs(w, columns=1:3, ggplot2::aes(colour=w[,4]))
A través de estas gráficas, observamos que la autocorrelación entre las variables regresoras es mínima. La distribución de las variables parece estar sesgada hacia la derecha, siendo este sesgo más marcado en las variables de duración e índice de masa corporal.
Selección de variables
Definimos un modelo lineal generalizado de la familia binomial con función de enlace logit donde seincluyen todas las posibles interacciones. Se utiliza la metología fordward y backward para la selección de las variables signfifcativas utilizando el criterio AIC por un lado y por otro lado el criterio de BIC.
attach(w)
m <- ret ~ (bmi + gly + dur)**2
fit <- glm(m, family=binomial(logit), data=wesdr)
stepCriterion(fit, direction="forward", criterion="aic", test="wald")##
## 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>)(*)
## + bmi:gly 1 785.23 803.25 0.1404 0.04303
## <none> 787.55 801.06 0.1369
## + dur 1 788.98 807.01 0.1362 0.45405
## - gly 1 909.45 918.46 0.0016 < 2e-16
##
## Step 3 : + bmi:gly
##
## df AIC BIC adj.R-squared P(Chisq>)(*)
## <none> 785.23 803.25 0.1404
## + dur 1 786.58 809.11 0.1398 0.4219
##
## Final model:
## ~ gly + bmi + gly:bmi
##
## ****************************************************************************
## (*) p-values of the Wald teststepCriterion(fit, direction="backward", criterion="aic", test="wald")##
## Family: binomial
## Link function: logit
## Criterion: AIC
##
## Initial model:
## ~ (bmi + gly + dur)^2
##
##
## Step 0 :
## df AIC BIC adj.R-squared P(Chisq>)(*)
## - bmi:dur 1 788.08 815.11 0.1391 0.75000
## - dur:gly 1 788.53 815.57 0.1386 0.44948
## <none> 789.98 821.52 0.1379
## - bmi:gly 1 792.88 819.92 0.1338 0.03163
##
## Step 1 : - bmi:dur
##
## df AIC BIC adj.R-squared P(Chisq>)(*)
## - dur:gly 1 786.58 809.11 0.1398 0.47333
## <none> 788.08 815.11 0.1391
## - bmi:gly 1 790.89 813.41 0.1351 0.03343
##
## Step 2 : - dur:gly
##
## df AIC BIC adj.R-squared P(Chisq>)(*)
## - dur 1 785.23 803.25 0.1404 0.42192
## <none> 786.58 809.11 0.1398
## + bmi:dur 1 788.53 815.57 0.1386 0.82757
## - bmi:gly 1 788.98 807.01 0.1362 0.04087
##
## Step 3 : - dur
##
## df AIC BIC adj.R-squared P(Chisq>)(*)
## <none> 785.23 803.25 0.1404
## - bmi:gly 1 787.55 801.06 0.1369 0.04303
##
## Final model:
## ~ bmi + gly + bmi:gly
##
## ****************************************************************************
## (*) p-values of the Wald testConforme a los resultados obtenidos, podemos evidenciar que el modelo final obtenido incluye las regresoras indice de masa corporal, glicemia y la interacción entre ellas.
Conforme a ello, se procede a realizar el ajuste del modelo final obtenido:
Ajuste del modelo
mod <- ret ~ bmi + gly + bmi:gly
# Ajustar el modelo de regresión logística
fit <- glm(mod, family=binomial(logit), data=w)
# Resumen del modelo ajustado
summary(fit)##
## Call:
## glm(formula = mod, family = binomial(logit), data = w)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.83653 2.97888 -0.281 0.779
## bmi -0.19762 0.13110 -1.507 0.132
## gly -0.08565 0.23478 -0.365 0.715
## bmi:gly 0.02116 0.01046 2.023 0.043 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 908.25 on 668 degrees of freedom
## Residual deviance: 777.23 on 665 degrees of freedom
## AIC: 785.23
##
## Number of Fisher Scoring iterations: 4MODELO ADITIVO GENERALIZADO (GAM)
De acuerdo con lo observado, no hay relación lineal entre los pares de variables, por lo tanto, se emplea un modelo GAM:
fit2 <- gam(ret ~ s(gly)+s(bmi)+ti(gly,bmi),
family = binomial(link = logit), data =wesdr, method = "ML")
summary(fit2)##
## Family: binomial
## Link function: logit
##
## Formula:
## ret ~ s(gly) + s(bmi) + ti(gly, bmi)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3971 0.0894 -4.441 8.94e-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(gly) 1.000 1.000 89.037 < 2e-16 ***
## s(bmi) 2.956 3.754 17.676 0.00129 **
## ti(gly,bmi) 1.959 2.425 6.196 0.07559 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.2 Deviance explained = 16.6%
## -ML = 385.57 Scale est. = 1 n = 669fit3 <- gam(ret ~ gly+s(bmi)+ti(gly,bmi),
family = binomial(link = logit), data =wesdr, method = "ML")
summary(fit3)##
## Family: binomial
## Link function: logit
##
## Formula:
## ret ~ gly + s(bmi) + ti(gly, bmi)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.31585 0.53268 -9.979 <2e-16 ***
## gly 0.39370 0.04172 9.436 <2e-16 ***
## ---
## 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(bmi) 2.956 3.754 17.676 0.00129 **
## ti(gly,bmi) 1.959 2.426 6.197 0.07562 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.2 Deviance explained = 16.6%
## -ML = 385.57 Scale est. = 1 n = 669MODELO GAM (libreria gamlss)
Ahora realizaremos el ajuste del modelo GAM pero usando la libreria “gamlss”.
suppressMessages(suppressWarnings(library(gamlss)))
fit4<-gamlss(formula = ret~pb(bmi)+pb(gly)+bmi:gly,family = BI, data = wesdr, trace = FALSE)
summary(fit4)## ******************************************************************
## Family: c("BI", "Binomial")
##
## Call: gamlss(formula = ret ~ pb(bmi) + pb(gly) + bmi:gly,
## family = BI, data = wesdr, trace = FALSE)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: logit
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.47813 3.02317 -0.820 0.413
## pb(bmi) -0.12731 0.13146 -0.968 0.333
## pb(gly) 0.04931 0.23219 0.212 0.832
## bmi:gly 0.01531 0.01026 1.493 0.136
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 669
## Degrees of Freedom for the fit: 6.082421
## Residual Deg. of Freedom: 662.9176
## at cycle: 3
##
## Global Deviance: 763.7219
## AIC: 775.8867
## SBC: 803.2928
## ******************************************************************Análisis de residuales
par(mfrow=c(2,2))
gam.check(fit2, type = "deviance", old.style = FALSE)
##
## Method: ML Optimizer: outer newton
## full convergence after 9 iterations.
## Gradient range [-9.826147e-05,4.904481e-06]
## (score 385.5746 & scale 1).
## Hessian positive definite, eigenvalue range [2.27553e-05,0.8915548].
## Model rank = 35 / 35
##
## 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.13
## s(bmi) 9.00 2.96 1.00 0.50
## ti(gly,bmi) 16.00 1.96 0.96 0.10plot(fit4)
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = 0.002338991
## variance = 1.046788
## coef. of skewness = 0.02119696
## coef. of kurtosis = 2.998924
## Filliben correlation coefficient = 0.9990314
## ******************************************************************COMPARACIÓN DE MODELOS
Comparamos los dos modelos empleando el AIC y BIC, con lo cual se observa que el modelo GAM se ajusta mejor:
data.frame(Model = c("fit", "fit2", "fit3", "fit4"), AIC = AIC(fit, fit2, fit3, fit4)[,2], BIC = BIC(fit, fit2, fit3, fit4)[,2])## Model AIC BIC
## 1 fit 785.2287 803.2518
## 2 fit2 774.2140 811.0650
## 3 fit3 774.2151 811.0698
## 4 fit4 775.8867 803.2928