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 razón de verosimilitud para responder la pregunta, y en la siguiente sección se muestra como realizar la prueba.

2 Prueba razón de verosimilitudes

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}\]

Denotemos con \(k_0\) el número de \(\beta\)´s en M0 y con \(k_1\) el número de \(\beta\)´s en M1. Asímismo supongamos que el deviance de M0 se denota por \(D_0\) y que el deviance de M1 se denota para \(D_1\).

2.1 Prueba razón de verosimilitud para glm’s sin \(\phi\)

Recordemos que los glm Poisson y binomial tienen \(\phi=1\). En este caso el estadístico está dado por:

\[lrt = D_0 - D_1,\] y suponiendo que \(H_0\) es verdadera, el estadístico \(lrt\) tendrá una distribución \(\chi^2_{k_1-k_0}\) (Dunn and Smyth (2018)).

2.2 Prueba razón de verosimilitud para glm’s con \(\phi\)

Recordemos que los glm normal, gamma y inversa gaussian tienen parámetro de dispersión \(\phi\) que se puede estimar de 4 formas diferentes. En este caso el estadístico está dado por:

\[F_0 = \frac{D_0 - D_1}{(k_1-k_0) \, \tilde{\phi}_1},\] y suponiendo que \(H_0\) es verdadera, el estadístico \(F_0\) tendrá una distribución \(F_{k1-k0, n-k1}\) (Dunn and Smyth (2018)). El elemento \(\tilde{\phi}_1\) es una estimación para el parámetro de dispersión \(\phi\) del modelo M1 y \(n\) representa el número de observaciones usadas para ajustar el modelo.

3 Ejemplo con glm Binomial

En este ejemplo se desea modelar la variable \(Y\) usando un modelo \(Binomial(\mu, m=1)\). La variable \(Y=0\) si la cangreja NO tiene satelites pegados al caparazón y \(Y=1\) si la cangreja SI tiene satelites pegados al caparazón.

El objetivo es comparar dos modelos anidades, uno sin covariables (M0) y otro con solo la covariable width (M1), las hipótesis de interés se pueden escribir así:

\[\begin{align} H_0 &: \text{un modelo binomial sin covariables es apropiado para explicar } \mu \\ H_1 &: \text{es mejor un modelo binomial con width para explicar } \mu \end{align}\]

Lo primero que se debe hacer es ajustar los modelos de la siguiente manera.

url <- "http://www.stat.ufl.edu/~aa/cat/data/Crabs.dat"
Crabs <- read.table(url, header=TRUE)
fit0 <- glm(y ~ 1, family=binomial, data=Crabs)
fit1 <- glm(y ~ width, family=binomial, data=Crabs)

Para realizar la prueba razón de verosimilitud manual se puede usar el siguiente código.

lrt <- deviance(fit0) - deviance(fit1)
lrt
## [1] 31.30586
pchisq(q=lrt, df=2-1, lower.tail=FALSE) # Valor-p
## [1] 2.204134e-08

De la salida anterior vemos que el valor-P es 2.204134e-08, lo cual es menor que cualquier nivel de significancia, por lo tanto podemos decir que hay evidencias en favor de M1, es decir, es mejor un modelo con la covariable width.

Es posible obtener los mismos resultados anteriores usando la función anova pero especificando el tipo de test deseado así:

anova(fit0, fit1, test="Chisq", dispersion=1)
## Analysis of Deviance Table
## 
## Model 1: y ~ 1
## Model 2: y ~ width
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       172     225.76                          
## 2       171     194.45  1   31.306 2.204e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4 Ejemplo con glm inverse Gaussian

En este ejemplo se desea modelar la variable \(Y\) que representa el precio de venta de unas casas en función de algunas covariables.

El objetivo es comparar dos modelos anidades, uno con 2 covariables (M0) y otro con 4 covariables (M1), las hipótesis de interés se pueden escribir así:

\[\begin{align} H_0 &: \text{un modelo inverse gaussian con size y taxes es apropiado para explicar } \mu \\ H_1 &: \text{es mejor un modelo inverse gaussian con size, taxes, new y beds para explicar } \mu \end{align}\]

En este ejemplo se usará como estimación de \(\phi\) el estimador Mean Deviance.

Lo primero que se debe hacer es ajustar los modelos de la siguiente manera.

url <- 'http://users.stat.ufl.edu/~aa/glm/data/Houses.dat'
datos <- read.table(url, header=TRUE)
fit0 <- glm(price ~ size + taxes,
            data=datos, family=inverse.gaussian(link=log))
fit1 <- glm(price ~ size + new + beds + taxes,
            data=datos, family=inverse.gaussian(link=log))

A continuación se exploran las estimaciones de los parámetros de los modelo M0 y M1.

coef(fit0) # estimaciones de los betas
##  (Intercept)         size        taxes 
## 3.8225775433 0.0003504572 0.0002936544
fit0$deviance / fit0$df.residual  # estimacion de phi
## [1] 0.0009512286
coef(fit1) # estimaciones de los betas
##   (Intercept)          size           new          beds         taxes 
##  4.0019724531  0.0004222138  0.0435855974 -0.1016675494  0.0002983505
fit1$deviance / fit1$df.residual  # estimacion de phi
## [1] 0.0009483615

Para realizar la prueba razón de verosimilitud manual se puede usar el siguiente código.

F0 <- (deviance(fit0)-deviance(fit1)) / ((5-3) * 0.0009483615)
F0
## [1] 1.146627
pf(q=F0, df1=5-3, df2=100-5, lower.tail=FALSE) # Valor-p
## [1] 0.3220637

De la salida anterior vemos que el valor-P es 0.3220637, lo cual es mayor que cualquier nivel de significancia usual, por lo tanto podemos decir que NO hay evidencias en favor de M1, es decir, es mejor el modelo M0 que sólo tiene las covariables size y taxes.

Es posible obtener los mismo resultados anteriores usando la función anova pero especificando el tipo de test deseado y el valor de \(\phi\) así:

anova(fit0, fit1, test="F", dispersion=0.0009483615)
## Analysis of Deviance Table
## 
## Model 1: price ~ size + taxes
## Model 2: price ~ size + new + beds + taxes
##   Resid. Df Resid. Dev Df  Deviance      F Pr(>F)
## 1        97   0.092269                           
## 2        95   0.090094  2 0.0021748 1.1466 0.3221



References

Dunn, Peter K, and Gordon K Smyth. 2018. Generalized Linear Models with Examples in R. Springer.