Step_04_Weekly_Transition_Probabilities.qmd

---
title: "Step 04 - Derive Weekly Transition Probabilities from Weibull Fits"
author: "Jason van Leuven"
date: "2025-10-25"
format:
  html:
    toc: true
execute:
  echo: true
  warning: false
  message: false
---

# Step 04. Derive Weekly Transition Probabilities from Weibull Fits

This step translates the continuous Weibull survival curves estimated in Step 03 into discrete-time transition probabilities, enabling their use in Markov cohort models.\
While the Weibull distribution describes healing as a smooth continuous process, health-economic models require state transitions at fixed intervals, here, weekly steps that capture the likelihood of moving from an active wound or infection state to a healed or cleared one.

Conceptually, the survival function $(S(t))$ represents the probability that healing has *not yet occurred* by time $(t)$.\
To obtain a usable transition probability between discrete intervals, the model computes the relative decline in survival from week $(t)$ to $(t!+!1)$.\
This transformation:

$$
p_t = 1 - \frac{S(t+1)}{S(t)}
$$

defines the weekly probability of event occurrence (e.g., wound closure or bacterial clearance).\
The result is a set of group-specific probability vectors, one for each cohort x treatment pair, that together form the empirical backbone of the Markov transition matrix used in later steps.

From a modelling perspective, this conversion achieves three things.\
First, it aligns continuous survival data with the discrete-time logic of state-based economic evaluation.\
Second, it preserves the shape and scale characteristics of the Weibull distribution, maintaining biological realism while enabling weekly cycle simulation.\
Third, it provides a transparent audit trail linking experimental data to the probabilistic structure used in Step 05 (scaling) and Step 07 (calibration).\
Each exported `.csv` file therefore serves as both a traceable data product and a reproducible modelling artefact.

### Mathematical Note - Discretisation Logic

In continuous time, the Weibull survival curve $S(t) = e{-(t/\lambda)k}$ expresses the probability that a wound remains unhealed at time (t).\
To translate this into a weekly Markov model, the change in survival between two consecutive time points reflects the probability that healing occurred within that week — conditional on still being unhealed the week before.

Formally, the conditional weekly probability is:

$$
p_t = 1 - \frac{S(t+1)}{S(t)} = 1 - e{-((t+1)/\lambda)k + (t/\lambda)^k} 
$$

This ratio structure ensures the total probability remains consistent with the original continuous hazard.\
In effect, $(S(t+1)/S(t))$ represents the fraction of the population that *survives one more week without healing*, so its complement $(1 - S(t+1)/S(t))$ gives the weekly transition probability to a healed or cleared state.

This discretisation approach preserves the essential shape of the underlying Weibull hazard, whether increasing $((k>1))$, constant $((k=1))$, or decreasing $((k<1))$ ; while making the results directly compatible with cycle-based health-economic models.


```{r setup_infrastructure}
# Load V2 Infrastructure
source("config/paths.R")
source("src/R/utils.R")
source("src/R/data_validation.R")
source("src/R/logging.R")

# Get standardized paths
paths <- get_project_paths()

# Setup logging
setup_logging(log_dir = paths$logs, log_level = "INFO")

# Load required packages
suppressPackageStartupMessages({
  library(tidyverse)
  library(glue)
})

log_info("=", paste(rep("=", 69), collapse = ""))
log_info("Step_04_Weekly_Transition_Probabilities")
log_info("=", paste(rep("=", 69), collapse = ""))
```

------------------------------------------------------------------------

```{r step04a_load}
# STEP 04A - Load Weibull fit tables
# Purpose: Import fitted parameters from Step 03 for closure and clearance cohorts.

library(tidyverse)
library(glue)

# Using paths from get_project_paths()
# output_path replaced by paths$outputs
weibull_closure_file   <- file.path(paths$outputs, "weibull_fits_closure.csv")
weibull_clearance_file <- file.path(paths$outputs, "weibull_fits_clearance.csv")

closure_fits   <- read_csv(weibull_closure_file, show_col_types = FALSE)
clearance_fits <- read_csv(weibull_clearance_file, show_col_types = FALSE)

cat(glue("Loaded {nrow(closure_fits)} closure fits and {nrow(clearance_fits)} clearance fits.\n"))
```

------------------------------------------------------------------------

## 2. Define Weibull - Transition Function

The Weibull survival function is expressed as:

