Modelo logit
Orlando Moscote Flórez
2022-04-06
1. El modelo logit.
El modelo logit se utiliza cuando se tiene una variable dependiente o de respuesta binaria, en la cual se supone que solo asume dos valores posibles \(0\) y \(1\), que son asignaciones arbitrarias a una respuesta cualitativa. Tal respuesta puede corresponder a no y sí, a incumplimiento y cumplimiento, en general, a fracaso y éxito.
Se puede considerar el modelo de la forma
\[\begin{equation*} y_{i}=X_{i}^{\prime}\beta+\varepsilon_{i} \end{equation*}\]
En donde:
\(X_{i}^{\prime}=[1, x_{i1}, x_{i2},x_{i3},\cdots, x_{ik}]\)
\(\beta^{\prime}=[\beta_{0},\beta_{1},\beta_{2},\cdots,\beta_{k}]\)
Y la variable respuesta \(Y\) solo toma los valores \(0\) y \(1\), o sea que esta es una variable aleatoria tipo Bernoulli, que tiene distribución de probabilidad:
| \(Y_{i}\) | Probabilidad |
|---|---|
| 0 | \(1-\pi\) |
| 1 | \(\pi\) |
Como se parte del supuesto de que \(E(\varepsilon_{i})=0\), entonces el valor esperado de la variable dependiente será:
\(E(y_{i})=0\cdot (1-\pi_{i})+1\cdot \pi_{i}\)
\(E(y_{i})=\pi_{i}\)
o sea:
\(E(y_{i})=E(X_{i}^{\prime }\beta)=\pi_{i}\)
Esto implica que la respuesta esperada es simplemente la probabilidad de que la variable dependiente tome el valor \(1\).
Una restricción importante del modelo es que se debe cumplir que \(0\leq E(y_{i})=\pi_{i} \leq 1\). Esta restricción puede causar problemas en términos de seleccionar la función de respuesta lineal según el modelo especificado; por lo tanto, la función de respuesta debería ser no lineal. Una función monótonamente creciente o decreciente, en forma de \(S\) o de \(S\) invertida es la que suele utilizarse. Esta función es la función logística o logit que tiene la forma
\(E(y)=\dfrac{e^{(X^{\prime}\beta)}}{1+e^{X^{\prime}\beta)}}\)
O en forma
\(\pi _{i}=\dfrac{1}{1+e^{(-X^{\prime}\beta)}}\)
Que, a su vez, es equivalente a:
\(1-\pi_{i}=\dfrac{1}{e^{(X^{\prime}\beta)}}\)
Con lo cual se tiene:
\(\dfrac{\pi_{i}}{1-\pi_{i}}=\dfrac{1+e^{(X^{\prime}\beta)}}{1+e^{(-X^{\prime}\beta)}}\)
A esta transformación se le conoce como transformación logit de la probabilidad \(\pi_{i}\) y la relación \(\dfrac{\pi_{i}}{1-\pi_{i}}\) es una razón de probabilidades o ventaja (odds ratio).
Si se toma logaritmo, se obtiene:
\(L_{i}=ln\,\left( \dfrac{\pi_{i}}{1-\pi_{i}}\right) =X^{\prime}\beta\)
O sea que el logaritmo de la razón de probabilidades es lineal en las variables y en los parámetros.
2. Estimación de parámetros del modelo logit. Método de máxima verosimilitud
Cada observación muestral tiene una distribución Bernoulli, cuya distribución tiene la forma:
\(f(y_{i})=\pi_{i}^{y_{i}}(1-\pi_{i})^{1-y_{i}} para i=1,2,3,\cdots,n\)
Donde cada observación \(y_{i}\) toma valor \(0\) o \(1\), y las observaciones son independientes, con lo cual la función de verosimilitud es:
\(L(y_{1}, y_{2}, y_{3},\cdots,y_{n},\beta)=\displaystyle\prod_{i=1}^{n} \pi_{i}^{y_{i}}(1-\pi_{i})^{1-y_{i}}\)
Tomando el logaritmo de la verosimilitud,
\(ln[ L(y_{1}, y_{2}, y_{3},\cdots,y_{n},\beta)]=Ln\left[ \displaystyle\prod_{i=1}^{n} \pi_{i}^{y_{i}}(1-\pi_{i})^{1-y_{i}}\right]\)
\(ln[ L(y_{1}, y_{2}, y_{3},\cdots,y_{n},\beta)]=\displaystyle\sum_{i=1}^{n}y_{i}\,Ln \left( \dfrac{\pi_{i}}{1-\pi_{i}} \right) +\displaystyle\sum_{i=1}^{n}y_{i}Ln\,(1-\pi_{i})\)
Como
\(1+\pi_{i}=[1+ e^{(X^{\prime}\beta)}]^{-1}\)
y,
\(ln\,\left( \dfrac{\pi_{i}}{1-\pi_{i}}\right) =X^{\prime}\beta\)
El logaritmo de la función de verosimilitud se puede expresar como:
\(ln(Y,\beta)=\displaystyle\sum_{i=1}^{n}y_{i}X^{\prime}\beta+\displaystyle\sum_{i=1}^{n}Ln[1+e^{(X^{\prime}\beta)}\)
3. Medida de bondad de ajuste.
La prueba de Hosmer-Lemeshow (1989), básicamente consiste en dividir el recorrido de la probabilidad en deciles, \(< 0.1\), \(< 0.2\) y así hasta \(<1\), y calcular tanto la distribución de éxitos como de fracasos, prevista por el modelo y los valores realmente observados. Ambas distribuciones, esperada y observada, se contrastan mediante una prueba chi-cuadrado.
Entonces, se tienen 10 grupos diferentes y se utiliza el estadístico:
\(G_{H-L}=\displaystyle\sum_{i=1}^{10}\dfrac{(O_{i}-E_{i})^{2}}{E_{i}\left( 1-\dfrac{E_{i}}{n_{i}}\right) }\)
Que tiene una distribución chi-cuadrado con 8 grados de libertad.
Donde:
\(n_{i}\) es el número de observaciones en cada grupo.
\(O_{i}\) es el número de casos en el grupo.
\(E_{i}\) es el número de casos esperados en cada grupo.
plot(logitp,p,t="l",xlab="log(p/1-p)",xaxt="n")
abline(h=c(0,0.5,1),v=0,col="grey")
4. Ejemplo 1.
Ejemplo tomado del libro Introducción al análisis de regresión lineal de Montgomery-Peck-Vining. Página 432.
Datos de 20 familias.
x: ingreso.
y: propiedad de la vivienda (1= sí, 0= no)
x<-c(38000,51200,39600,43400,47700,53000,41500,40800,45400,52400,38700,40100,
49500,38000,42000,54000,51700,39400,40900,52800)
y<-c(0,1,0,1,0,0,1,0,1,1,1,0,1,0,1,1,1,0,0,1)4.1. Gráfico de los datos.
plot(x,y,pch=16)
4.2. Modelo.
modelo1=glm(y~x,family=binomial(link="logit"))
summary(modelo1)##
## Call:
## glm(formula = y ~ x, family = binomial(link = "logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0232 -0.8766 0.5072 0.7980 1.6046
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.7395139 4.4394326 -1.969 0.0490 *
## x 0.0002009 0.0001006 1.998 0.0458 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.526 on 19 degrees of freedom
## Residual deviance: 22.435 on 18 degrees of freedom
## AIC: 26.435
##
## Number of Fisher Scoring iterations: 44.3. Coeficientes del modelo.
com<-coefficients(modelo1)
c1<-com[1];c1## (Intercept)
## -8.739514c2<-com[2];c2## x
## 0.0002009056exp(coefficients(modelo1))## (Intercept) x
## 0.0001601317 1.00020092584.4. Gráfico del modelo.
z<-function(x){
1/(1+exp(-(c1+c2*x)))
}
curve(z,xlim=c(30000,60000),lwd=2)
4.5. Predichos.
pred<-predict(modelo1)
pred## 1 2 3 4 5 6
## -1.10509968 1.54685473 -0.78365066 -0.02020924 0.84368500 1.90848488
## 7 8 9 10 11 12
## -0.40192995 -0.54256390 0.38160203 1.78794150 -0.96446574 -0.68319784
## 13 14 15 16 17 18
## 1.20531515 -1.10509968 -0.30147713 2.10939052 1.64730755 -0.82383179
## 19 20
## -0.52247333 1.86830375exp(pred)## 1 2 3 4 5 6 7 8
## 0.3311779 4.6966746 0.4567356 0.9799936 2.3249185 6.7428648 0.6690276 0.5812561
## 9 10 11 12 13 14 15 16
## 1.4646291 5.9771358 0.3811868 0.5049995 3.3378108 0.3311779 0.7397247 8.2432156
## 17 18 19 20
## 5.1929792 0.4387472 0.5930519 6.4773000pre1<-predict(modelo1,type="response");pre1## 1 2 3 4 5 6 7 8
## 0.2487856 0.8244590 0.3135336 0.4949479 0.6992408 0.8708488 0.4008487 0.3675914
## 9 10 11 12 13 14 15 16
## 0.5942594 0.8566747 0.2759850 0.3355480 0.7694690 0.2487856 0.4251964 0.8918125
## 17 18 19 20
## 0.8385268 0.3049509 0.3722741 0.8662619#probabilidad para un ingreso de 38000
odds<- predict(modelo1,data.frame(x=38000))
odds## 1
## -1.1051exp(odds)/(1+exp(odds))## 1
## 0.24878561/(1+exp(-(c1+c2*x)))## [1] 0.2487856 0.8244590 0.3135336 0.4949479 0.6992408 0.8708488 0.4008487
## [8] 0.3675914 0.5942594 0.8566747 0.2759850 0.3355480 0.7694690 0.2487856
## [15] 0.4251964 0.8918125 0.8385268 0.3049509 0.3722741 0.86626194.6. Intervalo de confianza para los parámetros.
confint(modelo1,level=0.95)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -1.89408e+01 -0.8444742832
## x 2.44334e-05 0.00043657214.7. Curva ROC.
roc(y,pre1,plot = TRUE, legacy.axes = TRUE,
percent = TRUE, xlab = "% Falsos positivos",
ylab = "% verdaderos postivios", col = "red", lwd = 2,
print.auc = TRUE)## Setting levels: control = 0, case = 1## Setting direction: controls < cases
##
## Call:
## roc.default(response = y, predictor = pre1, percent = TRUE, plot = TRUE, legacy.axes = TRUE, xlab = "% Falsos positivos", ylab = "% verdaderos postivios", col = "red", lwd = 2, print.auc = TRUE)
##
## Data: pre1 in 9 controls (y 0) < 11 cases (y 1).
## Area under the curve: 79.8%4.8. Matriz de confusión.
tabla1<-table(true=y,pred=round(fitted(modelo1)))
sum(diag(tabla1))/sum(tabla1)## [1] 0.7matriz1<-confusionMatrix(tabla1);matriz1## Confusion Matrix and Statistics
##
## pred
## true 0 1
## 0 7 2
## 1 4 7
##
## Accuracy : 0.7
## 95% CI : (0.4572, 0.8811)
## No Information Rate : 0.55
## P-Value [Acc > NIR] : 0.1299
##
## Kappa : 0.4059
##
## Mcnemar's Test P-Value : 0.6831
##
## Sensitivity : 0.6364
## Specificity : 0.7778
## Pos Pred Value : 0.7778
## Neg Pred Value : 0.6364
## Prevalence : 0.5500
## Detection Rate : 0.3500
## Detection Prevalence : 0.4500
## Balanced Accuracy : 0.7071
##
## 'Positive' Class : 0
## 4.9. Prueba de Hosmer-Lemeshow.
hl<-hoslem.test(modelo1$y,fitted(modelo1),g=10)
hl##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: modelo1$y, fitted(modelo1)
## X-squared = 11.809, df = 8, p-value = 0.1599cbind(x,y,pred,pre1)## x y pred pre1
## 1 38000 0 -1.10509968 0.2487856
## 2 51200 1 1.54685473 0.8244590
## 3 39600 0 -0.78365066 0.3135336
## 4 43400 1 -0.02020924 0.4949479
## 5 47700 0 0.84368500 0.6992408
## 6 53000 0 1.90848488 0.8708488
## 7 41500 1 -0.40192995 0.4008487
## 8 40800 0 -0.54256390 0.3675914
## 9 45400 1 0.38160203 0.5942594
## 10 52400 1 1.78794150 0.8566747
## 11 38700 1 -0.96446574 0.2759850
## 12 40100 0 -0.68319784 0.3355480
## 13 49500 1 1.20531515 0.7694690
## 14 38000 0 -1.10509968 0.2487856
## 15 42000 1 -0.30147713 0.4251964
## 16 54000 1 2.10939052 0.8918125
## 17 51700 1 1.64730755 0.8385268
## 18 39400 0 -0.82383179 0.3049509
## 19 40900 0 -0.52247333 0.3722741
## 20 52800 1 1.86830375 0.86626195. Ejemplo 2.
Ejemplo tomado del libro Introducción al análisis de regresión lineal de Montgomery-Peck-Vining. Página 434.
Datos de 20 familias.
x1: ingreso.
x2: edad (años) del coche mas viejo.
y: si seis meses después habian comprado otro coche (1= sí, 0= no)
x1<-c(45000,40000,60000,50000,55000,50000,35000,65000,53000,48000,37000,31000,
40000,75000,43000,49000,37500,71000,34000,27000)
x2<-c(2,4,3,2,2,5,7,2,2,1,5,7,4,2,9,2,4,1,5,6)
y2<-c(0,0,1,1,0,1,1,1,0,0,1,1,1,0,1,0,1,0,0,0)5.1. Modelo.
modelo2=glm(y2~x1+x2,family=binomial(link="logit"))
summary(modelo2)##
## Call:
## glm(formula = y2 ~ x1 + x2, family = binomial(link = "logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5635 -0.8045 -0.1397 0.9535 1.7915
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.047e+00 4.674e+00 -1.508 0.132
## x1 7.382e-05 6.371e-05 1.159 0.247
## x2 9.879e-01 5.274e-01 1.873 0.061 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.726 on 19 degrees of freedom
## Residual deviance: 21.082 on 17 degrees of freedom
## AIC: 27.082
##
## Number of Fisher Scoring iterations: 55.2. Predichos.
pred2<-predict(modelo2)
pred2## 1 2 3 4 5 6 7
## -1.7495337 -0.1428455 0.3456042 -1.3804498 -1.0113658 1.5832084 2.4517287
## 8 9 10 11 12 13 14
## -0.2731979 -1.1589994 -2.5159694 0.6235902 2.1564616 -0.1428455 0.4649700
## 15 16 17 18 19 20
## 5.0180352 -1.4542665 -0.3273875 -0.8181833 0.4021398 0.8733083exp(pred2)## 1 2 3 4 5 6
## 0.17385499 0.86688797 1.41284330 0.25146543 0.36372186 4.87055764
## 7 8 9 10 11 12
## 11.60839715 0.76094217 0.31380001 0.08078456 1.86561392 8.64050968
## 13 14 15 16 17 18
## 0.86688797 1.59196638 151.11409898 0.23357161 0.72080438 0.44123253
## 19 20
## 1.49502033 2.394820675.3. Coeficientes.
coma<-coefficients(modelo2); coma## (Intercept) x1 x2
## -7.047061e+00 7.381679e-05 9.878861e-01d1<-coma[1]
d2<-coma[2]
d3<-coma[3]5.4. Probabilidad para un ingreso de 45000 y 2 años.
odds2<- predict(modelo2,data.frame(x1=45000,x2=2))
odds2## 1
## -1.749534exp(odds2)/(1+exp(odds2))## 1
## 0.148106pre2<-predict(modelo2,type="response");pre2## 1 2 3 4 5 6 7
## 0.14810602 0.46434922 0.58555120 0.20093678 0.26671264 0.82965843 0.92068778
## 8 9 10 11 12 13 14
## 0.43212218 0.23884915 0.07474622 0.65103464 0.89627104 0.46434922 0.61419253
## 15 16 17 18 19 20
## 0.99342599 0.18934581 0.41887642 0.30614944 0.59920166 0.705433631/(1+exp(-(d1+d2*x1*d3*x2)))## [1] 0.38130110 0.99024866 0.99771477 0.56099382 0.72600005 0.99998610
## [7] 0.99997998 0.91929998 0.66434815 0.02801017 0.99841188 0.99984580
## [13] 0.99024866 0.97999030 1.00000000 0.52481740 0.97999030 0.13359283
## [19] 0.99527324 0.991560835.5. Intervalo de confianza.
confint(modelo2,level=0.95)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -1.805544e+01 1.0275430082
## x1 -4.361540e-05 0.0002184223
## x2 1.544228e-01 2.28721278555.6. Curva ROC.
roc(y2,pre2,plot = TRUE, legacy.axes = TRUE,
percent = TRUE, xlab = "% Falsos positivos",
ylab = "% verdaderos postivios", col = "red", lwd = 2,
print.auc = TRUE)## Setting levels: control = 0, case = 1## Setting direction: controls < cases
##
## Call:
## roc.default(response = y2, predictor = pre2, percent = TRUE, plot = TRUE, legacy.axes = TRUE, xlab = "% Falsos positivos", ylab = "% verdaderos postivios", col = "red", lwd = 2, print.auc = TRUE)
##
## Data: pre2 in 10 controls (y2 0) < 10 cases (y2 1).
## Area under the curve: 77.5%5.7. Matriz de confusión.
tabla2<-table(true=y2,pred=round(fitted(modelo2)))
sum(diag(tabla2))/sum(tabla2)## [1] 0.65matriz2<-confusionMatrix(tabla2);matriz2## Confusion Matrix and Statistics
##
## pred
## true 0 1
## 0 7 3
## 1 4 6
##
## Accuracy : 0.65
## 95% CI : (0.4078, 0.8461)
## No Information Rate : 0.55
## P-Value [Acc > NIR] : 0.252
##
## Kappa : 0.3
##
## Mcnemar's Test P-Value : 1.000
##
## Sensitivity : 0.6364
## Specificity : 0.6667
## Pos Pred Value : 0.7000
## Neg Pred Value : 0.6000
## Prevalence : 0.5500
## Detection Rate : 0.3500
## Detection Prevalence : 0.5000
## Balanced Accuracy : 0.6515
##
## 'Positive' Class : 0
## 5.8. Prueba de Hosmer-Lemeshow.
h2<-hoslem.test(modelo2$y,fitted(modelo2),g=10)
h2##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: modelo2$y, fitted(modelo2)
## X-squared = 7.3153, df = 8, p-value = 0.503