---
title: "Exploratory Data Analysis in anumaan"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Exploratory Data Analysis in anumaan}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 11,
  fig.height = 6,
  out.width = "100%"
)
```

```{r libraries}
library(anumaan)
library(dplyr)
```

# Overview

This vignette covers the **exploratory data analysis (EDA) layer**. It sits
between the two other pipelines: it is run on data that has already been
through [Preprocessing](preprocessing-workflow.html), and it is what a
stewardship team typically reviews *before* committing to
[DALY burden estimation](daly-pipeline.html). It does not transform data (like
Pipeline 1) or estimate burden (like Pipeline 3) -- it only visualises what is
already there, to catch data-quality problems and characterise the cohort
before further modelling.

All functions in this layer are named `plot_*` (with the shared styling
helper `eda_theme()`), which is deliberately distinct from the `prep_*`
(Pipeline 1) and `daly_*` (Pipeline 3) naming used elsewhere in the package.

Most `plot_*` functions here share a common `mode` argument:

- **`"faceted"`** -- one panel per centre (the default for most functions)
- **`"overall"`** -- all centres pooled into a single chart
- **`"single"`** -- one specific centre, passed via `center = "..."`

and a common `syndrome_col`/`syndrome_name` pair for restricting the plot to
one infectious syndrome (e.g. `syndrome_name = "BSI"`) when that column is
available.

## Connecting from Pipeline 1

The `plot_*` functions' *default* column names (`organism_name`,
`antibiotic_value`, `PatientInformation_id`, ...) are convenience defaults for
raw-ish data -- they do **not** automatically match the standardised column
names [Preprocessing](preprocessing-workflow.html) actually produces. When
passing real preprocessed output into these functions, map the arguments to
Pipeline 1's real output columns explicitly:

| EDA argument | Pipeline 1's actual output column |
|---|---|
| `patient_col` | `patient_id` |
| `organism_col` | `organism_normalized` (from `prep_standardize_organisms()`) |
| `antibiotic_col` | `antibiotic_normalized` (from `prep_standardize_antibiotics()`) |
| `value_col` | `ast_value_harmonized` (from `prep_harmonize_ast()`) |
| `sample_col` | `specimen_normalized` (from `prep_standardize_specimens()`) |
| `agebin_col` | `Age_bin` (from `prep_assign_age_bins()`) |
| `age_col` | `Age` (from `prep_fill_age()`) |
| `admission_col` | `date_of_admission` |
| `culture_col` | `date_of_culture` |
| `discharge_col` | `date_of_final_outcome` |
| `center_col` | `center_name` (a user-supplied raw column -- Pipeline 1 does not invent hospital identity, it just carries it through) |
| `outcome_col` | `final_outcome` (the raw outcome text persists unchanged; `prep_standardize_final_outcome()` *adds* a separate `outcome_std` column alongside it rather than replacing it) |

Every example below uses these Pipeline-1-aligned names, with the arguments
passed explicitly rather than relying on defaults.

# A Small Synthetic Dataset

The examples below use a small synthetic long-format dataset (one row per
isolate x antibiotic) with the same column names Pipeline 1 would hand off --
as if it were the output of the [Preprocessing](preprocessing-workflow.html)
vignette's `prepped` object, scaled up to two hospitals so the `mode =
"faceted"` examples have something to facet over.

```{r synthetic-data}
set.seed(42)
n_pts <- 60
centers <- c("General Hospital", "City Medical Center")
organisms <- c("Escherichia coli", "Klebsiella pneumoniae",
                "Staphylococcus aureus", "Pseudomonas aeruginosa")
antibiotics <- c("Amikacin", "Ciprofloxacin", "Meropenem", "Ceftriaxone")

eda_data <- data.frame(
  patient_id             = paste0("pt_", sprintf("%03d", rep(1:n_pts, each = 3))),
  center_name            = rep(sample(centers, n_pts, replace = TRUE), each = 3),
  organism_normalized    = rep(sample(organisms, n_pts, replace = TRUE), each = 3),
  antibiotic_normalized  = rep(antibiotics, length.out = n_pts * 3),
  ast_value_harmonized   = sample(c("R", "S"), n_pts * 3, replace = TRUE, prob = c(0.4, 0.6)),
  specimen_normalized    = rep(sample(c("Blood", "Urine", "Sputum"), n_pts, replace = TRUE), each = 3),
  final_outcome          = rep(sample(c("Death", "Discharged", "LAMA", "Transferred to other hospital"),
                                       n_pts, replace = TRUE, prob = c(0.15, 0.65, 0.1, 0.1)), each = 3),
  date_of_final_outcome  = rep(as.Date("2024-01-01") + sample(1:300, n_pts, replace = TRUE), each = 3),
  is_polymicrobial       = rep(sample(c(0, 1), n_pts, replace = TRUE, prob = c(0.75, 0.25)), each = 3),
  date_of_admission      = rep(as.Date("2024-01-01") + sample(1:300, n_pts, replace = TRUE), each = 3),
  location               = rep(sample(c("ICU", "Ward"), n_pts, replace = TRUE), each = 3),
  infection_type         = rep(sample(c("HAI", "CAI"), n_pts, replace = TRUE), each = 3),
  Age_bin                = rep(sample(c("<1", "1-5", "18-45", "45-65", "65+"), n_pts, replace = TRUE), each = 3),
  Age                    = rep(sample(0:85, n_pts, replace = TRUE), each = 3),
  infectious_syndrome    = rep(sample(c("BSI", "UTI", "RTI"), n_pts, replace = TRUE), each = 3),
  stringsAsFactors = FALSE
)
eda_data$date_of_culture <- eda_data$date_of_admission + sample(0:5, nrow(eda_data), replace = TRUE)

