Para estudiar la estimación de los parámetros en un GLM partimos de una sola observación \(y\sim\operatorname{EDM}(\mu,\phi/w)\). Al expresar la densidad en su forma canónica se observa que la derivada del logaritmo respecto al parámetro canónico \(\theta\) es
\[ \frac{\partial\log\mathcal{P}(y;\mu,\phi/w)}{\partial\theta} \;=\;\frac{w\,(y-\mu)}{\phi}, \]
donde se emplea la relación estándar \(\mu=\mathrm{d}\kappa(\theta)/\mathrm{d}\theta\). Al cambiar a la media obtenemos
\[ \frac{\partial\log\mathcal{P}(y;\mu,\phi/w)}{\partial\mu} \;=\;\frac{w\,(y-\mu)}{\phi\,V(\mu)}, \]
pues \(\mathrm{d}\mu/\mathrm{d}\theta=V(\mu)\) por definición de la función de varianza. Esta expresión resalta que el “score” individual es simplemente la discrepancia \(y-\mu\) re-escalada por la varianza y por el peso \(w\).
En el modelo completo se especifica un predictor
\[ g(\mu)=\eta \;=\; o + \sum_{j=0}^{p}\beta_j x_{j}, \]
donde \(o\) es un posible offset y \(x_{0}=1\) provee el intercepto. Aplicando la regla de la cadena resulta
\[ \frac{\partial\log\mathcal{P}(y;\mu,\phi/w)}{\partial\beta_{j}} \;=\frac{\partial \log \mathcal{P}(y ; \mu, \phi / w)}{\partial \mu} \frac{\partial \mu}{\partial \beta_{j}}= \; (y-\mu)\, \frac{w\,x_{j}}{\phi\,V(\mu)\,\mathrm{d}\eta/\mathrm{d}\mu}. \]
Esta fórmula muestra cómo cada covariable contribuye al vector score: la discrepancia \(y-\mu\) se pondera por el peso \(w\), se normaliza por la varianza \(V(\mu)\) y por la derivada de la función enlace.
Para los segundos derivados basta diferenciar de nuevo y tomar esperanza. El término que contiene \(y-\mu\) se anula y queda
\[ \mathrm{E}\!\left[ \frac{\partial^{2}\log\mathcal{P}(y;\mu,\phi/w)} {\partial\beta_{k}\,\partial\beta_{j}} \right] \;=\; -\frac{w}{\phi\,V(\mu)} \frac{x_{j}\,x_{k}}{(\mathrm{d}\eta/\mathrm{d}\mu)^{2}}. \]
Así se obtiene una expresión compacta para la matriz de información de Fisher: los elementos diagonales y cruzados dependen solo de los pesos, de la función de varianza y de la sensibilidad del enlace, pero no del dato en sí.
Al pasar de una sola observación al conjunto completo, consideramos \(n\) respuestas independientes con \(y_{i}\sim\operatorname{EDM}(\mu_{i},\phi/w_{i})\) y predictor lineal
\[ \eta_{i}=o_{i}+\beta_{0}+\sum_{j=1}^{p}\beta_{j}x_{ji}, \]
donde los \(w_{i}\) son pesos previos conocidos y \(o_{i}\) es un posible offset. El objetivo es estimar los \(p^{\prime}=p+1\) coeficientes \(\beta\) por máxima verosimilitud, para lo cual se requieren las derivadas de la log-verosimilitud.
La log-verosimilitud es la suma de los logaritmos de las densidades individuales. Aplicando las derivadas obtenidas en la sección anterior, el score asociado a cada \(\beta_{j}\) resulta
\[ U(\beta_{j})=\frac{1}{\phi}\sum_{i=1}^{n}W_{i}\, \frac{\mathrm{d}\eta_{i}}{\mathrm{d}\mu_{i}}\,(y_{i}-\mu_{i})\,x_{ji}, \qquad j=0,\dots,p, \]
donde se define el peso de trabajo
\[ W_{i}=\frac{w_{i}}{V(\mu_{i})\,[\mathrm{d}\eta_{i}/\mathrm{d}\mu_{i}]^{2}}. \]
Estos pesos dependen de la varianza teórica \(V(\mu)\) del modelo y de la sensibilidad del enlace; su aparición es la clave de los algoritmos iterativos tipo Fisher scoring.
Tomando esperanza de las segundas derivadas (y notando que \(\mathrm{E}[y_{i}-\mu_{i}]=0\)) se obtiene la matriz de información
\[ \mathcal{I}_{jk}(\beta)=\frac{1}{\phi}\sum_{i=1}^{n}W_{i}\,x_{ji}\,x_{ki}, \qquad j,k=0,\dots,p. \]
Esta forma cuadrática sencilla muestra que, dentro de un GLM, la precisión de los estimadores depende solo de los pesos \(W_{i}\) y de los valores de las covariables.
Para un GLM de Poisson con enlace logarítmico y pesos \(w_{i}=1\), se tiene \(V(\mu)=\mu\) y \(\mathrm{d}\eta/\mathrm{d}\mu=1/\mu\). Por tanto, \(W_{i}=\mu_{i}\) y las expresiones anteriores se reducen a
\[ U(\beta_{j})=\sum_{i=1}^{n}(y_{i}-\mu_{i})\,x_{ji}, \qquad \mathcal{I}_{jk}(\beta)=\sum_{i=1}^{n}\mu_{i}\,x_{ji}\,x_{ki}. \]
En palabras, el “score” es la covarianza entre las covariables y los residuos \((y-\mu)\), mientras que la información pesa los productos de covariables por los valores ajustados \(\mu_{i}\). Esta estructura algebraica hace que la búsqueda numérica de los \(\beta\) sea eficiente y estable, una de las principales ventajas prácticas de los GLM.
El algoritmo de Fisher scoring—equivalente en la práctica a la técnica de mínimos cuadrados iterativamente ponderados (IRLS)—permite obtener las estimaciones de máxima verosimilitud de los coeficientes \(\beta\) en cualquier GLM.
La idea central es reemplazar el problema no lineal original por una sucesión de regresiones lineales –cada una usa respuestas de trabajo y pesos de trabajo que se actualizan con las medias ajustadas de la iteración previa.
Para cada observación se define
\[ z_i=\eta_i+\frac{\mathrm d\eta_i}{\mathrm d\mu_i}(y_i-\mu_i) \]
y los pesos
\[ W_i=\frac{w_i}{V(\mu_i)\,[\mathrm d\eta_i/\mathrm d\mu_i]^2} \]
donde \(V(\mu)\) es la función de varianza de la familia EDM y \(w_i\) son los pesos previos (habitualmente \(w_i=1\)).
En una iteración dada se ajusta el modelo lineal
z ~ x0 + x1 + … + xp # usando pesos = Wobteniéndose nuevos coeficientes \(\hat\beta^{(r)}\), de los que derivan \(\eta_i\), \(\mu_i\) y, con ellos, nuevos \(z_i\) y \(W_i\). El proceso se repite hasta que los cambios en \(\hat\beta\) son numéricamente despreciables.
Obsérvese que \(\phi\) no interviene en estas actualizaciones; su estimación puede posponerse.
Nota: Elección de valores iniciales
Como \(z_i\) y \(W_i\) dependen solo de las medias \(\mu_i\), basta partir con \(\mu_i^{(0)}\) próximos a \(y_i\). Si aparecen ceros que provocan logaritmos o cocientes indeterminados, se añade un pequeño desplazamiento (p. ej. \(y_i+0.1\)). En modelos binomiales es común usar la corrección \((m\,y+0.5)/(m+1)\).
En los datos de aves y eucaliptos (\(x=\) Eucs) se propone
\[ \log\mu_i=\beta_0+\beta_1x_i. \]
Para Poisson con enlace log se cumple \(V(\mu)=\mu\) y \(\mathrm d\eta/\mathrm d\mu=1/\mu\), de modo que \(W_i=\mu_i\) y
\[ z_i=\log\hat\mu_i+\frac{y_i-\hat\mu_i}{\hat\mu_i}. \tag{6.10} \]
El algoritmo arranca con \(\hat\mu_i^{(0)}=y_i+0.1\); las primeras dos iteraciones reducen drásticamente la desviación y, tras seis iteraciones, las estimaciones se estabilizan en
\[ \hat\beta_0=-0.8762,\qquad \hat\beta_1=0.1140, \]
con una desviación total \(D=63.32\). El modelo final es
\[ \log\hat\mu=-0.8762+0.1140\,x. \tag{6.11} \]
En R todo este proceso se resume en una llamada a glm() especificando family = poisson(link = "log"); internamente, el motor usa exactamente el procedimiento IRLS descrito:
Ejemplo Ajuste del GLM de Poisson propuesto anteriormente (conjunto nminer):
library(GLMsData); data(nminer)
nm.m1 <- glm(Minerab ~ Eucs,
data = nminer,
family = poisson(link = "log"),
control = list(trace = TRUE)) # Muestra la deviancia en cada iteraciónDeviance = 82.14668 Iterations - 1
Deviance = 64.49515 Iterations - 2
Deviance = 63.32603 Iterations - 3
Deviance = 63.31798 Iterations - 4
Deviance = 63.31798 Iterations - 5nm.m1
Call: glm(formula = Minerab ~ Eucs, family = poisson(link = "log"),
data = nminer, control = list(trace = TRUE))
Coefficients:
(Intercept) Eucs
-0.8762 0.1140
Degrees of Freedom: 30 Total (i.e. Null); 29 Residual
Null Deviance: 150.5
Residual Deviance: 63.32 AIC: 121.5La devianza surge como la versión, para modelos lineales generalizados, de la suma de cuadrados residual que ya conocíamos del ajuste lineal normal. Partimos de la devianza unitaria \(d(y,\mu)\), introducida en la sección anterior, que cuantifica la distancia entre cada observación \(y\) y su media teórica \(\mu\) dentro de una distribución EDM. Al sumar esas distancias ponderadas en toda la muestra obtenemos la devianza total
\[ D(y,\mu)=\sum_{i=1}^{n} w_i\,d\!\left(y_i,\mu_i\right), \]
donde los pesos \(w_i\) son conocidos y, con frecuencia, todos iguales a 1. En la práctica el interés está en la devianza residual, que se define evaluando la expresión anterior en los valores ajustados \(\hat{\mu}_i\) obtenidos tras maximizar la verosimilitud:
\[ D(y,\hat{\mu})=\sum_{i=1}^{n} w_i\,d\!\left(y_i,\hat{\mu}_i\right) \]
Minimizar \(D(y,\mu)\) equivale a maximizar la log-verosimilitud, de modo que el algoritmo IRLS puede monitorizar el progreso observando cuánta devianza resta por explicar en cada iteración. En R, el criterio de convergencia compara la diferencia relativa entre dos devianzas sucesivas y se detiene cuando esa fracción es menor que \(10^{-8}\).
A menudo se considera la devianza escalada \(D^{*}(y,\hat{\mu})=D(y,\hat{\mu})/\phi\); su distribución asintótica es \(\chi^{2}\) con \(n\) grados de libertad siempre que el saddlepoint sea una buena aproximación. Para el modelo normal, el resultado es exacto y se recupera la antigua suma de cuadrados
\[
D(y,\hat{\mu})=\sum_{i=1}^{n} w_i\,(y_i-\hat{\mu}_i)^2=\text{RSS}.
\]
En un GLM de Poisson, donde \(V(\mu)=\mu\) y \(\phi=1\), la devianza residual adopta la forma
\[ D(y,\hat{\mu})=2\sum_{i=1}^{n}\Bigl\{y_i\log\!\frac{y_i}{\hat{\mu}_i}-(y_i-\hat{\mu}_i)\Bigr\}, \]
que coincide con la devianza escalada. El ajuste al conjunto noisy miner ilustrado en la sección anterior produce \(D(y,\hat{\mu})\approx63.32\), valor que R muestra vía deviance(fit):
deviance(nm.m1)[1] 63.31798Este número mide la variabilidad restante tras incorporar la relación log-lineal propuesta entre la abundancia de aves y el número de eucaliptos; una devianza menor indicaría un modelo más capaz de explicar los datos.
Una vez obtenidos los máximos verosímiles \(\hat{\beta}_j\), el siguiente paso es cuantificar su precisión. En un GLM esto se hace a través de la matriz de información de Fisher presentada en la siguiente ecuación . Recordemos que dicha matriz resume la curvatura (negativa) de la log-verosimilitud respecto de los parámetros y, por tanto, su inversa aproxima la varianza-covarianza de los estimadores.
Sea
\[
\mathcal I(\hat{\boldsymbol{\beta}})=\frac1\phi
\begin{pmatrix}
\sum_i W_i x_{i0}x_{i0} & \cdots & \sum_i W_i x_{i0}x_{ip}\\[2pt]
\vdots & \ddots & \vdots\\[2pt]
\sum_i W_i x_{ip}x_{i0} & \cdots & \sum_i W_i x_{ip}x_{ip}
\end{pmatrix},
\] donde \(W_i\) son los pesos de trabajo. Invertir esa matriz y tomar sus elementos diagonales nos da las varianzas asintóticas de cada \(\hat{\beta}_j\). En consecuencia, el error estándar asociado se calcula como
\[ \operatorname{se}\bigl(\hat\beta_j\bigr)=\sqrt{\phi}\,v_j, \]
con \(v_j\) la raíz cuadrada del elemento \((j,j)\) de \(\mathcal I^{-1}\) construido sin el factor \(1/\phi\). Si el parámetro de dispersión \(\phi\) es conocido (por ejemplo, \(\phi=1\) en Poisson o binomial), se introduce directamente; de lo contrario, se reemplaza por su estimación.
Ejemplo Considere el Modelo Poisson ajustado a los datos de noisy miner (conjunto de datos: nminer). El resumen del GLM en R muestra los EMV de los dos coeficientes y sus errores estándar correspondientes:
coef(summary(nm.m1)) Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.8762114 0.28279293 -3.098421 1.945551e-03
Eucs 0.1139813 0.01243104 9.169092 4.770189e-20En este contexto, el vector de scores se expresa como
\[ \mathbf U=\frac{1}{\phi}\,\mathbf X^{\mathsf T}\mathbf W\mathbf M\bigl(\mathbf y-\boldsymbol\mu\bigr), \]
donde \(\mathbf X\) es la matriz de diseño, \(\mathbf W=\operatorname{diag}(W_i)\) contiene los pesos de trabajo \(W_i=w_i\!\bigl[V(\mu_i)\{d\eta_i/d\mu_i\}^2\bigr]^{-1}\) y \(\mathbf M=\operatorname{diag}(d\eta_i/d\mu_i)\) incorpora la derivada del enlace.
La matriz de información de Fisher para \(\boldsymbol\beta\) adopta la forma
\[ \mathcal I=\frac{1}{\phi}\,\mathbf X^{\mathsf T}\mathbf W\mathbf X, \]
lo que anticipa que la precisión final de los estimadores dependerá de la estructura de \(\mathbf X\) y de los pesos \(W_i\).
A partir de estas expresiones, la iteración de Fisher-scoring se describe como
\[ \hat{\boldsymbol\beta}^{(r+1)}=\hat{\boldsymbol\beta}^{(r)} +\bigl(\mathbf X^{\mathsf T}\mathbf W\mathbf X\bigr)^{-1} \mathbf X^{\mathsf T}\mathbf W\mathbf M \bigl(\mathbf y-\hat{\boldsymbol\mu}^{(r)}\bigr), \]
y puede reordenarse como un problema de mínimos cuadrados ponderados:
\[ \hat{\boldsymbol\beta}^{(r+1)}=\bigl(\mathbf X^{\mathsf T}\mathbf W\mathbf X\bigr)^{-1} \mathbf X^{\mathsf T}\mathbf W\mathbf z, \]
con respuestas de trabajo
\[ \mathbf z=\hat{\boldsymbol\eta}^{(r)}+\mathbf M \bigl(\mathbf y-\hat{\boldsymbol\mu}^{(r)}\bigr). \]
En cada ciclo se actualizan \(\hat{\boldsymbol\eta}^{(r+1)}=\mathbf o+\mathbf X\hat{\boldsymbol\beta}^{(r+1)}\) y \(\hat{\boldsymbol\mu}^{(r+1)}=g^{-1}\!\bigl(\hat{\boldsymbol\eta}^{(r+1)}\bigr)\), hasta que los cambios sucesivos en la devianza sean despreciables. Este procedimiento recibe el nombre de Iteratively Reweighted Least Squares (IRLS) porque los pesos \(\mathbf W\) cambian iteración tras iteración.
Una vez alcanzada la convergencia, la matriz de covarianzas de los estimadores se aproxima por la inversa de la información:
\[ \operatorname{var}[\hat{\boldsymbol\beta}] \;=\;\phi\bigl(\mathbf X^{\mathsf T}\mathbf W\mathbf X\bigr)^{-1}, \]
de modo que el error estándar de cada coeficiente se obtiene mediante
\[ \operatorname{se}\bigl(\hat\beta_j\bigr)=\sqrt{\phi}\,v_j, \]
donde \(v_j\) es la raíz cuadrada del \(j\)-ésimo elemento diagonal de \(\bigl(\mathbf X^{\mathsf T}\mathbf W\mathbf X\bigr)^{-1}\). Si \(\phi\) es desconocido, se reemplaza por un estimador adecuado.
El hecho de que la estimación por máxima verosimilitud en modelos lineales generalizados se implemente mediante Iteratively Reweighted Least Squares (IRLS) no es un mero truco numérico: pone de manifiesto una analogía profunda con la regresión lineal ordinaria. Cuando el algoritmo ha convergido, cada paso de la IRLS equivale, en primera aproximación, a ajustar un modelo lineal sobre unas respuestas de trabajo \(z_i\) y con pesos de trabajo \(W_i\) ya estabilizados.
El residuo de trabajo queda definido por
\[ e_i = z_i - \hat\eta_i, \]
donde \(\hat\eta_i\) es el predictor lineal final. Tanto \(e_i\) como \(W_i\) se guardan en el objeto que devuelve la función glm() (fit$residuals, fit$weights), de modo que todos los diagnósticos desarrollados para la regresión lineal—varianzas de \(\hat\beta_j\), leverage \(h_i\), residuos crudos, distancias de Cook, DFFITs, DFBETAS, etc.—pueden calcularse reutilizando esas cantidades como si se tratara de un ajuste por mínimos cuadrados ponderados. En esencia, cada GLM puede interpretarse como un modelo lineal “local” en el espacio de los \(z_i\), lo que amplía las herramientas gráficas y analíticas disponibles para evaluar la calidad del ajuste y detectar observaciones influyentes.
La fase de ajuste del modelo por IRLS no exige conocer \(\phi\); sin embargo, este parámetro resulta indispensable cuando se desea cuantificar la variabilidad de los estimadores, construir intervalos de confianza o efectuar contrastes de hipótesis. Así, salvo que \(\phi\) sea conocido de antemano, conviene estimarlo después de haber obtenido los coeficientes \(\hat\beta_j\).
En los GLM binomial y Poisson el valor teórico es \(\phi=1\), pero más adelante veremos que relajar esa suposición—permitiendo sobre- o subdispersión—puede mejorar el ajuste en datos reales.
Podríamos aplicar de forma directa el principio de máxima verosimilitud; no obstante, el MLE de \(\phi\) suele presentar un sesgo importante cuando el tamaño muestral \(n\) no es muy grande comparado con el número de parámetros \(p'\).
La regresión normal convencional ilustra bien el problema. Si \(y_i\sim N(\mu_i,\sigma^2)\), el MLE de \(\sigma^2\) es
\[ \hat\sigma^{2}= \frac1n\sum_{i=1}^{n} w_i\,(y_i-\hat\mu_i)^2, \]
pero en la práctica se prefiere el estimador insesgado:
\[ s^{2}= \frac1{\,n-p'\,}\sum_{i=1}^{n} w_i\,(y_i-\hat\mu_i)^2, \]
El objetivo de las secciones siguientes será extender este razonamiento a los distintos GLM de modo que, al especializarse en el caso normal, se recupere exactamente la fórmula precedente.
En los GLM, la verosimilitud perfilada se construye fijando \(\phi\) y maximizando respecto a los coeficientes \(\beta\). Una vez obtenidos los estimadores \(\hat\beta_0,\dots,\hat\beta_p\) para ese valor de \(\phi\), la función resultante
\[ \ell(\phi)=\ell\!\bigl(\hat\beta_0,\dots,\hat\beta_p,\phi; y\bigr) \]
suele producir un estimador de máxima verosimilitud con sesgo cuando \(n\) no es muy grande frente al número de parámetros \(p'=p+1\).
El sesgo se corrige añadiendo un término penalizador que refuerza valores de \(\phi\) más grandes. Así nace la verosimilitud de perfil modificada (MPL):
\[ \ell^{0}(\phi)=\frac{p'}{2}\,\log\phi +\ell\!\bigl(\hat\beta_0,\dots,\hat\beta_p,\phi; y\bigr). \]
El valor que maximiza esta expresión,
\[ \hat\phi^{0}=\arg\max_{\phi}\,\ell^{0}(\phi), \]
se denomina estimador MPL. Sus principales virtudes son la consistencia y un sesgo prácticamente nulo incluso con tamaños muestrales moderados.
Sin embargo, el procedimiento presenta dos inconvenientes prácticos:
Ejemplo: Caso normal
Para una GLM normal, la verosimilitud perfilada respecto a \(\sigma^{2}\) es
\[ \ell\!\left(\sigma^{2}\right)= -\tfrac12\sum_{i=1}^{n}\log\bigl(2\pi\sigma^{2}/w_{i}\bigr)\; -\;\frac{1}{2\sigma^{2}} \sum_{i=1}^{n} w_{i}\bigl(y_{i}-\hat\mu_{i}\bigr)^{2}. \]
Al derivar, igualar a cero y despejar se recupera el MLE clásico. El perfil modificado añade un término de penalización por complejidad:
\[ \ell^{0}\!\left(\sigma^{2}\right)= \frac{p^{\prime}}{2}\log\sigma^{2} -\tfrac12\sum_{i=1}^{n}\log\bigl(2\pi\sigma^{2}/w_{i}\bigr) -\frac{1}{2\sigma^{2}} \sum_{i=1}^{n} w_{i}\bigl(y_{i}-\hat\mu_{i}\bigr)^{2}. \]
Repetir el proceso (derivar, igualar a cero, resolver) produce el estimador MPL
\[ \bigl(\hat\sigma^{2}\bigr)^{0} =\frac{1}{\,n-p^{\prime}\,}\, \sum_{i=1}^{n} w_{i}\bigl(y_{i}-\hat\mu_{i}\bigr)^{2}, \]
que coincide exactamente con \(s^{2}\): es insesgado y consistente, corrigiendo el sesgo del MLE cuando el tamaño muestral no es muy grande frente al número de parámetros ajustados.
Cuando la aproximación saddlepoint es exacta, el estimador de máxima verosimilitud (MLE) de la dispersión sería el promedio simple de la devianza residual, \(D(y,\hat\mu)/n\). Sin embargo, ese cociente no toma en cuenta los \(p^{\prime}\) parámetros que podrían afectar los grados de libertad. Inspirándonos en la regresión lineal, donde se corrige dividiendo por \(n-p^{\prime}\), se define el estimador de devianza media
\[ \tilde{\phi}\;=\;\frac{D\!\left(y,\hat\mu\right)}{\,n-p^{\prime}\,}. \]
Este valor compensa la pérdida de grados de libertad y resulta prácticamente insesgado incluso con muestras moderadas.
Ejemplo
En un GLM normal la devianza residual coincide con la suma de cuadrados de los residuos (RSS). Por tanto,
\[ \tilde{\phi} \;=\;\frac{\text{RSS}}{\,n-p^{\prime}\,}=s^{2}, \]
que es justamente el estimador insesgado clásico de la varianza \(\sigma^{2}\) en la regresión lineal.
Cuando tratamos un GLM como si fuera, localmente, una regresión por mínimos cuadrados ponderados, resulta natural recurrir al estadístico de Pearson para cuantificar la variabilidad residual. Partiendo de las respuestas de trabajo \(z_i\) y de los pesos de trabajo \(W_i\) que aparecen en el algoritmo IRLS, el residuo ponderado se resume en
\[ X^{2}\;=\;\sum_{i=1}^{n} W_i\bigl(z_i-\hat\eta_i\bigr)^2 \;=\;\sum_{i=1}^{n}\frac{w_i\bigl(y_i-\hat\mu_i\bigr)^2}{V\!\bigl(\hat\mu_i\bigr)}. \]
Cada término \(\displaystyle \frac{w_i\bigl(y_i-\hat\mu_i\bigr)^2}{V(\hat\mu_i)}\) se denomina unidad Pearson y mide la discrepancia observación – ajuste descontando la varianza teórica \(V(\mu)\). Si el modelo es correcto, \(X^{2}\) debería ser comparable con los grados de libertad residuales \(n-p'\). Por ello se define el estimador Pearson de la dispersión
\[ \bar{\phi}\;=\;\frac{X^{2}}{\,n-p'\,}, \]
que, al igual que su análogo en regresión lineal, es casi insesgado y fácil de calcular.
En el caso normal (\(V(\mu)=1\)) el estadístico Pearson coincide con la suma de cuadrados residual (RSS), de modo que \(\bar\phi=s^{2}\) y coincide con el estimador de devianza media. Para la distribución Poisson (\(V(\mu)=\mu\)) resulta
\[ X^{2}\;=\;\sum_{i=1}^{n}\frac{w_i\bigl(y_i-\hat\mu_i\bigr)^2}{\hat\mu_i}, \]
lo que permite detectar sobre-o sub-dispersión comparando \(X^{2}\) con \(n-p'\).
Para el ejemplo de los cerezos se propone un GLM gamma, con componente sistemático \(\log \mu=\beta_{0}+\beta_{1}\log d+\beta_{2}\log h\). Para ajustarlo en R escribimos:
data(trees)
cherry.m1 <- glm(Volume ~ log(Height) + log(Girth),
data = trees,
family = Gamma(link = "log"))Los coeficientes estimados son
coef(cherry.m1)(Intercept) log(Height) log(Girth)
-6.691109 1.132878 1.980412 El estimador Pearson de \(\phi\) se obtiene explícitamente con:
w <- weights(cherry.m1, type = "working")
e <- residuals(cherry.m1, type = "working")
sum(w * e^2) / df.residual(cherry.m1)[1] 0.006427286Como R adopta por defecto este estimador, podemos verificarlo con:
summary(cherry.m1)$dispersion[1] 0.006427286El estimador de devianza media resulta:
deviance(cherry.m1) / df.residual(cherry.m1)[1] 0.006554117Ambas estimaciones de \(\phi\) son prácticamente idénticas.