Pregunta 1

Suponga que está interesado en una variable aleatoria que tiene una distribución Poisson con parámetro \(\lambda\). En particular: \[P(X=x)=\frac{\lambda^x e^{-\lambda}}{x!}\] Suponga que tiene una muestra de n independientes e idénticamente distribuidas.

a. [2 puntos] Plantee la función de log verosimilitud del problema.

Primero, tenemos que las \(x_i\) son \(iid\), entonces se cumple que la función de densidad conjunta: \[ L_n(\lambda|x_i,X_i):=f(x|X,\lambda) = \prod_i^{n} f(x_i | X_i, \lambda) = \prod_{1}^{n}\frac{\lambda^{x_i} \exp{(-\lambda)}}{x_i!} \] Donde: \[ \begin{align*} \mathcal{L}_n(\lambda) = & ln(L_n(\lambda)) \nonumber\\ = & \sum_{1}^{n} \left[ -\lambda + x_i ln(\lambda) - ln(x_i !) \right] \end{align*} \]

En consecuencia, la función de log verosimilitud se define como: \[ \begin{align} Q_n(\lambda) = & \frac{1}{n} \mathcal{L}_n(\lambda) \nonumber \\ = & \frac{1}{n} \sum_{1}^{n} \left[ -\lambda + x_i ln(\lambda) - ln(x_i !) \right] \end{align}\]

b. [3 puntos] Obtenga las condiciones de primer orden y resuelva para \(\lambda\).

\[ \begin{align} \frac{\partial Q_n (\lambda)}{\partial \lambda} \vert_{\hat{\lambda}} = & 0 \nonumber \\ \Rightarrow 0 = & \frac{1}{n} \sum_{1}^{n} \left[ \frac{x_i}{\lambda} -1\right] |_{\hat{ \lambda }} \nonumber \\ \Rightarrow 0 = & \frac{1}{\hat{ \lambda }} \sum_{1}^{n} \left[ x_i \right] -n \nonumber \\ \Rightarrow \hat{ \lambda } & =\frac{ \sum_{1}^{n} x_i }{n} \end{align} \]

c. [2 puntos] ¿Cuál es la media y la varianza del estimador de máxima verosimilitud que ha encontrado?

Calculando el valor esperado:

\[\begin{align} \mathbb{E}(\hat{\lambda}) = & \mathbb{E} \left[ \frac{\sum_{1}^{n} x_i}{n} \right] \nonumber \\ = & \frac{\mathbb{E} \left( \sum_{1}^{n} x_i \right)}{n} \nonumber \\ = & \frac{\sum_{1}^{n} \mathbb{E}(x_i)}{n} \nonumber \\ = & \frac{\sum_{1}^{n} \lambda}{n} \nonumber \\ = & \frac{n \lambda}{n} \nonumber \\ = & \lambda \end{align}\] Ahora, la varianza, dado que las \(x_i\) son i.i.d, se puede calcular como: \[\begin{align} Var(\hat{\lambda}) = & Var \left( \frac{\sum_{1}^{n} x_i}{n} \right) \nonumber \\ = & \frac{\sum_{1}^{n} Var \left( x_i \right)}{n^2} \nonumber \\ = & \frac{n \lambda}{n^2} \quad \nonumber \\ = & \frac{\lambda}{n} \end{align}\]

d. [1 punto] Suponga que se las \(x_i\) son tales que se puede aplicar un teorema de límite central a la media muestral. Entonces, se puede hacer la siguiente afirmación sobre el parámetro poblacional \(\lambda\):

\[P\left(-1.96\leq\left(\frac{\lambda-a}{b}\right)\leq 1.96\right)=0.95\]

Buscamos \(a\) y \(b\) que satisfacen lo anterior. Ocupando que bajo i.i.d., con \(\mathbb{E}(X_i)= \mu\) y \(Var(X_i)=\sigma^2\), se cumple: \[\frac{\bar{X}_N-\mu}{\sigma / \sqrt{N}}\xrightarrow{d}\mathcal{N}(0,1)\] Luego recordemos que si, \[Z \sim \mathcal{N}(0,1) \; \Rightarrow \; \mathbb{P} \left( -1.96 \leq Z \leq 1.96 \right) = 0.95 \] Se tiene por tanto: \[\begin{align} a = & \lambda & b= & \frac{\sqrt{\lambda}}{\sqrt{N}} \end{align}\]

Pregunta 2 (Wooldridge, 2002)

Suponga que \(y_i|\mathbf{x}_i\sim\mathcal{N}(m(\mathbf{x}_i,\mathbf{\beta}_0),\sigma_0^2)\), donde \(m(\mathbf{x},\mathbf{\beta})\) es una función del vector de variables explicativas \(\mathbf{x}\) y del vector de parámetros \(\mathbf{\beta}\) de dimensión \((k\times 1)\). Entonces, \(E(y_i|\mathbf{x}_i)=m(\mathbf{x}_i,\mathbf{\beta}_0)\) y \(V(y_i|\mathbf{x}_i)=\sigma^2_0\).

a. [2 puntos] Escriba la función de log verosimilitud condicional para la observación \(i\). Muestre que el estimador de máxima verosimilitud \(\hat{\mathbf{\beta}}\) resuelve el problema de minimización \(\min_\mathbf{\beta}\sum_i(y_i-m(\mathbf{x}_i,\mathbf{\beta}))^2\).

Sabemos que la función de log-verosimilitud es: \[\mathcal{L}(\cdot) = \sum_{i=1}^{N} \ln f(y_i | x_i, \cdot) \] Por lo tanto, para la \(i\)-íesima observación tenemos: \[\begin{equation*} \ln f(y_i | x_i, \theta) = \frac{1}{2} \ln(2 \pi) - \frac{1}{2} \ln(\sigma^2) - \frac{1}{2 \sigma^2 } \left[ y_i - m(x_i, \beta)]^2 \right] \end{equation*}\] Donde \(\theta\) es el vector de parámetros \(\theta \equiv \begin{pmatrix} \beta^{T} & \sigma^2\\ \end{pmatrix}^{T}\)

Sabemos que el estimador es aquel que maximiza la función de log-verosimilitud, es decir: \[\begin{equation*} \hat{\theta}_{MV} =\underset{\theta}{\text{arg max}} \left\lbrace - \frac{N}{2} \ln(2 \pi) - \frac{N}{2} \ln(\sigma^2) - \frac{1}{2 \sigma^2} \sum_{i=1}^{N} \left[y_i - m(x_i, \beta) \right]^2 \right\rbrace \end{equation*}\] Asi, tenemos que maximizar dicha función respecto a \(\beta\), lo cual equivale a minimizar la última sumatoria (con \(\sigma^2 >0\)) . \[ \hat{\beta}_{MV} =\underset{\beta}{\text{arg min}} \sum_{i=1}^{N} \left[y_i - m(x_i, \beta) \right]^2 \] Esto significa que el estimador de maxima verosimilitud \(\hat{\beta}\) es el estimador de mínimos cuadrados no lineales.

