24 jul 2019

Modelos de probabilidad (Logit y Probit) binomiales y multinomiales en R


> library(wooldridge)
> library(stargazer)
> library(pscl) #PSCL = Bondad de ajuste: McFadden
> library(mfx, betareg) #efectos marginales
> library(car) #Estadistica, recodificar
> library(dplyr) #Gestión de bases de datos
> #Siempre es importante, realizar un análisis descriptivo preliminar. 
> attach(wage2)

> "WAGE2.DES
+ wage  hours     IQ        KWW       educ      exper     tenure    age      
+ married   black     south     urban   sibs   brthord  meduc     feduc    
+ lwage     
+ Obs:   935
+ 1. wage                     monthly earnings
+ 2. hours                    average weekly hours
+ 3. IQ                       IQ score
+ 4. KWW                      knowledge of world work score
+ 5. educ                     years of education
+ 6. exper                    years of work experience
+ 7. tenure                   years with current employer
+ 8. age                      age in years
+ 9. married                  =1 if married
+ 10. black                    =1 if black
+ 11. south                    =1 if live in south
+ 12. urban                    =1 if live in SMSA
+ 13. sibs                     number of siblings
+ 14. brthord                  birth order
+ 15. meduc                    mother's education
+ 16. feduc                    father's education
+ 17. lwage                    natural log of wage"

> # Modelo de probabilidad de que una persona presente un salario por encima del promedio
> wageHrs <- wage/(hours*4.2)
> SalarioALto <- as.numeric(wage>mean(wage))
> all(as.numeric(wage>mean(wage))==ifelse(wage>mean(wage),1,0))
[1] TRUE
> # Agregar una variable a mi base de datos
> wage2$salarioalto <- SalarioALto
> View(wage2)
> # Estimar los modelos de probabilidad

> "En R las funciones lm y glm, permite obtener los modelos de probabilidad lineal y modelos Probit y Logit. Posteriormente, utilizamos la función Stargazer para imprimir los resultados de los modelos. Tener claro, que salvo en el modelo de probabilidad (modeloPro), estos coeficientes no son interpretables de forma intuitiva. En el caso de los probit y logit, el signo me continúa indicando dirección de impacto en la probabilidad, note que coincide con + el signo del modelo de probabilidad."
> # Modelo de probabilidad
> modeloPro <- lm(SalarioALto~IQ+educ+exper+I(exper^2)+urban+black, data=wage2)
> # Modelo Probit
> modelo.probit <- glm(SalarioALto~IQ+educ+exper+I(exper^2)+urban+black,
+                     data=wage2, family=binomial(link="probit"))
> # Modelo Logit
> modelo.logit <- glm(SalarioALto~IQ+educ+exper+I(exper^2)+urban+black,
+                   data=wage2, family=binomial(link="logit"))
> # names(modelo.logit)
> stargazer(modeloPro, modelo.probit, modelo.logit, type="text",
+           title="Tabla 1. Modelos de probabilidad",
+           digits=2,
+           single.row=FALSE,
+           intercept.bottom=TRUE,
+           df = FALSE)

Tabla 1. Modelos de probabilidad
=================================================
                         Dependent variable:     
                    -----------------------------
                             SalarioALto         
                       OLS     probit   logistic 
                       (1)       (2)       (3)   
-------------------------------------------------
IQ                  0.004***   0.01***   0.02*** 
                     (0.001)   (0.004)   (0.01)  
                                                 
educ                 0.05***   0.14***   0.22*** 
                     (0.01)    (0.03)    (0.04)  
                                                 
exper                 0.02      0.07      0.11   
                     (0.02)    (0.05)    (0.08)  
                                                 
I(exper2)            -0.0000   -0.0002   -0.0004 
                     (0.001)   (0.002)   (0.003) 
                                                 
urban                0.23***   0.68***   1.11*** 
                     (0.03)    (0.10)    (0.17)  
                                                 
black               -0.17***  -0.52***  -0.87*** 
                     (0.05)    (0.15)    (0.25)  
                                                 
Constant            -1.01***  -4.42***  -7.23*** 
                     (0.17)    (0.51)    (0.87)  
                                                 
