---
title: "Sensitivity analyses for Pf and Pv"
output:
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
    code_folding: show
    toc: true
    number_sections: false
pkgdown:
  as_is: true
vignette: >
  %\VignetteIndexEntry{Sensitivity analyses for Pf and Pv}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 7,
  fig.height = 12,
  fig.asp = NULL,
  out.width = "95%"
)
```

## Overview

We defined triangular distributions for 18 models parameters, with bounds of \(\pm 20\%\) about the baseline values.
We used Latin hypercube sampling (LHS) to draw 1000 parameter combinations from these distributions, and simulated the model for each combination, to assess the model sensitivity.

The uninfected and infected red blood cell (RBC) populations in the circulation and spleen were compared to data collected from asymptomatic Papuan patients who underwent splenectomy.
Circulating RBC counts were also compared to data collected from asymptomatic Papuan patients who had not experienced trauma that required splenectomy.

```{r setup, class.source = 'fold-hide'}
library(spleenrbc, warn.conflicts = FALSE)
library(dplyr, warn.conflicts = FALSE)
library(tidyr)
library(ggplot2)
library(gghalves)
suppressPackageStartupMessages(library(here))

output_dir <- here("outputs")

# Load Pf results.
pf_results <- readRDS(
  file.path(output_dir, "sensitivity-analysis-pf-rbc-cints-all.rds")
) |>
  mutate(day = time / 24)

# Load Pv results.
pv_results <- readRDS(
  file.path(output_dir, "sensitivity-analysis-pv-rbc-cints-all.rds")
) |>
  mutate(day = time / 24)

load_clinical_data <- function(species = c("Pf", "Pv")) {
  species <- match.arg(species)

  # Load HHS total RBC counts for asymptomatic Pf infections.
  data(hhs_asymp, envir = environment())
  df_asymp <- hhs_asymp |>
    filter(species == !!species) |>
    mutate(kind = "Uninfected RBCs: Circulation")

  # Load uRBC and iRBC counts for splenectomised Pf patients.
  data(splenectomised, envir = environment())
  df_splen <- splenectomised |>
    filter(infection == !!species) |>
    pivot_longer(! c("patient", "infection")) |>
    filter(! is.na(value)) |>
    mutate(
      cell_type = case_when(
        grepl("^urbc", name) ~ "Uninfected RBCs",
        grepl("^irbc", name) ~ "Infected RBCs"
      ),
      location = case_when(
        grepl("_spleen", name) ~ "Spleen",
        grepl("_peri", name) ~ "Circulation"
      ),
      kind = paste0(cell_type, ": ", location),
      cell_type = NULL,
      location = NULL
    )

  # Calculate the ratio of iRBC % in spleen vs circulation.
  df_ratio <- df_splen |>
    select(patient, infection, name, value) |>
    pivot_wider(names_from = "name", values_from = "value") |>
    mutate(
      ratio = (irbc_spleen / urbc_spleen) / (irbc_peri / urbc_peri),
      day = 152
    ) |>
    select(day, ratio) |>
    filter(is.finite(ratio))

  list(asymp = df_asymp, splen = df_splen, ratio = df_ratio)
}

pf_clinical <- load_clinical_data("Pf")
pv_clinical <- load_clinical_data("Pv")
```

```{r plot-results-vs-data, class.source = 'fold-hide'}
plot_results_vs_data <- function(df, clinical) {
  # Plot uninfected and infected RBC populations against clinical data.
  quantiles <- c(0.05, 0.25, 0.5, 0.75, 0.95)

  labels <- c(
    "U_t" = "Uninfected RBCs: Circulation",
    "Ur_t" = "Uninfected RBCs: Spleen",
    "I_t" = "Infected RBCs: Circulation",
    "Ir_t" = "Infected RBCs: Spleen",
    "iRBC_Ratio" = "iRBC: (Spleen % / Circulation %)"
  )

  df <- df |>
    filter(name %in% names(labels)) |>
    mutate(kind = labels[name])

  df_splen <- bind_rows(
    clinical$splen |>
      select(value, kind),
    clinical$ratio |> mutate(
      value = ratio,
      kind = "iRBC: (Spleen % / Circulation %)"
    ) |>
      select(value, kind)
  )

  ggplot(df, aes()) +
    geom_half_violin(
      data = clinical$asymp,
      mapping = aes(2.5, total_rbc),
      linetype = "dashed",
      draw_quantiles = quantiles,
      side = "l",
      width = 5
    ) +
    geom_ribbon(
      mapping = aes(day, ymin = ci_min, ymax = ci_max,
                    fill = cred_int, colour = cred_int),
      linewidth = 0.5
    ) +
    geom_half_violin(
      data = df_splen,
      mapping = aes(max(df$day) + 2, value),
      position = position_identity(),
      draw_quantiles = quantiles,
      side = "r",
      width = 5
    ) +
    geom_jitter(
      data = df_splen,
      mapping = aes(max(df$day) + 7, value),
      colour = "#af3f3f",
      shape = 8,
      width = 2
    ) +
    scale_x_continuous(
      "Time (days)",
      breaks = scales::breaks_width(10)
    ) +
    scale_y_log10("Number of cells") +
    scale_colour_brewer(NULL) +
    scale_fill_brewer(NULL) +
    facet_wrap(~ kind, ncol = 1, scales = "free_y") +
    guides(colour = guide_none())
}
```

## Plasmodium falciparum

```{r plot-pf}
#| class.source = "fold-hide",
#| fig.cap = "Uninfected and infected RBC populations in the circulation and
#|     spleen (blue ribbons). These are compared to data from asymptomatic
#|     Papuans (violins) who underwent splenectomy (solid lines) and who
#|     retained their spleens (dashed lines)."
plot_results_vs_data(pf_results, pf_clinical)
```

## Plasmodium vivax

```{r plot-pv}
#| class.source = "fold-hide",
#| fig.cap = "Uninfected and infected RBC populations in the circulation and
#|     spleen (blue ribbons). These are compared to data from asymptomatic
#|     Papuans (violins) who underwent splenectomy (solid lines) and who
#|     retained their spleens (dashed lines)."
plot_results_vs_data(pv_results, pv_clinical)
```