ECNII Tarea 1
Rafael Martínez Martínez
9/13/2020
1 Problema 6
## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## ── Attaching packages ─────────────────────────────── tidyverse 1.3.0 ──## ✓ tibble 3.0.1 ✓ dplyr 1.0.0
## ✓ tidyr 1.1.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ✓ purrr 0.3.4## ── Conflicts ────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## some## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: betareg##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazer1.1 Inciso a)
Genere una muestra de \(10,000\) observaciones llamadas \(x\) tales que \(x\sim\mathcal{N}(1,1)\). Posteriormente, genere \(\lambda=exp(\beta_1+\beta_2x)\), donde \((\beta_1\;\beta_2)=(2\;-1)\). Note que \(1/\lambda\) es conocida como la tasa en la distribución exponencial. En R, requiere especificar como parámetro a la tasa en lugar de \(\lambda\)
1.2 Inciso b)
Reporte una tabla con la media, la desviación estándar, el mínimo y el máximo de \(x\), \(\lambda\) y \(y\)
## [1] "media"## x lambda y
## 0.9983458 4.5012090 0.6098114## [1] "desviación estandar"## x lambda y
## 1.004260 6.144595 1.280866## [1] "mínimo"## x lambda y
## -3.170064e+00 3.973818e-02 2.606608e-06## [1] "máximo"## x lambda y
## 5.225443 175.926037 37.9054741.3 Inciso c)
Reporte una gráfica donde muestre la relación entre \(x\) y \(\lambda\) en el plano \((x,\lambda)\). Realice otra gráfica similar, ahora para \((x,1/\lambda)\).
1.4 Inciso d)
Estime por MCO una regresión entre \(y\) y \(x\). Deberá obtener coeficientes parecidos a los de la primera columna de la Tabla 5.7 en CT.
##
## Call:
## lm(formula = y ~ x, data = var)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.995 -0.481 -0.116 0.250 35.853
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.005149 0.015816 -0.326 0.745
## x 0.615979 0.011169 55.149 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.122 on 9998 degrees of freedom
## Multiple R-squared: 0.2332, Adjusted R-squared: 0.2332
## F-statistic: 3041 on 1 and 9998 DF, p-value: < 2.2e-161.5 Inciso e)
¿Por qué difieren los coeficientes que obtuvo y los que se presentan en la Tabla 5.7 de CT?
Hoy diferentes motivos, uno es que las muestras de \(x\) son aleatorias, y con esto se genera \(\lambda\). Otro factor podría ser la forma en como se generan las muestras con distribución exponencial \(y\), podría ser que el método que utiliza CT en Stata no coincida con el de rexp en R
1.6 Inciso f)
Obtenga los errores robustos. En R, una librería que será muy útil es sandwich.
Se configura para obtener los errorres robustos tipo Stata
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0051487 0.0173073 -0.2975 0.7661
## x 0.6159791 0.0257321 23.9381 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Las estimaciones son 0.0173073 y 0.0257321
1.7 Inciso g)
¿El estimador de MCO es consistente? ¿Por qué?
No es consistente pues para que esto suceda \(\hat\beta_{mc}\xrightarrow{p}\beta_0\), lo cual no podría ser pues ajusta una recta a \(y\) cuando la distribución es exponencial
1.8 Inciso h)
Plantee la función de log verosimilitud.
La función de densidad para \(y\)
\[f(y)=\lambda e^{-\lambda y}\] Se tiene que
\[\lambda=e^{\beta_1+\beta_2 x}\] Así para cada \(y_i\)
\[f(y_i)=e^{\beta_1+\beta_2 x_i}e^{-e^{(\beta_1+\beta_2 x_i)}y_i}\]
Entonces, suponiendo muestras iid y aplicando ln
\[\mathcal{L}(\beta)=\sum_i \left(\beta_1+\beta_2 x_i - e^{(\beta_1+\beta_2 x_i)}y_i\right)\]
1.9 Inciso i)
Obtenga el estimador de máxima verosimilitud de \(\beta_1\) y \(\beta_2\) obteniendo la solución al negativo del problema de log verosimilitud. En R, puede emplear, por ejemplo nlm.
El estimador de máxima verosimilitud resuelve el problema
\[\hat\beta=\max_\beta{\frac{1}{N}}\mathcal{L}(\beta)=\max_\beta Q_N(\beta)\] Así
\[\frac{\partial Q_N}{\partial\beta_1}=\frac{1}{N}\sum_i\left(1-e^{(\beta_1+\beta_2 x_i)}y_i\right)=0\] \[\frac{\partial Q_N}{\partial\beta_2}=\frac{1}{N}\sum_i\left(x_i-e^{(\beta_1+\beta_2 x_i)}x_iy_i\right)=0\] Resolvemos este problema no lineal con nlm, pero ajutamos el signo, pues la herramienta de R resuelve minimizaciónes, no hacemos el escalamiento de \(1/N\) pues la matriz de varianza robusta está en términos de la matriz hesiana de \(\mathcal{L}\)
f <- function(Beta) {
beta1 <- Beta[1]
beta2 <- Beta[2]
-sum(beta1 + beta2 * x - exp(beta1 + beta2 * x) * y)
}
estnum <- nlm(f, c(2,-1), hessian = TRUE)
estnum## $minimum
## [1] 18.32851
##
## $estimate
## [1] 1.9884504 -0.9919251
##
## $gradient
## [1] 9.290707e-08 1.776357e-08
##
## $hessian
## [,1] [,2]
## [1,] 10001.979 9985.449
## [2,] 9985.449 20162.871
##
## $code
## [1] 1
##
## $iterations
## [1] 4Los coeficientes estimados son 1.9884504, -0.9919251
1.10 Inciso j)
Usando la matriz hesiana obtenida en la solución del problema de optimización, encuentre los errores estándar robustos de los coeficientes estimados.
#El menos de la matriz hessiana es porque se minimizo y el menos fuera de la #inversa es porque la matriz de varianzas y covarianzas es el negativo, al
# al final los signos se cancelan
A<- -solve(-estnum$hessian)
diag(A)^(1/2)## [1] 0.014062499 0.0099044251.11 Inciso k)
El modelo antes descrito puede expresarse como una regresión no lineal de la forma \(y=exp(-x'\beta)+u\). Encuentre la solución para \(\beta_1\) y \(\beta_2\). Reporte los errores estándar no robustos. ¿Son consistentes estos errores? ¿Por qué?
##
## Formula: y ~ exp(-(beta1 + beta2 * x))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## beta1 2.091460 0.029598 70.66 <2e-16 ***
## beta2 -1.037482 0.008995 -115.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.967 on 9998 degrees of freedom
##
## Number of iterations to convergence: 2
## Achieved convergence tolerance: 3.93e-06Los errores estandar no robustos \(0.02959\) y \(0.0089\) no son consistentes, pues los datos se optubieron emdiante el muestreo aleatorio y en este la varianza no fue constante
1.12 Inciso l)
Ahora implementará la matriz de varianzas y covarianzas robusta en la ecuación 5.81 de CT. Dé una expresión para \(\hat{D}\) en este problema
\[\hat D=\frac{\partial e^{-x'\beta}}{\partial \beta}\mid_{\hat\beta}=\begin{bmatrix}-e^{-\beta_1-\beta_2x}& -e^{-\beta_1-\beta_2x}x\end{bmatrix}\mid_{\hat\beta}\]
1.13 Inciso m)
Calcule el error estándar robusto definido como en la ecuación 5.81. En este caso \(\hat{\Omega}=Diad(\hat{u}_i^2)\)
uhat=y-exp(-beta1estnls-beta2estnls*x)
uhat2= uhat^2
omega= diag(uhat2)
#corección práctica N/(N-2)
Vr = (10000/9998)*solve(t(D)%*%D)%*%t(D)%*%omega%*%D%*%solve(t(D)%*%D)
diag(Vr)^(1/2)## [1] 0.09884300 0.040986721.14 Inciso n)
Calcule una versión alternativa de errores estándar (entre corchetes en Tabla 5.7), esta vez con \(\hat{\Omega}=Diag((exp(-x_i'\beta))^2)\)
uotro=(exp(-beta1estnls-beta2estnls*x))^2
omega2= diag(uotro)
#corección práctica N/(N-2)
Vr2 = (10000/9998)*solve(t(D)%*%D)%*%t(D)%*%omega2%*%D%*%solve(t(D)%*%D)
diag(Vr2)^(1/2)## [1] 0.3857799 0.14923141.15 Inciso o)
En este experimento, ¿qué estimador tiene las mejores propiedades?
EL MLE es mejor que el NLS y que el OLS
2 Problema 7
2.1 a.
plin <- lm(arr86 ~ pcnv + avgsen + tottime +
ptime86 + inc86 + black + hispan + born60,
data=data)
summary(plin)##
## Call:
## lm(formula = arr86 ~ pcnv + avgsen + tottime + ptime86 + inc86 +
## black + hispan + born60, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5665 -0.3152 -0.1942 0.5148 1.1544
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.360983 0.016093 22.432 < 2e-16 ***
## pcnv -0.154380 0.020934 -7.375 2.18e-13 ***
## avgsen 0.003502 0.006342 0.552 0.581
## tottime -0.002061 0.004888 -0.422 0.673
## ptime86 -0.021595 0.004468 -4.833 1.42e-06 ***
## inc86 -0.001225 0.000127 -9.648 < 2e-16 ***
## black 0.161718 0.023504 6.880 7.38e-12 ***
## hispan 0.089259 0.020559 4.342 1.47e-05 ***
## born60 0.002870 0.017199 0.167 0.867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4294 on 2716 degrees of freedom
## Multiple R-squared: 0.08239, Adjusted R-squared: 0.07969
## F-statistic: 30.48 on 8 and 2716 DF, p-value: < 2.2e-16##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.36098314 0.01670805 21.6053 < 2.2e-16 ***
## pcnv -0.15438020 0.01896403 -8.1407 5.921e-16 ***
## avgsen 0.00350240 0.00588763 0.5949 0.5520
## tottime -0.00206130 0.00422557 -0.4878 0.6257
## ptime86 -0.02159526 0.00275318 -7.8437 6.233e-15 ***
## inc86 -0.00122484 0.00011414 -10.7310 < 2.2e-16 ***
## black 0.16171828 0.02552788 6.3350 2.769e-10 ***
## hispan 0.08925863 0.02106893 4.2365 2.346e-05 ***
## born60 0.00286982 0.01715962 0.1672 0.8672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.2 Inciso a)
¿Cuál es el efecto en la probabilidad de arresto si pcnv pasa de \(0.25\) a \(0.75\)?
## pcnv
## -0.07719012.3 Inciso b)
Realice una prueba de significancia conjunta de avgsen y tottime. ¿Qué concluye?
## Linear hypothesis test
##
## Hypothesis:
## avgsen = 0
## tottime = 0
##
## Model 1: restricted model
## Model 2: arr86 ~ pcnv + avgsen + tottime + ptime86 + inc86 + black + hispan +
## born60
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2718 500.91
## 2 2716 500.84 2 0.06606 0.1791 0.836No hay evidencia estadística para rechazar que ambos coeficentes sean cero, es decir, la variable de salida no esta explicada por estas variables, pensando en una relación lineal.
2.4 Inciso c)
Estime un modelo probit relacionando las mismas variables. ¿Cuál es el efecto en la probabilidad de arresto cuando pcnv pasa de \(0.25\) a \(0.75\) para los valores promedio de avgsen, tottime, inc86 y ptime86 y cuando los individuos son de raza negra, no hispánicos y nacidos en 1960 (born60 igual a 1). ¿Cómo se comparan estos resultados con lo que encontró con el modelo de probabilidad lineal?
probit <- glm( arr86 ~ pcnv + avgsen + tottime +
ptime86 + inc86 + black + hispan + born60,
family = binomial(link = "probit"),
data = data)
coeftest(probit, vcov. = vcovHC, type = "HC1")##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.31383330 0.05192016 -6.0445 1.498e-09 ***
## pcnv -0.55292513 0.06892627 -8.0220 1.041e-15 ***
## avgsen 0.01273947 0.01903441 0.6693 0.5033
## tottime -0.00764851 0.01469718 -0.5204 0.6028
## ptime86 -0.08119959 0.01218345 -6.6647 2.651e-11 ***
## inc86 -0.00463461 0.00053825 -8.6104 < 2.2e-16 ***
## black 0.46660774 0.07138053 6.5369 6.280e-11 ***
## hispan 0.29110055 0.06540227 4.4509 8.550e-06 ***
## born60 0.01120673 0.05592090 0.2004 0.8412
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1efecto <- coef(probit)%*%c(1,0.75,mean(data$avgsen), mean(data$tottime), mean(data$inc86), mean(data$ptime86),1,0,1)-coef(probit)%*%c(1,0.25,mean(data$avgsen), mean(data$tottime), mean(data$inc86), mean(data$ptime86),1,0,1)
efecto## [,1]
## [1,] -0.2764626El cambio calculado es cuatro veces más grande, en la misma dirección, pero en ambos casos no hay evidencia estadística para rechazar que la variable de interés no dependa de estas variables.
3 Problema 8
3.1 Inciso a)
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. Se sugiere que posteriormente etiquete los valores de ahorro 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”.
## Parsed with column specification:
## cols(
## .default = col_double(),
## nom = col_character(),
## obs_caratu = col_character(),
## p12_des = col_character(),
## p14_des = col_character(),
## p16_des = col_character(),
## p19_1_des = col_character(),
## p19_2_des = col_character(),
## p23_des = col_character(),
## obs = col_character(),
## cs_ad_mot = col_logical(),
## cs_p20_des = col_logical(),
## cs_ad_des = col_logical(),
## cs_nr_mot = col_logical(),
## cs_p22_des = col_logical(),
## cs_nr_ori = col_logical()
## )## See spec(...) for full column specifications.3.2 Inciso b)
Estime un modelo logit multinomial (regresores invariantes a la alternativa) con la opción de ahorro como variable dependiente y con la edad (eda) y la condición como jefe de hogar jefe_hog como variables independientes. ¿Qué puede decir sobre el coeficiente de edad en la alternativa “Casa”?
## # weights: 12 (6 variable)
## initial value 2727.854313
## iter 10 value 2606.707110
## iter 10 value 2606.707107
## iter 10 value 2606.707107
## final value 2606.707107
## converged##
## ==============================================
## Dependent variable:
## ----------------------------
## Casa Otro
## (1) (2)
## ----------------------------------------------
## eda -0.042*** -0.003
## (0.005) (0.005)
##
## jefe_hog -0.352*** 0.008
## (0.114) (0.110)
##
## Constant 1.104*** -0.423**
## (0.172) (0.180)
##
## ----------------------------------------------
## Akaike Inf. Crit. 5,225.414 5,225.414
## ==============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Un incremento de la edad genera un descenso de la probabilidad del ahorro para casa respecto al ahorro en banco.
3.3 Inciso c)
Calcule los efectos marginales sobre la probabilidad de ahorrar en el banco. Al considerar el cambio de no ser jefe de hogar a serlo, ¿de qué tamaño es el efecto predicho en la probabilidad de ahorrar en el banco?
## factor AME SE z p lower upper
## eda 0.0056 NA NA NA NA NA
## jefe_hog 0.0434 NA NA NA NA NAEn promedio la probabilidad de ahorar en banco es 4.4% más alta cuando se trata de jefes de hogar que cundo no, pero el nivel de edad es el mismo
3.4 Inciso d)
Calcule los cocientes de riesgo relativo (relative risk ratios, RRR). ¿Qué significa el hecho de que el RRR de jefe de hogar sea mayor que 1 en la alternativa “Otro”?
## (Intercept) eda jefe_hog
## Casa 3.0167164 0.9591874 0.7031443
## Otro 0.6548042 0.9971362 1.0081656La tasa de riesgo relativo de ahorrar en la opción “otro” si se es jefe de hogar es 1.008 respecto a no ser jefe de hogar.
3.5 Inciso e)
Estime nuevamente el modelo, pero ahora, especifique que la alternativa “Casa” sea la alternativa base. ¿Cómo es el coeficiente de la edad en la alternativa “Banco”? ¿Es esto congruente con lo que noto en la parte d. de este problema?
## # weights: 12 (6 variable)
## initial value 2727.854313
## iter 10 value 2606.707293
## final value 2606.707107
## converged##
## =========================================================
## Dependent variable:
## ---------------------------------------
## Banco Otro Casa Otro
## (1) (2) (3) (4)
## ---------------------------------------------------------
## eda 0.042*** 0.039*** -0.042*** -0.003
## (0.005) (0.006) (0.005) (0.005)
##
## jefe_hog 0.352*** 0.360*** -0.352*** 0.008
## (0.114) (0.127) (0.114) (0.110)
##
## Constant -1.104*** -1.528*** 1.104*** -0.423**
## (0.172) (0.194) (0.172) (0.180)
##
## ---------------------------------------------------------
## Akaike Inf. Crit. 5,225.414 5,225.414 5,225.414 5,225.414
## =========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01## (Intercept) eda jefe_hog
## Banco 0.3314884 1.042549 1.422183
## Otro 0.2170632 1.039563 1.433776Hay un cambio de signo en edad. Pues un aumento de edad en una unidad afecta negativamente la tasa del riesgo relativo en ahorro en casa.
4 Problema 9
## Parsed with column specification:
## cols(
## caseid = col_double(),
## art = col_double(),
## female = col_character(),
## married = col_character(),
## kid5 = col_double(),
## phd = col_double(),
## mentor = col_double(),
## arttrunc = col_double(),
## kid5_0 = col_character(),
## kid5_1 = col_character(),
## kid5_2plus = col_character(),
## kid5cat = col_character(),
## male = col_character(),
## phdadeq = col_character(),
## phdcat = col_character(),
## phdelite = col_character(),
## phdgood = col_character(),
## phdstrong = col_character(),
## profage = col_double()
## )4.1 Inciso a)
¿Hay evidencia de sobredispersión en la variable ?
## [1] 3.709742## [1] 1.692896Observamos que hay evidencia de sobredispersión pues la varianza es más del doble que la media.
4.2 Inciso b)
Independientemente de si hay sobredispersión o no, estime un modelo Poisson que incluya variables dicotómicas para estudiantes mujeres y para estudiantes casadas o casados, la cantidad de hijos mejores de cinco años, el ranking de prestigio del doctorado () y el número de artículos publicados por su mentor. Interprete los coeficientes estimados.
data <- datos %>%
mutate(mujeres=ifelse(female=="Female",1,0))%>%
mutate(casados=ifelse(married=="Married",1,0))
poissa <- glm(art ~ mujeres + casados + kid5 + phd + mentor, family=poisson,
data=data)
summary(poissa)##
## Call:
## glm(formula = art ~ mujeres + casados + kid5 + phd + mentor,
## family = poisson, data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.5672 -1.5398 -0.3660 0.5722 5.4467
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.304617 0.102981 2.958 0.0031 **
## mujeres -0.224594 0.054613 -4.112 3.92e-05 ***
## casados 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: 5##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3046168 0.1470022 2.0722 0.038248 *
## mujeres -0.2245942 0.0718983 -3.1238 0.001785 **
## casados 0.1552434 0.0821991 1.8886 0.058942 .
## kid5 -0.1848827 0.0561477 -3.2928 0.000992 ***
## phd 0.0128226 0.0421024 0.3046 0.760704
## mentor 0.0255427 0.0038303 6.6685 2.584e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El efecto depende de \(\beta_jexp(x_i'\beta)\) y sera diferente para cada individuo, pero como \(exp(x_i'\beta)>0\) el signo de cada coeficente nos indicará si el efecto positivo o negativo si se se aumenta en una unidad \(x_j\), por ejemplo el echo de ser mujer o tener hijos menores de 5 años afecta de manera negativa la cantidad de articulos, y el echo de estar casado, tener un mentor porductivo y estar en un posgrado de exigencia, aumentara el numero de artículos
4.3 Inciso c)
Obtenga la razón de tasas de incidencia (IRR) para los coeficientes e interprete los resultados.
## Call:
## poissonirr(formula = art ~ mujeres + casados + kid5 + phd + mentor,
## data = data)
##
## Incidence-Rate Ratio:
## IRR Std. Err. z P>|z|
## mujeres 0.798840 0.043627 -4.1124 3.915e-05 ***
## casados 1.167942 0.071682 2.5294 0.01142 *
## kid5 0.831202 0.033354 -4.6075 4.076e-06 ***
## phd 1.012905 0.026738 0.4858 0.62714
## mentor 1.025872 0.002058 12.7327 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Las tasas de incidencias para coeficientes positivos (previamente obtenidos) son mayores a 1; mientras que als tasas de incidencia para coefientes negativos son menores a 1.
4.4 Inciso d)
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 en R.
poissd <- glm(art ~ mujeres + casados + kid5 + phd + mentor,
offset = profage,
family = poisson,
data = data )
summary(poissd)##
## Call:
## glm(formula = art ~ mujeres + casados + kid5 + phd + mentor,
## family = poisson, data = data, offset = profage)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -23.7550 -0.1226 2.2364 4.0564 17.9577
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -35.835084 0.092001 -389.506 < 2e-16 ***
## mujeres 15.259449 0.058453 261.056 < 2e-16 ***
## casados -0.414288 0.065193 -6.355 2.09e-10 ***
## kid5 -0.318704 0.040167 -7.935 2.11e-15 ***
## phd 0.666198 0.026137 25.489 < 2e-16 ***
## mentor 0.051771 0.002517 20.572 < 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: 33647 on 914 degrees of freedom
## Residual deviance: 14026 on 909 degrees of freedom
## AIC: 15706
##
## Number of Fisher Scoring iterations: 11## Call:
## poissonirr(formula = poissd, data = data)
##
## Incidence-Rate Ratio:
## IRR Std. Err. z P>|z|
## mujeres 4.2374e+06 2.4768e+05 261.0564 < 2.2e-16 ***
## casados 6.6081e-01 4.3080e-02 -6.3548 2.086e-10 ***
## kid5 7.2709e-01 2.9205e-02 -7.9345 2.113e-15 ***
## phd 1.9468e+00 5.0884e-02 25.4889 < 2.2e-16 ***
## mentor 1.0531e+00 2.6504e-03 20.5717 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Si se toma en cuenta que las mujeres estuvieron menos expuestas al evento publicar; las mujeres publicaron más articulos que los hombres.
El efecto marginal de ser mujer en el numero de articulos publicados es mucho mayor respecto a ser hombre. Hay un cambio de signo, estar casado ahora afecta de manera negativa. En las demás variables la dirección es la misma con mayor magnitud.
4.5 Inciso e)
Emplee ahora un modelo negativo binomial con sobredispersión constante 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.
bne <- glm(art ~ mujeres + casados + kid5 + phd + mentor,
family = quasipoisson,
data = data)
stargazer(bne, poissa, type= "text")##
## ==============================================
## Dependent variable:
## ----------------------------
## art
## glm: quasipoisson Poisson
## link = log
## (1) (2)
## ----------------------------------------------
## mujeres -0.225*** -0.225***
## (0.074) (0.055)
##
## casados 0.155* 0.155**
## (0.083) (0.061)
##
## kid5 -0.185*** -0.185***
## (0.054) (0.040)
##
## phd 0.013 0.013
## (0.036) (0.026)
##
## mentor 0.026*** 0.026***
## (0.003) (0.002)
##
## Constant 0.305** 0.305***
## (0.139) (0.103)
##
## ----------------------------------------------
## Observations 915 915
## Log Likelihood -1,651.056
## Akaike Inf. Crit. 3,314.113
## ==============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Nuevamente se tiene la misma interpretación que el inciso b), los coeficientes incluso son los mismos que en el inciso mencionado. Los errores estándar ahora son más grandes y el grado de significancia se mantuvo en el mismo nivel.
4.6 Inciso f)
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?
##
## Call:
## glm.nb(formula = art ~ mujeres + casados + kid5 + phd + mentor,
## data = data, init.theta = 2.264387695, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1678 -1.3617 -0.2806 0.4476 3.4524
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.256144 0.137348 1.865 0.062191 .
## mujeres -0.216418 0.072636 -2.979 0.002887 **
## casados 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.917##
## ===============================================================
## Dependent variable:
## ---------------------------------------------
## art
## glm: quasipoisson negative Poisson
## link = log binomial
## (1) (2) (3)
## ---------------------------------------------------------------
## mujeres -0.225*** -0.216*** -0.225***
## (0.074) (0.073) (0.055)
##
## casados 0.155* 0.150* 0.155**
## (0.083) (0.082) (0.061)
##
## kid5 -0.185*** -0.176*** -0.185***
## (0.054) (0.053) (0.040)
##
## phd 0.013 0.015 0.013
## (0.036) (0.036) (0.026)
##
## mentor 0.026*** 0.029*** 0.026***
## (0.003) (0.003) (0.002)
##
## Constant 0.305** 0.256* 0.305***
## (0.139) (0.137) (0.103)
##
## ---------------------------------------------------------------
## Observations 915 915 915
## Log Likelihood -1,561.958 -1,651.056
## theta 2.264*** (0.271)
## Akaike Inf. Crit. 3,135.917 3,314.113
## ===============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01La interpretación se mantiene ahora los errores son un poco menores a los del sobredispersión constante pero siguen muy alejados de lo errores de poisson.
La estimación del parámetro es 2.26 con el 1% de significancia, así se rechaza la hipótesis nula \(H_0: \alpha=0\) a favor de utilizar un modelo binomial negativo.