Logit, Probit y MLP con R

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.

Función aggregate: Ejemplos

En R es fácil resumir información a partir de apliación de funciones sobre los subgrupos creados según las variables de agrupación (factor o numérica) por las cuales se desean hacer las agrupaciones. Es posible hacer esto con una gran variedad de funciones (aggregate, split, by, tapply y con otras funciones del paquete plyr), sin embargo este post sólo trata de la función aggregate porque está en el paquete base de R y es muy manejable.

Según la documentación de la función aggregate, esta función divide la muestra en subconjuntos, calcula los estadísticos deseados por cada subconjunto y devuelve el resultado en una forma conveniente, esto significa que el resutado puede ser un vector, un data.frame o una lista según sea más conveniente.

A continuación algunos ejemplos usando la base datos iris.

data(iris)  # cargando la base de datos `iris`
head(iris)  # así lucen los datos.

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

Nota En todo este post se usará la notación de fórmula dentro de la función aggregate donde la lógica es escribir al lado izquierdo del símbolo ~ la variable que se quiere agrupar y al lado derecho, la variable de agrupación: variable_a_ser_agrupada ~ variable_de_agrupación. Así la sintaxis entera será:

aggregate(variable_a_ser_agrupada ~ variable_de_agrupación, FUN=función_deseada, data=base_de_datos)

Ejemplo 1

Calculando la suma de Sepal.Length por Species, como ha de esperarse, el resultado final será un data.frame con la suma total de la variable Sepal.Length por cada especie contenida en la variable Species.

aggregate(Sepal.Length ~ Species, FUN = sum, data = iris)

##      Species Sepal.Length
## 1     setosa        250.3
## 2 versicolor        296.8
## 3  virginica        329.4

Ejemplo 2

Una facilidad que permite la función aggregate es que se pueden resumir más de una variable, veamos la suma de todas las demás varibles según la especie.

aggregate(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, 
    FUN = sum, data = iris)

##      Species Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1     setosa        250.3       171.4         73.1        12.3
## 2 versicolor        296.8       138.5        213.0        66.3
## 3  virginica        329.4       148.7        277.6       101.3

Nótese que cuando se quieren “agregar” muchas variables es necesario usar la función cbind para que funcione. Sin embargo, no es muy práctico escribir el nombre de cada variable como se muestra en el ejemplo anterior, pero como se está utilizando el estilo de “fórmula” para escribir dentro de la función aggregate esta permite reemplazar lo anterior por este nuevo código:

aggregate(. ~ Species, FUN = sum, data = iris)

##      Species Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1     setosa        250.3       171.4         73.1        12.3
## 2 versicolor        296.8       138.5        213.0        66.3
## 3  virginica        329.4       148.7        277.6       101.3

y como ha de esperarse, los resultados son los mismos.

Ejemplo 3

¿Qué tal si tuviéramos más de 1 variable de agregación, por ejemplo Species y Petal.Size?, Esto no representa ningún problema para aggregate, veamos un ejemplo.

Primero se ha de crear la nueva variable Petal.Size que indica si la longitud del pétalo (Petal.Length) es mayor que su valor mediano (4.350 cm) tendrá el valor Big y si es menor o igual, será Small.

iris$Petal.Size <- with(iris, ifelse(Petal.Length > median(Petal.Length), "Big", 
    "Small"))
aggregate(Petal.Length ~ Species + Petal.Size, FUN = sum, data = iris)  # resultado sólo para Petal.Length

##      Species Petal.Size Petal.Length
## 1 versicolor        Big        115.9
## 2  virginica        Big        277.6
## 3     setosa      Small         73.1
## 4 versicolor      Small         97.1

aggregate(. ~ Species + Petal.Size, FUN = sum, data = iris)  # resultado para todos

##      Species Petal.Size Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 versicolor        Big        156.4        72.8        115.9        36.5
## 2  virginica        Big        329.4       148.7        277.6       101.3
## 3     setosa      Small        250.3       171.4         73.1        12.3
## 4 versicolor      Small        140.4        65.7         97.1        29.8

Ejemplo 4 (un truco sencillo)

¿Y si en lugar de sólo querer la suma, quiero además la media y desviación estándar de los datos de Petal.Length agrupados por Species?, Esto se puede hacer sin ningún problema.

resultado <- aggregate(Petal.Length ~ Species, FUN = function(x) c(Suma = sum(x), 
    Media = mean(x), SD = sd(x)), data = iris)
resultado

