Logit, Probit y MLP con R

Logit

En este ocasión, se presenta una reproducción de los resultados de las páginas 605 y 613 del libro Basic Econometrics 4ta ed. de Gujarati que tratan sobre la estimación de modelos logits probit y mlp (modelo lineal de probabilidades).

Para la estimación de tanto el logit como el probit es necesario usar la función glm, sin embargo para la estimación del mlp se puede usar tanto lm como glm, reemplazando familia de probabilidad por la identidad, estas funciones dan resultados idénticos.

Los datos con los que se trabajará son los de la Tabla 15.7 que se importarán directamente de la web.

table_15.7 <- read.table("http://www.auburn.edu/~sabaric/forcasting/dataSets/Gujarati/Table%2015.7.txt", 
    skip = 16, header = TRUE)
head(table_15.7)  # un vistazo a las primeras observaciones de la tabla
##   OBS  GPA TUCE PSI GRADE LETTER
## 1   1 2.66   20   0     0      C
## 2   2 2.89   22   0     0      B
## 3   3 3.28   24   0     0      B
## 4   4 2.92   12   0     0      B
## 5   5 4.00   21   0     1      A
## 6   6 2.86   17   0     0      B

Logit

Una vez que se tienen los datos en R, vamos a estimar el logit que se presenta en la tabla 15.8 de la página 605 del libro.

logit <- glm(GRADE ~ GPA + TUCE + PSI, data = table_15.7, family = "binomial")
summary(logit)  # así de fácil.
## 
## Call:
## glm(formula = GRADE ~ GPA + TUCE + PSI, family = "binomial", 
##     data = table_15.7)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.955  -0.645  -0.257   0.589   2.097  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -13.0213     4.9313   -2.64   0.0083 **
## GPA           2.8261     1.2629    2.24   0.0252 * 
## TUCE          0.0952     0.1416    0.67   0.5014   
## PSI           2.3787     1.0646    2.23   0.0255 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 41.183  on 31  degrees of freedom
## Residual deviance: 25.779  on 28  degrees of freedom
## AIC: 33.78
## 
## Number of Fisher Scoring iterations: 5

La tabla 15.8 es muy fácil de construir si ya se tiene estimado el modelo:

# Creando la tabla 15.9 pagina 607
actual <- table_15.7[, "GRADE"]  # valores reales
fitted.prob <- plogis(predict(logit, type = "link"))  # valores estimados


table_15.9 <- data.frame(Actual = actual, Fitted = fitted.prob, Residual = actual - 
    fitted.prob)
table_15.9
##    Actual  Fitted Residual
## 1       0 0.02658 -0.02658
## 2       0 0.05950 -0.05950
## 3       0 0.18726 -0.18726
## 4       0 0.02590 -0.02590
## 5       1 0.56989  0.43011
## 6       0 0.03486 -0.03486
## 7       0 0.02650 -0.02650
## 8       0 0.05156 -0.05156
## 9       0 0.11113 -0.11113
## 10      1 0.69351  0.30649
## 11      0 0.02447 -0.02447
## 12      0 0.19000 -0.19000
## 13      0 0.32224 -0.32224
## 14      1 0.19321  0.80679
## 15      0 0.36099 -0.36099
## 16      0 0.03018 -0.03018
## 17      0 0.05363 -0.05363
## 18      0 0.03859 -0.03859
## 19      0 0.58987 -0.58987
## 20      1 0.66079  0.33921
## 21      0 0.06138 -0.06138
## 22      1 0.90485  0.09515
## 23      0 0.24177 -0.24177
## 24      0 0.85209 -0.85209
## 25      1 0.83829  0.16171
## 26      1 0.48113  0.51887
## 27      1 0.63542  0.36458
## 28      0 0.30722 -0.30722
## 29      1 0.84170  0.15830
## 30      1 0.94534  0.05466
## 31      0 0.52912 -0.52912
## 32      1 0.11103  0.88897

Para comparar mejor los resultados, podríamos crear una tabla que nos diga cuál es el valor real y el valor estimado en términos de los ceros y unos que son los que realmente se tienen en la tabla original de datos, así pues, se puede hacer los siguiente

data.frame(Actual = actual, Fitted = round(fitted.prob), OK_Wrong = ifelse(actual == 
    round(fitted.prob), "OK", "Wrong"))
##    Actual Fitted OK_Wrong
## 1       0      0       OK
## 2       0      0       OK
## 3       0      0       OK
## 4       0      0       OK
## 5       1      1       OK
## 6       0      0       OK
## 7       0      0       OK
## 8       0      0       OK
## 9       0      0       OK
## 10      1      1       OK
## 11      0      0       OK
## 12      0      0       OK
## 13      0      0       OK
## 14      1      0    Wrong
## 15      0      0       OK
## 16      0      0       OK
## 17      0      0       OK
## 18      0      0       OK
## 19      0      1    Wrong
## 20      1      1       OK
## 21      0      0       OK
## 22      1      1       OK
## 23      0      0       OK
## 24      0      1    Wrong
## 25      1      1       OK
## 26      1      0    Wrong
## 27      1      1       OK
## 28      0      0       OK
## 29      1      1       OK
## 30      1      1       OK
## 31      0      1    Wrong
## 32      1      0    Wrong

