Ejemplo: Número de muertes por cancer de pulmón -Regresión de Poisson y Binomial Negativa-

Especialización en Estadística Aplicada

Roberto Trespalacios

Vamos a usar el conjunto de datos correspondiente al número de muertes por cáncer de pulmón (cander.dat) el cual puede ser descargado de la siguiente dirección: http://stat.ufl.edu/~aa/glm/data/ del autor [@agresti2015foundations].

La base de datos contiene 3 filas y 5 columnas(variables):

• Conteo (count) del número de pacientes que murieron clasificados por el estado de la enfermedad (stage)
• Tiempo de seguimiento (time).
• Histología (histology), es el tipo de cancer que padece el paciente.
  
• Etapa (stage) en la que se encuentra la enfermedad del paciente.
• Risktime (risktime) corresponde al número de pacientes en riesgo en el período de seguimiento.

• Note que la devianza residual es \( G^2 = 79.989(gl = 52) \). Entonces \( G^2 /gl = 79.989/ 52=1.5382 \)
• \( G^2 /gl \) es una cantidad no muy lejana de 1, lo cual indica un ajuste razonable del modelo.
• En el caso de conteos con distribución Poisson este resultado es señal de sobredispersión, es decir, varianza mayor que la media.
• La devianza nula corresponde a la razón de verosimilitud entre el modelo saturado y el modelo nulo (con intercepto únicamente). Algunos autores usan la cantidad

\[ R^2_{pseudo}=\frac{Null\ deviance-Residual\ deviance}{Null\ deviance} \]

como una medida de ajuste similar al coeficiente de determinación en modelos lineales.

• Esta cantidad se conoce como el pseudo \( R^2 \). En este caso \( R^2_{pseudo} = 0.81 \), lo cual indicaría que las variables explican aproximadamente el 82\% de la variabilidad del número de muertes.

• En términos de interpretación de los coeficientes podemos decir, por ejemplo, que el número promedio de muertes en la etapa 3 es aproximadamente \( \exp(1.3286) = 3.776 \) veces el número promedio de muertes de la etapa 1, después de ajustar por tiempo de seguimiento y histología (intervalo del 95\% de confianza (2.848, 5.087)).

• Según la salida, todas las variables son significativas para explicar el número promedio de muertes por cáncer. Este resultado se debe revisar ya que estamos ignorando el número de pacientes en riesgo.

• Si queremos determinar si hay diferencia significativa entre el número promedio de muertes entre la etapa 2 y la etapa 1 (referencia) es suficiente con llevar a cabo una prueba de hipótesis para el coeficiente de la etapa 2.

• Como el \( p-valor \) es menor que 0.05 (\( p = 0.121084 \)) entonces no rechazamos \( H_0: \beta_2^S = 0 \).

• Es decir, podemos concluir que el número promedio de muertes en la etapa 2 y 1 no son estadísticamente diferentes.

• Note que la estimación puntual de esta razón es igual a \( \exp(0.2703) = 1.31 \) con un intervalo de 95\% de confianza (0.933, 1.851).

• El hecho que el intervalo incluya el 1 respalda el resultado de la prueba de hipótesis (razón de promedio de muertes se puede considerar 1:1).

• Si el investigador está interesado en estimar la razón del promedio de muertes entre la etapa 2 y 3 entonces se puede recurrir a dos estrategias.

• La primera consiste en ajustar nuevamente el modelo pero considerando el tercer nivel de la variable stage como referencia.

• La seguna estrategia es escribir la razón del promedio de muertes entre la etapa 2 y 3 en términos del modelo original.

• La segunda estrategia se usa la función estimable de la librería gmodels con el modelo original que tiene el primer nivel de la etapa como referencia. Antes de usar el programa veamos como escribimos la razón promedio de muertes entre la etapa 2 y 3 en términos del modelo. Se deja para el lector.

# Usando la combinación lineal de coeficientes
fit.cancer = glm( count ~ factor(histology) + factor(stage) + factor(time),
family = poisson(link = log), data = cancer)

library(gmodels)
contrast = estimable(fit.cancer, c("factor(stage)2"=1, "factor(stage)3"=-1),
conf.int=.95)
contrast

                         Estimate Std. Error X^2 value DF Pr(>|X^2|)