##      Species Petal.Length.Suma Petal.Length.Media Petal.Length.SD
## 1     setosa           73.1000             1.4620          0.1737
## 2 versicolor          213.0000             4.2600          0.4699
## 3  virginica          277.6000             5.5520          0.5519

A pesar que el resultado mostrado es el adecuado, en realidad, no es muy útil si se desea guardar para luego volverlo a usar, porque en realidad el output de este último ejemplo es un data.frame que debería tener la información calculada, pero si intentamos acceder a esa info a través de su nombre o su posición nos dice que no existe tal información!!! :(

resultado[, "Petal.Length.Suma"]  # Error en `[.data.frame`(resultado, , 3) : undefined columns selected  :(
resultado[, "Petal.Length.Media"]  # Error en `[.data.frame`(resultado, , 3) : undefined columns selected  :(
resultado[, 3]  # Error en `[.data.frame`(resultado, , 3) : undefined columns selected  :(

Pero tranquilos,“calma! que no panda el cúnico”, hay una sencilla solución para ello, basta con convertir esto a un data.frame usado la función do.call

resultado2 <- do.call(data.frame, resultado)

Aparentemente no ha habido ningún cambio respecto a lo anterior, sin embargo, ahora sí se pueden acceder a los datos de las columnas, porque ahora sí existen!!! :D

resultado2[, "Petal.Length.Suma", drop = FALSE]  # Sí funciona!!!
resultado2[, "Petal.Length.Media", drop = FALSE]  # Sí funciona!!!
resultado2[, 3, drop = FALSE]  # Sí funciona!!!

Histogramas con etiquetas de frecuencias

Este sencillo post nace de una pregunta que fue hecha en Stackoverflow. En su momento me pareció interesante para mostrar el uso de ... dentro la definición de nuevas funciones que sólo requieren fijar un par de parámetros y mantener la funcionalidad del resto, en este caso se redefine la función hist, utilizada para estimar y dibujar histogramas, que a partir de ella se crea histL que en esencia hace lo mismo pero adicionalmente crea y pinta las etiquetas (labels) de los conteos de frecuencia o probabilidad, según sea el caso, sobre cada barra del histograma, algo así parecido a lo que hace SPSS, pero sin encasillarlas dentro de cuadros. El nombre histL es debido a que se hacen histogramas con Labels.

Para resumir (lo que ya era corto), histL estima un histograma a partir de datos tal y cual lo hace hist, la diferencia es que esta nueva función añade etiquetas sobre cada barra que indica el conteo o la frecuencia asociada a esa barra particular, se imprime un indicador por cada barra del histograma. Todos los demás parámetros de la función hist son aplicables a histL de manera que la documentación de R disponible para hist también es válida para histL esto es cierto porque los parámetros son pasados a través del operator ... como se muestra a continuación.

Definición de la función histL

# ============================================================================
# Nombre : histL
# Autor : J. Urbina
# Version : 0.01
# Descripción : Histogramas con etiquetas de frecuencias/probabilidades
# ============================================================================


histL <- function(x, ylim = NULL, freq = TRUE, ...) {

    y <- hist(x, plot = FALSE)

    if (freq == TRUE) {
        z <- y$counts
    } else {
        z <- y$density
    }

    if (is.null(ylim)) {
        ylim <- c(0, max(z) + ifelse(is.integer(z), 3, 0.1))
    }


    plot(y, ylim = ylim, freq = freq, ...)
    # parámetros perdeterminados más `...` para aprovechar la funcionalidad ya
    # implementada de `hist` a través de `plot.histogram`

    text(y$mids, z + ifelse(is.integer(z), 1, 0.03), z, cex = 0.75)

}

Ejemplos

Algunos ejemplos que muestran que se pueden pasar todos los argumentos de hist a histL.

# Generando una muestra pseudo-aleatoria a partir de una normal estándar
set.seed(1)
x <- rnorm(100)
# Simple
hist(x)  # sin etiquetas  

histL(x)  # con etiquetas

# se pueden usar todos los demás argumentos tal cual si fuera `hist`
histL(x, main = "Histograma de x con título")

hist(x, freq = FALSE)  # histograma de la densidad de probabilidades sin etiquetas

histL(x, freq = FALSE)  # histograma de la densidad de probabilidades con etiquetas

histL(x, freq = FALSE, main = "Histograma con líndea de densidad \n título en cursiva", 
    cex.main = 1.25, font.main = 3, xlab = "", las = 1, cex.axis = 0.8)
lines(density(x), col = "red")