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

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: \beta_j = \beta_j^{(0)},\] donde \(\beta_j^{(0)}\) es el valor de referencia del parámetro \(\beta_j\).

En esta situación pueden darse una de las siguientes hipótesis alterna

\[ H_1: \beta_j < \beta_j^{(0)}, \quad H_1: \beta_j \neq \beta_j^{(0)}, \quad H_1: \beta_j > \beta_j^{(0)}. \]

2.1 Prueba de Wald con \(\phi\) conocido

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

\[Z = \frac{\hat{\beta}_j - \beta_j^{(0)}}{se(\hat{\beta}_j)},\] y suponiendo que \(H_0\) es verdadera, el estadístico \(Z\) tendrá una distribución \(N(0, 1)\) (Dunn and Smyth (2018), sección 7.2.1).

2.2 Prueba de Wald 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:

\[T = \frac{\hat{\beta}_j - \beta_j^{(0)}}{se(\hat{\beta}_j)},\] y suponiendo que \(H_0\) es verdadera, el estadístico \(T\) tendrá una distribución \(t_{n-k-1}\) (Dunn and Smyth (2018), sección 7.6.1). 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} = 0.1 \\ H_1 &: \beta_{weight} > 0.1 \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)

La tabla de resumen del modelo ajustado se obtiene de la siguiente manera.

summary(fit)
## 
## Call:
## glm(formula = y ~ weight + color + width + spine, family = poisson(link = log), 
##     data = Crabs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0290  -1.8630  -0.5988   0.9331   4.9446  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.36180    0.96655  -0.374  0.70817   
## weight       0.49647    0.16626   2.986  0.00283 **
## color2      -0.26485    0.16811  -1.575  0.11515   
## color3      -0.51371    0.19536  -2.629  0.00855 **
## color4      -0.53086    0.22692  -2.339  0.01931 * 
## width        0.01675    0.04892   0.342  0.73207   
## spine2      -0.15037    0.21358  -0.704  0.48139   
## spine3       0.08728    0.11993   0.728  0.46674   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 632.79  on 172  degrees of freedom
## Residual deviance: 549.59  on 165  degrees of freedom
## AIC: 920.88
## 
## Number of Fisher Scoring iterations: 6

Para realizar la prueba Wald se puede usar el siguiente código.

z <- (0.49647 - 0.1) / 0.16626
z
## [1] 2.384639
pnorm(q=z, lower.tail=FALSE) # valor-P
## [1] 0.008547955

De la salida anterior vemos que el valor-P es 0.008547955, lo cual es menor que cualquier nivel de significancia, por lo tanto podemos decir que hay evidencias para afirmar que \(\beta_{weight} > 0.1\).

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_{size} = -0.000001 \\ H_1 &: \beta_{size} < -0.000001 \end{align}\]

Suponga que vamos a usar como estimador de \(\phi\) el Mean Deviance 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))

La tabla de resumen del modelo ajustado se obtiene de la siguiente manera.

summary(fit)
## 
## Call:
## glm(formula = price ~ size + new + beds, family = Gamma(link = inverse), 
##     data = datos)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.28183  -0.25692   0.02211   0.14277   0.87938  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.183e-02  9.237e-04  12.804  < 2e-16 ***
## size        -2.179e-06  2.636e-07  -8.266 7.74e-13 ***
## new         -1.229e-03  4.705e-04  -2.612   0.0105 *  
## beds        -2.540e-04  3.455e-04  -0.735   0.4640    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1286359)
## 
##     Null deviance: 31.940  on 99  degrees of freedom
## Residual deviance: 13.581  on 96  degrees of freedom
## AIC: 1073.6
## 
## Number of Fisher Scoring iterations: 6

Para obtener el Mean Deviance estimator de \(\phi\) hacemos lo siguiente:

deviance(fit) / df.residual(fit)
## [1] 0.1414655

Con el estimador \(\hat{\phi}\) anterior debemos actualizar la tabla de resumen para obtener los errores estándar apropiados para la prueba.

printCoefmat(coef(summary(fit, dispersion=0.1414655)))
##                Estimate  Std. Error z value  Pr(>|z|)    
## (Intercept)  1.1828e-02  9.6869e-04 12.2099 < 2.2e-16 ***
## size        -2.1790e-06  2.7644e-07 -7.8824 3.212e-15 ***
## new         -1.2289e-03  4.9343e-04 -2.4905   0.01276 *  
## beds        -2.5404e-04  3.6237e-04 -0.7011   0.48326    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Para realizar la prueba Wald se puede usar el siguiente código.

T <- (-2.1790e-06 - (-0.000001)) / 2.7644e-07
T
## [1] -4.26494
pt(q=T, df=100-4, lower.tail=TRUE) # valor-P
## [1] 2.341904e-05

De la salida anterior vemos que el valor-P es 2.341904e-05, por lo cual se concluye que SI hay evidencias para rechazar \(H_0\) y concluir que \(\beta_{size} < -0.000001\).



References

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