1 Introducción

Al construir un modelo lineal generalizado para explicar la media (\(\mu\)) de una variable respuesta (\(Y\)), es posible que nos encontremos con la duda que se muestra en la siguiente figura.

¿Cómo se puede proceder en estos casos?

Una de las formas es utilizar la prueba Score para responder la pregunta y en la siguiente sección se muestra como realizar la prueba.

2 Prueba Score

Supongamos que se tienen los siguientes los modelos glm anidados M0 y M1 para explicar \(\mu\) en función varias covariables así:

\[\begin{align*} \text{M0: } g(\mu) &= \beta_0 + \beta_1 x_1 + \ldots + \beta_{A} x_{A} \\ \text{M1: } g(\mu) &= \beta_0 + \beta_1 x_1 + \ldots + \beta_{A} x_{A} + \beta_{A+1} x_{A+1} + \beta_{A+2} x_{A+2} + \ldots + \beta_{B} x_{B} \\ \end{align*}\]

El modelo M0 es un caso especial del modelo M1, porque M1 tiene las covariables de M0 más otras covariables, en este caso es claro que \(A < B\).

Las hipótesis para comparar los modelos M0 y M1 se pueden escribir de dos formas:

\[\begin{align} H_0 &: \text{el modelo M0 es apropiado para explicar } \mu \\ H_1 &: \text{es mejor el modelo M1 para explicar } \mu \end{align}\]

o como:

\[\begin{align} H_0 &: \beta_{A+1}=\beta_{A+2}=\ldots=\beta_B=0 \\ H_1 &: \text{al menos uno de los} \, \beta_j \neq0 \, \text{con} \, j=A+1,A+2,\ldots,B, \end{align}\]

El estadístico de la prueba Score está dado por:

\[ S = \bf{U}^\top(\hat{\boldsymbol{\theta}}_0) \, \bf{I}^{-1}(\hat{\boldsymbol{\theta}}_0) \, \bf{U}(\hat{\boldsymbol{\theta}}_0), \] donde \(\bf{U}\) es el vector Score evaluado en \(\hat{\boldsymbol{\theta}}_0\) y \(\bf{I}\) es la matriz de información observada evaluada en \(\hat{\boldsymbol{\theta}}_0\). Suponiendo que \(H_0\) es verdadera, el estadístico \(S\) tendrá una distribución \(\chi^2_{k_1-k_0}\) donde \(k_0\) es el número de \(\beta\)´s en M0 y con \(k_1\) el número de \(\beta\)´s en M1.

3 Ejemplo con glm Gamma

Este ejemplo está basado en la viñeta del paquete enrichwith de Kosmidis (2020). En este ejemplo se quiere modelar la media \(\mu\) de la variable respuesta time en función de las variables conc y lot asumiendo una distribución gamma para la variable respuesta time. Las dos opciones para modelar \(\mu\) se resumen a continuación.

\[\begin{align*} \text{M0: } \frac{1}{\mu} &= \beta_0 + \beta_1 \log(conc) \\ \text{M1: } \frac{1}{\mu} &= \beta_0 + \beta_1 \log(conc) + \beta_2 lot2 + \beta_3 \log(conc) \times lot2 \\ \end{align*}\]

Los datos para el ejemplo se muestran a continuación.

clotting <- data.frame(conc = c(5,10,15,20,30,40,60,80,100,
                                5,10,15,20,30,40,60,80,100),
                       time = c(118, 58, 42, 35, 27, 25, 21, 19, 18,
                                69, 35, 26, 21, 18, 16, 13, 12, 12),
                       lot = factor(c(rep(1, 9), rep(2, 9))))

Lo primero que se debe hacer es ajustar los modelos de interés.

fit0 <- glm(time ~ log(conc), family=Gamma, data=clotting)
fit1 <- glm(time ~ log(conc) * lot, family=Gamma, data=clotting)

Lo siguiente que se debe hacer para aplicar la prueba Score es determinar el vector \(\hat{\boldsymbol{\theta}}_0\), el subíndice 0 indica que el \(\hat{\boldsymbol{\theta}}\) es asumiendo \(H_0\) verdadera.

# Parametros estimados en fit0
coef(fit0) # betas
## (Intercept)   log(conc) 
## -0.01963451  0.01860894
MASS::gamma.dispersion(fit0) # Maximum likelihood estimate
## [1] 0.05604664
summary(fit0)$dispersion # Pearson Estimator
## [1] 0.06219859
fit0$deviance / fit0$df.residual # Mean Deviance Estimator
## [1] 0.06364126