(0 0 0 1 -1 0 0 0 0 0 0)    -1.06      0.133      63.2  1   1.89e-15
                         Lower.CI Upper.CI
(0 0 0 1 -1 0 0 0 0 0 0)    -1.32   -0.794

Ejemplo 1. Los administradores escolares estudian el comportamiento de ausencia de los estudiantes de tercer año de secundaria en dos escuelas. Los predictores de la cantidad de días de ausencia incluyen el tipo de programa en el que se inscribe el alumno y una prueba estandarizada en matemáticas.

Ejemplo 2. Un investigador relacionado con la salud está estudiando el número de visitas hospitalarias realizadas en los últimos 12 meses por personas de la tercera edad en una comunidad según las características de las personas y los tipos de planes de salud según cada uno de ellos.

• R primero muestra la llamada y los residuos de la desviación.
• A continuación, vemos los coeficientes de regresión para cada una de las variables, junto con los errores estándar, los puntajes z y los valores p.
• La variable matemática tiene un coeficiente de -0.006, que es estadísticamente significativo.
• Esto significa que por cada aumento de una unidad en matemáticas, el recuento de registros esperado del número de días de ausencia disminuye en 0.006.
• La variable indicadora mostrada como progAcademic es la diferencia esperada en el recuento de registros entre el grupo 2 y el grupo de referencia (prog = 1).
• El recuento de registros esperado para el nivel 2 de prog es 0.44 menor que el recuento de registros esperado para el nivel 1.

• La variable indicadora para progVocational es la diferencia esperada en el recuento de registros entre el grupo 3 y el grupo de referencia(prog = 1).
• El recuento de registros esperado para el nivel 3 de prog es 1.28 más bajo que el recuento de registros esperado para el nivel 1.
• Para determinar si el progreso en sí mismo, en general, es estadísticamente significativo, podemos comparar un modelo con y sin prog.
• La razón por la cual es importante ajustar modelos separados es que, a menos que lo hagamos, el parámetro de sobredispersión se mantiene constante.

Comparemos los modelos mod.nb y mod.nb_sin_prog con la proeba LRT (anova)

anova(mod.nb, mod.nb_sin_prog)

Likelihood ratio tests of Negative Binomial Models

Response: daysabs
        Model theta Resid. df    2 x log-lik.   Test    df LR stat.
1        math 0.856       312           -1776                      
2 math + prog 1.033       310           -1731 1 vs 2     2       45
         Pr(Chi)
1               
2 0.000000000165

• La prueba de chi-cuadrado con dos grados de libertad indica que el progreso es un predictor estadísticamente significativo de daysabs.

• La desviación nula se calcula a partir de un modelo de solo interceptación con 313 grados de libertad. Luego vemos la desviación residual, la desviación del modelo completo. También se nos muestra el AIC y la probabilidad de \( 2*log \).

• Como mencionamos anteriormente, los modelos binomiales negativos suponen que las medias condicionales no son iguales a las varianzas condicionales.
• Esta desigualdad se captura al estimar un parámetro de dispersión (no mostrado en la salida) que se mantiene constante en un modelo de Poisson.
• Por lo tanto, el modelo de Poisson está realmente anidado en el modelo binomial negativo. Luego podemos usar una prueba de razón de verosimilitud para comparar estos dos y probar la suposición de este modelo.
• Para hacer esto, ejecutaremos nuestro modelo con función enlace (link) Poisson.

La forma de la ecuación modelo para la regresión binomial negativa es la misma que para la regresión de Poisson. El registro del resultado se predice con una combinación lineal de los predictores:

\[ ln(\widehat{daysabs_i}) = \hat{\beta}_0 + \hat{\beta}_1(prog_i = 2) + \hat{\beta}(prog_i = 3) + \hat{\beta}math_i \]

\[ \begin{align*} \widehat{daysabs_i} &= e^{\hat{\beta}_0 + \hat{\beta}_1(prog_i = 2) + \hat{\beta}(prog_i = 3) + \hat{\beta}math_i} \\ &= e^{\hat{\beta}_0}e^{\hat{\beta}_1(prog_i = 2)}e^{b_2(prog_i = 3)}e^{\hat{\beta}math_i} \end{align*} \]