Analysis.Rmd

---
title: "Acidentes_ofidicos"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r}
#Load specific function created for this 
source("HelpFunctions.R")
```


```{r}
#read data
dados <- read_csv("data/2019_06_25_Snakebite.csv")
```


```{r}
#Genral characteristics of accidents
desc_vars <- c("CS_SEXO", "ANT_ZONA", "EVOLUCAO", "ANT_LOCA_1", "ANI_SERPEN", "ANT_TEMPO_", "TRA_CLASSI")
for(desc_var in desc_vars){
    print(knitr::kable(compute_prop(dados, !! rlang::sym(desc_var))))
}

#compute descriptive analysis in each city
compute_desc(dados, AGE, CITY, "median") %>% spread(key = "CITY", value = "Median")
for(desc_var in desc_vars){
 print("Proportion values")
 print(knitr::kable(compute_desc(dados, !! rlang::sym(desc_var), CITY) %>% 
    select(-numb) %>% 
    spread(key = "CITY", value = "Proportion")))
 
  print("N values")
  print(knitr::kable(compute_desc(dados, !! rlang::sym(desc_var), CITY) %>% 
    select(-Proportion) %>% 
    spread(key = "CITY", value = "numb")))
}
```


```{r}
#Spliting data for each city
pvh <-  split_city(data = dados, city = "Porto Velho")
cacoal <- split_city(data = dados, city = "Cacoal")
ariquemes <- split_city(data = dados, city = "Ariquemes")
vilhena <- split_city(data = dados, city = "Vilhena")
```


```{r}
#read climate datasets
max_temp <- make_df(path = "data/2019_06_24_LSTD.csv") %>% rename(max_temp = var_sum)
min_temp <- make_df(path = "data/2019_06_24_LSTN.csv")  %>% rename(min_temp = var_sum)
rainfall <- make_df(path = "data/2019_06_24_Rainfall.csv")  %>% rename(rain = var_sum)
humidity <- make_df(path = "data/2019_06_24_Humidity.csv") %>% rename(hum = var_sum)
humidity <- humidity %>% 
   mutate(hum = hum * 1000)

#Join all datasets
clima <- list(max_temp, min_temp, rainfall, humidity)
clima_df <- clima %>% reduce(left_join, by = c("ym", "city"))
clima_df <- clima_df %>% 
      mutate(amplitude = max_temp - min_temp)


climate_vars <- c("rain", "hum", "max_temp", "min_temp")

for(clima_var in climate_vars){
  print(clima_var)
  print(knitr::kable(clima_df %>% 
  group_by(city) %>% 
  summarise(IQ = quantile(!! rlang::sym(clima_var), probs = c(0.25), na.rm = TRUE),
            mean = mean(!! rlang::sym(clima_var), na.rm = TRUE),
            median = median(!! rlang::sym(clima_var), na.rm = TRUE),
            IIQ = quantile(!! rlang::sym(clima_var), probs = c(0.75), na.rm = TRUE),
            sd = sd(!! rlang::sym(clima_var), na.rm = TRUE))))
  
}
```


```{r}
#read population data
df_pop <-  read_csv("data/2019_pop.csv")
df_pop <- df_pop %>% 
  filter(cod_ibg %in% list(110020, 110004, 110002, 110030))
df_pop <- df_pop %>% 
  mutate(city = as.character(factor(cod_ibg, levels = c(110020, 110002, 110004, 110030),
                                    labels = c("Porto Velho", "Ariquemes", "Cacoal", "Vilhena"))))

clima_df <- clima_df %>% 
  mutate(ano = as.numeric(substr(ym, start = 1, stop = 4)))

clima_df <- left_join(clima_df,df_pop)


#Remove data from memory
rm(list = c("clima", "min_temp", "max_temp", "rainfall", "df_pop"))
```


```{r}
#Plot incidence by year
pvh <- left_join(pvh, clima_df, by = c("ym", "city"))
cacoal <- left_join(cacoal, clima_df, by = c("ym", "city"))
ariquemes <- left_join(ariquemes, clima_df, by = c("ym", "city"))
vilhena <- left_join(vilhena, clima_df, by = c("ym", "city"))
```


```{r}
#Descriptive analysis
desc_df <-  rbind(pvh, cacoal, ariquemes, vilhena)
desc_df <- desc_df %>%
  mutate(incidence = total/pop * 100000) %>%
  rename(time = data)

