El modelo

Sea $p$ la probabilidad de que $y_i=1$, es decir, de que la persona i-ésima (seleccionada al azar) presente rezago escolar:

$\mathbb{P}(y_i=1)=p,$ y, naturalmente:

$\mathbb{P}(y_i=0)=1-p$

En este caso, es fácil deducir que la función de probabilidad para $y_i$ es (Bernoulli):

$f(y_i) =p^{y_i}(1-p)^{1-y_i}; \ y_i =0,1,$ y que: $\begin{aligned} \mathbb{E}(y_i)&=p \\ \mathbb{V}(y_i)&=1-p \end{aligned}$

El objetivo del modelo es estimar la probabilidad de que un individuo presente rezago escolar. De manera general, el modelo que utilizaremos tiene la forma:

$y_i = F(z_i) + u_i,$ tal que la probabilidad de que un individuo presente rezago puede estimarse como:

$\begin{aligned} \widehat{y}_i &= \mathbb{E}(y_i|x_i)\\ &=\widehat{p}_i \\ &=F(\widehat{z}_i), \end{aligned}$ donde $z_i$ no es más que una combinación lineal de las variables explicativas:

$\begin{aligned} z_i &= \mathbf{x}\boldsymbol{\beta}\\ &=\beta_1 + \beta_2 x_{i2} + \dots + \beta_k x_{ik} \end{aligned}$

El problema al que nos enfrentamos en este punto es que no sabemos quién es la función $F(\cdot)$. Aquí veremos tres posibilidades.

El modelo de probabilidad lineal

La opción más simple para modelar la función $F(\cdot)$ es suponer que esta es lineal. Es decir:

$\begin{aligned} F(z_i)&=z_i\\ &= \mathbf{x}\boldsymbol{\beta}\\ &= \beta_1 + \beta_2 x_{i2} + \dots + \beta_k x_{ik} \end{aligned}$ En este caso, el modelo a estimar se vuelve:

$y_i = \beta_1 + \beta_2 x_{i2} + \dots + \beta_k x_{ik} + u_i$ Esta especificación es idéntica a la del modelo de regresión lineal. Por lo tanto, si se cumplen los supuestos de Gauss-Markov sobre el error, los coeficientes podrían estimarse utilizando el procedimiento de mínimos cuadrados ordinarios. En R, podemos estimar un modelo por el método de los mínimos cuadrados ordinarios ya sea con la función lm o con la función glm.

m.ols <- lm(d_asistencia_cr~anemia+d_mujer+d_indigena+d_vpadres+edad_d, data=db)
summary(m.ols)

## 
## Call:
## lm(formula = d_asistencia_cr ~ anemia + d_mujer + d_indigena + 
##     d_vpadres + edad_d, data = db)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5502 -0.3441 -0.2819  0.6080  0.7900 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.557335   0.031373  17.765  < 2e-16 ***
## anemia       0.063650   0.018155   3.506 0.000457 ***
## d_mujer     -0.048067   0.008757  -5.489 4.14e-08 ***
## d_indigena   0.094351   0.010484   9.000  < 2e-16 ***
## d_vpadres    0.023807   0.009144   2.603 0.009243 ** 
## edad_d      -0.015749   0.001931  -8.156 3.83e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4672 on 11495 degrees of freedom
##   (21088 observations deleted due to missingness)
## Multiple R-squared:  0.01774,    Adjusted R-squared:  0.01731 
## F-statistic: 41.52 on 5 and 11495 DF,  p-value: < 2.2e-16

Los resultados se intepretan de la manera usual. Por ejemplo, la probabilidad de que un individuo presente rezago se incrementa bajo la presencia de anemia. En particular, la probabilidad de presentar rezago es, en promedio, 6.8% mayor para los individuos con anemia que para los individuos sanos. Tenemos una interpretación similar para las variables d_indigena y d_vpadres. En el caso del género, las mujeres tienen un 4.89% menos de probabilidades de presentar rezago que los hombres. Finalmente, la probabilidad de rezagarse disminuye con la edad a una tasa de 1.57% por cada año. Todas las variables del modelo son estadísticamente significativas.

Existen al menos dos problemas asociados a la estimación de los coeficientes por mínimos cuadrados ordinarios:

1. Los errores no son continuos, ni normales, ni homoscedásticos. De hecho, es fácil demostrar que la varianza de los errores depende de las variables explicativas del modelo: $\mathbb{V}(u_i) = z_i(1-z_i)$.

