1 Problema 6

library(ggplot2)
library(estimatr)
library(sandwich)
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(tidyverse)

## ── 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()

library(dplyr)
library(car)

## 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

library(mfx)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: betareg

library(stargazer)

## 
## 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=stargazer

1.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\)

set.seed(2003)
N <- 10000
mean <- 1
sd <- 1
beta <- c(2,-1)
# numeros con distribución normal
x <- rnorm(N,mean,sd)
# Calculo de lambda con los valores anteriores
lambda <- exp(beta[1]+beta[2]*x)
# generamos un numero aleatorio con ditribucion
y <- (1/lambda)*rexp(1/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\)

var <- data.frame(x=x, lambda=lambda,
                  y=y)
print("media")

## [1] "media"

sapply(var,mean)

##         x    lambda         y 
## 0.9983458 4.5012090 0.6098114

print("desviación estandar")

## [1] "desviación estandar"

sapply(var,sd)

##        x   lambda        y 
## 1.004260 6.144595 1.280866

print("mínimo")

## [1] "mínimo"

sapply(var,min)

##             x        lambda             y 
## -3.170064e+00  3.973818e-02  2.606608e-06

print("máximo")

## [1] "máximo"

sapply(var,max)

##          x     lambda          y 
##   5.225443 175.926037  37.905474

1.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)\).

ggplot(data=var, aes(x=x, y=lambda, group=1)) +
  geom_line(color="red")+
  geom_point()

ggplot(data=var, aes(x=x, y=1/lambda, group=1)) +
  geom_line(color="green")+
  geom_point()

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.

mco <- lm(y ~ x, data = var)
summary(mco)

## 
## 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-16

1.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

coeftest(mco, vcov = vcovHC(mco, "HC1"))

## 
## 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 ' ' 1

Las 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] 4

Los 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.009904425

1.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é?

estnls<-nls(y~exp(-(beta1+beta2*x)),start=list(beta1=2,beta2=-1), data=var)
summary(estnls)

## 
## 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-06

Los 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}\]

beta1estnls = coef(estnls)[1]
beta2estnls = coef(estnls)[2]


D = cbind(-exp(-beta1estnls-beta2estnls*x),-exp(-beta1estnls-beta2estnls*x)*x)

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.04098672

1.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.1492314

1.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 9

datos <- read_csv("phd_articulos.csv")

## 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()
## )

datos$married

