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