2. El efecto marginal (o incremental, en el caso de las dummies) de $x_k$ es constante e igual a $\beta_k$. Matemáticamente, esto permite que el valor de la probabilidad estimada, $\widehat{p}_i$, pueda exceder 1 o ser inferior que $0$. Esto es, por supuesto, totalmente contraintuitivo.

Para el primer problema, tenemos una solución “simple”: estimar los coeficientes por el método de los mínimos cuadrados generalizados factibles o bien por el método de máxima verosimilitud. Sin embargo, para el segundo problema la solución es más elaborada. Esto implica buscar una especificación para $F(z_i)$ que nos permita acotar la probabilidad estimada a valores entre 0 y 1. Esto nos conduce al siguiente par de modelos.

El modelo PROBIT

Recuerda que, si $z$ es una variable aleatoria normal estándar, su función de densidad de probabilidad se calcula como:

$\phi_z = \frac{1}{\sqrt{2\pi}}\exp\left\{-\frac{1}{2}z^2\right\}$ Y su función de densidad acumulada es:

$\Phi(z)=\mathbb{P}(Z\leq z) = \int_{-\infty}^{z} \frac{1}{\sqrt{2\pi}}\exp\left\{-\frac{1}{2}u^2 du\right\}$

Al ser esta una función de densidad acumulada, los valores de $\mathbb{P}(Z\leq z)$ están acotados entre 0 y 1. En el modelo PROBIT suponemos que:

$F(z_i) = \Phi(z_i)$ Naturalmente, dado que el modelo resultante deja de ser lineal, sus coeficientes ya no puede estimarse por el método de los mínimos cuadrados ordinarios. En cambio, el modelo se estima por máxima verosimilitud. Además, es importante destacar que el modelo puede interpretarse como un modelo lineal generalizado en el que la función de conexión una función normal estándar y la distribución del error es binomial.

m.probit <- glm(d_asistencia_cr~anemia+d_mujer+d_indigena+d_vpadres+edad_d, family=binomial(link="probit"), 
                   data=db)
summary(m.probit)

## 
## Call:
## glm(formula = d_asistencia_cr ~ anemia + d_mujer + d_indigena + 
##     d_vpadres + edad_d, family = binomial(link = "probit"), data = db)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2850  -0.9174  -0.8139   1.3722   1.7466  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.18978    0.08745   2.170 0.030003 *  
## anemia       0.17116    0.04978   3.439 0.000585 ***
## d_mujer     -0.13396    0.02446  -5.478 4.31e-08 ***
## d_indigena   0.25399    0.02872   8.843  < 2e-16 ***
## d_vpadres    0.06938    0.02564   2.706 0.006804 ** 
## edad_d      -0.04401    0.00541  -8.135 4.10e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 14635  on 11500  degrees of freedom
## Residual deviance: 14431  on 11495  degrees of freedom
##   (21088 observations deleted due to missingness)
## AIC: 14443
## 
## Number of Fisher Scoring iterations: 4

En este punto, es importante destacar que los valores de los coeficientes estimados nos permiten solo recuperar el valor estimado de $z_i$, $\widehat{z}_i$. Por lo tanto, estos no tienen una interpretación directa. Para poder interpretar los resultados aun tenemos que evaluar los valores de $\widehat{z}_i$ en la función $F(\widehat{z}_i) = \Phi(\widehat{z}_i)=\widehat{p}$. Para hacer esto en R, podemos utilizar la función predict(). Esta implica crear una base de datos con los valores que deseamos utilizar para pronosticar utilizando el modelo PROBIT. Por ejemplo, creamos un data frame llamado means que contiene la edad promedio de nuestra muestra y todas las variables dummy fijadas en 1. Nota que conservamos los nombres y la estructura de las variables en la base de datos original:

means <- data.frame(anemia=1, d_mujer=1, d_indigena=1, d_vpadres=1, edad_d = mean(db$edad_d))

Ahora usamos la función predict(). Esta tiene como primer argumento el modelo estimado. En nuestro ejemplo, el modelo estimado usando una función PROBIT ha sido guardado en un objeto llamado m.probit. El segundo argumento, newdata, corresponde a la base de datos con los valores que deseamos usar para pronosticar. El argumento type=response indica que el pronóstico se hace en función la escala de la variable dependiente; en nuestro caso esto corresponde a probabilidades. El último argumento indica que deseamos que R también calcule el error estándar del pronóstico.

predict(m.probit, newdata=means, type="response", se.fit=TRUE)

## $fit
##         1 
## 0.4469431 
## 
## $se.fit
##          1 
## 0.02133842 
## 
## $residual.scale
## [1] 1