nrow(eda_data)
```

# Theme

`eda_theme()` is the shared visual style behind every plot in this vignette
-- based on `ggpubr::theme_pubr()` when `ggpubr` is installed, falling back
to `ggplot2::theme_bw()` otherwise. You do not normally call it directly; it
is applied internally by every `plot_*` function.

```{r theme, eval = FALSE}
eda_theme(base_size = 14)
```

# Organism and Susceptibility Overview

## Top organisms

```{r top-organisms, fig.width = 14, fig.height = 5.5}
plot_top_organisms(
  eda_data, n = 3, mode = "faceted",
  patient_col = "patient_id", organism_col = "organism_normalized", center_col = "center_name"
)
```

Pool across all centres with `mode = "overall"`, or restrict to one centre
with `mode = "single", center = "General Hospital"`.

## Susceptibility pattern (stacked R/S bars)

`plot_abx_susceptibility()` applies a **worst-phenotype rule**: if a patient
has multiple records for the same organism-antibiotic combination, a single R
result marks the whole episode as R, so repeat cultures from one infection
episode are not double-counted.

```{r abx-susceptibility}
plot_abx_susceptibility(
  eda_data, n = 3, mode = "overall",
  patient_col = "patient_id", antibiotic_col = "antibiotic_normalized",
  value_col = "ast_value_harmonized", organism_col = "organism_normalized",
  center_col = "center_name"
)
```

## Resistance heatmap

```{r abx-heatmap, fig.width = 11, fig.height = 4.5}
plot_abx_heatmap(
  eda_data, mode = "all",
  patient_col = "patient_id", antibiotic_col = "antibiotic_normalized",
  value_col = "ast_value_harmonized", organism_col = "organism_normalized",
  center_col = "center_name"
)
```

## Resistance by specimen type

```{r resistance-by-sample}
plot_resistance_by_sample(
  eda_data, n = 3, mode = "overall",
  patient_col = "patient_id", sample_col = "specimen_normalized",
  organism_col = "organism_normalized", antibiotic_col = "antibiotic_normalized",
  value_col = "ast_value_harmonized", center_col = "center_name"
)
```

## Resistance by age bin and by organism

```{r resistance-by-agebin, fig.width = 12, fig.height = 5}
plot_resistance_by_agebin(
  eda_data, mode = "overall",
  patient_col = "patient_id", organism_col = "organism_normalized",
  antibiotic_col = "antibiotic_normalized", value_col = "ast_value_harmonized",
  agebin_col = "Age_bin", center_col = "center_name"
)
```

```{r resistance-by-organism, fig.width = 12, fig.height = 5}
plot_resistance_by_organism(
  eda_data, n = 3, mode = "overall",
  patient_col = "patient_id", organism_col = "organism_normalized",
  antibiotic_col = "antibiotic_normalized", value_col = "ast_value_harmonized",
  center_col = "center_name"
)
```

# Outcomes

## Overall distribution

By default `merge_referred = TRUE` recodes `"Transferred to other hospital"`
to `"Referred"`, since both represent the same clinical event across centres.

```{r outcome-distribution}
plot_outcome_distribution(
  eda_data, mode = "overall",
  patient_col = "patient_id", outcome_col = "final_outcome", center_col = "center_name"
)
```

## Outcome by organism, split by resistance status

`resistance_filter = "R"` shows patients with *at least one* R result for
that organism; `"S"` shows patients with only S results. Because these groups
overlap, calling the function once for each gives the fuller picture.

```{r outcome-by-organism}
plot_outcome_by_organism(
  eda_data, n = 3, resistance_filter = "R", mode = "overall",
  patient_col = "patient_id", organism_col = "organism_normalized",
  antibiotic_col = "antibiotic_normalized", value_col = "ast_value_harmonized",
  outcome_col = "final_outcome", center_col = "center_name"
)
```

## Death vs. discharged

```{r death-discharged}
plot_death_discharged(
  eda_data, n = 3, mode = "overall",
  patient_col = "patient_id", organism_col = "organism_normalized",
  outcome_col = "final_outcome", center_col = "center_name",
  antibiotic_col = "antibiotic_normalized", value_col = "ast_value_harmonized"
)
```

## Outcome by age bin and by year

```{r outcome-by-agebin}
plot_outcome_by_agebin(
  eda_data, mode = "overall",
  patient_col = "patient_id", agebin_col = "Age_bin",
  outcome_col = "final_outcome", center_col = "center_name"
)
```

```{r outcome-by-year}
plot_outcome_by_year(
  eda_data, mode = "overall",
  patient_col = "patient_id", outcome_col = "final_outcome",
  date_col = "date_of_final_outcome", center_col = "center_name"
)
```

# Case-Mix and Facility Patterns

## Patients per hospital

```{r patients-by-hospital}
plot_patients_by_hospital(eda_data, patient_col = "patient_id", center_col = "center_name")
```

## Monomicrobial vs. polymicrobial, by facility

```{r mono-poly}
plot_mono_poly_by_facility(
  eda_data, patient_col = "patient_id", poly_col = "is_polymicrobial", center_col = "center_name"
)
```

## HAI vs. CAI, by facility

Computed from the gap between `admission_col` and `culture_col`, using
`hai_threshold_hours` (default 48h) when an explicit infection-type column
is not already available.

```{r hai-cai}
plot_hai_cai_by_facility(
  eda_data, patient_col = "patient_id", center_col = "center_name",
  admission_col = "date_of_admission", culture_col = "date_of_culture"
)
```

## Location (ICU vs. Ward), by facility

```{r location-facility}
plot_location_by_facility(
  eda_data, patient_col = "patient_id", location_col = "location", center_col = "center_name"
)
```

## Infection type by location

```{r infection-type-location}
plot_infection_type_by_location(
  eda_data, style = "bar", patient_col = "patient_id", center_col = "center_name",
  location_col = "location", infection_col = "infection_type"
)
```

## Syndrome distribution

```{r syndrome-distribution}
plot_syndrome_distribution(
  eda_data, mode = "overall", patient_col = "patient_id",
  syndrome_col = "infectious_syndrome", center_col = "center_name"
)
```

# Length of Stay and Age Distributions

## LOS ridge plot

```{r los-ridge, fig.width = 11, fig.height = 5.5}
plot_los_ridge(
  eda_data, mode = "all", patient_col = "patient_id", center_col = "center_name",
  admission_col = "date_of_admission", discharge_col = "date_of_final_outcome"
)
```

## Age ridge plot

```{r age-ridge, fig.width = 11, fig.height = 5.5}
plot_age_ridge(eda_data, mode = "all", patient_col = "patient_id", age_col = "Age", center_col = "center_name")
```

## LOS by age bin

```{r los-by-agebin, fig.width = 12, fig.height = 5.5}
plot_los_by_agebin(
  eda_data, mode = "overall", patient_col = "patient_id", agebin_col = "Age_bin",
  age_col = "Age", center_col = "center_name",
  admission_col = "date_of_admission", discharge_col = "date_of_final_outcome"
)
```

---

# Pipeline at a Glance

```
Preprocessing output (prep_* pipeline, analysis-ready long-format data)
          │
          ▼
  plot_top_organisms() / plot_abx_susceptibility() / plot_abx_heatmap()
  plot_resistance_by_sample() / plot_resistance_by_agebin() / plot_resistance_by_organism()
          │            (organism & susceptibility overview)
          ▼
  plot_outcome_distribution() / plot_outcome_by_organism() / plot_death_discharged()
  plot_outcome_by_agebin() / plot_outcome_by_year()
          │            (outcome patterns)
          ▼
  plot_patients_by_hospital() / plot_mono_poly_by_facility() / plot_hai_cai_by_facility()
  plot_location_by_facility() / plot_infection_type_by_location() / plot_syndrome_distribution()
          │            (case-mix and facility patterns)
          ▼
  plot_los_ridge() / plot_age_ridge() / plot_los_by_agebin()
          │            (LOS / age distributions)
          ▼
  Reviewed cohort -> proceed to DALY Burden Estimation (Pipeline 3)
```

# Summary

The EDA layer is meant to be run interactively, plot by plot, while deciding
whether the preprocessed cohort is ready for burden estimation -- not wrapped
into a single automated entry point the way `run_preprocess()` is for
Pipeline 1. A typical review pass:

1. Check organism and susceptibility patterns for anything unexpected
   (`plot_top_organisms()`, `plot_abx_heatmap()`)
2. Check outcome distributions and case-mix balance across centres
3. Check LOS/age distributions for implausible values before they reach
   burden estimation
4. Only then move on to the
   [DALY Burden Estimation vignette](daly-pipeline.html)

# Session Info

```{r session}
sessionInfo()
```