data.financiero <- read_csv("../files/motral2012.csv",
locale = locale(encoding = "latin1")) %>%
clean_names() %>%
mutate(ahorra_inf = case_when(p14 %in% c(4,5,6,7) ~ 1,
.default = 0),
mujer=ifelse(sex==2,1,0))Respuestas a la tarea 2
Respuestas
Pregunta 1
Use los datos en el archivo motral2012.csv, que incluye una muestra de individuos con sus características socioeconómicas. Nos interesa conocer los factores que afectan la probabilidad de que los individuos tengan ahorros formales. Considere lo siguiente sobre las opciones de ahorro de los entrevistados, contenida en la variable p14:
- p14 igual a 1 significa cuentas de ahorro bancarias
- p14 igual a 2 significa cuenta de inversión bancaria
- p14 igual a 3 significa inversiones en bienes raíces
- p14 igual a 4 significa caja de ahorro en su trabajo
- p14 igual a 5 significa caja de ahorro con sus amigos
- p14 igual a 6 significa tandas
- p14 igual a 7 significa que ahorra en su casa o alcancías
- p14 igual a 8 significa otro lugar
- p14 NA significa que no ahorra
[2 puntos] Comience generando una variable binaria ahorra_inf que tome el valor de 1 para las personas que ahorran en instrumentos informales y 0 en otro caso. Se consideran instrumentos informales las cajas de ahorro en el trabajo o amigos, las tandas y el ahorro en casa o alcancías . Construya también la variable mujer que tome el valor de 1 cuando sex toma el valor de 2 y 0 en otro caso. Luego, estime un modelo de probabilidad lineal que relacione ahorra_inf como variable dependiente con eda (edad), anios_esc (años de escolaridad) y mujer. Reporte los errores que asumen homocedasticidad y los errores robustos a heteroscedasticidad. ¿Qué observa respecto a los errores y por qué sucede?
Generamos variables:
Estimamos el modelo lineal y obtenemos la matriz de varianzas robusta usando vcovHC:
summary(reg.lineal <- lm(ahorra_inf ~ eda + anios_esc + mujer, data = data.financiero))Call: lm(formula = ahorra_inf ~ eda + anios_esc + mujer, data = data.financiero) Residuals: Min 1Q Median 3Q Max -0.3365 -0.2580 -0.1844 -0.1092 0.9467 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.4519266 0.0254686 17.744 <2e-16 *** eda -0.0063753 0.0005491 -11.610 <2e-16 *** anios_esc -0.0006780 0.0012186 -0.556 0.578 mujer -0.0006674 0.0113305 -0.059 0.953 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.4107 on 5260 degrees of freedom Multiple R-squared: 0.02507, Adjusted R-squared: 0.02452 F-statistic: 45.09 on 3 and 5260 DF, p-value: < 2.2e-16#Matriz robusta v_rob <- vcovHC(reg.lineal, type = "HC0") se_rob <- sqrt(diag(v_rob))Presentamos usando modelsummary.
modelsummary(list(reg.lineal, reg.lineal), vcov = c("iid", "HC0"), stars = c('*'=.1, '**'=.05, '***'=.01))tinytable_fp0cpcl5vevz02j6o6gl (1) (2) * p < 0.1, ** p < 0.05, *** p < 0.01 (Intercept) 0.452*** 0.452*** (0.025) (0.029) eda -0.006*** -0.006*** (0.001) (0.001) anios_esc -0.001 -0.001 (0.001) (0.002) mujer -0.001 -0.001 (0.011) (0.011) Num.Obs. 5264 5264 R2 0.025 0.025 R2 Adj. 0.025 0.025 AIC 5575.3 5575.3 BIC 5608.1 5608.1 Log.Lik. -2782.626 -2782.626 F 45.091 46.854 RMSE 0.41 0.41 Std.Errors IID HC0 En en problema binario, la media condicional está bien planteada. Y dado que la densidad pertenece a la familia lineal exponencial, basta con que la media condicional esté bien planteada para la consistencia de \(\beta\). Sin embargo, en el caso de los modelos binarios, la varianza condicional también está siempre bien planteada, pues \(V(y)=p(1-p)\). Esto implica que no hay ninguna ganancia en usar la matriz de varianzas robustas. Ver CT, p. 469 para una discusión.
[3 puntos] ¿Cuál es el efecto en la probabilidad de ahorrar informalmente si los años de educación se incrementan en una unidad, pasando de 4 a 5 años de educación?
Una reducción de 0.1 puntos porcentuales (0.001/100), estadísticamente no significativa.
[2 puntos] Realice una prueba de significancia conjunta de eda y anios_esc. ¿Qué concluye?
Podemos usar la función linearHypothesis:
car::linearHypothesis(reg.lineal, c("eda=0", "anios_esc=0"))Linear hypothesis test Hypothesis: eda = 0 anios_esc = 0 Model 1: restricted model Model 2: ahorra_inf ~ eda + anios_esc + mujer Res.Df RSS Df Sum of Sq F Pr(>F) 1 5262 909.91 2 5260 887.14 2 22.773 67.512 < 2.2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Concluimos que no hay evidencia para afirmar que \(\beta_{eda}=\beta_{anios\_esc}=0\).
[3 puntos] Estime un modelo logit relacionando las mismas variables. Use la función avg_slopes del paquete marginaleffects para obtener los efectos marginales promedio de un cambio en cada uno de los regresores. ¿Por qué difiere la magnitud de este efecto marginal con respecto a la parte b.?
Estimamos el modelo probit:
reg.logit <- glm(ahorra_inf ~ eda + anios_esc + mujer, family = binomial(link = "logit"), data = data.financiero) summary(reg.logit)$coefEstimate Std. Error z value Pr(>|z|) (Intercept) 0.081021587 0.150793629 0.53730113 5.910596e-01 eda -0.038339619 0.003385824 -11.32357193 1.003001e-29 anios_esc -0.003787198 0.007627194 -0.49653889 6.195143e-01 mujer -0.001539617 0.067250054 -0.02289392 9.817349e-01Noten que el signo de los coeficientes coinciden con el promedio de los efectos marginales:
avg_slopes(reg.logit)Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % anios_esc dY/dX -0.000638 0.001285 -0.4966 0.619 0.7 -0.00316 eda dY/dX -0.006459 0.000554 -11.6533 <0.001 101.8 -0.00755 mujer 1 - 0 -0.000259 0.011331 -0.0229 0.982 0.0 -0.02247 97.5 % 0.00188 -0.00537 0.02195 Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high Type: responseEl promedio del efecto marginal de un cambio en los años de educación es de una reducción de 0.06 puntos porcentuales.
[2 puntos] Ahora estime el efecto marginal en la media para eda y anios_esc y para los hombres, usando la función slopes. ¿Por qué difiere la magnitud de este efecto marginal respecto a la parte b. y la d.?
Para obtener los efectos marginales evaluados en algún valor \(X_i\) de los covariables, debemos especificar estos valores usando datagrid:
avg_slopes(reg.logit, newdata = datagrid(eda = mean(data.financiero$eda), anios_esc = mean(data.financiero$anios_esc), mujer = 1))Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % anios_esc dY/dX -0.000639 0.001287 -0.4963 0.620 0.7 -0.00316 eda dY/dX -0.006465 0.000575 -11.2472 <0.001 95.1 -0.00759 mujer 1 - 0 -0.000260 0.011346 -0.0229 0.982 0.0 -0.02250 97.5 % 0.00188 -0.00534 0.02198 Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high Type: responseEl efecto marginal de un cambio en los años de escolaridad evaluados en la media de los años de educación y edad, para las mujeres, es de -0.06 puntos. Esto difiere del modelo lineal porque en el modelo lineal los efectos marginales son constantes, mientras que los efectos marginal del modelo no lineal dependen del punto de evaluación. También difiere de los efectos marginales promedio pues aquí solo hemos calculado el efecto marginal una sola vez, para un valor \(X_i\), mientras que el promedio de efectos marginales implica calcular el efecto marginal para cada individuo y luego obtener el promedio.
En clase les pregunté cómo estimarían el error estándar de los cambios marginales y brevemente mencioné que una forma muy usada es el Método Delta, el cual se basa en que los efectos marginales son funciones no lineales de los parámetros. Esto es lo que efectivamente se usa en la función avg_slopes para obtener los errores estándar y los intervalos de confianza. Aquí pueden leer al respecto.
[3 puntos] Provea una expresión para la maginitud de:
\[\frac{\frac{\partial P(y=1)}{\partial \; anios\_esc}}{\frac{\partial P(y=1)}{\partial \; eda}}\]
La razón de efectos marginales es la razón de coeficientes:
\[\frac{\frac{\partial P(y=1)}{\partial \; anios\_esc}}{\frac{\partial P(y=1)}{\partial \; eda}}=\frac{\beta_{anios\_esc}}{\beta_{eda}}=\frac{0.0038 }{0.0383}\]
Pregunta 2
Ahora estimará un modelo multinomial empleando los mismos datos en motral2012.csv. El propósito será ahora estudiar los factores relevantes para predecir la forma de ahorro que tienen las personas que ahorran.
[2 puntos] Genere una variable categórica llamada ahorro que sea igual a 1 cuando p14 sea igual a 1 o 2, igual a 2 cuando p14 sea igual a 7, e igual a 3 cuando p14 sea igual a 3, 4, 5, 6 u 8. Haga que esa variable sea missing cuando p14 sea missing. Posteriormente, convierta esta nueva variable en una de factores de forma que el valor 1 tenga la etiqueta “Banco”, el valor 2 tenga la etiqueta “Casa” y el valor 3 tenga la etiqueta “Otro”.
Construimos la variable dependiente:
data.financiero <- read_csv("../files/motral2012.csv", locale = locale(encoding = "latin1")) %>% clean_names() %>% mutate(ahorro=NA) %>% mutate(ahorro=ifelse(p14%in%c(1,2),1,ahorro)) %>% mutate(ahorro=ifelse(p14==7,2,ahorro)) %>% mutate(ahorro=ifelse(p14%in%c(3,4,5,6,8),3,ahorro)) %>% mutate(ahorro=factor(ahorro, levels=c(1,2,3), labels=c("Banco","Casa","Otro"))) %>% mutate(mujer=ifelse(sex==2,1,0))[4 puntos] Estime un modelo logit multinomial (regresores invariantes a la alternativa) con la opción de ahorro como variable dependiente y los mismos regresores de la pregunta 1. Hay varios paquetes para hacer esto, pero recomiendo usar la función multinom del paquete nnet. ¿Qué puede decir sobre el coeficiente de años de educación en la alternativa “Casa”?
Usamos multinom para estimar el modelo logit multinomial:
multilogit <- nnet::multinom(ahorro~ eda + anios_esc + mujer, data=data.financiero)# weights: 15 (8 variable) initial value 2727.854313 iter 10 value 2546.085070 final value 2545.712541 convergedsummary(multilogit)Call: nnet::multinom(formula = ahorro ~ eda + anios_esc + mujer, data = data.financiero) Coefficients: (Intercept) eda anios_esc mujer Casa 3.026880 -0.052310196 -0.14829719 0.09265405 Otro 0.206704 -0.003501367 -0.04628175 -0.06459305 Std. Errors: (Intercept) eda anios_esc mujer Casa 0.2487107 0.005207439 0.01355319 0.09956544 Otro 0.2476498 0.004964123 0.01285100 0.09995715 Residual Deviance: 5091.425 AIC: 5107.425En el logit multinominal (regresores invariantes) el coeficiente se interpreta con respecto a una categoría base. En este caso, la categoría base es Banco. El modelo implica que la probabilidad de ahorrar en casa disminuye con un año más de educación, en comparación con la probabilidad de ahorrar en el banco. En particular, sabemos que podemos escribir el log del cociente de la probabilidad de las categorías \(j\) y \(k\) sean escogidas, normalizando \(k\) a ser la base, como:
\[\ln\left(\frac{P(y=Casa)}{P(y=Banco)}\right)=x'\beta=\beta_0+\beta_1 edad + \beta_2 educación + \beta_3 mujer \]
Es decir, un año más de educación se asocia con una reducción en el log de la razón de momios de 0.15.
[4 puntos] Calcule los efectos marginales promedio sobre la probabilidad de ahorrar en el banco. Al considerar el cambio en la probabilidad para el caso de las mujeres (cuando la variable mujer pasa de 0 a 1), ¿de qué tamaño es el efecto predicho en la probabilidad de ahorrar en el banco?
Usamos avg_slopes:
avg_slopes(multilogit)Group Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % Banco anios_esc dY/dX 0.02285 0.002381 9.595 <0.001 70.0 0.01818 Banco eda dY/dX 0.00657 0.000935 7.027 <0.001 38.8 0.00474 Banco mujer 1 - 0 -0.00343 0.019432 -0.177 0.860 0.2 -0.04152 Casa anios_esc dY/dX -0.02512 0.002231 -11.262 <0.001 95.3 -0.02949 Casa eda dY/dX -0.00982 0.000860 -11.422 <0.001 98.0 -0.01151 Casa mujer 1 - 0 0.02271 0.017662 1.286 0.199 2.3 -0.01191 Otro anios_esc dY/dX 0.00228 0.002174 1.046 0.295 1.8 -0.00199 Otro eda dY/dX 0.00325 0.000842 3.863 <0.001 13.1 0.00160 Otro mujer 1 - 0 -0.01927 0.017549 -1.098 0.272 1.9 -0.05367 97.5 % 0.02751 0.00840 0.03465 -0.02075 -0.00814 0.05732 0.00654 0.00490 0.01512 Columns: term, group, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high Type: probsEl efecto de ser mujer es de una reducción de 0.3 puntos en la probabilidad de ahorrar en el banco al estimar el promedio de los efectos marginales.
[3 puntos] Calcule los cocientes de riesgo relativo (relative risk ratios o RRR). ¿Qué significa el hecho de que el RRR asociado a ser mujer sea mayor que 1 en la alternativa “Casa”?
(multilogit_rrr = exp(coef(multilogit)))(Intercept) eda anios_esc mujer Casa 20.632752 0.9490344 0.8621748 1.0970821 Otro 1.229619 0.9965048 0.9547729 0.9374489Los coeficientes en forma de RRR tienen la interpretación del cambio en el riesgo relativo que una categoría sea elegida con relación al riesgo de escoger la categoría base. En este caso, el ser mujer está asociado con una probabilidad de ahorrar en “Casa” 1.097 veces mayor de que la de ahorrar en “Banco”.
[2 puntos] Estime nuevamente el modelo, pero ahora, especifique que la alternativa “Casa” sea la alternativa base. ¿Cómo es el RRR de la edad en la alternativa “Banco”? ¿Es esto congruente con lo que obtuvo en la parte d. de esta pregunta?
Primero tenemos que cambiar la base. Para esto hacemos uso de que ahorro es una variable de factores. Luego estimamos:
data.financiero <- data.financiero %>% mutate(ahorro = relevel(ahorro, ref = "Casa")) multilogit2 <- nnet::multinom(ahorro ~ eda + anios_esc + mujer, data=data.financiero)# weights: 15 (8 variable) initial value 2727.854313 iter 10 value 2552.696634 final value 2545.712541 convergedObtenemos el RRR:
(multilogit2_rrr = exp(coef(multilogit2)))(Intercept) eda anios_esc mujer Banco 0.04846638 1.053703 1.159857 0.9115111 Otro 0.05959489 1.050020 1.107401 0.8544952Al cambiar la categoría base a Casa solo se modifica la interpretación relativa. En la parte d. el RRR de la edad para la opción de Casa era 0.949, es decir, si la edad se incrementa en una unidad, la probabilidad de ahorrar en Casa es 0.949 veces la de ahorrar en Banco. Con la nueva categoría base, el RRR de la edad para ahorrar en Banco es 1.054, es decir, si la edad se incrementa en un año, la probabilidad de ahorrar en Banco es 1.054 veces más probable que la probabilidad de ahorrar en Casa. La parte d. implica que \(P(Casa)=0.949(Banco)\). Mientras que estimando el modelo con la nueva categoría, \(P(Banco)=1.054(Casa)\), o \(P(Casa)=1/1.054(Banco)\). Empleando todos los decimales en R se puede notar que 1/1.054≅0.949 Ambos resultados son consistentes.
Pregunta 3: modelo Poisson inflado en cero
Otra manera de resolver el problema del exceso de ceros que a veces nos molesta en los modelos Poisson es usar un modelo Poisson inflado en cero (CT, p. 681). La idea es introducir un proceso binario con densidad \(f_1(\cdot)\) para modelar la probabilidad de que \(y=0\) y luego una densidad de conteo \(f_2(\cdot)\). Si el proceso binario toma el valor de 0, con probabilidad \(f_1(0)\), entonces \(y=0\), pero si el proceso binario toma el valor de 1, entonces \(y={0,1,2,\ldots}\). Note que podemos entonces observar ceros por dos razones, por el proceso binomial o por el conteo.
Un modelo inflado en cero tendrá como densidad:
\[ g(y)= \begin{cases} f_1(0)+(1-f_1(0))f_2(0) & \text{si }y=0 \\ (1-f_1(0))f_2(y)& \text{si }y\geq 1 \end{cases} \]
Considere la variable aleatoria \(Y\) con observaciones iid que sigue una distribución Poisson con parámetro \(\lambda\). Y considere una variable un proceso binomial tal que \(\pi\) es la probabilidad de que el conteo no se realice. Entonces:
\[ g(y)= \begin{cases} \pi+(1-\pi)f_2(0) & \text{si }y=0 \\ (1-\pi)f_2(y)& \text{si }y\geq 1 \end{cases} \]
[4 puntos] Termine de especializar la expresión anterior unsando la distribución Poisson para \(f_2(\cdot)\) para obtener la función de masa de probabilidad del modelo Poisson inflado en cero \(g(y|\lambda, \pi)\).
En el caso particular de un modelo Poisson, sabemos que \(f_2(0)=P(Y=0)=exp(-\lambda)\). Definiendo la probabilidad de observar un conteo cero como \(\pi\), la función de masa de probabilidad del modelo inflado en cero es:
\[g(y)= \begin{cases} \pi+(1-\pi)exp(-\lambda) \quad\text{si }y=0 \\ (1-\pi)\frac{\lambda^y exp(-\lambda)}{y!} \quad\text{si } y \geq1 \\ \end{cases} \]
[5 puntos] Provea una expresión para la función de verosimilitud \(L(\lambda,\pi)=\prod_{i=1}^N g(y_i|\lambda, \pi)\). Una sugerencia para simplificar sus cálculos es definir una variable \(X\) igual al numero de veces que \(Y_i\) que toma el valor de cero.
La función de verosimilitud del problema es:
\[L(\pi,\lambda,y_i)=\prod_i P(Y_i=y_i)\]
Con las formas específicas para el modelo Poisson inflado en cero:
\[L(\pi,\lambda,y_i)=\prod_{i\in y_i=0}\left(\pi+(1-\pi)exp(-\lambda) \right) \prod_{i\in y_i>0}\left((1-\pi)\frac{\lambda^{y_i} exp(-\lambda)}{y!}\right)\]
Haciendo \(X\) el número de veces que \(y_i\) toma el valor de cero, el primer producto es \(\left(\pi+(1-\pi)exp(-\lambda) \right)\) elevado a la potencia \(X\).
¿Cuántos términos distintos de cero quedan? Quedan \(n-X\). El segundo producto en la verosimilitud es:
\[\left((1-\pi)exp(-\lambda)\right)^{n-X}\frac{\lambda^{\sum_i y_i}}{\prod_{i\in y_i>0} y!}\]
La verosimilitud es por tanto:
\[L(\pi,\lambda,y_i)=\left(\pi+(1-\pi)exp(-\lambda) \right)^X \left((1-\pi)exp(-\lambda)\right)^{n-X}\frac{\lambda^{\sum_i y_i}}{\prod_{i\in y_i>0} y!}\]
[3 puntos] Provea una expresión para la log verosimilitud del problema, \(\mathcal{L}(\lambda,\pi)\).
Dada la verosimilitud planteada en la parte anterior, la log verosimilitud es:
\[\mathcal{L}(\pi,\lambda,y_i)=X\ln \left(\pi+(1-\pi)exp(-\lambda) \right)+(n-X)\ln(1-\pi)-(n-X)\lambda+n\bar{Y}\ln (\lambda)- \ln\left(\prod_{i\in y_i>0} y! \right)\]
[3 puntos] Obtenga las condiciones de primer orden que caracterizan la solución del problema de máxima verosimilitud, derivando la log verosimilitud con respecto a \(\lambda\) y a \(\pi\).
Tenemos dos parámetros, así que tenemos dos condiciones de primer orden. Derivando la log verosimilitud con respecto a \(\pi\) obtenemos:
\[\frac{\partial \mathcal{L}}{\partial \pi}=\frac{X}{\pi+(1-\pi)exp(-\lambda)}(1-exp(-\lambda))-\frac{n-X}{1-\pi}=0\]
La primer condición (A) es:
\[\frac{X(1-exp(-\lambda))(1-\pi)}{\pi+(1-\pi)exp(-\lambda)}=n-X\quad\quad\ldots(A)\]
Ahora derivando la log verosimilitud con respecto a \(\lambda\):
\[\frac{\partial \mathcal{L}}{\partial \lambda}=-\frac{X}{\pi+(1-\pi)exp(-\lambda)}(1-\pi)exp(-\lambda)-(n-X)+\frac{n\bar{Y}}{\lambda}=0\]
La segunda condición (B) es:
\[\frac{X(1-\pi)exp(-\lambda)}{\pi+(1-\pi)exp(-\lambda)}+(n-X)=\frac{n\bar{Y}}{\lambda}\quad\quad\ldots(B)\]
\((\hat{\pi}_{MV},\hat{\lambda}_{MV})\) son los valores de los parámetros que resulven el sistema dado por (A) y (B).
Pregunta 4
Use los datos phd_articulos.csv, los cuales contienen información sobre el número de artículos publicados para una muestra de entonces estudiantes de doctorado. Nuestra variable de interés será el número de artículos art.
[4 puntos] Estime un modelo Poisson que incluya variables dicotómicas para estudiantes mujeres (female) y para estudiantes casadas o casados (married), el número de hijos mejores de cinco años (kid5), el ranking de prestigio del doctorado (phd) y el número de artículos publicados por su mentor (mentor). Realice la estimación de la matriz de varianzas primero a partir de la varianza teórica que resulta de la igualdad de la matriz de información y luego usando una matriz de sándwich. Interprete los coeficientes estimados.
data.phd<-read_csv("../files/phd_articulos.csv", locale = locale(encoding = "latin1")) data.phd <- data.phd %>% mutate(female=factor(female, levels=c('Male','Female'))) mpoisson <- glm(art ~ factor(female) + factor(married) + kid5 + phd + mentor, family="poisson", data=data.phd) summary(mpoisson)Call: glm(formula = art ~ factor(female) + factor(married) + kid5 + phd + mentor, family = "poisson", data = data.phd) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.459860 0.093335 4.927 8.35e-07 *** factor(female)Female -0.224594 0.054613 -4.112 3.92e-05 *** factor(married)Single -0.155243 0.061374 -2.529 0.0114 * kid5 -0.184883 0.040127 -4.607 4.08e-06 *** phd 0.012823 0.026397 0.486 0.6271 mentor 0.025543 0.002006 12.733 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 1817.4 on 914 degrees of freedom Residual deviance: 1634.4 on 909 degrees of freedom AIC: 3314.1 Number of Fisher Scoring iterations: 5Presentamos errores heterocedásticos y robustos a la heterocedasticidad. Aquí les muestro otro paquete que puede servirles para presentar resultados en trabajos y tesinas, alterntivo a stargazer, modelsummary:
modelsummary(list(mpoisson, mpoisson), vcov = list('classical', 'robust'), stars = c('***' = 0.01, '**' = 0.05, '*' = 0.1))tinytable_9o8r9e5l8tfbm2rxphd5 (1) (2) * p < 0.1, ** p < 0.05, *** p < 0.01 (Intercept) 0.460*** 0.460*** (0.093) (0.151) factor(female)Female -0.225*** -0.225*** (0.055) (0.072) factor(married)Single -0.155** -0.155* (0.061) (0.083) kid5 -0.185*** -0.185*** (0.040) (0.057) phd 0.013 0.013 (0.026) (0.044) mentor 0.026*** 0.026*** (0.002) (0.004) Num.Obs. 915 915 AIC 3314.1 3314.1 BIC 3343.0 3343.0 Log.Lik. -1651.056 -1651.056 F 43.333 14.855 RMSE 1.84 1.84 Std.Errors IID HC3 Para las variables continuas, como el número de artículos publicados por el mentor, la interpretación es el cambio en el log conteo esperado. En este caso, un artículo más publicado por el mentor incrementa el log conteo esperado en 0.026. También sabemos que los coeficientes tienen una interpretación de semielasticidad; en este caso, la semielasticidad del conteo con respecto al número de artículos publicados es 0.026. Para las variables dicotómicas, por ejemplo female, la interpretación es la diferencia entre el log conteo esperado entre mujeres y la categoría base (hombres).
[3 puntos] Obtenga la razón de tasas de incidencia (IRR) para los coeficientes e interprete los resultados.
exp(summary(mpoisson)$coef)Estimate Std. Error z value Pr(>|z|) (Intercept) 1.5838526 1.097829 1.379638e+02 1.000001 factor(female)Female 0.7988403 1.056132 1.636793e-02 1.000039 factor(married)Single 0.8562068 1.063297 7.970295e-02 1.011490 kid5 0.8312018 1.040943 9.977222e-03 1.000004 phd 1.0129051 1.026749 1.625407e+00 1.872246 mentor 1.0258718 1.002008 3.386456e+05 1.000000Aunque esto también puede hacerse directamente en modelsummary:
modelsummary(list(mpoisson, mpoisson), exponentiate = TRUE, vcov = list('classical', 'robust'), stars = c('***' = 0.01, '**' = 0.05, '*' = 0.1))tinytable_t66btvk6013qfhtg2axp (1) (2) * p < 0.1, ** p < 0.05, *** p < 0.01 (Intercept) 1.584*** 1.584*** (0.148) (0.240) factor(female)Female 0.799*** 0.799*** (0.044) (0.058) factor(married)Single 0.856** 0.856* (0.053) (0.071) kid5 0.831*** 0.831*** (0.033) (0.047) phd 1.013 1.013 (0.027) (0.045) mentor 1.026*** 1.026*** (0.002) (0.004) Num.Obs. 915 915 AIC 3314.1 3314.1 BIC 3343.0 3343.0 Log.Lik. -1651.056 -1651.056 F 43.333 14.855 RMSE 1.84 1.84 Std.Errors IID HC3 La interpretación de los coeficientes se vuelve más sencilla usando irr. Para la variable continua mentor, un artículo más publicado por el mentor está asociado con 1.026 veces más artículos publicados por el estudiante, es decir, un 2.6% más artículos. En cambio, la variable dicotómica para mujeres indica que las mujeres publican 0.8 veces el número de artículos que los hombres.
[2 puntos] Considere ahora que las mujeres han tenido carreras profesionales más cortas que los hombres, es decir, han estado menos expuestas a la ocurrencia de los eventos publicar. Incorpore esto al análisis y reinterprete los resultados. Pista: explore la opción offeset en glm de R. La columna profage mide la duración efectiva de las carreras profesionales de cada individuo.
El razonamiento es que ahora queremos conocer cuál es la tasa de publicación, es decir, \(art/profage\). Pero como nuestro podemos Poisson solo puede manejar conteos, podemos modificar el modelo para pasar la edad de la carrera del lado derecho:
\[\begin{aligned}ln(art/profage)&=x'\beta \\ ln(art)&=x'\beta+\ln(profage) \end{aligned}\]
mpoisson_duracion <- glm(art ~ factor(female) + factor(married) + kid5 + phd + mentor, offset = log(profage), family="poisson", data=data.phd) summary(mpoisson_duracion)$coefEstimate Std. Error z value Pr(>|z|) (Intercept) -2.95404558 0.093812104 -31.4889600 1.230266e-217 factor(female)Female 0.45874678 0.054721432 8.3833109 5.145931e-17 factor(married)Single -0.15598278 0.061347334 -2.5426171 1.100257e-02 kid5 -0.18643454 0.040135522 -4.6451256 3.398696e-06 phd 0.01801602 0.026428953 0.6816773 4.954430e-01 mentor 0.02573493 0.002001731 12.8563329 7.924799e-38Hasta ahora hemos asumido que cada individuo ha estado “en riesgo” de publicar por el mismo periodo de tiempo, lo cual puede ser no cierto si, por ejemplo, algunos estudiantes se graduaron antes, o si otros han tenido pausas en sus carreras. Al controlar por el hecho de que las mujeres han tenido carreras más cortas, la variable female deja de ser negativa y se convierte en positiva. Las mujeres publican más que los hombres al tomar en cuenta la duración de las carreras.
Comparando los tres modelos:
modelsummary(list(mpoisson, mpoisson, mpoisson_duracion), vcov = list('classical', 'robust', 'robust'), stars = c('***' = 0.01, '**' = 0.05, '*' = 0.1))tinytable_f8i8yavwm2trbbd1swxy (1) (2) (3) * p < 0.1, ** p < 0.05, *** p < 0.01 (Intercept) 0.460*** 0.460*** -2.954*** (0.093) (0.151) (0.155) factor(female)Female -0.225*** -0.225*** 0.459*** (0.055) (0.072) (0.073) factor(married)Single -0.155** -0.155* -0.156* (0.061) (0.083) (0.083) kid5 -0.185*** -0.185*** -0.186*** (0.040) (0.057) (0.057) phd 0.013 0.013 0.018 (0.026) (0.044) (0.045) mentor 0.026*** 0.026*** 0.026*** (0.002) (0.004) (0.005) Num.Obs. 915 915 915 AIC 3314.1 3314.1 3322.8 BIC 3343.0 3343.0 3351.7 Log.Lik. -1651.056 -1651.056 -1655.393 F 43.333 14.855 20.311 RMSE 1.84 1.84 1.85 Std.Errors IID HC3 HC3 [2 puntos] Implemente la prueba de dispersión de Cameron y Trivedi (1990) usando una regresión auxiliar y los coeficientes estimados en la parte a. ¿Qué concluye?
Seguimos a CT, p671 y construimos:
\[\frac{(y_i-\hat{\mu}_i)^2}{\hat{\mu}_i}=\alpha\frac{ g(\hat{\mu}_i)}{\hat{\mu}_i}+u_i\] Creamos el lado izquierdo:
data.phd <- data.phd %>% mutate(xb_hat = predict(mpoisson), mu_hat = exp(xb_hat), lado_izq = (art-mu_hat)^2/mu_hat)Noten que si especificamos \(g(\hat{\mu}_i)=\hat{\mu}^2_i\), el lado derecho simplemente es \(\alpha \hat{\mu}_i+u_i\). Estimamos entonces la regresión, sin constante:
Corremos la regresión:
summary(lm(lado_izq ~ -1 + mu_hat, data = data.phd))Call: lm(formula = lado_izq ~ -1 + mu_hat, data = data.phd) Residuals: Min 1Q Median 3Q Max -5.570 -1.431 -0.798 -0.027 69.967 Coefficients: Estimate Std. Error t value Pr(>|t|) mu_hat 1.02020 0.08866 11.51 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 4.881 on 914 degrees of freedom Multiple R-squared: 0.1265, Adjusted R-squared: 0.1256 F-statistic: 132.4 on 1 and 914 DF, p-value: < 2.2e-16El coeficiente sobre \(\alpha\) es estadísticamente significativo, sugiriendo una relación entre la media y la varianza.
[5 puntos] Emplee ahora un modelo negativo binomial con sobredispersión cuadrática en la media para estimar la relación entre el número de artículos publicados y las variables explicativas antes enumeradas. Interprete el coeficiente asociado al número de hijos y a la variable dicotómica para estudiantes mujeres. ¿Qué puede decir sobre la significancia del \(\alpha\) estimado?
mnb2 <- MASS::glm.nb(art ~ factor(female) + factor(married) + kid5 + phd + mentor, data = data.phd) summary(mnb2)Call: MASS::glm.nb(formula = art ~ factor(female) + factor(married) + kid5 + phd + mentor, data = data.phd, init.theta = 2.264387695, link = log) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.406633 0.125778 3.233 0.001225 ** factor(female)Female -0.216418 0.072636 -2.979 0.002887 ** factor(married)Single -0.150489 0.082097 -1.833 0.066791 . kid5 -0.176415 0.052813 -3.340 0.000837 *** phd 0.015271 0.035873 0.426 0.670326 mentor 0.029082 0.003214 9.048 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for Negative Binomial(2.2644) family taken to be 1) Null deviance: 1109.0 on 914 degrees of freedom Residual deviance: 1004.3 on 909 degrees of freedom AIC: 3135.9 Number of Fisher Scoring iterations: 1 Theta: 2.264 Std. Err.: 0.271 2 x log-likelihood: -3121.917Ponemos todo junto:
modelsummary(list(mpoisson, mpoisson, mpoisson_duracion, mnb2), vcov = list('classical', 'robust', 'robust', 'robust'), stars = c('***' = 0.01, '**' = 0.05, '*' = 0.1))tinytable_8r6ak5uj72zjcl881ed0 (1) (2) (3) (4) * p < 0.1, ** p < 0.05, *** p < 0.01 (Intercept) 0.460*** 0.460*** -2.954*** 0.407*** (0.093) (0.151) (0.155) (0.135) factor(female)Female -0.225*** -0.225*** 0.459*** -0.216*** (0.055) (0.072) (0.073) (0.071) factor(married)Single -0.155** -0.155* -0.156* -0.150* (0.061) (0.083) (0.083) (0.081) kid5 -0.185*** -0.185*** -0.186*** -0.176*** (0.040) (0.057) (0.057) (0.053) phd 0.013 0.013 0.018 0.015 (0.026) (0.044) (0.045) (0.038) mentor 0.026*** 0.026*** 0.026*** 0.029*** (0.002) (0.004) (0.005) (0.003) Num.Obs. 915 915 915 915 AIC 3314.1 3314.1 3322.8 3135.9 BIC 3343.0 3343.0 3351.7 3169.6 Log.Lik. -1651.056 -1651.056 -1655.393 -1560.958 F 43.333 14.855 20.311 20.935 RMSE 1.84 1.84 1.85 1.86 Std.Errors IID HC3 HC3 HC3 A diferencia de otros paquetes, glm.nb reporta \(\theta=1/\alpha\):
(alpha <- 1/summary(mnb2)$theta)[1] 0.4416205Este es el modelo NB2 visto en clase y la forma más usada para implementar un modelo negativo binomial. Se asume una sobredispersión cuadrática en la media, con la varianza parametrizada usando \(\alpha\). La interpretación de los coeficientes se mantiene con respecto al modelo Poisson. Los coeficientes tienen magnitudes similares, pero se prefiere el modelo NB2 si el propósito es pronóstico pues toma en cuenta la sobredispersión y le da suficiente flexibilidad a la varianza para depender de manera cuadrática de la media.
Un poco más cuidado hay que poner en \(\alpha\). En este caso, \(\hat{\alpha}=0.44\). Pero noten que lo que se reporta es el error estándar de \(\theta\). Como platicamos en clase, con un estadístico podemos hacer un test y obtener un valor \(p\), pero una función no lineal del mismo puede que no tenga el mismo valor \(p\). Esto ocurre aquí, deberíamos recurrir al método delta para calcular el error estándar de \(\alpha\).
Pregunta 5
Retome los datos del archivo motral2012.csv usado en la Tarea 1. Estimará un modelo Tobit para explicar los factores que afectan la oferta laboral femenina. En este archivo de datos la variable hrsocup registra las horas trabajadas a la semana.
[2 punto] ¿Qué proporción de la muestra femenina reporta horas trabajadas iguales a cero?
Si hacemos una dummy de horas positivas, al sacarle la media obtenemos la proporción.
data.salarios<-read_csv("../files/motral2012.csv", locale = locale(encoding = "latin1")) data.salarios <- data.salarios %>% filter(sex==2) %>% mutate(zerohrs=ifelse(hrsocup==0,1,0)) summary(data.salarios$zerohrs)Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0000 0.0000 0.0000 0.3528 1.0000 1.0000El 35% de las observaciones tienen cero horas trabajadas.
[3 puntos] Se desea estimar el efecto de los años de educación (anios_esc) sobre la oferta laboral femenina controlando por el estado marital (casada), la edad (eda) y el número de hijos (n_hij) como una variable continua. En la base, e_con toma el valor de 5 para las personas casadas. Genere la variable dummy casada que tome el valor de 1 para las mujeres casadas y cero en otro caso. Estime un modelo de MCO para hrsocup mayor que cero, usando solo la población femenina. Reporte errores robustos. ¿Cuál es la interpretación sobre el coeficiente de los años de escolaridad?
El estimar por MCO, un año más de escolaridad se asocia con 0.17 horas trabajadas más a la semana. Sin embargo, este efecto no es estadísticamente significativo.
data.salarios <- data.salarios %>% mutate(casada=ifelse(e_con==5,1,0)) reg_mco <- lm(hrsocup ~ anios_esc+casada+eda+n_hij, data=filter(data.salarios,hrsocup>0)) coeftest(reg_mco, vcov = vcovHC(reg_mco, "HC1"))[1:4,]Estimate Std. Error t value Pr(>|t|) (Intercept) 36.70129720 1.99116828 18.432042 2.742336e-69 anios_esc 0.17465627 0.10353350 1.686954 9.179628e-02 casada -3.52571327 0.89724706 -3.929479 8.855253e-05 eda 0.06949593 0.04914655 1.414055 1.575295e-01Con modelsummary podemos hacer pedir la tabla de coeficientes. Podemos especificar qué tipo de errores robustos queremos en la opción vcov:
modelsummary(list(reg_mco, reg_mco), vcov = list('classical', 'HC1'), stars = c('***' = 0.01, '**' = 0.05, '*' = 0.1))tinytable_wtd585jx2rmv8jidwuq4 (1) (2) * p < 0.1, ** p < 0.05, *** p < 0.01 (Intercept) 36.701*** 36.701*** (1.959) (1.991) anios_esc 0.175* 0.175* (0.104) (0.104) casada -3.526*** -3.526*** (0.862) (0.897) eda 0.069 0.069 (0.047) (0.049) n_hij -1.149*** -1.149*** (0.336) (0.372) Num.Obs. 1699 1699 R2 0.030 0.030 R2 Adj. 0.028 0.028 AIC 14234.8 14234.8 BIC 14267.4 14267.4 Log.Lik. -7111.383 -7111.383 F 13.171 12.526 RMSE 15.91 15.91 Std.Errors IID HC1 [3 puntos] ¿Qué problema existe con el modelo planteado en el punto anterior en términos de la selección? ¿Considera que se trata de un caso de censura o de truncamiento?
Podemos racionalizar las horas trabajadas en un modelo microeconómico de oferta laboral. Las horas trabajadas observadas son positivas cuando la solución óptima es una cantidad positiva de horas. Sin embargo, si la solución óptima implicara horas negativas, las horas observadas serían cero. En este caso tenemos datos censurados en cero. Si existe una relación positiva entre educación y horas trabajadas, al estimar un modelo por MCO usando solo los datos con horas positivas estamos sobreestimando la media condicional pues se habrán omitido del análisis aquellas mujeres cuya solución a su problema de optimización eran horas iguales a cero o negativas.
[8 puntos] Estime un modelo Tobit de datos censurados. ¿Qué resuelve el modelo Tobit en este caso? Interprete nuevamente el coeficiente sobre los años de escolaridad.
La función tobit permite hacer esto muy fácilmente. Noten que left especifica dónde está la censura. La opción gaussian pone explícito uno de los supuestos críticos del modelo tobit visto en clase: errores normales. Además, se asume homocedasticidad.
reg_tobit <- tobit(hrsocup ~ anios_esc+casada+eda+n_hij, left = 0, right = Inf, dist = "gaussian", data = data.salarios) summary(reg_tobit)Call: tobit(formula = hrsocup ~ anios_esc + casada + eda + n_hij, left = 0, right = Inf, dist = "gaussian", data = data.salarios) Observations: Total Left-censored Uncensored Right-censored 2625 926 1699 0 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.88236 3.19905 0.276 0.78269 anios_esc 0.85530 0.17509 4.885 1.04e-06 *** casada -10.99515 1.43025 -7.688 1.50e-14 *** eda 0.41621 0.07665 5.430 5.64e-08 *** n_hij -1.73840 0.55887 -3.111 0.00187 ** Log(scale) 3.44512 0.01887 182.608 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Scale: 31.35 Gaussian distribution Number of Newton-Raphson Iterations: 3 Log-likelihood: -9086 on 6 Df Wald-statistic: 127.9 on 4 Df, p-value: < 2.22e-16El modelo tobit para datos censurados toma en cuenta que hay una masa de ceros en las horas trabajadas para individuos para los que disponemos de sus características en la base de datos. El modelo tobit ajusta la probabilidad de observar esta masa de ceros. El coeficiente estimado será ahora consistente si el modelo está bien especificado, es decir, si el proceso subyacente es lineal en los parámetros y con un error normal homoscedástico (los supuestos de tobit básico). En este caso, un año más de educación se asocia con 0.85 más horas semanales trabajadas, un efecto estadísticamente significativo. Usar MCO subestimaba el efecto de la escolaridad.
modelsummary acepta la salida de la función tobit:
modelsummary(list(reg_mco, reg_mco, reg_tobit), vcov = list('classical', 'HC1', 'classical'), stars = c('***' = 0.01, '**' = 0.05, '*' = 0.1))tinytable_6pw2w4wa6i7rviiggdof (1) (2) (3) * p < 0.1, ** p < 0.05, *** p < 0.01 (Intercept) 36.701*** 36.701*** 0.882 (1.959) (1.991) (3.199) anios_esc 0.175* 0.175* 0.855*** (0.104) (0.104) (0.175) casada -3.526*** -3.526*** -10.995*** (0.862) (0.897) (1.430) eda 0.069 0.069 0.416*** (0.047) (0.049) (0.077) n_hij -1.149*** -1.149*** -1.738*** (0.336) (0.372) (0.559) Num.Obs. 1699 1699 2625 R2 0.030 0.030 R2 Adj. 0.028 0.028 AIC 14234.8 14234.8 18184.9 BIC 14267.4 14267.4 18220.1 Log.Lik. -7111.383 -7111.383 F 13.171 12.526 RMSE 15.91 15.91 23.02 Std.Errors IID HC1 IID [4 puntos] ¿Cuál es el efecto marginal de un incremento de un año de educación en la oferta laboral? ¿Cómo cambia su respuesta si, en lugar de considerar la variable latente, considera la variable censurada?
El efecto marginal en la variable latente es directamente el coficiente estimado en la parte d., es decir 0.855.
El efecto marginal en la media censurada está dado por:
\[\frac{\partial E(y|x)}{\partial x_j}=\beta_j\Phi(x_i'\beta)\]
Lo que hice aquí fue calcular este efecto marginal para cada individuo y luego obtener el promedio de los efectos marginales en aquellos individuos con horas positivas.
data.salarios <- data.salarios %>% mutate(index1=predict(reg_tobit,.)) %>% mutate(phi=pnorm(index1/reg_tobit$scale)) %>% mutate(mfx_anis_esc=reg_tobit$coefficients[2]*phi, mfx_eda=reg_tobit$coefficients[4]*phi, mfx_n_hij=reg_tobit$coefficients[5]*phi) data.salarios %>% filter(hrsocup>0) %>% summarise(mfx_anis_esc=mean(mfx_anis_esc))# A tibble: 1 × 1 mfx_anis_esc <dbl> 1 0.612
Pregunta 6
Usando los mismos datos del archivo motral2012.csv implementará un ejercicio en el mismo espíritu del famoso estudio de Mroz (1987)1 sobre la oferta laboral femenina. El propósito es estimar la relación entre el salario y el número de horas trabajadas, concentrándonos en la muestra de mujeres.
[5 puntos] El primer problema al que nos enfrentamos es que el salario no se observa para las mujeres que no trabajan. Estime un modelo lineal para el log del salario por hora, ing_x_hrs, usando las variables anios_esc, eda, n_hij, el cuadrado de n_hij, busqueda y casada, usando la submuestra de mujeres con salario por hora positivo. Dichas variables representan los años de escolaridad, la edad, el número de hijos, el cuadrado del número de hijos, si la persona buscó trabajo recientemente y si la persona es casada, respectivamente. Use los coeficientes estimados para imputar el ingreso por hora, faltante para las mujeres que reportan 0 en las horas trabajadas.
Imputamos el salario faltante:
data.salarios<-read_csv("../files/motral2012.csv", locale = locale(encoding = "latin1")) %>% filter(sex==2) %>% mutate(casada=ifelse(e_con==5,1,0)) data.salarios <- data.salarios %>% mutate(ly=ifelse(ing_x_hrs>0,log(ing_x_hrs),NA)) reg_imput <- lm(ly ~ anios_esc+casada+eda+n_hij+n_hij^2+busqueda, data = data.salarios) data.salarios <- data.salarios %>% mutate(lyhat = predict(reg_imput, .)) %>% mutate(ly = ifelse(is.na(ly), lyhat, ly))Aquí tomé en cuenta que hay personas con horas trabajadas positivas e ingreso cero. En ese caso puse un NA al log del salario. Luego, en la imputación, le asigné el valor ajustado a estas observaciones junto con todas las que tienen el log del salario faltante.
[5 puntos] Use heckit de la librería sampleSelection para estimar por máxima verosimilitud un heckit para las horas trabajadas hrsocup. En la ecuación de selección (si la persona trabaja o no) incluya como variable explicativa el salario por hora (imputado para las mujeres que no trabajan), además de anios_esc, eda, n_hij, el cuadrado de n_hij, casada y busqueda (esta última es un indicador de si se buscó trabajo en la última semana). En la ecuación de horas, incluya los mismos regresores, excepto n_hij, su cuadrado y busqueda.
La función heckit permite estimar el modelo de Heckman por máxima verosimilitud de manera muy simple. Hay que especificar method=“ml” para que la estimación sea por máxima verosimilitud:
data.salarios <- data.salarios %>% mutate(trabaja = ifelse(hrsocup>0,1,0)) %>% mutate(trabaja = factor(trabaja,levels=c(0,1))) reg_heckit_mv <- heckit(trabaja ~ anios_esc+casada+eda+ly+n_hij+n_hij^2+busqueda, hrsocup ~ anios_esc+casada+eda+ly, method="ml", data = data.salarios) summary(reg_heckit_mv)-------------------------------------------- Tobit 2 model (sample selection model) Maximum Likelihood estimation Newton-Raphson maximisation, 3 iterations Return code 8: successive function values within relative tolerance limit (reltol) Log-Likelihood: -7181.675 2625 observations (926 censored and 1699 observed) 14 free parameters (df = 2611) Probit selection equation: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.583614 0.448320 -5.763 9.24e-09 *** anios_esc 0.005346 0.020341 0.263 0.793 casada -0.213125 0.145135 -1.468 0.142 eda -0.003391 0.008137 -0.417 0.677 ly -0.004236 0.133344 -0.032 0.975 n_hij 0.023985 0.058900 0.407 0.684 busqueda 2.406669 0.104595 23.009 < 2e-16 *** Outcome equation: Estimate Std. Error t value Pr(>|t|) (Intercept) 55.62469 2.17656 25.556 < 2e-16 *** anios_esc 1.04819 0.09995 10.487 < 2e-16 *** casada -3.58856 0.77967 -4.603 4.37e-06 *** eda 0.11614 0.03902 2.977 0.00294 ** ly -9.83418 0.60389 -16.285 < 2e-16 *** Error terms: Estimate Std. Error t value Pr(>|t|) sigma 14.8579 0.2591 57.350 <2e-16 *** rho -0.1606 0.1964 -0.818 0.414 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 --------------------------------------------Podemos reportar con modelsummary, aunque realmente lo hace muy mal:
modelsummary(list(reg_heckit_mv), stars = c('***' = 0.01, '**' = 0.05, '*' = 0.1))tinytable_vsm7u8tt6ve87piqbswr (1) * p < 0.1, ** p < 0.05, *** p < 0.01 (Intercept) -2.584*** 55.625*** (0.448) (2.177) anios_esc 0.005 1.048*** (0.020) (0.100) casada -0.213 -3.589*** (0.145) (0.780) eda -0.003 0.116*** (0.008) (0.039) ly -0.004 -9.834*** (0.133) (0.604) n_hij 0.024 (0.059) busqueda 2.407*** (0.105) sigma 14.858*** (0.259) rho -0.161 (0.196) Num.Obs. 2625 AIC 14391.4 BIC 14473.6 RMSE 14.84 [10 puntos] Estime ahora el heckit en dos pasos, a mano. Es decir, siga los siguientes pasos: i) estime un probit para la ecuación de selección y obtenga el índice \(x_i'\hat{\beta}\); ii) calcule el inverso de la razón de Mills \(\lambda_i(x_i'\hat{\beta})\); y iii) estime por MCO la ecuación para las horas trabajadas con la submuestra que tiene horas trabajadas positivas, incluyendo como regresor el inverso de la razón de Mills estimado y el resto de los regresores. Compare los coeficientes y los errores estándar obtenidos en esta parte con los de la parte b. ¿Por qué son iguales o por qué difieren?
Estimamos ahora el heckit a mano. Estimamos el probit y obtenemos el valor ajustado del IMR:
reg_heckit_pe <- glm(trabaja ~ anios_esc+casada+eda+ly+n_hij+n_hij^2+busqueda, family = binomial(link = "probit"), data = data.salarios) data.salarios <- data.salarios %>% mutate(index = predict(reg_heckit_pe, .)) %>% mutate(imr = dnorm(index)/pnorm(index))reg_heckit_se <- lm(hrsocup ~ anios_esc+casada+eda+ly+imr, data=filter(data.salarios,trabaja==1)) summary(reg_heckit_se)Call: lm(formula = hrsocup ~ anios_esc + casada + eda + ly + imr, data = filter(data.salarios, trabaja == 1)) Residuals: Min 1Q Median 3Q Max -46.172 -10.085 1.915 9.253 57.689 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 55.68676 2.17948 25.550 < 2e-16 *** anios_esc 1.04814 0.10000 10.481 < 2e-16 *** casada -3.56927 0.78057 -4.573 5.17e-06 *** eda 0.11621 0.03904 2.977 0.00296 ** ly -9.82971 0.60406 -16.273 < 2e-16 *** imr -3.94669 3.62684 -1.088 0.27667 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 14.86 on 1693 degrees of freedom Multiple R-squared: 0.1563, Adjusted R-squared: 0.1539 F-statistic: 62.75 on 5 and 1693 DF, p-value: < 2.2e-16Comparamos coeficientes (aquí stargazer lo hace mejor):
stargazer(reg_heckit_mv, reg_heckit_se, type="html", df=FALSE, digits=4)Dependent variable:
hrsocup
Heckman
OLS
selection
(1)
(2)
anios_esc
1.0482***
1.0481***
(0.0999)
(0.1000)
casada
-3.5886***
-3.5693***
(0.7797)
(0.7806)
eda
0.1161***
0.1162***
(0.0390)
(0.0390)
ly
-9.8342***
-9.8297***
(0.6039)
(0.6041)
imr
-3.9467
(3.6268)
Constant
55.6247***
55.6868***
(2.1766)
(2.1795)
Observations
2,625
1,699
R2
0.1563
Adjusted R2
0.1539
Residual Std. Error
14.8614
F Statistic
62.7490***
Note:
p<0.1; p<0.05; p<0.01
La magnitud de los coeficientes es práctiamente la misma entre el modelo estimado por máxima verosimilitud y con un procedimiento en dos etapas a mano. En este ejemplo las diferencias son sutiles, aunque recordemos que en general la estimación por MV es más eficiente si la verosimilitud está bien planteada.
Notas
Mroz, T. A. (1987). The sensitivity of an empirical model of married women’s hours of work to economic and statistical assumptions. Econometrica: Journal of the econometric society, 765-799.↩︎