Entonces, la probabilidad de rezagarse para un individuo que presenta anemia, es mujer, es indígena, vive con ambos padres y tiene actualmente 15.53 años (la edad promedio) es de 44.69%.

Otra complicación en la lectura de los resultados es que no podemos calcular de forma simple los efectos marginales de las variables explicativas sobre la variable dependiente. Para un modelo PROBIT, el efecto marginal de la variable $x_k$ sobre $p=1$ se calcula como:

$\frac{\partial p}{\partial x_k}=\frac{\partial \Phi(z)}{\partial x_k} = \phi(z)\beta_k$ Para realizar este cálculo, en R podemos utilizar la función probitmfx(). Los argumentos más importantes de esta función son la fórmula del modelo y la base de datos con la cual hacemos la estimación. En el fondo, la función probitmfx() vuelve a estimar el modelo utilizando una especificación PROBIT, pero ahora nos devuelve los efectos marginales.

probitmfx(formula=d_asistencia_cr~anemia+d_mujer+d_indigena+d_vpadres+edad_d, data = db)

## Call:
## probitmfx(formula = d_asistencia_cr ~ anemia + d_mujer + d_indigena + 
##     d_vpadres + edad_d, data = db)
## 
## Marginal Effects:
##                 dF/dx  Std. Err.       z     P>|z|    
## anemia      0.0638614  0.0190359  3.3548 0.0007942 ***
## d_mujer    -0.0485331  0.0088492 -5.4845 4.147e-08 ***
## d_indigena  0.0944929  0.0109073  8.6633 < 2.2e-16 ***
## d_vpadres   0.0250327  0.0092055  2.7193 0.0065417 ** 
## edad_d     -0.0159536  0.0019601 -8.1393 3.976e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
## [1] "anemia"     "d_mujer"    "d_indigena" "d_vpadres"

Así pues, cada año adicional de vida disminuye, en promedio, la probabilidad de rezagarse en un 1.59% y la probabilidad de rezagarse es, en promedio, 6.38% veces mayor para los individuos con anemia.

El modelo LOGIT

El modelo LOGIT tiene exactamente la misma lógica que el modelo PROBIT, excepto que ahora la función de densidad acumulada tiene la forma:

$\Lambda(z)=\mathbb{P}(Z\leq z) = p = \frac{1}{1+e^{-z}}$

Y, entonces:

$F(z_i) = \Lambda(z_i)$

Además, la función de densidad de probabilidad es:

$\frac{\partial p}{\partial z} = \lambda(z) = \frac{e^{-z}}{(1+e^{-z})^2}$ Al igual que en el caso del PROBIT, el modelo LOGIT se estima por máxima verosimilitud. El modelo también puede interpretarse como un modelo lineal generalizado en el que la función de conexión una función logística y la distribución del error es binomial.

m.logit <- glm(d_asistencia_cr~anemia+d_mujer+d_indigena+d_vpadres+edad_d, family=binomial(link="logit"), 
                   data=db)
summary(m.logit)

## 
## Call:
## glm(formula = d_asistencia_cr ~ anemia + d_mujer + d_indigena + 
##     d_vpadres + edad_d, family = binomial(link = "logit"), data = db)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2921  -0.9162  -0.8131   1.3732   1.7412  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.326036   0.143446   2.273 0.023033 *  
## anemia       0.281071   0.080684   3.484 0.000495 ***
## d_mujer     -0.220130   0.040164  -5.481 4.24e-08 ***
## d_indigena   0.415407   0.046547   8.925  < 2e-16 ***
## d_vpadres    0.110717   0.042201   2.624 0.008702 ** 
## edad_d      -0.072307   0.008905  -8.119 4.69e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 14635  on 11500  degrees of freedom
## Residual deviance: 14431  on 11495  degrees of freedom
##   (21088 observations deleted due to missingness)
## AIC: 14443
## 
## Number of Fisher Scoring iterations: 4

Para poder interpretar los resultados del modelo LOGIT, hacemos uso nuevamente de la función predict(). Esta tiene exactamente la misma estructura de antes, solo que ahora hacemos referencia al modelo almacenado en el objeto m.logit:

predict(m.logit, newdata=means, type="response", se.fit=TRUE)

## $fit
##         1 
## 0.4476621 
## 
## $se.fit
##          1 
## 0.02160377 
## 
## $residual.scale
## [1] 1

Nota que los resultados son muy parecidos que en el caso del modelo PROBIT: la probabilidad de rezagarse para un individuo que presenta anemia, es mujer, es indígena, vive con ambos padres y tiene actualmente 15.53 años (la edad promedio) es de 44.76%.

