CO y Mortalidad. Análisis con R usando datos cargados de Internet

CO y Mortalidad. Análisis con R usando datos cargados de Internet

En Octubre del 2009, en las V Jornadas Latinoamericanas de Química y Física Ambiental, presenté un trabajo que compara el efecto del monóxido de carbono ambiental sobre la mortalidad de los habitantes de Santiago, Chile y Buenos Aires, Argentina. Para este trabajo se tomaron datos de mortalidad total, masculina, femenina, de menores de 14 años, de 14 a 65 y mayores a 65 años de las ciudades de Santiago y Buenos Aires. Para que fuesen comparables utilicé la misma metodología que consiste en un Modelo Semi-Paramétricode regresión quasi-poisson de la mortalidad esperada sobre los días de la semana, mes, monóxido de carbono, temperatura, humedad relativa y tiempo (variables estadísticamente significativas). Las variables de confusión fueron controladas usando los suavizadores cubic splines. Se consideró sólamente el CO pues esos eran los datos disponibles de monitoreo en Buenos Aires.

En Marzo del 2011 se realizó el curso R-Básico a cargo de AESPRO. Durante ese curso presenté un breve resumen de mi presentación  del CO, y en ella mostré como cargar datos desde Internet y realizar algunas estadísticas básicas usando ggplot 2.

Ahora repetiré las instrucciones:

Primero Abrir R en Windows o escribir R en una consola de Linux, yo uso Ubuntu-Linux.

# Establecer el working directory

setwd("/home/sergio/r") #Mi carpeta "r" dentro de mi home "sergio"

datos <- read.csv("http://aire.cenma.cl/sibarra/datos/saba.csv", header=T)

attach(datos)

names(datos)

#tenemos todas las variables para operar con ellas

range(dia) #vemos que contamos con 665 días de observación

#luego cargamos la librería ggplot2 para ver histogramas a color

library(ggplot2)

png("histograma.png", width = 480, height = 680,units="px")

#Con esto creamos el archivo "histograma.png" que contendrá la imagen

#que se especifica a continuación. Este archivo se encuentra en el Working Directory

m <- ggplot(datos, aes(x=m_s))

m + geom_histogram(aes(fill = ..count..)) +

scale_fill_gradient("Count", low = "yellow", high = "red") +

opts(title = "Histograma de la mortalidad en Santiago")

dev.off()

#la imagen está posteada al comienzo

#Se aprecia ligeramente una distribución no normal, de contar con mas datos 

#de apreciaría mejor

mean(m_s)

var(m_s)

#se aprecia que la varianza es mayor que la media, otro indicador que no es normal,

# y además es poisson

#Para hacer estas pruebas uso el software STATA, resultados que aquí no incluiré

#Ahora, algo de regresiones

library(mgcv)

m1 <- gam(m_s~lu+ma+mi+ju+vi+sa+do, family=quasipoisson)

summary(m1) #se ve que el día lunes es estadísticamente significativo

m2 <- gam(m_male_s~lu+ma+mi+ju+vi+sa+do, family=quasipoisson)

summary(m2)

#se ve que el día lunes, sábado y domingo son estadísticamente significativos

m3 <- gam(m_female_s~lu+ma+mi+ju+vi+sa+do, family=quasipoisson)

summary(m3)

#No se aprecian efectos del día estadísticamente significativo
s