02-time-vaying-coef-models.Rmd

---
title: "Time-varying coefficients model using GAMs"
output:
  html_document:
    toc: yes
  word_document:
    toc: yes
editor_options:
  chunk_output_type: inline
---

```{r setup, echo=FALSE}
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
```

We are going to use *generalized additive models* (GAMs) to estimate time-varying coefficients for climate variables in order to quantify their effects on malaria incidence by specie thought the time period of the study.

Lets load the packages we are going to use for data reading, wrangling and analysis.

```{r message=FALSE, warning=FALSE}
# Utilities
library(htmltools)
library(tictoc)
library(fs)

# Data reading
library(readr)
library(skimr)

# Data manipulation
library(tibble)
library(lubridate)
library(dplyr)
library(scales)
library(DT)
library(stringr)
library(tidyr)

# Plots
library(ggplot2)
library(ggsci)
library(GGally)
library(DHARMa)
library(grid)
library(gridExtra)
library(graphics)

# GAM
library(mgcv)
library(gratia)
```

Setting up the theme for our plots.

```{r ggplot-setup}
ggplot2::theme_set(theme_bw())
ggplot2::theme_update(panel.grid = element_blank())
```

# Data reading

The cleaned up data can be found on the `data/processed/` directory on this repository. We are going to define the file path using the function `path()` of the `fs` packages. 

```{r}
data_path <- path("data", "processed")
file_name <- "pro_malaria-monthly-cases-district-loreto.csv"
file_path <- path(data_path, file_name)
file_path
```

Create an object with the column names.

```{r message=FALSE, warning=FALSE}
col_names <- c(
    "district", "year", "month", "falciparum", "vivax", "aet", "prcp", "q", 
    "soilm", "tmax", "tmin", "pop2015", "province", "region", "dttm", "time"
)
```

Create an object with the column types. We exclude `soilm` and `region` from the analysis.

```{r}
col_types <-
  cols_only(
    district = readr::col_factor(), year = readr::col_factor(), 
    month = readr::col_factor(),
    falciparum = col_integer(), vivax = col_integer(), aet = col_double(),
    prcp = col_double(), q = col_double(), tmax = col_double(),
    tmin = col_double(), pop2015 = col_integer(), province = readr::col_factor(),
    dttm = col_datetime(), time = col_integer()
  )
```

Reading data file into a data frame.

```{r message=FALSE, warning=FALSE}
monthly_cases <- 
  read_csv(
    file = file_path, col_names = col_names, col_types = col_types, skip = 1,
    locale = readr::locale(encoding = "UTF-8")
  )
```

Display the data.

```{r}
monthly_cases
```

# Data wrangling

The following wrangling pipe consists on these steps:

- Re-scale the numeric predictors (`aet`, `prcp`, `q`, `tmax` and `tmin`) from 0 to 1.
- Create a column with the population values on the natural logarithmic scale.
- For each district, create columns for the 1, 3, 6 and 12 months-lagged values of the predictors.
- Remove `NA`s from the lagged values.

```{r}
monthly_cases %>% 
  mutate(
    across(c(aet, prcp, q, tmax, tmin), scales::rescale, to = c(0, 1)),
    pop2015_ln = log(pop2015)
  ) %>% 
  group_by(district) %>% 
  arrange(dttm) %>% 
  mutate(
    across(
      aet:tmin, 
      list(
        lag1 = ~ lag(.x, n = 1L), lag3 = ~ lag(.x, n = 3L), 
        lag6 = ~ lag(.x, n = 6L), lag12 = ~ lag(.x, n = 12L)
      ), 
      .names = "{.col}_{.fn}"
    )
  ) %>% 
  na.omit() %>% 
  ungroup() -> 
  monthly_cases_analysis
```

Display the final data for analysis.

```{r}
monthly_cases_analysis
```

# Time-varying coefficients model building

## P. vivax

We start with the GAM with a time-varying effects configuration for P. vivax cases. The basis dimension `k` for the `s(district)` (random effect) term was set to 49 because there are 49 districts in the data. The remaining basis dimensions were selected by checking the test for adequacy of basis dimensions produced by the function `gam.check` iteratively until the simulated p-values were big enough to have a strong likelihood that they were adequate. 

```{r message=FALSE}
tic("GAM fitting")
vivax_tvcm <- bam(
  vivax ~ offset(pop2015_ln) + s(district, bs = "re", k = 49) +
    s(time, by = aet,  bs = "tp", k = 50) +
    s(time, by = prcp, bs = "tp", k = 30) +
    s(time, by = tmin, bs = "tp", k = 100),
  family = nb(), data = monthly_cases_analysis, method = "fREML", discrete = TRUE
)
toc()

cat("\nConverged?", vivax_tvcm$converged, "\n")

summary(vivax_tvcm)
```

Now we check the model using the default function in the `mgcv` package for model diagnosis.

```{r include=TRUE}
par(mfrow = c(2, 2))
gam.check(vivax_tvcm)
```

