En este ejercicio se busca replicar los resultados obtenidos por Favi, et. al. sobre los casos de rabia en Chile de 1989 a 2005, publicados en la revista [Revista Chilena de Infectologia][http://ref.scielo.org/jws58f].
Básicamente la idea es replicar la siguiente tabla:
Sin embargo está tabla no se pudo reproducir con los mismos resultados, adicionalmente cabe señalar que algo extraño que tiene el artículo es que utiliza como variable respuesta la proporción de murcielagos con rabia, sin embargo la distribución de probabilidad utilizada corresponde a datos de conte. Por lo cual en este ejercicio se propone una metodología diferente. Con los mismos datos:
año <- 1989:2005
p_sos <- c(3, 5, 7, 7, 12, 15, 16, 48, 22, 9, 38, 67, 63, 108, 73, 81, 103)
t_sos <- c(15, 267, 92, 74, 55, 98, 270, 467, 301, 265, 76, 848, 741, 878, 744,
896, 1042)
p_vig <- c(0, 4, 1, 1, 5, 2, 2, 9, 8, 0, 0, 0, 0, 3, 1, 5, 1)
t_vig <- c(291, 164, 314, 305, 206, 61, 347, 756, 837, 625, 1510, 509, 337,
637, 694, 922, 747)
Supongamos que el interes de estudio consiste en observar si la proporción de murcielagos con rabia a aumentado, o cambiado, en el tiempo:
plot(año, (p_sos + p_vig)/(t_sos + t_vig), type = "b", pch = 16, col = 2, ylab = "Proporción de murcielagos con rabia")
Sin considerar autocorrelación en los datos, la primera opción consiste en utilizar un modelo de regresión logística, sin embargo a la manera en como se encuentra agrupada la información esto no es posible por lo cual se utiliza la función expand.table en el paquete epitools para generar un “data.frame” con un registro para cada observación:
library(epitools)
datos <- matrix(nrow = 17, ncol = 2, dimnames = list(año = año, rabia = c(1,
0)))
datos[, 1] <- p_sos + p_vig
datos[, 2] <- t_sos + t_vig - (p_sos + p_vig)
datos #esta es la tabla con los datos agrupados
## rabia
## año 1 0
## 1989 3 303
## 1990 9 422
## 1991 8 398
## 1992 8 371
## 1993 17 244
## 1994 17 142
## 1995 18 599
## 1996 57 1166
## 1997 30 1108
## 1998 9 881
## 1999 38 1548
## 2000 67 1290
## 2001 63 1015
## 2002 111 1404
## 2003 74 1364
## 2004 86 1732
## 2005 104 1685
library(epitools)
expandidos <- expand.table(datos) #con esta función expandimos la tabla de contingencia a una versión completa por registro
# mostramos los primeros diez registros
expandidos[1:10, ]
## año rabia
## 1 1989 1
## 2 1989 1
## 3 1989 1
## 4 1989 0
## 5 1989 0
## 6 1989 0
## 7 1989 0
## 8 1989 0
## 9 1989 0
## 10 1989 0
A partir de los datos expandidos ajustamos el modelo de regresión logística:
m_log_exp <- glm(rabia ~ as.integer(año), family = binomial(logit), data = expandidos)
summary(m_log_exp)
##
## Call:
## glm(formula = rabia ~ as.integer(año), family = binomial(logit),
## data = expandidos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.769 0.263 0.299 0.329 0.351
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.87683 0.13023 29.77 < 2e-16 ***
## as.integer(año) -0.06587 0.00992 -6.64 3.1e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5902.1 on 16390 degrees of freedom
## Residual deviance: 5854.6 on 16389 degrees of freedom
## AIC: 5859
##
## Number of Fisher Scoring iterations: 6
Otra alternativa consiste en utilizar la opción de pesos en la función de modelos lineales generalizados, pero para esto es necesario presentar de otra manera los datos, de manera similar como una tabla de contingencia sin embargo ahora ponemos la variable rabia tambien en los renglones, y utilizamos las frecuencias como pesos
datos_2 <- matrix(nrow = 17 * 2, ncol = 3, dimnames = list(rep(1:17, 2), c("año",
"rabia", "Freq")))
datos_2[, 1] <- c(año, año)
datos_2[, 2] <- c(rep(1, 17), rep(0, 17))
datos_2[, 3] <- c(datos[, 1], datos[, 2])
datos_2 #esta tabla es la que se utiliza para ajustar el modelo.
## año rabia Freq
## 1 1989 1 3
## 2 1990 1 9
## 3 1991 1 8
## 4 1992 1 8
## 5 1993 1 17
## 6 1994 1 17
## 7 1995 1 18
## 8 1996 1 57
## 9 1997 1 30
## 10 1998 1 9
## 11 1999 1 38
## 12 2000 1 67
## 13 2001 1 63
## 14 2002 1 111
## 15 2003 1 74
## 16 2004 1 86
## 17 2005 1 104
## 1 1989 0 303
## 2 1990 0 422
## 3 1991 0 398
## 4 1992 0 371
## 5 1993 0 244
## 6 1994 0 142
## 7 1995 0 599
## 8 1996 0 1166
## 9 1997 0 1108
## 10 1998 0 881
## 11 1999 0 1548
## 12 2000 0 1290
## 13 2001 0 1015
## 14 2002 0 1404
## 15 2003 0 1364
## 16 2004 0 1732
## 17 2005 0 1685
Y el modelo es el siguiente
m_log_tab <- glm(rabia ~ año, family = binomial(logit), weights = Freq, data = as.data.frame(datos_2))
summary(m_log_tab)
##
## Call:
## glm(formula = rabia ~ año, family = binomial(logit), data = as.data.frame(datos_2),
## weights = Freq)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.401 -8.805 0.931 13.372 25.831
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.35e+02 1.98e+01 -6.79 1.1e-11 ***
## año 6.59e-02 9.92e-03 6.64 3.1e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5902.1 on 33 degrees of freedom
## Residual deviance: 5854.6 on 32 degrees of freedom
## AIC: 5859
##
## Number of Fisher Scoring iterations: 6
Contrario a lo que esperabamos, los resultados no son los mismos, es decir los esperados, de manera análoga a lo que sucedio con la opción de pesos para la funcion de modelo lineal en un ejercicio anterior. Sería necesario revisar la documentación, y probablemente directamente el código de la función, sin embargo no es el objetivo de este ejercicio.
Ahora supongamos que el objetivo del estudio no consiste en saber si hay una mayor proporción de murcielagos con rabia, sino saber si hay más murcielagos con rabia, lo cual puede ser una pregunta más relevante en términos de salud pública, sin considerar los tamaños de muestra tenemos los siguientes conteos de murcielagos con rabia:
plot(año, (p_sos + p_vig), type = "b", pch = 16, col = 4, ylab = "Murcielagos con rabia")
De entrada notamos que la tendencia en este caso es muy clara. Ajustemos pues el modelo sin considerar el tamaño de muestra:
library(MASS)
m_bn_t <- glm.nb(datos[, 1] ~ año)
summary(m_bn_t)
##
## Call:
## glm.nb(formula = datos[, 1] ~ año, init.theta = 6.994900521,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.569 -0.583 -0.278 0.227 2.339
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -368.9320 44.0678 -8.37 <2e-16 ***
## año 0.1864 0.0221 8.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(6.995) family taken to be 1)
##
## Null deviance: 86.097 on 16 degrees of freedom
## Residual deviance: 17.088 on 15 degrees of freedom
## AIC: 138.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 6.99
## Std. Err.: 3.04
##
## 2 x log-likelihood: -132.88
…Sin embargo quizá vale la pena considerar los tamaños de muestra, ¿qué pasa cuando vemos el número de conteos frente al tamaño de muestra?..
plot(año, (t_sos + t_vig), type = "b", pch = 16, col = 3, ylab = "Murcielagos",
ylim = c(0, max(t_sos + t_vig)))
lines(año, (p_sos + p_vig), type = "b", pch = 16, col = 4)
legend("topleft", bty = "n", legend = c("Con rabia", "Totales"), lty = 1, pch = 16,
col = c(4, 3))
Ok… ya no se que pensar…
Ajustemos pues considerando el tamaño de muestra como pesos para el modelo:
m_bn_w <- glm.nb(datos[, 1] ~ año, weights = datos[, 1] + datos[, 2])
summary(m_bn_t)
##
## Call:
## glm.nb(formula = datos[, 1] ~ año, init.theta = 6.994900521,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.569 -0.583 -0.278 0.227 2.339
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -368.9320 44.0678 -8.37 <2e-16 ***
## año 0.1864 0.0221 8.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(6.995) family taken to be 1)
##
## Null deviance: 86.097 on 16 degrees of freedom
## Residual deviance: 17.088 on 15 degrees of freedom
## AIC: 138.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 6.99
## Std. Err.: 3.04
##
## 2 x log-likelihood: -132.88
Pero queda la pregunta ¿Esto que significa?