r examples.Rmd

---
title: "Quantitative Methods in Infectious Disease"
subtitle: "Time series examples"
author: "Dan Weinberger"
date: 'November 15, 2023'
runtime: shiny
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(shiny)
source('./R/WaveletPkg.R')
source('./R/plot.wave2.R')
source('./R/biennial.func.R')
library(dplyr)
library(dplyr)
library(ggplot2)
library(MASS)
```

## Data

For these examples, we will use data from the WHO surveillance database for RSV (from the Flu datamart)

Downloads the latest data from WHO

```{r}
#code Courtesy of DeusThindwa; data from WHO
#library(curl)
#library(dplyr)
#library(ggplot2)
#rsv <- read.csv(curl("https://frontdoor-l4uikgap6gz3m.azurefd.net/FLUMART/VIW_FNT?$format=csv"))
#saveRDS(rsv,'./Data/rsv_who.rds')
```

Import and clean the data. Fill in 0s to create a time series. data are not collected year round; assume 0 cases during weeks when there is no data

```{r}
rsv <- readRDS('./Data/rsv_who.rds')

#Clean the data, fill in 0s if not observed
rsvds <-
  rsv %>%
  dplyr::filter(!is.na(RSV)) %>%
  dplyr::select(WHOREGION, FLUSEASON, ORIGIN_SOURCE,HEMISPHERE, COUNTRY_AREA_TERRITORY, MMWR_WEEKSTARTDATE, MMWR_YEAR, MMWR_WEEK, RSV) %>%
  arrange(WHOREGION, COUNTRY_AREA_TERRITORY, MMWR_WEEKSTARTDATE,ORIGIN_SOURCE) %>%
  filter(COUNTRY_AREA_TERRITORY %in% c('France', 'Sweden')) %>%
  mutate(date=as.Date(MMWR_WEEKSTARTDATE)) %>%
  filter(date>='2015-01-01' & date <= '2023-03-01') %>%
  rename(country=COUNTRY_AREA_TERRITORY) %>%
  tidyr::complete(country,ORIGIN_SOURCE,date = seq.Date(from=as.Date('2015-01-04') ,to=as.Date('2022-11-06'), by='week'), fill=list(RSV=0)) %>% #fill time series
  arrange(country,ORIGIN_SOURCE, date ) %>%
  group_by(country, ORIGIN_SOURCE) %>%
  mutate(t= row_number()) %>%  #index for time
  ungroup()
```

Plot the time series

```{r}
rsvds %>% filter(ORIGIN_SOURCE=='NONSENTINEL') %>%
  ggplot( aes(x=date, y=RSV, group=country ,color=country)) +
  geom_line() +
  theme_classic()
```

## Examples

Here are some worked examples of what we are covering in class

## Generating harmonic curves

-First create an index variable for time (1,2,3...t) -Then use sine and cosine functions to generate a wave with specified period -52 week harmonic (sine and cosine waves)

```{r harmonics, echo=TRUE}

rsvds <- rsvds %>%
  mutate( sin52 = sin(2*pi*t/52.143),
          cos52 = cos(2*pi*t/52.143)
  )
```

```{r}
ggplot(rsvds)+
  geom_line( aes(x=date, y=sin52),color='red') +
  geom_line(data=rsvds,aes(x=date, y=cos52),color='blue') +
  theme_classic()
```

## Vary the period

```{r period.slider, echo=FALSE}
#library(shiny)
shinyApp(

  ui = fluidPage(
    sliderInput("period", "Period:",
               min=6, max=208, value=52),
    plotOutput("periodPlot")
  ),

  server = function(input, output) {
    output$periodPlot = renderPlot({
      t=1:520
      sin52=sin(2*pi*t/52.143)
      cos52=cos(2*pi*t/52.143)
      plot(sin52, type='l', col='gray')
      points(sin(2*pi*t/input$period),type='l', col='red', xlab='time' )
      title(paste0('sin(2*pi*t/',input$period,')'))
    },width = "auto", height = "auto")
  }
)
```

## Harmonic Regression

### Generate harmonics with specified periods (Fourier series), fit to data

## Create the needed harmonic variables with 52,26, 17, 104 week periods

focus on Sweden in pre-pandemic period

```{r harm1, echo=TRUE}
mod_ds <- rsvds %>%
  mutate( sin52 = sin(2*pi*t/52.143),
          cos52 = cos(2*pi*t/52.143),
          
          sin26 = sin(2*pi*t*2/52.143),
          cos26 = cos(2*pi*t*2/52.143),
        
          sin17 = sin(2*pi*t*3/52.143),
          cos17 = cos(2*pi*t*3/52.143),
          
          sin104 = sin(2*pi*t*0.5/52.143),
          cos104 = cos(2*pi*t*0.5/52.143),
        ) %>%
  filter(country=='Sweden' &   ORIGIN_SOURCE=='NONSENTINEL' & date<='2020-03-01')


