CLASE 4 - CURSO DE CLIMA Y SALUD
Cargamos las librerias
library(tidyverse)
library(gridExtra)
library(WaveletComp)Cargamos las bases de datos
salud_clima_lamb<-readRDS(file = "data_cruda/salud_clima_lamb.rds")Procedimiento
1) Conceptos basicos de periodicidad y frecuencia
1.1) Definiendo una onda conocida
Vamos a definir 3 ondas seno con amplitud de 1 y periodos de 2, 4 y 5 y una onda suma de estas
max<-60
min<-0
dx<-0.01
onda_t2<-data.frame(t=seq(min,max,dx),valor_t2=sin(seq(min,max,dx)*2*pi/2)) #perido de 2
onda_t4<-data.frame(t=seq(min,max,dx),valor_t4=sin(seq(min,max,dx)*2*pi/4)) #periodo de 4
onda_t5<-data.frame(t=seq(min,max,dx),valor_t5=sin(seq(min,max,dx)*2*pi/5)) #periodo de 5
ondas<-onda_t2%>%
merge(onda_t4,by="t")%>%
merge(onda_t5,by="t")%>%
mutate(onda_suma=valor_t2+valor_t4+valor_t5)
Resultamos teniendo una serie de periĆ³dica que se repite cada 30 unidades de tiempo
1.2) Temperatura como onda
Vamos a trabajar con los datos de la provincia de lambayeque y vamos a eliminar la semana 53 para tener periodos anuales de 52 semanas
lamb_prov<-salud_clima_lamb%>%
filter(provincia=="LAMBAYEQUE",semana!=53)%>%
rename(pp=pp_sem_prov,tmax=tmax_sem_prov,tmin=tmin_sem_prov,
tmean=tmean_sem_prov)%>%
ungroup()
ggplot(data=lamb_prov,aes(x=fecha))+
geom_line(aes(y=tmean),color="steelblue4", lwd=0.8)+
scale_x_date(expand = c(0,0),date_labels = "%W-%Y",
breaks = "12 weeks") +
scale_y_continuous(expand = c(0,0))+
geom_vline(xintercept=as.Date(seq(from = as.Date("2010-01-01"),
to = as.Date("2021-01-01"),by = "1 year")), linetype=4)+
theme_bw()+
theme(axis.text.x = element_text(angle=90,size = 8),
legend.position = "bottom")+
labs(x="Fecha",y="Temperatura (ĀŗC)",
caption = "Fuente: ElaboraciĆ³n propia",
title = "Temperatura semanal de Lambayeque, periodo 2010-2020")
En este caso tenemos a la temperatura, la cual, asemeja una onda periĆ³dica con una periodicidad de un aƱo (52 semanas), a esto lo conocemos como la estacionalidad.
2) Serie de fourier (periodograma)
2.1) Periodograma con nuestras ondas conocidas
Para calcular el periodograma usaremos el comando spectrum, el cual calcula el periodograma haciendo uso de la tranformada rapida de foulier, en este grƔfico vamos a notar a que frecuencias se encuentra distribuidas nuestras funciones.
#' periodograma para onda con t=2
mspect <- spectrum(ondas$valor_t2, log="no", plot=F)
specx <- mspect$freq/dx
specy <- 2*mspect$spec
spectro<-data.frame(x=specx,y=specy)
a1<-ggplot(spectro,aes(x=x,y=y))+
geom_line()+
scale_x_continuous(limits=c(0,2),breaks = seq(0,2,0.1))+
theme_bw()+
labs(title = "Periodograma de onda seno de t=2 /f=0.5",
y="Densidad espectral",x="")
#' periodograma para onda con t=4
mspect <- spectrum(ondas$valor_t4, log="no", plot=F)
specx <- mspect$freq/dx
specy <- 2*mspect$spec
spectro<-data.frame(x=specx,y=specy)
a2<-ggplot(spectro,aes(x=x,y=y))+
geom_line()+
scale_x_continuous(limits=c(0,2),breaks = seq(0,2,0.1))+
theme_bw()+
labs(title = "Periodograma de onda seno de t=4 /f=0.25",
y="Densidad espectral",x="")
#' periodograma para onda con t=5
mspect <- spectrum(ondas$valor_t5, log="no", plot=F)
specx <- mspect$freq/dx
specy <- 2*mspect$spec
spectro<-data.frame(x=specx,y=specy)
a3<-ggplot(spectro,aes(x=x,y=y))+
geom_line()+
scale_x_continuous(limits=c(0,2),breaks = seq(0,2,0.1))+
theme_bw()+
labs(title = "Periodograma de onda seno de t=5 /f=0.2",
y="Densidad espectral",x="")
mspect <- spectrum(ondas$onda_suma, log="no", plot=F)
specx <- mspect$freq/dx
specy <- 2*mspect$spec
spectro<-data.frame(x=specx,y=specy)
a4<-ggplot(spectro,aes(x=x,y=y))+
geom_line()+
scale_x_continuous(limits=c(0,2),breaks = seq(0,2,0.1))+
theme_bw()+
labs(title = "Periodograma de onda suma de t=10 /f=0.5",
y="Densidad espectral",x="Frecuencias")
grid.arrange(a1,a2,a3,a4, nrow = 4)
En el grafico vemos en el eje x las frecuencias y en el eje Y su densidad espectral vemos que las densidades espectrales coinciden con la frecuencia de nuestras ondas, asimismo, vemos las frecuencias de la que fue formada la onda de suma
2.2) Periodograma con datos de temperatura
mspect <- spectrum(lamb_prov$tmean, log="no", plot=F)
dx<-1
specx <- mspect$freq/dx
specy <- 2*mspect$spec
spectro<-data.frame(x=specx,y=specy)
spectro_tmean
ggplot(spectro,aes(x=x,y=y))+
geom_line(na.rm = T)+
scale_x_continuous(limits=c(0,0.2),breaks = seq(0,0.2,0.005))+
geom_vline(xintercept=1/52, linetype=4)+
geom_label(aes(label="Frecuencia de 1/52",x=1/52,y=max(specy)*0.8),hjust="left")+
theme_bw()+
theme(axis.text.x = element_text(angle = 90))+
labs(title = "Periodograma de temperatura media",
y="Densidad espectral",x="Frecuencias")
Vemos que es tiene gran densidad espectral en la frecuencia anual, ademƔs de otras frecuencias espectrales de menor frecuencia.
3) Analisis de Wavelet power spectrum
Con el anƔlisis del periodograma podemos comprender la periodicidad de las ondas, el problema o debilidad de esta tƩcnica es que no sabemos si esta periodicidad cambia en el tiempo, para eso podemos usar el anƔlisis de wavelet, para lo cual vamos a definir una onda seno con amplitud de 1 y periodo variable:
x<-periodic.series(start.period = 5, end.period = 100, phase = 0, length = 1000, make.plot = F)
onda_t_var<-data.frame(t=seq(0,60,length.out = 1000),valor=x)
ggplot(onda_t_var,aes(x=t))+
geom_line(aes(y=valor))+
scale_x_continuous(n.breaks = 60,expand = c(0,0))+
theme_bw()+
labs(x="",y="",title="Onda periodo variable (t<-5,100) (f <-0.1,0.01)")
El periodograma para estas ondas es muy difĆcil de interpretar ya que no nos da informaciĆ³n de en que momento cambia el periodo.
mspect <- spectrum(onda_t_var$valor, log="no", plot=F)
specx <- mspect$freq/1
specy <- 2*mspect$spec
spectro<-data.frame(x=specx,y=specy)
ggplot(spectro,aes(x=x,y=y))+
geom_line(na.rm = T)+
scale_x_continuous(limits=c(0,0.5),breaks = seq(0,0.5,0.02))+
geom_vline(xintercept=0.1, linetype=4)+
geom_vline(xintercept=0.01, linetype=4)+
theme_bw()+
labs(title = "Periodograma de onda de periodo variable (t<-10,100) (f <-0.1,0.01)",
y="Densidad espectral",x="Frecuencias")
Por eso es mejor usar el anĆ”lisis de Wavelet Power Spectrum para ondas en las que el periodo puede que sea variable, este nos da una mejor comprensiĆ³n de la periodicidad de la data.
data <- onda_t_var %>%
select(x=valor)
my.w<-analyze.wavelet(my.data = data,
my.series = "x",
loess.span = 0,
dt = 1,
dj = 1/104,
lowerPeriod = 16,
upperPeriod = 128,
make.pval = TRUE,
n.sim = 10)
wt.image(my.w,
periodlab = "Periodo",
timelab = 'X',
legend.params = list(lab = "Power"),
color.palette = "rainbow(n.levels, start = 0, end = .7)",
main = "Onda periodo variable")
Al hacer el anƔlisis de wavelet power spectrum podemos notar como la periodicidad varia en el tiempo, en este caso va a aumentando
4) Aplicado a datos de Salud y clima
4.1) Preparando la data de salud
Vamos a ver los datos de casos del departamento de Lambayeque, eliminaremos la semana 53 para tener periodos anuales de 52 semanas
analisis<-salud_clima_lamb%>%
group_by(departamento,fecha,semana)%>%
summarise(casos=sum(casos),
pp=mean(pp_sem_prov),
tmax=mean(tmax_sem_prov),
tmin=mean(tmin_sem_prov),
tmean=mean(tmean_sem_prov))%>%
filter(semana!=53)%>%
select(-semana)%>%
ungroup()
ggplot(analisis,aes(x=fecha,y=casos))+
geom_line()+
scale_x_date(expand = c(0,0),date_labels = "%Y",
breaks = "1 year") +
scale_y_continuous(expand = c(0,0))+
geom_vline(xintercept=as.Date(seq(from = as.Date("2010-01-03"),
to = as.Date("2020-01-03"),by = "1 year")), linetype=4)+
theme_bw()+
labs(x="AƱo",y="Casos")
4.2) Periodograma
Con el Periodograma podremos explorar de manera general de que frecuencias espectrales esta conformada nuestra data.
mspect <- spectrum(analisis$casos, log="no", plot=F)
delta <- 1
specx <- mspect$freq/delta
specy <- 2*mspect$spec
spectro<-data.frame(x=specx,y=specy)
ggplot(spectro,aes(x=x,y=y))+
geom_line(na.rm = T)+
scale_x_continuous(limits=c(0,0.2),breaks = seq(0,0.2,0.005))+
geom_vline(xintercept=1/52, linetype=4)+
geom_label(aes(label="Frecuencia de 1/52",x=1/52,y=max(specy)*0.9),hjust="left")+
theme_bw()+
theme(axis.text.x = element_text(angle = 90))+
labs(title = "Periodograma de Casos",
y="Densidad espectral",x="Frecuencias")
En este grafico identificar una gran densidad espectral para la frecuencia alrededor de 0.019, lo cual equivale a un periodo de 52 semanas o un aƱo
4.3) Wavelet power spectrum
Como mencionamos anteriormente, con el anƔlisis de Wavelet Power Spectrum podremos explorar si nuestra frecuencia anual (0.019) encontrada en el periodograma se mantiene durante todo momento.
data <- analisis %>%
select(x= casos)
wavelet_ps<-analyze.wavelet(my.data = analisis,
my.series = "casos",
loess.span = 0,
dt = 1/52,
dj = 1/104,
lowerPeriod = 0.1,
upperPeriod = 3,
make.pval = TRUE,
n.sim = 10)
index.ticks <- seq(1, nrow(analisis), by = 52) #comando para poner aƱos en eje X
index.labels <- ((index.ticks-1)/52)+2010 #comando para poner aƱos en eje X
png(file = "Wavelet_casos.png",width = 900,height = 450)
wps_imagen<-wt.image(wavelet_ps,
periodlab = "Periodo (aƱos)",
timelab = 'Tiempo (semanas)',
legend.params = list(lab = "Power"),
color.palette = "rainbow(n.levels, start = 0, end = .7)",
main = "Wavelet Power Spectrum casos - Dpt de Lambayeque",
spec.time.axis = list(at = index.ticks, labels = index.labels))
dev.off()
Al hacer el anƔlisis de wavelet power spectrum, vemos que el periodo anual que encontramos anteriormente se pierde durante los aƱos 2010 a 2015.
Si hacemos el mismo anĆ”lisis con los datos de precipitaciĆ³n encontramos que se mantiene una frecuencia anual durante casi todo del 2010-2020 (tener en cuenta que no hay datos de pp para el final del 2020)
4.4) Wavelet Coherence
Para este ejemplo veremos si existe una coherencia entre casos de dengue y valores de precipitaciĆ³n, la ventaja de usar este anĆ”lisis es que no es necesario usar un retraso entre alguno de los datos ya que analizamos la data a nivel de frecuencias:
wavelet_cohe<-analyze.coherency(
as.data.frame(analisis),my.pair=c("casos","pp"),
loess.span = 0.75,
lowerPeriod = 0.125,
upperPeriod = 4,
dt = 1/52, dj = 1/104,
n.sim = 50,
window.type.t = 3, window.type.s = 3,
window.size.t = 1/4, window.size.s = 1)
png(file = "Wavelet_coherencia_casos_pp.png",width = 900,height = 450)
wc.image(wavelet_cohe,which.image = "wc",
periodlab = "Periodo (aƱos)",
timelab = 'Tiempo (semanas)',
main = "Coherencia Casos Vs PP - Dpt de Lambayeque",
color.palette = "rainbow(n.levels, start = 0, end = .7)",
legend.params = list(lab = "Power"),
spec.time.axis = list(at = index.ticks, labels = index.labels))
dev.off()
Las Ć”reas con las flechas se grafican en las zonas donde existe significancia entre nuestras variables, segĆŗn este grafico la significancia para el periodo anual entre la PP y casos se pierde durante los aƱos 2013 a 2015 y 2019 a 2020, tambiĆ©n vemos significancias a periodos mĆ”s pequeƱos.
La orientaciĆ³n de las flechas indican si los datos se encuentran en fase o fuera de fase. Si estĆ”n en fase significa que, si los valores de uno aumentan, los del otro tambiĆ©n aumentan, y si los valores de uno disminuyen los del otro tambiĆ©n disminuyen.
En el presente caso, las flechas en su mayorĆa estĆ”n apuntando hacia abajo a la derecha (aproximadamente -pi/4), eso representa que los valores de precipitaciĆ³n y casos estĆ”n en fase y los de precipitaciĆ³n se encuentran adelantados por aproximadamente un octavo de ciclo de sobre los valores de los casos.
Fuente: WaveletComp:A guided tour through the R-package. A. Rosch, H. Schmidbauer. 2014
Para entender los parĆ”metros de las funciones del paquete WaveletComp a continuaciĆ³n presentamos una tabla resumen explicando la utilidad de cada parĆ”metro