We can use the `gratia` package to produce tidier outputs for GAM models using `mgcv`. Here, we use the  `appraise` function to show diagnostics plots.

```{r fig.width=12, fig.height=8}
appraise(vivax_tvcm, point_col = "steelblue", point_alpha = 0.4)
```

We plot the fitted smooth functions for the random effect term and the covariates.

```{r fig.width=12, fig.height=8}
draw(vivax_tvcm, residuals = TRUE)
```

Finlly, we check for concurvity. 

```{r}
model_concurvity(vivax_tvcm) %>% 
  filter(type == "estimate")
```

In the Supplementary Information section of the article there is a further discussion why the concurvity is high in our models.

### P. falciparum

Now, we repeat the GAM fitting and diagnosis for P. falciparum.

```{r message=FALSE}
tic("GAM fitting")
falciparum_tvcm <- bam(
  falciparum ~ offset(pop2015_ln) + s(district, bs = "re", k = 49) +
    s(time, by = aet,  bs = "tp", k = 50) +
    s(time, by = prcp, bs = "tp", k = 30) +
    s(time, by = tmin, bs = "tp", k = 100),
  family = nb(), data = monthly_cases_analysis, method = "fREML", discrete = TRUE
)
toc()

cat("\nConverged?", falciparum_tvcm$converged, "\n")

summary(falciparum_tvcm)
```
```{r fig.width=12, fig.height=8}
appraise(falciparum_tvcm, point_col = "steelblue", point_alpha = 0.4)
```

```{r fig.width=12, fig.height=8}
draw(falciparum_tvcm, residuals = TRUE)
```

```{r}
model_concurvity(falciparum_tvcm) %>% 
  filter(type == "estimate")
```

# Figure 2

To construct Figure 2 on the article, we first create a data frame with the start and end dates of the PAMAFRO intervention.

```{r message=FALSE, warning=FALSE}
pamafro_period <- 
  tibble(
    start = make_datetime(year = 2005, month = 10, day = 1L), 
    end = make_datetime(year = 2010, month = 09, day = 1L)
  )
```

Then we define a function called `get_est` to extract the smooth estimates evaluated over a range of covariate values in `data` for each model using the `smooth_estimates` of the `gratia` package. 

```{r}
get_est <- function(vivax_model, falciparum_model, data) {
  vivax_model %>% 
  smooth_estimates(
    smooth = c("s(time):aet", "s(time):prcp", "s(time):tmin"), data = data
  ) -> vivax_smooth_est
  
  falciparum_model %>% 
  smooth_estimates(
    smooth = c("s(time):aet", "s(time):prcp", "s(time):tmin"), data = data
  ) -> falciparum_smooth_est
  
  smooth_est_list <- list(
    "vivax" = vivax_smooth_est, "falciparum" = falciparum_smooth_est
  )
  
  full_smooth_est <- 
  bind_rows(smooth_est_list, .id = "specie") %>% 
  mutate(dttm = as_datetime(time))
  
  full_smooth_est
}
```