#compute overall incidence among cities
desc_df %>% 
  filter(city == "Vilhena") %>% 
  group_by(ano) %>% 
  summarise(cases = sum(total),
            av = sum(incidence)) %>% 
  mutate(N = sum(cases))
    

desc_df %>%
  mutate(city = factor(city, levels = c("Porto Velho", "Ariquemes", "Cacoal", "Vilhena"))) %>% 
  ggplot(aes(x = time, y = incidence)) +
          geom_line() +
          theme_bw(base_size = 28) +
          labs(x = NULL, y = "Snakebite incidence per 100,000 persons") +
          facet_wrap(~city) +
          ggsave("results/Fig_1.png", device = "png", width = 18, height = 12)


cities <- c("Cacoal", "Ariquemes", "Porto Velho", "Vilhena")

for(c in cities){
  print(paste("Yearly Descriptive analsys for", c))
  print(desc_df  %>% 
    filter(city == c) %>% 
    group_by(ano, pop) %>% 
    summarise(cases = sum(total)) %>%
    mutate(incidence = cases/pop * 100000) %>% 
    arrange(desc(incidence)))
}



#Compute overall monthly incidence
desc_df %>% 
  group_by(city) %>% 
  summarise(av = mean(incidence),
            LB = mean(incidence) - (1.96 * sd(incidence)/sqrt(144)),
            UB = mean(incidence) + (1.96 * sd(incidence)/sqrt(144)))

for(c in cities){
  print(desc_df %>% 
    filter(city == c) %>% 
    mutate(month = month(time)) %>% 
    group_by(month) %>% 
    summarise(average = mean(incidence),
              sd = sd(incidence)) %>% 
    arrange(desc(average))) 
}
```


```{r}
#Time series plot
clima_df <- clima_df %>% 
  mutate(time = ymd(paste0(ym, "-01")))

plot_seas_ts(clima_df, rain, expression(paste("Average rainfall in ", mm^3))) +
  ggsave("results/Supp_rain.png", device = "png", width = 18, height = 12)

plot_seas_ts(clima_df, max_temp, expression(paste("Average maximum temperature", ~degree, "C"))) +
  ggsave("results/Supp_maxTemp.png", device = "png", width = 18, height = 12)

plot_seas_ts(clima_df, min_temp, expression(paste("Average minimum temperature", ~degree, "C"))) +
  ggsave("results/Supp_minTemp.png", device = "png", width = 18, height = 12)

plot_seas_ts(clima_df, hum, expression(paste("Average humidity ", "kg/kg"))) +
  ggsave("results/Supp_hum.png", device = "png", width = 18, height = 12)

plot_seas_ts(desc_df, incidence, "Seasonality of snakebite incidence per 100,000 persons") +
  ggsave("results/Figure_1.png", device = "png", width = 18, height = 12)
```


```{r}
#Plot monthly ts with covariates
ggplot(desc_df) +
  geom_line(aes(x = time, y = total/pop * 100000)) +
  geom_line(aes(x = time, y = max_temp/2), color = "blue") +
  geom_line(aes(x = time, y = min_temp/2), color = "purple") +
  scale_y_continuous(sec.axis = sec_axis(~.*2, 
                                         name = expression(paste("Average temperature", ~degree, "C")))) +
  facet_wrap(~city) +
  labs(y = NULL, x = NULL) + 
  theme_bw(base_size = 28) +
  ggsave("results/Figure_1_Temp.tiff", device = "tiff", width = 20, height = 8)

ggplot(desc_df) +
  geom_line(aes(x = time, y = total/pop * 100000)) +
  geom_line(aes(x = time, y = hum -3), color = "grey") +
  scale_y_continuous(
                     sec.axis = sec_axis(~.+ 3,
                        name = expression(paste("Average humidity ", "kg/kg")))) +
  labs(y = NULL, x = NULL) +
  facet_wrap(~city) +
  theme_bw(base_size = 28) +
  ggsave("results/Figure_1_Hum.tiff", device = "tiff", width = 20, height = 8)

