6.9 Modelos lineales generalizados
Como ya se comentó, los modelos lineales generalizados son una extensión de los modelos lineales para el caso de que la distribución condicional de la variable respuesta no sea normal, introduciendo una función de enlace (o link) \(g\) de forma que \[g\left(E(Y | \mathbf{X} )\right) = \beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{p}X_{p}\] y su ajuste en la práctica se realiza empleando el método de máxima verosimilitud (habrá que especificar también una familia de distribuciones para la respuesta).
La función link debe ser invertible, de forma que se pueda volver a transformar el modelo ajustado (en la escala lineal de las puntuaciones) a la escala original. Por ejemplo, como se comentó al final de la Sección 1.2.1, para modelar una variable indicadora, con distribución de Bernouilli (caso particular de la Binomial) donde \(E(Y | \mathbf{X} ) = p(\mathbf{X})\) es la probabilidad de éxito, podemos considerar la función logit \[\operatorname{logit}(p(\mathbf{X}))=\log\left( \frac{p(\mathbf{X})}{1-p(\mathbf{X})} \right) = \beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{p}X_{p}\] (que proyecta el intervalo \([0, 1]\) en \(\mathbb{R}\)), siendo su inversa la función logística \[p(\mathbf{X}) = \frac{e^{\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{p}X_{p}}}{1 + e^{\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{p}X_{p}}}\] Esto da lugar al modelo de regresión logística (múltiple), que será el que utilizaremos como ejemplo en esta sección. Para un tratamiento más completo de los métodos de regresión lineal generalizada se recomienda consultar alguno de los libros de referencia, como Faraway (2016), T. J. Hastie y Pregibon (2017), Dunn y Smyth (2018) o McCullagh (2019).
Para el ajuste (estimación de los parámetros) de un modelo lineal generalizado a un conjunto de datos (por máxima verosimilitud) se emplea la función glm() (la mayoría de los principales parámetros coinciden con los de la función lm()):
ajuste <- glm(formula, family = gaussian, data, weights, subset, na.action, ...)El parámetro family especifica la distribución y opcionalmente la función de enlace.
Por ejemplo:
gaussian(link = "identity"),gaussian(link = "log")binomial(link = "logit"),binomial(link = "probit")poisson(link = "log")Gamma(link = "inverse")
Para cada distribución se toma por defecto una función de enlace (el denominado enlace canónico, mostrada en primer lugar en la lista anterior; ver help(family) para más detalles).
Por ejemplo, en el caso del modelo logístico bastará con establecer family = binomial.
También se podría emplear la función bigglm() del paquete biglm para ajustar modelos lineales generalizados a grandes conjuntos de datos, aunque en este caso los requerimientos computacionales pueden ser mayores.
Como se comentó en la Sección 6.1, muchas de las herramientas y funciones genéricas disponibles para los modelos lineales son válidas también para este tipo de modelos, como por ejemplo las mostradas en las tablas 6.1 y 6.2.
Como ejemplo continuaremos con los datos de clientes de la compañía de distribución industrial HBAT, pero consideraremos como respuesta la variable alianza y como predictores las percepciones de HBAT (al igual que en las secciones anteriores consideraremos únicamente variables explicativas continuas, sin interacciones, por comodidad).
df <- hbat[c(24, 7:19)]
set.seed(1)
nobs <- nrow(df)
itrain <- sample(nobs, 0.8 * nobs)
train <- df[itrain, ]
test <- df[-itrain, ]En primer lugar se suele realizar un análisis descriptivo. Por ejemplo, si el número de variables no es muy grande, podemos generar un gráfico de dispersión matricial, diferenciando las observaciones pertenecientes a las distintas clases (ver Figura 6.21).
plot(train[-1], pch = as.numeric(train$alianza), col = as.numeric(train$alianza))
Figura 6.21: Gráfico de dispersión matricial, con colores y símbolos dependiendo de alianza.
Para ajustar un modelo de regresión logística bastaría con establecer el argumento family = binomial en la llamada a glm() (por defecto utiliza link = "logit"):
modelo <- glm(alianza ~ velocida + calidadp, family = binomial, data = train)
modelo##
## Call: glm(formula = alianza ~ velocida + calidadp, family = binomial,
## data = train)
##
## Coefficients:
## (Intercept) velocida calidadp
## -12.5218 1.6475 0.7207
##
## Degrees of Freedom: 159 Total (i.e. Null); 157 Residual
## Null Deviance: 218.2
## Residual Deviance: 160.5 AIC: 166.5La razón de ventajas (OR) permite cuantificar el efecto de las variables explicativas en la respuesta (incremento proporcional en la razón entre la probabilidad de éxito y la de fracaso, al aumentar una unidad la variable manteniendo las demás fijas):
exp(coef(modelo)) # Razones de ventajas ("odds ratios")## (Intercept) velocida calidadp
## 3.646214e-06 5.194162e+00 2.055887e+00exp(confint(modelo))## 2.5 % 97.5 %
## (Intercept) 4.465945e-08 1.593277e-04
## velocida 2.766629e+00 1.068554e+01
## calidadp 1.557441e+00 2.789897e+00Para obtener un resumen más completo del ajuste también se utiliza summary():
summary(modelo)##
## Call:
## glm(formula = alianza ~ velocida + calidadp, family = binomial,
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8273 -0.7622 -0.2998 0.7837 1.8375
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.5218 2.0758 -6.032 1.62e-09 ***
## velocida 1.6475 0.3426 4.809 1.52e-06 ***
## calidadp 0.7207 0.1479 4.872 1.11e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 218.19 on 159 degrees of freedom
## Residual deviance: 160.55 on 157 degrees of freedom
## AIC: 166.55
##
## Number of Fisher Scoring iterations: 5La desvianza (deviance) es una medida de la bondad del ajuste de un modelo lineal generalizado (sería equivalente a la suma de cuadrados residual de un modelo lineal; valores más altos indican peor ajuste). La Null deviance se correspondería con un modelo solo con la constante y la Residual deviance con el modelo ajustado. En este caso hay una reducción de 57.65 con una pérdida de 2 grados de libertad (una reducción significativa).
Para contrastar globalmente el efecto de las covariables también podemos emplear:
modelo.null <- glm(alianza ~ 1, binomial, train)
anova(modelo.null, modelo, test = "Chi")## Analysis of Deviance Table
##
## Model 1: alianza ~ 1
## Model 2: alianza ~ velocida + calidadp
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 159 218.19
## 2 157 160.55 2 57.646 3.036e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.9.1 Selección de variables explicativas
El objetivo sería conseguir un buen ajuste con el menor número de variables explicativas posible.
Al igual que en el caso del modelo de regresión lineal múltiple, se podría seguir un proceso interactivo, eliminando o añadiendo variables con la función update(), aunque también están disponibles métodos automáticos de selección de variables.
Para obtener el modelo “óptimo” lo ideal sería evaluar todos los modelos posibles.
En este caso no se puede emplear la función regsubsets del paquete leaps (sólo para modelos lineales),
pero por ejemplo el paquete
bestglm
proporciona una herramienta equivalente (bestglm()).
También se podría emplear la función stepwise() del paquete RcmdrMisc para seleccionar un modelo por pasos según criterio AIC o BIC:
# library(RcmdrMisc)
modelo.completo <- glm(alianza ~ ., family = binomial, data = train)
modelo <- stepwise(modelo.completo, direction='forward/backward', criterion='BIC')##
## Direction: forward/backward
## Criterion: BIC
##
## Start: AIC=223.27
## alianza ~ 1
##
## Df Deviance AIC
## + velocida 1 189.38 199.53
## + calidadp 1 192.15 202.30
## + facturac 1 193.45 203.60
## + producto 1 196.91 207.06
## + quejas 1 198.10 208.25
## + imgfvent 1 198.80 208.95
## + web 1 204.40 214.55
## + publi 1 209.28 219.43
## + precio 1 211.97 222.12
## + garantia 1 212.37 222.52
## <none> 218.19 223.27
## + nprod 1 213.97 224.12
## + soporte 1 216.50 226.65
## + flexprec 1 216.99 227.14
##
## Step: AIC=199.53
## alianza ~ velocida
##
## Df Deviance AIC
## + calidadp 1 160.55 175.77
## + imgfvent 1 178.43 193.65
## + web 1 181.52 196.74
## + precio 1 183.38 198.61
## <none> 189.38 199.53
## + flexprec 1 185.47 200.69
## + producto 1 185.54 200.77
## + nprod 1 186.92 202.14
## + facturac 1 186.93 202.15
## + garantia 1 187.02 202.25
## + publi 1 187.44 202.67
## + soporte 1 189.06 204.29
## + quejas 1 189.27 204.50
## - velocida 1 218.19 223.27
##
## Step: AIC=175.77
## alianza ~ velocida + calidadp
##
## Df Deviance AIC
## + imgfvent 1 137.05 157.35
## + web 1 145.63 165.93
## <none> 160.55 175.77
## + publi 1 156.03 176.33
## + facturac 1 157.35 177.66
## + flexprec 1 157.44 177.74
## + producto 1 157.75 178.06
## + garantia 1 158.97 179.27
## + nprod 1 160.18 180.47
## + quejas 1 160.20 180.50
## + soporte 1 160.37 180.67
## + precio 1 160.37 180.67
## - calidadp 1 189.38 199.53
## - velocida 1 192.15 202.30
##
## Step: AIC=157.35
## alianza ~ velocida + calidadp + imgfvent
##
## Df Deviance AIC
## <none> 137.05 157.35
## + precio 1 134.97 160.35
## + flexprec 1 135.39 160.77
## + publi 1 135.65 161.03
## + producto 1 135.72 161.09
## + facturac 1 135.81 161.19
## + garantia 1 136.31 161.69
## + nprod 1 136.63 162.00
## + soporte 1 136.79 162.16
## + quejas 1 136.96 162.33
## + web 1 137.03 162.40
## - velocida 1 160.34 175.57
## - imgfvent 1 160.55 175.77
## - calidadp 1 178.43 193.65summary(modelo)##
## Call:
## glm(formula = alianza ~ velocida + calidadp + imgfvent, family = binomial,
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6527 -0.6822 -0.1482 0.7631 2.0561
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.5164 3.4593 -5.931 3.02e-09 ***
## velocida 1.6631 0.3981 4.177 2.95e-05 ***
## calidadp 1.0469 0.2014 5.197 2.02e-07 ***
## imgfvent 1.0085 0.2398 4.205 2.61e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 218.19 on 159 degrees of freedom
## Residual deviance: 137.05 on 156 degrees of freedom
## AIC: 145.05
##
## Number of Fisher Scoring iterations: 56.9.2 Análisis e interpretación del modelo
Las hipótesis estructurales del modelo son similares al caso de regresión lineal (aunque algunas como la linealidad se suponen en la escala transformada). Si no se verifican, los resultados basados en la teoría estadística pueden no ser fiables o totalmente erróneos.
Con el método plot() se pueden generar gráficos de interés para la diagnosis del modelo (ver Figura 6.22):
oldpar <- par( mfrow=c(2,2))
plot(modelo)
Figura 6.22: Gráficos de diagnóstico del ajuste lineal generalizado.
par(oldpar)Su interpretación es similar a la de los modelos lineales (consultar las referencias incluidas al principio de la sección para más detalles).
En este caso destacan tres posibles datos atípicos, aunque aparentemente no son muy influyentes a posteriori (en el modelo ajustado).
Adicionalmente se pueden generar gráficos parciales de residuos, por ejemplo con la función crPlots() del paquete car (ver Figura 6.23):
# library(car)
id <- list(method = which(abs(residuals(modelo, type = "pearson")) > 3), col = 2)
crPlots(modelo, id = id, main = "")
Figura 6.23: Gráficos parciales de residuos del ajuste generalizado.
Se pueden emplear las mismas funciones vistas en los modelos lineales para obtener medidas de diagnosis de interés (tablas 6.1 y 6.2). Por ejemplo residuals(model, type = "deviance") proporcionará los residuos deviance.
Por supuesto también pueden aparecer problemas de colinealidad, y podemos emplear las mismas herramientas para detectarla:
# library(car)
vif(modelo)## velocida calidadp imgfvent
## 1.193557 1.656649 1.451237Si no se satisfacen los supuestos básicos también se pueden intentar distintas alternativas (se puede cambiar la función de enlace y la familia de distribuciones, que puede incluir parámetros para modelar dispersión, además de las descritas en la Sección 6.4), incluyendo emplear modelos más flexibles o técnicas de aprendizaje estadístico que no dependan de ellas (sustancialmente).
6.9.3 Evaluación de la precisión
Para evaluar la calidad de la predicción en nuevas observaciones podemos seguir los pasos mostrados en la Sección 1.3.5.
Obteniendo las estimaciones de la probabilidad (de la segunda categoría) empleando predict() con type = "response":
p.est <- predict(modelo, type = "response", newdata = test)
pred <- factor(p.est > 0.5, labels = c("No", "Si")) # levels = c('FALSE', 'TRUE')y las medidas de precisión de la predicción (además de los criterios AIC o BIC tradicionales):
caret::confusionMatrix(pred, test$alianza, positive = "Si", mode = "everything")## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 19 5
## Si 3 13
##
## Accuracy : 0.8
## 95% CI : (0.6435, 0.9095)
## No Information Rate : 0.55
## P-Value [Acc > NIR] : 0.0008833
##
## Kappa : 0.5918
##
## Mcnemar's Test P-Value : 0.7236736
##
## Sensitivity : 0.7222
## Specificity : 0.8636
## Pos Pred Value : 0.8125
## Neg Pred Value : 0.7917
## Precision : 0.8125
## Recall : 0.7222
## F1 : 0.7647
## Prevalence : 0.4500
## Detection Rate : 0.3250
## Detection Prevalence : 0.4000
## Balanced Accuracy : 0.7929
##
## 'Positive' Class : Si
## o de las estimaciones de las probabilidades (como el AUC).
6.9.4 Extensiones
Se pueden imponer restricciones a las estimaciones de los parámetros de modo análogo al caso de modelos lineales (secciones 6.7 y 6.8). Por ejemplo, en los métodos de regularización (ridge, LASSO o elastic net; Sección 6.7) bastaría con cambiar en la función de pérdidas la suma residual de cuadrados por el logaritmo negativo de la función de verosimilitud.
Ejercicio 6.1 Emplear el paquete glmnet para ajustar modelos logísticos con penalización ridge y LASSO a la muestra de entrenamiento de los datos de clientes de la compañía de distribución industrial HBAT, considerando como respuesta la variable alianza y seleccionando un valor “óptimo” del hiperparámetro \(\lambda\).
Ajustar también un modelo con penalización elastic net empleando caret (seleccionando los valores óptimos de los hiperparámetros).
El método PCR (Sección 6.8) se extendería de forma inmediata al caso de modelos generalizados, simplemente cambiando el modelo ajustado. También están disponibles métodos PLSR para modelos generalizados.
Ejercicio 6.2 Emplear el paquete caret para ajustar modelos logísticos con reducción de la dimensión a los datos de clientes de la compañía de distribución industrial HBAT. Comparar el modelo obtenido con preprocesado "pca" y el método "glmStepAIC", con el obtenido empleando el método "plsRglm".