$$ S(t) = e^{-(t/\lambda)^k} $$

The discrete-time probability of transition between week *t-1* and *t* is defined as:

$$ p_t = 1 - \frac{S(t)}{S(t-1)} = 1 - e^{-((t/\lambda)^k - ((t-1)/\lambda)^k)} $$

This transformation provides the weekly probability of closure or clearance.

```{r step04b_define_function}
# STEP 04B - Define conversion from Weibull to weekly transitions
# Reference: Adapted from parametric survival to discrete Markov conversion (Briggs et al., 2006).

weibull_to_transition <- function(shape, scale, horizon_weeks = 12) {
  t <- seq_len(horizon_weeks)
  S <- exp(-(t / scale)^shape)
  S_prev <- c(1, S[-length(S)])   # Define S(0) = 1 baseline
  p <- 1 - (S / S_prev)
  tibble(week = t, S = S, p = p)
}

cat("Weibull -> Transition conversion function defined.\n")
```

------------------------------------------------------------------------

## 3. Generate Weekly Probabilities per Group

```{r step04c_generate}
# STEP 04C - Generate transition probability series per cohort x treatment pair

cat("\nSTEP 04C - GENERATE WEEKLY PROBABILITIES\n============================\n")

# Closure transitions
closure_transitions <- closure_fits %>%
  filter(converged) %>%
  mutate(data = map2(shape, scale, weibull_to_transition)) %>%
  unnest(data) %>%
  mutate(type = "closure")

# Clearance transitions
clearance_transitions <- clearance_fits %>%
  filter(converged) %>%
  mutate(data = map2(shape, scale, weibull_to_transition)) %>%
  unnest(data) %>%
  mutate(type = "clearance")

cat(glue("Generated transitions for {n_distinct(closure_transitions$treatment_group)} closure and {n_distinct(clearance_transitions$treatment_group)} clearance treatments.\n"))
```

------------------------------------------------------------------------

## 4. Visualise Transition Curves