b. [2 puntos] Sea \(\mathbf{\theta}\equiv(\mathbf{\beta}'\;\sigma^2)'\) un vector de parámetros de dimensión \((k+1)\times 1\). Encuentre el vector score para la observación \(i\). Muestre que \(E(\mathbf{s}_i(\mathbf{\theta}_0)|\mathbf{x}_i)=\mathbf{0}\).

Sea \(\theta= ( \overbrace{(\beta_1 , \beta_2, \ldots, \beta_K)}^{\beta^T} \; \sigma^2 )^{T}\) un vector de dimensión \((K+1) \times 1\) Se define al vector score como el gradiente de la función de log-verosimilitud \[ \nabla \mathcal{L}_{N}(\theta) = \nabla \mathcal{L}_{N}(\beta^{T} \; \sigma^2)^{T} \] Calculando las derivadas tenemos que para la \(i\)-ésima observación: \[\begin{align*} \nabla_{\beta} \ln f(y_i | x_i , \beta^{T}, \sigma^2 ) = & \frac{1}{\sigma^2} \left[y_i - m(x_i, \beta) \right] \nabla_{\beta} m(x_i, \beta) \\ \frac{\partial \ln f(y_i | x_i , \beta^{T}, \sigma^2 )}{\partial \sigma^2} = & -\frac{1}{2 \sigma^2} + \frac{1}{2 \sigma^4}\left[y_i - m(x_i, \beta) \right]^2 \end{align*}\] Así, el vector score para la \(i\)-ésima observación es: \[\begin{equation*} S_i(\theta) = \begin{pmatrix} \frac{\left[ y_i - m(x_i, \beta) \right]}{\sigma^2}) \nabla_{\beta}(m(x_i, \beta))^{T} \\ -\frac{1}{2 \sigma^2} + \frac{1}{2 \sigma^4}\left[y_i - m(x_i, \beta) \right]^2 \end{pmatrix} \end{equation*}\] Evaluando el score en \(\theta_0 = (\beta_{0}^{T} \; \sigma_{0}^{2})\), se tiene \[ \begin{equation*} S_i(\theta_{0}) = \begin{pmatrix} \frac{\left[ y_i - m(x_i, \beta_0) \right]}{\sigma_{0}^{2}}) \nabla_{\beta}(m(x_i, \beta_0))^{T} \\ -\frac{1}{2 \sigma_{0}^{2}} + \frac{1}{2 \sigma_{0}^{4}}\left[y_i - m(x_i, \beta_0) \right]^2 \end{pmatrix} \end{equation*}\] \[\begin{equation*} \mathbb{E}[S_i(\theta_0) | x_i ] = \begin{pmatrix} \frac{1}{\sigma_{0}^{2}} \nabla_{\beta}(m(x_i, \beta_0))^{T} \left( \mathbb{E}[y_i |x_i] -m(x_i, \beta_0) \right) \\ -\frac{1}{2 \sigma_{0}^{2}} + \frac{1}{2 \sigma_{0}^{4}} \mathbb{E} \left[(y_i - m(x_i, \beta_0))^2|x_i \right] \end{pmatrix} \end{equation*}\] Pero sabemos que \(\mathbb{E}[y_i| x_i] = m(x_i, \beta_0)\) y \(V(y_i| x_i) =\sigma_{0}^{2}\) donde la varianza condicional se define como: \(V(y_i| x_i) = \mathbb{E}[y_i - m(x_i, \beta_0))^2|x_i]\). Entonces: \[\begin{equation*} \mathbb{E}[S_i(\theta_0) | x_i ] = \begin{pmatrix} \frac{1}{\sigma_{0}^{2}} \nabla_{\beta}(m(x_i, \beta_0))^{T} \cdot 0 \\ -\frac{1}{2 \sigma_{0}^{2}} + \frac{1}{2 \sigma_{0}^{4}} \sigma_{0}^{2} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} = \bar{0} \end{equation*}\]

c. [2 puntos] Use las condiciones de primer orden, encuentre \(\hat{\sigma}^2\) en términos de \(\hat{\mathbf{\beta}}\).

La condición de primer orden para \(\sigma_{MV}^{2}\) es: \[-\frac{N}{2 \hat{\sigma}^2} + \frac{1}{2 \sigma^{4}} \sum_{i=1}^{N} \hat{u}_i^{2} =0\] Donde \(\hat{u}_i = y_i -m(x_i, \hat{\beta})\) \[\begin{equation*} \Rightarrow \hat{\sigma}^{2} = \frac{1}{N} \sum_{i=1}^{N} \hat{u}_i^{2} \end{equation*} \]

d. [4 puntos] Encuentre la matriz hesiana de la función de log verosimilitud con respecto a \(\mathbf{\theta}\).

Dado que no conocemos la forma de \(m(\cdot)\), calculamos la Hessina para la observación \(i\)-ésima partiendo del vector score \[\begin{equation*} S_i(\theta) = \begin{bmatrix} \frac{u_i}{\sigma^2} \nabla_{\beta}(m_i(\beta))^{T} \\ -\frac{1}{2 \sigma^2} + \frac{1}{2 \sigma^4} u_{i}^{2} \end{bmatrix} \end{equation*}\]
Donde \(u_i \equiv y_i - m_i(\beta)\) y \(m_i(\beta) \equiv m_i(x_i, \beta)\) \[\begin{equation*} H_i(\theta)= \begin{bmatrix} \frac{\partial^2}{\partial \beta \partial \beta} \ln f(y_i | x_i, \theta) & \frac{\partial^2}{\partial \beta \partial \sigma^2} \ln f(y_i | x_i, \theta) \\ \frac{\partial^2}{\partial \sigma^2 \partial \beta} \ln f(y_i | x_i, \theta) & \frac{\partial^2}{\partial \sigma^2 \partial \sigma^2} \ln f(y_i | x_i, \theta) \end{bmatrix} \end{equation*}\] Usando las reglas de diferenciación: \[\begin{equation*} H_i(\theta)= \begin{bmatrix} - \nabla_{\beta}m_i(\beta)^{T} \frac{\nabla_{\beta} m_i(\beta)}{\sigma^2} + \frac{\nabla_{\beta}^{2} m_i(\beta) u_i}{\sigma^2} & - \frac{u_i \nabla_{\beta} \left(m_i(\beta) \right)^{T}}{\sigma^4} \\ - \frac{u_i \nabla_{\beta} \left(m_i(\beta) \right)^{T}}{\sigma^4} & \frac{1}{2 \sigma^4} - \frac{1}{\sigma^6}u_{i}^{2} \end{bmatrix} \end{equation*}\] Donde \(\nabla_{\beta}^{2} m_i(\beta)\) representa la matriz de segundas derivadas de \(m_i (\beta)\) respecto a \(\beta\).