-------------------------------------------------
Observations           935       935       935   
R2                    0.15                       
Adjusted R2           0.15                       
Log Likelihood                 -565.43   -565.81 
Akaike Inf. Crit.             1,144.86  1,145.63 
Residual Std. Error   0.46                       
F Statistic         27.63***                     
=================================================
Note:                 *p<0.1; **p<0.05; ***p<0.01
> "Intervalos de confianza, sobre los coeficientes estimados"
[1] "Intervalos de confianza, sobre los coeficientes estimados"
> confint(modelo.logit, level = 0.95)
Waiting for profiling to be done...
                  2.5 %       97.5 %
(Intercept) -8.96827199 -5.555087814
IQ           0.00862557  0.032245906
educ         0.14065842  0.305938128
exper       -0.04003188  0.264492650
I(exper^2)  -0.00670212  0.005887174
urban        0.78776839  1.445851370
black       -1.37765309 -0.378928881
> "Medida de bondad de ajuste, a partir de una tabla de clasificación: Aquí se estima la variable independiente a partir del modelo logit, para verificar, que porcentaje de las observaciones han sido predichas correctamente por el modelo. "
> "Agrega la variable yhat a la base de datos: cuando la probabilidad estimada por el Logit es mayor a 0.5, yhat=1, en caso contrario es igual a cero. Luego utilizo la función mean y la propiedad de las variables binarias, para preguntarme en que porcentaje, coinciden el valor estimados con el observado."

> wage2$salarioaltoHat <- round(fitted(modelo.logit))
>  mean(wage2$salarioalto==wage2$salarioaltoHat)
[1] 0.6609626
>  
> # o Cambiar el valor (DEFAULT) a partir de la proporción de 1 en yobs
>  # Promedio de mi población con salario por encima de la media.
>  
> p<- mean(wage2$salarioalto) 
> wage2$salarioaltoHat2 <- ifelse(fitted(modelo.logit) >= (1-p), 1, 0)
> mean(wage2$salarioalto==wage2$salarioaltoHat2)
[1] 0.659893
> "Este análisis nunca es recomendable verlo a nivel agregado, porque el modelo puede tener diversos niveles de eficiencia, en los grupos que se estan considerando."

> #hacemos una tabla
> table(     Obs = wage2$salarioalto,
+       hatLogit = wage2$salarioaltoHat2) ->table1
> "Porcentaje de casos predichos correctamente, cuando yobs=0"

> table1[1,1]/sum(wage2$salarioalto==0)
[1] 0.7825243
> "Porcentaje de casos predichos correctamente, cuando yobs=1. Se verifica, que el modelo predice un porcentaje menor de los datos."

> table1[2,2]/sum(wage2$salarioalto==1)
[1] 0.5095238
> "Pseudo R2"
[1] "Pseudo R2"
> pR2(modelo.logit)
         llh      llhNull           G2     McFadden         r2ML         r2CU 
-565.8149986 -643.2580723  154.8861474    0.1203919    0.1526603    0.2042548 
> # Interpretación
> "Estos se utilizan para facilitar la compresión de los coeficientes estimados. Su interpretación se realiza en torno a la unidad. Con respecto al nivel educativo, el odds-ratio, indica que cada año adicional de educación, se espera un aumento en las probabilidades de tener un salario por encima de la media, igual a 24%, manteniendo el resto de variables constantes."

> coef(modelo.logit)   # Coeficientes del modelo Logit

> exp(coef(modelo.logit)) # odds-ratios

> cbind(coeff=coef(modelo.logit), 
+       odd=exp(coef(modelo.logit)))
                    coeff          odd
(Intercept) -7.2281717932 0.0007258466
IQ           0.0203578903 1.0205665255
educ         0.2225410899 1.2492471501
exper        0.1105312180 1.1168712150
I(exper^2)  -0.0003661957 0.9996338714
urban        1.1125439758 3.0420875560
black       -0.8665315393 0.4204071892
> 
> "Verificar el nivel de significancia del odds-ratio: para ver significancia, verificó si el odds-ratio, es diferente o no de 1."

Waiting for profiling to be done...
                   2.5 %      97.5 %
(Intercept) 0.0001273881 0.003867729
IQ          1.0086628769 1.032771439
educ        1.1510314105 1.357898288
exper       0.9607588131 1.302769848
I(exper^2)  0.9933202891 1.005904538
urban       2.1984847963 4.245465068
black       0.2521696782 0.684594299

> # Presentar odds en tabla de ecuaciones

> logit.odds = exp(coef(modelo.logit))
> stargazer(modelo.logit, coef=list(logit.odds), type = "text")

=============================================
                      Dependent variable:    
                  ---------------------------
                          SalarioALto        
