04 Modelos Lineales Generalizados
Eduardo Villarreal
4/4/2021
1 Modelos Lineales Generalizados (GLM)
1.1 Caso 1: \(y \in [0, 1]\) es Binaria. Regresión Logística
Hasta ahora sólo hemos revisado el caso en el que \(y\) es una variable continua. Pero hay casos en los que nos interesa que la variable sea binaria: Por ejemplo, en calidad si un producto es rechazado o aceptado. En este caso aceptado = 1y rechazado = 0.
Para ilustrar el modelo que vamos a revisar, vamos a usar la libreria ISLR que contiene el dataset Smarket que se refiere al retorno diario del S&P500 entre 2001 y 2005. No olvides instalar la libreria con install.packages('ISLR')
require(ISLR)
require(tidyverse)#Explorando los datos
head(Smarket)La variable que nos interesa predecir es Direction que tiene valores Up cuando el retorno del dia es positivo y Down cuando el retorno es negativo.
1.1.1 Usado Regresión Lineal
Vamos a utilizar la herramienta que hemos visto hasta ahora y haremos una regresión lineal simple entre Direction ~ Lag1. Pero antes, tenemos que codificar la variable Direction. Para facilitar el manejo, vamos a codificar Up = 1 y Down = 0 (La variable Lag1 se refiere al rendimiento del día anterior)
#Codificación de la variable Direction
Smarket$Direction_coded = ifelse(Smarket$Direction == 'Up', 0, 1)
#Visualizacion de los datos
ggplot(Smarket, aes(x = Lag1, y = Direction_coded)) +
geom_point() +
geom_smooth(method = 'lm') +
labs(title = 'Regresion Lineal')## `geom_smooth()` using formula 'y ~ x'
#Análisis de Regresión
m.lineal = lm(Direction_coded ~ Lag1, data = Smarket)
summary(m.lineal)##
## Call:
## lm(formula = Direction_coded ~ Lag1, data = Smarket)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5761 -0.4814 -0.4469 0.5161 0.6046
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.48153 0.01413 34.072 <2e-16 ***
## Lag1 0.01749 0.01244 1.406 0.16
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4997 on 1248 degrees of freedom
## Multiple R-squared: 0.001581, Adjusted R-squared: 0.0007806
## F-statistic: 1.976 on 1 and 1248 DF, p-value: 0.1601Después de revisar el resumen, vemos que la variable Lag1 es estadíticamente no significativa y que la \(R_{adj}^2 = 0.00\) lo que demuestra que este modelo quizás no es aecuado. Ahora vamos a revisar los residuales:
par(mfrow = c(2,2))
residuals = m.lineal$residuals
#Normalidad
plot(density(residuals), col = 'red', main = 'Desityplot de residuales') # density plot for 'speed'
polygon(density(residuals), col="red")
qqnorm(residuals)
qqline(residuals, col = 'red')
#Varianza Constante
plot(residuals, main = 'Constant Variance', type = 'l')
#Residuales independientes
acf(residuals)
Una imagen dice más que mil palabras!!! el modelo claramente viola los supuestos de normalidad y aunque la media del error es \(0\), la distribución parece ser bimodal. Esto es porque el modelo no está ayudando a clasificar correctamente si el mercado sube o baja. La mejor estimación es quedarse a la mitad del camino y decir: no se que va a suceder.
1.1.2 La Distribución Binomial
En estadística, la distribución binomial una distribución de probabilidad discreta que cuenta el número de éxitos en una secuencia \(n\) nensayos de Bernoulli independientes entre sí con una probabilidad fija \(p\) de ocurrencia de éxito entre los ensayos.
Un experimento de Bernoulli se caracteriza por ser dicotómico, esto es, solo dos resultados son posibles, a uno de estos se le denomina éxito y tiene una probabilidad de ocurrencia \(p\) y al otro se le denomina fracaso y tiene una probabilidad \(q = 1-p\).
En general, una variable aleatoria discreta \(X\) tiene una distribución binomial con parámetros \(n\) y \(0\leq p \leq 1\) y escribimos \(x\backsim\text{Binom}(n, p)\) con función de probabilidades:
\(f(x) = {n\choose x}p^x(1-p)^{n-x}\)
con la función rbino(n, x, p) podemos simular algunas distribuciones binomiales asumiendo que en 10 intentos queremos el 30%, 50% y 95% de éxito del evento de interés.
#1,000 numeros aleatorios con distribucion binomial con exito de 30%
x1 = rbinom(1000, 10, 0.3)
#1,000 numeros aleatorios con distribucion binomial con exito de 50%
x2 = rbinom(1000, 10, 0.5)
#1,000 numeros aleatorios con distribucion binomial con exito de 95%
x3 = rbinom(1000, 10, 0.95)
#hacer 3 histigramas en 1 panel
par(mfrow = c(1, 3))
hist(x1)
hist(x2)
hist(x3)
1.1.3 Modelo Logistico
Una forma de modelar esto es determinar \(P(Direccion = Up|Lag1)\) (la probabilidad de que la dirección sea Updado Lag1). De una forma mas general, podemos escribirlo como \(P(Y = 1 | X)\)
Pariendo de la base de un modelo lineal:
\(log(\frac{p(x)}{1-p(x)}) = \beta_0 + \beta_1x\)
A la fracción \(\frac{p(x)}{1-p(x)}\) se le conoce como razón de momios o en inglés odds ratio.
Si elevamos todo al número \(e\):
\(\frac{p(x)}{1-p(x)} = e^{\beta_0 + \beta_1x}\)
Despejando \(p(x)\):
\(p(x) = \frac{e^{\beta_0 + \beta_1x}}{1+e^{\beta_0 + \beta_1x}}\)
1.1.4 Estimación por Máxima Verosimilitud
En la regresión lineal vimos que podemos estimar los parámetros utilizando el método de mínimos cuadrados. sin embargo, en la regresión logística, no es posible llegar a una forma cerrada como lo hicimos en el caso de la regresión lineal.
La función de verosimilitud en el caso de la regresión logística es:
\(\ell(\beta_0,\beta_1) = \prod_{i=1}^np(x_i)^{y_{i}}(1-p(x_i))^{1-y_i}\)
Tomando logaritmo natural:
\(log(\ell(\beta_0,\beta_1)) = \sum_{i=1}^n[y_i(\beta_0+\beta_1x_i) - log(1+e^{\beta_0 + \beta_1x_i})]\)
Esta es la función que tenemos que optimizar para encontrar los parámetros \(\beta_0,\beta_1\)
Solo como demostración vamos a optimizar la función de maxima verosimilitud from the cratch!!
Vamos a usar el dataset Default en la libreria ISLR que contiene créditos vencidos:
head(Default)Vamos a hacer la regresión default ~ balance. Este proceso requiere 2 pasos:
Computar la máxima verosimilitud dados
x,yy un vector de parámetrosOptimizar lal función. Vamos a usar la librería
optimxque pidepar = vector de valores iniciales,fn = funcion a optimizar. Recordemos que queremos encontrar \(\beta_0, \beta_1\); es decir, 2 parámetros!
x2 = model.matrix(default ~ balance, data = Default)
y2 = ifelse(Default$default == 'No', 0, 1)
#beta_op = as.numeric(c(-3.09, -8.353e-06)) #Como referencia este es el valor optimo
logVerosimilitud = function(beta){
#beta = vector de parámetros
#x = matriz variable explicativa con la primer columna con 1's
#y = vector de respuestas entre 0 y 1
x = model.matrix(default ~ balance, data = Default)
y = ifelse(Default$default == 'No', 0, 1)
linearTransf = x %*% beta
parcial1 = y * linearTransf
parcial2 = log(1 + exp(linearTransf))
likelihood = -sum(parcial1 - parcial2)
return (likelihood)
}
require(optimx)## Loading required package: optimxbeta_min = optimx(par = c(0, 0), fn = logVerosimilitud)
beta_min$p1[1]## [1] -10.65206beta_min$p2[1]## [1] 0.005499188En R no es necesario hacer optimizaciones para estimar los valores de la regresión. Podemos usar la función glm(formula, data, family = 'binomial') para estimar nuestra regresión logística:
m.logistica = glm(default ~ balance, data = Default, family = 'binomial')
summary(m.logistica)##
## Call:
## glm(formula = default ~ balance, family = "binomial", data = Default)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2697 -0.1465 -0.0589 -0.0221 3.7589
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.065e+01 3.612e-01 -29.49 <2e-16 ***
## balance 5.499e-03 2.204e-04 24.95 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2920.6 on 9999 degrees of freedom
## Residual deviance: 1596.5 on 9998 degrees of freedom
## AIC: 1600.5
##
## Number of Fisher Scoring iterations: 8Gráficamente, la solución se ve así:
Default$y_code = ifelse(Default$default == 'No', 0, 1)
ggplot(Default, aes(x = balance, y = y_code)) +
geom_smooth(method = 'glm', method.args = list(family = "binomial"))## `geom_smooth()` using formula 'y ~ x'
1.1.5 Interpretación
La interpretación de los coeficientes de regresión en el caso de la Regresión Logística no es tan directa. Para poder obtener algo de significado, tenemos que calcular la razón de momios u odds ratio:
\(\text{odds} = \frac{p(y|x)}{1-p(y|x)} = e^{\beta_i}\)
En donde \(e\approx2.718\)
De esta forma, para la variable balance el odd ratio es:
odd_beta1 = exp(coef(m.logistica)[2])
print(paste("El odd ratio es: ", odd_beta1, sep=' '))## [1] "El odd ratio es: 1.00551406372555"La interpretación práctica es que por cada dolar de balance en deuda la probabilidad de pasar a cartera vencida se incrementa 1.005. Entonces, si tenemos 1,000 dolares de deuda, la probabilidad de pasar a cartera vencida es:
riesgo = exp(coef(m.logistica)[2] * 1000)
print(paste("Si tenemos $1000 usd en deuda, la probabilidad de pasar a cartera vencida se incrementa: ", riesgo, sep = ' '))## [1] "Si tenemos $1000 usd en deuda, la probabilidad de pasar a cartera vencida se incrementa: 244.42705746094"En este caso, no existe un valor de \(R^2\) como en la regresión lineal, pero sí podemos computar un equivalente: McFadden´s Pseudo R-squared
\(^2_{Mc} = 1 - \frac{ln \mathscr{L_c}}{ln \mathscr{L}_{null}}\)
En donde \(ln \mathscr{L_c}\) es el velor de la función de log-verosimilitud del modelo y \(ln \mathscr{L}_{null}\) es el valor de la función de log-verosimilitud del modelo nulo \(y_i = \beta_0\) (un modelo que sólo tiene la constante o intercepto):
#Computo del modelo nulo
m.null = glm(default ~ 1, data = Default, family = 'binomial')
#Computo de la pseudo r2
r2 = 1 -(logLik(m.logistica) / logLik(m.null))
print(paste("La pseudo-R2 es:", r2, sep = ' '))## [1] "La pseudo-R2 es: 0.453391593901766"Esto quiere decir que el modelo default ~ balance explica el 45% de la variación.
1.1.6 Predicción de Diabetes
Ahora vamos a intentar predecir casos de Diabetes en mujeres en el dataset PimaIndiansDiabetes2 dentro de la libreria mlbench. Para ayudar a visualizar el modelo vamos a requerir las librerias visreg, margins, rcompanion
library(mlbench) # for PimaIndiansDiabetes2 dataset
library(visreg) # for potting logodds and probability
library(margins) # to calculate Average Marginal Effects
library(rcompanion) # to calculate pseudo R2
#Loa Data
data(PimaIndiansDiabetes2)
# See the data strcuture
glimpse(PimaIndiansDiabetes2, width = 50)## Rows: 768
## Columns: 9
## $ pregnant <dbl> 6, 1, 8, 1, 0, 5, 3, 10, 2, 8, ~
## $ glucose <dbl> 148, 85, 183, 89, 137, 116, 78,~
## $ pressure <dbl> 72, 66, 64, 66, 40, 74, 50, NA,~
## $ triceps <dbl> 35, 29, NA, 23, 35, NA, 32, NA,~
## $ insulin <dbl> NA, NA, NA, 94, 168, NA, 88, NA~
## $ mass <dbl> 33.6, 26.6, 23.3, 28.1, 43.1, 2~
## $ pedigree <dbl> 0.627, 0.351, 0.672, 0.167, 2.2~
## $ age <dbl> 50, 31, 32, 21, 33, 30, 26, 29,~
## $ diabetes <fct> pos, neg, pos, neg, pos, neg, p~Haciendo una inspección rápida con glimpse podemos ver que tenermos 768 casos o renglones y 9 columnas. La columna diabetes es nuestra variable de respuesta. Las variables explicativas son pregnant, glucose, pressure, triceps, insulin, pedigree, age. También podemos ver que hay iformación faltante o con NA.
Em R podemos eliminar los renglones que al menos tienen 1 valor de NA o faltante (missing data) con la función na.omit(). Ya que los datos no tienen NA´s, vamos a correr la regresión con glm():
#Eliminar los datos faltantes
data_full = na.omit(PimaIndiansDiabetes2)
#Correr el modelo
m.diabetes = glm(diabetes ~., data = data_full, family = 'binomial')
summary(m.diabetes)##
## Call:
## glm(formula = diabetes ~ ., family = "binomial", data = data_full)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7823 -0.6603 -0.3642 0.6409 2.5612
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.004e+01 1.218e+00 -8.246 < 2e-16 ***
## pregnant 8.216e-02 5.543e-02 1.482 0.13825
## glucose 3.827e-02 5.768e-03 6.635 3.24e-11 ***
## pressure -1.420e-03 1.183e-02 -0.120 0.90446
## triceps 1.122e-02 1.708e-02 0.657 0.51128
## insulin -8.253e-04 1.306e-03 -0.632 0.52757
## mass 7.054e-02 2.734e-02 2.580 0.00989 **
## pedigree 1.141e+00 4.274e-01 2.669 0.00760 **
## age 3.395e-02 1.838e-02 1.847 0.06474 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.10 on 391 degrees of freedom
## Residual deviance: 344.02 on 383 degrees of freedom
## AIC: 362.02
##
## Number of Fisher Scoring iterations: 5Después de correr el modelo notamos que hay algunas variables que nos significativas por ejemplo preassure cuyo Pr(>z) = 0.904 que es mayor de 0.05. Por lo tanto, la variable no es estadíticamente significativa.
Antes de continuar, ahora es un buen momento para revisar el conceto de stepwise selection
1.1.7 Feature Selection (seleccion de variables)
La selección de vriables o Feature Selection se refiere al procedimiento para descartar variables no significtivas en los modelos. Lo que buscamos es Parsimonia es decir, tener un modelo lo más simple posible y que explique la mayor cantidad de variación.
Existen varios algoritmos para lograr parsimonia pero ahora vamos a revisar el algoritmo de Stepwise. El objetivo de este algoritmo es introducir o eliminar variables de tal forma que la menor cantidad de variables presenten un impacto significativo para el problema que estamos modelando.
Podemos correr el algoritmo con la función step(modelo)
m.diabetes2 = step(m.diabetes)## Start: AIC=362.02
## diabetes ~ pregnant + glucose + pressure + triceps + insulin +
## mass + pedigree + age
##
## Df Deviance AIC
## - pressure 1 344.04 360.04
## - insulin 1 344.42 360.42
## - triceps 1 344.45 360.45
## <none> 344.02 362.02
## - pregnant 1 346.24 362.24
## - age 1 347.55 363.55
## - mass 1 350.89 366.89
## - pedigree 1 351.58 367.58
## - glucose 1 396.95 412.95
##
## Step: AIC=360.04
## diabetes ~ pregnant + glucose + triceps + insulin + mass + pedigree +
## age
##
## Df Deviance AIC
## - insulin 1 344.42 358.42
## - triceps 1 344.46 358.46
## <none> 344.04 360.04
## - pregnant 1 346.24 360.24
## - age 1 347.60 361.60
## - mass 1 351.28 365.28
## - pedigree 1 351.67 365.67
## - glucose 1 397.31 411.31
##
## Step: AIC=358.42
## diabetes ~ pregnant + glucose + triceps + mass + pedigree + age
##
## Df Deviance AIC
## - triceps 1 344.89 356.89
## <none> 344.42 358.42
## - pregnant 1 346.74 358.74
## - age 1 347.87 359.87
## - mass 1 351.32 363.32
## - pedigree 1 351.90 363.90
## - glucose 1 411.11 423.11
##
## Step: AIC=356.89
## diabetes ~ pregnant + glucose + mass + pedigree + age
##
## Df Deviance AIC
## <none> 344.89 356.89
## - pregnant 1 347.23 357.23
## - age 1 348.72 358.72
## - pedigree 1 352.72 362.72
## - mass 1 360.44 370.44
## - glucose 1 411.85 421.85summary(m.diabetes2)##
## Call:
## glm(formula = diabetes ~ pregnant + glucose + mass + pedigree +
## age, family = "binomial", data = data_full)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8827 -0.6535 -0.3694 0.6521 2.5814
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.992080 1.086866 -9.193 < 2e-16 ***
## pregnant 0.083953 0.055031 1.526 0.127117
## glucose 0.036458 0.004978 7.324 2.41e-13 ***
## mass 0.078139 0.020605 3.792 0.000149 ***
## pedigree 1.150913 0.424242 2.713 0.006670 **
## age 0.034360 0.017810 1.929 0.053692 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.10 on 391 degrees of freedom
## Residual deviance: 344.89 on 386 degrees of freedom
## AIC: 356.89
##
## Number of Fisher Scoring iterations: 5Ahora tenemos un modelo con 5 variables vs 8 variables del modelo anterior. Si calculamos la \(R^2\) tenemos:
#Computo del modelo nulo
m.null = glm(diabetes ~ 1, data = data_full, family = 'binomial')
#Computo de la pseudo r2
r2_m1 = 1 -(logLik(m.diabetes) / logLik(m.null))
r2_m2 = 1 -(logLik(m.diabetes2) / logLik(m.null))
print(paste("La pseudo-R2 del modelo de 8 var es:", r2_m1, sep = ' '))## [1] "La pseudo-R2 del modelo de 8 var es: 0.309329958576541"print(paste("La pseudo-R2 del modelo de 5 var es:", r2_m2, sep = ' '))## [1] "La pseudo-R2 del modelo de 5 var es: 0.307595560186559"Y lo que observamos es que nuestro 2do modelo es mas simple (tiene menos variables) y tiene prácticamente el mismo desempeño.
1.2 Caso 2: \(y\) es un conteo. Regresión de Poisson
1.2.1 La distribución de Poisson
En teoría de probabilidad y estadística, la distribución de Poisson es una distribución de probabilidad discreta que expresa, a partir de una frecuencia de ocurrencia media, la probabilidad de que ocurra un determinado número de eventos durante cierto período de tiempo.
Fue propuesta por Siméon-Denis Poisson, que la dio a conocer en 1838 en su trabajo Recherches sur la probabilité des jugements en matières criminelles et matière civile (Investigación sobre la probabilidad de los juicios en materias criminales y civiles).
Sea \(\lambda > 0\) un parámetro y \(X\) una variable aleatoria discreta, se dice que la variable aleatoria \(X\) tiene una distribución de Poisson con parámetro \(\lambda\) y escribimos \(X \sim \text{Poisson}(\lambda)\) con función de probabilidades:
\(p(k) = P[X=k] = \frac{e^{-\lambda}\lambda^k}{k!}\)
en donde \(k=0,1,2,..,n\) es el número de ocurrencias de un evento. Esta distribución tiene la eculariedad de:
\(\mu = \lambda\)
\(\sigma^2 = \lambda\)
https://es.wikipedia.org/wiki/Distribuci%C3%B3n_de_Poisson
Supongamos que una enfermedad se propaga una tasa de 2.3 casos por dia. Si asumimos que el número de casos por dia sigue una distribución de Poisson entonces, la distribución de probabilidad es:
lambda = 2.3
x = rpois(1000, lambda)
hist(x)
Qué pasa si asumimos que la distribución de casos por día siguiera una distibución normal?
#simulamos 1000 datos con distribución normal
normal = rnorm(1000, lambda, lambda)
#Ponemos todo en un dataframe
hist(normal)
abline(v = 0, col = 'red', lty = 2)
Lo que vemos, es que asumiendo de manera equivocada este comportamiento, tendríamos una probabilidad de generar eventos negativos, cuando sabemos que no podemos tener casos diarios negativos. En este caso la distribución de Poisson es un supuesto más adecuado.
1.2.2 Estimación de la Regresión de Poisson
Si se cumple que \(y_i \sim \text{Poisson}(\lambda)\) entonces:
\(\hat{y_i} = \text{exp}(\sum_{j=1}^p \beta_jx_{ij})\)
Es común hacer \(log(y_i)\), entonces:
\(\text{log}(y_i) = \sum_{j=1}^p \beta_jx_{ij} + \epsilon_i\)
Con el cual podemos usar mínimos cuadrados para resolver y encontrar los valores de \(\beta_j\).
Vamos a retomar los datos del precio de diamantes de diamonds y vamos a generar una regresión de poisson de price~carat:
m.poisson = glm(price ~ carat, data = diamonds, family = 'poisson')
summary(m.poisson)##
## Call:
## glm(formula = price ~ carat, family = "poisson", data = diamonds)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1200.42 -22.44 -8.66 8.53 162.42
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.965e+00 1.379e-04 50513 <2e-16 ***
## carat 1.324e+00 9.636e-05 13745 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 183127012 on 53939 degrees of freedom
## Residual deviance: 42213813 on 53938 degrees of freedom
## AIC: 42732976
##
## Number of Fisher Scoring iterations: 5Computamos la medida de desempeño \(R^2\)
#El modelo nulo es:
m.null = glm(price ~ 1, data = diamonds, family = 'poisson')
#Computo de R2
LogLik.model = logLik(m.poisson)
LogLik.null = logLik(m.null)
r2_m1 = 1 -(LogLik.model / LogLik.null)
r2_m1## 'log Lik.' 0.7673081 (df=2)Usando la libreria visreg para visualizar:
visreg(m.poisson, 'carat', scale = 'response', rug = 2)