e. [3 puntos] Muestre que \(-E(\mathbf{H}_i(\mathbf{\theta}_0)|\mathbf{x}_i)=E(\mathbf{s}_i(\mathbf{\theta}_0)\mathbf{s}_i(\mathbf{\theta}_0)'|\mathbf{x}_i)\).

Por un lado \(\mathbb{E}[H_i(\theta_0)|X_i]\) es \[\begin{equation*} \mathbb{E}[H_i(\theta_0)|X_i] = \begin{bmatrix} - \frac{\nabla_{\beta}m_i(\beta_0)^{T} \nabla_{\beta}m_i(\beta_0)}{\sigma_{0}^{2}} & 0 \\ 0 & -\frac{1}{2 \sigma_{0}^{4}} \end{bmatrix} \end{equation*}\] Esto porque \(\mathbb{E}[u_{i0}|x_i]=0\) y \(\mathbb{E}[u_{i0}^{2}|x_i]= \sigma_{0}^{2} \quad (\ast)\) Donde \(u_{i0}= y_i - m(x_i, \beta_0)\) \[\begin{align*} S_i(\theta_0) \cdot S_i(\theta_0)^{T} = & \begin{bmatrix} \frac{u_{i0}}{\sigma_{0}^{2}} \nabla_{\beta} (m_i(\beta_0))^{T} \\ -\frac{1}{2 \sigma_{0}^{2}} - \frac{1}{2 \sigma_{0}^{4}} u_{i0}^{2} \end{bmatrix} \begin{bmatrix} \frac{u_{i0}}{\sigma_{0}^{2}} \nabla_{\beta} (m_i(\beta_0))^{T} & -\frac{1}{2 \sigma_{0}^{2}} - \frac{1}{2 \sigma_{0}^{4}} u_{i0}^{2} \end{bmatrix}\\ = & \begin{bmatrix} \frac{u_{i0}^{2}}{\sigma_{0}^{4}} \nabla_{\beta} (m_i(\beta_0))^{T} \nabla_{\beta} (m_i(\beta_0)) & \frac{ u_{i0}}{ \sigma_{0}^{2}} \nabla_{\beta} (m_i(\beta_0))^{T} \left(-\frac{1}{2 \sigma_{0}^{2}} - \frac{1}{2 \sigma_{0}^{4}} u_{i0}^{2} \right) \\ \frac{ u_{i0}}{ \sigma_{0}^{2}} \nabla_{\beta} (m_i(\beta_0))^{T} \left(-\frac{1}{2 \sigma_{0}^{2}} - \frac{1}{2 \sigma_{0}^{4}} u_{i0}^{2} \right) & \left( \frac{1}{2 \sigma_{0}^{2}} - \frac{1}{2\sigma_{0}^{4}} u_{i0}^{2} \right)^2 \end{bmatrix} \end{align*}\] Calculando la esperanza condicional de \(S_i(\theta_0) \cdot S_i(\theta_0)^{T}\), sabiendo que se satisface \((\ast)\) entonces:

\[\begin{align*} \mathbb{E}[S_i(\theta_0) \cdot S_i(\theta_0)^{T}] = & \begin{bmatrix} \frac{\nabla_{\beta}m_{i}(\beta_{0})^{T}\nabla_{\beta}m_{i}(\beta_{0})}{\sigma^{2}} & 0\\ 0 & \mathbb{E}[(-\frac{1}{2\sigma_{0}^{2}}+\frac{1}{\sigma_{0}^{4}}u_{i0}|x_{i})] \end{bmatrix} \end{align*}\] Ahora, dado que \(y_{i}|X_{i}\sim N(m_{i}(x_{i},\beta_{0}),\sigma_{0}^{2})\), se tiene que \(\frac{1}{\sigma_{0}^{2}}u_{i0}|x_{i}\sim N(o,1)\). Llamemos \(z_{i0}=\frac{1}{\sigma_{0}^{2}}u_{i0}|x_{i}\). Como \(z_{i0}\) es normal estándar, sabemos que sus momentos impares son cero, y sus momentos pares vienen dados por \[\begin{equation*} \mathbb{E}[z_{i0}^{2n}]=\frac{2^{1-n}(2n-1)!}{(n-1)!}. \end{equation*}\] En especial, \(\mathbb{E}[z_{i0}^{4}]=3\), de donde \(\mathbb{E}[u_{i0}^{4}]=3\sigma_{0}^{4}\). Desarrollando: \[\begin{align*} \mathbb{E}\left[\frac{1}{4\sigma_{0}^{4}}-\frac{u_{i0}^{2}}{2\sigma_{0}^{6}}+\frac{u_{i0}^{4}}{4\sigma_{0}^{6}}|x_{i}\right] & = \frac{1}{4\sigma_{0}^{4}}-\frac{1}{2\sigma_{0}^{6}}\mathbb{E}[u_{i0}^{2}|x_{i}]+\frac{1}{4\sigma_{0}^{6}}\mathbb{E}[u_{i0}^{4}|x_{i}]\\ & = \frac{1}{4\sigma_{0}^{4}}-\frac{\sigma_{0}^{2}}{2\sigma_{0}^{6}}+\frac{3\sigma_{0}^{4}}{4\sigma_{0}^{6}}\\ & = \frac{1}{2\sigma_{0}^{4}} \end{align*}\] Con esto vemos que \(-\mathbb{E}[H_{i}(\theta_{0})|x_{i}]=\mathbb{E}[S_{i}(\theta_{0})S_{i}(\theta_{0})^{T}|x_{i}]\)

f. [2 puntos] Encuentre la varianza asintótica estimada de \(\hat{\mathbf{\beta}}\) y explique cómo obtendría los errores estándar.

Sabemos que el estimador de maxima versomilitud es consistente para \(\theta_0\) y \[\sqrt{N}(\hat{\theta_{MV}}- \theta_0)\xrightarrow[]{d}\mathcal{N}(0,-A_0^{-1})\] Lo que nos lleva a:

\[ \hat{\theta_{MV}} \xrightarrow[]{a}\mathcal{N}(\theta_0,-{\mathbb{E}[H(\theta_0)]}^{-1}) \] Este resultado es el límite inferior para la matr'iz de varianza de estimadores consistentes asint'oticamente normales en convergencia a la normalidad de \(\sqrt{N}(\hat{\theta}-\theta_0)\) uniforme en intervalos compactos de \(\theta_0\), es decir:

\[AVar(\sqrt{N}(\hat{\beta}-\beta_0))=\sigma_0^2{\mathbb{E}[\nabla_{\beta} (m_i(\beta_0))^{T}\nabla_{\beta} (m_i(\beta_0))]}^{-1} \] Esto es resultado de que \(\mathbb{E}[H(\theta_0)]\) es una matríz diagonal.

Así, el estimador consistente es:

\[{\hat{\sigma_0}^2}\Big\{ {\frac{1}{N}\sum_{i=1}^{N}[\nabla_{\beta} (m_i(\hat{\beta})^{T}\nabla_{\beta} (m_i(\hat{\beta}))]}\Big\}^{-1}\]

\[\widehat{AVar}(\hat{\beta})= {\hat{\sigma_0}^2}\Big\{ \sum_{i=1}^{N}[\nabla_{\beta} (m_i(\hat{\beta})^{T}\nabla_{\beta} (m_i(\hat{\beta}))]\Big\}^{-1}\]

Para hallar los errores estándar, bastaría con calcular la raíz cuadrada de cada elemento en la diagonal de la matriz de varianza asintótica estimada.

Pregunta 3

Suponga una variable aleatoria \(X_i\) con distribución desconocida. Sin embargo, sí conocemos que \(E(X)=\mu=54\) y que \(\sqrt{V(X)}=\sigma=6\). Suponga que se recolecta una muestra de 50 observaciones.

a. [1 punto] ¿Cuál es la distribución asintótica de la media muestral \(\bar{X}\)?

En nuestro caso tenemos que: \[\begin{equation*} \bar{X}=\frac{1}{N} \sum_{i=1}^{N} X_{i}. \end{equation*}\] Luego, \[\begin{equation*} \begin{aligned} \mathbb{E}\left[\bar{X}\right]&=\frac{1}{N} \sum_{i=1}^{N}\mathbb{E}\left[X_{i}\right]\\ &=\frac{N}{N} \mu=\mu. \end{aligned} \end{equation*}\] Es decir que \(\bar{X}\) es un estimador insesgedo de \(\mu\). Por otro lado \(V(\bar{X})=V\left(\frac{1}{N} \sum_{i=1}^{N} X_{i}\right)\) y como la varianza de una suma de observaciones estadisticamente independientes, cada una escogida aleatoriamente de la misma población, es la suma de sus varianzas podemos escribir: \[\begin{equation*} V(\bar{X})=\frac{1}{N^{2}} \sum_{i=1}^{n} \sigma_{X_{i}}^{2} \end{equation*}\] Más aún, como las \(X_{i}\) son identicamente distribuidas tienen la misma varianza. Así, \[\begin{equation*} \begin{aligned} V(\bar{X}) &=\frac{N \sigma^{2}}{N^{2}} \\ &=\frac{\sigma^{2}}{N} \end{aligned} \end{equation*}\] Si se satisfacen las condiciones del TLC, es decir, si el muestreo es aleatorio como hemos venido suponiendo, tenemos que: \[\begin{equation*} \sqrt{n}\left(\bar{X}_{n}-\mu\right) \stackrel{d}{\rightarrow} N\left(0, \sigma^{2}\right) \end{equation*}\] Y que: \[\begin{equation*} \bar{X}_{n} \stackrel{a}{\rightarrow} N\left(\mu, \frac{\sigma^{2}}{n}\right) \end{equation*}\]

b. [2 punto] ¿Cuál es la probabilidad de que \(\bar{X}>60\)?

Tenemos por el TLC que: \[\begin{equation*} \frac{\sqrt{n}}{\sigma}\left(\bar{X}_{n}-\mu\right) \stackrel{d}{\rightarrow} N\left(0,1\right) \end{equation*}\] Así, \[\begin{equation*} \frac{\sqrt{50}}{6}\left(\bar{X}_{n}-54\right) \stackrel{d}{\rightarrow} N\left(0,1\right) \end{equation*}\] Para \(\bar{X}=60\), el \(z-value\) es: \[\begin{equation*} \frac{\sqrt{50}}{6}\left(\bar{X}_{n}-54\right)=7.07 \end{equation*}\] Sabemos que \(99 \%\) de las medias muestrales se ubican a \(3\) desviaciones estándar de la media. Por lo tanto la probabilidad de que \(\bar{X}>60\) es menor a \(1 \%\).

c. [1 punto]¿Cuál es la probabilidad de que una observación elegida al azar sea tal que \(X_i<50\)?

Por TLC tenemos que: \[\begin{equation*} \frac{\sqrt{1}}{6}\left(X_{i}-54\right) \stackrel{d}{\rightarrow} N\left(0,1\right) \end{equation*}\] Para \(\bar{X}=50\), el \(z-value\) es: \[\begin{equation*} \frac{\sqrt{1}}{6}\left(50-54\right)=-0.66 \end{equation*}\] \(P(Z<z)= 0.2546\), la probabilidad de obtener una observación al azar y que sea menor a 50 es de \(25.26\%\).

d. [1 punto] Provea un intervalo de confianza de 90% para la media muestral.

El \(z-value\) asociado a \(5 \%\) de probabilidad a la izquierda de \(z\) es \(-1.645\). Es decir que \(90 \%\) de las observaciones aterrizan en el intervalo [-1.645, 1.645]. \[\begin{equation*} \frac{\sqrt{50}}{6}\left(\bar{x}_{2}-54\right)=-1.645 \Rightarrow \bar{X}_{1} \approx 52.60 \end{equation*}\] \[\begin{equation*} \frac{\sqrt{50}}{6}\left(\bar{x}_{2}-54\right)=1.645 \Rightarrow \bar{X}_{2} \approx 55.39 \end{equation*}\] Por lo tanto, la media muestral de una muestra aleatoria de 50 observaciones se encuentra, con \(90\%\) de confianza en el intervalo [52.60,55.39].

Pregunta 4

Sea \(x_1\) un vector de variables continuas, \(x_2\) una variable continua y \(d_1\) una variable dicotómica. Considere el siguiente modelo probit: \[P(y=1│x_1,x_2 )=\Phi(x_1'\alpha+\beta x_2+\gamma x_2^2 )\]

a. [2 punto] Provea una expresión para el efecto marginal de \(x_2\) en la probabilidad. ¿Cómo estimaría este efecto marginal?

\(x_{2}\) es una variable continua y tenemos que su efecto marginal sobre \(Pr(y=1)\) es: \[\begin{equation*} \begin{aligned} \frac{\partial \operatorname{Pr}\left(Y_{i}=1\right)}{\partial x_{i 2}}&=\frac{\partial \Phi\left(X_{i}^{\top} \beta\right)}{\partial X_{i 2}} \\ &=\frac{\partial \Phi\left(x_{i 1}^{\top} \alpha+\beta x_{i 2}+\gamma x_{i 2}^{2}\right)}{\partial x_{i 2}} \end{aligned} \end{equation*}\] Aplicando la regla de la cadena tenemo que el efecto marginal de \(x_{2}\) es: \[\begin{equation*} \frac{\partial \Phi\left(x_{i 1}^{\top} \alpha+\beta x_{i 2}+\gamma x_{i 2}^{2}\right)}{\partial\left(x_{i 1}^{\top} \alpha+\beta x_{i 2}+\partial x_{i 2}^{2}\right)} \frac{\partial\left(x_{i 1}^{\top} \alpha+\beta x_{i 2}+\gamma x_{i 2}^{2}\right)}{\partial x_{i 2}} \end{equation*}\] O bien, \[\begin{equation*} \phi\left(x_{i 1}^{\top} \alpha+\beta x_{i 2}+\gamma x_{i 2}^{2}\right)\left(\beta+2 \gamma x_{i 2}\right) \end{equation*}\] Donde \(\phi\left(x_{i 1}^{\top} \alpha+\beta x_{i 2}+\gamma x_{i 2}^{2}\right)\) es el valor de la p.d.f normal estándar evaluado en \(x_{i 1}^{\top} \alpha+\beta x_{i 2}+\gamma x_{i 2}^{2}\).

b. [2 punto] Considere ahora el modelo:

\[P(y=1│x_1 ,x_2 ,d_1)=\Phi(x_1 '\delta+\pi x_2+\rho d_1+\nu x_2 d_1 )\] Provea la nueva expresión para el efecto marginal de \(x_2\) Para evaluar el efecto marginal podriamos computar las medias muestrales de las variables independientes para con ellas evaluar \(\phi()\) y \((\beta+2 \gamma x_{i 2})\).

El efecto marginal de \(X_{2}\) sobre \(Pr(y=1)\) es:

\[\begin{equation*} \phi\left(x_{1}^{\top} \delta+\pi x_{2}+\rho d_{1}+v x_{2} d_{1}\right) \frac{\partial\left(x_{1}^{\top} \delta+\pi x_{2}+\rho d_{1}+v x_{2} d_{1}\right)}{\partial x_{12}} \end{equation*}\]

\[\begin{equation*} \begin{aligned} &=\phi\left(x_{1}^{\top} \delta+\pi x_{2}+\rho d_{1}+v x_{2} d_{1}\right)\left(\pi+v d_{1}\right)\\ &=\phi\left(x_{1}^{\top} \delta+\pi x_{2}+\rho d_{1}+v x_{2} d_{1}\right)(\pi+v), \quad si\quad d_{1}=1\\ &=\phi\left(x_{1}^{\top} \delta+\pi x_{2}+p d_{1}+v x_{2} d_{1}\right) (\pi), \quad si \quad d_{1}=0 \end{aligned} \end{equation*}\]

c. [2 punto] En el modelo de la parte b., ¿cómo evaluaría el efecto de un cambio en \(d_1\) en la probabilidad? Provea una expresión para este efecto.

El efecto marginal de \(d_{1}\) es igual al valor de \(\Phi()\) cuando \(d_{1 i}\)\(=1\) y las demás variables toman valores definidos menos el valor de \(\Phi()\) cuando \(d_{1 i}\)\(=0\) y el resto de las variables permanecen fijas en los mismos valores. Es decir: \[\begin{equation*} \Phi\left(x_{1}^{\top} \delta+\pi x_{2}+\rho+v x_{2}\right)-\Phi\left(x_{1}^{\top} \delta+\pi x_{2}\right) \end{equation*}\] Para evaluarlo se deben elegir valores de referencia para las variables \(x_{1}\) y la variable \(x_{2}\).

Pregunta 5

Considere el modelo Poisson visto en clase y un vector de variables explicativas \(x\), todas continuas, usadas para parametrizar la media.

a. [1 puntos] ¿Cuál es el efecto de un cambio en el \(j\)ésimo regresor sobre \(E(y│x)\)?

El modelo de regresión Poisson es: \[\begin{equation*} E\left(y_{i} \mid x_{i}, \beta\right)=\lambda_{i}=e^{x_{i}^{\top} \beta} \end{equation*}\] Así, el efecto marginal es: \[\begin{equation*} \frac{\partial E\left(y_{i} \mid x_{i}, \beta\right)}{\partial x_{i j}}=\frac{\partial e^{x_{i}^{\top} \beta}}{\partial x_{i j}}=\beta_{j} e^{x_{i}^{\top} \beta} \end{equation*}\] es el efecto del cambio en el j-ésimo regresor sobre \(E(y \mid x)\).

b. [2 puntos] Usando esta expresión, muestre que si el \(j\)ésimo regresor es \(\beta_j\), entonces \(100 \beta_j\) es la semielasticidad de \(E(y│x)\) con respecto a \(x_j\). Nota: Este punto es muy útil para la interpretación de los coeficientes de un modelo Poisson.

Tenemos que \[E\left(y_{i} \mid x_{i}, \beta\right)=e^{x_{i}^{\top}\beta}\] Aplicando logaritmo a ambos lados obtenemos, \[\begin{equation*} \log \left[E\left(y_{i} \mid x_{i}, \beta\right)\right]=x_{i}^{\top} \beta \end{equation*}\] Calculando la razón de cambio de \(ln E[yjxi]\), ante un cambio en el j-ésimo regresor tenemos: \[\begin{equation*} \frac{\partial \log \left[E\left(y_{i} \mid x_{i}, \beta\right)\right]}{\partial x_{i j}}=\beta_{j} \end{equation*}\] \[\begin{equation*} \Rightarrow \frac{\partial E\left(y_{i} \mid x_{i}, \beta\right)}{\partial x_{i j}} \frac{1}{E\left(y_{i} \mid x_{i}, \beta\right)}=\beta_{j} \end{equation*}\] Si \(x_{i j}\) aumenta en una unidad, \(\lambda_{i}\) aumenta en \(100^{*}\) \(\beta_{j}\%\)

c. [2 puntos] ¿Cómo se interpreta \(\beta_j\) si reemplazamos \(x_j\) por \(\log(x_j)\))?

En dicho caso tendriamos que: \[\begin{equation*} E\left(y_{i} \mid x_{i}, \beta\right)=e^{x_{-j}^{\top} \beta_{-j}+\beta_{j} \log \left(x_{i j}\right)} \end{equation*}\] \[\begin{equation*} \Rightarrow \log E(\cdot)=x_{-j}^{\top} \beta_{-j}+\beta_{j} \log \left(x_{i j}\right) \end{equation*}\] Calculando la razón de cambio de \(ln E[y_j|x_i]\), ante un cambio en el j-ésimo regresor tenemos: \[\begin{equation*} \Rightarrow \frac{\partial \log E(\cdot)}{\partial x_{i j}}=\frac{\partial E(\cdot)}{\partial x_{i j}} \frac{x_{i j}}{E(\cdot)}=\beta_{j} \end{equation*}\]

Si \(x_{ij}\) aumenta \(1\%\), \(\lambda_{i}\) cambia en \(\beta_{j}\%\). Lo cual se interpreta como la elasticidad de la esperanza condicional.

Pregunta 6 (Cameron & Trivedi, 2005)

En esta pregunta comparará el estimador de MCO, de MV y de MCNL. Antes de comenzar, recuerde fijar una semilla en R (o el software que utilice) para poder replicar sus cálculos. Se recomienda repasar la sección 5.9 en CT.

Cameron y Trivedi proveen pistas para replicar esta tabla aquí y aquí, aunque ellos trabajan en Stata. La idea es entender la anatomía de los distintos estimadores estudiados en clase.

pacman::p_load(stargazer, psych, ggplot2, data.table, MASS, sandwich, 
               lmtest, rgdal, modelsummary)

a. [2 puntos] 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, rexp requiere especificar como parámetro a la tasa en lugar de \(\lambda\).

Primero, seleccionamos una semilla, para ello elijan el número que más les guste.

set.seed(1234)

Vamos a generar la muestra

x <- rnorm(10000,1,1) #10000 obs con ~N(1,1) 
#Beta:
b1 <- 2
b2 <- -1
#Lambda:
lam <- exp(b1 + b2*x)
Ey <- 1/lam

Entonces la muestra de la muestra de distribución exponencial es:

y <- rexp(10000,lam)

b. [2 puntos] Reporte una tabla con la media, la desviación estándar, el mínimo y el máximo de \(x\), \(\lambda\) y \(y\).

Generamos las tablas de datos

datos <- as.data.table(cbind(x,y,lam))
describe(datos, skew = F)  #estadística descriptiva

c. [2 puntos] 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)\).

datos1 <- as.data.table(cbind(x,y,lam,Ey))
ggplot(datos) + geom_point(aes(x,lam))

ggplot(datos1) + geom_point(aes(x,Ey))

d. [2 puntos] 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.

reg_mco <- lm(y ~ x, datos)
stargazer(reg_mco, type = "text")
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                  y             
## -----------------------------------------------
## x                            0.590***          
##                               (0.010)          
##                                                
## Constant                      -0.0001          
##                               (0.014)          
##                                                
## -----------------------------------------------
## Observations                  10,000           
## R2                             0.252           
## Adjusted R2                    0.252           
## Residual Std. Error      1.003 (df = 9998)     
## F Statistic         3,373.201*** (df = 1; 9998)
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01
Tabla Cameron y Trivedi

Tabla Cameron y Trivedi

e. [1 punto] ¿Por qué difieren los coeficientes que obtuvo y los que se presentan en la Tabla 5.7 de CT?

Esto sucede por que obtenemos una muestra aleatoria, que no coincide precisamente con la muestra que se usa en libro.

f. [2 puntos] Obtenga los errores robustos. En R, una librería que será muy útil es sandwich.

reg_mcor <- coeftest(reg_mco, vcov =  vcovHC(reg_mco, type = "HC1"))
stargazer(reg_mcor,reg_mco, type = "text")
## 
## ===========================================================
##                               Dependent variable:          
##                     ---------------------------------------
##                                              y             
##                     coefficient             OLS            
##                        test                                
##                         (1)                 (2)            
## -----------------------------------------------------------
## x                    0.590***            0.590***          
##                       (0.021)             (0.010)          
##                                                            
## Constant              -0.0001             -0.0001          
##                       (0.014)             (0.014)          
##                                                            
## -----------------------------------------------------------
## Observations                              10,000           
## R2                                         0.252           
## Adjusted R2                                0.252           
## Residual Std. Error                  1.003 (df = 9998)     
## F Statistic                     3,373.201*** (df = 1; 9998)
## ===========================================================
## Note:                           *p<0.1; **p<0.05; ***p<0.01

g. [1 punto] ¿El estimador de MCO es consistente? ¿Por qué?

Recordemos que para que el estimador sea consistente se requiere que el proceso generador de datos este bien especificado, lo que no ocurre en este caso, ya que MCO asume linealidad y nuestra pgd es exponencial.

h. [2 puntos] Plantee la función de log verosimilitud.

Si tenemos que \(y_i\sim \exp{(\lambda_i)}\), entonces su función de densidad viene dada por: \[f(y_i)=\lambda_i e^{-\lambda y_i}\] Donde el parámetro de la distribución esta parametrizado por la muestra como:
\[ \lambda_i=\exp{(\beta_1+\beta_2x_i)} \] Entonces para una observación, tenemos que la función de log-verosimilitud es: \[ \ln f(y_i \mid \mathbf{x})=\beta_1+\beta_2x_i-\exp{(\beta_1+\beta_2x_i)}y_i \]

i. [4 puntos] 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.

La función de log-verosimilitud esta dada por: \[ \mathcal{L}_N(\cdot)=\sum{i=1}^N \ln f(y_i \mid \mathbf{x})=N\beta_1+N\beta_2\bar{x}+e^{\beta_1} \sum_{i=1}^N e^{\beta_2x_i}y_i \] Para obtener el estimador de MV, debemos maximizar la función anterior o bien minimizar su negativo. Como vemos el problema resulta particularmente complicado de resolver de forma analítica. Por lo que hacemos uso del siguiente cálculo numérico:

f <- function(B) sum(-(B[1] + x*B[2] - y*exp(B[1] + x*B[2])))
reg_ml <- nlm(f, B <- c(2,-1), hessian = T)
reg_ml
## $minimum
## [1] -80.43464
## 
## $estimate
## [1]  2.0130271 -0.9988757
## 
## $gradient
## [1] 1.115392e-06 2.572165e-06
## 
## $hessian
##          [,1]     [,2]
## [1,] 10002.00 10063.16
## [2,] 10063.16 19980.90
## 
## $code
## [1] 1
## 
## $iterations
## [1] 5

j. [3 puntos] Usando la matriz hesiana obtenida en la solución del problema de optimización, encuentre los errores estándar robustos de los coeficientes estimados.

X <- cbind(rep(1,10000),datos$x) 
index <- X%*%t(t(reg_ml$estimate))
g <- matrix(c(1-datos$y*exp(index) , datos$x - datos$y*datos$x*exp(index)),ncol = 2) #gradiente betas
B_1 <- t(g)%*%g               #Matriz B
A_1 <- solve(reg_ml$hessian)  #Matriz A
VBr_mle <- A_1%*%B_1%*%A_1    #A^(-1)*B*A^(-1)
sqrt(diag(VBr_mle))
## [1] 0.01455311 0.01013300

k. [4 puntos] 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é?

Obtenemos la estimación de MínimoS Cuadrados No Lineales

reg_nls <- nls(y ~ exp(-(b1+b2*x)), datos1,start = list(b1 = 2,b2 = -1))
summary(reg_nls)
## 
## Formula: y ~ exp(-(b1 + b2 * x))
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## b1  1.858173   0.029785   62.39   <2e-16 ***
## b2 -0.934542   0.009959  -93.84   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9046 on 9998 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 1.943e-06

Para MCNL, lo errores estándar no robustos no son válidos, porque el proceso generador de datos tiene errores heterocedásticos y exagera en gran medida la precisión de las estimaciones de MCNL.

l. [3 puntos] 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.

X <- cbind(rep(1,10000),x)
y_hat <- fitted(reg_nls,X)
D <- matrix(c(y_hat,datos1$x*y_hat), ncol = 2) # exp(b1 + b2*x) = (exp(b1+b2*x), x*exp(b1+b2*x) )

m. [3 puntos] Calcule el error estándar robusto definido como en la ecuación 5.81. En este caso \(\hat{\Omega}=Diag(\hat{u}_i^2)\).

res_nls <- residuals(reg_nls)
u_2 <- diag(res_nls%*%t(res_nls))
ue_2 <- diag(u_2,10000,10000)  #creamos Omega
B <- t(D)%*%ue_2%*%D #Matriz B
A <- t(D)%*%D        #Matriz A
VB <- solve(A)%*%B%*%solve(A)  # A^(-1)*B*A^(-1)
sqrt(diag(VB))
## [1] 0.10708159 0.04910057

n. [3 puntos] 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)\).

yhat_2 <- diag(y_hat%*%t(y_hat))
yhate_2 <- diag(yhat_2,10000,10000) #creamos nuestra nueva Omega
B <- t(D)%*%yhate_2%*%D  #Matriz B
A <- t(D)%*%D            #Matriz A
VB <- solve(A)%*%B%*%solve(A)  # A^(-1)*B*A^(-1)
sqrt(diag(VB))
## [1] 0.1500731 0.0671277

o. [1 puntos] En este experimento, ¿qué estimador tiene las mejores propiedades?

modelos <- list()
modelos[["MCO"]]<- reg_mco
modelos[["NLS"]]<- reg_nls
gg <- msummary(modelos, "html")
gg
MCO NLS
(Intercept) -0.000
(0.014)
x 0.590
(0.010)
b1 1.858
(0.030)
b2 -0.935
(0.010)
Num.Obs. 10000 10000
R2 0.252
R2 Adj. 0.252
AIC 28451.2 26378.4
BIC 28472.8 26400.1
Log.Lik. -14222.579 -13186.210
F 3373.201
finTol 0.000
isConv TRUE

EL estimador MCO es inconsistente, lo que produce estimaciones no relacionadas con beta en la función exponencial. Los estimadores de ML y NLS son consistentes y están dentro de dos errores estándar de los valores de los parámetros verdaderos \((2,-1)\), donde los errores estándar robustos deben usarse para NLS.
Para MLE los errores estándar robustos y no robustos son bastante similares. Esto se espera ya que son asintóticamente equivalentes, dado que la igualdad de la matriz de información se mantiene si el MLE se basa en la densidad real, y el tamano de la muestra aquí es grande.
Para NLS, los errores estándar no robustos no son válidos, porque el dgp tiene errores heterocedasticos y exagera en gran medida la precisión de las estimaciones de NLS.

Pregunta 7

Use la base grogger.csv para esta pregunta. Esta base contiene información sobre arrestos y características socioeconómicas de individuos arrestados.

rm(list = ls())
pacman::p_load(ggplot2,tidyverse,knitr,wooldridge,stargazer, lmtest, sandwich,
               xtable,stringr,tibble,sjPlot,sjmisc,jtools,interactions, mfx)

grogger <- read.csv("C:/Users/migue/Downloads/grogger.csv")
str(grogger)
## 'data.frame':    2725 obs. of  17 variables:
##  $ narr86 : int  0 2 1 2 1 0 2 5 0 0 ...
##  $ nfarr86: int  0 2 1 2 1 0 2 3 0 0 ...
##  $ nparr86: int  0 0 0 1 0 0 1 5 0 0 ...
##  $ pcnv   : num  0.38 0.44 0.33 0.25 0 ...
##  $ avgsen : num  17.6 0 22.8 0 0 ...
##  $ tottime: num  35.2 0 22.8 0 0 ...
##  $ ptime86: int  12 0 0 5 0 0 0 0 9 0 ...
##  $ qemp86 : num  0 1 0 2 2 4 0 0 0 3 ...
##  $ inc86  : num  0 0.8 0 8.8 8.1 ...
##  $ durat  : num  0 0 11 0 1 ...
##  $ black  : int  0 0 1 0 0 0 1 0 1 0 ...
##  $ hispan : int  0 1 0 1 0 0 0 0 0 1 ...
##  $ born60 : int  1 0 1 1 0 1 1 1 1 1 ...
##  $ pcnvsq : num  0.1444 0.1936 0.1089 0.0625 0 ...
##  $ pt86sq : int  144 0 0 25 0 0 0 0 81 0 ...
##  $ inc86sq: num  0 0.64 0 77.44 65.61 ...
##  $ arr86  : int  0 1 1 1 1 0 1 1 0 0 ...

a. [1 punto] Estime un modelo de probabilidad lineal que relacione arr86 (haber si arrestado al menos una vez en 1986) con pcnv, avgsen, tottime, ptime86, inc86, black, hispan y born60. Reporte los errores que asumen homocedasticidad y los errores robustos a heteroscedasticidad.

El ejercicio pide la siguiente regresión:

\[arr86 = \alpha_0 + \alpha_1*pcnv + \alpha_2 *avgsen + \alpha_3 *tottime +\alpha_4*inc86 + \alpha_5 *ptime86 +\alpha_6*inc86 + \alpha_7*black + \alpha_8 * hispan + \alpha_9 * born60 + e_i \] De ahí que, los resultados cuando se asumen homocedasticidad son:

reg_a <- lm(arr86 ~ 
              pcnv + avgsen + tottime + ptime86  + inc86 + black +
              hispan + born60, 
            data = grogger)
stargazer(reg_a, type = "text")
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                arr86           
## -----------------------------------------------
## pcnv                         -0.154***         
##                               (0.021)          
##                                                
## avgsen                         0.004           
##                               (0.006)          
##                                                
## tottime                       -0.002           
##                               (0.005)          
##                                                
## ptime86                      -0.022***         
##                               (0.004)          
##                                                
## inc86                        -0.001***         
##                              (0.0001)          
##                                                
## black                        0.162***          
##                               (0.024)          
##                                                
## hispan                       0.089***          
##                               (0.021)          
##                                                
## born60                         0.003           
##                               (0.017)          
##                                                
## Constant                     0.361***          
##                               (0.016)          
##                                                
## -----------------------------------------------
## Observations                   2,725           
## R2                             0.082           
## Adjusted R2                    0.080           
## Residual Std. Error      0.429 (df = 2716)     
## F Statistic          30.485*** (df = 8; 2716)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Luego, los errores robustos a heterocedasticidad son

coeftest(reg_a, vcov = vcovHC(reg_a, type = "HC0"))
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)  0.36098314  0.01668044  21.6411 < 2.2e-16 ***
## pcnv        -0.15438020  0.01893269  -8.1542 5.311e-16 ***
## avgsen       0.00350240  0.00587790   0.5959    0.5513    
## tottime     -0.00206130  0.00421859  -0.4886    0.6251    
## ptime86     -0.02159526  0.00274863  -7.8567 5.633e-15 ***
## inc86       -0.00122484  0.00011395 -10.7487 < 2.2e-16 ***
## black        0.16171828  0.02548569   6.3455 2.590e-10 ***
## hispan       0.08925863  0.02103411   4.2435 2.274e-05 ***
## born60       0.00286982  0.01713126   0.1675    0.8670    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

b. [2 punto] ¿Cuál es el efecto en la probabilidad de arresto si pcnv pasa de 0.25 a 0.75?

Sabemos que en los modelos de probabilidad lineal la intepretación de los parametros \(\alpha_i\) es directa. En consecuencía, el efecto en la probabilidad de ser arretado si pcnv pasa de 0.25 a 0.75 es igual al estimador asociado a este último por el cambio. Es decir, \[\alpha_{pcnv}*(0.75-0.25) = 0.5*\alpha_{pcnv} = 0.5(-0.15438020) = -0.0771901\]

c. [2 punto] Realice una prueba de significancia conjunta de avgsen y tottime. ¿Qué concluye?

Presentamos las pruebas de significancia conjunta sencilla y robusta a heterocedasticidad.

car::linearHypothesis(reg_a, c("avgsen = 0", "tottime = 0"))
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
car::linearHypothesis(reg_a, c("avgsen = 0", "tottime = 0"), vcov = vcovHC(reg_a, type = "HC0"))

Concluimos así, que no es posible rechazar la hipótesis nula. Es decir, las variables no sen estadísticamente significativas.

d. [2 punto] 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?

Primero, corremos el modelo probit.

probit <- glm( arr86 ~ 
              pcnv + avgsen + tottime + ptime86  + inc86 + black +
              hispan + born60, 
              family = binomial(link = "probit"),
              data = grogger)
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 ' ' 1

Luego, buscamos el efecto en la probailidad de ser arrestado 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). Entonces, tomando como base los ejemplos del laboratorio

summarise(grogger, mean(avgsen), mean(tottime), mean(inc86), mean(ptime86))
predict <- predict(probit, 
                       newdata = data.frame("pcnv" = c(0.25, 0.75),
                       "black" = c(1, 1),
                       "hispan" = c(0, 0),
                       "born60" = c(1, 1),
                       "avgsen" = c(0.6322936, 0.6322936),
                       "tottime" = c(0.8387523, 0.8387523),
                       "inc86" = c(54.96705, 54.96705),
                       "ptime86" = c(0.387156, 0.387156)),
                        type = "response")
diff(predict)
##          2 
## -0.1016607

El cambio en la probablidad es del -10.10%. Esto es, que la probabilidad disminuye en aproximadamente 10 puntos porcentuales 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.

Pregunta 8

Ahora estimará un modelo multinomial empleando la misma base motral2012.csv. El propósito será estudiar los factores relevantes para predecir la forma de ahorro que tienen las personas. 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

rm(list = ls())
pacman::p_load(Hmisc, haven, car, readr, dplyr, foreign, sandwich, lmtest,
       data.table, filehash, nnet, mlogit, margins)

motral <- read_csv("C:/Users/migue/Downloads/motral2012.csv")
## 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.

a. [1 punto] 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”.

motral <- motral %>%
  mutate(motral, ahorro = ifelse(p14 %in% 1:2, 1,
                                            ifelse(p14 %in% 7,2,
                                              ifelse(p14 %in% 3:6, 3,
                                                  ifelse(p14 %in% 8,3,NA)))))

motral1 <-motral %>%
  mutate(ahorro=as.factor(case_when(ahorro=="1" ~ "Banco",
                                                     ahorro=="2" ~ "Casa",
                                                     ahorro=="3" ~ "Otro")))

b. [2 puntos] 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”?

lgm <- multinom(ahorro ~ 
                      eda + jefe_hog,
                    data = motral1)
## # 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
stargazer(lgm, type = "text")
## 
## ==============================================
##                       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.01

Vemos que un incremento de una unidad en edad está asociado con una disminución de la probabilidad del ahorro tipo casa respecto al ahorro en el banco.

c. [2 puntos] 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?

summary(margins(lgm))

Entonces, la probabilidad de ahorrar en el banco es, en promedio, 4.36% más alto para los jefe de hogares que para los que no lo son bajo un mismo nivel de edad.

d. [4 puntos] 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”?

exp(coef(lgm))
##      (Intercept)       eda  jefe_hog
## Casa   3.0167164 0.9591874 0.7031443
## Otro   0.6548042 0.9971362 1.0081656

Esto es que, la tasa de riesgo relativo asociada a ahorrar en la opción otro dado que es jefe de hogar es 1.008 veces con respecto a no ser jefe de hogar.

e. [1 puntos] 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?

lgm2 <- multinom(relevel(ahorro, ref=2) ~ 
                   eda + jefe_hog, 
                 data = motral1)
## # weights:  12 (6 variable)
## initial  value 2727.854313 
## iter  10 value 2606.707293
## final  value 2606.707107 
## converged
stargazer(lgm2, lgm, type = "text")
## 
## =========================================================
##                             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
exp(coef(lgm2))
##       (Intercept)      eda jefe_hog
## Banco   0.3314884 1.042549 1.422183
## Otro    0.2170632 1.039563 1.433776

El coeficiente de la edad en la alternativa Banco para este último modelo tiene la misma magnitud del coeficiente en el inciso d, pero signo contrario. Esto es congruente, ya que el aumento en una unidad de edad c.r.a grupo de referencia banco, conduce a una disminución de la tasa del riesgo relativo en la alternativa de ahorro en casa.

Pregunta 9

Use la base de datos phd_articulos.csv, la cual contiene información sobre el número de artículos publicados en los últimos tres años del doctorado para una muestra de entonces estudiantes. Nuestra variable de interés será entonces art.

#rm(list = ls())
pacman::p_load(ggplot2,tidyverse,knitr,wooldridge,stargazer, lmtest, sandwich,
               xtable,stringr,tibble,sjPlot,sjmisc,jtools,interactions, MASS,
               margins, mfx)

phda <- read.csv("C:/Users/migue/Downloads/phd_articulos.csv")
str(phda)
## 'data.frame':    915 obs. of  19 variables:
##  $ caseid    : int  221 612 779 518 906 498 348 575 403 635 ...
##  $ art       : int  0 2 3 1 7 1 1 2 1 2 ...
##  $ female    : chr  "Male" "Male" "Male" "Male" ...
##  $ married   : chr  "Married" "Married" "Married" "Married" ...
##  $ kid5      : int  0 0 2 2 0 1 1 0 0 3 ...
##  $ phd       : num  2.51 2.96 1.38 3.19 3.69 ...
##  $ mentor    : int  3 2 8 5 19 5 16 1 11 8 ...
##  $ arttrunc  : int  NA 2 3 1 7 1 1 2 1 2 ...
##  $ kid5_0    : chr  "Yes" "Yes" "No" "No" ...
##  $ kid5_1    : chr  "No" "No" "No" "No" ...
##  $ kid5_2plus: chr  "No" "No" "Yes" "Yes" ...
##  $ kid5cat   : chr  "0-kids" "0-kids" "2-3-kids" "2-3-kids" ...
##  $ male      : chr  "Male" "Male" "Male" "Male" ...
##  $ phdadeq   : chr  "No" "No" "Yes" "No" ...
##  $ phdcat    : chr  "Good" "Good" "Adequate" "Strong" ...
##  $ phdelite  : chr  "No" "No" "No" "No" ...
##  $ phdgood   : chr  "Yes" "Yes" "No" "No" ...
##  $ phdstrong : chr  "No" "No" "No" "Yes" ...
##  $ profage   : num  25.4 26.6 26.4 33.6 33.4 ...

a. [1 punto] ¿Hay evidencia de sobredispersión en la variable art?

La sobredispersión está definida sobre una distribución Poisson si su media y varianza no coinciden. Entonces, calculamos la media y varianza de la variable art.

phda %>%
  summarise(Media=mean(art), Varianza= var(art))

En este caso la varianza es mayor que la media, de ahí que hay sobredispersión en la variable art.

b. [2 puntos] 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 (phd) y el número de artículos publicados por su mentor. Interprete los coeficientes estimados.

El modelo a estimar con errores robustos sería

poi <- glm(art ~ 
             factor(female) + factor(married) + kid5 + phd + mentor,
           family = poisson, data = phda )
poir <- coeftest(poi, vcov. = vcovHC, type = "HC1")
poir
## 
## z test of coefficients:
## 
##                         Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)            0.2352660  0.1347301  1.7462  0.080776 .  
## factor(female)Male     0.2245942  0.0718983  3.1238  0.001785 ** 
## factor(married)Single -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

Recordemos que la interpretación de los estimadores no puede ser directa. Sin embargo, el signo de estos últimos determinan el sentido del efecto. Vease por ejemplo que las varibales married y kid5 tienen signo negativo, lo cual indica que estar casada o la cantidad de hijos menores de 5 años disminuye la probabilidad de publicar en articulos. Ahora, calculando los efectos marginales promedio para dar una mejor interpretación se tiene que:

marg <- margins(poi)
summary(marg)

En promedio, el tener un hijo adicional menor de 5 años disminuye en 0.313 unidades la cantidad de artículos publicados manteniendo lo demás constante.

c. [2 puntos] Obtenga la razón de tasas de incidencia (IRR) para los coeficientes e interprete los resultados.

Podemos calcular la tasa de incidencia mediante una función en R o simplemente calculando la exponencial de los coeficientes.

poissonirr(poi, data = phda)
## Call:
## poissonirr(formula = poi, data = phda)
## 
## Incidence-Rate Ratio:
##                            IRR Std. Err.       z     P>|z|    
## factor(female)Male    1.251815  0.068366  4.1124 3.915e-05 ***
## factor(married)Single 0.856207  0.052549 -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
exp(poi$coefficients)
##           (Intercept)    factor(female)Male factor(married)Single 
##             1.2652453             1.2518147             0.8562068 
##                  kid5                   phd                mentor 
##             0.8312018             1.0129051             1.0258718

Ahora, los hombres tienen 1.25 más probabilidad de publicar un artículo en comparación con los hombres. En contraste, vemos que estar soltera o tener hijos disminuye la probabilidad de públicar.

d. [1 punto] 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 R. Note que la variable profage mide la duración efectiva de las carreras profesionales de cada individuo.

El nuevo modelo a estimar es el siguiente. Asimismo, calculamos la tasa de incidencia.

poi2 <- update(poi, .  ~ . + offset(profage))

phda %>%
  group_by(Sexo = female) %>%
  summarise( mean(profage), mean(art))
## `summarise()` ungrouping output (override with `.groups` argument)
poi2r <- coeftest(poi2, vcov. = vcovHC, type = "HC1")
stargazer(poi, poi2, poi2r, type = "text")
## 
## =======================================================
##                              Dependent variable:       
##                       ---------------------------------
##                                art                     
##                              Poisson        coefficient
##                                                test    
##                          (1)        (2)         (3)    
## -------------------------------------------------------
## factor(female)Male     0.225***  -15.259*** -15.259*** 
##                        (0.055)    (0.058)     (0.462)  
##                                                        
## factor(married)Single  -0.155**   0.414***     0.414   
##                        (0.061)    (0.065)     (0.524)  
##                                                        
## kid5                  -0.185***  -0.319***    -0.319   
##                        (0.040)    (0.040)     (0.420)  
##                                                        
## phd                     0.013     0.666***   0.666***  
##                        (0.026)    (0.026)     (0.231)  
##                                                        
## mentor                 0.026***   0.052***   0.052***  
##                        (0.002)    (0.003)     (0.017)  
##                                                        
## Constant               0.235**   -20.990*** -20.990*** 
##                        (0.093)    (0.093)     (1.011)  
##                                                        
## -------------------------------------------------------
## Observations             915        915                
## Log Likelihood        -1,651.056 -7,847.028            
## Akaike Inf. Crit.     3,314.113  15,706.060            
## =======================================================
## Note:                       *p<0.1; **p<0.05; ***p<0.01
poissonirr(poi2, data= phda, robust = TRUE)
## Call:
## poissonirr(formula = poi2, data = phda, robust = TRUE)
## 
## Incidence-Rate Ratio:
##                                 IRR     Std. Err.        z
## factor(female)Male    0.00000023600 0.00000010866 -33.1414
## factor(married)Single 1.51329362189 0.79017689896   0.7934
## kid5                  0.72709107279 0.30425770757  -0.7616
## phd                   1.94682123056 0.44756597499   2.8978
## mentor                1.05313504201 0.01734931548   3.1426
##                                       P>|z|    
## factor(female)Male    < 0.00000000000000022 ***
## factor(married)Single              0.427535    
## kid5                               0.446291    
## phd                                0.003758 ** 
## mentor                             0.001674 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Es visible que al tener en cuenta la hetereogeneidad de la duración efectiva de las carreras profesionales de cada individuo, las hombres tienen menos probabilidad de publicar artículos. Vale mencionar, que las diferencias son muy grandes y se deben en gran medida a que las mujeres publicaron solamente un 22% menos a pesar de estar la mitad del tiempo expuestas a la posibilidad de públicar en comparación con los hombres. Finalmente, observamos que con errores robustos, el efecto en casado/soltero desaparece.

e. [2 puntos] 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.

mbn <-  glm(art ~ 
             factor(female) + factor(married) + kid5 + phd + mentor,
           family = quasipoisson, data = phda )
mbnr <- coeftest(mbn, vcov. = vcovHC, type = "HC1")

stargazer(poi, mbn, type = "text")
## 
## ==================================================
##                           Dependent variable:     
##                       ----------------------------
##                                   art             
##                        Poisson   glm: quasipoisson
##                                     link = log    
##                          (1)            (2)       
## --------------------------------------------------
## factor(female)Male     0.225***      0.225***     
##                        (0.055)        (0.074)     
##                                                   
## factor(married)Single  -0.155**       -0.155*     
##                        (0.061)        (0.083)     
##                                                   
## kid5                  -0.185***      -0.185***    
##                        (0.040)        (0.054)     
##                                                   
## phd                     0.013          0.013      
##                        (0.026)        (0.036)     
##                                                   
## mentor                 0.026***      0.026***     
##                        (0.002)        (0.003)     
##                                                   
## Constant               0.235**        0.235*      
##                        (0.093)        (0.126)     
##                                                   
## --------------------------------------------------
## Observations             915            915       
## Log Likelihood        -1,651.056                  
## Akaike Inf. Crit.     3,314.113                   
## ==================================================
## Note:                  *p<0.1; **p<0.05; ***p<0.01

El modelo de Negativo Binomial con sobredispersión constante presenta los mismos coeficientes que el modelo Poisson. Sin embargo, aumentan los errores estándar los cuales erán altos en el modelo Poisson. Primero, vemos que ser hombre aumenta la probabilidad de publicar un artículo; en términos de su tasa de incidencia, ser hombre aumenta la probabbilidad en 1.25 veces de publicar en comparación con una mujer. Segundo, tener un hijo adicional menor de 5 años reduce la probabilidad de públicar, especificamente 0.83 veces.

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

mbn2 <- glm.nb(art ~ 
                 factor(female) + factor(married) + kid5 + phd + mentor,
               data = phda)
stargazer(poi, mbn, mbn2 ,type = "text")
## 
## ===================================================================
##                                    Dependent variable:             
##                       ---------------------------------------------
##                                            art                     
##                        Poisson   glm: quasipoisson     negative    
##                                     link = log         binomial    
##                          (1)            (2)              (3)       
## -------------------------------------------------------------------
## factor(female)Male     0.225***      0.225***          0.216***    
##                        (0.055)        (0.074)          (0.073)     
##                                                                    
## factor(married)Single  -0.155**       -0.155*          -0.150*     
##                        (0.061)        (0.083)          (0.082)     
##                                                                    
## kid5                  -0.185***      -0.185***        -0.176***    
##                        (0.040)        (0.054)          (0.053)     
##                                                                    
## phd                     0.013          0.013            0.015      
##                        (0.026)        (0.036)          (0.036)     
##                                                                    
## mentor                 0.026***      0.026***          0.029***    
##                        (0.002)        (0.003)          (0.003)     
##                                                                    
## Constant               0.235**        0.235*            0.190      
##                        (0.093)        (0.126)          (0.125)     
##                                                                    
## -------------------------------------------------------------------
## Observations             915            915              915       
## Log Likelihood        -1,651.056                      -1,561.958   
## theta                                              2.264*** (0.271)
## Akaike Inf. Crit.     3,314.113                       3,135.917    
## ===================================================================
## Note:                                   *p<0.1; **p<0.05; ***p<0.01

Incialmente, observamos que el valor de \(\hat{\alpha}\) es de 2.26 y significativo al 1% de confiabilidad. Es decir, que se rechaza la hipótesis nula de \(\alpha = 0\). Concluyendo así que es mejor utilizar un Modelo de Binomial Negativo a uno de Poisson.

Luego, para comparar los resultados de una manera más didactica comparamos las tasas de incidencia para cada modelo.

label <- c("Intercept", "factor(female)Male", "factor(married)Single", "kid5", "phd", "mentor")
a <- as_tibble_col(exp(poi$coefficients), column_name = "Poisson")
b <- as_tibble_col(exp(mbn$coefficients), column_name = "MBN Constante")
c <- as_tibble_col(exp(mbn2$coefficients), column_name = "MBN Cuadratico")
cbind(label, a , b, c )

Recordemos que para los modelos de Poisson y MBN Constante, los coeficientes no cambian. En consecuencia, las tasas de incidencia entre estos dos no cambian. No obstante, para el MBN Cuadratico la tasa de incidencia de ser hombre disminuye. En cambio, para la variable de kid5 se interpreta que el efecto disminuye un poco.