ggplot(desc_df) +
  geom_line(aes(x = time, y = total/pop * 100000)) +
  geom_line(aes(x = time, y = rain/27), color = "steelblue") +
  scale_y_continuous("Snakebite incidence per 100,000 inhab.",
                     sec.axis = sec_axis(~.*27,
                        name = expression(paste("Average rainfall in ", mm^3)))) +
  labs(y = NULL, x = NULL) +
  facet_wrap(~city) +
  theme_bw(base_size = 28) +
  ggsave("results/Figure_1_Rain.tiff", device = "tiff", width = 20, height = 8)
```


```{r}
#Run models
model_ariquemes <- run_inla(ariquemes, clima_df, TRUE)
model_cacoal <- run_inla(cacoal, clima_df, TRUE)
model_pvh <-  run_inla(pvh, clima_df, TRUE)
model_vilhena <-  run_inla(vilhena, clima_df, TRUE)
```



```{r}
summary(model_pvh[[3]])
summary(model_ariquemes[[3]])
summary(model_cacoal[[3]])
summary(model_vilhena[[3]])
```

```{r}
ariquemes <- ariquemes %>% 
  mutate(fitted = model_ariquemes[[3]]$summary.fitted.values$mean * expected(pop,total,1),
         lb = model_ariquemes[[3]]$summary.fitted.values$`0.025quant` * expected(pop,total,1),
         ub = model_ariquemes[[3]]$summary.fitted.values$`0.975quant` * expected(pop,total,1),
         residuals = total - fitted)

cacoal <- cacoal %>% 
  mutate(fitted = model_cacoal[[3]]$summary.fitted.values$mean * expected(pop,total,1),
         lb = model_cacoal[[3]]$summary.fitted.values$`0.025quant` * expected(pop,total,1),
         ub = model_cacoal[[3]]$summary.fitted.values$`0.975quant` * expected(pop,total,1),
         residuals = total - fitted)

pvh <- pvh %>% 
  mutate(fitted = model_pvh[[3]]$summary.fitted.values$mean * expected(pop,total,1),
         lb = model_pvh[[3]]$summary.fitted.values$`0.025quant` * expected(pop,total,1),
         ub = model_pvh[[3]]$summary.fitted.values$`0.975quant` * expected(pop,total,1),
         residuals = total - fitted)

vilhena <- vilhena %>% 
  mutate(fitted = model_vilhena[[3]]$summary.fitted.values$mean * expected(pop,total,1),
         lb = model_vilhena[[3]]$summary.fitted.values$`0.025quant` * expected(pop,total,1),
         ub = model_vilhena[[3]]$summary.fitted.values$`0.975quant` * expected(pop,total,1),
         residuals = total - fitted)

dfs <- list(ariquemes, cacoal, pvh, vilhena)
check_df <- purrr::reduce(dfs, function(x,y){rbind(x,y)})
```


```{r}
ggplot(check_df) +
  geom_ribbon(aes(x = data, ymin = lb, ymax = ub), fill = "grey70" ) +
  geom_line(aes(x = data, y = total), color = "black") +
  geom_line(aes(x = data, y = fitted), color = "orange") +
  theme_bw(base_size = 28) +
  facet_wrap(~city, scales = "free_y") +
  labs(x = NULL, y = "Cases") +
  ggsave("results/Figure_S3.png", device = "png", width = 18, height = 12)


ggplot(check_df) +
  geom_point(aes(x = time_str, y = residuals), color = "black") +
  theme_bw(base_size = 28) +
  facet_wrap(~city, scales = "free_y") +
  labs(x = NULL, y = "Cases") +
  ggsave("results/Figure_S4.png", device = "png", width = 18, height = 12)
```


```{r}
ggarrange(model_pvh[[1]],
          model_ariquemes[[1]],
          model_cacoal[[1]], 
          model_vilhena[[1]],
          labels = c("(A)", "(B)", "(C)", "(D)"),
          font.label = list(size = 20)) +
  ggsave("results/Figure_2.png", device = "png", width = 18, height = 12)
```