Finalmente, calculamos el efecto marginal de las variables explicativas sobre la probabilidad de presentar rezago escolar. Para el caso de un modelo LOGIT tenemos que:

$\frac{\partial p}{\partial x_k}=\frac{\partial \Phi(z)}{\partial x_k} = \lambda(z)\beta_k$

Para realizar este cálculo en R utilizar la función logitmfx(). Esta se utiliza e interpreta igual que la función probitmfx() que vimos en la sección anterior.

logitmfx(formula=d_asistencia_cr~anemia+d_mujer+d_indigena+d_vpadres+edad_d, data = db)

## Call:
## logitmfx(formula = d_asistencia_cr ~ anemia + d_mujer + d_indigena + 
##     d_vpadres + edad_d, data = db)
## 
## Marginal Effects:
##                 dF/dx  Std. Err.       z     P>|z|    
## anemia      0.0645098  0.0191137  3.3750  0.000738 ***
## d_mujer    -0.0486564  0.0088619 -5.4905 4.008e-08 ***
## d_indigena  0.0949932  0.0109239  8.6959 < 2.2e-16 ***
## d_vpadres   0.0243474  0.0092251  2.6393  0.008309 ** 
## edad_d     -0.0159895  0.0019663 -8.1317 4.233e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
## [1] "anemia"     "d_mujer"    "d_indigena" "d_vpadres"

Cada año adicional de vida disminuye, en promedio, la probabilidad de rezagarse en 1.59% y la probabilidad de rezagarse es, en promedio, 6.45% veces mayor para los individuos con anemia.

PROBIT vs LOGIT

Antes vimos que los dos modelos, PROBIT y LOGIT, arrojan resultados similares. Esto no es casualidad. De hecho, dado que la única diferencia entre ambos modelos es la función de densidad, esperamos que las probabilidades y los efectos marginales (o incrementales, en el caso de las dummies) sean virtualmente los mismos.

¿Entonces, cuál modelo preferimos? En la práctica no existe ninguna razón suficientemente importante para preferir uno u otro modelo, fuera de la tradición. Por ejemplo, en Economía suele prefirse el modelo PROBIT, mientras en el resto de las ciencias sociales y en las ciencias de la salud se acostumbra más el modelo LOGIT.

Quizás una de las razones que podrían inclinar la balanza en favor del modelo LOGIT es su fácil interpretación. Esto es porque los resultados pueden exprarse fácilmente en términos del cociente de probabilidades, también conocido como cociente de momios o riesgo relativo:

$\frac{p}{1-p}=\frac{\frac{1}{1+e^{-z}}}{1-\frac{1}{1+e^{-z}}}=e^z = \exp\{\mathbf{x}\boldsymbol{\beta}\}$ Por ejemplo, supongamos que en un estudio para verificar la eficacia de la vacuna contra el COVID-19, $y=1$ significa que el paciente sobrevivió y $y=0$ que el paciente murió, y entre las variables explicativas se encuentra alguna medida para la dosis tomada. Un cociente de probabilidades de 2 significa que la probabilidad de sobrevivir es 2 veces más alta que la probabilidad de morir.

En el caso del modelo de rezago escolar tenemos:

logit.odds <- cbind(Coeficiente=round(coef(m.logit),4), OR=round(exp(coef(m.logit)),4))
print(logit.odds)

##             Coeficiente     OR
## (Intercept)      0.3260 1.3855
## anemia           0.2811 1.3245
## d_mujer         -0.2201 0.8024
## d_indigena       0.4154 1.5150
## d_vpadres        0.1107 1.1171
## edad_d          -0.0723 0.9302

Entonces, por ejemplo, el cociente de probabilidades es 1.32 para la variable anemia. Esto quiere decir que, en promedio, la probabilidad de rezarse es 1.32 veces más grande para los individuos con anemia. De igual forma, la probabilidad de rezagarse es 1.51 veces más grande para los indígenas y 1.11 veces más grande para las personas que viven con ambos padres. En cambio, la probabilidad de rezagarse es 1.25 (1/0.80) veces más grandre para los hombres que para las mujeres. En el caso de la variable edad, tenemos un cociente de probabilidades de 0.93. Entonces, con cada incremento unitario en la variable edad_d esperamos ver una caída de alrededor del 7% (0.93-1) en el cociente de probabilidades de presentar rezago.

Bondad de ajuste

