Modelos Lineales Generalizados
¶
Estimación ¶
1 Estimación
Métodos de estimación por intervalos y puntual.
Este capítulo aborda la estimación de Máxima Verosimilitud \((\mathrm{EMV})\)
Algunos EMV tienen expresiones matematicas explícitas, pero por lo general se necesitan métodos númericos.
Los métodos númericos son iterativos (Newton-Raphson o Fisher scoring). \[ x_{j+1}=x_j-\frac{f\left(x_j\right)}{f^{\prime}\left(x_j\right)} \]
2 Estimación de máxima verosimilitud
Considere las variables aleatorias independientes \(Y_1, \ldots Y_N\) que satisfacen las propiedades de los modelos lineales generalizados. Deseamos estimar los parámetros \(\boldsymbol{\beta}\) que están relacionados con las \(Y_i\) ’s a través de \(E\left(Y_i\right)=\mu_i\) y \(g\left(\mu_i\right)=\mathbf{x}_i^T \boldsymbol{\beta}\). Para cada \(Y_i\), la función de log-verosimilitud es \[ l_i=y_i b\left(\theta_i\right)+c\left(\theta_i\right)+d\left(y_i\right) \] donde las funciones \(b, c\) y \(d\) son definidas. Tambien \[ \begin{aligned} E\left(Y_i\right) & =\mu_i=-\frac{c^{\prime}\left(\theta_i\right)}{ b^{\prime}\left(\theta_i\right)} \\ \operatorname{var}\left(Y_i\right) & = \frac{\left[b^{\prime \prime}\left(\theta_i\right) c^{\prime}\left(\theta_i\right)-c^{\prime \prime}\left(\theta_i\right) b^{\prime}\left(\theta_i\right)\right]}{\left[b^{\prime}\left(\theta_i\right)\right]^3} \\ g\left(\mu_i\right) & =\mathbf{x}_i^T \boldsymbol{\beta}=\eta_i \end{aligned} \] donde \(\mathbf{x}_i\) es un vector con elementos \(x_{i j}, j=1, \ldots, p\).
La función de log-verosimilitud de todas \(Y_i\) ’s es \[ l=\sum_{i=1}^N l_i=\sum y_i b\left(\theta_i\right)+\sum c\left(\theta_i\right)+\sum d\left(y_i\right) . \] Para obtener el estimador de máxima verosimilitud para el parámetro \(\beta_j\) necesitamos \[ \begin{aligned} \frac{\partial l}{\partial \beta_j}=U_j & =\sum_{i=1}^N\left[\frac{\partial l_i}{\partial \beta_j}\right]=\sum_{i=1}^N\left[\frac{\partial l_i}{\partial \theta_i} \cdot \frac{\partial \theta_i}{\partial \mu_i} \cdot \frac{\partial \mu_i}{\partial \beta_j}\right] \\ &=\sum_{i=1}^N\left[\frac{\partial l_i}{\partial \beta_j}\right] \\ &=\sum_{i=1}^N\left[\frac{\partial (y_i b\left(\theta_i\right)+c\left(\theta_i\right)+d\left(y_i\right))}{\partial \beta_j}\right] \\ &=\sum_{i=1}^N\left[y_i b'\left(\theta_i\right)+c'\left(\theta_i\right)\right] \\ &=\sum_{i=1}^N\left[y_i b'\left(\theta_i\right)- \mu_i b'\left(\theta_i\right)\right] \\ \frac{\partial l}{\partial \beta_j}=U_j&=\sum_{i=1}^N b'\left(\theta_i\right)\left(y_i- \mu_i \right) \end{aligned} \]
- Por diferenciación.
\[ \begin{aligned} \frac{\partial \theta_i}{\partial \mu_i}&= \frac{1}{\left(\frac{\partial \mu_i}{\partial \theta_i}\right)} \\ &= \frac{1}{\left(\frac{\partial (-\frac{c^{\prime}\left(\theta_i\right)}{ b^{\prime}\left(\theta_i\right)})}{\partial \theta_i}\right)} \\ \frac{\partial \theta_i}{\partial \mu_i}&= \frac{1}{-\frac{c''(\theta)}{b'(\theta)}+\frac{c'(\theta) b''(\theta)}{(b'(\theta))^2}} \end{aligned} \] De lo cual tenemos
\[ \begin{aligned} \frac{\partial \mu_i}{\partial \theta_i} & =\frac{-c^{\prime \prime}\left(\theta_i\right)}{b^{\prime}\left(\theta_i\right)}+\frac{c^{\prime}\left(\theta_i\right) b^{\prime \prime}\left(\theta_i\right)}{\left[b^{\prime}\left(\theta_i\right)\right]^2} \\ & =b^{\prime}\left(\theta_i\right) \operatorname{var}\left(Y_i\right) \end{aligned} \]
De \(g\left(\mu_i\right) =\mathbf{x}_i^T \boldsymbol{\beta}=\eta_i\) tenemos \[ \frac{\partial \mu_i}{\partial \beta_j}=\frac{\partial \mu_i}{\partial \eta_i} \cdot \frac{\partial \eta_i}{\partial \beta_i}=\frac{\partial \mu_i}{\partial \eta_i} x_{i j} \] Por lo tanto la puntuación es \[ U_j=\sum_{i=1}^N\left[\frac{\left(y_i-\mu_i\right)}{\operatorname{var}\left(Y_i\right)} x_{i j}\left(\frac{\partial \mu_i}{\partial \eta_i}\right)\right] . \]
La matriz de varinzas-covarianzas de los \(U_j\) ’s tiene los terminos \[ \beth_{j k}=E\left[U_j U_k\right] \] que forman la matriz de información. \[ \begin{aligned} \beth_{j k} & =E\left\{\sum_{i=1}^N\left[\frac{\left(y_i-\mu_i\right)}{\operatorname{var}\left(Y_i\right)} x_{i j}\left(\frac{\partial \mu_i}{\partial \eta_i}\right)\right] \sum_{l=1}^N\left[\frac{\left(y_l-\mu_l\right)}{\operatorname{var}\left(Y_l\right)} x_{l j}\left(\frac{\partial \mu_l}{\partial \eta_l}\right)\right]\right\} \\ & =\sum \frac{E\left[\left(y_i-\mu_i\right)^2\right] x_{i j} x_{i k}}{\left[\operatorname{var}\left(Y_i\right)\right]^2}\left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 \end{aligned} \] porque \(E\left[\left(Y_i-\mu_i\right)\left(Y_l-\mu_l\right)\right]=0\) para \(i \neq l\) ya que las \(Y_i\) ’s son independientes. Usando que \(E\left[\left(y_i-\mu_i\right)^2\right]=\operatorname{var}\left(Y_i\right),(4.19)\) se puede simplificar a \[ \boxed{ \beth_{j k}=\sum_{i=1}^N \frac{x_{i j} x_{i k}}{\operatorname{var}\left(Y_i\right)}\left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2} \] Matricialmente: \[ \begin{gathered} \mathfrak{I}=\mathbf{X}^T \mathbf{W X} \\ \operatorname{donde} \mathbf{W}_{N \times N}=\operatorname{diag}\left\{w_{11}, \ldots, w_{N N}\right\} \qquad \operatorname{con} \qquad w_{i i}=\frac{1}{\operatorname{var}\left(Y_i\right)}\left\{\frac{\partial \mu_i}{\partial \eta_i}\right\}^2 \end{gathered} \] ## Estimación por método de puntuación
El método de puntuación se generaliza a
\[ \mathbf{b}^{(m)}=\mathbf{b}^{(m-1)}+\left[\beth^{(m-1)}\right]^{-1} \mathbf{U}^{(m-1)} \]
donde \(\mathbf{b}^{(m)}\) es el vector de estimaciones de los parámetros \(\beta_1, \ldots, \beta_p\) en la \(m\)-ésima iteración. En la ecuación \(\left[\beth^{(m-1)}\right]^{-1}\) es la inversa de la matriz de información con los elementos \(\beth_{j k}\) y \(\mathbf{U}^{(m-1)}\) es el vector de elementos dados. todo evaluado en \(\mathbf{b}^{(m-1)}\). Si ambos lados de la ecuación (4.21)son multiplicados por \(\beth^{(m-1)}\) obtenemos \[ \beth^{(m-1)} \mathbf{b}^{(m)}=\beth^{(m-1)} \mathbf{b}^{(m-1)}+\mathbf{U}^{(m-1)} \] \(\beth\) puede escribirse como
\[ \beth=\mathbf{X}^T \mathbf{W X} \] donde \(\mathbf{W}\) es una matriz diagonal de \(N \times N\) con elememtos \[ w_{i i}=\frac{1}{\operatorname{var}\left(Y_i\right)}\left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 \] La expresión en la parte derecha de (4.22) es el vector con elementos \[ \sum_{k=1}^p \sum_{i=1}^N \frac{x_{i j} x_{i k}}{\operatorname{var}\left(Y_i\right)}\left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 b_k^{(m-1)}+\sum_{i=1}^N \frac{\left(y_i-\mu_i\right) x_{i j}}{\operatorname{var}\left(Y_i\right)}\left(\frac{\partial \mu_i}{\partial \eta_i}\right) \] evaluado en \(\mathbf{b}^{(m-1)}\); esto se sigue de las ecuaciones (4.20) y (4.18). Así, el lado derecho de la ecuación (4.22) se puede escribir como \[ \mathbf{X}^T \mathbf{W} \mathbf{z} \] donde \(\mathbf{z}\) tiene los elementos
\[ z_i=\sum_{k=1}^p x_{i k} b_k^{(m-1)}+\left(y_i-\mu_i\right)\left(\frac{\partial \eta_i}{\partial \mu_i}\right) \] con \(\mu_i\) y \(\partial \eta_i / \partial \mu_i\) evaluadas en \(\mathbf{b}^{(m-1)}\) Por lo tanto la ecuación iterativa se puede escribir como \[ \mathbf{X}^T \mathbf{W} \mathbf{X} \mathbf{b}^{(m)}=\mathbf{X}^T \mathbf{W} \mathbf{z} \]
2.1 Ejemplo de regresión de Poisson.
Supongamos que las respuestas \(Y_i\) son variables aleatorias Poisson . En la práctica, tal supuesto se haría ya sea por motivos sustantivos o de darse cuenta de que en la figura. los variabilidad incrementa con \(Y\). Los datos de ejemplo de regresión de Poisson.
\[ \begin{array}{cccccccccc} \hline y_i & 2 & 3 & 6 & 7 & 8 & 9 & 10 & 12 & 15 \\ x_i & -1 & -1 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \\ \hline \end{array} \] datos
y=c(2,3,6,7,8,9,10,12,15)
x=c(-1,-1,0,0,0,0,1,1,1)
X=cbind(1, x)plot(x, y)
abline(v=0, col = "grey")
abline(a=7, b=5, col = "OliveDrab", lwd=2)
Esta observación apoya el uso de la distribución de Poisson, que tiene la propiedad de que el valor esperado y la varianza de \(Y_i\) son iguales \[ E\left(Y_i\right)=\operatorname{Var}\left(Y_i\right) \] Vamos a modelar la relación entre \(Y_i\) y \(x_i\) por la línea recta \[ \begin{aligned} E\left(Y_i\right) & =\mu_i=\beta_1+\beta_2 x_i \\ & =\mathbf{x}_i^T \boldsymbol{\beta} \end{aligned} \] donde \(\boldsymbol{\beta}=\left[\begin{array}{l}\beta_1 \\ \beta_2\end{array}\right]\) y \(\mathbf{x}_i=\left[\begin{array}{c}1 \\ x_i\end{array}\right]\) para \(i=1, \ldots, N\). Por lo tanto que tomamos la función liga \(g\left(\mu_i\right)\) como la función identidad \[ g\left(\mu_i\right)=\mu_i=\mathbf{x}_i^T \boldsymbol{\beta}=\eta_i \] Por lo tanto \(\partial \mu_i / \partial \eta_i=1\) lo que simplifica las ecuaciones. \[ w_{i i}=\frac{1}{\operatorname{var}\left(Y_i\right)}=\frac{1}{\beta_1+\beta_2 x_i} . \] Usando la estimación de \(\mathbf{b}=\left[\begin{array}{l}b_1 \\ b_2\end{array}\right]\), para \(\boldsymbol{\beta}\), la ecuación \(y_i\) se convierte en \[ z_i=b_1+b_2 x_i+\left(y_i-b_1+b_2 x_i\right)=y_i \] Tambien \[ \mathbf{I}=\mathbf{X}^T \mathbf{W} \mathbf{X}=\left[\begin{array}{cc} \sum_{i=1}^N \frac{1}{b_1+b_2 x_i} & \sum_{i=1}^N \frac{x_i}{b_1+b_2 x_i} \\ \sum_{i=1}^N \frac{x_i}{b_1+b_2 x_i} & \sum_{i=1}^N \frac{x_i^2}{b_1+b_2 x_i} \end{array}\right] \] y \[ \mathbf{X}^T \mathbf{W} \mathbf{z}=\left[\begin{array}{l} \sum_{i=1}^N \frac{y_i}{b_1+b_2 x_i} \\ \sum_{i=1}^N \frac{y_i x_i}{b_1+b_2 x_i} \end{array}\right] \] Las estimaciones de máxima verosimilitud se obtienen de forma iterativa a partir de las ecuaciones \[ \left(\mathbf{X}^T \mathbf{W X}\right)^{(m-1)} \mathbf{b}^{(m)}=\mathbf{X}^T \mathbf{W} \mathbf{z}^{m-1} \] donde el superindice \({ }_{(m-1)}\) denota la evaluación en \(\mathbf{b}^{(m-1)}\). Para estos datos \(N=9\) \[ \mathbf{y}=\mathbf{z}=\left[\begin{array}{c} 2 \\ 3 \\ 6 \\ 7 \\ 8 \\ 9 \\ 10 \\ 12 \\ 15 \end{array}\right] \quad y \quad \mathbf{X}=\left[\begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_9 \end{array}\right]=\left[\begin{array}{cc} 1 & -1 \\ 1 & -1 \\ 1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 1 \\ 1 & 1 \\ 1 & 1 \end{array}\right] \] De la figura 4.5 podemos obtener estimaciones iniciales \(b_1^{(1)}=7 \mathrm{y} b_2^{(1)}=5\). Por lo tanto \[ \left(\mathbf{X}^T \mathbf{W X}\right)^{(1)}=\left[\begin{array}{cc} 1.821429 & -0.75 \\ -0.75 & 1.25 \end{array}\right], \quad\left(\mathbf{X}^T \mathbf{W z}\right)^{(1)}=\left[\begin{array}{c} 9.869048 \\ 0.583333 \end{array}\right] \] así \[ \begin{aligned} \mathbf{b}^{(2)} & =\left[\left(\mathbf{X}^T \mathbf{W} \mathbf{X}\right)^{(1)}\right]^{-1}\left(\mathbf{X}^T \mathbf{W} \mathbf{z}\right)^{(1)} \\ & =\left[\begin{array}{cc} 0.729167 & 0.4375 \\ 0.4375 & 1.0625 \end{array}\right]\left[\begin{array}{c} 9.869048 \\ 0.583333 \end{array}\right] \\ & =\left[\begin{array}{c} 7.4514 \\ 4.9375 \end{array}\right] . \end{aligned} \] Este proceso iterativo se continúa hasta que converge.
\[ \begin{aligned} \mathbf{b}^{(m+1)} & =\left[\left(\mathbf{X}^T \mathbf{W} \mathbf{X}\right)^{(m)}\right]^{-1}\left(\mathbf{X}^T \mathbf{W} \mathbf{z}\right)^{(m)} \\. \end{aligned} \] Vector inicial
Beta1=c(7, 5); Beta1[1] 7 5\[ w_{i i}=\frac{1}{\operatorname{var}\left(Y_i\right)}=\frac{1}{\beta_1+\beta_2 x_i} . \] La matriz \(W\) se obtiene medianate:
W= solve( diag(as.vector(X %*% Beta1)) ); round(W, 3) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.5 0.0 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[2,] 0.0 0.5 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[3,] 0.0 0.0 0.143 0.000 0.000 0.000 0.000 0.000 0.000
[4,] 0.0 0.0 0.000 0.143 0.000 0.000 0.000 0.000 0.000
[5,] 0.0 0.0 0.000 0.000 0.143 0.000 0.000 0.000 0.000
[6,] 0.0 0.0 0.000 0.000 0.000 0.143 0.000 0.000 0.000
[7,] 0.0 0.0 0.000 0.000 0.000 0.000 0.083 0.000 0.000
[8,] 0.0 0.0 0.000 0.000 0.000 0.000 0.000 0.083 0.000
[9,] 0.0 0.0 0.000 0.000 0.000 0.000 0.000 0.000 0.083si \[ \begin{aligned} A= &\left[\left(\mathbf{X}^T \mathbf{W} \mathbf{X}\right)^{(m)}\right] \\ B= & \left(\mathbf{X}^T \mathbf{W} \mathbf{z}\right) \end{aligned} \] Entonces nuestra b toma la forma:
\[ \begin{aligned} \mathbf{b}^{(m+1)} & =A^{-1}B \end{aligned} \]
A= t(X) %*% W %*% X; A x
1.821429 -0.75
x -0.750000 1.25B= t(X)%*%W%*%y; B [,1]
9.8690476
x 0.5833333Beta2 = qr.solve(A, B); round(Beta2, 5) [,1]
7.45139
x 4.93750Todo junto seria:
Beta2 = qr.solve(t(X) %*% W %*% X, t(X)%*%W%*%y); round(Beta2, 5) [,1]
7.45139
x 4.93750abs(Beta2-Beta1) [,1]
0.4513889
x 0.0625000max(abs(Beta2-Beta1) )[1] 0.4513889# Paso 3
W= solve( diag(as.vector(X%*% Beta2)) ); round(W, 3) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.398 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[2,] 0.000 0.398 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[3,] 0.000 0.000 0.134 0.000 0.000 0.000 0.000 0.000 0.000
[4,] 0.000 0.000 0.000 0.134 0.000 0.000 0.000 0.000 0.000
[5,] 0.000 0.000 0.000 0.000 0.134 0.000 0.000 0.000 0.000
[6,] 0.000 0.000 0.000 0.000 0.000 0.134 0.000 0.000 0.000
[7,] 0.000 0.000 0.000 0.000 0.000 0.000 0.081 0.000 0.000
[8,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.081 0.000
[9,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.081A = t(X)%*%W%*%X; A x
1.5745453 -0.5534276
x -0.5534276 1.0377326B = t(X)%*%W%*%y; B [,1]
9.0015924
x 0.9975968Beta3 = qr.solve(A, B); round(Beta3, 5) [,1]
7.45163
x 4.93531abs(Beta3-Beta2) [,1]
0.0002429178
x 0.0021862604max(abs(Beta3-Beta2) )[1] 0.00218626# Paso 4
W= solve( diag(as.vector(X%*% Beta3)) ); round(W, 3) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.397 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[2,] 0.000 0.397 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[3,] 0.000 0.000 0.134 0.000 0.000 0.000 0.000 0.000 0.000
[4,] 0.000 0.000 0.000 0.134 0.000 0.000 0.000 0.000 0.000
[5,] 0.000 0.000 0.000 0.000 0.134 0.000 0.000 0.000 0.000
[6,] 0.000 0.000 0.000 0.000 0.000 0.134 0.000 0.000 0.000
[7,] 0.000 0.000 0.000 0.000 0.000 0.000 0.081 0.000 0.000
[8,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.081 0.000
[9,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.081A= t(X)%*%W%*%X; A x
1.5737977 -0.5526216
x -0.5526216 1.0370025B= t(X)%*%W%*%y; B [,1]
9.0000096
x 0.9999854Beta4 = qr.solve(A, B); round(Beta4, 5) [,1]
7.45163
x 4.93530abs(Beta4-Beta3) [,1]
1.473804e-06
x 1.326424e-05max(abs(Beta4-Beta3) )[1] 1.326424e-05\[ \begin{array}{ccccc} \hline m & 1 & 2 & 3 & 4 \\ b_1^{(m)} & 7 & 7.45139 & 7.45163 & 7.45163 \\ b_2^{(m)} & 5 & 4.93750 & 4.93531 & 4.93530 \\ \hline \end{array} \]
Utilizando el paquete GLM.
res.p=glm(y~x,family=poisson(link="identity"))
summary(res.p)
Call:
glm(formula = y ~ x, family = poisson(link = "identity"))
Deviance Residuals:
Min 1Q Median 3Q Max
-0.7019 -0.3377 -0.1105 0.2958 0.7184
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.4516 0.8841 8.428 < 2e-16 ***
x 4.9353 1.0892 4.531 5.86e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 18.4206 on 8 degrees of freedom
Residual deviance: 1.8947 on 7 degrees of freedom
AIC: 40.008
Number of Fisher Scoring iterations: 3