02-sensitivity_analyses.Rmd

---
title: "Sensitivity analyses"
author: "Paloma Cárcamo"
output: html_document
---

## Load packages
```{r}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, broom, spdep, stats)
```

## Load data
```{r}
data(Peru, package = "innovar")

loreto <- Peru |> 
  filter(dep == "LORETO")

data <- read_rds("data/dengue-malaria.rds")
```

## Function for TLCC
```{r}
cross_corr <- function(df, lag = 104) {
    tidy(ccf(x = df$cases_d_var,
             y = df$cases_m_var,
             lag.max = lag,
             plot = FALSE))
}
```

## Create neighors list and calculate spatial weights
```{r}
nb <- poly2nb(loreto, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
```

## Run TLCCs with random permutations of date
```{r, eval = FALSE}
# Set number of simulations to run
n = 999

# Create empty df for results of loop
sens_db <- vector("list", length = n)

for (i in 1:n) {
  set.seed(n)
  # Reshuffle week start
  by_distr <- data |> 
    group_by(distr) |> 
    mutate(week_start = sample(week_start)) |> 
    arrange(distr, week_start) |> 
    ungroup() |> 
    select(distr, week_start, 
           cases_m, cases_d, 
           cases_m_var, cases_d_var,
           cases_m_var_abs, cases_d_var_abs) |> 
    group_by(distr) |> 
    nest()
  
  # Calculate TLCC coefficients
  by_distr2 <- by_distr |> 
    mutate(crosscorr = purrr::map(data, cross_corr)) 
  
  ccfs <- unnest(by_distr2, crosscorr)
  
  # Extract max ccfs
  max_ccfs <- ccfs |> 
    select(distr, lag, acf)|> 
    group_by(distr) |> 
    mutate(lag = if_else(is.na(acf), NA, lag)) |> 
    slice_max(order_by = abs(acf)) |> 
    mutate(acf = round(acf,2)) |> 
    unique()
  
  # Join max ccfs to map
  map <- loreto |> 
    left_join(max_ccfs, by = "distr") |> 
    replace_na(list(lag = 0, acf = 0))
  
  # Run Local Moran's
  loc_mor <- localmoran(map$acf, lw)
  loc_mor_full <- cbind(max_ccfs, loc_mor) |> 
    mutate(stat = if_else(`Pr(z != E(Ii))` < 0.05, Ii, NA),
           run = i) |> 
    select(run, distr, lag, acf, stat, stat2 = `Pr(z != E(Ii))`)
  
  sens_db[[i]] <- loc_mor_full
  
}

sens_db_full <- do.call(rbind, sens_db)

# write_rds(sens_db_full, "data/sens_results_full.rds")
```

## Rerun original TLCC
```{r}
# Reference: TLCC with original data
by_distr <- data |> 
    ungroup() |> 
    select(distr, week_start, 
           cases_m, cases_d, 
           cases_m_var, cases_d_var,
           cases_m_var_abs, cases_d_var_abs) |> 
    group_by(distr) |> 
    nest()
  
# Calculate TLCC coefficients
by_distr2 <- by_distr |> 
  mutate(crosscorr = purrr::map(data, cross_corr)) 

ccfs <- unnest(by_distr2, crosscorr)

# Extract max ccfs
max_ccfs <- ccfs |> 
  select(distr, lag, acf)|> 
  group_by(distr) |> 
  mutate(lag = if_else(is.na(acf), NA, lag)) |> 
  slice_max(order_by = abs(acf)) |> 
  mutate(acf = round(acf,2)) |> 
  unique()

# Join max ccfs to map
map <- loreto |> 
  left_join(max_ccfs, by = "distr") |> 
  replace_na(list(lag = 0, acf = 0))

# Run Local Moran's
loc_mor <- localmoran(map$acf, lw)

reference <- cbind(max_ccfs, loc_mor) |> 
  mutate(stat = if_else(`Pr(z != E(Ii))` < 0.05, Ii, NA)) |> 
  select(distr, lag, acf, stat, stat2 = `Pr(z != E(Ii))`)
```

## Compare distributions of lags, TLCC coefficients and Moran's test statistics in original vs. randomly reshuffled samples
```{r}
sens_db_full <- read_rds("data/sens_results_full.rds")

# Lags
hist(reference$lag, breaks = 12)
hist(sens_db_full$lag, breaks = 12)

ks.test(reference$lag, sens_db_full$lag)

# TLCC coefficients
hist(reference$acf, breaks = 12, freq = FALSE)
hist(sens_db_full$acf, breaks = 12, freq = FALSE)

ks.test(reference$acf, sens_db_full$acf)

# Moran's test statistics (only significant stats)
# hist(reference$stat, breaks = 12, freq = FALSE)
# hist(sens_db_full$stat, breaks = 12, freq = FALSE)
# 
# ks.test(reference$stat, sens_db_full$stat)

# Moran's test statistics (all stats)
hist(reference$stat2, breaks = 12, freq = FALSE)
hist(sens_db_full$stat2, breaks = 12, freq = FALSE)

ks.test(reference$stat2, sens_db_full$stat2)
```