---------------------------------------------
IQ                         1.021***          
                            (0.006)          
                                             
educ                       1.249***          
                            (0.042)          
                                             
exper                      1.117***          
                            (0.077)          
                                             
I(exper2)                  1.000***          
                            (0.003)          
                                             
urban                      3.042***          
                            (0.168)          
                                             
black                       0.420*           
                            (0.254)          
                                             
Constant                     0.001           
                            (0.870)          
                                             
---------------------------------------------
Observations                  935            
Log Likelihood             -565.815          
Akaike Inf. Crit.          1,145.630         
=============================================
Note:             *p<0.1; **p<0.05; ***p<0.01
> # Efectos marginales
> "Estos efectos marginales, me permiten ver el efecto esperado en la probabilidad de Y, cuando cambia una de mis variables X."

> logitMfx<-logitmfx(SalarioALto~IQ+educ+exper+I(exper^2)+urban+black,
+                         data=wage2)
> options(scipen=999) # How to disable scientific notation? 
> logitMfx

Call:
logitmfx(formula = SalarioALto ~ IQ + educ + exper + I(exper^2) + 
    urban + black, data = wage2)

Marginal Effects:
                  dF/dx    Std. Err.       z              P>|z|    
IQ          0.005004668  0.001477740  3.3867          0.0007074 ***
educ        0.054708239  0.010360957  5.2802 0.0000001290215583 ***
exper       0.027172368  0.019048727  1.4265          0.1537338    
I(exper^2) -0.000090023  0.000787632 -0.1143          0.9090028    
urban       0.255938991  0.034690150  7.3779 0.0000000000001609 ***
black      -0.196951878  0.051113564 -3.8532          0.0001166 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dF/dx is for discrete change for the following variables:

[1] "urban" "black"

# Estima probabilidades

> allmeans$predprob <- predict(modelo.logit, newdata=allmeans, type="response")
> allmeans
        IQ     educ    exper I.exper.2. urban black  predprob
1 101.2824 13.46845 11.56364   152.8342     1     1 0.3331872
2 101.2824 13.46845 11.56364   152.8342     0     1 0.1410800
> #Ejemplo 2
> "La probabilidad que Y=1, es 0.33 si el individuo vive en la zona urbana y 0.14 si vive en la zona rural."
> allmeans1 <- data.frame(IQ = c(mean(IQ),mean(IQ)),        
+                        educ = rep(mean(educ),2),  
+                        exper = rep(mean(exper),2), 
+                        "I(exper^2)" = rep(mean(exper^2),2),
+                        urban  =c(1,0),    
+                        black = rep(1,2) )
> allmeans1$predprob <- predict(modelo.logit, newdata=allmeans1, type="response")
> allmeans1
        IQ     educ    exper I.exper.2. urban black  predprob
1 101.2824 13.46845 11.56364   152.8342     1     1 0.3331872
2 101.2824 13.46845 11.56364   152.8342     0     1 0.1410800
> #Ejemplo 3
> "Verificar como cambia la probabilidad de ganar por encima del promedio, para un individuo con IQ promedio, experiencia laboral promedio, que vive en la zona urbana y no es negro. El efecto del cambio en el año educativo no es lineal."
> n<-30
> allmeans2 <- data.frame(IQ = rep(mean(IQ),n),        
+                         educ = 0:(n-1),  
+                         exper = rep(mean(exper),n), 
+                         "I(exper^2)" = rep(mean(exper^2),n),
+                         urban  =rep(1,n),    
+                         black = rep(0,n) )
> allmeans2$predprob <- predict(modelo.logit, newdata=allmeans2, type="response")
> allmeans2
         IQ educ    exper I.exper.2. urban black   predprob