```

## Fit a simple poisson regression with just 52 week period

```{r harm_reg1}
fit1a <- glm.nb(RSV ~ sin52 + cos52, data=mod_ds )
mod_ds$pred1a<- fitted(fit1a)
```

```{r}

ggplot(mod_ds, aes(x=date, y=RSV)) +
  geom_line( color='gray') +
  geom_line(aes(x=date, y=pred1a), color='red') +
  theme_classic()

```

## Why does this work?

This is effectively a simple regression, relating X (sin and cos) to log(Y)(cases)-- but we are plotting both as a function of time. If we plot compared to each other, the relationship is more obvious:

```{r}
ggplot(mod_ds, aes(x=sin52, y=log(RSV+0.5))) +
  geom_point( color='gray') +
  theme_classic()

ggplot(mod_ds, aes(x=cos52, y=log(RSV+0.5))) +
  geom_point( color='gray') +
  theme_classic()

ggplot(mod_ds, aes(x=(2*cos52 + 2*sin52), y=log(RSV+0.5))) +
  geom_point( color='gray') +
  theme_classic()
```

## Add in 2 year periodicity

```{r harm_reg2}
fit2a <- glm.nb(RSV ~ sin52+cos52 + sin104 + cos104,  data=mod_ds )
mod_ds$pred2a<- fitted(fit2a)
```

```{r}
ggplot(mod_ds, aes(x=date, y=RSV)) +
  geom_line( color='gray') +
  geom_line(aes(x=date, y=pred2a), color='red') +
  theme_classic()

```

## Add in 26 week period

```{r harm_reg3a}
fit3a <- glm.nb(RSV ~ sin52+cos52 + sin104 + cos104 + sin26+cos26,  data=mod_ds )
mod_ds$pred3a<- fitted(fit3a)

```

```{r}
ggplot(mod_ds, aes(x=date, y=RSV)) +
  geom_line( color='gray') +
  geom_line(aes(x=date, y=pred3a), color='red') +
  theme_classic()

```

```{r harm_reg3}
fit4a <- glm.nb(RSV ~ sin52+cos52 + sin104 + cos104 + sin26+cos26 + sin17 + cos17 , data=mod_ds )
mod_ds$pred4a<- fitted(fit4a)

```

```{r}
ggplot(mod_ds, aes(x=date, y=RSV)) +
  geom_line( color='black') +
  geom_line(aes(x=date, y=pred2a), color='#1b9e77') +
    geom_line(aes(x=date, y=pred3a), color='#d95f02', lty=2) +
    geom_line(aes(x=date, y=pred4a), color='#7570b3', lty=3) +
    theme_classic()

```

## Determine best model with AIC (smaller=better)

(AIC might not be appropriate here due t residual autocorrelation in the data; should probably have some sort of autoregressive random effect)

```{r aic}
 AIC(fit1a)
 AIC(fit2a)
 AIC(fit3a)
 AIC(fit4a)

```

## Calculate amplitude of the 52 and 104 month periods and their ratio

First for the 12 month period. Amplitude = sqrt(b1\^2 + b2\^2)

```{r amp}
beta_sin52<-coef(fit2a)['sin52']
beta_cos52<-coef(fit2a)['cos52']

amp52<-sqrt(beta_sin52^2 +beta_cos52^2)

amp52 
```

Same for the 24 month period

```{r}
beta_sin104<-coef(fit4a)['sin104']
beta_cos104<-coef(fit4a)['cos104']
amp104<-sqrt(beta_sin104^2 +beta_cos104^2)

amp104 
```

Ratio of 24 to 12 month gives estimate of importance of biennial vs annual cycles

```{r}
ratio104_52<-amp104/amp52
ratio104_52
```

## Fast Fourier Transform

If we don't already know what periods should be present, can we use first use FFT to detect the correct period? Can sqrt or log transform data prior to FFT

```{r fft1, echo=TRUE}
# Do a FFT using the spectrum function 
mspect <- spectrum(sqrt(mod_ds$RSV),  spans=c(2,2), plot=FALSE)
```

This returns estimates of spectral density for each frequency. We will convert to period=1/frequency

```{r}
periodx <- 1/mspect$freq 
specy <- 2*mspect$spec
plot(periodx, specy, xlab="Period (weeks)",
     ylab="Spectral Density", type="l", bty='l')