## Run TLCCs with random numbers of cases
```{r, eval = FALSE}
# Set number of simulations to run
n = 999

# Create empty df for results of loop
sens_db2 <- vector("list", length = n)

# run nb regression for each district, extract parameters and save in summary dataset
# (remove Rosa Panduro, has no dengue cases)
distr_stats <- data |>
  filter(distr != "ROSA PANDURO") |> 
  group_by(distr) |> 
  do({
    model_m <- MASS::glm.nb(cases_m ~ 1, data = .)
    parameters_m <- coefficients(summary(model_m))
    mean_cases_m <- exp(parameters_m[1, 1])
    size_param_m <- parameters_m[1, 2]^(-1)
    model_d <- MASS::glm.nb(cases_d ~ 1, data = .)
    parameters_d <- coefficients(summary(model_d))
    mean_cases_d <- exp(parameters_d[1, 1])
    size_param_d <- parameters_d[1, 2]^(-1)
    data.frame(distr = unique(.$distr), 
               mean_cases_m = mean_cases_m,
               size_param_m = size_param_m,
               mean_cases_d = mean_cases_d,
               size_param_d = size_param_d)
  })

# create new nb list and spatial weights without Rosa Panduro
nb <- poly2nb(loreto[loreto$distr != "ROSA PANDURO",], queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

for (i in 1:n) {
  set.seed(n)
  # Assign random case numbers
  by_distr <- data |>
    filter(distr != "ROSA PANDURO") |> 
    ungroup() |> 
    left_join(distr_stats, by = "distr") |> 
    mutate(cases_m = rnbinom(n(), size = size_param_m, mu = mean_cases_m),
           cases_d = rnbinom(n(), size = size_param_d, mu = mean_cases_d)) |> 
    mutate(cases_m_var = (cases_m/lag(cases_m) - 1) * 100,
           cases_d_var = (cases_d/lag(cases_d) - 1) * 100) |> 
    mutate(cases_m_var = if_else(is.finite(cases_m_var), cases_m_var, 0),
           cases_d_var = if_else(is.finite(cases_d_var), cases_d_var, 0),
           cases_m_var = if_else(is.nan(cases_m_var), 0, cases_m_var),
           cases_d_var = if_else(is.nan(cases_d_var), 0, cases_d_var)) |> 
    mutate(cases_m_var_abs = (cases_m - lag(cases_m)),
           cases_d_var_abs = (cases_d - lag(cases_d))) |> 
    mutate(cases_m_var_abs = if_else(is.finite(cases_m_var_abs), cases_m_var_abs, 0),
           cases_d_var_abs = if_else(is.finite(cases_d_var_abs), cases_d_var_abs, 0),
           cases_m_var_abs = if_else(is.nan(cases_m_var_abs), 0, cases_m_var_abs),
           cases_d_var_abs = if_else(is.nan(cases_d_var_abs), 0, cases_d_var_abs)) |> 
    filter(!is.na(cases_m_var)) |> 
    filter(!is.na(cases_d_var)) |> 
    filter(!is.na(cases_m_var_abs)) |> 
    filter(!is.na(cases_d_var_abs)) |> 
    select(distr, week_start, 
           cases_m, cases_d, 
           cases_m_var, cases_d_var,
           cases_m_var_abs, cases_d_var_abs) |> 
    group_by(distr) |> 
    nest()
  
  # Calculate TLCC coefficients
  by_distr2 <- by_distr |> 
    mutate(crosscorr = purrr::map(data, cross_corr)) 
  
  ccfs <- unnest(by_distr2, crosscorr)
  
  # Extract max ccfs
  max_ccfs <- ccfs |> 
    select(distr, lag, acf)|> 
    group_by(distr) |> 
    mutate(lag = if_else(is.na(acf), NA, lag)) |> 
    slice_max(order_by = abs(acf)) |> 
    mutate(acf = round(acf,2)) |> 
    unique()
  
  # Join max ccfs to map
  map <- loreto |> 
    left_join(max_ccfs, by = "distr") |> 
    filter(distr != "ROSA PANDURO") |> 
    replace_na(list(lag = 0, acf = 0))
  
  # Run Local Moran's
  loc_mor <- localmoran(map$acf, lw)
  loc_mor_full <- cbind(max_ccfs, loc_mor) |> 
    mutate(stat = if_else(`Pr(z != E(Ii))` < 0.05, Ii, NA),
           run = i) |> 
    select(run, distr, lag, acf, stat, stat2 = `Pr(z != E(Ii))`)
  
  sens_db2[[i]] <- loc_mor_full
  
}

sens_db2_full <- do.call(rbind, sens_db2)

# write_rds(sens_db2_full, "data/sens_results_full_2.rds")
```

## Compare distributions of lags, TLCC coefficients and Moran's test statistics in original vs. randomly reshuffled samples
```{r}
sens_db2_full <- read_rds("data/sens_results_full_2.rds")

# Lags
hist(reference$lag, breaks = 12)
hist(sens_db2_full$lag, breaks = 12)

ks.test(reference$lag, sens_db2_full$lag)

# TLCC coefficients
hist(reference$acf, breaks = 12, freq = FALSE)
hist(sens_db2_full$acf, breaks = 12, freq = FALSE)

ks.test(reference$acf, sens_db2_full$acf)

# Moran's test statistics (only significant stats)
# hist(reference$stat, breaks = 12, freq = FALSE)
# hist(sens_db2_full$stat, breaks = 12, freq = FALSE)
# 
# ks.test(reference$stat, sens_db2_full$stat)

# Moran's test statistics (all stats)
hist(reference$stat2, breaks = 12, freq = FALSE)
hist(sens_db2_full$stat2, breaks = 12, freq = FALSE)

ks.test(reference$stat2, sens_db2_full$stat2)
```