> 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
+ "
