#install.packages(c("AER", "robust", "gcc"))CLASE 13. Regresión de Poisson y Binomial Negativa en R
Necesitamos instalar la paquetería AER, robust, gcc
REGRESIÓN DE POISSON
Visualizar la base de datos
data(breslow.dat, package="robust")
names(breslow.dat) [1] "ID" "Y1" "Y2" "Y3" "Y4" "Base" "Age" "Trt" "Ysum"
[10] "sumY" "Age10" "Base4"summary(breslow.dat[c(6, 7, 8, 10)]) Base Age Trt sumY
Min. : 6.00 Min. :18.00 placebo :28 Min. : 0.00
1st Qu.: 12.00 1st Qu.:23.00 progabide:31 1st Qu.: 11.50
Median : 22.00 Median :28.00 Median : 16.00
Mean : 31.22 Mean :28.34 Mean : 33.05
3rd Qu.: 41.00 3rd Qu.:32.00 3rd Qu.: 36.00
Max. :151.00 Max. :42.00 Max. :302.00 head(breslow.dat) ID Y1 Y2 Y3 Y4 Base Age Trt Ysum sumY Age10 Base4
1 104 5 3 3 3 11 31 placebo 14 14 3.1 2.75
2 106 3 5 3 3 11 30 placebo 14 14 3.0 2.75
3 107 2 4 0 5 6 25 placebo 11 11 2.5 1.50
4 114 4 4 1 4 8 36 placebo 13 13 3.6 2.00
5 116 7 18 9 21 66 22 placebo 55 55 2.2 16.50
6 118 5 2 8 7 27 29 placebo 22 22 2.9 6.75- Variables Independiente: Base, Age, Trt
- Variable dependiente: sumY
Distribución gráfica de los recuentos de convulsiones posteriores al tratamiento
opar <- par(no.readonly=TRUE)
par(mfrow=c(1, 2))
attach(breslow.dat)
hist(sumY, breaks=20, xlab="Seizure Count",
main="Distribution of Seizures")
boxplot(sumY ~ Trt, xlab="Treatment", main="Group Comparisons")
par(opar)Ajuste de Regresión (Pronósticos)
fit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson())
summary(fit)
Call:
glm(formula = sumY ~ Base + Age + Trt, family = poisson(), data = breslow.dat)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.9488259 0.1356191 14.370 < 2e-16 ***
Base 0.0226517 0.0005093 44.476 < 2e-16 ***
Age 0.0227401 0.0040240 5.651 1.59e-08 ***
Trtprogabide -0.1527009 0.0478051 -3.194 0.0014 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 2122.73 on 58 degrees of freedom
Residual deviance: 559.44 on 55 degrees of freedom
AIC: 850.71
Number of Fisher Scoring iterations: 5Interpretar los parámetros del modelo
coef(fit) # logarítmica (Intercept) Base Age Trtprogabide
1.94882593 0.02265174 0.02274013 -0.15270095 exp(coef(fit)) # exponencial (Intercept) Base Age Trtprogabide
7.0204403 1.0229102 1.0230007 0.8583864 Evaluar la sobredispersión
deviance(fit)/df.residual(fit)[1] 10.1717library(qcc)Package 'qcc' version 2.7Type 'citation("qcc")' for citing this R package in publications.qcc.overdispersion.test(breslow.dat$sumY, type="poisson")
Overdispersion test Obs.Var/Theor.Var Statistic p-value
poisson data 62.87013 3646.468 0REGRESIÓN CUASI-POISSON
Con los datos existentes, hacemos la aplicación del modelo, para ajustar a la estrucutra de un pronóstico con Cuasi-Poisson
Ajustar el modelo con Cuasipoisson (Pronóstico)
fit.od <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,
family=quasipoisson())
summary(fit.od)
Call:
glm(formula = sumY ~ Base + Age + Trt, family = quasipoisson(),
data = breslow.dat)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.948826 0.465091 4.190 0.000102 ***
Base 0.022652 0.001747 12.969 < 2e-16 ***
Age 0.022740 0.013800 1.648 0.105085
Trtprogabide -0.152701 0.163943 -0.931 0.355702
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 11.76075)
Null deviance: 2122.73 on 58 degrees of freedom
Residual deviance: 559.44 on 55 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5anova(fit,fit.od,test='LR')Analysis of Deviance Table
Model 1: sumY ~ Base + Age + Trt
Model 2: sumY ~ Base + Age + Trt
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 55 559.44
2 55 559.44 0 0
Nota
El paso siguiente es intentar ver como se adecua el modelo, en vez de a una variable de Poisson, a una variable binomial negativa. La definición de variable binomial negativa aparece en cualquier manual de estadística básica. A nosotros nos interesa que es una variable que, como la de Poisson, en principio toma valores desde 0 en adelante, pero la varianza puede ser bastante más grande que la media (que es lo que nos ocurre).
REGRESIÓN BINOMIAL NEGATIVA
Para estimar el modelo, debemos cambiar el comando glm por glm.nb.
library(MASS)
fm_nbin <- glm.nb(sumY ~ Base + Age + Trt, data=breslow.dat)
summary(fm_nbin)
Call:
glm.nb(formula = sumY ~ Base + Age + Trt, data = breslow.dat,
init.theta = 3.361505307, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.931832 0.408372 4.731 2.24e-06 ***
Base 0.027736 0.002817 9.846 < 2e-16 ***
Age 0.016815 0.012535 1.341 0.18
Trtprogabide -0.193518 0.154271 -1.254 0.21
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(3.3615) family taken to be 1)
Null deviance: 181.856 on 58 degrees of freedom
Residual deviance: 63.667 on 55 degrees of freedom
AIC: 476.46
Number of Fisher Scoring iterations: 1
Theta: 3.362
Std. Err.: 0.710
2 x log-likelihood: -466.457 De entrada, comparando las devianzas y el AIC de este modelo con los términos correspondientes en el modelo de Poisson, se reducen bastante los números. El AIC pasa de 850.71 a 476.46, con lo que este último modelo es mejor que el anterior.
Observamos que no aparecen como variables significativas ni Age ni Trtprogabide. Vamos entonces a ver si mejora el modelo quitando estas variables (ahora debemos escribirlas, puesto que no son todas las del conjunto de datos que estamos usando):
fm_nbin2<- glm.nb(sumY ~ Base + Age , data=breslow.dat)
summary(fm_nbin2)
Call:
glm.nb(formula = sumY ~ Base + Age, data = breslow.dat, init.theta = 3.266828461,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.778076 0.394931 4.502 6.72e-06 ***
Base 0.027770 0.002856 9.723 < 2e-16 ***
Age 0.018761 0.012602 1.489 0.137
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(3.2668) family taken to be 1)
Null deviance: 177.407 on 58 degrees of freedom
Residual deviance: 63.771 on 56 degrees of freedom
AIC: 476.02
Number of Fisher Scoring iterations: 1
Theta: 3.267
Std. Err.: 0.687
2 x log-likelihood: -468.017 fm_nbin3<- glm.nb(sumY ~ Base , data=breslow.dat)
summary(fm_nbin3)
Call:
glm.nb(formula = sumY ~ Base, data = breslow.dat, init.theta = 3.12053628,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.327085 0.122420 19.009 <2e-16 ***
Base 0.027418 0.002858 9.595 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(3.1205) family taken to be 1)
Null deviance: 170.474 on 58 degrees of freedom
Residual deviance: 63.519 on 57 degrees of freedom
AIC: 476.12
Number of Fisher Scoring iterations: 1
Theta: 3.121
Std. Err.: 0.647
2 x log-likelihood: -470.116 Es mejor el modelo fm_nbin2
anova(fm_nbin2,fm_nbin3)Likelihood ratio tests of Negative Binomial Models
Response: sumY
Model theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
1 Base 3.120536 57 -470.1163
2 Base + Age 3.266828 56 -468.0168 1 vs 2 1 2.099421 0.1473549Los resultados de AIC y devianzas son muy similares al modelo anterior. Vamos a hacer un test ANOVA para comparar entre los dos modelos
anova(fm_nbin2,fm_nbin3,test = 'Chisq')Likelihood ratio tests of Negative Binomial Models
Response: sumY
Model theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
1 Base 3.120536 57 -470.1163
2 Base + Age 3.266828 56 -468.0168 1 vs 2 1 2.099421 0.1473549El p_{value} nos dice que los modelos no pueden considerarse diferentes.
Modelos con exceso de ceros
El siguiente paso que se debe dar es considerar los llamados modelos que tienen “excesos de ceros”. Se utilizan cuando, como en nuestro caso, la frecuencia de aparición de ceros es bastante más grande que la frecuencia de aparición del resto de valores. Sin embargo, hay ocasiones en que tienen mucho mayor interés. Por ejemplo, el número de partes a un seguro al año por sus asegurados (la mayoría son cero y las probabilidades de valores mayores van disminuyendo). En seguros de salud, las visitas al centro hospitalario. El número de muertos por una enfermedad al año entre quienes la padecieron.
Este tipo de situaciones se modelizan con los llamados “Zero inflated” models (modelos inflados de ceros) y los “Hurdle models” (modelos con salto). Un modelo con salto es un modelo con exceso de ceros donde además puede existir una diferencia sustancial entre el valor cero y el resto de valores. Por ejemplo, si se realiza una encuesta en una población sobre el número de cigarrillos consumidos en una semana, hay dos subpoblaciones: los que fuman y los que no. Los que no fuman siempre van a dar el valor cero, y los que fuman pueden dar el valor cero o mayor. Si consideramos un grupo de personas con discapacidad y otro sin ella, el valor cero de caidas sería posiblemente diferente para cada subpoblación y habría que estudiarlo.
Precaución
Los modelos con exceso de ceros se pueden analizar de manera equivalente a la que hemos hecho. Al estudiarlo con estos datos, no se han observado mejoras con respecto al modelo que utiliza la variable binomial negativa.
Dado que, por un lado, estos modelos son matemáticamente más complejos, suele ser necesario un número mayor de datos del que tenemos en esta ocasión para poder realizar estimaciones fiables. Por otro lado, al no haber obtenido mejoras con el modelo de variable binomial negativa, el principio de parsimonia (si se tienen varios modelos equivalentes, nos quedamos con el más sencillo) nos dice que nos quedemos con este último.
Importante
No sería necesario, puesto que hemos visto ya que el modelo de binomial negativa es mejor que el de Poisson, pero vamos a realizar una comparación entre ellos (Poisson, Quasi-Poisson, fm_nbin y fm_nbin2):
zero-countsb
round(c("Obs" = sum(breslow.dat$sumY < 1),
"ML-Pois" = sum(dpois(0, fitted(fit))),
"NB" = sum(dnbinom(0, mu = fitted(fm_nbin),
size = fm_nbin$theta)),
"NB_2" = sum(dnbinom(0, mu = fitted(fm_nbin2),
size = fm_nbin$theta))
)) Obs ML-Pois NB NB_2
1 0 0 0 Esto simplemente para comparar los AIC de los modelos y elegir el menor
tmp <- c( AIC(fit), AIC(fit.od), AIC(fm_nbin), AIC(fm_nbin3), AIC(fm_nbin2))
( BestAIC<- which.min(tmp) )[1] 5Considerando como mejor modelo de pronóstico fm\_nbin2, debido a que tiene el menor AIC y no existen ceros inflados, entonces ejecutamos lo siguiente:
#DATOS REALES
head(breslow.dat) ID Y1 Y2 Y3 Y4 Base Age Trt Ysum sumY Age10 Base4
1 104 5 3 3 3 11 31 placebo 14 14 3.1 2.75
2 106 3 5 3 3 11 30 placebo 14 14 3.0 2.75
3 107 2 4 0 5 6 25 placebo 11 11 2.5 1.50
4 114 4 4 1 4 8 36 placebo 13 13 3.6 2.00
5 116 7 18 9 21 66 22 placebo 55 55 2.2 16.50
6 118 5 2 8 7 27 29 placebo 22 22 2.9 6.75Prueba con modelo fm\_nbin2
prediccion=predict(fm_nbin2,data.frame(Base=11,Age=31))
prediccion #sumY (logarítmica) 1
2.665121 exp(2.665121)[1] 14.36969prediccion=predict(fm_nbin2,data.frame(Base=6,Age=25))
prediccion #sumY (logarítmica) 1
2.413709 exp(2.413709)[1] 11.17533# PRUEBA
p_pron=predict(fm_nbin2,data.frame(Base=14,Age=36))
exp(p_pron) 1
17.15403 Prueba con modelo fit
pred=predict(fit,data.frame(Base=11,Age=31,Trt="placebo"))
pred #sumY (logarítmica) 1
2.902939 exp(2.750238)[1] 15.64636pred=predict(fit,data.frame(Base=6,Age=25,Trt="placebo"))
pred #sumY (logarítmica) 1
2.65324 exp(2.500539)[1] 12.18906# PRUEBA
p_pronos=predict(fit,data.frame(Base=14,Age=36,Trt="progabide"))
exp(p_pronos) 1
18.76314