Como se aprecia en este último resultado, se tienen 6 resultados incorrectos que corresponden a los alumnos 14, 19, 24, 26, 31 y 32, estos son los mismos que se señalan en el libro.

Probit

Para la esiimación del probit en R se modifica ligeramente el comando anterior y sólo se cambia family = "binomial" a family = binomial(link = "probit") y se tiene el probit.

probit <- glm(GRADE ~ GPA + TUCE + PSI, data = table_15.7, family = binomial(link = "probit"))
summary(probit)
## 
## Call:
## glm(formula = GRADE ~ GPA + TUCE + PSI, family = binomial(link = "probit"), 
##     data = table_15.7)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.939  -0.651  -0.223   0.593   2.045  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -7.4523     2.5715   -2.90   0.0038 **
## GPA           1.6258     0.6897    2.36   0.0184 * 
## TUCE          0.0517     0.0812    0.64   0.5241   
## PSI           1.4263     0.5870    2.43   0.0151 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 41.183  on 31  degrees of freedom
## Residual deviance: 25.638  on 28  degrees of freedom
## AIC: 33.64
## 
## Number of Fisher Scoring iterations: 6

MLP

Por último, se estima el mlp. Como se mencionó anteriormente, hay dos formas de hacerlo: usando glm y la otra es usando lm.

# Estimando el linear probability model (mlp) con glm
mlp <- glm(GRADE ~ GPA + TUCE + PSI, data = table_15.7, family = gaussian(link = "identity"))
summary(mlp)
## 
## Call:
## glm(formula = GRADE ~ GPA + TUCE + PSI, family = gaussian(link = "identity"), 
##     data = table_15.7)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7815  -0.2773   0.0053   0.2109   0.8114  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -1.4980     0.5239   -2.86   0.0079 **
## GPA           0.4639     0.1620    2.86   0.0078 **
## TUCE          0.0105     0.0195    0.54   0.5944   
## PSI           0.3786     0.1392    2.72   0.0111 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1506)
## 
##     Null deviance: 7.2188  on 31  degrees of freedom
## Residual deviance: 4.2165  on 28  degrees of freedom
## AIC: 35.96
## 
## Number of Fisher Scoring iterations: 2
# Otra forma de estimar el mlp con lm
mlp2 <- lm(GRADE ~ GPA + TUCE + PSI, data = table_15.7)
summary(mlp2)  # es el mismo resultado que antes.

Por último, comparando los resultados.

library(memisc)
mtable(probit, logit, mlp)
## 
## Calls:
## probit: glm(formula = GRADE ~ GPA + TUCE + PSI, family = binomial(link = "probit"), 
##     data = table_15.7)
## logit: glm(formula = GRADE ~ GPA + TUCE + PSI, family = "binomial", 
##     data = table_15.7)
## mlp: glm(formula = GRADE ~ GPA + TUCE + PSI, family = gaussian(link = "identity"), 
##     data = table_15.7)
## 
## ===================================================
##                        probit     logit      mlp   
## ---------------------------------------------------
## (Intercept)           -7.452**  -13.021** -1.498** 
##                       (2.572)    (4.931)  (0.524)  
## GPA                    1.626*     2.826*   0.464** 
##                       (0.690)    (1.263)  (0.162)  
## TUCE                   0.052      0.095    0.010   
##                       (0.081)    (0.142)  (0.019)  
## PSI                    1.426*     2.379*   0.379*  
##                       (0.587)    (1.065)  (0.139)  
## ---------------------------------------------------
## Aldrich-Nelson R-sq.     0.327     0.325     0.086 
## McFadden R-sq.           0.377     0.374     0.416 
## Cox-Snell R-sq.          0.385     0.382     0.090 
## Nagelkerke R-sq.         0.532     0.528     0.443 
## phi                      1.000     1.000     0.151 
## Likelihood-ratio        15.546    15.404     3.002 
## p                        0.001     0.002     0.391 
## Log-likelihood         -12.819   -12.890   -12.978 
## Deviance                25.638    25.779     4.216 
## AIC                     33.638    33.779    35.956 
## BIC                     39.501    39.642    43.285 
## N                       32        32        32     
## ===================================================

Una breve presentación de qué hace mtable se hace en este post: mtable el equivalente en R de estimates table de stata.

No hay comentarios:

Publicar un comentario