##   [1] "Married" "Married" "Married" "Married" "Single"  "Married" "Married"
##   [8] "Married" "Single"  "Married" "Married" "Single"  "Married" "Married"
##  [15] "Married" "Married" "Married" "Married" "Married" "Married" "Single" 
##  [22] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
##  [29] "Married" "Single"  "Single"  "Married" "Single"  "Married" "Married"
##  [36] "Married" "Married" "Single"  "Married" "Married" "Married" "Married"
##  [43] "Single"  "Married" "Married" "Married" "Married" "Married" "Married"
##  [50] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
##  [57] "Single"  "Single"  "Single"  "Married" "Single"  "Married" "Married"
##  [64] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
##  [71] "Married" "Single"  "Married" "Married" "Married" "Married" "Married"
##  [78] "Married" "Single"  "Married" "Married" "Married" "Married" "Single" 
##  [85] "Married" "Married" "Single"  "Married" "Married" "Married" "Married"
##  [92] "Married" "Married" "Single"  "Married" "Single"  "Married" "Married"
##  [99] "Married" "Married" "Married" "Married" "Married" "Married" "Single" 
## [106] "Married" "Married" "Married" "Single"  "Single"  "Married" "Married"
## [113] "Married" "Married" "Single"  "Single"  "Married" "Married" "Single" 
## [120] "Married" "Single"  "Single"  "Married" "Married" "Married" "Married"
## [127] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
## [134] "Married" "Married" "Married" "Married" "Married" "Single"  "Married"
## [141] "Married" "Married" "Single"  "Married" "Single"  "Married" "Married"
## [148] "Married" "Single"  "Married" "Married" "Married" "Married" "Married"
## [155] "Single"  "Married" "Married" "Married" "Married" "Married" "Married"
## [162] "Married" "Married" "Married" "Married" "Single"  "Married" "Married"
## [169] "Married" "Single"  "Married" "Married" "Married" "Married" "Single" 
## [176] "Single"  "Married" "Married" "Married" "Married" "Married" "Married"
## [183] "Married" "Married" "Married" "Married" "Married" "Married" "Single" 
## [190] "Married" "Married" "Married" "Single"  "Married" "Married" "Married"
## [197] "Married" "Married" "Married" "Married" "Single"  "Single"  "Married"
## [204] "Married" "Single"  "Single"  "Married" "Married" "Married" "Married"
## [211] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
## [218] "Married" "Single"  "Single"  "Married" "Married" "Married" "Single" 
## [225] "Married" "Single"  "Married" "Married" "Married" "Single"  "Married"
## [232] "Married" "Married" "Married" "Single"  "Married" "Married" "Married"
## [239] "Married" "Single"  "Married" "Married" "Married" "Single"  "Single" 
## [246] "Married" "Married" "Single"  "Married" "Married" "Married" "Married"
## [253] "Married" "Married" "Single"  "Married" "Married" "Married" "Married"
## [260] "Married" "Single"  "Married" "Married" "Married" "Married" "Married"
## [267] "Married" "Married" "Married" "Married" "Single"  "Married" "Married"
## [274] "Married" "Married" "Married" "Married" "Married" "Single"  "Married"
## [281] "Married" "Single"  "Married" "Single"  "Married" "Single"  "Single" 
## [288] "Married" "Married" "Married" "Single"  "Married" "Married" "Married"
## [295] "Married" "Married" "Married" "Married" "Married" "Single"  "Married"
## [302] "Married" "Single"  "Married" "Married" "Single"  "Married" "Married"
## [309] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
## [316] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
## [323] "Single"  "Married" "Single"  "Married" "Married" "Married" "Married"
## [330] "Married" "Single"  "Married" "Married" "Single"  "Single"  "Married"
## [337] "Married" "Married" "Single"  "Single"  "Single"  "Single"  "Married"
## [344] "Married" "Married" "Married" "Single"  "Married" "Married" "Married"
## [351] "Single"  "Married" "Single"  "Single"  "Married" "Married" "Married"
## [358] "Single"  "Married" "Married" "Married" "Married" "Married" "Married"
## [365] "Single"  "Single"  "Single"  "Married" "Married" "Married" "Married"
## [372] "Single"  "Married" "Single"  "Married" "Single"  "Married" "Married"
## [379] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
## [386] "Married" "Married" "Married" "Married" "Married" "Married" "Single" 
## [393] "Married" "Married" "Single"  "Married" "Married" "Married" "Married"
## [400] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
## [407] "Married" "Married" "Single"  "Single"  "Married" "Married" "Married"
## [414] "Single"  "Married" "Married" "Married" "Married" "Single"  "Single" 
## [421] "Single"  "Married" "Married" "Single"  "Married" "Single"  "Married"
## [428] "Married" "Married" "Married" "Married" "Single"  "Married" "Married"
## [435] "Single"  "Married" "Married" "Married" "Single"  "Married" "Married"
## [442] "Married" "Married" "Single"  "Single"  "Single"  "Single"  "Married"
## [449] "Married" "Married" "Single"  "Single"  "Married" "Married" "Married"
## [456] "Married" "Single"  "Single"  "Married" "Married" "Married" "Married"
## [463] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
## [470] "Single"  "Married" "Married" "Married" "Married" "Single"  "Single" 
## [477] "Married" "Married" "Married" "Married" "Married" "Single"  "Married"
## [484] "Single"  "Single"  "Married" "Married" "Married" "Married" "Single" 
## [491] "Married" "Single"  "Married" "Married" "Married" "Married" "Single" 
## [498] "Married" "Married" "Married" "Married" "Single"  "Single"  "Married"
## [505] "Single"  "Married" "Married" "Married" "Married" "Married" "Married"
## [512] "Single"  "Single"  "Single"  "Married" "Married" "Single"  "Single" 
## [519] "Single"  "Married" "Single"  "Single"  "Single"  "Married" "Single" 
## [526] "Married" "Single"  "Single"  "Married" "Married" "Single"  "Married"
## [533] "Single"  "Married" "Single"  "Married" "Married" "Married" "Single" 
## [540] "Single"  "Single"  "Married" "Married" "Married" "Married" "Single" 
## [547] "Married" "Single"  "Married" "Single"  "Married" "Married" "Married"
## [554] "Single"  "Married" "Married" "Married" "Single"  "Single"  "Married"
## [561] "Married" "Single"  "Single"  "Single"  "Married" "Single"  "Single" 
## [568] "Single"  "Single"  "Married" "Single"  "Single"  "Married" "Married"
## [575] "Married" "Single"  "Single"  "Single"  "Married" "Single"  "Married"
## [582] "Single"  "Single"  "Single"  "Single"  "Married" "Married" "Married"
## [589] "Married" "Single"  "Married" "Single"  "Single"  "Single"  "Married"
## [596] "Single"  "Single"  "Married" "Married" "Married" "Married" "Single" 
## [603] "Married" "Single"  "Married" "Married" "Single"  "Married" "Married"
## [610] "Married" "Married" "Single"  "Single"  "Married" "Single"  "Single" 
## [617] "Married" "Single"  "Married" "Single"  "Married" "Married" "Single" 
## [624] "Married" "Married" "Single"  "Married" "Married" "Married" "Married"
## [631] "Single"  "Married" "Single"  "Single"  "Married" "Married" "Married"
## [638] "Married" "Married" "Single"  "Married" "Single"  "Married" "Married"
## [645] "Single"  "Married" "Single"  "Single"  "Married" "Married" "Married"
## [652] "Married" "Married" "Single"  "Married" "Single"  "Married" "Single" 
## [659] "Married" "Married" "Married" "Married" "Single"  "Single"  "Single" 
## [666] "Married" "Married" "Single"  "Single"  "Single"  "Married" "Married"
## [673] "Single"  "Single"  "Single"  "Single"  "Married" "Married" "Married"
## [680] "Single"  "Married" "Married" "Single"  "Married" "Married" "Single" 
## [687] "Married" "Single"  "Married" "Single"  "Single"  "Married" "Married"
## [694] "Married" "Married" "Single"  "Married" "Single"  "Married" "Married"
## [701] "Single"  "Married" "Married" "Married" "Single"  "Single"  "Married"
## [708] "Single"  "Single"  "Single"  "Married" "Single"  "Single"  "Married"
## [715] "Single"  "Married" "Married" "Married" "Married" "Married" "Married"
## [722] "Married" "Single"  "Single"  "Single"  "Single"  "Single"  "Single" 
## [729] "Married" "Married" "Single"  "Married" "Married" "Single"  "Single" 
## [736] "Married" "Single"  "Single"  "Married" "Married" "Married" "Married"
## [743] "Single"  "Married" "Single"  "Married" "Single"  "Married" "Married"
## [750] "Married" "Married" "Married" "Single"  "Married" "Married" "Single" 
## [757] "Married" "Married" "Married" "Married" "Married" "Married" "Married"
## [764] "Married" "Married" "Single"  "Married" "Married" "Married" "Married"
## [771] "Married" "Single"  "Married" "Married" "Married" "Single"  "Married"
## [778] "Married" "Single"  "Married" "Single"  "Married" "Married" "Married"
## [785] "Single"  "Single"  "Single"  "Married" "Single"  "Single"  "Married"
## [792] "Single"  "Married" "Married" "Single"  "Single"  "Married" "Single" 
## [799] "Single"  "Married" "Married" "Married" "Married" "Married" "Married"
## [806] "Single"  "Married" "Married" "Married" "Married" "Single"  "Married"
## [813] "Married" "Single"  "Married" "Married" "Single"  "Single"  "Single" 
## [820] "Married" "Married" "Single"  "Married" "Single"  "Married" "Single" 
## [827] "Married" "Married" "Single"  "Single"  "Married" "Single"  "Single" 
## [834] "Single"  "Single"  "Married" "Married" "Single"  "Single"  "Single" 
## [841] "Single"  "Married" "Married" "Single"  "Single"  "Single"  "Single" 
## [848] "Married" "Single"  "Single"  "Single"  "Single"  "Single"  "Married"
## [855] "Married" "Single"  "Married" "Single"  "Married" "Single"  "Married"
## [862] "Single"  "Single"  "Married" "Single"  "Married" "Single"  "Single" 
## [869] "Married" "Married" "Single"  "Single"  "Single"  "Single"  "Single" 
## [876] "Married" "Married" "Married" "Married" "Single"  "Single"  "Married"
## [883] "Single"  "Single"  "Married" "Married" "Married" "Single"  "Married"
## [890] "Married" "Married" "Single"  "Married" "Single"  "Single"  "Single" 
## [897] "Married" "Single"  "Married" "Single"  "Single"  "Single"  "Single" 
## [904] "Married" "Single"  "Single"  "Single"  "Single"  "Single"  "Single" 
## [911] "Married" "Single"  "Single"  "Single"  "Single"

2.1 Inciso a)

¿Hay evidencia de sobredispersión en la variable ?

var(datos$art)

## [1] 3.709742

mean(datos$art)

## [1] 1.692896

Observamos que hay evidencia de sobredispersión pues la varianza es más del doble que la media.

2.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

coeftest(poissa, vcov. = vcovHC, type = "HC1")

## 
## 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 ' ' 1

El 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

2.3 Inciso c)

Obtenga la razón de tasas de incidencia (IRR) para los coeficientes e interprete los resultados.

poissonirr(art ~ mujeres + casados + kid5 + phd + mentor, data=data)

## 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 ' ' 1

Las 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.

2.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

poissonirr(poissd, data=data)

## 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 ' ' 1

Si 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.

2.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.01

Nuevamente 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.

2.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?

bnf <- glm.nb(art ~ mujeres + casados + kid5 + phd + mentor, 
              data = data)
summary(bnf)

## 
## 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

stargazer(bne, bnf, poissa, type = "text")

## 
## ===============================================================
##                                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.01

La 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.

ECNII Tarea 1

Rafael Martínez Martínez

9/13/2020