Prueba de Wald multivariada en GLM’s
Freddy Hernández
22 julio, 2020
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 de Wald para responder la pregunta, y en la siguiente sección se muestra como realizar la prueba.
2 Prueba de Wald multivariada
Supongamos que se ajustó un modelo glm para explicar \(\mu\) en función de \(k\) covariables así:
\[ g(\mu) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_{k} x_{k} \] La hipótesis nula en la prueba de Wald es
\[H_0: \underset{q \times (k+1)}{\bf{L}} \underset{(k+1) \times 1}{\boldsymbol{\beta}} = \underset{q \times 1}{\bf{C}},\] donde \(\bf{L}\) es una matriz de ceros y unos para extraer los elementos de interés del vector \(\boldsymbol{\beta}\), la matriz \(\bf{C}\) contiene los valores de referencia de cada uno de los \(\beta\)’s.
La hipótesis alterna en esta prueba es
\[H_1: \underset{q \times (k+1)}{\bf{L}} \underset{(k+1) \times 1}{\boldsymbol{\beta}} \neq \underset{q \times 1}{\bf{C}},\]
2.1 Prueba de Wald multivariada con \(\phi\) conocido
Recordemos que los glm Poisson y binomial tienen \(\phi=1\). En este caso el estadístico está dado por:
\[Z_0^2 = \left( \bf{L} \hat{\boldsymbol{\beta}} - \bf{C} \right)^{\top} \left( \bf{L} \widehat{Var}(\hat{\boldsymbol{\beta}}) \bf{L}^{\top} \right)^{-1} \left( \bf{L} \hat{\boldsymbol{\beta}} - \bf{C} \right),\] y suponiendo que \(H_0\) es verdadera, el estadístico \(Z_0^2\) tendrá una distribución \(\chi^2_{q}\) (Fox (2015), sección 15.3.3).
2.2 Prueba de Wald multivariada con \(\phi\) desconocido
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{\left( \bf{L} \hat{\boldsymbol{\beta}} - \bf{C} \right)^{\top} \left( \bf{L} \widehat{Var}(\hat{\boldsymbol{\beta}}) \bf{L}^{\top} \right)^{-1} \left( \bf{L} \hat{\boldsymbol{\beta}} - \bf{C} \right)}{q},\] y suponiendo que \(H_0\) es verdadera, el estadístico \(F_0\) tendrá una distribución \(F_{q, n-k-1}\) (Fox (2015), sección 15.3.3). El valor de \(n\) corresponde al número de observaciones usadas en el modelo.
3 Ejemplo con glm Poisson
En este ejemplo se desea modelar la variable \(Y\) que representa el número de cangrejos satelites (machos) pegados al caparazón de las hembras en funcion de las covariables weight, color, width y spine.
Se va a ajustar un modelo Poisson(log) y se desea estudiar las siguientes hipótesis.
\[\begin{align} H_0 &: (\beta_{weight}, \beta_{width})^\top = (0, 0)^\top \\ H_1 &: (\beta_{weight}, \beta_{width})^\top \neq (0, 0)^\top \end{align}\]
Lo primero que se debe hacer es ajustar el modelo de la siguiente manera.
url <- "http://users.stat.ufl.edu/~aa/glm/data/Crabs.dat"
Crabs <- read.table(url, header=TRUE)
Crabs$color <- as.factor(Crabs$color)
Crabs$spine <- as.factor(Crabs$spine)
fit <- glm(y ~ weight + color + width + spine,
family=poisson(link=log), data=Crabs)El vector con los \(\beta\)’s estimados se obtiene así:
beta_hat <- matrix(coef(fit), ncol=1,
dimnames=list(names(coef(fit)), NULL))
beta_hat## [,1]
## (Intercept) -0.36180032
## weight 0.49647004
## color2 -0.26485119
## color3 -0.51370507
## color4 -0.53086011
## width 0.01674870
## spine2 -0.15037184
## spine3 0.08728258Como la hipótesis nula involucra sólo \(\beta_{weight}\) y \(\beta_{width}\) entonces la matriz \(\bf{L}\) debe ser la siguiente.
L <- matrix(c(0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0), ncol=8, byrow=TRUE)
L## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 1 0 0 0 0 0 0
## [2,] 0 0 0 0 0 1 0 0Para verificar si la matriz \(\bf{L}\) anterior es la correcta se debe multiplicar por \(\hat{\boldsymbol{\beta}}\) para obtener solo los \(\beta\)’s involucrados en la hipótesis nula.
L %*% beta_hat # subvector con los betas en H0## [,1]
## [1,] 0.4964700
## [2,] 0.0167487El vector anterior de dimensión \(2 \times 1\) si contiene la información de los \(\beta\)’s involucrados en la hipótesis nula.
La matriz \(\bf{C}\) con los lados derechos de la hipótesis nula es la siguiente.
C <- matrix(c(0, 0), ncol=1, byrow=TRUE)
C## [,1]
## [1,] 0
## [2,] 0Para realizar la prueba Wald se puede usar el siguiente código.
aux <- L %*% beta_hat - C
z2 <- t(aux) %*% solve(L %*% vcov(fit) %*% t(L)) %*% aux
z2 <- as.numeric(z2)
z2## [1] 55.86176pchisq(q=z2, df=2, lower.tail=FALSE) # valor-P## [1] 7.409242e-13De la salida anterior vemos que el valor-P es 7.409242e-13, lo cual es menor que cualquier nivel de significancia, por lo tanto podemos decir que hay evidencias para afirmar que el vector \((\beta_{weight}, \beta_{width})^\top\) no es igual al vector nulo, en otras palabras, que al menos uno de los \(\beta\)’s involucrados es diferente de cero.
4 Ejemplo con glm Gamma
En este ejemplo se desea modelar la variable \(Y\) que representa el precio de venta de unas casas en funcion de las covariables size, new y beds.
Se va a ajustar un modelo Gamma(inv) y se desea estudiar las siguientes hipótesis.
\[\begin{align} H_0 &: (\beta_{new}, \beta_{beds})^\top = (0, 0)^\top \\ H_1 &: (\beta_{new}, \beta_{beds})^\top \neq (0, 0)^\top \end{align}\]
Suponga que vamos a usar como estimador de \(\phi\) el Maximum Likelihood estimator.
Lo primero que se debe hacer es ajustar el modelo de la siguiente manera.
url <- 'http://users.stat.ufl.edu/~aa/glm/data/Houses.dat'
datos <- read.table(url, header=TRUE)
fit <- glm(price ~ size + new + beds,
data=datos, family=Gamma(link=inverse))El vector con los \(\beta\)’s estimados se obtiene así:
beta_hat <- matrix(coef(fit), ncol=1,
dimnames=list(names(coef(fit)), NULL))
beta_hat## [,1]
## (Intercept) 1.182760e-02
## size -2.179041e-06
## new -1.228912e-03
## beds -2.540445e-04Como la hipótesis nula involucra sólo \(\beta_{new}\) y \(\beta_{beds}\) entonces la matriz \(\bf{L}\) debe ser la siguiente.
L <- matrix(c(0, 0, 1, 0,
0, 0, 0, 1), ncol=4, byrow=TRUE)
L## [,1] [,2] [,3] [,4]
## [1,] 0 0 1 0
## [2,] 0 0 0 1Para verificar si la matriz \(\bf{L}\) anterior es la correcta se debe multiplicar por \(\hat{\boldsymbol{\beta}}\) para obtener solo los \(\beta\)’s involucrados en la hipótesis nula.
L %*% beta_hat # subvector con los betas en H0## [,1]
## [1,] -0.0012289122
## [2,] -0.0002540445El vector anterior de dimensión \(2 \times 1\) si contiene la información de los \(\beta\)’s involucrados en la hipótesis nula.
La matriz \(\bf{C}\) con los lados derechos de la hipótesis nula es la siguiente.
C <- matrix(c(0, 0), ncol=1, byrow=TRUE)
C## [,1]
## [1,] 0
## [2,] 0Para obtener el Maximum Likelihood estimator de \(\phi\) se usa lo siguiente.
MASS::gamma.dispersion(fit)## [1] 0.1328696La matriz \(\widehat{Var}(\hat{\boldsymbol{\beta}})\) se debe actualizar usando el \(\hat{\phi}\) obtenido antes.
var_cov_hat <- vcov(summary(fit, dispersion=0.1328696))
var_cov_hat## (Intercept) size new beds
## (Intercept) 8.813446e-07 5.156914e-11 -4.086979e-08 -2.761805e-07
## size 5.156914e-11 7.177816e-14 -7.109259e-11 -6.263171e-11
## new -4.086979e-08 -7.109259e-11 2.286833e-07 3.935481e-08
## beds -2.761805e-07 -6.263171e-11 3.935481e-08 1.233321e-07Para realizar la prueba Wald se puede usar el siguiente código.
aux <- L %*% beta_hat - C
F0 <- t(aux) %*% solve(L %*% var_cov_hat %*% t(L)) %*% aux
F0 <- as.numeric(F0) / 2
F0## [1] 3.309771pf(q=F0, df1=2, df2=100-4, lower.tail=FALSE) # valor-P## [1] 0.04073577De la salida anterior vemos que el valor-P es 0.04073577, un poco menor que el nivel de significancia usual de 5%, por lo tanto podemos decir que hay evidencias para afirmar que el vector \((\beta_{new}, \beta_{beds})^\top\) no es igual al vector nulo, en otras palabras, que al menos uno de los \(\beta\)’s involucrados es diferente de cero.
Nota: Para conocer otras publicaciones relacionadas con glm visite este enlace.
References
Fox, John. 2015. Applied Regression Analysis and Generalized Linear Models. Sage Publications.