log(media)
1.079 Modelos Lineales Generalizados (GLM): Definición y fundamentos
EP7120-Modelos Lineales Generalizados Aplicados
Enver Gerald Tarazona Vargas
Universidad Nacional Agraria La Molina (UNALM), Perú
Introducción
El problema de partida
El modelo lineal normal fue usado como herramienta básica para estudiar relaciones entre una respuesta y un conjunto de covariables.
Sin embargo, en muchas aplicaciones esta suposición puede resultar forzada:
\[ Y_i \sim N(\mu_i,\sigma^2) \]
El problema no es solo ajustar un modelo, sino representar adecuadamente la naturaleza de la variable respuesta.
La respuesta no siempre es normal
En la práctica, la variable respuesta puede tener restricciones naturales:
\[ Y_i \geq 0 \]
\[ 0 \leq Y_i \leq 1 \]
o puede tomar solo valores discretos.
Ejemplos:
- número de eventos;
- presencia o ausencia de una característica;
- proporción de éxitos;
- tiempos o costos positivos.
El modelo lineal normal en tres componentes
El modelo lineal normal puede separarse en tres partes porque distingue:
- la distribución asumida para la respuesta;
- la forma en que entran las covariables;
- la relación entre la media y el predictor lineal.
En el modelo lineal normal, la función de enlace es la identidad.
La alternativa clásica: transformar la respuesta
Una alternativa usual fue transformar la variable respuesta para aproximar las condiciones del modelo normal.
Por ejemplo:
\[ Z_i = h(Y_i) \]
con transformaciones como:
\[ \log(Y_i), \qquad \sqrt{Y_i} \]
La idea era ajustar un modelo lineal normal sobre la respuesta transformada.
Transformación de Box-Cox
Una familia de transformaciones ampliamente usada fue propuesta por Box y Cox:
\[ Z_i = \begin{cases} \dfrac{Y_i^\lambda-1}{\lambda}, & \lambda \neq 0, \\[3mm] \log(Y_i), & \lambda = 0. \end{cases} \]
donde \(\lambda\) es un parámetro que debe estimarse.
Qué se buscaba con la transformación
La transformación buscaba que la nueva respuesta (Z_i) cumpliera aproximadamente:
\[ Z_i \sim N \]
\[ Var(Z_i)=\text{constante} \]
\[ E(Z_i)=\eta_i \]
Es decir, se intentaba recuperar normalidad, homogeneidad de varianza y linealidad.
Limitación de esta estrategia
En la práctica, una misma transformación rara vez logra simultáneamente:
- normalidad;
- varianza constante;
- linealidad en la media.
Además, el modelo pasa a explicar:
\[ E(Z_i)=E\{h(Y_i)\} \]
y no directamente:
\[ E(Y_i)=\mu_i \]
Cambio de enfoque
Los modelos lineales generalizados surgen como una alternativa a esta estrategia.
En lugar de transformar la respuesta para forzar normalidad, se busca modelar directamente:
\[ \mu_i=E(Y_i) \]
respetando la naturaleza de la variable respuesta.
Idea que abre la siguiente sección
La generalización consiste en mantener una estructura de regresión, pero permitir:
- una distribución más adecuada para \(Y_i\);
- una relación más flexible entre \(\mu_i\) y el predictor lineal.
Con esto se da paso a la definición formal de un modelo lineal generalizado.
Definición y estructura de los GLM
Definición
Componente aleatorio: Sean \(Y_1,\ldots,Y_n\) v.a. independientes tal que
\[ Y_i \sim FE(\mu_i,\phi) \]
esto es
\[ E(Y_i)=\mu_i,\qquad i=1,\ldots,n. \]
Componente sistemático:
\[ \eta_i=\beta_1x_{i1}+\cdots+\beta_px_{ip}=x_i^\top\beta \]
donde \(\eta_i\) es denominado predictor lineal.
Función de enlace:
\[ g(\mu_i)=\eta_i \]
donde \(g(\cdot)\) es una función monótona y diferenciable.
Componente aleatorio
El componente aleatorio especifica la distribución de la variable respuesta.
En los GLM se asume que \(Y_i\) pertenece a la familia exponencial:
\[ f(y;\theta,\phi) = \exp\left\{ \phi[y\theta-b(\theta)]+c(y,\phi) \right\}. \]
donde:
- \(\theta\): parámetro canónico;
- \(\phi>0\): parámetro de precisión;
- \(b(\theta)\): función cumulante;
- \(c(y,\phi)\): función que no depende de \(\theta\).
Media, varianza y función de varianza
Si \(Y\) pertenece a la familia exponencial, entonces:
\[ E(Y)=\mu=b'(\theta) \]
y
\[ Var(Y)=\phi^{-1}b''(\theta). \]
Como:
\[ b''(\theta)=V(\mu), \]
se escribe:
\[ Var(Y)=\phi^{-1}V(\mu). \]
La función \(V(\mu)\) se denomina función de varianza y describe cómo cambia la varianza de \(Y\) cuando cambia su media.
Función de varianza
En la notación de Paula, para una variable de la familia exponencial:
\[ Var(Y)=\phi^{-1}V(\mu). \]
La varianza se descompone en dos partes:
\[ Var(Y) = \underbrace{\phi^{-1}}_{\text{dispersión}} \cdot \underbrace{V(\mu)}_{\text{relación media-varianza}}. \]
Aquí \(\phi>0\) es el parámetro de precisión, por lo que \(\phi^{-1}\) puede leerse como factor de dispersión.
¿Qué representa \(V(\mu)\)?
La función \(V(\mu)\) describe cómo cambia la variabilidad de la respuesta cuando cambia su media.
Por ejemplo:
\[ V(\mu)=1 \]
indica varianza constante.
\[ V(\mu)=\mu \]
indica que la varianza crece proporcionalmente con la media.
\[ V(\mu)=\mu^2 \]
indica que la varianza crece cuadráticamente con la media.
Propiedad asintótica
Si \(Y\) pertenece a la familia exponencial, entonces:
\[ \sqrt{\phi}(Y-\mu) \xrightarrow{d} N(0,V(\mu)), \qquad \phi\to\infty. \]
Equivalente:
\[ Y \approx N\left(\mu,\phi^{-1}V(\mu)\right) \]
para valores grandes de \(\phi\).
Esta propiedad muestra el papel de \(\phi\) como parámetro de precisión.
Interpretación de la función de varianza
| Función de varianza | Interpretación |
|---|---|
| \(V(\mu)=1\) | Varianza constante |
| \(V(\mu)=\mu\) | La varianza crece proporcionalmente con la media |
| \(V(\mu)=\mu(1-\mu)\) | La varianza depende de una media acotada entre 0 y 1 |
| \(V(\mu)=\mu+\mu^2/k\) | La varianza crece más que en Poisson |
| \(V(\mu)=\mu^2\) | La varianza crece cuadráticamente con la media |
| \(V(\mu)=\mu^3\) | La varianza crece más rápidamente con la media |
Función de varianza y elección del modelo
En aplicaciones, la relación media-varianza orienta la elección del componente aleatorio.
- Si la varianza es aproximadamente constante: modelo normal.
- Si la varianza crece con la media: modelo Poisson.
- Si la varianza crece más que la media: binomial negativa u otro modelo con sobredispersión.
- Si la respuesta es positiva y la varianza crece con una potencia de la media: Gamma o normal inversa.
Dentro de la familia exponencial, \(V(\mu)\) no solo describe variabilidad: también ayuda a caracterizar la distribución asociada.
La función de varianza conecta:
\[ \text{naturaleza de }Y_i \quad\longrightarrow\quad \text{distribución} \quad\longrightarrow\quad V(\mu). \]
Función de varianza: exploración empírica
Una forma exploratoria de estudiar la relación media-varianza es agrupar los datos y calcular, para cada grupo:
\[ \bar{y}_g \qquad \text{y} \qquad s_g^2. \]
Luego se grafica:
\[ \log(s_g^2) \quad \text{vs.} \quad \log(\bar{y}_g). \]
Si se observa una relación aproximadamente lineal:
\[ \log(s_g^2)\approx a+b\log(\bar{y}_g), \]
entonces:
\[ Var(Y)\propto \mu^b. \]
Así, el gráfico ayuda a identificar una posible función de varianza.
Familia exponencial: casos principales
| Distribución | Tipo de respuesta | Soporte de \(Y\) | Rango de \(\mu\) | \(\theta\) | \(b(\theta)\) | \(\phi\) | \(V(\mu)\) |
|---|---|---|---|---|---|---|---|
| Normal | Continua simétrica | \(\mathbb{R}\) | \(\mathbb{R}\) | \(\mu\) | \(\theta^2/2\) | \(\sigma^{-2}\) | \(1\) |
| Poisson | Conteo | \(\mathbb{N}_0\) | \(\mathbb{R}^+\) | \(\log(\mu)\) | \(e^\theta\) | \(1\) | \(\mu\) |
| Binomial | Binaria / proporción | \(\{0,1,\ldots,m\}/m\) | \((0,1)\) | \(\log\left(\dfrac{\mu}{1-\mu}\right)\) | \(\log(1+e^\theta)\) | \(m\) | \(\mu(1-\mu)\) |
| Binomial negativa | Conteo con sobredispersión | \(\mathbb{N}_0\) | \(\mathbb{R}^+\) | \(\log\left(\dfrac{\mu}{\mu+k}\right)\) | \(-k\log(1-e^\theta)\) | \(1\) | \(\mu+\dfrac{\mu^2}{k}\) |
| Gamma | Positiva asimétrica | \(\mathbb{R}^+\) | \(\mathbb{R}^+\) | \(-\dfrac{1}{\mu}\) | \(-\log(-\theta)\) | \(\dfrac{1}{CV^2}\) | \(\mu^2\) |
| Normal inversa | Positiva asimétrica | \(\mathbb{R}^+\) | \(\mathbb{R}^+\) | \(-\dfrac{1}{2\mu^2}\) | \(-\sqrt{-2\theta}\) | \(\phi\) | \(\mu^3\) |
Nota. \(\mathbb{N}_0=\{0,1,2,\ldots\}\); \(\mathbb{R}^+=(0,\infty)\); en la binomial, \(m\) representa el número de ensayos y la respuesta se expresa como proporción \(Y^*\) de éxitos; en la binomial negativa, \(k>0\) es el parámetro de dispersión; \(CV\) representa el coeficiente de variación. Además, para la notación de la tabla: \(Y^*\sim B(m,\mu)\), donde \(Y^*\) es la proporción de éxitos en \(m\) ensayos independientes; \(Y\sim BN(\mu,k)\), con \(\mu>0\); y \(Y\sim G(\mu,\phi)\), donde \(\phi^{-1}\) corresponde al coeficiente de variación. La función \(V(\mu)\) corresponde a la función de varianza.
Ejemplo: contexto de simulación
Se simula un estudio de eventos de plaga en parcelas agrícolas.
Las parcelas se agrupan según su nivel de exposición al riesgo. Para el grupo \(g\), se asume que el número esperado de eventos es:
\[ \lambda_g \in [1,20]. \]
Para cada parcela \(i\) del grupo \(g\) se genera un conteo:
\[ Y_{ig}\sim Poisson(\lambda_g). \]
Como en una distribución Poisson:
\[ E(Y)=Var(Y)=\mu, \]
se espera observar:
\[ V(\mu)=\mu. \]
Ejemplo: gráfico media-varianza
Ejemplo: lectura del gráfico
Para cada grupo se calculó:
\[ \bar{y}_g \qquad \text{y} \qquad s_g^2. \]
Luego se graficó:
\[ \log(s_g^2) \quad \text{frente a} \quad \log(\bar{y}_g). \]
Si:
\[ Var(Y)\propto \mu^b, \]
entonces:
\[ \log{Var(Y)}=a+b\log(\mu). \]
En esta simulación, la pendiente estimada es:
Ejemplo: conclusión
La pendiente estimada es cercana a 1.
Por tanto:
\[ Var(Y)\propto \mu. \]
Esto coincide con lo esperado para una respuesta Poisson:
\[ V(\mu)=\mu. \]
El ejemplo muestra que la relación media-varianza puede explorarse antes de elegir el componente aleatorio del GLM.
Ver código en R
set.seed(123)
n_grupos <- 30
tam_grupo <- 40
lambda_grupo <- seq(1, 20, length.out = n_grupos)
datos_var <- data.frame(
grupo = rep(1:n_grupos, each = tam_grupo),
y = unlist(lapply(lambda_grupo, function(lambda) {
rpois(tam_grupo, lambda = lambda)
}))
)
resumen_grupo <- aggregate(
y ~ grupo,
data = datos_var,
FUN = function(z) c(media = mean(z), varianza = var(z))
)
resumen_grupo <- data.frame(
grupo = resumen_grupo$grupo,
media = resumen_grupo$y[, "media"],
varianza = resumen_grupo$y[, "varianza"]
)
ajuste_mv <- lm(log(varianza) ~ log(media), data = resumen_grupo)
plot(
log(resumen_grupo$media),
log(resumen_grupo$varianza),
pch = 19,
xlab = "log(media por grupo)",
ylab = "log(varianza por grupo)",
main = "Relación empírica media-varianza"
)
abline(ajuste_mv, lwd = 2)
coef(ajuste_mv)Componente sistemático
El componente sistemático define cómo ingresan las covariables al modelo mediante el predictor lineal:
\[ \eta_i=x_i^\top\beta. \]
Es decir:
\[ \eta_i = \beta_1x_{i1} + \beta_2x_{i2} + \cdots + \beta_px_{ip}. \]
Aquí pueden incorporarse variables cuantitativas, variables cualitativas codificadas, interacciones o transformaciones de covariables.
Función de enlace
La función de enlace conecta la media de la respuesta con el predictor lineal:
\[ g(\mu_i)=\eta_i. \]
La función \(g(\cdot)\) debe ser:
- monótona;
- diferenciable;
- invertible en el rango de \(\mu_i\).
Por tanto:
\[ \mu_i=g^{-1}(\eta_i). \]
Su papel es permitir que \(\mu_i\) permanezca en el rango adecuado para cada tipo de respuesta:
\[ \mu_i>0 \]
para conteos o respuestas positivas, y
\[ 0<\mu_i<1 \]
para probabilidades o proporciones.
Enlaces comunes
| Enlace | Forma | Transforma | Uso frecuente |
|---|---|---|---|
| Identidad | \(\mu_i=\eta_i\) | \(\mu_i\in\mathbb{R}\) hacia \(\eta_i\in\mathbb{R}\) | Respuestas continuas aproximadamente normales |
| Log | \(\log(\mu_i)=\eta_i\) | \(\mu_i>0\) hacia \(\eta_i\in\mathbb{R}\) | Conteos, tasas y respuestas positivas |
| Logit | \(\log\left(\dfrac{\mu_i}{1-\mu_i}\right)=\eta_i\) | \(0<\mu_i<1\) hacia \(\eta_i\in\mathbb{R}\) | Respuestas binarias o proporciones |
| Probit | \(\Phi^{-1}(\mu_i)=\eta_i\) | \(0<\mu_i<1\) hacia \(\eta_i\in\mathbb{R}\) | Respuestas binarias o proporciones |
| Complemento log-log | \(\log\{-\log(1-\mu_i)\}=\eta_i\) | \(0<\mu_i<1\) hacia \(\eta_i\in\mathbb{R}\) | Respuestas binarias con efecto asimétrico |
| Inverso | \(\mu_i^{-1}=\eta_i\) | \(\mu_i>0\) hacia \(\eta_i\) restringido | Respuestas positivas; enlace canónico Gamma |
| Inverso cuadrático | \(\mu_i^{-2}=\eta_i\) | \(\mu_i>0\) hacia \(\eta_i>0\) | Respuestas positivas; enlace canónico normal inversa |
Nota. La misma función de enlace puede usarse con más de una distribución. La elección del enlace debe considerar el rango de \(\mu_i\), la interpretación de los coeficientes y, en algunos casos, si se desea usar el enlace canónico.
El modelo lineal normal como caso particular
El modelo lineal normal se recupera tomando:
\[ Y_i\sim N(\mu_i,\sigma^2), \]
\[ \eta_i=x_i^\top\beta, \]
y enlace identidad:
\[ g(\mu_i)=\mu_i. \]
Por tanto:
\[ \mu_i=\eta_i=x_i^\top\beta. \]
Ejemplo: respuesta de conteo
Si \(Y_i\) representa un conteo, una opción natural es:
\[ Y_i\sim Poisson(\mu_i). \]
Como \(\mu_i>0\), se puede usar enlace log:
\[ \log(\mu_i)=x_i^\top\beta. \]
Equivalente a:
\[ \mu_i=\exp(x_i^\top\beta)>0. \]
Ejemplo: respuesta binaria
Si \(Y_i\) representa éxito o fracaso:
\[ Y_i\sim Bernoulli(\mu_i), \]
donde:
\[ \mu_i=P(Y_i=1). \]
Como \(0<\mu_i<1\), se puede usar enlace logit:
\[ \log\left(\dfrac{\mu_i}{1-\mu_i}\right)=x_i^\top\beta. \]
La cantidad \(\dfrac{\mu_i}{1-\mu_i}\) representa la razón de chances.
Ejemplo: respuesta positiva continua
Si \(Y_i\) es una variable positiva y asimétrica:
\[ Y_i>0, \]
una opción posible es:
\[ Y_i\sim Gamma(\mu_i,\phi). \]
Se puede usar enlace inverso:
\[ \mu_i^{-1}=x_i^\top\beta, \]
o enlace log:
\[ \log(\mu_i)=x_i^\top\beta, \]
lo que implica:
\[ \mu_i=\exp(x_i^\top\beta)>0. \]
Estructura de elección del modelo
Para definir un GLM se deben elegir tres elementos:
\[\begin{array}{lll} 1. & \text{Distribución de }Y_i, & \text{según la naturaleza de la respuesta.}\\[1mm] 2. & \text{Predictor lineal}, & \text{según las covariables disponibles.}\\[1mm] 3. & \text{Función de enlace}, & \text{según el rango de }\mu_i. \end{array}\]Síntesis de la sección
Un GLM permite escribir, en una misma estructura:
\[ Y_i\sim FE(\mu_i,\phi), \]
\[ \eta_i=x_i^\top\beta, \]
\[ g(\mu_i)=\eta_i. \]
Así, el modelo lineal normal, la regresión de Poisson, la regresión logística y otros modelos se estudian bajo un marco común.
Estimación e inferencia
Punto de partida
Se considera:
\[ Y_i\sim FE(\mu_i,\phi), \qquad i=1,\ldots,n, \]
con función de probabilidad o densidad:
\[ f(y_i;\theta_i,\phi) = \exp\left\{ \phi[y_i\theta_i-b(\theta_i)]+c(y_i,\phi) \right\}. \]
Además:
\[ E(Y_i)=\mu_i=b'(\theta_i) \]
y
\[ Var(Y_i)=\phi^{-1}V(\mu_i). \]
El predictor lineal está dado por:
\[ \eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}. \]
Verosimilitud
Para respuestas independientes, el logaritmo de la función de verosimilitud es:
\[ L(\boldsymbol{\beta},\phi) = \sum_{i=1}^n \left[ \phi\{y_i\theta_i-b(\theta_i)\} + c(y_i,\phi) \right]. \]
La relación con \(\boldsymbol{\beta}\) se obtiene a través de:
\[ \eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}, \]
\[ \mu_i=g^{-1}(\eta_i) = g^{-1}(\boldsymbol{x}_i^\top\boldsymbol{\beta}), \]
y
\[ \mu_i=b'(\theta_i). \]
Por tanto, al cambiar \(\boldsymbol{\beta}\) cambia \(\eta_i\), cambia \(\mu_i\) y también cambia el parámetro canónico \(\theta_i\) asociado a esa media.
Así, la verosimilitud depende de \(\boldsymbol{\beta}\) mediante la secuencia:
\[ \boldsymbol{\beta} \longrightarrow \eta_i \longrightarrow \mu_i \longrightarrow \theta_i \longrightarrow L(\boldsymbol{\beta},\phi). \]
Enlace canónico
La función de enlace relaciona la media de cada observación con el predictor lineal:
\[ g(\mu_i)=\eta_i, \qquad i=1,\ldots,n. \]
Se dice que el enlace es canónico cuando el predictor lineal coincide con el parámetro canónico:
\[ \theta_i=\eta_i. \]
Como:
\[ \eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}, \]
bajo enlace canónico se tiene:
\[ \theta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}. \]
Verosimilitud con enlace canónico
La log-verosimilitud general es:
\[ L(\boldsymbol{\beta},\phi) = \sum_{i=1}^n \left[ \phi\{y_i\theta_i-b(\theta_i)\} + c(y_i,\phi) \right]. \]
Si el enlace es canónico:
\[ \theta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}. \]
Entonces:
\[ L(\boldsymbol{\beta},\phi) = \phi \sum_{i=1}^n \left[ y_i\boldsymbol{x}_i^\top\boldsymbol{\beta} - b(\boldsymbol{x}_i^\top\boldsymbol{\beta}) \right] + \sum_{i=1}^n c(y_i,\phi). \]
Forma equivalente de la verosimilitud
Como:
\[ \boldsymbol{x}_i^\top\boldsymbol{\beta} = \sum_{j=1}^p x_{ij}\beta_j, \]
se tiene:
\[ \sum_{i=1}^n y_i\boldsymbol{x}_i^\top\boldsymbol{\beta} = \sum_{j=1}^p \beta_j \sum_{i=1}^n y_i x_{ij}. \]
Definiendo:
\[ S_j=\sum_{i=1}^n y_i x_{ij}, \qquad j=1,\ldots,p, \]
la log-verosimilitud queda:
\[ L(\boldsymbol{\beta},\phi) = \phi \left[ \sum_{j=1}^p \beta_j S_j - \sum_{i=1}^n b(\boldsymbol{x}_i^\top\boldsymbol{\beta}) \right] + \sum_{i=1}^n c(y_i,\phi). \]
Enlaces canónicos principales
| Distribución | Parámetro canónico | Enlace canónico |
|---|---|---|
| Normal | \(\theta=\mu\) | \(\mu_i=\eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}\) |
| Poisson | \(\theta=\log(\mu)\) | \(\log(\mu_i)=\eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}\) |
| Binomial | \(\theta=\log\left(\dfrac{\mu}{1-\mu}\right)\) | \(\log\left(\dfrac{\mu_i}{1-\mu_i}\right)=\eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}\) |
| Gamma | \(\theta=-\dfrac{1}{\mu}\) | \(\mu_i^{-1}=\eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}\) |
| Normal inversa | \(\theta=-\dfrac{1}{2\mu^2}\) | \(\mu_i^{-2}=\eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}\) |
Ventaja del enlace canónico
Una ventaja de usar enlaces canónicos es que la log-verosimilitud adopta una forma más simple:
\[ L(\boldsymbol{\beta},\phi) = \phi \left[ \sum_{j=1}^p \beta_j S_j - \sum_{i=1}^n b(\boldsymbol{x}_i^\top\boldsymbol{\beta}) \right] + \sum_{i=1}^n c(y_i,\phi). \]
Esta forma facilita:
- la derivación de la función score;
- la obtención de la información de Fisher;
- el estudio de la concavidad de \(L(\boldsymbol{\beta},\phi)\);
- el análisis de la existencia del estimador de máxima verosimilitud.
Por eso el enlace canónico no es solo una elección de escala: también simplifica la inferencia.
Función de desvío
Considere un GLM con matriz de modelo:
\[ \boldsymbol{X}_{n\times p}, \qquad p<n. \]
En este modelo, las medias no son libres, sino que están determinadas por:
\[ \mu_i=g^{-1}(\eta_i), \qquad \eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}. \]
Por tanto, el vector de medias:
\[ \boldsymbol{\mu}=(\mu_1,\ldots,\mu_n)^\top \]
queda restringido por la estructura del modelo.
La función de desvío compara este modelo con el modelo saturado.
Modelo saturado
El modelo saturado es aquel que asigna una media propia a cada observación.
En este caso:
\[ \tilde{\mu}_i=y_i, \qquad i=1,\ldots,n. \]
Por tanto, el modelo saturado reproduce exactamente los datos observados.
Como se requiere una media para cada observación, este modelo usa típicamente \(n\) parámetros para \(n\) datos.
El modelo saturado representa el ajuste más complejo posible, pero no necesariamente aporta una explicación del fenómeno: reproduce los datos, pero no resume una estructura más simple.
Si se escribe la log-verosimilitud como función de las medias:
\[ L(\boldsymbol{\mu};\boldsymbol{y}) = \sum_{i=1}^n L(\mu_i;y_i), \]
entonces, para el modelo saturado:
\[ L(\boldsymbol{y};\boldsymbol{y}) = \sum_{i=1}^n L(y_i;y_i). \]
Modelo nulo
El modelo nulo es el modelo más simple dentro de la clase considerada.
En este modelo todas las observaciones tienen la misma media:
\[ \mu_i=\mu, \qquad i=1,\ldots,n. \]
Por tanto, el predictor lineal contiene solo una constante:
\[ \eta_i=\beta_0. \]
El modelo nulo no incorpora covariables explicativas.
En ese sentido, representa la situación en la que no se modela una relación entre los predictores y la variable respuesta.
Así, el modelo propuesto se ubica entre dos referencias:
\[ \text{modelo nulo} \quad \longrightarrow \quad \text{modelo propuesto} \quad \longrightarrow \quad \text{modelo saturado}. \]
El modelo nulo concentra la explicación en una media común; el modelo saturado reproduce cada observación; el modelo propuesto busca una estructura intermedia y más informativa.
Log-verosimilitud del modelo propuesto
Para el modelo propuesto:
\[ \mu_i=g^{-1}(\boldsymbol{x}_i^\top\boldsymbol{\beta}). \]
Entonces:
\[ L(\boldsymbol{\mu};\boldsymbol{y}) = \sum_{i=1}^n L(\mu_i;y_i). \]
Una vez estimado el modelo, se obtiene:
\[ \hat{\mu}_i = g^{-1}(\boldsymbol{x}_i^\top\hat{\boldsymbol{\beta}}). \]
Por tanto, la log-verosimilitud maximizada bajo el modelo propuesto es:
\[ L(\hat{\boldsymbol{\mu}};\boldsymbol{y}) = \sum_{i=1}^n L(\hat{\mu}_i;y_i). \]
El modelo propuesto separa los datos en una parte sistemática, representada por \(\boldsymbol{x}_i^\top\boldsymbol{\beta}\), y una parte aleatoria, representada por la distribución de \(Y_i\) alrededor de \(\mu_i\).
Log-verosimilitudes de referencia
Para el modelo saturado:
\[ \tilde{\mu}_i=y_i. \]
Por tanto:
\[ L(\boldsymbol{y};\boldsymbol{y}) = \sum_{i=1}^n L(y_i;y_i). \]
Para el modelo nulo:
\[ \mu_i=\hat{\mu}_0. \]
Por tanto:
\[ L(\hat{\boldsymbol{\mu}}_0;\boldsymbol{y}) = \sum_{i=1}^n L(\hat{\mu}_0;y_i). \]
Estas tres log-verosimilitudes permiten comparar el ajuste del modelo propuesto con el modelo nulo y con el modelo saturado.
El modelo saturado funciona como referencia superior de ajuste, mientras que el modelo nulo funciona como referencia básica.
Función de desvío
La función de desvío, también denominada deviance escalada, compara la log-verosimilitud del modelo saturado con la log-verosimilitud de un modelo con medias \(\boldsymbol{\mu}\).
Se define como:
\[ D^*(\boldsymbol{y};\boldsymbol{\mu}) = 2 \left[ L(\boldsymbol{y};\boldsymbol{y}) - L(\boldsymbol{\mu};\boldsymbol{y}) \right]. \]
Cuando se evalúa en el modelo ajustado:
\[ D^*(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = 2 \left[ L(\boldsymbol{y};\boldsymbol{y}) - L(\hat{\boldsymbol{\mu}};\boldsymbol{y}) \right]. \]
Esta cantidad mide la pérdida de ajuste al reemplazar el modelo saturado por el modelo propuesto.
Desvío
Como \(\phi\) es el parámetro de precisión:
\[ D^*(\boldsymbol{y};\boldsymbol{\mu}) = \phi D(\boldsymbol{y};\boldsymbol{\mu}). \]
Por tanto:
\[ D(\boldsymbol{y};\boldsymbol{\mu}) = \frac{1}{\phi} D^*(\boldsymbol{y};\boldsymbol{\mu}). \]
La cantidad:
\[ D(\boldsymbol{y};\boldsymbol{\mu}) \]
se denomina desvío, también llamado deviance.
Al evaluar en el modelo ajustado se obtiene:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}). \]
Forma alternativa del desvío
Sea \(\tilde{\theta}_i\) el parámetro canónico asociado al modelo saturado:
\[ \tilde{\mu}_i=y_i. \]
Sea \(\theta_i\) el parámetro canónico asociado a una media \(\mu_i\).
Entonces:
\[ D(\boldsymbol{y};\boldsymbol{\mu}) = 2 \sum_{i=1}^n \left[ y_i(\tilde{\theta}_i-\theta_i) - b(\tilde{\theta}_i) + b(\theta_i) \right]. \]
Equivalente:
\[ D(\boldsymbol{y};\boldsymbol{\mu}) = 2 \sum_{i=1}^n \left[ y_i(\tilde{\theta}_i-\theta_i) + b(\theta_i) - b(\tilde{\theta}_i) \right]. \]
Componentes del desvío
El desvío puede escribirse como suma de componentes:
\[ D(\boldsymbol{y};\boldsymbol{\mu}) = \sum_{i=1}^n d_i^2(y_i;\mu_i). \]
Cada componente:
\[ d_i^2(y_i;\mu_i) \]
mide la contribución de la observación \(i\) al desvío total.
Al evaluar en el modelo ajustado:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = \sum_{i=1}^n d_i^2(y_i;\hat{\mu}_i). \]
Esta descomposición será importante para definir residuos de devianza.
Interpretación del desvío
El desvío no mide una distancia ordinaria entre \(y_i\) y \(\mu_i\).
Mide una diferencia de log-verosimilitudes:
\[ L(\boldsymbol{y};\boldsymbol{y}) - L(\boldsymbol{\mu};\boldsymbol{y}). \]
Por tanto, cuantifica cuánto ajuste se pierde al pasar del modelo saturado a un modelo con medias restringidas por:
\[ \mu_i=g^{-1}(\boldsymbol{x}_i^\top\boldsymbol{\beta}). \]
Si el modelo ajusta bien, entonces:
\[ L(\hat{\boldsymbol{\mu}};\boldsymbol{y}) \]
será cercano a:
\[ L(\boldsymbol{y};\boldsymbol{y}). \]
En este sentido, el modelo bajo estudio busca explicar parte de la estructura sistemática de los datos sin llegar al extremo no informativo del modelo saturado.
Nomenclatura
En esta presentación se usa la parametrización con \(\phi\) como parámetro de precisión:
\[ Var(Y)=\phi^{-1}V(\mu). \]
Por ello:
\[ D^*(\boldsymbol{y};\boldsymbol{\mu}) = \phi D(\boldsymbol{y};\boldsymbol{\mu}). \]
Entonces:
| Notación | Nombre usado |
|---|---|
| \(D^*(\boldsymbol{y};\boldsymbol{\mu})\) | función de desvío |
| \(D(\boldsymbol{y};\boldsymbol{\mu})\) | desvío |
| \(D^*(\boldsymbol{y};\hat{\boldsymbol{\mu}})\) | función de desvío del modelo ajustado |
| \(D(\boldsymbol{y};\hat{\boldsymbol{\mu}})\) | desvío del modelo ajustado |
En otros textos, donde \(\phi\) se usa como parámetro de dispersión, la relación puede escribirse con una división por \(\phi\).
Desvío: caso normal (I)
Para el caso normal:
\[ \theta=\mu, \qquad b(\theta)=\frac{\theta^2}{2}. \]
La forma general del desvío es:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = 2 \sum_{i=1}^n \left[ y_i(\tilde{\theta}_i-\hat{\theta}_i) - b(\tilde{\theta}_i) + b(\hat{\theta}_i) \right]. \]
En el modelo saturado:
\[ \tilde{\mu}_i=y_i \quad \Rightarrow \quad \tilde{\theta}_i=y_i. \]
En el modelo ajustado:
\[ \hat{\mu}_i = g^{-1}(\boldsymbol{x}_i^\top\hat{\boldsymbol{\beta}}) \quad \Rightarrow \quad \hat{\theta}_i=\hat{\mu}_i. \]
Desvío: caso normal (II)
Sustituyendo en la forma general:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = 2 \sum_{i=1}^n \left[ y_i(y_i-\hat{\mu}_i) - \frac{y_i^2}{2} + \frac{\hat{\mu}_i^2}{2} \right]. \]
Como:
\[ y_i(y_i-\hat{\mu}_i) = y_i^2-y_i\hat{\mu}_i, \]
entonces:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = 2 \sum_{i=1}^n \left[ y_i^2-y_i\hat{\mu}_i - \frac{y_i^2}{2} + \frac{\hat{\mu}_i^2}{2} \right]. \]
Desvío: caso normal (III)
Agrupando términos:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = 2 \sum_{i=1}^n \left[ \frac{y_i^2}{2} - y_i\hat{\mu}_i + \frac{\hat{\mu}_i^2}{2} \right]. \]
Por tanto:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = \sum_{i=1}^n \left[ y_i^2 - 2y_i\hat{\mu}_i + \hat{\mu}_i^2 \right]. \]
Finalmente:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = \sum_{i=1}^n (y_i-\hat{\mu}_i)^2. \]
que coincide con la suma de cuadrados de residuos.
Desvío: caso Poisson
Para el caso Poisson:
\[ \theta=\log(\mu), \qquad b(\theta)=e^\theta=\mu. \]
En el modelo ajustado:
\[ \hat{\theta}_i=\log(\hat{\mu}_i). \]
En el modelo saturado, si \(y_i>0\):
\[ \tilde{\theta}_i=\log(y_i). \]
Entonces:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = 2 \sum_{i=1}^n \left[ y_i\log\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i-\hat{\mu}_i) \right]. \]
Si \(y_i=0\), el término correspondiente se define por continuidad y queda:
\[ d_i^2(0;\hat{\mu}_i)=2\hat{\mu}_i. \]
Por tanto, para el modelo Poisson:
\[ d_i^2(y_i;\hat{\mu}_i) = \begin{cases} 2\left[ y_i\log\left(\dfrac{y_i}{\hat{\mu}_i}\right) - (y_i-\hat{\mu}_i) \right], & y_i>0,\\[6pt] 2\hat{\mu}_i, & y_i=0. \end{cases} \]
Desvío: caso Gamma
Para el caso Gamma:
\[ \theta=-\frac{1}{\mu}, \qquad b(\theta)=-\log(-\theta). \]
En el modelo ajustado:
\[ \hat{\theta}_i=-\frac{1}{\hat{\mu}_i}. \]
En el modelo saturado:
\[ \tilde{\theta}_i=-\frac{1}{y_i}. \]
Cuando \(y_i>0\) y \(\hat{\mu}_i>0\):
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = 2 \sum_{i=1}^n \left[ -\log\left(\frac{y_i}{\hat{\mu}_i}\right) + \frac{y_i-\hat{\mu}_i}{\hat{\mu}_i} \right]. \]
Papel de la función de desvío
La función de desvío se define antes de estimar el modelo porque nace de la comparación entre dos log-verosimilitudes:
\[ L(\boldsymbol{y};\boldsymbol{y}) \qquad \text{y} \qquad L(\boldsymbol{\mu};\boldsymbol{y}). \]
Sin embargo, la estimación de \(\boldsymbol{\beta}\) no se obtiene minimizando directamente el desvío.
La estimación se obtiene maximizando la log-verosimilitud, o equivalentemente resolviendo la función score.
Después de estimar el modelo, el desvío se evalúa en:
\[ \hat{\boldsymbol{\mu}}. \]
Ejercicio: componentes del desvío para el modelo binomial
Suponga que:
\[ Y_i \sim Binomial(m_i,\pi_i), \qquad i=1,\ldots,n. \]
Sea:
\[ \mu_i=E(Y_i)=m_i\pi_i. \]
Obtenga el componente individual del desvío:
\[ d_i^2(y_i;\hat{\mu}_i), \]
donde \(\hat{\mu}_i=m_i\hat{\pi}_i\).
Particularice el resultado para el caso:
\[ m_i=1, \]
es decir, para una variable Bernoulli.
Solución: forma binomial
La función de probabilidad binomial es:
\[ P(Y_i=y_i) = \binom{m_i}{y_i} \pi_i^{y_i} (1-\pi_i)^{m_i-y_i}. \]
Como:
\[ \mu_i=m_i\pi_i, \]
se tiene:
\[ \pi_i=\frac{\mu_i}{m_i}. \]
Entonces:
\[ P(Y_i=y_i) = \binom{m_i}{y_i} \left(\frac{\mu_i}{m_i}\right)^{y_i} \left(1-\frac{\mu_i}{m_i}\right)^{m_i-y_i}. \]
Solución: parámetro canónico
Para la distribución binomial, el parámetro canónico es:
\[ \theta_i = \log\left( \frac{\pi_i}{1-\pi_i} \right). \]
Como:
\[ \pi_i=\frac{\mu_i}{m_i}, \]
entonces:
\[ \theta_i = \log\left( \frac{\mu_i/m_i}{1-\mu_i/m_i} \right). \]
Por tanto:
\[ \theta_i = \log\left( \frac{\mu_i}{m_i-\mu_i} \right). \]
Además:
\[ b(\theta_i) = m_i\log(1+e^{\theta_i}). \]
Solución: modelo saturado y ajustado
En el modelo saturado:
\[ \tilde{\mu}_i=y_i. \]
Como la respuesta es un conteo binomial:
\[ 0\leq y_i\leq m_i. \]
El parámetro canónico para una media \(\mu_i\) es:
\[ \theta_i = \log\left( \frac{\mu_i}{m_i-\mu_i} \right). \]
Por tanto, en el modelo saturado:
\[ \tilde{\theta}_i = \log\left( \frac{y_i}{m_i-y_i} \right). \]
En el modelo ajustado:
\[ \hat{\theta}_i = \log\left( \frac{\hat{\mu}_i}{m_i-\hat{\mu}_i} \right). \]
Solución: punto de partida
La forma general del desvío es:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = 2 \sum_{i=1}^n \left[ y_i(\tilde{\theta}_i-\hat{\theta}_i) - b(\tilde{\theta}_i) + b(\hat{\theta}_i) \right]. \]
Como:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = \sum_{i=1}^n d_i^2(y_i;\hat{\mu}_i), \]
el componente individual es:
\[ d_i^2(y_i;\hat{\mu}_i) = 2 \left[ y_i(\tilde{\theta}_i-\hat{\theta}_i) - b(\tilde{\theta}_i) + b(\hat{\theta}_i) \right]. \]
Solución: término con \(\theta_i\)
Primero se calcula:
\[ \tilde{\theta}_i-\hat{\theta}_i = \log\left( \frac{y_i}{m_i-y_i} \right) - \log\left( \frac{\hat{\mu}_i}{m_i-\hat{\mu}_i} \right). \]
Usando propiedades de logaritmos:
\[ \tilde{\theta}_i-\hat{\theta}_i = \log\left[ \frac{y_i}{m_i-y_i} \cdot \frac{m_i-\hat{\mu}_i}{\hat{\mu}_i} \right]. \]
Entonces:
\[ \tilde{\theta}_i-\hat{\theta}_i = \log\left( \frac{y_i(m_i-\hat{\mu}_i)} {\hat{\mu}_i(m_i-y_i)} \right). \]
Por tanto:
\[ y_i(\tilde{\theta}_i-\hat{\theta}_i) = y_i \log\left( \frac{y_i(m_i-\hat{\mu}_i)} {\hat{\mu}_i(m_i-y_i)} \right). \]
Solución: función \(b(\theta_i)\)
Para el conteo binomial:
\[ b(\theta_i) = m_i\log(1+e^{\theta_i}). \]
Como:
\[ e^{\theta_i} = \frac{\mu_i}{m_i-\mu_i}, \]
entonces:
\[ 1+e^{\theta_i} = 1+ \frac{\mu_i}{m_i-\mu_i}. \]
Por tanto:
\[ 1+e^{\theta_i} = \frac{m_i-\mu_i+\mu_i}{m_i-\mu_i} = \frac{m_i}{m_i-\mu_i}. \]
Así:
\[ b(\theta_i) = m_i \log\left( \frac{m_i}{m_i-\mu_i} \right). \]
Solución: \(b(\tilde{\theta}_i)\) y \(b(\hat{\theta}_i)\)
En el modelo saturado:
\[ \tilde{\mu}_i=y_i. \]
Entonces:
\[ b(\tilde{\theta}_i) = m_i \log\left( \frac{m_i}{m_i-y_i} \right). \]
En el modelo ajustado:
\[ \hat{\mu}_i. \]
Entonces:
\[ b(\hat{\theta}_i) = m_i \log\left( \frac{m_i}{m_i-\hat{\mu}_i} \right). \]
Por tanto:
\[ -b(\tilde{\theta}_i)+b(\hat{\theta}_i) = -m_i \log\left( \frac{m_i}{m_i-y_i} \right) + m_i \log\left( \frac{m_i}{m_i-\hat{\mu}_i} \right). \]
Solución: simplificación de \(b(\theta_i)\)
Agrupando logaritmos:
\[ -b(\tilde{\theta}_i)+b(\hat{\theta}_i) = m_i \log\left[ \frac{ m_i/(m_i-\hat{\mu}_i) }{ m_i/(m_i-y_i) } \right]. \]
Entonces:
\[ -b(\tilde{\theta}_i)+b(\hat{\theta}_i) = m_i \log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right). \]
Así, el componente individual queda:
\[ d_i^2(y_i;\hat{\mu}_i) = 2 \left[ y_i \log\left( \frac{y_i(m_i-\hat{\mu}_i)} {\hat{\mu}_i(m_i-y_i)} \right) + m_i \log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right) \right]. \]
Solución: separación de términos
Usando:
\[ \log\left( \frac{y_i(m_i-\hat{\mu}_i)} {\hat{\mu}_i(m_i-y_i)} \right) = \log\left(\frac{y_i}{\hat{\mu}_i}\right) + \log\left( \frac{m_i-\hat{\mu}_i}{m_i-y_i} \right), \]
se obtiene:
\[ d_i^2(y_i;\hat{\mu}_i) = 2 \left[ y_i\log\left(\frac{y_i}{\hat{\mu}_i}\right) + y_i\log\left( \frac{m_i-\hat{\mu}_i}{m_i-y_i} \right) + m_i\log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right) \right]. \]
Solución: agrupación final
Los dos últimos términos son:
\[ y_i\log\left( \frac{m_i-\hat{\mu}_i}{m_i-y_i} \right) + m_i\log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right). \]
Como:
\[ \log\left( \frac{m_i-\hat{\mu}_i}{m_i-y_i} \right) = -\log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right), \]
entonces:
\[ -y_i\log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right) + m_i\log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right) \]
queda como:
\[ (m_i-y_i) \log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right). \]
Solución: componente binomial
Por tanto:
\[ d_i^2(y_i;\hat{\mu}_i) = 2 \left[ y_i\log\left(\frac{y_i}{\hat{\mu}_i}\right) + (m_i-y_i) \log\left( \frac{m_i-y_i}{m_i-\hat{\mu}_i} \right) \right]. \]
Esta es la contribución de la observación \(i\) al desvío binomial cuando la respuesta es el conteo de éxitos.
Usando la convención usual:
\[ 0\log(0)=0, \]
esta expresión también cubre los casos de borde \(y_i=0\) y \(y_i=m_i\).
Caso particular: Bernoulli
Si:
\[ m_i=1, \]
entonces:
\[ Y_i\sim Bernoulli(\pi_i), \qquad \mu_i=\pi_i. \]
El componente del desvío queda:
\[ d_i^2(y_i;\hat{\mu}_i) = 2 \left[ y_i\log\left(\frac{y_i}{\hat{\mu}_i}\right) + (1-y_i) \log\left( \frac{1-y_i}{1-\hat{\mu}_i} \right) \right]. \]
Caso particular: Bernoulli por valores
Si \(y_i=1\):
\[ d_i^2(1;\hat{\mu}_i) = 2\log\left(\frac{1}{\hat{\mu}_i}\right). \]
Si \(y_i=0\):
\[ d_i^2(0;\hat{\mu}_i) = 2\log\left(\frac{1}{1-\hat{\mu}_i}\right). \]
Por tanto:
\[ d_i^2(y_i;\hat{\mu}_i) = \begin{cases} 2\log\left(\dfrac{1}{\hat{\mu}_i}\right), & y_i=1,\\[8pt] 2\log\left(\dfrac{1}{1-\hat{\mu}_i}\right), & y_i=0. \end{cases} \]
Resultados asintóticos del desvío I
Aunque es usual comparar el desvío del modelo ajustado con cuantiles de una distribución chi-cuadrado, en general:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \not\sim \chi^2_{n-p}. \]
Por tanto, la aproximación:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \approx \chi^2_{n-p} \]
no debe usarse automáticamente.
Su validez depende de la distribución considerada y del régimen asintótico bajo el cual se trabaja.
Resultados asintóticos del desvío II
En algunos casos particulares sí se obtiene una aproximación chi-cuadrado.
Para el caso binomial, si el número de observaciones permanece fijo y los tamaños binomiales crecen:
\[ m_i\to\infty, \qquad i=1,\ldots,n, \]
entonces, bajo el modelo verdadero:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \sim \chi^2_{n-p}. \]
Para el caso Poisson, si las medias crecen:
\[ \mu_i\to\infty, \qquad i=1,\ldots,n, \]
entonces:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \sim \chi^2_{n-p}. \]
Resultados asintóticos del desvío III
Para el caso normal, con \(\sigma^2\) fijo:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \sim \sigma^2\chi^2_{n-p}. \]
Equivalente:
\[ \frac{ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) }{ \sigma^2 } \sim \chi^2_{n-p}. \]
Como, en el modelo normal:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = \sum_{i=1}^n (y_i-\hat{\mu}_i)^2, \]
se obtiene:
\[ \frac{ \sum_{i=1}^n (y_i-\hat{\mu}_i)^2 }{ \sigma^2 } \sim \chi^2_{n-p}. \]
Resultados asintóticos del desvío IV
En general, cuando la función de desvío depende del parámetro de precisión \(\phi\), se considera:
\[ D^*(\boldsymbol{y};\hat{\boldsymbol{\mu}}) = \phi D(\boldsymbol{y};\hat{\boldsymbol{\mu}}). \]
Bajo condiciones regulares:
\[ D^*(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \sim \chi^2_{n-p}, \qquad \phi\to\infty. \]
Como \(\phi\) es un parámetro de precisión, la condición:
\[ \phi\to\infty \]
equivale a considerar dispersión pequeña.
Por tanto, cuando la dispersión es pequeña, puede ser razonable comparar la función de desvío con cuantiles de:
\[ \chi^2_{n-p}. \]
Resultados asintóticos del desvío V
En el caso Gamma, el desvío se aproxima a una distribución chi-cuadrado con \(n-p\) grados de libertad cuando el coeficiente de variación es cercano a cero.
En síntesis:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \sim \chi^2_{n-p} \]
no es un resultado general.
La aproximación debe justificarse según:
- la distribución de la respuesta;
- el tamaño de los parámetros involucrados;
- el comportamiento de la dispersión;
- el régimen asintótico considerado.
Análisis del desvío
El análisis del desvío permite comparar modelos GLM anidados.
Se considera una partición del vector de parámetros:
\[ \boldsymbol{\beta} = (\boldsymbol{\beta}_1^\top,\boldsymbol{\beta}_2^\top)^\top, \]
donde \(\boldsymbol{\beta}_1\) tiene dimensión \(q\) y \(\boldsymbol{\beta}_2\) tiene dimensión \(p-q\).
El interés es probar:
\[ H_0:\boldsymbol{\beta}_1=\boldsymbol{0} \]
contra:
\[ H_1:\boldsymbol{\beta}_1\neq\boldsymbol{0}. \]
Modelos bajo \(H_0\) y bajo \(H_1\)
Bajo \(H_0\), se impone:
\[ \boldsymbol{\beta}_1=\boldsymbol{0}. \]
Por tanto, se obtiene un modelo reducido.
Sea:
\[ \hat{\boldsymbol{\mu}}^0 \]
el estimador de máxima verosimilitud de \(\boldsymbol{\mu}\) bajo \(H_0\).
Bajo \(H_1\), se ajusta el modelo completo.
Sea:
\[ \hat{\boldsymbol{\mu}} \]
el estimador de máxima verosimilitud de \(\boldsymbol{\mu}\) bajo \(H_1\).
Desvíos bajo \(H_0\) y bajo \(H_1\)
El desvío bajo el modelo reducido es:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0). \]
El desvío bajo el modelo completo es:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}). \]
Como el modelo completo contiene al modelo reducido, su ajuste no puede ser peor.
Por tanto:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \leq D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0). \]
Reducción del desvío
La reducción del desvío al pasar del modelo reducido al modelo completo es:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}). \]
Si esta reducción es grande, los términos asociados a \(\boldsymbol{\beta}_1\) mejoran de manera importante el ajuste del modelo.
Si esta reducción es pequeña, el modelo reducido puede ser suficiente.
La reducción del desvío se evalúa considerando el número de restricciones:
\[ q. \]
Estadística de razón de verosimilitud
La estadística de razón de verosimilitud queda dada por:
\[ \xi_{RV} = \phi \left[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \right]. \]
Bajo \(H_0\):
\[ \xi_{RV} \sim \chi^2_q, \qquad n\to\infty. \]
Por tanto, valores grandes de \(\xi_{RV}\) indican evidencia contra \(H_0\).
En ese caso, los términos asociados a \(\boldsymbol{\beta}_1\) contribuyen al modelo.
Caso con \(\phi=1\)
En modelos como Poisson y binomial, usualmente:
\[ \phi=1. \]
Entonces:
\[ \xi_{RV} = D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}). \]
Bajo \(H_0\):
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \sim \chi^2_q. \]
Así, en estos casos la comparación se basa directamente en la reducción del desvío.
Estadística tipo \(F\)
De forma similar, se puede definir la estadística:
\[ F = \frac{ \left[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \right]/q }{ D(\boldsymbol{y};\hat{\boldsymbol{\mu}})/(n-p) }. \]
Bajo \(H_0\):
\[ F \sim F_{q,n-p}. \]
Esta forma es útil porque no depende directamente del parámetro de precisión \(\phi\).
Lectura de la estadística tipo \(F\)
El numerador mide la reducción media del desvío al pasar del modelo reducido al modelo completo:
\[ \frac{ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) }{ q }. \]
El denominador mide el desvío residual medio del modelo completo:
\[ \frac{ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) }{ n-p }. \]
Por tanto, la estadística \(F\) compara la mejora obtenida con el desvío residual que permanece en el modelo completo.
Tabla de análisis del desvío
El análisis del desvío puede organizarse en una tabla:
| Modelo | gl residuales | Desvío |
|---|---|---|
| Bajo \(H_0\) | \(n-(p-q)\) | \(D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0)\) |
| Bajo \(H_1\) | \(n-p\) | \(D(\boldsymbol{y};\hat{\boldsymbol{\mu}})\) |
La diferencia de desvíos es:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}). \]
La diferencia de grados de libertad es:
\[ q. \]
Interpretación del análisis del desvío
El análisis del desvío permite evaluar si un conjunto de términos debe permanecer en el modelo.
Si el valor-p es pequeño, se rechaza:
\[ H_0:\boldsymbol{\beta}_1=\boldsymbol{0}. \]
En ese caso, los términos asociados a \(\boldsymbol{\beta}_1\) aportan información para explicar la respuesta.
Si el valor-p no es pequeño, el modelo reducido puede ser suficiente.
Análisis del desvío en R
En R, la comparación entre modelos GLM anidados puede realizarse con:
Función score para \(\boldsymbol{\beta}\) I
Sea el vector de parámetros del modelo:
\[ \boldsymbol{\theta} = (\boldsymbol{\beta}^\top,\phi)^\top. \]
El logaritmo de la verosimilitud es:
\[ L(\boldsymbol{\theta}) = \sum_{i=1}^n \left[ \phi\{y_i\theta_i-b(\theta_i)\} + c(y_i,\phi) \right]. \]
Para estimar \(\boldsymbol{\beta}\) se estudia la función score:
\[ \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\theta}) = \frac{\partial L(\boldsymbol{\theta})} {\partial \boldsymbol{\beta}}. \]
Como:
\[ \boldsymbol{\beta} \longrightarrow \eta_i \longrightarrow \mu_i \longrightarrow \theta_i \longrightarrow L_i, \]
la derivada se obtiene usando la regla de la cadena.
Función score para \(\boldsymbol{\beta}\) II
Para cada parámetro \(\beta_j\):
\[ \frac{\partial L(\boldsymbol{\theta})}{\partial \beta_j} = \sum_{i=1}^n \phi \left[ y_i \frac{d\theta_i}{d\mu_i} \frac{d\mu_i}{d\eta_i} \frac{\partial \eta_i}{\partial \beta_j} - \frac{db(\theta_i)}{d\theta_i} \frac{d\theta_i}{d\mu_i} \frac{d\mu_i}{d\eta_i} \frac{\partial \eta_i}{\partial \beta_j} \right]. \]
Como:
\[ \frac{db(\theta_i)}{d\theta_i} = \mu_i, \]
entonces:
\[ \frac{\partial L(\boldsymbol{\theta})}{\partial \beta_j} = \sum_{i=1}^n \phi (y_i-\mu_i) \frac{d\theta_i}{d\mu_i} \frac{d\mu_i}{d\eta_i} x_{ij}. \]
Función score para \(\boldsymbol{\beta}\) III
Además:
\[ V_i=V(\mu_i)=\frac{d\mu_i}{d\theta_i}. \]
Por tanto:
\[ \frac{d\theta_i}{d\mu_i} = \frac{1}{V_i}. \]
Así:
\[ \frac{\partial L(\boldsymbol{\theta})}{\partial \beta_j} = \phi \sum_{i=1}^n \frac{1}{V_i} \left( \frac{d\mu_i}{d\eta_i} \right) (y_i-\mu_i)x_{ij}. \]
Esta expresión muestra que el score depende de:
- los residuos \(y_i-\mu_i\);
- la función de varianza \(V(\mu_i)\);
- la función de enlace mediante \(d\mu_i/d\eta_i\).
Matriz de pesos
Se define:
\[ w_i = \frac{ \left( \dfrac{d\mu_i}{d\eta_i} \right)^2 }{ V_i }. \]
Entonces:
\[ \boldsymbol{W} = diag(w_1,\ldots,w_n). \]
También se define:
\[ \boldsymbol{V} = diag(V_1,\ldots,V_n), \]
donde:
\[ V_i=V(\mu_i). \]
La matriz \(\boldsymbol{W}\) contiene los pesos que cambian según la media y según la función de enlace.
Ejercicio guiado: Poisson con enlace log I
Suponga que:
\[ Y_i\sim Poisson(\lambda_i), \qquad i=1,\ldots,n. \]
Considere el modelo:
\[ \log(\lambda_i)=\beta_1+\beta_2x_i. \]
Determine:
- el vector de medias \(\boldsymbol{\mu}\);
- la matriz de diseño \(\boldsymbol{X}\);
- la función de varianza \(V(\mu_i)\);
- la matriz de pesos \(\boldsymbol{W}\).
Solución: media de la respuesta
Como:
\[ Y_i\sim Poisson(\lambda_i), \]
entonces:
\[ E(Y_i)=\lambda_i. \]
En un GLM se denota:
\[ \mu_i=E(Y_i). \]
Por tanto:
\[ \mu_i=\lambda_i. \]
Así, el modelo dado:
\[ \log(\lambda_i)=\beta_1+\beta_2x_i \]
puede escribirse como:
\[ \log(\mu_i)=\beta_1+\beta_2x_i. \]
Solución: predictor lineal
En un GLM:
\[ g(\mu_i)=\eta_i. \]
Para el modelo Poisson con enlace log:
\[ g(\mu_i)=\log(\mu_i). \]
Entonces:
\[ \eta_i=\log(\mu_i). \]
Como el modelo propone:
\[ \log(\mu_i)=\beta_1+\beta_2x_i, \]
se obtiene:
\[ \eta_i=\beta_1+\beta_2x_i. \]
Por tanto, el predictor lineal es:
\[ \eta_i=\beta_1+\beta_2x_i. \]
Solución: media en función de los parámetros
A partir de:
\[ \log(\mu_i)=\beta_1+\beta_2x_i, \]
aplicamos exponencial en ambos lados:
\[ \exp\{\log(\mu_i)\} = \exp(\beta_1+\beta_2x_i). \]
Como:
\[ \exp\{\log(\mu_i)\}=\mu_i, \]
entonces:
\[ \mu_i=\exp(\beta_1+\beta_2x_i). \]
Por tanto, para cada observación:
\[ \mu_i=e^{\beta_1+\beta_2x_i}. \]
Solución: vector de medias
El vector de medias es:
\[ \boldsymbol{\mu} = \begin{pmatrix} \mu_1\\ \mu_2\\ \vdots\\ \mu_n \end{pmatrix}. \]
Como:
\[ \mu_i=\exp(\beta_1+\beta_2x_i), \]
se obtiene:
\[ \boldsymbol{\mu} = \begin{pmatrix} \exp(\beta_1+\beta_2x_1)\\ \exp(\beta_1+\beta_2x_2)\\ \vdots\\ \exp(\beta_1+\beta_2x_n) \end{pmatrix}. \]
También puede escribirse como:
\[ \boldsymbol{\mu} = \exp(\boldsymbol{\eta}), \]
donde la exponencial se aplica componente a componente.
Solución: matriz de diseño I
El predictor lineal es:
\[ \eta_i=\beta_1+\beta_2x_i. \]
Esto puede escribirse como producto fila por columna:
\[ \eta_i = \begin{pmatrix} 1 & x_i \end{pmatrix} \begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}. \]
Es decir:
\[ \eta_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}, \]
donde:
\[ \boldsymbol{x}_i^\top= \begin{pmatrix} 1 & x_i \end{pmatrix} \]
y:
\[ \boldsymbol{\beta} = \begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}. \]
Solución: matriz de diseño II
Apilando todas las observaciones:
\[ \boldsymbol{\eta} = \begin{pmatrix} \eta_1\\ \eta_2\\ \vdots\\ \eta_n \end{pmatrix}. \]
Como:
\[ \eta_i=\beta_1+\beta_2x_i, \]
se tiene:
\[ \boldsymbol{\eta} = \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{pmatrix} \begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}. \]
Por tanto:
\[ \boldsymbol{\eta} = \boldsymbol{X}\boldsymbol{\beta}, \]
donde:
\[ \boldsymbol{X} = \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{pmatrix}. \]
Solución: función de varianza
Para una variable Poisson:
\[ E(Y_i)=\lambda_i, \qquad Var(Y_i)=\lambda_i. \]
Como:
\[ \mu_i=\lambda_i, \]
entonces:
\[ Var(Y_i)=\mu_i. \]
En la notación del GLM:
\[ Var(Y_i)=\phi^{-1}V(\mu_i). \]
Para el modelo Poisson:
\[ \phi=1. \]
Por tanto:
\[ V(\mu_i)=\mu_i. \]
Solución: pesos I
La matriz de pesos se define como:
\[ \boldsymbol{W} = diag(w_1,\ldots,w_n), \]
donde:
\[ w_i = \frac{ \left( \dfrac{d\mu_i}{d\eta_i} \right)^2 }{ V(\mu_i) }. \]
Ya se obtuvo que:
\[ V(\mu_i)=\mu_i. \]
Ahora se necesita calcular:
\[ \frac{d\mu_i}{d\eta_i}. \]
Solución: derivada de la media
Como el enlace es log:
\[ \eta_i=\log(\mu_i). \]
Entonces:
\[ \mu_i=e^{\eta_i}. \]
Derivando respecto de \(\eta_i\):
\[ \frac{d\mu_i}{d\eta_i} = \frac{d}{d\eta_i}e^{\eta_i}. \]
Por tanto:
\[ \frac{d\mu_i}{d\eta_i} = e^{\eta_i}. \]
Pero:
\[ e^{\eta_i}=\mu_i. \]
Así:
\[ \frac{d\mu_i}{d\eta_i} = \mu_i. \]
Solución: pesos II
Sustituyendo en la definición de \(w_i\):
\[ w_i = \frac{ \left( \dfrac{d\mu_i}{d\eta_i} \right)^2 }{ V(\mu_i) }. \]
Como:
\[ \frac{d\mu_i}{d\eta_i}=\mu_i \]
y:
\[ V(\mu_i)=\mu_i, \]
entonces:
\[ w_i = \frac{\mu_i^2}{\mu_i}. \]
Por tanto:
\[ w_i=\mu_i. \]
Solución: matriz de pesos final
Como:
\[ w_i=\mu_i, \qquad i=1,\ldots,n, \]
la matriz de pesos es:
\[ \boldsymbol{W} = diag(\mu_1,\ldots,\mu_n). \]
Es decir:
\[ \boldsymbol{W} = \begin{pmatrix} \mu_1 & 0 & \cdots & 0\\ 0 & \mu_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \mu_n \end{pmatrix}. \]
Como:
\[ \mu_i=\exp(\beta_1+\beta_2x_i), \]
también puede escribirse como:
\[ \boldsymbol{W} = diag\left\{ \exp(\beta_1+\beta_2x_1), \ldots, \exp(\beta_1+\beta_2x_n) \right\}. \]
Solución: matriz de pesos ajustada
En la práctica, la matriz de pesos se evalúa en el modelo ajustado.
Es decir, se reemplaza:
\[ \beta_1,\beta_2 \]
por:
\[ \hat{\beta}_1,\hat{\beta}_2. \]
Entonces:
\[ \hat{\mu}_i = \exp(\hat{\beta}_1+\hat{\beta}_2x_i). \]
Por tanto:
\[ \hat{\boldsymbol{W}} = diag(\hat{\mu}_1,\ldots,\hat{\mu}_n). \]
Equivalentemente:
\[ \hat{\boldsymbol{W}} = diag\left\{ \exp(\hat{\beta}_1+\hat{\beta}_2x_1), \ldots, \exp(\hat{\beta}_1+\hat{\beta}_2x_n) \right\}. \]
Resumen
Para el modelo:
\[ Y_i\sim Poisson(\lambda_i), \qquad \log(\lambda_i)=\beta_1+\beta_2x_i, \]
se tiene:
\[ \mu_i=\lambda_i = \exp(\beta_1+\beta_2x_i). \]
La matriz de diseño es:
\[ \boldsymbol{X} = \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{pmatrix}. \]
La función de varianza es:
\[ V(\mu_i)=\mu_i. \]
La matriz de pesos es:
\[ \boldsymbol{W} = diag(\mu_1,\ldots,\mu_n). \]
Función score en forma matricial
Sea:
\[ \boldsymbol{y} = (y_1,\ldots,y_n)^\top, \qquad \boldsymbol{\mu} = (\mu_1,\ldots,\mu_n)^\top. \]
Sea \(\boldsymbol{X}\) la matriz de diseño de dimensión \(n\times p\).
Entonces, la función score para \(\boldsymbol{\beta}\) puede escribirse como:
\[ \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top \boldsymbol{W}^{1/2} \boldsymbol{V}^{-1/2} (\boldsymbol{y}-\boldsymbol{\mu}). \]
Equivalentemente, usando la derivada del enlace inverso:
\[ \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top \boldsymbol{\Delta} \boldsymbol{V}^{-1} (\boldsymbol{y}-\boldsymbol{\mu}), \]
donde:
\[ \boldsymbol{\Delta} = diag\left( \frac{d\mu_1}{d\eta_1}, \ldots, \frac{d\mu_n}{d\eta_n} \right). \]
Información de Fisher para \(\boldsymbol{\beta}\) I
Para obtener la información de Fisher se parte de la segunda derivada:
\[ \frac{\partial^2 L(\boldsymbol{\theta})} {\partial \beta_j\partial \beta_l}. \]
Después de derivar la función score, aparecen tres tipos de términos:
\[ \phi \sum_{i=1}^n (y_i-\mu_i) \frac{d^2\theta_i}{d\mu_i^2} \left( \frac{d\mu_i}{d\eta_i} \right)^2 x_{ij}x_{il}, \]
\[ \phi \sum_{i=1}^n (y_i-\mu_i) \frac{d\theta_i}{d\mu_i} \frac{d^2\mu_i}{d\eta_i^2} x_{ij}x_{il}, \]
y
\[ -\phi \sum_{i=1}^n \frac{d\theta_i}{d\mu_i} \left( \frac{d\mu_i}{d\eta_i} \right)^2 x_{ij}x_{il}. \]
Información de Fisher para \(\boldsymbol{\beta}\) II
Como:
\[ E(Y_i-\mu_i)=0, \]
los términos que contienen \((Y_i-\mu_i)\) desaparecen al tomar esperanza.
Por tanto:
\[ E\left[ \frac{\partial^2 L(\boldsymbol{\theta})} {\partial \beta_j\partial \beta_l} \right] = -\phi \sum_{i=1}^n \frac{d\theta_i}{d\mu_i} \left( \frac{d\mu_i}{d\eta_i} \right)^2 x_{ij}x_{il}. \]
Como:
\[ \frac{d\theta_i}{d\mu_i} = \frac{1}{V_i}, \]
se obtiene:
\[ E\left[ \frac{\partial^2 L(\boldsymbol{\theta})} {\partial \beta_j\partial \beta_l} \right] = -\phi \sum_{i=1}^n w_i x_{ij}x_{il}. \]
Información de Fisher para \(\boldsymbol{\beta}\) III
La información de Fisher para \(\boldsymbol{\beta}\) se define como:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}(\boldsymbol{\theta}) = E\left[ - \frac{\partial^2 L(\boldsymbol{\theta})} {\partial \boldsymbol{\beta}\partial\boldsymbol{\beta}^\top} \right]. \]
Usando la matriz de pesos \(\boldsymbol{W}\):
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top \boldsymbol{W} \boldsymbol{X}. \]
Esta matriz resume la información disponible para estimar \(\boldsymbol{\beta}\).
Depende de:
- la matriz de diseño \(\boldsymbol{X}\);
- la función de varianza;
- la función de enlace;
- el parámetro de precisión \(\phi\).
Casos particulares de \(U_{\boldsymbol{\beta}}\) y \(K_{\boldsymbol{\beta}\boldsymbol{\beta}}\)
A partir de la forma general:
\[ \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top \boldsymbol{W}^{1/2} \boldsymbol{V}^{-1/2} (\boldsymbol{y}-\boldsymbol{\mu}), \]
y
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top \boldsymbol{W} \boldsymbol{X}, \]
se obtienen las siguientes expresiones para algunos modelos usuales:
| Modelo | \(\boldsymbol{U}_{\boldsymbol{\beta}}\) | \(\boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}\) |
|---|---|---|
| Normal | \(\sigma^{-2}\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\sigma^{-2}\boldsymbol{X}^\top\boldsymbol{X}\) |
| Poisson | \(\boldsymbol{X}^\top\boldsymbol{V}^{-1/2}(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\boldsymbol{X}^\top\boldsymbol{X}\) |
| Binomial | \(\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\boldsymbol{X}^\top\boldsymbol{V}\boldsymbol{X}\) |
| Gamma | \(\phi\,\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\phi\,\boldsymbol{X}^\top\boldsymbol{X}\) |
| Normal inversa | \(\phi\,\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\phi\,\boldsymbol{X}^\top\boldsymbol{V}\boldsymbol{X}\) |
donde:
\[ \boldsymbol{V} = diag\{V(\mu_1),\ldots,V(\mu_n)\}. \]
Enlace canónico: simplificación
Cuando el enlace es canónico:
\[ \theta_i=\eta_i. \]
Entonces:
\[ \frac{d\mu_i}{d\eta_i} = \frac{d\mu_i}{d\theta_i} = V_i. \]
Por tanto:
\[ w_i = \frac{V_i^2}{V_i} = V_i. \]
Así:
\[ \boldsymbol{W}=\boldsymbol{V}. \]
Score e información con enlace canónico
Bajo enlace canónico, la función score se simplifica a:
\[ \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top (\boldsymbol{y}-\boldsymbol{\mu}). \]
Además, la información de Fisher queda:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top \boldsymbol{V} \boldsymbol{X}. \]
Esta simplificación es una de las razones por las que los enlaces canónicos son importantes en la inferencia de los GLM.
Casos con enlace canónico: resumen
| Modelo | Enlace canónico | \(\phi\) | \(V(\mu)\) | \(\boldsymbol{U}_{\boldsymbol{\beta}}\) | \(\boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}\) |
|---|---|---|---|---|---|
| Normal | \(\mu\) | \(\sigma^{-2}\) | \(1\) | \(\sigma^{-2}\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\sigma^{-2}\boldsymbol{X}^\top\boldsymbol{X}\) |
| Poisson | \(\log(\mu)\) | \(1\) | \(\mu\) | \(\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\boldsymbol{X}^\top\boldsymbol{V}\boldsymbol{X}\) |
| Binomial | \(\log\left(\dfrac{\mu}{1-\mu}\right)\) | \(m\) | \(\mu(1-\mu)\) | \(m\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(m\boldsymbol{X}^\top\boldsymbol{V}\boldsymbol{X}\) |
| Gamma | \(\mu^{-1}\) | \(\phi\) | \(\mu^2\) | \(\phi\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\phi\boldsymbol{X}^\top\boldsymbol{V}\boldsymbol{X}\) |
| Normal inversa | \(\mu^{-2}\) | \(\phi\) | \(\mu^3\) | \(\phi\boldsymbol{X}^\top(\boldsymbol{y}-\boldsymbol{\mu})\) | \(\phi\boldsymbol{X}^\top\boldsymbol{V}\boldsymbol{X}\) |
Esta tabla muestra la simplificación que se obtiene cuando \(\theta_i=\eta_i\).
Partición de parámetros
Si se particiona:
\[ \boldsymbol{\beta} = (\boldsymbol{\beta}_1^\top,\boldsymbol{\beta}_2^\top)^\top, \]
también se puede particionar la matriz de diseño:
\[ \boldsymbol{X} = [\boldsymbol{X}_1 \ \boldsymbol{X}_2]. \]
La función score para \(\boldsymbol{\beta}_1\) es:
\[ \boldsymbol{U}_{\boldsymbol{\beta}_1}(\boldsymbol{\theta}) = \phi \boldsymbol{X}_1^\top \boldsymbol{W}^{1/2} \boldsymbol{V}^{-1/2} (\boldsymbol{y}-\boldsymbol{\mu}). \]
La información de Fisher correspondiente es:
\[ \boldsymbol{K}_{\boldsymbol{\beta}_1\boldsymbol{\beta}_1}(\boldsymbol{\theta}) = \phi \boldsymbol{X}_1^\top \boldsymbol{W} \boldsymbol{X}_1. \]
Función score para \(\phi\)
El parámetro \(\phi\) es el parámetro de precisión.
La función score para \(\phi\) se obtiene derivando la log-verosimilitud respecto de \(\phi\):
\[ U_\phi(\boldsymbol{\theta}) = \frac{\partial L(\boldsymbol{\theta})}{\partial \phi}. \]
Como:
\[ L(\boldsymbol{\theta}) = \sum_{i=1}^n \left[ \phi\{y_i\theta_i-b(\theta_i)\} + c(y_i,\phi) \right], \]
se tiene:
\[ U_\phi(\boldsymbol{\theta}) = \sum_{i=1}^n \{y_i\theta_i-b(\theta_i)\} + \sum_{i=1}^n c'(y_i,\phi), \]
donde:
\[ c'(y_i,\phi)=\frac{\partial c(y_i,\phi)}{\partial \phi}. \]
Información de Fisher para \(\phi\)
La información de Fisher para \(\phi\) está dada por:
\[ K_{\phi\phi}(\boldsymbol{\theta}) = E\left[ - \frac{\partial^2 L(\boldsymbol{\theta})}{\partial \phi^2} \right]. \]
Como la parte que depende de \(\phi\) en la segunda derivada proviene de \(c(y_i,\phi)\):
\[ K_{\phi\phi}(\boldsymbol{\theta}) = - \sum_{i=1}^n E\left[ c''(Y_i,\phi) \right], \]
donde:
\[ c''(Y_i,\phi) = \frac{\partial^2 c(Y_i,\phi)}{\partial \phi^2}. \]
Ortogonalidad entre \(\boldsymbol{\beta}\) y \(\phi\)
La derivada cruzada entre \(\boldsymbol{\beta}\) y \(\phi\) es proporcional a:
\[ \sum_{i=1}^n \frac{1}{V_i} \left( \frac{d\mu_i}{d\eta_i} \right) (y_i-\mu_i)\boldsymbol{x}_i. \]
Como:
\[ E(Y_i-\mu_i)=0, \]
entonces:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\phi}(\boldsymbol{\theta}) = \boldsymbol{0}. \]
Por tanto, \(\boldsymbol{\beta}\) y \(\phi\) son ortogonales en el sentido de la información de Fisher.
Matriz de información completa
Debido a la ortogonalidad entre \(\boldsymbol{\beta}\) y \(\phi\), la matriz de información de Fisher tiene forma bloque diagonal:
\[ \boldsymbol{K}_{\boldsymbol{\theta}\boldsymbol{\theta}} = \begin{pmatrix} \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}} & \boldsymbol{0}\\ \boldsymbol{0} & K_{\phi\phi} \end{pmatrix}. \]
La función score completa puede escribirse como:
\[ \boldsymbol{U}_{\boldsymbol{\theta}} = \begin{pmatrix} \boldsymbol{U}_{\boldsymbol{\beta}}\\ U_\phi \end{pmatrix}. \]
Esta separación facilita el estudio de la inferencia para \(\boldsymbol{\beta}\) y para \(\phi\).
Estimación de \(\boldsymbol{\beta}\) I
La estimación de \(\boldsymbol{\beta}\) se realiza resolviendo:
\[ \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\theta})=\boldsymbol{0}. \]
En general, estas ecuaciones no tienen solución cerrada.
Por eso se usa un procedimiento iterativo.
En los GLM, el método usual es el de mínimos cuadrados reponderados iterativos.
También se conoce como:
\[ IRWLS. \]
Estimación de \(\boldsymbol{\beta}\) II
En la iteración \(m\), se actualiza la estimación de \(\boldsymbol{\beta}\) mediante:
\[ \boldsymbol{\beta}^{(m+1)} = \left[ \boldsymbol{X}^\top \boldsymbol{W}^{(m)} \boldsymbol{X} \right]^{-1} \boldsymbol{X}^\top \boldsymbol{W}^{(m)} \boldsymbol{z}^{(m)}. \]
Aquí:
\[ \boldsymbol{W}^{(m)} \]
es la matriz de pesos evaluada con las medias de la iteración \(m\).
Por tanto, los pesos cambian en cada paso del algoritmo.
Variable dependiente modificada
La variable dependiente modificada se define como:
\[ \boldsymbol{z}^{(m)} = \boldsymbol{\eta}^{(m)} + \left[ \boldsymbol{W}^{(m)} \right]^{-1/2} \left[ \boldsymbol{V}^{(m)} \right]^{-1/2} (\boldsymbol{y}-\boldsymbol{\mu}^{(m)}). \]
Esta variable cumple el rol de una respuesta auxiliar en cada iteración.
En cada paso, el problema se aproxima por un modelo lineal ponderado.
Por eso el procedimiento puede interpretarse como una sucesión de ajustes de mínimos cuadrados ponderados.
Lectura del algoritmo IRWLS
El algoritmo IRWLS puede resumirse así:
- partir de un valor inicial \(\boldsymbol{\beta}^{(0)}\);
- calcular \(\boldsymbol{\eta}^{(m)}\) y \(\boldsymbol{\mu}^{(m)}\);
- construir \(\boldsymbol{W}^{(m)}\) y \(\boldsymbol{z}^{(m)}\);
- actualizar \(\boldsymbol{\beta}^{(m+1)}\);
- repetir hasta convergencia.
La convergencia se alcanza cuando los cambios en \(\boldsymbol{\beta}\), en la log-verosimilitud o en el desvío son suficientemente pequeños.
Estimación de \(\phi\) I
El parámetro \(\phi\) es el parámetro de precisión.
Para algunos modelos, su estimación depende del desvío.
En el caso normal y en el caso normal inversa, el estimador de máxima verosimilitud puede expresarse como:
\[ \hat{\phi} = \frac{n}{D(\boldsymbol{y};\hat{\boldsymbol{\mu}})}. \]
Por tanto, si el desvío es grande, la precisión estimada disminuye.
Si el desvío es pequeño, la precisión estimada aumenta.
Estimación de \(\phi\) II
En el caso Gamma, el estimador de máxima verosimilitud de \(\phi\) se obtiene resolviendo:
\[ 2n \left[ \log(\hat{\phi}) - \psi(\hat{\phi}) \right] = D(\boldsymbol{y};\hat{\boldsymbol{\mu}}), \]
donde \(\psi(\cdot)\) es la función digamma.
Esta ecuación no tiene, en general, una solución cerrada simple.
Por ello, se resuelve numéricamente.
Estimación de \(\phi\) III
Un estimador alternativo de \(\phi\) se basa en la estadística de Pearson.
Como:
\[ Var(Y_i)=\phi^{-1}V(\mu_i), \]
se considera:
\[ \sum_{i=1}^n \frac{(y_i-\hat{\mu}_i)^2} {V(\hat{\mu}_i)}. \]
Un estimador de momentos para \(\phi\) es:
\[ \hat{\phi} = \frac{ n-p }{ \sum_{i=1}^n \dfrac{(y_i-\hat{\mu}_i)^2}{V(\hat{\mu}_i)} }. \]
Este estimador usa la variabilidad residual estandarizada por la función de varianza.
Estimación de \(\phi\) IV
También aparece una estimación basada en el desvío:
\[ \hat{\phi} = \frac{n-p} {D(\boldsymbol{y};\hat{\boldsymbol{\mu}})}. \]
Equivalente:
\[ \hat{\phi}^{-1} = \frac{ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) }{ n-p }. \]
Esta forma usa el desvío residual medio como estimación del factor de dispersión.
Sin embargo, debe distinguirse del estimador de máxima verosimilitud, que en algunos casos usa \(n\) en lugar de \(n-p\).
Síntesis de estimación
La estimación en GLM se apoya en tres elementos:
\[ \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top \boldsymbol{W}^{1/2} \boldsymbol{V}^{-1/2} (\boldsymbol{y}-\boldsymbol{\mu}), \]
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^\top \boldsymbol{W} \boldsymbol{X}, \]
y el algoritmo IRWLS:
\[ \boldsymbol{\beta}^{(m+1)} = \left[ \boldsymbol{X}^\top \boldsymbol{W}^{(m)} \boldsymbol{X} \right]^{-1} \boldsymbol{X}^\top \boldsymbol{W}^{(m)} \boldsymbol{z}^{(m)}. \]
Con esto se obtiene \(\hat{\boldsymbol{\beta}}\) y, posteriormente, se estima \(\phi\) cuando corresponde.
Distribución asintótica de \(\hat{\boldsymbol{\beta}}\)
Bajo condiciones regulares, el estimador de máxima verosimilitud de \(\boldsymbol{\beta}\) es consistente y asintóticamente normal.
Como:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}(\boldsymbol{\theta}) = \phi \boldsymbol{X}^{\top} \boldsymbol{W} \boldsymbol{X}, \]
se tiene:
\[ \hat{\boldsymbol{\beta}} \approx N_p \left( \boldsymbol{\beta}, \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}^{-1}(\boldsymbol{\theta}) \right). \]
En la práctica, se reemplaza \(\boldsymbol{\theta}\) por su estimador:
\[ \hat{\boldsymbol{\theta}} = (\hat{\boldsymbol{\beta}}^\top,\hat{\phi})^\top. \]
Por tanto:
\[ \hat{\boldsymbol{\beta}} \approx N_p \left( \boldsymbol{\beta}, \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}^{-1}(\hat{\boldsymbol{\theta}}) \right). \]
Matriz de varianzas y covarianzas
La matriz de varianzas y covarianzas aproximada de \(\hat{\boldsymbol{\beta}}\) es:
\[ \widehat{Var}(\hat{\boldsymbol{\beta}}) = \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}^{-1}(\hat{\boldsymbol{\theta}}). \]
Como:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}(\hat{\boldsymbol{\theta}}) = \hat{\phi} \boldsymbol{X}^{\top} \hat{\boldsymbol{W}} \boldsymbol{X}, \]
entonces:
\[ \widehat{Var}(\hat{\boldsymbol{\beta}}) = \left[ \hat{\phi} \boldsymbol{X}^{\top} \hat{\boldsymbol{W}} \boldsymbol{X} \right]^{-1}. \]
Los errores estándar se obtienen como:
\[ EP(\hat{\beta}_j) = \sqrt{ \left[ \widehat{Var}(\hat{\boldsymbol{\beta}}) \right]_{jj} }. \]
Distribución asintótica componente a componente
Para cada parámetro \(\beta_j\):
\[ \hat{\beta}_j \approx N \left( \beta_j, EP^2(\hat{\beta}_j) \right). \]
Por tanto:
\[ Z_j = \frac{ \hat{\beta}_j-\beta_j }{ EP(\hat{\beta}_j) } \approx N(0,1). \]
Esta aproximación permite construir intervalos de confianza y contrastes de hipótesis para coeficientes individuales.
Un intervalo aproximado de confianza para \(\beta_j\) es:
\[ \hat{\beta}_j \pm z_{1-\alpha/2} EP(\hat{\beta}_j). \]
Ejercicio guiado: score, Fisher y varianza asintótica
Suponga que:
\[ Y_i\sim Poisson(\lambda_i), \qquad i=1,\ldots,n, \]
y considere el modelo:
\[ \log(\lambda_i)=\beta_1+\beta_2x_i. \]
Sea:
\[ \boldsymbol{\beta} = \begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}. \]
Obtenga:
La función score \(\boldsymbol{U}_{\boldsymbol{\beta}}\).
La matriz de información de Fisher \(\boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}\).
La matriz de varianzas y covarianzas asintótica de \(\hat{\boldsymbol{\beta}}\).
Solución: resumen de elementos ya obtenidos
Para el modelo Poisson:
\[ \mu_i=E(Y_i)=\lambda_i. \]
Como:
\[ \log(\lambda_i)=\beta_1+\beta_2x_i, \]
entonces:
\[ \log(\mu_i)=\beta_1+\beta_2x_i. \]
Por tanto:
\[ \eta_i=\beta_1+\beta_2x_i, \]
y:
\[ \mu_i=\exp(\beta_1+\beta_2x_i). \]
La matriz de diseño es:
\[ \boldsymbol{X} = \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{pmatrix}. \]
Solución: resumen de pesos
Para Poisson:
\[ V(\mu_i)=\mu_i, \qquad \phi=1. \]
Con enlace log:
\[ \eta_i=\log(\mu_i), \qquad \mu_i=e^{\eta_i}. \]
Entonces:
\[ \frac{d\mu_i}{d\eta_i}=\mu_i. \]
Luego:
\[ w_i = \frac{ \left( \dfrac{d\mu_i}{d\eta_i} \right)^2 }{ V(\mu_i) } = \frac{\mu_i^2}{\mu_i} = \mu_i. \]
Por tanto:
\[ \boldsymbol{W}=diag(\mu_1,\ldots,\mu_n). \]
Además:
\[ \boldsymbol{V}=diag(V(\mu_1),\ldots,V(\mu_n)) = diag(\mu_1,\ldots,\mu_n). \]
Solución: función score general
En la notación de Paula, la función score para \(\boldsymbol{\beta}\) es:
\[ \boldsymbol{U}_{\boldsymbol{\beta}} = \phi \boldsymbol{X}^{\top} \boldsymbol{W}^{1/2} \boldsymbol{V}^{-1/2} (\boldsymbol{y}-\boldsymbol{\mu}). \]
Para este modelo:
\[ \phi=1, \]
\[ \boldsymbol{W}=diag(\mu_1,\ldots,\mu_n), \]
y:
\[ \boldsymbol{V}=diag(\mu_1,\ldots,\mu_n). \]
Por tanto:
\[ \boldsymbol{W}^{1/2} \boldsymbol{V}^{-1/2} = \boldsymbol{I}_n. \]
Solución: función score simplificada
Sustituyendo en la expresión general:
\[ \boldsymbol{U}_{\boldsymbol{\beta}} = \phi \boldsymbol{X}^{\top} \boldsymbol{W}^{1/2} \boldsymbol{V}^{-1/2} (\boldsymbol{y}-\boldsymbol{\mu}), \]
se obtiene:
\[ \boldsymbol{U}_{\boldsymbol{\beta}} = 1\cdot \boldsymbol{X}^{\top} \boldsymbol{I}_n (\boldsymbol{y}-\boldsymbol{\mu}). \]
Entonces:
\[ \boldsymbol{U}_{\boldsymbol{\beta}} = \boldsymbol{X}^{\top} (\boldsymbol{y}-\boldsymbol{\mu}). \]
Solución: score explícito I
Como:
\[ \boldsymbol{X}^{\top} = \begin{pmatrix} 1 & 1 & \cdots & 1\\ x_1 & x_2 & \cdots & x_n \end{pmatrix}, \]
y:
\[ \boldsymbol{y}-\boldsymbol{\mu} = \begin{pmatrix} y_1-\mu_1\\ y_2-\mu_2\\ \vdots\\ y_n-\mu_n \end{pmatrix}, \]
entonces:
\[ \boldsymbol{U}_{\boldsymbol{\beta}} = \begin{pmatrix} 1 & 1 & \cdots & 1\\ x_1 & x_2 & \cdots & x_n \end{pmatrix} \begin{pmatrix} y_1-\mu_1\\ y_2-\mu_2\\ \vdots\\ y_n-\mu_n \end{pmatrix}. \]
Solución: score explícito II
Realizando el producto matricial:
\[ \boldsymbol{U}_{\boldsymbol{\beta}} = \begin{pmatrix} \sum_{i=1}^n 1(y_i-\mu_i)\\[6pt] \sum_{i=1}^n x_i(y_i-\mu_i) \end{pmatrix}. \]
Por tanto:
\[ \boldsymbol{U}_{\boldsymbol{\beta}} = \begin{pmatrix} \sum_{i=1}^n (y_i-\mu_i)\\[6pt] \sum_{i=1}^n x_i(y_i-\mu_i) \end{pmatrix}. \]
Solución: información de Fisher general
En la notación de Paula, la matriz de información de Fisher para \(\boldsymbol{\beta}\) es:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}} = \phi \boldsymbol{X}^{\top} \boldsymbol{W} \boldsymbol{X}. \]
Para Poisson:
\[ \phi=1. \]
Entonces:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}} = \boldsymbol{X}^{\top} \boldsymbol{W} \boldsymbol{X}. \]
Como:
\[ \boldsymbol{W}=diag(\mu_1,\ldots,\mu_n), \]
se debe calcular:
\[ \boldsymbol{X}^{\top} diag(\mu_1,\ldots,\mu_n) \boldsymbol{X}. \]
Solución: cálculo de \(\boldsymbol{W}\boldsymbol{X}\)
Se tiene:
\[ \boldsymbol{W} = \begin{pmatrix} \mu_1 & 0 & \cdots & 0\\ 0 & \mu_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \mu_n \end{pmatrix}, \]
y:
\[ \boldsymbol{X} = \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{pmatrix}. \]
Entonces:
\[ \boldsymbol{W}\boldsymbol{X} = \begin{pmatrix} \mu_1 & \mu_1x_1\\ \mu_2 & \mu_2x_2\\ \vdots & \vdots\\ \mu_n & \mu_nx_n \end{pmatrix}. \]
Solución: cálculo de \(\boldsymbol{X}^{\top}\boldsymbol{W}\boldsymbol{X}\) I
Ahora:
\[ \boldsymbol{X}^{\top} = \begin{pmatrix} 1 & 1 & \cdots & 1\\ x_1 & x_2 & \cdots & x_n \end{pmatrix}. \]
Por tanto:
\[ \boldsymbol{X}^{\top}\boldsymbol{W}\boldsymbol{X} = \begin{pmatrix} 1 & 1 & \cdots & 1\\ x_1 & x_2 & \cdots & x_n \end{pmatrix} \begin{pmatrix} \mu_1 & \mu_1x_1\\ \mu_2 & \mu_2x_2\\ \vdots & \vdots\\ \mu_n & \mu_nx_n \end{pmatrix}. \]
Solución: cálculo de \(\boldsymbol{X}^{\top}\boldsymbol{W}\boldsymbol{X}\) II
El elemento \((1,1)\) es:
\[ K_{11} = \sum_{i=1}^n 1\cdot \mu_i = \sum_{i=1}^n \mu_i. \]
El elemento \((1,2)\) es:
\[ K_{12} = \sum_{i=1}^n 1\cdot \mu_i x_i = \sum_{i=1}^n \mu_i x_i. \]
El elemento \((2,1)\) es:
\[ K_{21} = \sum_{i=1}^n x_i\mu_i = \sum_{i=1}^n \mu_i x_i. \]
Solución: cálculo de \(\boldsymbol{X}^{\top}\boldsymbol{W}\boldsymbol{X}\) III
El elemento \((2,2)\) es:
\[ K_{22} = \sum_{i=1}^n x_i(\mu_i x_i). \]
Por tanto:
\[ K_{22} = \sum_{i=1}^n \mu_i x_i^2. \]
Así:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}} = \begin{pmatrix} \sum_{i=1}^n \mu_i & \sum_{i=1}^n \mu_i x_i\\[6pt] \sum_{i=1}^n \mu_i x_i & \sum_{i=1}^n \mu_i x_i^2 \end{pmatrix}. \]
Solución: notación resumida
Defina:
\[ S_0=\sum_{i=1}^n\mu_i, \]
\[ S_1=\sum_{i=1}^n\mu_i x_i, \]
y:
\[ S_2=\sum_{i=1}^n\mu_i x_i^2. \]
Entonces:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}} = \begin{pmatrix} S_0 & S_1\\ S_1 & S_2 \end{pmatrix}. \]
Solución: varianza asintótica I
Bajo condiciones regulares:
\[ \hat{\boldsymbol{\beta}} \approx N_2 \left( \boldsymbol{\beta}, \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}^{-1} \right). \]
Por tanto:
\[ Var(\hat{\boldsymbol{\beta}}) \approx \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}^{-1}. \]
Como:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}} = \begin{pmatrix} S_0 & S_1\\ S_1 & S_2 \end{pmatrix}, \]
se debe invertir esta matriz.
Solución: varianza asintótica II
Para una matriz:
\[ \boldsymbol{A} = \begin{pmatrix} a & b\\ b & c \end{pmatrix}, \]
su inversa es:
\[ \boldsymbol{A}^{-1} = \frac{1}{ac-b^2} \begin{pmatrix} c & -b\\ -b & a \end{pmatrix}. \]
En este caso:
\[ a=S_0, \qquad b=S_1, \qquad c=S_2. \]
Entonces:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}^{-1} = \frac{1}{S_0S_2-S_1^2} \begin{pmatrix} S_2 & -S_1\\ -S_1 & S_0 \end{pmatrix}. \]
Solución: matriz de varianza asintótica
Por tanto:
\[ Var(\hat{\boldsymbol{\beta}}) \approx \frac{1}{S_0S_2-S_1^2} \begin{pmatrix} S_2 & -S_1\\ -S_1 & S_0 \end{pmatrix}. \]
Es decir:
\[ Var(\hat{\beta}_1) \approx \frac{S_2}{S_0S_2-S_1^2}, \]
\[ Var(\hat{\beta}_2) \approx \frac{S_0}{S_0S_2-S_1^2}, \]
y:
\[ Cov(\hat{\beta}_1,\hat{\beta}_2) \approx -\frac{S_1}{S_0S_2-S_1^2}. \]
Solución: evaluación en el modelo ajustado
En la práctica, la matriz se evalúa reemplazando \(\mu_i\) por \(\hat{\mu}_i\).
Como:
\[ \hat{\mu}_i = \exp(\hat{\beta}_1+\hat{\beta}_2x_i), \]
defina:
\[ \hat{S}_0=\sum_{i=1}^n\hat{\mu}_i, \qquad \hat{S}_1=\sum_{i=1}^n\hat{\mu}_i x_i, \qquad \hat{S}_2=\sum_{i=1}^n\hat{\mu}_i x_i^2. \]
Entonces:
\[ \widehat{Var}(\hat{\boldsymbol{\beta}}) = \frac{1}{\hat{S}_0\hat{S}_2-\hat{S}_1^2} \begin{pmatrix} \hat{S}_2 & -\hat{S}_1\\ -\hat{S}_1 & \hat{S}_0 \end{pmatrix}. \]
Resumen del ejercicio
Para el modelo:
\[ Y_i\sim Poisson(\lambda_i), \qquad \log(\lambda_i)=\beta_1+\beta_2x_i, \]
se obtuvo:
\[ \boldsymbol{U}_{\boldsymbol{\beta}} = \begin{pmatrix} \sum_{i=1}^n (y_i-\mu_i)\\[6pt] \sum_{i=1}^n x_i(y_i-\mu_i) \end{pmatrix}. \]
Además:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}} = \begin{pmatrix} S_0 & S_1\\ S_1 & S_2 \end{pmatrix}, \]
donde:
\[ S_0=\sum_{i=1}^n\mu_i, \qquad S_1=\sum_{i=1}^n\mu_i x_i, \qquad S_2=\sum_{i=1}^n\mu_i x_i^2. \]
Finalmente:
\[ Var(\hat{\boldsymbol{\beta}}) \approx \frac{1}{S_0S_2-S_1^2} \begin{pmatrix} S_2 & -S_1\\ -S_1 & S_0 \end{pmatrix}. \]
Contrastes de hipótesis
Los contrastes de hipótesis permiten evaluar restricciones sobre los parámetros del GLM.
Una hipótesis simple sobre todo el vector de parámetros puede escribirse como:
\[ H_0: \boldsymbol{\beta} = \boldsymbol{\beta}^{0} \]
frente a:
\[ H_1: \boldsymbol{\beta} \neq \boldsymbol{\beta}^{0}. \]
También puede interesar contrastar solo un subconjunto de parámetros.
Para ello se usa una partición:
\[ \boldsymbol{\beta} = (\boldsymbol{\beta}_1^\top,\boldsymbol{\beta}_2^\top)^\top, \]
donde \(\boldsymbol{\beta}_1\) tiene dimensión \(q\).
Hipótesis sobre un subconjunto de parámetros
La hipótesis sobre un subconjunto de coeficientes puede escribirse como:
\[ H_0: \boldsymbol{\beta}_1 = \boldsymbol{\beta}_1^0 \]
frente a:
\[ H_1: \boldsymbol{\beta}_1 \neq \boldsymbol{\beta}_1^0. \]
Un caso frecuente es:
\[ H_0: \boldsymbol{\beta}_1 = \boldsymbol{0}. \]
Esta hipótesis evalúa si un bloque de variables puede retirarse del modelo.
Los tres contrastes clásicos son:
\[ \text{razón de verosimilitud}, \qquad \text{Wald}, \qquad \text{score}. \]
Contraste de razón de verosimilitud I
Para la hipótesis simple:
\[ H_0: \boldsymbol{\beta} = \boldsymbol{\beta}^{0}, \]
el contraste de razón de verosimilitud compara la log-verosimilitud bajo el modelo restringido con la log-verosimilitud maximizada.
La estadística es:
\[ \xi_{RV} = 2 \left[ L(\hat{\boldsymbol{\beta}}) - L(\boldsymbol{\beta}^{0}) \right]. \]
Bajo \(H_0\):
\[ \xi_{RV} \sim \chi^2_p. \]
Valores grandes de \(\xi_{RV}\) indican evidencia contra \(H_0\).
Contraste de razón de verosimilitud II
En GLM, esta estadística puede escribirse usando el desvío.
Sea:
\[ \hat{\boldsymbol{\mu}}^0 = g^{-1}(\boldsymbol{X}\boldsymbol{\beta}^{0}) \]
el vector de medias bajo \(H_0\).
Sea:
\[ \hat{\boldsymbol{\mu}} = g^{-1}(\boldsymbol{X}\hat{\boldsymbol{\beta}}) \]
el vector de medias bajo el modelo ajustado.
Entonces:
\[ \xi_{RV} = \phi \left[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \right]. \]
Contraste de razón de verosimilitud III
Para una hipótesis sobre \(q\) parámetros:
\[ H_0: \boldsymbol{\beta}_1 = \boldsymbol{\beta}_1^0, \]
la estadística de razón de verosimilitud es:
\[ \xi_{RV} = \phi \left[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}) \right], \]
donde \(\hat{\boldsymbol{\mu}}^0\) se obtiene bajo la restricción impuesta por \(H_0\).
Bajo \(H_0\):
\[ \xi_{RV} \sim \chi^2_q. \]
Si \(\phi=1\), como en Poisson o binomial, la comparación se basa directamente en:
\[ D(\boldsymbol{y};\hat{\boldsymbol{\mu}}^0) - D(\boldsymbol{y};\hat{\boldsymbol{\mu}}). \]
Contraste de Wald I
El contraste de Wald usa la distancia entre el estimador y el valor propuesto por la hipótesis nula.
Para:
\[ H_0: \boldsymbol{\beta} = \boldsymbol{\beta}^{0}, \]
se define:
\[ \xi_W = (\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})^\top \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}(\hat{\boldsymbol{\theta}}) (\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0}). \]
Bajo \(H_0\):
\[ \xi_W \sim \chi^2_p. \]
Este contraste se basa en la distribución asintótica de \(\hat{\boldsymbol{\beta}}\).
Contraste de Wald II
Para contrastar un solo coeficiente:
\[ H_0: \beta_j=\beta_j^0, \]
la estadística de Wald es:
\[ Z_W = \frac{ \hat{\beta}_j-\beta_j^0 }{ EP(\hat{\beta}_j) }. \]
Bajo \(H_0\):
\[ Z_W \approx N(0,1). \]
Equivalentemente:
\[ Z_W^2 \approx \chi^2_1. \]
Este es el contraste que usualmente aparece en la tabla de coeficientes del modelo.
Contraste de Wald III
Para un subconjunto de parámetros:
\[ H_0: \boldsymbol{\beta}_1 = \boldsymbol{\beta}_1^0, \]
sea:
\[ \widehat{Var}(\hat{\boldsymbol{\beta}}_1) \]
la submatriz correspondiente dentro de:
\[ \widehat{Var}(\hat{\boldsymbol{\beta}}). \]
Entonces:
\[ \xi_W = (\hat{\boldsymbol{\beta}}_1-\boldsymbol{\beta}_1^0)^\top \left[ \widehat{Var}(\hat{\boldsymbol{\beta}}_1) \right]^{-1} (\hat{\boldsymbol{\beta}}_1-\boldsymbol{\beta}_1^0). \]
Bajo \(H_0\):
\[ \xi_W \sim \chi^2_q. \]
Contraste score I
El contraste score evalúa la pendiente de la log-verosimilitud bajo el modelo restringido.
Para:
\[ H_0: \boldsymbol{\beta} = \boldsymbol{\beta}^{0}, \]
se calcula la función score en el valor impuesto por la hipótesis nula:
\[ \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\beta}^{0}). \]
La estadística score es:
\[ \xi_S = \boldsymbol{U}_{\boldsymbol{\beta}}^\top(\boldsymbol{\beta}^{0}) \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}^{-1}(\boldsymbol{\beta}^{0}) \boldsymbol{U}_{\boldsymbol{\beta}}(\boldsymbol{\beta}^{0}). \]
Bajo \(H_0\):
\[ \xi_S \sim \chi^2_p. \]
Contraste score II
Para una hipótesis sobre un subconjunto de parámetros:
\[ H_0: \boldsymbol{\beta}_1 = \boldsymbol{\beta}_1^0, \]
se ajusta primero el modelo restringido.
Luego se evalúa si la función score asociada a los parámetros restringidos es suficientemente grande.
La idea es:
- si \(H_0\) es adecuada, el score bajo el modelo restringido no debe ser grande;
- si el score es grande, hay evidencia de que liberar esos parámetros mejoraría el modelo.
El contraste score es útil porque no requiere ajustar el modelo completo.
Comparación de los tres contrastes
Los tres contrastes se basan en la misma log-verosimilitud, pero usan información diferente.
| Contraste | Usa principalmente | Se evalúa en |
|---|---|---|
| Razón de verosimilitud | Diferencia de log-verosimilitudes o desvíos | Modelo restringido y modelo completo |
| Wald | Distancia entre \(\hat{\boldsymbol{\beta}}\) y el valor nulo | Modelo completo |
| Score | Pendiente de la log-verosimilitud bajo \(H_0\) | Modelo restringido |
Asintóticamente, bajo \(H_0\), los tres contrastes siguen una distribución chi-cuadrado con grados de libertad iguales al número de restricciones.
Decisión estadística
Para cualquiera de los tres contrastes, se calcula un valor-p usando:
\[ \chi^2_q, \]
donde \(q\) es el número de restricciones.
La regla de decisión es:
\[ \text{rechazar } H_0 \quad \text{si} \quad p\text{-valor}<\alpha. \]
Si se rechaza \(H_0\), se concluye que los parámetros evaluados aportan evidencia estadística dentro del modelo.
Si no se rechaza \(H_0\), no hay evidencia suficiente para descartar la restricción propuesta.
Lectura en términos del GLM
En un GLM, contrastar parámetros equivale a evaluar si ciertos predictores aportan a la explicación de la media:
\[ g(\mu_i) = \eta_i = \boldsymbol{x}_i^\top\boldsymbol{\beta}. \]
Por ejemplo, si:
\[ H_0: \beta_j=0, \]
entonces se evalúa si el predictor \(x_j\) contribuye al predictor lineal, manteniendo los demás términos del modelo.
Si se rechaza \(H_0\), el predictor se considera relevante dentro de la estructura del GLM ajustado.
Síntesis de contrastes
La inferencia en GLM usa la aproximación asintótica de:
\[ \hat{\boldsymbol{\beta}} \]
y la información de Fisher:
\[ \boldsymbol{K}_{\boldsymbol{\beta}\boldsymbol{\beta}}. \]
A partir de ellas se construyen tres contrastes principales:
\[ \xi_{RV}, \qquad \xi_W, \qquad \xi_S. \]
Bajo \(H_0\), todos tienen distribución aproximada:
\[ \chi^2_q. \]
Estos contrastes permiten evaluar hipótesis sobre coeficientes individuales, bloques de variables o restricciones generales del modelo.
Ejercicio guiado: contraste de Wald
Suponga que:
\[ Y_i\sim Poisson(\lambda_i), \qquad i=1,\ldots,n, \]
y considere el modelo:
\[ \log(\lambda_i)=\beta_1+\beta_2x_i. \]
Se desea contrastar:
\[ H_0:\beta_2=0 \]
frente a:
\[ H_1:\beta_2\neq 0. \]
A partir de la matriz estimada de varianzas y covarianzas de \(\hat{\boldsymbol{\beta}}\), construya el contraste de Wald para \(\beta_2\).
Solución: estimador y varianza requerida
Sea:
\[ \hat{\boldsymbol{\beta}} = \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2 \end{pmatrix}. \]
Además, sea:
\[ \widehat{Var}(\hat{\boldsymbol{\beta}}) = \begin{pmatrix} \widehat{Var}(\hat{\beta}_1) & \widehat{Cov}(\hat{\beta}_1,\hat{\beta}_2)\\ \widehat{Cov}(\hat{\beta}_1,\hat{\beta}_2) & \widehat{Var}(\hat{\beta}_2) \end{pmatrix}. \]
Para contrastar \(H_0:\beta_2=0\), solo se necesita:
\[ \widehat{Var}(\hat{\beta}_2). \]
El error estándar estimado de \(\hat{\beta}_2\) es:
\[ EE(\hat{\beta}_2) = \sqrt{\widehat{Var}(\hat{\beta}_2)}. \]
Solución: estadístico de Wald
El estadístico de Wald para contrastar:
\[ H_0:\beta_2=0 \]
es:
\[ Z_W = \frac{\hat{\beta}_2-0} {EE(\hat{\beta}_2)}. \]
Por tanto:
\[ Z_W = \frac{\hat{\beta}_2} {\sqrt{\widehat{Var}(\hat{\beta}_2)}}. \]
Bajo \(H_0\):
\[ Z_W \overset{aprox}{\sim} N(0,1). \]
Solución: regla de decisión bilateral
Para una prueba bilateral con nivel de significancia \(\alpha\):
\[ H_0:\beta_2=0 \]
frente a:
\[ H_1:\beta_2\neq 0, \]
se rechaza \(H_0\) si:
\[ |Z_W|>z_{1-\alpha/2}. \]
Equivalentemente, usando el valor-p:
\[ p\text{-valor} = 2P(Z>|Z_W|). \]
Se rechaza \(H_0\) si:
\[ p\text{-valor}<\alpha. \]
Solución: forma chi-cuadrado
También puede escribirse el contraste de Wald como:
\[ W = Z_W^2. \]
Entonces:
\[ W = \left( \frac{\hat{\beta}_2} {\sqrt{\widehat{Var}(\hat{\beta}_2)}} \right)^2. \]
Por tanto:
\[ W = \frac{\hat{\beta}_2^2} {\widehat{Var}(\hat{\beta}_2)}. \]
Bajo \(H_0\):
\[ W \overset{aprox}{\sim} \chi^2_1. \]
Solución: interpretación en el modelo
El modelo es:
\[ \log(\lambda_i)=\beta_1+\beta_2x_i. \]
Si:
\[ \beta_2=0, \]
entonces:
\[ \log(\lambda_i)=\beta_1. \]
Por tanto, bajo \(H_0\), la media de la respuesta no cambia con \(x_i\).
En cambio, si:
\[ \beta_2\neq 0, \]
entonces \(x_i\) tiene efecto sobre la media esperada de la respuesta.
Resumen del contraste
Para contrastar:
\[ H_0:\beta_2=0 \]
se usa:
\[ Z_W = \frac{\hat{\beta}_2} {\sqrt{\widehat{Var}(\hat{\beta}_2)}}. \]
Bajo \(H_0\):
\[ Z_W \overset{aprox}{\sim} N(0,1). \]
Equivalente:
\[ W = \frac{\hat{\beta}_2^2} {\widehat{Var}(\hat{\beta}_2)} \overset{aprox}{\sim} \chi^2_1. \]
Si se rechaza \(H_0\), se concluye que \(x_i\) está asociado con la media esperada de la respuesta en el modelo Poisson.
Referencias
Agresti, A. (2015). Foundations of linear and generalized linear models. Wiley.
Blitzstein, J. K., & Hwang, J. (2019). Introduction to probability (2nd ed.). Chapman; Hall/CRC.
Casella, G., & Berger, R. L. (2002). Statistical inference (2nd ed.). Duxbury.
DeGroot, M. H., & Schervish, M. J. (2012). Probability and statistics (4th ed.). Pearson.
Dobson, A. J., & Barnett, A. G. (2018). An introduction to generalized linear models (4th ed.). Chapman; Hall/CRC.
Faraway, J. J. (2016). Extending the linear model with R: Generalized linear, mixed effects and nonparametric regression models (2nd ed.). Chapman; Hall/CRC.
Hogg, R. V., McKean, J. W., & Craig, A. T. (2019). Introduction to mathematical statistics (8th ed.). Pearson.
Larsen, R. J., & Marx, M. L. (2008). An introduction to mathematical statistics and its applications (4th ed.). Pearson.
McCullagh, P., & Nelder, J. A. (1989). Generalized linear models (2nd ed.). Chapman; Hall.
Pawitan, Y. (2001). In all likelihood: Statistical modelling and inference using likelihood. Oxford University Press.
Pitman, J. (1993). Probability. Springer.
Rice, J. A. (2006). Mathematical statistics and data analysis (3rd ed.). Duxbury Press.
Ross, S. (2014). A first course in probability (9th ed.). Pearson.
Wackerly, D. D., Mendenhall, W., & Scheaffer, R. L. (2008). Mathematical statistics with applications (7th ed.). Thomson Brooks/Cole.
Weisberg, S. (2014). Applied linear regression (4th ed.). Wiley.