```{r step04d_visualise, fig.width=14, fig.height=12}
# STEP 04D - Visualise transition probability curves
# Purpose: Verify smooth decay patterns and between-group separation.

library(ggplot2)
library(patchwork)

cat("\nSTEP 04D - VISUAL DIAGNOSTIC PLOTS\n============================\n")

log_info("Generating enhanced transition probability visualizations")

# Enhanced Closure transition probabilities
p_closure_trans <- ggplot(closure_transitions, aes(x = week, y = p, colour = treatment_group)) +
  geom_line(linewidth = 1.3) +
  geom_point(size = 2.5) +
  facet_wrap(~ cohort + treatment_group, scales = "free_y", ncol = 3) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 11),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "none",
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "A. Weekly Closure Transition Probabilities",
    subtitle = "Probability of wound healing in each week, conditional on not yet being healed",
    x = "Week (Human Timescale After Scaling)", 
    y = "Probability (Infected -> Healing/Healed)"
  )

# Enhanced Clearance transition probabilities
p_clearance_trans <- ggplot(clearance_transitions, aes(x = week, y = p, colour = treatment_group)) +
  geom_line(linewidth = 1.3) +
  geom_point(size = 2.5) +
  facet_wrap(~ cohort + treatment_group, scales = "free_y", ncol = 3) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 11),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "none",
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "B. Weekly Clearance Transition Probabilities",
    subtitle = "Probability of bacterial clearance in each week, conditional on not yet being cleared",
    x = "Week (Human Timescale After Scaling)", 
    y = "Probability (Infected -> Cleared)"
  )

p_closure_trans
p_clearance_trans

# NEW: Cumulative transition probability (shows eventual healing/clearance)
closure_transitions <- closure_transitions %>%
  group_by(cohort, treatment_group) %>%
  arrange(week) %>%
  mutate(cum_prob = cumsum(p * cumprod(c(1, 1 - p[-n()])))) %>%
  ungroup()

clearance_transitions <- clearance_transitions %>%
  group_by(cohort, treatment_group) %>%
  arrange(week) %>%
  mutate(cum_prob = cumsum(p * cumprod(c(1, 1 - p[-n()])))) %>%
  ungroup()

p_closure_cum <- ggplot(closure_transitions, aes(x = week, y = cum_prob, colour = treatment_group)) +
  geom_line(linewidth = 1.3) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0.5, linetype = "dashed", colour = "gray40", linewidth = 0.8) +
  geom_hline(yintercept = 0.9, linetype = "dotted", colour = "gray40", linewidth = 0.8) +
  facet_wrap(~ cohort + treatment_group, ncol = 3) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 11),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "none",
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "C. Cumulative Wound Healing Probability",
    subtitle = "Dashed line = 50% healed (median), dotted line = 90% healed",
    x = "Week",
    y = "Cumulative Probability of Healing"
  )

p_clearance_cum <- ggplot(clearance_transitions, aes(x = week, y = cum_prob, colour = treatment_group)) +
  geom_line(linewidth = 1.3) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0.5, linetype = "dashed", colour = "gray40", linewidth = 0.8) +
  geom_hline(yintercept = 0.9, linetype = "dotted", colour = "gray40", linewidth = 0.8) +
  facet_wrap(~ cohort + treatment_group, ncol = 3) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 11),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "none",
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "D. Cumulative Bacterial Clearance Probability",
    subtitle = "Dashed line = 50% cleared (median), dotted line = 90% cleared",
    x = "Week",
    y = "Cumulative Probability of Clearance"
  )

p_closure_cum
p_clearance_cum

# NEW: Hazard rate visualization
closure_transitions <- closure_transitions %>%
  mutate(hazard = -log(1 - p))

clearance_transitions <- clearance_transitions %>%
  mutate(hazard = -log(1 - p))

p_closure_hazard <- ggplot(closure_transitions, aes(x = week, y = hazard, colour = treatment_group)) +
  geom_line(linewidth = 1.3) +
  geom_point(size = 2) +
  facet_wrap(~ cohort + treatment_group, scales = "free_y", ncol = 3) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 11),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "none",
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "E. Weekly Hazard Rate (Wound Closure)",
    subtitle = "Instantaneous risk of healing; upward trend indicates accelerating healing",
    x = "Week",
    y = "Hazard Rate h(t)"
  )

p_clearance_hazard <- ggplot(clearance_transitions, aes(x = week, y = hazard, colour = treatment_group)) +
  geom_line(linewidth = 1.3) +
  geom_point(size = 2) +
  facet_wrap(~ cohort + treatment_group, scales = "free_y", ncol = 3) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 11),
    strip.text = element_text(face = "bold", size = 10),
    legend.position = "none",
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "F. Weekly Hazard Rate (Bacterial Clearance)",
    subtitle = "Instantaneous risk of clearance; upward trend indicates accelerating clearance",
    x = "Week",
    y = "Hazard Rate h(t)"
  )

p_closure_hazard
p_clearance_hazard

log_info("Enhanced transition probability visualizations generated")

cat("Enhanced diagnostic plots generated.\n")
```

------------------------------------------------------------------------

## 5. Save Transition Tables

```{r step04e_export}
# STEP 04E - Export transition probabilities
# Purpose: Store weekly transition probabilities for later Markov analysis (Step 05).

cat("\nSTEP 04E - EXPORT RESULTS\n============================\n")

# Ensure the output directory exists
if (!dir.exists(output_path)) {
  dir.create(output_path, recursive = TRUE)
  cat(glue("Created output directory -> {output_path}\n"))
}

closure_outfile   <- file.path(paths$outputs, "weekly_transitions_closure.csv")
clearance_outfile <- file.path(paths$outputs, "weekly_transitions_clearance.csv")

write_csv(closure_transitions, closure_outfile)
write_csv(clearance_transitions, clearance_outfile)

if (file.exists(closure_outfile) && file.exists(clearance_outfile)) {
  cat(glue("Saved closure transitions -> {closure_outfile}\n"))
  cat(glue("Saved clearance transitions -> {clearance_outfile}\n"))
  cat("Step 04 complete - weekly transition probabilities exported.\n")
} else {
  cat("Warning: Files not found after export - check write permissions or path.\n")
}
```

------------------------------------------------------------------------

*References:*

-   Briggs, A., Claxton, K., & Sculpher, M. (2006). *Decision Modelling for Health Economic Evaluation.* Oxford University Press.\
-   Cukjati, D., et al. (2001). *A reliable method of determining wound healing rate.* *Medical and Biological Engineering and Computing, 39*(2), 263-271.\
-   Wei, X., et al. (2017). *Allometric scaling of skin thickness, elasticity, viscoelasticity to mass for micro-medical device translation.* *Journal of the Mechanical Behavior of Biomedical Materials, 68*, 303-310.