yes

Ejemplo: Número de muertes por cancer de pulmón (Agresti 2015)

Vamos a usar el conjunto de datos correspondiente al número de muertes por cáncer de pulmón (cander.dat). La base de datos contiene el número de pacientes que murieron clasificados por el estado de la enfermedad (stage), tiempo de seguimiento (time) e histología (histology). La base de datos contiene 63 filas y 5 columnas (variables). La variable risktime corresponde al número de pacientes en riesgo en el período de seguimiento.

# Lee base de datos del sitio web
cancer <- read.csv("http://stat.ufl.edu/~aa/glm/data/Cancer.dat", sep="")
str(cancer)

## 'data.frame':    63 obs. of  5 variables:
##  $ time     : int  1 1 1 2 2 2 3 3 3 4 ...
##  $ histology: int  1 2 3 1 2 3 1 2 3 1 ...
##  $ stage    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ count    : int  9 5 1 2 2 1 9 3 1 10 ...
##  $ risktime : int  157 77 21 139 68 17 126 63 14 102 ...

head(cancer)

Veamos la distribución de los conteos por estadío (stage) y la media y desviación estándar de los conteos en cada estadío.

# Histogram of counts by stage
library(ggplot2)
ggplot(cancer, aes(count, fill = stage)) +
  geom_histogram(binwidth=1) +
  facet_grid(stage ~ ., margins=TRUE, scales="free")

# Mean and SD by stage
with(cancer, tapply(count, stage, function(x) {
  sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))
}))

##                        1                        2                        3 
##   "M (SD) = 2.76 (2.96)"   "M (SD) = 3.62 (3.06)" "M (SD) = 10.43 (10.77)"

Por ahora vamos a modelar el número de pacientes que murieron sin tener en cuenta el número de pacientes en riesgo. Asumamos que \(Y_{i} \sim\) Poisson \(\left(\mu_{i}\right)\) donde \[ \log \left(\mu_{i}\right)=\beta_{0}+\beta_{j}^{H}+\beta_{k}^{S}+\beta_{l}^{T} \] Note que todas las variables explicativas son factores. Por ejemplo, para histología hay tres coeficientes \(\beta_{1}^{H}\) \(\beta_{2}^{H} y \beta_{3}^{H}\) pero como es costumbre en \(\mathrm{R}\) el primer coeficiente es igual a 0 (primer nivel de referencia)

# MLG para conteo de muertes de cáncer
fit.cancer = glm( count ~ factor(histology) + factor(stage) + factor(time), 
                  family = poisson(link = log), data = cancer)
summary(fit.cancer)

## 
## Call:
## glm(formula = count ~ factor(histology) + factor(stage) + factor(time), 
##     family = poisson(link = log), data = cancer)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5404  -0.9275  -0.2889   0.7252   2.4529  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          2.3163     0.1594  14.535  < 2e-16 ***
## factor(histology)2  -0.5108     0.1217  -4.197 2.71e-05 ***
## factor(histology)3  -1.0186     0.1447  -7.039 1.94e-12 ***
## factor(stage)2       0.2703     0.1744   1.550 0.121084    
## factor(stage)3       1.3286     0.1477   8.997  < 2e-16 ***
## factor(time)2       -0.5191     0.1488  -3.488 0.000487 ***
## factor(time)3       -0.7885     0.1626  -4.848 1.24e-06 ***
## factor(time)4       -0.9040     0.1693  -5.339 9.37e-08 ***
## factor(time)5       -1.9626     0.2590  -7.577 3.53e-14 ***
## factor(time)6       -1.8001     0.2414  -7.457 8.83e-14 ***
## factor(time)7       -1.8514     0.2468  -7.502 6.27e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 419.708  on 62  degrees of freedom
## Residual deviance:  79.989  on 52  degrees of freedom
## AIC: 287.81
## 
## Number of Fisher Scoring iterations: 5

# MLG para conteo de muertes de cáncer
library(car)

## Loading required package: carData

Anova(fit.cancer, type=3)

exp(confint(fit.cancer)) # intervalo confianza razón

## Waiting for profiling to be done...

##                         2.5 %     97.5 %
## (Intercept)        7.34476142 13.7260360
## factor(histology)2 0.47139765  0.7600624
## factor(histology)3 0.27008568  0.4767558
## factor(stage)2     0.93271862  1.8505531
## factor(stage)3     2.84843861  5.0874067
## factor(time)2      0.44264355  0.7941470
## factor(time)3      0.32815278  0.6216442
## factor(time)4      0.28809458  0.5604103
## factor(time)5      0.08158672  0.2266697
## factor(time)6      0.09999561  0.2589176
## factor(time)7      0.09381927  0.2482086

La salida de este modelo merece varios comentarios:

Note que la devianza residual es \(G^{2}=79.989(g l=52)\). El cociente de estas dos cantidades es igual a 1.5382, una cantidad no muy lejana de 1 , lo cual indica un ajuste razonable del modelo (recuerde que valores lejanos de 1 indican que el modelo no ajusta bien a los datos. 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_{\text {pseudo }}^{2}=\frac{\text { Null deviance-Residual deviance }}{\text { 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_{\text {pseudo }}^{2}=0.81\), lo cual indicaria 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 \mathrm{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 \(.05(p=.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 (.2703)=1.31\) con un intervalo de \(95 \%\) de confianza \([0.9327,1.8505]\). 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 esta 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.

# Modelo con stage=3 como nivel de referencia
f.stage = factor(cancer$stage)
fit.cancer2 = glm( count ~ factor(histology) + f.stage + factor(time), 
                   contrasts=list(f.stage="contr.SAS"), 
                   family = poisson(link = log), data = cancer)
coef(fit.cancer2)

##        (Intercept) factor(histology)2 factor(histology)3           f.stage1 
##          3.6448830         -0.5108256         -1.0185696         -1.3286287 
##           f.stage2      factor(time)2      factor(time)3      factor(time)4 
##         -1.0583384         -0.5191244         -0.7884574         -0.9039702 
##      factor(time)5      factor(time)6      factor(time)7 
##         -1.9625772         -1.8000583         -1.8513516

coef(fit.cancer2)

##        (Intercept) factor(histology)2 factor(histology)3           f.stage1 
##          3.6448830         -0.5108256         -1.0185696         -1.3286287 
##           f.stage2      factor(time)2      factor(time)3      factor(time)4 
##         -1.0583384         -0.5191244         -0.7884574         -0.9039702 
##      factor(time)5      factor(time)6      factor(time)7 
##         -1.9625772         -1.8000583         -1.8513516

exp( confint(fit.cancer2) )

## Waiting for profiling to be done...

##                          2.5 %     97.5 %
## (Intercept)        30.52331596 47.5029762
## factor(histology)2  0.47139765  0.7600624
## factor(histology)3  0.27008568  0.4767558
## f.stage1            0.19656380  0.3510695
## f.stage2            0.26575247  0.4481972
## factor(time)2       0.44264355  0.7941470
## factor(time)3       0.32815278  0.6216442
## factor(time)4       0.28809458  0.5604103
## factor(time)5       0.08158672  0.2266697
## factor(time)6       0.09999561  0.2589176
## factor(time)7       0.09381927  0.2482086

Usando la primera estrategia tenemos que la razón de muertes promedio de la etapa 2 es \(\exp (-1.0583)=0.347\) veces el promedio de muertes de la etapa \(3 .\) Es decir, por cada 3 muertes en la etapa 3 ocurren en promedio una muerte en la etapa \(2(\mathrm{IC} 95 \%[0.2657,0.4482])\) Para la segunda estrategia usaremos 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. Con un poco de álgebra podemos ver que esta razón se puede escribir como \(\exp \left(\beta_{2}^{S}-\beta_{3}^{S}\right)\). Este resultado se obtiene diviendo la razón entre la etapa \(2 y 1\) sobre la razón entre la etapa 3 y 1 (Ejercicio).

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

exp(contrast$Estimate) # estimación puntual

## [1] 0.347032

# Intervalo de confianza - dos formas de calcularlo
exp(contrast$Lower.CI) # límite inferior 95% CI

## [1] 0.2664767

exp(contrast$Estimate-1.96*contrast$Std.) # límite inferior 95% CI

## [1] 0.2673285

exp(contrast$Upper.CI) # límite superior 95% CI

## [1] 0.4519388

exp(contrast$Estimate+1.96*contrast$Std.) # límite superior 95% CI

## [1] 0.4504988

Finalmente, vamos a graficar los conteos observados versus las predicciones y los residuales para el modelo ajustado (figura ).

# Observed vs. predicted counts
par( mfrow=c(1,2) )
plot(fitted(fit.cancer), resid(fit.cancer))
plot(cancer$count, fitted(fit.cancer))

Modelos para Tasas

En el ejemplo de las muertes por cáncer el número de muertes se debe dividir por el número de sujetos en riesgo para obtener la tasa de muertes por cáncer. En este caso, en lugar de modelar los conteos definimos una tasa (conteos/total) en la respuesta. Es decir. \[ \log \left(\mu_{i} / t_{i}\right)=\mathbf{x}_{i}^{t} \beta_{i} \] donde \(t_{i}\) es el número de sujetos en riesgo. Note que si \(t_{i}\) es conocida y usando propiedades de la función logaritimica la expresión anterior se puede escribir como \[ \log \left(\mu_{i}\right)=\mathbf{x}_{i}^{t} \beta+\log \left(t_{i}\right) \] La cantidad \(\log \left(t_{i}\right)\) se le conoce como el of fset . De esta manera, el modelo anterior se ajusta fijando en uno el coeficiente para esta variable. Note que en este caso \[ \mu_{i}=t_{i} e^{\mathrm{x}^{\mathrm{\beta}}, \boldsymbol{\beta}} \]

Vamos a modelar los datos de cáncer usando el offset igual al logaritmo del número de pacientes en riesgo risktime Según la salida, únicamente la variable stage es significativa para explicar la tasa promedio de muertes. Este resultado contrasta con el modelo de la sección anterior el cual no tuvo en cuenta el total de pacientes en riesgo. Para este modelo la devianza residual es \(G^{2}=43.923(g l=52) .\) El cociente de estas dos cantidades es igual a \(0.84\), indicando que el modelo tiene un buen ajuste (cercano a 1) En términos de interpretación para el coeficiente de stage 3 podemos decir la tasa promedio de muertes en la etapa 3 es aproximadamente exp \((1.324)=3.76\) veces la tasa promedio de muertes de la etapa 1 . después de ajustar por tiempo de seguimiento y histología (intervalo del \(95 \%\) de confianza \([2.810,5.106]\}\) Esta cantidad \(3.759\) se le conoce como la razón de tasas (Rate Ratio)

# MLG para tasa de muertes de cáncer
logrisk = log(cancer$risktime)
fit.cancer1 = glm( count ~ factor(histology) + factor(stage) + factor(time), 
                  family = poisson(link = log), offset=logrisk, data = cancer)
summary(fit.cancer1)

## 
## Call:
## glm(formula = count ~ factor(histology) + factor(stage) + factor(time), 
##     family = poisson(link = log), data = cancer, offset = logrisk)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.00333  -0.74769  -0.03194   0.46468   1.70832  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -3.00928    0.16651 -18.073  < 2e-16 ***
## factor(histology)2  0.16244    0.12195   1.332  0.18285    
## factor(histology)3  0.10754    0.14745   0.729  0.46580    
## factor(stage)2      0.47001    0.17444   2.694  0.00705 ** 
## factor(stage)3      1.32431    0.15205   8.709  < 2e-16 ***
## factor(time)2      -0.12745    0.14908  -0.855  0.39259    
## factor(time)3      -0.07973    0.16352  -0.488  0.62585    
## factor(time)4       0.11892    0.17107   0.695  0.48694    
## factor(time)5      -0.66511    0.26061  -2.552  0.01071 *  
## factor(time)6      -0.35015    0.24348  -1.438  0.15040    
## factor(time)7      -0.17518    0.24985  -0.701  0.48321    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 175.718  on 62  degrees of freedom
## Residual deviance:  43.923  on 52  degrees of freedom
## AIC: 251.74
## 
## Number of Fisher Scoring iterations: 5

Anova(fit.cancer1, type=3)

exp(confint(fit.cancer1)) # exponencial de coeficientes

## Waiting for profiling to be done...

##                         2.5 %     97.5 %
## (Intercept)        0.03524982 0.06774327
## factor(histology)2 0.92382408 1.49090908
## factor(histology)3 0.82864241 1.47851437
## factor(stage)2     1.13871258 2.26002502
## factor(stage)3     2.81058090 5.10653503
## factor(time)2      0.65456940 1.17546240
## factor(time)3      0.66548164 1.26510936
## factor(time)4      0.79866879 1.56415206
## factor(time)5      0.29779345 0.83253478
## factor(time)6      0.42469720 1.10872589
## factor(time)7      0.49880474 1.33560856

# Pearson Chi-square statistic
chisq = sum( residuals(fit.cancer1, type="pearson")^2 )
chisq

## [1] 43.05481

# Deviance statistic
dev = sum( residuals(fit.cancer1, type="deviance")^2 )
dev

## [1] 43.92253

# Otra forma de calcular chisq
chisq = sum( ((cancer$count-fitted(fit.cancer1)) / sqrt(fitted(fit.cancer1)))^2)
chisq

## [1] 43.05481

Note que en este caso \(X^{2}=43.05\) y \(G^{2}=43.92\). En ambos casos, el cociente entre \(X^{2} / d f \circ G^{2} / d f\) es aproximadamente igual a \(0.83\) (en este modelo \(d f=52\) ). Esto indica que el ajuste del modelo es razonable dado que este cociente es cercano a 1. Cuando este cociente se es mucho menor de 1 (e.g., \(0.50\) ) se dice que hay subdispersión en los conteos (underdispersion), es decir, \(\operatorname{Var}\left(Y_{i}\right)<E\left(Y_{i}\right)\). Este escenario no es muy común en la práctica pero existen algunas distribuciones tales como la Poisson Generalizada que se pueden usar en estos casos (Hilbe 2014). Un ejemplo de regresión usando la Poisson Generalizada se encuentra en la ayuda de PROC GLIMMIX de \(5 A 5\). Otro tipo de modelo para subdispersión es la Regresión de Poisson Conway-Maxwell (COM-Poisson). Se puede usar la librería compoissonReg de \(\mathrm{R}\) para ajustar este tipo de modelo

En \(\mathrm{R}\) existe la función \(\mathrm{R}\) para llevar a cabo la prueba de hipótesis si un término es significativo o no en el i oy elo. Este proceso le da al investigador argumentos para llevar a cabo la selección de variables. Claro esta, siempre hay que ser conscientes de las comparaciones múltiples y su consecuencia en la inflación de los errores en las pruebas de hipótesis (esta función no permite tomar decisiones sobre la significancia de dos o más términos). Por ejemplo, en este la variable se histology se podría eliminar del modelo \(p>0.05\). La idea es reajustar el modelo sin ese término y volver a correr la función drop1 para el modelo nuevo. De la misma manera, existe la función add1 para agregar un término a la vez en el modelo y step para selección de modelos usando AIC

# Eliminando un término a la vez en el modelo
drop1(fit.cancer1, test="Chi")

Sobredispersión en Conteos.

Para ciertas situaciones, los conteos exhiben alta variabilidad tal que \(Var\left(Y_{i}\right)>E\left(Y_{i}\right)\), lo cual no satisface los supuestos de la distribución Poisson. Este fenómeno, bastante común en datos de conteos, se conoce como sobredispersión (overdispersion). Usualmente valores de \(X^{2} /(n-p) \circ G^{2} /(n-p)\) mucho mayores que 1 son señal de sobredispersión (algunos autores consideran valores de \(1.5\) en adelante)

Sobredispersión puede suceder por varias razones

La estructura de la media omite variables explicativas o predictores relevantes
Los datos contienen valores atípicos (outliers)
El modelo omite términos de interacción.
Un predictor necesita transformación (e.g., log)
La función de enlace no se especifica correctamente. Es decir, la relación entre la media transformada y los predictores no es lineal (poco frecuente en conteos)

Una consecuencia de ajustar una regresión con datos sobredispersos es que los errores estándar de los coeficientes se pueden subestimar. Esto trae como consecuencia que se declaren predictores significativos cuando en realidad no lo son (falsos positivos). Un mensaje importante es que se debe revisar detenidamente los datos para descartar posibles causas de sobredispersión antes de recurrir a las estrategias que vamos a mencionar para modelar datos de conteos sobredispersos.

Vamos a ilustrar el problema de sobredispersión en la regresión para conteos usando la base de datos del número de casos del virus de West Nile en equinos y aves en condados de Carolina del Sur (bajar base datos). Las variables explicativas son: PBR (Positive Bird Rate=# de aves/población) y densidad humana (HD). La variable respuesta es la tasa de casos en equinos (# de casos equinos/granjas) (ver descripción)

# Lee base de datos del sitio web
westnilesc <- read.csv("http://users.stat.ufl.edu/~winner/data/westnilesc.dat", sep="",header=TRUE)

names(westnilesc)<-c("count_y","state","bird","equine","farms","area",
                     "pop","humanden","pbr","per")
head(westnilesc)

str(westnilesc)

## 'data.frame':    45 obs. of  10 variables:
##  $ count_y : chr  "Aiken," "Allendale," "Anderson," "Bamberg," ...
##  $ state   : chr  "SC" "SC" "SC" "SC" ...
##  $ bird    : int  14 2 9 1 9 9 15 0 19 2 ...
##  $ equine  : int  6 0 1 2 1 0 2 2 2 0 ...
##  $ farms   : int  729 114 1271 254 325 99 292 261 266 412 ...
##  $ area    : int  1073 408 718 393 548 587 1098 380 919 393 ...
##  $ pop     : int  142557 11226 165719 16674 23472 120952 142704 15171 309997 52542 ...
##  $ humanden: num  132.9 27.5 230.8 42.4 42.8 ...
##  $ pbr     : num  9.82e-05 1.78e-04 5.43e-05 6.00e-05 3.83e-04 ...
##  $ per     : num  0.00823 0 0.000787 0.007874 0.003077 ...

Vamos a considerar un primer modelo con la variable densidad humana como la única predictora de la tasa promedio de equinos con el virus. Note que definimos el logaritmo del número de granjas como el offset para modelar una tasa

# Equine vs. PBR
ggplot(westnilesc, aes(x=humanden, y=equine))+geom_point()# Equine vs. PBR

# # Equine Cases =  Human Density with log(farms) as offset
westnilesc$log.farms = log(westnilesc$farms)
# Model without PBR
wn.fit1 = glm(equine ~  humanden, family=poisson, data=westnilesc, 
             offset=log.farms)
summary(wn.fit1)

## 
## Call:
## glm(formula = equine ~ humanden, family = poisson, data = westnilesc, 
##     offset = log.farms)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0761  -1.4395  -0.9995   0.5054   3.6755  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.8802118  0.2273985 -25.859   <2e-16 ***
## humanden    -0.0001143  0.0011885  -0.096    0.923    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 106.77  on 44  degrees of freedom
## Residual deviance: 106.76  on 43  degrees of freedom
## AIC: 163.29
## 
## Number of Fisher Scoring iterations: 6

# X^2 / df
chisq = sum( ((westnilesc$equine-fitted(wn.fit1)) / sqrt(fitted(wn.fit1)))^2)
chisq / 44

## [1] 2.752723

# G^2 / df
dev = sum( residuals(wn.fit1, type="deviance")^2 )
dev/44

## [1] 2.426345

Usando las dos cantidades \(X^{2} / d f \quad \text{y} \quad G^{2} /df\) vemos que ambas son mayores a \(2.5\), lo cual es una clara señal de sobredispersión. Antes de usar cualquiera de los métodos existentes para modelar sobredisperión es importante descartar algunas de las posibles causas antes de continuar con un modelo o estrategia más compleja. Primero, veamos los residuales del modelo

par(mfrow=c(2,2))
plot(fitted(wn.fit1), residuals(wn.fit1, type="deviance"))
plot(fitted(wn.fit1), residuals(wn.fit1, type="pearson"))
# Standardized Pearson residual
stpearson = resid(wn.fit1, type="pearson")/sqrt(1 - hatvalues(wn.fit1))
plot(fitted(wn.fit1), stpearson)
# Standardized deviance residual
stdeviance = rstandard(wn.fit1) # default
plot(fitted(wn.fit1), stdeviance)

Un segundo paso a seguir es verificar si se omitieron variables importantes en el modelo. Para este conjunto de datos vamos a incluir la variable PBR y su interacción con la densidad humana, tal como proponen los autores del artículo original.

# Equine vs. PBR
ggplot(westnilesc, aes(x=pbr, y=equine))+geom_point()

# Model with PBR and humanden and its interaction with log(farms) as offset
wn.fit2 = glm(equine ~  pbr*humanden + offset(log.farms), family=poisson, data=westnilesc)
summary(wn.fit2)

## 
## Call:
## glm(formula = equine ~ pbr * humanden + offset(log.farms), family = poisson, 
##     data = westnilesc)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6372  -1.0902  -0.7596   0.3453   2.4889  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -6.931e+00  3.697e-01 -18.748  < 2e-16 ***
## pbr           5.180e+03  1.608e+03   3.222  0.00128 ** 
## humanden      4.383e-04  1.687e-03   0.260  0.79503    
## pbr:humanden  2.630e+01  1.241e+01   2.119  0.03407 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 106.768  on 44  degrees of freedom
## Residual deviance:  61.086  on 41  degrees of freedom
## AIC: 121.62
## 
## Number of Fisher Scoring iterations: 6

# X^2 / df
chisq = sum( ((westnilesc$equine-fitted(wn.fit2)) / sqrt(fitted(wn.fit2)))^2)
chisq / 42

## [1] 1.663582

# G^2 / df
dev = sum( residuals(wn.fit2, type="deviance")^2 )
dev/42

## [1] 1.454429

Para este modelo los diagnósticos indican valores mayores a 1.5, lo cual podría ser un indicio de sobredispersión leve. Revisemos los residuales de este modelo.

par(mfrow=c(2,2))
plot(fitted(wn.fit2), residuals(wn.fit2, type="deviance"))
plot(fitted(wn.fit2), residuals(wn.fit2, type="pearson"))
# Standardized Pearson residual
stpearson = resid(wn.fit2, type="pearson")/sqrt(1 - hatvalues(wn.fit2))
plot(fitted(wn.fit2), stpearson)
# Standardized deviance residual
stdeviance = rstandard(wn.fit2) # default
plot(fitted(wn.fit2), stdeviance)

A pesar del ajuste razonable de este modelo vamos a explorar algunas alternativas para modelar datos de conteo sobredispersos.

Regresión Causi-Poisson

Una estrategia para modelar sobredispersión consiste en ajustar un modelo en el cual se asume que \(\operatorname{Var}\left(Y_{i}\right)=\phi \nu\left(\mu_{i}\right)\) para una constante \(\phi\). Cuando \(\phi>1\) entonces se dice que hay sobredispersión. Para la regresión con Poisson, \(\operatorname{Var}\left(Y_{i}\right)=\phi \mu_{i}>E\left(Y_{i}\right)=\mu_{i}\) cuando \(\phi>1\). La estimación de los parámetros \(\beta\) en la regresión se llevan cabo usando el método de cuasi-verosimilitud (quasi-likelihood). Los coeficientes en el modelo se estiman de igual manera que un GLM (verosimilitud) pero los errores estándar asintóticos de los coeficientes se multiplican por la cantidad \(\sqrt{\hat{\phi}}=\sqrt{X^{2} /(n-p)} .\) Es importante mencionar que este método corrige los errores estándar por sobredispersión pero no aborda el problema de especificación incorrecta de la componente lineal del modelo. Vamos a ajustar una regresión cuasi-Poisson o Poisson usando cuasi-verosimilitud para los datos del virus.

wn.fit3<-glm(formula = equine ~ pbr * humanden, family = quasipoisson, 
    data = westnilesc, offset = log.farms)
summary(wn.fit3)

## 
## Call:
## glm(formula = equine ~ pbr * humanden, family = quasipoisson, 
##     data = westnilesc, offset = log.farms)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6372  -1.0902  -0.7596   0.3453   2.4889  
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6.931e+00  4.826e-01 -14.361   <2e-16 ***
## pbr           5.180e+03  2.099e+03   2.468   0.0179 *  
## humanden      4.383e-04  2.203e-03   0.199   0.8432    
## pbr:humanden  2.630e+01  1.620e+01   1.623   0.1122    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.704157)
## 
##     Null deviance: 106.768  on 44  degrees of freedom
## Residual deviance:  61.086  on 41  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 6

Regresión Binomial Negativa (BN)

En los modelos anteriores asumimos que \(Y_{i} \sim\) Poisson \(\left(\mu_{i}\right)\). Como consecuencia de este supuesto distribucional tenemos que \(\operatorname{Var}\left(Y_{i}\right)=\mu_{i} .\) En esta sección vamos a modelar los conteos usando la distribución BN. La función de densi dad de esta distribución esta dada por \[ f\left(y_{i} ; \mu_{i}, k\right)=\frac{\Gamma(y+k)}{\Gamma(k) \Gamma(y+1)}\left(\frac{\mu_{i}}{\mu_{i}+k}\right)^{y}\left(\frac{k}{\mu_{i}+k}\right)^{k}, \quad y=0,1, \ldots \] Esta distribución resulta de la distribución margin: \(\underline{\alpha}_{\Rightarrow} Y_{i}\) si se asume que \(Y_{i} \mid \lambda_{i} \sim\) Porsson \(\left(\lambda_{i}\right) \mathrm{y}\) \(\lambda_{i} \sim Gama\left(k, k / \mu_{i}\right)\) (parametrización tradicional con \(E\left(\lambda_{i}\right)=\mu_{i} y Var\left(\lambda_{i}\right)=\mu_{i}^{2} / k\). Es importante recordar que la definición de las densidades, particularmente para la BN, cambia dependiendo de la referencia (la parametrización de la densidad de la BN en \(R\) es diferente a la de SAS. Inclusive dentro de SAS la parametrización de la densidad de la BN varía entre procedimientos) En la distribución BN se tiene entonces que \[ E\left(Y_{i}\right)=\mu_{i,} \quad var\left(Y_{i}\right)=\mu_{i}+\gamma \mu_{j^{2}}^{2} \] donde \(\gamma=1 / k\) se conoce como el parámetro de dispersión. Sin embargo, este parámetro no tiene la misma connotación del parámetro \(\phi\) en la definición de la familia exponencial de dispersión (se recomienda leer el siguiente documento de la ayuda de SAS ). Note que para \(\gamma>0\) se tiene que \(\operatorname{Var}\left(Y_{i}\right)>\mu_{i,}\) lo cual justifica el uso de la distribución BN para sobredispersión. Además, cuando \(\gamma\) tiende a 0 entonces la BN tiende a una Poisson. Es decir, si \(\phi\) no es significativamente diferente de 0 entonces este resultado sugiere que el modelo Poisson es mejor para el conjunto de datos. De la misma manera, se pueden llevar a cabo una prueba LRT para comparar el modelo Poisson y el BN, pero el valor p se debe ajustar debido a las características de los modelos En el caso de la regresión BN seguimos usando al misma función de enlace La regresión BN se ajusta en \(\mathrm{R}\) usando la función glm.nb de la librería mass . Es decir, esta distribución no esta incluida en la función glm que hemos venido utilizando Para este modelo la devianza esta dada por \[ D(\mathbf{y}, \hat{\mu})=2 \sum_{i}\left[y_{i} \log \left(\frac{y_{i}}{\hat{\mu}_{i}}\right)-\left(y_{i}+\frac{1}{\hat{\gamma}}\right) \log \left(\frac{1+\hat{\gamma} y_{i}}{1+\hat{\gamma} \hat{\mu}_{i}}\right)\right] \] Note que cuando \(\hat{\gamma}\) es muy pequeño entonces la devianza de la BN se acerca a aquella de la Poisson.

library(MASS)
# Model with PBR and interaction, likelihood, Negative Binomial
wn.fit4 = glm.nb(equine ~  pbr*humanden + offset(log.farms), data=westnilesc,
                 control=list(maxit=100, trace=1, epsilon=1e-14)) # no necesario en general

## Theta(1) = 2.847760, 2(Ls - Lm) = 46.092300

## Theta(2) = 2.515400, 2(Ls - Lm) = 44.905400

## Theta(3) = 2.515400, 2(Ls - Lm) = 44.901600

## Theta(4) = 2.466850, 2(Ls - Lm) = 44.710400

## Theta(5) = 2.466850, 2(Ls - Lm) = 44.710300

## Theta(6) = 2.458950, 2(Ls - Lm) = 44.678700

## Theta(7) = 2.458950, 2(Ls - Lm) = 44.678700

## Theta(8) = 2.457650, 2(Ls - Lm) = 44.673500

## Theta(9) = 2.457650, 2(Ls - Lm) = 44.673500

## Theta(10) = 2.457430, 2(Ls - Lm) = 44.672700

## Theta(11) = 2.457430, 2(Ls - Lm) = 44.672700

## Theta(12) = 2.457400, 2(Ls - Lm) = 44.672500

## Theta(13) = 2.457400, 2(Ls - Lm) = 44.672500

## Theta(14) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(15) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(16) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(17) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(18) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(19) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(20) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(21) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(22) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(23) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(24) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(25) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(26) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(27) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(28) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(29) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(30) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(31) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(32) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(33) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(34) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(35) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(36) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(37) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(38) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(39) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(40) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(41) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(42) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(43) = 2.457390, 2(Ls - Lm) = 44.672500

## Theta(44) = 2.457390, 2(Ls - Lm) = 44.672500

summary(wn.fit4)

## 
## Call:
## glm.nb(formula = equine ~ pbr * humanden + offset(log.farms), 
##     data = westnilesc, control = list(maxit = 100, trace = 1, 
##         epsilon = 1e-14), init.theta = 2.457391009, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2886  -1.0068  -0.7581   0.3880   1.9213  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -6.860e+00  4.378e-01 -15.672   <2e-16 ***
## pbr           4.344e+03  2.403e+03   1.807   0.0707 .  
## humanden     -8.230e-04  2.300e-03  -0.358   0.7205    
## pbr:humanden  4.459e+01  2.274e+01   1.961   0.0499 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.4574) family taken to be 1)
## 
##     Null deviance: 70.781  on 44  degrees of freedom
## Residual deviance: 44.672  on 41  degrees of freedom
## AIC: 120.97
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.46 
##           Std. Err.:  1.98 
## 
##  2 x log-likelihood:  -110.973

# Deviance G^2 
dev = sum( residuals(wn.fit4, type="deviance")^2 )
dev # devianza no escalada

## [1] 44.67248

dev/wn.fit4$theta # devianza escalada

## [1] 18.17882

Hay varias comentarios con respecto a este modelo

El cociente \(G^{2} / d f=45.246 / 42=1.08\), lo cual indica un buen ajuste
ELAIC del modelo BN es \(121.82\) comparado a \(122.54\) del modelo Poisson. Es decir, el modelo BN es ligeramente mejor que el modelo Poisson.
Sin embargo, la estimación del parámetro de dispersión \(\gamma\) de la BN es igual \(a \hat{\gamma}=2.39(1.91)\). Note que el error estándar es grande comparado a la estimación puntual. Esto sugiere que hay nucha incertidumbre en la estimación del parámetro de dispersión, por lo que el modelo Poisson podría ser un buen candidato para modelar este conjunto de datos
Más aún, la prueba LRT entre el modelo Poisson \(\left(H_{0}\right)\) y la BN \(\left(H_{1}\right)\) sugiere el modelo BN es preferido (&) para este conjunto de datos \((p<0.0495)\). Esta conclusión se debe tomar con precaución debido a la cercanía del valor \(p\) al nivel de significancia \(.05\). Una recomendación personal es usar el modelo Poisson el cual es más sencillo y ajusta razonablemente al conjunto de datos

anova(wn.fit2, wn.fit4) # DO NOT use it!

# LRT test Poisson vs. Negative Binomial
lrt  = -2*(logLik(wn.fit2) - logLik(wn.fit4))
lrt

## 'log Lik.' 2.645644 (df=4)

# p-value is 0.5*p-value of the test with a df=1 due to 
# the fact that the Poisson is a boundary case of the NB distribution (gamma=0)
0.5*pchisq(lrt, df=1, lower.tail=F) # p-value

## 'log Lik.' 0.05191703 (df=4)

# Resultado similar usando función pdTest()  de la librería pscl
library(pscl)

## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis

odTest(wn.fit4)

## Likelihood ratio test of H0: Poisson, as restricted NB model:
## n.b., the distribution of the test-statistic under H0 is non-standard
## e.g., see help(odTest) for details/references
## 
## Critical value of test statistic at the alpha= 0.05 level: 2.7055 
## Chi-Square Test Statistic =  2.6456 p-value = 0.05192

A manera de conclusión, podemos decir que el modelo Poisson es un buen modelo para describir la tasa de casos del virus en función de las covariables. Note que a pesar de que la verosimilitud en el modelo BN es más compleja, la interpretación de los \(\beta^{\prime}\) sen ambos modelos (Poisson o BN) es la misma ya que esta depende de la función de enlace

El fenómeno de sobredispersión se puede abordar de otras maneras tales como distribuciones diferentes para la respuesta o modelos mixtos. En cuanto a otras distribuciones diferentes a la Poisson y BN no serán discutidas en este curso. Con base en la experiencia, la regresión BN suele ser muy versátil para manejar casos de sobredispersión bajo diferentes escenarios. Es decir, rara vez es necesario recurrir a distribuciones más complejas que podrían complicar la estimación de los parámetros del modelo.

Exceso de Ceros en los Datos

En muchas situaciones la proporción de ceros en los conteos es considerable. La presencia de demasiados ceros suele ocasionar sobredispersión que en muchas circunstancias se pueden modelar con la regresión BN. Sin embargo, si hay demasiados ceros en ocasiones es mejor considerar un modelo que tome en cuenta el proceso o los procesos que generan los ceros tales como el ZIP modelo Zero-inflated Poisson model (los ceros pueden resultar de dos procesos separados), o un modelo de hurdle (los ceros resultan de un único proceso). Tanto sas como \(\mathrm{R}\) permite el ajuste de este tipo de modelos. El uso de un modelo en particular dependerá de las caracteristicas del problema investigado. En cualquier caso, los modelos se definen usando los MLG que se darán en este curso En un modelo con exceso de ceros se asume que hay dos procesos: uno que genera solo ceros y otro que genera conteos de una distribución \(f\left(y_{i}\right)\). Este último también incluye ceros. Esto es, se modela \(y_{i}=0\) con probabilidad \(\psi_{i} y y_{i}>0\) con probabilidad \(1-\psi_{i}\). Por lo tanto, la función de probabilidad de \(Y_{i}\) dado los vectores de covariables \(\mathbf{x}_{i} y \mathbf{z}_{i}\) se puede escribir como \[ \begin{array}{c} P\left(Y_{i}=0 \mid \mathrm{x}_{i}, \mathrm{z}_{i}\right)=\psi_{i}+\left(1-\psi_{i}\right) f(0) \\ P\left(Y_{i}=y_{i} \mid \mathrm{x}_{i}, \mathrm{z}_{i}\right)=\left(1-\psi_{i}\right) f\left(y_{i}\right), \quad y_{i}>0 \end{array} \]

En el caso de un modelo ZIP \(f\left(y_{i}\right)\) es una Poisson y para un modelo ZINB \(f\left(y_{i}\right)\) es una binomial negativa Lenalmente, la probabilidad de tomar el valor cero \(\psi\) se modela usando una función de distribución \(F\) (por efiniplo, la función probit \(\Phi(t) \circ\) logit \(\log (t /(1-t)))\) \[ P\left(Y_{i}=0\right)=F\left(\mathrm{z}_{i} \gamma\right) \] En el caso del modelo \(Z \mid P\) \[ \mu_{i}=e^{\mathrm{x}_{4}} \] así que el modelo se define como \[ \begin{array}{c} P\left(Y_{i}=0 \mid \mathrm{x}_{i}, \mathrm{z}_{i}\right)=\psi_{i}+\left(1-\psi_{i}\right) e^{-\mu_{4}} \\ P\left(Y_{i}=y_{i} \mid \mathrm{x}_{i}, \mathrm{z}_{i}\right)=\left(1-\psi_{i}\right) \frac{e^{-\mu_{i}} \mu_{i}^{y_{i}}}{y_{i} !}, \quad y_{i}>0 . \end{array} \] De lo anterior, se tiene que \[ E\left(Y_{i} \mid \mathrm{x}_{i,} \mathrm{z}_{i}\right)=\mu_{i}\left(1-\psi_{i}\right) \] \[ Var\left(Y_{i} \mid \mathbf{x}_{i}, \mathrm{z}_{i}\right)=E\left(Y_{i} \mid \mathrm{x}_{i}, \mathbf{z}_{i}\right)\left(1+\mu_{i} \psi_{i}\right) \] El logaritmo de la verosimilitud se maximiza para obtener las estimaciones de los parámetros en cada modelo

El logaritmo de la verosimilitud se maximiza para obtener las estimaciones de los parámetros en cada modelo Q. Le otro tipo de modelos llamados Hurdle Models. Sin embargo, en este caso se modelan los ceros y los valores positivos por separado. Siguiendo con la distribución Poisson, un modelo de hurdle se define como \[ \begin{array}{c} P\left(Y_{i}=0 \mid \mathrm{x}_{i\rangle} \mathrm{z}_{i}\right)=\psi_{i} \\ P\left(Y_{i}=y_{i}\left|\mathbf{x}_{i}\right\rangle \mathbf{z}_{i}\right)=\left(1-\psi_{i}\right) \frac{e^{-\mu_{i}} \mu_{i}^{y_{i}}}{\left(1-e^{\mu_{4}}\right) y_{i} !}, \quad y_{i}>0 . \end{array} \] Note que en la segunda parte, a diferencia del modelo \(Z I P\), se define una distribución Poisson truncada ya que no se permiten ceros del segundo proceso \(f\left(y_{i}\right)\).

Por ejemplo, suponga que usted se enfrenta a una situación de compra en la cuál primero decide si compra o no un producto y luego decide cuántos. Si una vez toma decisión de comprar tiene que comprar al menos un producto entonces se puede usar el modelo de hurdle para modelar esta situación. Sin embargo, si después de decidir que va a comprar le permiten retractarse y comprar cero productos entonces debería usar un modelo similar al ZIP (los ceros pueden generarse de dos procesos)

Regresion Poisson Y Binomial Negativa: Numero de muertes por cancer de pulmón