Cuando el parámetro de dispersión \(\phi\) se asume conocido, las pruebas de hipótesis sobre un único coeficiente en un GLM pueden realizarse mediante el estadístico de Wald. Bajo muestras suficientemente grandes, el estimador \(\hat\beta_j\) es aproximadamente normal, de modo que para contrastar
\[ H_0:\;\beta_j=\beta_j^0 \]
se emplea
\[ Z=\frac{\hat{\beta}_{j}-\beta_{j}^{0}}{\mathrm{se}(\hat{\beta}_{j})}, \]
donde \(\mathrm{se}(\hat\beta_j)=\sqrt{\phi},v_j\) surge de la matriz de información. Si \(H_0\) es cierta, entonces \(Z\sim N(0,1)\) asintóticamente, y el \(p\)-valor de dos colas se calcula como \(2\Phi(-|Z|)\).
En R, tras ajustar el modelo con glm(), la función summary() muestra para cada coeficiente su estimación, error estándar, estadístico \(z\) y \(p\)-valor. Por ejemplo, usando los datos de los “noisy miners”:
library(GLMsData)
data(nminer)
nm.m1 <- glm(Minerab ~ Eucs, data = nminer, family = poisson(link = "log"))
printCoefmat(coef(summary(nm.m1))) Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.876211 0.282793 -3.0984 0.001946 **
Eucs 0.113981 0.012431 9.1691 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Estos resultados indican que tanto el intercepto como el coeficiente de Eucs son significativamente distintos de cero.
Para obtener intervalos de confianza de un solo coeficiente en un GLM con \(\phi\) conocido, se pueden emplear las estadísticas de Wald, de score o de razón de verosimilitud, tal como se describió en la secciones anteriores. En la práctica, el intervalo de Wald es el más común porque el ajuste ya proporciona la estimación \(\hat{\beta}_{j}\) y su error estándar sin cálculos adicionales.
El intervalo de confianza aproximado al \(100(1-\alpha)%\) para \(\beta_{j}\) es
\[ \hat{\beta}_{j} \pm z_{\alpha/2}^{*}\,\mathrm{se}(\hat{\beta}_{j}) \]
donde \(z_{\alpha/2}^{*}\) es el cuantíl de la normal estándar con área \(\alpha/2\) en cada cola, y \(\mathrm{se}(\hat{\beta}_{j})\) es el error estándar obtenido del inverso de la matriz de información. Este intervalo es simétrico en la escala de la predicción lineal \(\eta\).
En R, la función confint() calcula directamente estos intervalos de un objeto glm. Por ejemplo, para el modelo ajustado a los datos de noisy miners:
confint(nm.m1)Waiting for profiling to be done... 2.5 % 97.5 %
(Intercept) -1.45700887 -0.3465538
Eucs 0.08985068 0.1386685Estos resultados muestran que ambos coeficientes excluyen el cero al 95 %, reforzando la evidencia de su relevancia en el modelo.
Al ajustar un GLM obtenemos estimaciones de los parámetros \(\hat{\beta}_j\), y de ahí derivamos predicciones de la escala del predictor lineal \(\hat\eta = g(\hat\mu)\). Dado que \(\hat\eta\) tiene incertidumbre, sus transformaciones a \(\hat\mu = g^{-1}(\hat\eta)\) también la tendrán. Cuando \(\phi\) es conocido, un intervalo de confianza aproximado al \(100(1-\alpha)%\) para \(\eta\) viene dado por
\[ \hat{\eta}\pm z_{\alpha/2}^{*}\,\mathrm{se}(\hat{\eta}), \]
donde \(\mathrm{se}(\hat{\eta})=\sqrt{\operatorname{var}[\hat{\eta}]}\) y \(z_{\alpha/2}^{*}\) es el cuantil de la distribución normal estándar con área \(\alpha/2\) en cada cola. Después, aplicando la inversa del enlace \(\mu=g^{-1}(\eta)\) a los límites de este intervalo, se obtiene el intervalo de confianza para \(\mu\). Este intervalo es simétrico en la escala \(\eta\) pero, al transformar de vuelta, suele resultar asimétrico en la escala de \(\mu\).
En R, la función predict() permite calcular directamente estas predicciones y sus errores estándar. Por ejemplo, para el modelo de los noisy miners (datos nminer), la predicción de \(\hat\mu\) cuando Eucs = 10 y la construcción de su intervalo de confianza al 95 % se harían así:
# Predicción del predictor lineal con error estándar
out <- predict(nm.m1,
newdata = data.frame(Eucs = 10),
se.fit = TRUE)
# Predicción en la escala de la respuesta
out2 <- predict(nm.m1,
newdata = data.frame(Eucs = 10),
se.fit = TRUE,
type = "response")
# Confirmar que ambas predicciones coinciden
c(exp(out$fit), out2$fit) 1 1
1.30161 1.30161 # Construcción del intervalo de confianza para mu
zstar <- qnorm(0.975)
ci.lo <- exp(out$fit - zstar * out$se.fit)
ci.hi <- exp(out$fit + zstar * out$se.fit)
c(Lower = ci.lo, Estimate = exp(out$fit), Upper = ci.hi) Lower.1 Estimate.1 Upper.1
0.924013 1.301610 1.833512 Esto arroja una estimación \(\hat\mu\approx1.302\) con intervalo de confianza 95 % aproximadamente \([0.924,,1.834]\). Estos límites no son simétricos respecto a \(\hat\mu\) debido a la transformación inversa del enlace logarítmico y al hecho de que, para una Poisson, \(V(\mu)=\mu\), lo que hace que los intervalos se ensanchen conforme aumenta \(\hat\mu\).
Cuando dos GLM basados en la misma familia EDM difieren sólo en sus componentes sistemáticos, uno puede ser anidado en el otro si el modelo más simple (A) se obtiene imponiendo ciertas restricciones de cero sobre los coeficientes del modelo más completo (B). Específicamente, si
\[ \begin{aligned} \text{Modelo A:}\quad & g(\hat\mu_{A})=\hat\beta_{0}+\hat\beta_{1}x_{1}+\cdots+\hat\beta_{p_A}x_{p_A},\\ \text{Modelo B:}\quad & g(\hat\mu_{B})=\hat\beta_{0}+\hat\beta_{1}x_{1}+\cdots+\hat\beta_{p_A}x_{p_A}+\cdots+\hat\beta_{p_B}x_{p_B}, \end{aligned} \]
entonces el contraste
\[ H_{0}:\ \beta_{p_A+1}=\cdots=\beta_{p_B}=0 \]
se puede probar mediante la estadística de razón de verosimilitudes. Usando la devianza, si \(\phi\) es conocido, el estadístico es
\[ L \;=\; 2\bigl\{\ell_B - \ell_A\bigr\} \;=\; \frac{D\bigl(y,\hat\mu_{A}\bigr)\;-\;D\bigl(y,\hat\mu_{B}\bigr)}{\phi}, \]
donde \(D(y,\hat\mu)\) es la devianza residual del modelo. Bajo \(H_{0}\) y muestras grandes,
\[ L \sim \chi^2_{\,p'_B-p'_A}. \]
Cuando los modelos sólo difieren en un coeficiente (\(p'_B-p'_A=1\)), puede definirse la estadística “signada”:
\[ Z=\operatorname{sign}\bigl(\hat\beta_{p_B}\bigr)\,\sqrt{L}, \]
que satisface \(Z\sim N(0,1)\) bajo \(H_{0}\) y permite probar hipótesis unilaterales.
Ejemplo
# Modelo nulo (sólo intercepto) y modelo con Eucs
nm.m0 <- glm(Minerab ~ 1, data=nminer, family=poisson)
nm.m1 <- glm(Minerab ~ Eucs, data=nminer, family=poisson)
# Deviancias y grados de libertad residuales
dev0 <- deviance(nm.m0); df0 <- df.residual(nm.m0)
dev1 <- deviance(nm.m1); df1 <- df.residual(nm.m1)
# Estadística de razón de verosimilitudes y p-valor
L <- dev0 - dev1
p <- pchisq(L, df0 - df1, lower.tail=FALSE)
c(L = L, p.value = p) L p.value
8.722735e+01 9.673697e-21 El resultado muestra un valor de \(L\approx87.23\) y \(p\text{-value}\approx9.7\times10^{-21}\), indicando que la inclusión de Eucs mejora significativamente el ajuste frente al modelo sólo con intercepto.
Cuando comparamos una secuencia de modelos anidados añadiendo una variable explicativa tras otra, cada contraste mediante prueba de razón de verosimilitudes equivale a la diferencia en devianza residual entre dos modelos sucesivos. Estos cálculos se estructuran en una tabla de análisis de devianza, que extiende directamente las tablas de ANOVA de los modelos lineales al caso de los GLM.
En R, la función anova() aplicada a un objeto glm produce esta tabla. Para obtener los valores p basados en \(\chi^{2}\) hay que especificar test="Chisq". Si la dispersión \(\phi\) fuese distinta de 1, puede indicarse con el argumento dispersion=.
# Supongamos nm.m1 es el GLM que incluye Eucs y nm.m0 el modelo nulo
anova(nm.m1, test="Chisq")Analysis of Deviance Table
Model: poisson, link: log
Response: Minerab
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 30 150.545
Eucs 1 87.227 29 63.318 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En esta tabla, la fila Eucs muestra que añadir el número de eucaliptos reduce la devianza en \(L=87.23\) con 1 grado de libertad, lo que corresponde a un \(P<0.001\) según la distribución \(\chi^{2}_{1}\). El valor de la devianza residual y los grados de libertad residuales aparecen también para cada modelo, facilitando la comparación directa.
Las pruebas de score permiten evaluar si un nuevo predictor podría mejorar un modelo GLM ya ajustado, contrastando \(H_{0}\colon \beta_{p+1}=0\) sin necesidad de refit completo. Para ello se aprovecha la relación con los residuos de Pearson y la interpretación del GLM como regresión ponderada.
Sea \(e(y)_i\) el residuo de trabajo del modelo actual (definido en (6.20)) y \(e\left(x_{p+1}\right)_{i}\) el residuo de la regresión por mínimos cuadrados ponderados de \(x_{p+1}\) sobre los predictores existentes, con pesos \(W_{i}\). El estadístico de score es
\[ Z=\frac{\sum_{i=1}^{n} e(x_{p+1})_{i}\,e(y)_{i}} {\bigl\{\sum_{i=1}^{n}e(x_{p+1})_{i}^{2}\bigr\}^{1/2}}. \]
Bajo \(H_{0}\), \(Z\sim N(0,1)\) aproximadamente. En R, la función glm.scoretest() del paquete statmod calcula directamente este \(Z\) para un GLM ya ajustado.
library(statmod)
# Modelo nulo sin Eucs
nm.m0 <- glm(Minerab ~ 1, data = nminer, family = poisson)
# Estadístico de score para añadir Eucs
z.stat <- glm.scoretest(nm.m0, nminer$Eucs)
p.val <- 2 * pnorm(abs(z.stat), lower.tail = FALSE)
round(c(score.stat = z.stat, P = p.val), 4)score.stat P
9.7565 0.0000 En este ejemplo, el valor de \(Z\approx9.76\) y \(P\approx0\) indican que Eucs debería incorporarse al modelo.
Además, el tradicional test de independencia de Pearson en tablas de contingencia es un caso especial de prueba score para un log-lineal Poisson con interacción. Por ejemplo:
# Tabla 2×2 de conteos
Y <- matrix(c(10,20,20,10), 2, 2,
dimnames = list(c("A1","A2"), c("B1","B2")))
# Test de chi-cuadrado de Pearson
chisq.test(Y, correct = FALSE)$p.value[1] 0.009823275# Misma prueba como score test de interacción en GLM
y <- as.vector(Y)
A <- factor(rep(1:2, each = 2))
B <- factor(rep(1:2, times = 2))
fit <- glm(y ~ A + B, family = poisson)
# x2 codifica sólo el efecto A2:B2
z.stat2 <- glm.scoretest(fit, x2 = c(0,0,0,1))
2 * pnorm(-abs(z.stat2))[1] 0.009823231Aquí se confirma que la prueba de score reproduce el \(P\)-valor de Pearson para testar la interacción.
En inferencia basada en verosimilitud, los resultados de distribución para los estadísticos de prueba (Wald, score y razón de verosimilitudes) son aproximaciones asintóticas que mejoran a medida que el número de observaciones \(n\) crece. Estas asintóticas de muestra grande suelen ser fiables para las pruebas de score y de razón de verosimilitudes incluso con muestras moderadas. En cambio, las pruebas de Wald—y en particular para EDM binomiales con tamaños pequeños de \(m\) requieren muestras más grandes para no resultar demasiado conservadoras. Cuando \(n\) es suficientemente grande para que \(\mathrm{se}(\hat\beta_j)\) sea pequeña, las aproximaciones asintóticas suelen ser razonablemente precisas. Cabe recordar que, para GLM normales lineales, todos estos resultados son exactos.
Las pruebas de bondad de ajuste permiten verificar si el modelo actual capta toda la variabilidad sistemática de los datos o si faltan variables explicativas importantes. Para ello se compara:
La prueba utiliza el estadístico de razón de verosimilitudes, que en forma de devianza queda
\[ L=2\{\ell_B-\ell_A\}=\frac{D(y,\hat{\mu}_A)-D(y,\hat{\mu}_B)}{\phi}\,, \]
y como \(D(y,\hat{\mu}_B)=0\) (modelo saturado),
\[ L=\frac{D(y,\hat{\mu}_A)}{\phi}\sim \chi^2_{\,n-p'}\,. \]
Bajo la hipótesis nula de buen ajuste, \(L\) sigue aproximadamente una \(\chi^2_{n-p'}\). Rechazar \(H_0\) indica que el modelo actual no es suficiente para explicar todos los patrones sistemáticos y que conviene incorporar nuevas variables.
En la comparación del modelo actual frente al modelo saturado, la devianza residual del modelo saturado es cero, de modo que el estadístico de razón de verosimilitudes cobra la forma muy sencilla
\[ L \;=\; \frac{D(y,\hat{\mu}_\text{actual}) - D(y,\hat{\mu}_\text{saturado})}{\phi} \;=\; \frac{D(y,\hat{\mu})}{\phi}\,. \]
Como \(D(y,\hat{\mu}_\text{saturado})=0\), la prueba reduce a evaluar directamente \(D(y,\hat{\mu})\).
Aunque a primera vista podría parecer razonable comparar \(D(y,\hat{\mu})\) con una \(\chi^{2}_{n-p'}\), esta asintótica clásica no se cumple porque el número de parámetros del modelo saturado crece con \(n\). Para justificar la aproximación de la distribución de la devianza en este caso se recurre a la aproximación de saddlepoint (más detalles adelante).
Un ejemplo muy conocido de esta prueba de bondad de ajuste es el test \(G\) de independencia en tablas de contingencia, que utiliza precisamente la devianza como estadístico de contraste.
La prueba de bondad de ajuste basada en el estadístico de Pearson surge como un test de score para comparar el modelo actual con el modelo saturado. En lugar de la devianza, el estadístico es
\[ X^{2}=\sum_{i=1}^{n}W_{i}\bigl(z_{i}-\hat{\eta}_{i}\bigr)^{2} \quad\Longleftrightarrow\quad X^{2}=\sum_{i=1}^{n}\frac{w_{i}\,(y_{i}-\hat{\mu}_{i})^{2}}{V(\hat{\mu}_{i})}\,, \]
que en R se obtiene directamente con
pearson.gof <- sum(fit$weights * fit$residuals^2)Siguiendo la teoría de tests de score, cabría comparar \(X^{2}\) con una \(\chi^{2}_{n-p'}\), pero esa asintótica clásica no se cumple cuando el modelo saturado crece con \(n\) por el mismo motivo anterior. Para fundamentar la aproximación de la distribución de Pearson en este caso se recurre al Teorema Central del Límite (ver sección siguiente).
Ejemplos En la investigación genética molecular moderna, es habitual estudiar ratones transgénicos que presentan mutaciones en un gen específico pero que, por lo demás, son idénticos a los ratones normales. En un estudio del Walter and Eliza Hall Institute of Medical Research (Melbourne), se cruzaron varios ratones heterocigotos (con un alelo normal \(A\) y un alelo mutante \(a\) para el gen de interés). La herencia mendeliana simple implicaría que las genotipos \(AA\) (normal), \(Aa\) (heterocigoto mutante) y \(aa\) (homocigoto mutante) deberían aparecer en las crías en proporciones \(1/4\), \(1/2\) y \(1/4\), respectivamente. En un experimento concreto, se obtuvieron los números de crías que se muestran en la siguiente tabla.
¿Son estos recuentos compatibles con la herencia mendeliana? Respondemos a esta pregunta ajustando un GLM de Poisson cuyos valores ajustados sigan las proporciones mendelianas:
y <- c(15, 26, 4)
x <- c(1/4, 1/2, 1/4)
fit <- glm(y ~ 0 + x, family=poisson)
pearson.gof <- sum(fit$weights * fit$residuals^2)
deviance.gof <- fit$deviance
df <- fit$df.residual
p.values <- pchisq(c(deviance.gof, pearson.gof), df, lower.tail=FALSE)Tanto la devianza como el estadístico de Pearson arrojan \(p<0.01\), rechazando la adequación del modelo mendeliano y sugiriendo un déficit de homocigotos mutantes (posible letalidad prenatal).
En contraste con la teoría asintótica de grandes muestras, que requiere un \(n\) grande respecto al número de parámetros, las pruebas de bondad de ajuste exigen resultados que funcionen bien a nivel de cada observación. Para ello se emplean asintóticos de pequeña dispersión, basados en dos herramientas principales: la aproximación de saddlepoint para las estadísticas de devianza y el Teorema del Límite Central para las estadísticas de Pearson.
La aproximación de saddlepoint es sorprendentemente precisa incluso en colas extremas y da lugar a que el estadístico de devianza residual
\[ \frac{D(y,\hat{\mu})}{\phi}\sim \chi_{n-p^{\prime}}^{2} \]
sea una buena aproximación, siempre que se cumpla el criterio de que
\[ \tau = \frac{\phi\,V(y)}{(y-\text{límite})^{2}} \le 1/3 \]
para todas las observaciones y límite: borde del soporte \(S\) para \(y\). De este modo se garantiza que tanto los \(y_i\) como los \(\hat\mu_i\) estén en rango donde la chi-cuadrado aplique.
El Teorema del Límite Central, de convergencia más lenta \(O(\phi^{1/2})\), justifica que la estadística de Pearson
\[ \frac{X^{2}}{\phi} \sim \chi_{n-p^{\prime}}^{2},\quad X^{2}=\sum_{i=1}^n \frac{w_i\,(y_i-\hat\mu_i)^2}{V(\hat\mu_i)}, \]
también siga una chi-cuadrado asintótica siempre que se refuerce el criterio a \(\tau\le1/5\) para cada observación.
Para las familias EDM más comunes, estas reglas prácticas son:
En los casos extremos de \(\phi\to0\), devianza y Pearson se hacen casi idénticos. Para GLMs normales (lineales) ambas aproximaciones son exactas.
Ejemplo En el ejemplo de la herencia mendeliana, el recuento mínimo es \(y_i=4\), por lo que se satisface el criterio de saddlepoint pero no el de CLT. Así, la prueba de devianza es más fiable que la de Pearson.
Ejemplo En los datos de los “noisy miners”, hay múltiples ceros, por lo que no se cumple \(\min\{y_i\}\ge3\) y ninguna de las pruebas de bondad de ajuste (devianza o Pearson) es confiable para ese conjunto de datos.
Cuando la dispersión \(\phi\) no es conocida, los tests de Wald se construyen de forma semejante al caso con \(\phi\) conocida, pero incorporando un estimador \(s^2\) de \(\phi\) en el cómputo de los errores estándar. El estadístico de Wald para contrastar \(H_{0}: \beta_{j}=\beta_{j}^{0}\) pasa a ser
\[ T=\frac{\hat{\beta}_{j}-\beta_{j}^{0}}{\operatorname{se}(\hat{\beta}_{j})}, \]
donde ahora \(\operatorname{se}(\hat{\beta}_{j})=sv_{j}\), siendo los \(v_{j}\) las raíces de los elementos diagonales de la inversa de la información. En R, el estimador por defecto de \(\phi\) es el estadístico de Pearson \(\bar\phi=X^2/(n-p')\), de modo que el summary() de un objeto glm muestra un estadístico etiquetado como “t” y calcula los \(p\)-valores mediante la distribución \(t_{n-p'}\). Para muestras muy grandes, \(T\) se aproxima a una \(N(0,1)\), recuperándose en el caso normal-lineal la exactitud del \(t\)-test.
Ejemplo
Para los datos de los cerezos (trees), se ajusta un GLM gamma con enlace logarítmico:
data(trees)
tr.m2 <- glm( Volume ~ log(Girth) + log(Height),
family=Gamma(link="log"), data=trees )
printCoefmat(coef(summary(tr.m2))) Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.69111 0.78784 -8.4929 3.108e-09 ***
log(Girth) 1.98041 0.07389 26.8021 < 2.2e-16 ***
log(Height) 1.13288 0.20138 5.6255 5.037e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La salida muestra coeficientes, errores estándar y estadísticos \(t\) con \(n-p'=28\) grados de libertad, indicando que tanto \(\log(\mathrm{Girth})\) como \(\log(\mathrm{Height})\) son significativos. Los dos estimadores de dispersión son muy similares:
phi.meandev <- deviance(tr.m2) / df.residual(tr.m2)
phi.pearson <- summary(tr.m2)$dispersion
c(Mean.deviance=phi.meandev, Pearson=phi.pearson)Mean.deviance Pearson
0.006554117 0.006427286 Si se fuerza el uso del estimador de media de la devianza para los errores estándar:
printCoefmat(coef(summary(tr.m2, dispersion=phi.meandev))) Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.691109 0.795578 -8.4104 < 2.2e-16 ***
log(Girth) 1.980412 0.074616 26.5415 < 2.2e-16 ***
log(Height) 1.132878 0.203361 5.5708 2.536e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1R realiza entonces z-tests tratando la dispersión como conocida, lo que tiende a sobrestimar ligeramente la significación de los coeficientes.
Comentario: Cuando \(\phi\) es muy pequeña, los distintos estimadores (\(\bar\phi\), media de la devianza) coinciden y el efecto de referir a \(t\) frente a \(N(0,1)\) es mínimo; sin embargo, para muestras moderadas conviene mantener el \(t\)-test para reflejar la incertidumbre en \(\phi\).
Cuando la dispersión \(\phi\) es desconocida, los intervalos de confianza basados en Wald se construyen de forma análoga al caso con \(\phi\) conocida, pero empleando un estimador de \(\phi\) (por defecto, el estimador de Pearson) en el cálculo de los errores estándar. El intervalo de confianza del \(100(1-\alpha)%\) para \(\beta_{j}\) es
\[ \hat{\beta}_{j} \pm t_{\alpha / 2,\,n-p^{\prime}}^{*}\,\operatorname{se}\bigl(\hat{\beta}_{j}\bigr), \]
donde \(t_{\alpha / 2,\,n-p^{\prime}}^{*}\) es el cuantil de la \(t\) de Student con \(n-p^{\prime}\) grados de libertad que deja un área \(\alpha/2\) en cada cola, y \(\operatorname{se}(\hat{\beta}_{j})\) incorpora la estimación de la dispersión. Este enfoque coincide exactamente con el caso normal lineal, donde los intervalos son exactos.
En R, la función confint() aplicada a un objeto glm devuelve automáticamente estos intervalos basados en la distribución \(t\):
confint(tr.m2)Waiting for profiling to be done... 2.5 % 97.5 %
(Intercept) -8.2358004 -5.139294
log(Girth) 1.8359439 2.124974
log(Height) 0.7364235 1.528266En el ejemplo de los cerezos (trees), se observa que el intervalo para \(\log(\mathrm{Girth})\) contiene el valor teórico \(\beta_{1}=2\), mientras que el de \(\log(\mathrm{Height})\) contiene \(\beta_{2}=1\), corroborando las predicciones del modelo geométrico.
Cuando la dispersión \(\phi\) es desconocida, la incertidumbre en las predicciones de la media \(\hat{\mu}\) procede tanto de la estimación de los coeficientes como de \(\phi\). Para construir intervalos de confianza utilizamos primero la escala del predictor lineal \(\hat{\eta}=g(\hat{\mu})\), cuya varianza estimada \(\widehat{\operatorname{var}}[\hat{\eta}]\) se obtiene directamente de la matriz de información. El intervalo de confianza del \(100(1-\alpha)%\) para \(\eta\) es
\[ \hat{\eta}\pm t_{\alpha/2,\,n-p'}^{*}\,\operatorname{se}\bigl(\hat{\eta}\bigr), \]
donde \(t_{\alpha/2,,n-p'}^{*}\) es el cuantil de la \(t\) de Student con \(n-p'\) grados de libertad y \(\operatorname{se}(\hat{\eta})=\sqrt{\widehat{\operatorname{var}}[\hat{\eta}]}\). Para obtener el intervalo de confianza en la escala original de la media, se aplica la función inversa del enlace:
\[ \bigl[g^{-1}(\hat{\eta}-t^*\operatorname{se}),\,g^{-1}(\hat{\eta}+t^*\operatorname{se})\bigr]. \]
En R, predict(..., se.fit=TRUE) devuelve el ajuste y su error estándar en la escala de \(\eta\), mientras que type="response" transforma automáticamente ambos al nivel de \(\mu\).
Ejemplo
Para los datos de los cerezos (altura 70 ft, circunferencia 15 in), calculamos primero la predicción y su error estándar:
out <- predict(tr.m2,
newdata = data.frame(Height = 70, Girth = 15),
se.fit = TRUE)A continuación formamos el intervalo de confianza del 95 % en la escala de \(\mu\):
tstar <- qt(0.975, df = df.residual(tr.m2))
ci.lo <- exp(out$fit - tstar * out$se.fit)
ci.hi <- exp(out$fit + tstar * out$se.fit)
c(Lower = ci.lo, Estimate = exp(out$fit), Upper = ci.hi) Lower.1 Estimate.1 Upper.1
30.81902 32.62157 34.52955 Este resultado indica que el volumen medio estimado es \(\hat{\mu}=32.62\) ft³, con intervalo \((30.82,;34.53)\) ft³.
Podemos extenderlo a varias alturas manteniendo la circunferencia fija en 15 in:
newHt <- seq(min(trees$Height), max(trees$Height), by = 4)
newVol <- predict(tr.m2, se.fit = TRUE,
newdata = data.frame(Height = newHt, Girth = 15))
ci.lo <- exp(newVol$fit - tstar * newVol$se.fit)
ci.hi <- exp(newVol$fit + tstar * newVol$se.fit)
cbind(newHt, Lower = ci.lo, Estimate = exp(newVol$fit), Upper = ci.hi) newHt Lower Estimate Upper
1 63 26.33168 28.95124 31.83141
2 67 28.88896 31.04230 33.35614
3 71 31.45834 33.15002 34.93267
4 75 33.93192 35.27358 36.66829
5 79 36.10127 37.41225 38.77084
6 83 37.87594 39.56537 41.33016
7 87 39.40973 41.73232 44.19180De este modo obtenemos un conjunto de intervalos que reflejan cómo varía la incertidumbre de \(\mu\) con la altura, respetando siempre la simetría en la escala de \(\eta\).
Cuando la dispersión \(\phi\) no es conocida, la comparación de dos modelos anidados A y B (con \(p_{A}'\) y \(p_{B}'\) parámetros respectivamente) se basa en la diferencia de sus deviancias escaladas por un estimador de \(\phi\). El estadístico general es
\[ F=\frac{\{D(y,\hat\mu_A)-D(y,\hat\mu_B)\}/(p_B'-p_A')}{s^2}\,, \]
donde \(s^2\) es una estimación apropiada de \(\phi\) obtenida del modelo B. Este enfoque refleja exactamente el caso clásico de regresión lineal y, bajo el supuesto de normalidad y con \(s^2\) no sesgado, sigue una distribución \(F_{,p_B'-p_A',,n-p_B'}\).
Para otros GLM se usan tres variantes de \(F\) dependiendo del estimador de \(\phi\):
En la práctica, \(\tilde F\) suele dar resultados razonablemente aproximados a la \(F\) teórica.
\[ \begin{align*} \hat{F}^{0} & =\frac{\left\{D\left(y,\hat{\mu}_{A}\right)-D\left(y,\hat{\mu}_{B}\right)\right\}/\left(p_{B}^{\prime}-p_{A}^{\prime}\right)}{\hat{\phi}_{B}^{0}} \\ \bar{F} & =\frac{\left\{D\left(y,\hat{\mu}_{A}\right)-D\left(y,\hat{\mu}_{B}\right)\right\}/\left(p_{B}^{\prime}-p_{A}^{\prime}\right)}{\bar{\phi}_{B}} \\ \tilde{F} & =\frac{\left\{D\left(y,\hat{\mu}_{A}\right)-D\left(y,\hat{\mu}_{B}\right)\right\}/\left(p_{B}^{\prime}-p_{A}^{\prime}\right)}{\tilde{\phi}_{B}} \end{align*} \]
donde todas las estimaciones de \(\phi\) se basan en el modelo B. Como es habitual, las tres estadísticas \(F\) coinciden en el caso de regresión lineal, y en ese escenario siguen exactamente una distribución \(F\bigl(p_{B}^{\prime}-p_{A}^{\prime},,n-p_{B}^{\prime}\bigr)\) bajo la hipótesis nula de igualdad entre los modelos A y B. Para otros GLM, las estadísticas \(F\) se aproximan asintóticamente a esa misma distribución bajo la hipótesis nula. La aproximación suele ser buena siempre que el denominador del estadístico \(F\) siga aproximadamente una chi-cuadrado escalada.
Cuando sólo se añade un coeficiente (\(p_B'=p_A'+1\)), puede definirse un estadístico firmado
\[ t=\operatorname{sign}(\hat\beta_{p_B})\,\sqrt{F} \]
que, bajo \(H_0:\beta_{p_B}=0\), es aproximadamente \(t_{n-p_B'}\) y sirve para pruebas unilaterales.
Ejemplo (cerezos). Consideremos los datos de los cerezos (trees) y el modelo tr.m2 ajustado. Ajustamos secuencialmente las dos variables explicativas log(Girth) y log(Height) y anotamos la devianza residual y los grados de libertad residuales para cada modelo:
data(trees)
tr.m0 <- glm( Volume ~ 1, family=Gamma(link="log"), data=trees)
tr.m1 <- update(tr.m0, . ~ . + log(Girth) )
tr.m2 <- update(tr.m1, . ~ . + log(Height) )
c( deviance(tr.m0), deviance(tr.m1), deviance(tr.m2) )[1] 8.3172012 0.3840839 0.1835153c( df.residual(tr.m0), df.residual(tr.m1), df.residual(tr.m2) )[1] 30 29 28A continuación, calculamos las diferencias de devianza entre los modelos y los correspondientes cambios en los grados de libertad:
dev1 <- deviance(tr.m0) - deviance(tr.m1)
dev2 <- deviance(tr.m1) - deviance(tr.m2)
df1 <- df.residual(tr.m0) - df.residual(tr.m1)
df2 <- df.residual(tr.m1) - df.residual(tr.m2)
c( dev1, dev2 )[1] 7.9331173 0.2005686c( df1, df2 )[1] 1 1Para calcular las estadísticas de prueba \(F\) según , necesitamos primero estimar \(\phi\):
phi.meandev <- deviance(tr.m2) / df.residual(tr.m2) # Estimador media de devianza
phi.Pearson <- summary(tr.m2)$dispersion # Estimador de Pearson
c("Mean deviance" = phi.meandev, "Pearson" = phi.Pearson )Mean deviance Pearson
0.006554117 0.006427286 Los estimadores de Pearson y de media de devianza resultan muy similares. De igual modo, las estadísticas \(F\) y los \(P\)-valores correspondientes, calculados con ambos estimadores, son comparables:
F.Pearson <- c(dev1/df1, dev2/df2) / phi.Pearson
F.meandev <- c(dev1/df1, dev2/df2) / phi.meandev
P.Pearson <- pf(F.Pearson, df1, df.residual(tr.m2), lower.tail=FALSE)
P.meandev <- pf(F.meandev, df2, df.residual(tr.m2), lower.tail=FALSE)
tab <- data.frame(F.Pearson, P.Pearson, F.meandev, P.meandev)
rownames(tab) <- c("Girth", "Height")
print(tab, digits=3) F.Pearson P.Pearson F.meandev P.meandev
Girth 1234.3 1.05e-24 1210.4 1.38e-24
Height 31.2 5.60e-06 30.6 6.50e-06Estos resultados muestran que log(Girth) es significativo en el modelo y que, una vez ajustado por log(Girth), también log(Height) resulta significativo.
Cuando ajustamos una serie de GLMs anidados, registramos para cada modelo la devianza residual \(D(y,\hat\mu)=\sum_i w_i\,d(y_i,\hat\mu_i)\) y los grados de libertad residuales \(n-p'\). Las diferencias sucesivas de estas devianzas y de los grados de libertad se organizan en una tabla de análisis de devianza, que generaliza la tabla ANOVA de regresión lineal.
En R, se obtiene con
anova(fit, test="F")donde test="F" pide que las diferencias de devianza se evalúen frente a una distribución \(F\) con \((\Delta p',,n-p'_B)\) grados de libertad, usando por defecto el estimador de Pearson \(\bar\phi\). Si se prefiere otro estimador de la dispersión (por ejemplo la media de devianza \(\tilde\phi\)), basta con añadir
anova(fit, test="F", dispersion=phi_meandev)El orden de entrada de los términos en el modelo influye en los valores de \(F\) y sus \(P\)–valores, aunque la conclusión global (qué variables resultan significativas) suele mantenerse.
Ejemplo. Con los datos trees y el modelo
tr.m2 <- glm( Volume ~ log(Girth) + log(Height),
family=Gamma(link="log"), data=trees )la tabla anova(tr.m2, test="F") muestra la contribución de log(Girth) y log(Height) al ajuste, sus estadísticas \(F\) y los \(P\)–valores asociados:
anova(tr.m2, test="F")Analysis of Deviance Table
Model: Gamma, link: log
Response: Volume
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 30 8.3172
log(Girth) 1 7.9331 29 0.3841 1234.287 < 2.2e-16 ***
log(Height) 1 0.2006 28 0.1835 31.206 5.604e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Si invertimos el orden:
tr.rev <- glm( Volume ~ log(Height) + log(Girth),
family=Gamma(link="log"), data=trees )
anova(tr.rev, test="F")Analysis of Deviance Table
Model: Gamma, link: log
Response: Volume
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 30 8.3172
log(Height) 1 3.5345 29 4.7827 549.92 < 2.2e-16 ***
log(Girth) 1 4.5992 28 0.1835 715.57 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1vemos que aunque ambos términos siguen siendo altamente significativos, los valores numéricos de \(F\) y \(P\) cambian, subrayando la interpretación condicional de cada efecto.
Los tests de score se basan en la derivada de la verosimilitud y asumen estrictamente que \(\phi\) es conocido. Sin embargo, en la práctica pueden aplicarse aproximando \(\phi\) mediante un estimador, típicamente el estimador de Pearson \(\bar\phi\). En R, la función glm.scoretest() (paquete statmod) usa por defecto \(\bar\phi\), aunque puede recibir otro valor con el argumento dispersion=. Al sustituir un estimador de \(\phi\), la estadística de score se considera aproximadamente \(t_{n-p'}\) bajo la hipótesis nula.
En el ejemplo de los datos de los cerezos (trees), se usan tests de score para evaluar si conviene añadir cada covariable al modelo que ya contiene la otra. Primero, ajustando
mA <- glm(Volume ~ log(Girth), family=Gamma(link="log"), data=trees)
t.Ht <- glm.scoretest(mA, log(trees$Height))
p.Ht <- 2 * pt( -abs(t.Ht), df=df.residual(mA) ) # Two-tailed P-value
tab <- data.frame(Score.stat = t.Ht, P.Value=p.Ht )
print(tab, digits=3) Score.stat P.Value
1 3.83 0.00063se obtiene un estadístico de score y su \(P\)-valor (tratado como \(t\) con \(n-p'\) grados de libertad) que indican si \(\log(\text{Height})\) aporta información adicional dado \(\log(\text{Girth})\). De forma análoga,
mB <- glm(Volume ~ log(Height), family=Gamma(link="log"), data=trees)
t.Girth <- glm.scoretest(mB, log(trees$Girth))
p.Girth <- 2 * pt( -abs(t.Girth), df=df.residual(mB) )
tab <- data.frame(Score.stat = t.Girth, P.Value=p.Girth )
print(tab, digits=3) Score.stat P.Value
1 5.22 1.36e-05evalúa la relevancia de \(\log(\text{Girth})\) en presencia de \(\log(\text{Height})\).
Aunque los \(P\)-valores resultantes suelen ser algo más conservadores que los de los Wald tests —debido a que la dispersión se estima en el modelo reducido—, en ambos casos los tests de score confirman que tanto \(\log(\text{Girth})\) como \(\log(\text{Height})\) son necesarios para describir adecuadamente el volumen de los árboles.
Cuando los modelos no están anidados, las pruebas de razón de verosimilitud no aplican, pero es posible comparar su capacidad predictiva usando el AIC y el BIC. Para un GLM con \(n\) observaciones, \(p'\) parámetros de regresión y \(\phi\) conocido se definen
\[ \mathrm{AIC}=-2\times\ell\bigl(\hat\beta_0,\dots,\hat\beta_p,\phi;y\bigr)+2\,p' \]
\[ \mathrm{BIC}=-2\times\ell\bigl(\hat\beta_0,\dots,\hat\beta_p,\phi;y\bigr)+(\log n)\,p' \]
Si \(\phi\) es desconocido, se sustituye por su mLE \(\hat\phi\), pasando a
\[ \mathrm{AIC}=-2\times\ell\bigl(\hat\beta_0,\dots,\hat\beta_p,\hat\phi;y\bigr)+2\,(p'+1) \]
\[ \mathrm{BIC}=-2\times\ell\bigl(\hat\beta_0,\dots,\hat\beta_p,\hat\phi;y\bigr)+(\log n)\,(p'+1). \]
En R, las funciones AIC() y BIC() usan para \(\hat\phi\) el estimador de media de la devianza \(D(y, \hat{\mu}) / n\), que coincide con el mLE en los GLM normales y de inversa gaussiana, y es aproximado en gamma cuando la aproximación de saddlepoint es válida.
Queremos comparar dos modelos no anidados que fijan teóricamente el coeficiente de \(\log(\text{Girth})\) o de \(\log(\text{Height})\) mediante un offset:
tr.aic1 <- glm( Volume ~ offset(2*log(Girth)) + log(Height),
family=Gamma(link="log"), data=trees)
tr.aic2 <- glm( Volume ~ log(Girth) + offset(log(Height)),
family=Gamma(link="log"), data=trees)
AIC(tr.aic1, tr.aic2) df AIC
tr.aic1 3 137.9780
tr.aic2 3 138.3677Ambas casos, AIC() y extractAIC(), devuelven valores numéricos (por ejemplo, ~137.98 vs ~138.37). El modelo con AIC menor se considera preferible en términos de predicción. En este caso, el AIC sugiere elegir el modelo que fija el coeficiente de \(\log(\text{Girth})\) en 2 y estima libremente el coeficiente de \(\log(\text{Height})\).
En los GLM pueden emplearse las mismas funciones automáticas que en regresión lineal—drop1(), add1() y step()—para añadir o eliminar términos basándose en el AIC. Por defecto, R selecciona el modelo que minimiza el AIC. No obstante, todas las reservas sobre las técnicas de selección automática en regresión lineal (sobreajuste, dependencia del criterio, etc.) siguen vigentes en el contexto de los GLM.
Cuando la dispersión $$ se estima (mediante la media de la devianza), el AIC calculado para cada modelo es solo aproximado, pues la estimación de $$ varía con el modelo y no coincide con el mLE. Hay que interpretar los resultados con cautela y, de ser posible, contrastarlos con criterios teóricos o validación externa.
Para ilustrar el uso de step() en el conjunto de datos trees (modelo Gamma con enlace logarítmico), primero definimos un modelo mínimo y otro máximo:
min.model <- glm(Volume ~ 1,
data = trees,
family = Gamma(link="log"))
max.model <- glm(Volume ~ log(Girth) + log(Height),
data = trees,
family = Gamma(link="log"))A continuación, aplicamos un procedimiento de eliminación hacia atrás, de adelante y ambos:
m.b <- step(max.model,
scope = list(lower = min.model, upper = max.model),
direction = "backward")Start: AIC=139.9
Volume ~ log(Girth) + log(Height)
Df Deviance AIC
<none> 0.1835 139.90
- log(Height) 1 0.3841 169.11
- log(Girth) 1 4.7827 853.48m.f <- step(min.model,
scope = list(lower = min.model, upper = max.model),
direction = "forward")Start: AIC=255.48
Volume ~ 1
Df Deviance AIC
+ log(Girth) 1 0.3841 230.75
+ log(Height) 1 4.7827 245.57
<none> 8.3172 255.48
Step: AIC=160.83
Volume ~ log(Girth)
Df Deviance AIC
+ log(Height) 1 0.18352 147.82
<none> 0.38408 160.83
Step: AIC=139.9
Volume ~ log(Girth) + log(Height)m.s <- step(min.model,
scope = list(lower = min.model, upper = max.model),
direction = "both")Start: AIC=255.48
Volume ~ 1
Df Deviance AIC
+ log(Girth) 1 0.3841 230.75
+ log(Height) 1 4.7827 245.57
<none> 8.3172 255.48
Step: AIC=160.83
Volume ~ log(Girth)
Df Deviance AIC
+ log(Height) 1 0.1835 147.82
<none> 0.3841 160.83
- log(Girth) 1 8.3172 752.55
Step: AIC=139.9
Volume ~ log(Girth) + log(Height)
Df Deviance AIC
<none> 0.1835 139.90
- log(Height) 1 0.3841 169.11
- log(Girth) 1 4.7827 853.48En este caso todas las estrategias conducen al mismo modelo que la fundamentación teórica anticipaba:
coef(m.s)(Intercept) log(Girth) log(Height)
-6.691109 1.980412 1.132878