```

##FFT slider To understand how FFT works, let's simulate monthly time series data with different characteristics. The time series has two periods, and we can change the length of the periods and the amplitude wit these sliders

```{r fft.slider, echo=FALSE}
n=1000
shinyApp(

  ui = fluidPage(
    sidebarLayout(position = "left",
    sidebarPanel(
    sliderInput("period1", "Period1:",
               min=2, max=48, value=12),
     sliderInput("amplitude1", "Amplitude1:",
               min=0.2, max=2.5, value=0.75),
    sliderInput("period2", "Period2:",
               min=2, max=48, value=24),
      sliderInput("amplitude2", "Amplitude2:",
               min=0.2, max=2.5, value=1.5),
      sliderInput("int1", "Intercept:",
               min=0, max=5, value=0)
      ),
    mainPanel(
       plotOutput("periodPlot")
    )
  )
  ),

  server = function(input, output) {
    t=1:120
    xt1a=reactive({input$amplitude1*cos(2*3.14159*t/input$period1)}) #component 1
    xt1b=reactive({input$amplitude2*cos(2*3.14159*t/input$period2)}) #component 2
    intercept<-reactive({input$int1})
    lambda<- reactive ({exp(intercept()+xt1a()+xt1b())})
    xt.slide=reactive ({rpois(length(t),lambda())})
    fft.s=reactive ({spectrum(sqrt(xt.slide()), spans=c(2,2), plot=FALSE)})

    output$periodPlot = renderPlot({
      par(mfrow=c(2,1), mar=c(3,2,1,1))
   plot(t, xt.slide(),type="l", lwd=3, bty="l",ylab="Cases of Disease", xlab="Month(n)", main="Cases")

      periodx <- 1/fft.s()$freq 
      specy <- 2*fft.s()$spec
plot(periodx, specy, xlab="Period (months)",
     ylab="Spectral Density", type="l", bty='l')
      

    },width = "auto", height = "auto")
  })
```

## Wavelets

### RSV data from WHO surveillance. Compare Sweden and France

##Run wavelets on the sqrt(Cases) for France and Sweden

Note: these Wavelet functions originally developed by Johannson et al PLOS Medicine <https://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.1000168>

```{r wave2, echo=TRUE}
wave1<- rsvds %>% 
    filter(country=='Sweden' & ORIGIN_SOURCE=='NONSENTINEL' ) %>% 
    mutate(x=sqrt(RSV)) %>%
   dplyr::select(x, t)
  
cwt1 <- complete.cwt(wave1, dj=1/26, dt=1/52)

plot.cwt(cwt1)

plot.wave2(cwt1 )

```

```{r}
	wave2<- rsvds %>% 
    filter(country=='France' & ORIGIN_SOURCE=='NONSENTINEL' ) %>% 
    mutate(x=sqrt(RSV)) %>%
  dplyr::select(x, t)	

cwt2 <- complete.cwt(wave2, dj=1/26, dt=1/52)

plot.cwt(cwt2)
plot.wave2(cwt2 )

```

##Compare phase angle of the 2 series we generated

For the 52 week harmonic

```{r wave3, echo=FALSE}

	phase1<-phase(cwt1, s=c(0.8, 1.2))
	phase2<-phase(cwt2, s=c(0.8, 1.2))
	
	dates <- unique(rsvds$date)
	plot(dates, phase1,type='l', col="blue", lty=1) #Sweden
	points(dates, phase2,type='l', col="red", lty=1) #France
```

##Calculate and view the phase difference

```{r wave4, echo=TRUE}
#Take modulo (%%) of phase difference to account for edge effects as described by Grenfell. This corrects for problems at start and begining of cycles
	phase.diff2.1= phase2-phase1

	phasediff2.1_correct <- ((phase.diff2.1 + 3*pi) %% (2*pi)) - (pi)
	#shows corrected phase diff in modulo for the 12 month phase

	plot(dates,phasediff2.1_correct,type='l',col="blue", lty=1,ylim=c(-pi,pi))

```

##Convert phase difference to weeks

```{r wave5, echo=FALSE}
	phase.diff.wk <- phasediff2.1_correct * 52.143/(2*pi)
	plot(dates,phase.diff.wk,type='l',col="blue", lty=1,ylim=c(-26,26))
```

##Reconstruct series, filtering out noise

```{r wave.recon, echo=FALSE}

y = wave1$x #Original data

 recon1 <- reconstruct.series(cwt1)  #square to get back to original scale
 recon1a <- reconstruct.series(cwt1, sp1=0.5, sp2=100) 
 recon1b <- reconstruct.series(cwt1, sp1=0.8, sp2=100)
 
 plot(dates,recon1,type='l',col="#1b9e77", lty=1)
  points(dates,recon1a,type="l",col="#d95f02",lty=2)
 points(dates,recon1b,type="l",col="#7570b3",lty=3)

```

## Compare biennial vs annual strength in the 2 series

Confirms that biennial pattern is stronger in Sweden than France

```{r annual.biennial}

series1.biennial.ratio <- biennial.func(ds=cwt1,period.range.short=c(0.8,1.2),period.range.long=c(1.8,2.2))

series2.biennial.ratio <- biennial.func(ds=cwt2,period.range.short=c(0.8,1.2),period.range.long=c(1.8,2.2))

series1.biennial.ratio # Sweden

series2.biennial.ratio # France
```

## Coherence analysis Evaluate coherence between the 2 series. At which periodicity do we have correlation?

periods and which time points are the series in sync? should actually use more monte carlo draws, but takes too long for in class demo (note if set dt=1/52, gives same result but takes much longer)

```{r coherence}
#Coherence analysis
coh1.2 <- complete.coh(wave1,wave2, dj=1/26,mc.sim=25) #mc.sim=1000
plot.coh(coh1.2)
```