Hay 3 estimadores disponibles para \(\phi\) y en este ejemplo vamos a utilizar el Maximum likelihood estimate, de esta forma se tiene que

\[ \hat{\boldsymbol{\theta}}_0 = (\hat{\beta}_0=-0.01963451, \hat{\beta}_1=0.01860894, \hat{\phi}=0.05604664)^\top \] No es necesario obtener \(\hat{\boldsymbol{\theta}}_1\) porque el estadístico de la prueba no lo requiere.

Lo siguiente que se debe hacer es obtener el vector score \(\bf{U}\) y la matriz de información observada \(\bf{I}\), esto se hace “enriqueciendo” el modelo fit1 con la función enrich del paquete enrichwith de Kosmidis (2020), el resultado es un objeto de la clase “enriched_glm,” a continuación el procedimiento.

# Para enriquecer el modelo glm mas general fit1
library("enrichwith")
enriched_fit1 <- enrich(fit1, with="auxiliary functions")

Para obtener el vector score \(\bf{U}\) se procede así:

scores_fit1 <- get_score_function(fit1)

Para explorar el vector score \(\bf{U}\) evaluado en \(\hat{\boldsymbol{\theta}}_1\) se hace lo siguiente.

scores_fit1()
##    (Intercept)      log(conc)           lot2 log(conc):lot2     dispersion 
##   1.227818e-05   2.253249e-05   6.113431e-07   1.548145e-06   1.647695e-09 
## attr(,"coefficients")
##    (Intercept)      log(conc)           lot2 log(conc):lot2 
##   -0.016554382    0.015343115   -0.007354088    0.008256099 
## attr(,"dispersion")
##  dispersion 
## 0.001632971

De la salida anterior se observa que el resultado es un vector con 5 elementos con valores muy cercanos a cero porque son las pendientes de la función de log verosimilitud en \(\hat{\boldsymbol{\theta}}_1\). Adicionalmente se obtienen como atributos los valores de \(\hat{\boldsymbol{\theta}}_1\).

Para obtener la matriz de información observada \(\bf{I}\) se procede así:

info_fit1 <- get_information_function(fit1)

Para explorar el vector score \(\bf{I}\) evaluado en \(\hat{\boldsymbol{\theta}}_1\) se hace lo siguiente.

info_fit1()
##                (Intercept) log(conc)     lot2 log(conc):lot2 dispersion
## (Intercept)       19327145  40768701  5060148       10887183          0
## log(conc)         40768701  98046714 10887183       26762283          0
## lot2               5060148  10887183  5060148       10887183          0
## log(conc):lot2    10887183  26762283 10887183       26762283          0
## dispersion               0         0        0              0    3376930
## attr(,"coefficients")
##    (Intercept)      log(conc)           lot2 log(conc):lot2 
##   -0.016554382    0.015343115   -0.007354088    0.008256099 
## attr(,"dispersion")
##  dispersion 
## 0.001632971

Luego se debe obtener \(\bf{U}(\hat{\boldsymbol{\theta}}_0)\) y \(\bf{I}(\hat{\boldsymbol{\theta}}_0)\) así:

# Creando los vectores con beta y disp bajo H0
beta_hat_0 <- c(-0.01963451, 0.01860894, 0, 0)
disp_hat_0 <- 0.05604664 # Maximum likelihood estimate

# Creando el vector y matriz para la prueba
scores <- scores_fit1(beta_hat_0, disp_hat_0)
info <- info_fit1(beta_hat_0, disp_hat_0)

Por último calculamos el estadístico usando la expresión \(S = \bf{U}^\top(\hat{\boldsymbol{\theta}}_0) \, \bf{I}^{-1}(\hat{\boldsymbol{\theta}}_0) \, \bf{U}(\hat{\boldsymbol{\theta}}_0)\), computacionalmente los cálculos se harían así:

score_statistic <- drop(scores%*%solve(info)%*%scores)
score_statistic
## [1] 17.15989
pchisq(score_statistic, 2, lower.tail=FALSE) # valor-P
## [1] 0.0001878352

De la salida anterior se obtiene que el valor-P es 0.0001878352, más pequeño que cualquier nivel de significancia \(\alpha\), por lo tanto, hay evidencias para afirmar que es mejor el modelo M1, es decir, es mejor incluir en el modelo las variables \(\log(conc)\), \(lot\) y la interacción entre ellas.



References

Kosmidis, Ioannis. 2020. enrichwith: Methods to Enrich List-Like r Objects with Extra Components. https://github.com/ikosmidis/enrichwith.