We create a data frame for the `smooth_estimates` to evaluate over. We need to include all the terms inthe model, but we are going to extract just the estimates for the smooths associated with the covariates `aet`, `prcp` and `tmin`, so it does not matter which values we put in the other terms (`pop2015_ln` and `district`. And since the `smooth_estimates` function evaluates the smooth times the covariate value, that is $s\left(t\right)x$, we fill with ones the values for the covariates, since we just want $s\left(t\right)$.

```{r}
dummy_data <- 
  tibble(
    pop2015_ln = rep(1, 500), district = factor(rep(1, 500)), 
    time = seq_min_max(monthly_cases_analysis$time, 500), aet = rep(1, 500), 
    prcp = rep(1, 500), tmin = rep(1, 500)
)
```

Now we use the function we created to extract the smooth estimates.

```{r}
smooth_est <- get_est(vivax_tvcm, falciparum_tvcm, dummy_data)
```

Display the smooth estimates.

```{r}
smooth_est
```

Here, we define a function to plot the smooth estimates and add a PAMAFRO period layer and some labels for the final figure.

```{r}
plot_tve <- function(smooth_est, pamafro_period, labels) {
  tve_plot <- 
  smooth_est %>% 
  ggplot() +
  facet_wrap(
    ~ by, ncol = 1, scales = "free_y", strip = "left", 
    labeller = labeller(by = labels)
  ) +
  geom_line(aes(dttm, est, color = specie), size = 1) +
  geom_ribbon(
    aes(dttm, est, ymin = est + 2*se, ymax = est - 2*se, fill = specie), 
    alpha = 0.1
  ) +
  geom_rect(
    data = pamafro_period, aes(xmin = start, xmax = end), ymin = -100, ymax = 100, 
    alpha = 0.2
  ) +
  geom_vline(
    xintercept = c(pamafro_period$start, pamafro_period$end), linetype = "dashed", 
    alpha = 0.6
  ) +
  geom_hline(aes(yintercept = 0), alpha = 0.6) +
  scale_x_datetime(
    date_labels = "%Y", date_breaks = "1 year", expand = expansion(add = 0)
  ) +
  labs(y = NULL, x = NULL) +
  theme(
    legend.position = "top", legend.title = element_blank(), 
    strip.placement = "outside", strip.background = element_blank()
  ) +
  scale_color_npg(labels = c("P. vivax", "P. falciparum")) +
  scale_fill_npg(labels = c("P. vivax", "P. falciparum"))
  
  tve_plot
}
```

Define the term labels.

```{r}
covariate_labels <- c(
  aet = "Actual evapotranspiration\ntime-varying coefficients (linear predictor)",
  prcp = "Precipitation\ntime-varying coefficients (linear predictor)",
  tmin = "Minimum temperature\ntime-varying coefficients (linear predictor)"
)
```

Plot Figure 2.

```{r fig.width=10, fig.height=9}
plot_tve(smooth_est, pamafro_period, covariate_labels)
```

```{r}
# png(
#   filename = "figs/figure-2.png",
#   width = 18.5,
#   height = 22.5,
#   units = "cm",
#   pointsize = 12,
#   res = 300
# )
# 
# tve_plot
# 
# dev.off()
```

```{r}
# pdf("figs/figure-2.pdf", width = 6, height = 8)
# tve_plot
# dev.off()
```

# Table 2

```{r}
table_02_sketch = withTags(
  table(
    class = "display",
    thead(
      tr(
        th(rowspan = 3, "Climate variable"),
        th(colspan = 4, "Before"),
        th(colspan = 4, "During"),
        th(colspan = 4, "After")
      ),
      tr(
        lapply(rep(c("Min", "Max"), 3), th, colspan = 2)
      ),
      tr(
        lapply(rep(c("Estimate", "95% CI"), 6), th)
      )
    )
  )
)
```

```{r}
smooth_est %>% 
  mutate(
    period = case_when(
      dttm < pamafro_period$start ~ "before",
      dttm > pamafro_period$end ~ "after",
      TRUE ~ "during"
    )
  ) %>% 
  mutate(period = factor(period, c("before", "during", "after"))) ->
  smooth_est_period
```

```{r}
smooth_est_period %>% 
  group_by(by , period, specie) %>% 
  slice_max(est, n = 1) %>% 
  ungroup() -> smooth_est_max
```

```{r}
smooth_est_period %>% 
  group_by(by, period, specie) %>% 
  slice_min(est, n = 1) %>% 
  ungroup() -> smooth_est_min
```

```{r}
smooth_est_extrema_list <- list("max" = smooth_est_max, "min" = smooth_est_min)
```

```{r}
smooth_est_extrema <- bind_rows(smooth_est_extrema_list, .id = "extrema")
```

```{r}
smooth_est_extrema %>% 
  mutate(ci = sprintf("[%s, %s]", round(est - 2*se, 2), round(est + 2*se, 2))) -> 
  smooth_est_extrema
```

## P. vivax

```{r message=FALSE}
smooth_est_extrema %>% 
  filter(specie == "vivax") %>% 
  select(by, period, extrema, est, ci) %>%
  mutate(est = round(est, 4)) %>% 
  pivot_wider(
    by, names_from  = c(period, extrema), values_from = c(est, ci), names_sep = "_"
  ) %>% 
  select(by, ends_with("min"), ends_with("max")) -> smooth_est_extrema_vivax
```

```{r}
smooth_est_extrema_vivax %>% 
  column_to_rownames(var = "by") %>% 
  datatable(
    container = table_02_sketch, extensions = 'Buttons',
    options = list(
      dom = "Bt", 
      scrollX = TRUE,
      columnDefs = list(list(className = 'dt-center', targets = 6)),
      buttons = c('copy', 'csv')
    ), 
    caption = "P. vivax"
  )
```

## P. falciparum

```{r message=FALSE}
smooth_est_extrema %>% 
  filter(specie == "falciparum") %>% 
  select(by, period, extrema, est, ci) %>%
  mutate(est = round(est, 4)) %>% 
  pivot_wider(
    by, names_from  = c(period, extrema), values_from = c(est, ci), names_sep = "_"
  ) %>% 
  select(by, ends_with("min"), ends_with("max")) -> smooth_est_extrema_falciparum
```

```{r}
smooth_est_extrema_falciparum %>% 
  column_to_rownames(var = "by") %>% 
  datatable(
    container = table_02_sketch, extensions = 'Buttons',
    options = list(
      dom = "Bt", 
      scrollX = TRUE,
      columnDefs = list(list(className = 'dt-center', targets = 6)),
      buttons = c('copy', 'csv')
    ), 
    caption = "P. falciparum"
  )
```

```{r}
smooth_est_period %>% 
  group_by(period, specie, by) %>% 
  filter(est == min(est) | est == max(est))
```