Chapter 11 Modelado de regresión de Poisson utilizando datos de tasas
Hasta ahora, en este tutorial, hemos modelado datos de recuento, pero también podemos modelar datos de tasa que predicen el número de recuentos durante un período de tiempo o agrupación. La fórmula para modelar datos de tasa viene dada por: \[\log (\frac{X}{n})=\beta_{0}+\sum \beta_i X_{i}\] esto es equivalente a:
\[\log (X)-\log (n)=\beta_{0}+\sum \beta_{i} X_{i}\]
\[\log (X)=\log (n)+\beta_{0}+\sum \beta_i X_{i}\] Por lo tanto, los datos de tasa se pueden modelar incluyendo el término log(n) con un coeficiente de 1.
Esto se denomina compensación. Este desplazamiento se modela con offset()en R.
Usemos otro conjunto de datos llamado eba1977 del paquete ISwR para modelar el modelo de regresión de Poisson para datos de tasas. Primero, instalaremos el paquete:
# install.packages("ISwR")
library(ISwR)Ahora, echemos un vistazo a algunos detalles sobre los datos e imprimamos las primeras diez filas para tener una idea de lo que incluye el conjunto de datos.
data(eba1977)
cancer.data = eba1977
cancer.data[1:10, ]## city age pop cases
## 1 Fredericia 40-54 3059 11
## 2 Horsens 40-54 2879 13
## 3 Kolding 40-54 3142 4
## 4 Vejle 40-54 2520 5
## 5 Fredericia 55-59 800 11
## 6 Horsens 55-59 1083 6
## 7 Kolding 55-59 1050 8
## 8 Vejle 55-59 878 7
## 9 Fredericia 60-64 710 11
## 10 Horsens 60-64 923 15# Description
# Lung cancer incidence in four Danish cities 1968-1971
# Description:
# This data set contains counts of incident lung cancer cases and
# population size in four neighbouring Danish cities by age group.
# Format:
# A data frame with 24 observations on the following 4 variables:
# city a factor with levels Fredericia, Horsens, Kolding, and Vejle.
# age a factor with levels 40-54, 55-59, 60-64, 65-69,70-74, and 75+.
# pop a numeric vector, number of inhabitants.
# cases a numeric vector, number of lung cancer cases.Para modelar datos de tasa, usamos \(\dfrac{X}{n}\) donde X es el evento que sucederá y n es la agrupación. En este ejemplo, X = casos (el evento es un caso de cáncer) y n = pop (la población es la agrupación).
Como en la fórmula anterior, los datos de tasa se contabilizan mediante log (n) y en estos datos n es la población, por lo que primero encontraremos el log de la población. Podemos modelar para casos / población de la siguiente manera:
# find the log(n) of each value in 'pop' column. It is the third column
logpop = log(cancer.data[ ,3])
# add the log values to the dataframe using 'cbind()'
new.cancer.data = cbind(cancer.data, logpop)
# display new dataframe
new.cancer.data## city age pop cases logpop
## 1 Fredericia 40-54 3059 11 8.025843
## 2 Horsens 40-54 2879 13 7.965198
## 3 Kolding 40-54 3142 4 8.052615
## 4 Vejle 40-54 2520 5 7.832014
## 5 Fredericia 55-59 800 11 6.684612
## 6 Horsens 55-59 1083 6 6.987490
## 7 Kolding 55-59 1050 8 6.956545
## 8 Vejle 55-59 878 7 6.777647
## 9 Fredericia 60-64 710 11 6.565265
## 10 Horsens 60-64 923 15 6.827629
## 11 Kolding 60-64 895 7 6.796824
## 12 Vejle 60-64 839 10 6.732211
## 13 Fredericia 65-69 581 10 6.364751
## 14 Horsens 65-69 834 10 6.726233
## 15 Kolding 65-69 702 11 6.553933
## 16 Vejle 65-69 631 14 6.447306
## 17 Fredericia 70-74 509 11 6.232448
## 18 Horsens 70-74 634 12 6.452049
## 19 Kolding 70-74 535 9 6.282267
## 20 Vejle 70-74 539 8 6.289716
## 21 Fredericia 75+ 605 10 6.405228
## 22 Horsens 75+ 782 2 6.661855
## 23 Kolding 75+ 659 12 6.490724
## 24 Vejle 75+ 619 7 6.428105Ahora, modelemos los datos de la tasa con offset()
poisson.model.rate<-glm(cases ~ city + age+ offset(logpop), family = poisson(link = "log"), data = cancer.data)
#display summary
summary(poisson.model.rate)##
## Call:
## glm(formula = cases ~ city + age + offset(logpop), family = poisson(link = "log"),
## data = cancer.data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.63573 -0.67296 -0.03436 0.37258 1.85267
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.6321 0.2003 -28.125 < 2e-16 ***
## cityHorsens -0.3301 0.1815 -1.818 0.0690 .
## cityKolding -0.3715 0.1878 -1.978 0.0479 *
## cityVejle -0.2723 0.1879 -1.450 0.1472
## age55-59 1.1010 0.2483 4.434 9.23e-06 ***
## age60-64 1.5186 0.2316 6.556 5.53e-11 ***
## age65-69 1.7677 0.2294 7.704 1.31e-14 ***
## age70-74 1.8569 0.2353 7.891 3.00e-15 ***
## age75+ 1.4197 0.2503 5.672 1.41e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 129.908 on 23 degrees of freedom
## Residual deviance: 23.447 on 15 degrees of freedom
## AIC: 137.84
##
## Number of Fisher Scoring iterations: 5En este conjunto de datos, podemos ver que la desviación residual está cerca de los grados de libertad y el parámetro de dispersión es 1.5 (23.447 / 15) que es pequeño, por lo que el modelo es un buen ajuste.
Usamos fitted(model) para devolver valores ajustados por el modelo. Devuelve resultados utilizando los datos de entrenamiento sobre los que se construye el modelo. Hagamos un intento:
fitted(poisson.model.rate)## 1 2 3 4 5 6 7 8
## 10.954812 7.411803 7.760169 6.873215 8.615485 8.384458 7.798635 7.201421
## 9 10 11 12 13 14 15 16
## 11.609373 10.849479 10.092831 10.448316 12.187276 12.576313 10.155638 10.080773
## 17 18 19 20 21 22 23 24
## 11.672630 10.451942 8.461440 9.413988 8.960422 8.326004 6.731286 6.982287Usando este modelo, podemos predecir el número de casos por 1000 habitantes para un nuevo conjunto de datos, usando la función predict (), de manera muy similar a como lo hicimos para nuestro modelo de conteo de datos anteriormente:
# create a test dataframe containing new values of variables
test.data = data.frame(city = "Kolding", age = "40-54", pop = 1000, logpop = log(1000))
# predict outcomes (responses) using 'predict()'
predicted.value<-predict(poisson.model.rate, test.data, type = "response")
# show predicted value
predicted.value## 1
## 2.469818Entonces, para la ciudad de Kolding entre las personas en el grupo de edad de 40 a 54 años, podríamos esperar aproximadamente 2 o 3 casos de cáncer de pulmón por cada 1000 personas.