Dado que el modelo de probabilidad lineal se estima e interpreta de la misma forma que un modelo de regresión lineal, los primeros indicadores para la bondad de ajuste son la R cuadrada y la R cuadrada ajustada. En nuestro ejemplo tenemos una R cuadrada ajustada de 1.7%. Este pésimo ajuste es relativamente normal, puesto que nuestra especificación carece de muchas otras variables que podrían ayudar a explicar el rezago escolar.

summary(m.ols)$adj.r.squared

## [1] 0.01731254

En el caso de los modelos PROBIT y LOGIT no existe una medida de bondad de ajuste equivalente a la R cuadrada. Existen, sin embargo, varias otras medidas que nos permiten comparar entre diferentes especificaciones del modelo. Una de ellas es el número de casos correctamente predichos. Sea $\widehat{p}_i$ la probabilidad estimada a partir del modelo PROBIT o LOGIT. Si $\widehat{p}_i>0.5$ fijamos $\widehat{y}_i=1$; en caso contrario fijamos $\widehat{y}_i=0$. Un caso es correctamente predicho en tanto $\widehat{y}_i = y_i$. Este cálculo puede hacerse muy fácilmente en R mediante tablas de confusión.

Consideremos primero el modelo PROBIT. El primer paso es utilizar nuevamente la función predict(). A diferencia de nuestros ejemplos anteriores, ahora usamos como argumento la base de datos original para hacer la predicción. Esto es porque queremos calcular la probabilidad de presentar anemia para cada individio de la muestra dadas sus características específicas. Luego usamos la función confusionMatrix(). En su versión básica, esta función tiene como argumentos dos vectores: el vector para el cual calculamos la matriz de confusión, $\widehat{y}_i$, y un vector de referencia, $y_i$.

probit.probs <- predict(m.probit, newdata=db, type="response")

PROBIT.y <- as.factor(as.numeric(probit.probs>0.5))
PROBIT.ref <- as.factor(db$d_asistencia_cr)
probit.cmtx <- confusionMatrix(data = PROBIT.y, reference = PROBIT.ref)
print(probit.cmtx)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 7648 3798
##          1   24   31
##                                          
##                Accuracy : 0.6677         
##                  95% CI : (0.659, 0.6763)
##     No Information Rate : 0.6671         
##     P-Value [Acc > NIR] : 0.4493         
##                                          
##                   Kappa : 0.0066         
##                                          
##  Mcnemar's Test P-Value : <2e-16         
##                                          
##             Sensitivity : 0.996872       
##             Specificity : 0.008096       
##          Pos Pred Value : 0.668181       
##          Neg Pred Value : 0.563636       
##              Prevalence : 0.667072       
##          Detection Rate : 0.664986       
##    Detection Prevalence : 0.995218       
##       Balanced Accuracy : 0.502484       
##                                          
##        'Positive' Class : 0              
## 

Los resultados nos indican que, de los individuos en la muestra que no presentaban rezago escolar al momento de la encuesta, 7648 fueron predichos correctamente por el modelo y 24 fueron predichos incorrectamente. De la misma forma, de los individuos en la muestra con rezago escolar al momento de la encuesta, 31 fueron predichos correctamente y 3798 casos fueron predichos incorrectamente. Entonces, el porcentaje de efectividad del modelo es de 66% (7669 casos de los 11501 considerados). 66% es la medida de bondad de ajuste.

El análisis es fundamentalmente el mismo para un modelo LOGIT. De hecho, los resultados, como esperábamos, son virtualmente idénticos.

logit.probs <- predict(m.logit, newdata=db, type="response")

LOGIT.y <- as.factor(as.numeric(logit.probs>0.5))
LOGIT.ref <- as.factor(db$d_asistencia_cr)
logit.cmtx <- confusionMatrix(data = LOGIT.y, reference = LOGIT.ref)
print(logit.cmtx)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 7648 3798
##          1   24   31
##                                          
##                Accuracy : 0.6677         
##                  95% CI : (0.659, 0.6763)
##     No Information Rate : 0.6671         
##     P-Value [Acc > NIR] : 0.4493         
##                                          
##                   Kappa : 0.0066         
##                                          
##  Mcnemar's Test P-Value : <2e-16         
##                                          
##             Sensitivity : 0.996872       
##             Specificity : 0.008096       
##          Pos Pred Value : 0.668181       
##          Neg Pred Value : 0.563636       
##              Prevalence : 0.667072       
##          Detection Rate : 0.664986       
##    Detection Prevalence : 0.995218       
##       Balanced Accuracy : 0.502484       
##                                          
##        'Positive' Class : 0              
## 

Última actualización: 27-08-2020.