1  101.2824    0 11.56364   152.8342     1     0 0.05601151
2  101.2824    1 11.56364   152.8342     1     0 0.06900880
3  101.2824    2 11.56364   152.8342     1     0 0.08475131
4  101.2824    3 11.56364   152.8342     1     0 0.10368509
5  101.2824    4 11.56364   152.8342     1     0 0.12626520
6  101.2824    5 11.56364   152.8342     1     0 0.15292374
7  101.2824    6 11.56364   152.8342     1     0 0.18402527
8  101.2824    7 11.56364   152.8342     1     0 0.21981082
9  101.2824    8 11.56364   152.8342     1     0 0.26033501
10 101.2824    9 11.56364   152.8342     1     0 0.30540567
11 101.2824   10 11.56364   152.8342     1     0 0.35453912
12 101.2824   11 11.56364   152.8342     1     0 0.40694603
13 101.2824   12 11.56364   152.8342     1     0 0.46156007
14 101.2824   13 11.56364   152.8342     1     0 0.51711265
15 101.2824   14 11.56364   152.8342     1     0 0.57224545
16 101.2824   15 11.56364   152.8342     1     0 0.62564054
17 101.2824   16 11.56364   152.8342     1     0 0.67614258
18 101.2824   17 11.56364   152.8342     1     0 0.72284977
19 101.2824   18 11.56364   152.8342     1     0 0.76516041
20 101.2824   19 11.56364   152.8342     1     0 0.80277415
21 101.2824   20 11.56364   152.8342     1     0 0.83565734
22 101.2824   21 11.56364   152.8342     1     0 0.86398686
23 101.2824   22 11.56364   152.8342     1     0 0.88808692
24 101.2824   23 11.56364   152.8342     1     0 0.90836960
25 101.2824   24 11.56364   152.8342     1     0 0.92528558
26 101.2824   25 11.56364   152.8342     1     0 0.93928741
27 101.2824   26 11.56364   152.8342     1     0 0.95080474
28 101.2824   27 11.56364   152.8342     1     0 0.96022972
29 101.2824   28 11.56364   152.8342     1     0 0.96790997
30 101.2824   29 11.56364   152.8342     1     0 0.97414698
> plot(allmeans2$educ, allmeans2$predprob)
> # - Probabildad en la media para cada individuo de la muestra
> yhat<-modelo.logit$fitted.values
> prob<-exp(yhat)/(1+exp(yhat))
> # ::::::::::::::::::::::::::::::::::::::::::::::::::::
> "Logit-multinomiales: se recodifica la variable salario, para que la misma tenga tres categorías"

> summary(wage)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  115.0   669.0   905.0   957.9  1160.0  3078.0 
> q1<-summary(wage)[2]
> med<-summary(wage)[3]
> "Para recodificar, se utilizan un ifelse anidado: el primero pregunta si wage<q1, en caso afirmativo, el programa colocará un 1, en caso contrario, el programa entrará al segundo if, donde se verifica si el salario está en un rango: wage>q1 & wage<med, en caso de estarlo coloca un 2, y en caso contrario, que significa wage>mediana, coloca un tres"
> wage2$rWage2 <- ifelse(wage<q1, 1,ifelse(wage>q1 & wage<med,2,3))
> wage2$rWage2 <- factor(wage2$rWage2)
> levels(wage2$rWage2)
[1] "1" "2" "3"
> # Estadístico por grupo
> "Cuantas obs. estan en cada grupo"

> table(wage2$rWage2)

  1   2   3 
234 232 469 
> # Estadístico por grupo
> table(wage2$rWage2)

  1   2   3 
234 232 469 
> wage2 %>%  #requiere el dplyr
+   group_by(rWage2) %>%
+   summarize(no_obs = n(),
+             ing.medio=mean(wage),
+             ingRelativ= ing.medio/max(wage))
# A tibble: 3 x 4
  rWage2 no_obs ing.medio ingRelativ
  <fct>   <int>     <dbl>      <dbl>
1 1         234      519.      0.776
2 2         232      792.      0.876
3 3         469     1260.      0.409
> # Logit multinomial
> library(nnet)
> model <- multinom(rWage2 ~ educ, data=wage2)
# weights:  9 (4 variable)
initial  value 1027.202490 
final  value 944.306828 
converged
> summary(model)
Call:
multinom(formula = rWage2 ~ educ, data = wage2)

Coefficients:
  (Intercept)       educ
2  -0.8821114 0.06742855
3  -2.7213052 0.25555737

Std. Errors:
  (Intercept)       educ
2   0.6133235 0.04681189
3   0.5378879 0.04045376

Residual Deviance: 1888.614 
AIC: 1896.614 
> "
+ [,1]       [,2]
+ [1,] -0.8712192 0.06725081
+ [2,] -2.7374661 0.25640808
+ log(pr[y=2]/Pr[y=1])=-0.8712+0.067educ
+ log(pr[y=3]/Pr[y=1])=-2.7374+0.256educ
+ "

Recesión plot en R usando ggplot (recession plot in r)

En el siguiente ejemplo se simular y replica -parcialmente- un ejemplo usado por el FMI para ilustrar la importancia del uso de series de ti...