---
title: "Cancer Survival Analysis in Western Kenya"
description: "A comprehensive survival analysis of cancer patients at KEMRI, focusing on lost-to-follow-up patterns and survival outcomes"
date: "2025-10-19"
image: images/cancer-survival-preview.jpg
categories: ["cancer", "survival-analysis", "kenya", "kemri", "healthcare"]
execute: 
  enabled: true
  freeze: auto
format:
  html:
    code-tools: true
    code-copy: true
    code-line-numbers: true
    toc: true
    toc-depth: 3
    toc-location: left

---

```{r setup, include=FALSE}
# Global options
knitr::opts_chunk$set(
  echo = TRUE,           # Show code
  warning = FALSE,       # Suppress warnings
  message = FALSE,       # Suppress messages
  fig.width = 10,        # Default figure width
  fig.height = 6,        # Default figure height
  out.width = "100%"     # Make figures responsive
)

# Load required packages
library(tidyverse)      # Data manipulation
library(survival)       # Survival analysis
library(survminer)      # Survival visualizations
library(plotly)         # Interactive plots
library(DT)             # Interactive tables
library(lubridate)      # Date handling
library(flexsurv)       # Flexible parametric survival models
library(riskRegression) # Competing risks analysis
library(mstate)         # Multi-state models
library(cmprsk)         # Competing risks
library(prodlim)        # Product limit estimator
```

::: {.highlight-box}
### 🎯 Featured Blog Post: The Kisumu Cancer Dashboard Story

This survival analysis was the foundation for **East Africa's first Cancer Epidemiology Dashboard**, launched by Governor Anyang' Nyong'o in November 2024.

As **Research Data Manager** for this KEMRI project, I led:
- REDCap database design & implementation
- Research assistant training & supervision
- Chart abstraction across 3,438 patients
- Quality checks using Python, R, and Stata
- Survival analysis & epidemiological insights
- Interactive Tableau dashboard development

**Read the full story:** [Transforming Cancer Care: Building East Africa's First Cancer Epidemiology Dashboard →](../../posts/kisumu-cancer-dashboard/index.qmd)

**Impact:** Real-time patient tracking, 59% first-year LTFU identification, data-driven clinical interventions.
:::

## 🔬 Project Overview

This analysis examines cancer survival patterns among patients at Jaramogi Oginga Odinga Teaching and Referral Hospital (JOOTRH) in Kisumu County, focusing on:

- **Lost-to-follow-up (LTFU) patterns** (59% in first year)
- **Survival outcomes** by cancer type and stage (9% 5-year survival)
- **Time-to-event analysis** with competing risks
- **Predictive modeling** of LTFU risk factors

## 📊 Dataset

```{r load-data}
# Simulated dataset based on KEMRI's oncology registry
set.seed(2025)

# Generate patient data
n_patients <- 2000
cancer_types <- c("Esophageal", "Cervical", "Breast", "Prostate", "Kaposi's Sarcoma")
districts <- c("Kisumu", "Kakamega", "Bungoma", "Vihiga", "Busia", "Siaya")

patient_data <- tibble(
  patient_id = 1:n_patients,
  age = round(rnorm(n_patients, 52, 12)),
  sex = sample(c("Male", "Female"), n_patients, replace = TRUE, prob = c(0.45, 0.55)),
  hiv_status = sample(c("Positive", "Negative"), n_patients, replace = TRUE, 
                     prob = c(0.35, 0.65)),
  cancer_type = sample(cancer_types, n_patients, replace = TRUE, 
                      prob = c(0.30, 0.25, 0.20, 0.15, 0.10)),
  district = sample(districts, n_patients, replace = TRUE),
  stage = sample(c("I", "II", "III", "IV"), n_patients, replace = TRUE, 
                prob = c(0.1, 0.2, 0.3, 0.4)),
  date_diagnosis = sample(seq(as.Date('2018-01-01'), as.Date('2024-12-31'), by="day"), 
                         n_patients, replace = TRUE),
  # Simulate LTFU (60% rate)
  lost_to_followup = rbinom(n_patients, 1, 0.6),
  # Simulate death (competing risk)
  death = ifelse(lost_to_followup == 0, rbinom(sum(lost_to_followup == 0), 1, 0.7), 0)
)

# Calculate follow-up time (in months)
max_followup <- as.numeric(difftime(Sys.Date(), min(patient_data$date_diagnosis), units = "days")) / 30.44
patient_data <- patient_data %>%
  mutate(
    # Random follow-up time (0 to max possible)
    followup_months = if_else(
      lost_to_followup == 1,
      runif(n(), 1, 24),  # LTFU: 1-24 months
      pmin(runif(n(), 1, 60),  # Censored or death: up to 5 years
           as.numeric(difftime(Sys.Date(), date_diagnosis, units = "days")) / 30.44)
    ),
    # Event indicator (1=death, 0=censored, 2=LTFU)
    status = case_when(
      death == 1 ~ 1,
      lost_to_followup == 1 ~ 2,
      TRUE ~ 0
    ),
    # Format for survival analysis
    time = followup_months,
    event = factor(status, levels = 0:2, labels = c("Censored", "Death", "LTFU")),
    age_group = cut(age, breaks = c(0, 39, 59, Inf), labels = c("<40", "40-59", "60+"), right = TRUE),
    hiv_positive = hiv_status == "Positive",
    travel_cost_kes = round(runif(n(), 100, 1500), 0),
    income_quintile = sample(1:5, n(), replace = TRUE)
  )

# Save for download
write_csv(patient_data, "kemri_cancer_data.csv")
```

## 📈 1. Patient Characteristics

```{r patient-characteristics}
# Summary table
summary_stats <- patient_data %>%
  group_by(cancer_type) %>%
  summarise(
    Patients = n(),
    `Median Age` = median(age),
    `Female (%)` = round(mean(sex == "Female") * 100, 1),
    `HIV+ (%)` = round(mean(hiv_status == "Positive") * 100, 1),
    `LTFU Rate (%)` = round(mean(lost_to_followup) * 100, 1),
    .groups = 'drop'
  )

# Interactive table
DT::datatable(
  summary_stats,
  rownames = FALSE,
  options = list(
    pageLength = 5,
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel')
  ),
  caption = "Table 1: Patient Characteristics by Cancer Type"
) %>%
  DT::formatStyle(
    'LTFU Rate (%)',
    backgroundColor = styleInterval(
      c(40, 60), 
      c('#d9f3d9', '#fff2cc', '#f8d7da')
    )
  )
```

## 📉 2. Kaplan-Meier Survival Analysis

### Overall Survival

```{r km-overall}
# Fit survival model
fit_os <- survfit(Surv(time, status == 1) ~ 1, data = patient_data)

# Plot
km_plot <- ggsurvplot(
  fit_os,
  data = patient_data,
  risk.table = TRUE,
  pval = TRUE,
  conf.int = TRUE,
  xlab = "Months from Diagnosis",
  ylab = "Overall Survival Probability",
  title = "Figure 1: Overall Survival (All Patients)",
  legend = "none",
  palette = "#E5243B",
  ggtheme = theme_minimal()
)

# Convert to plotly for interactivity
ggplotly(km_plot$plot) %>%
  layout(
    hovermode = "x unified",
    xaxis = list(title = "Months from Diagnosis"),
    yaxis = list(title = "Survival Probability")
  )
```

### By Cancer Type

```{r km-by-type}
# Fit model by cancer type
fit_by_type <- survfit(Surv(time, status == 1) ~ cancer_type, data = patient_data)

# Plot
ggsurvplot(
  fit_by_type,
  data = patient_data,
  conf.int = TRUE,
  pval = TRUE,
  risk.table = TRUE,
  legend.title = "Cancer Type",
  legend.labs = unique(patient_data$cancer_type),
  xlab = "Months from Diagnosis",
  ylab = "Survival Probability",
  title = "Figure 2: Overall Survival by Cancer Type",
  palette = "jco",
  ggtheme = theme_minimal()
)
```

## 🔎 2b. Focused KM: Esophageal & HIV+

```{r km-focused}
# Esophageal cancer cohort
escohort <- patient_data %>% filter(cancer_type == "Esophageal")
fit_es <- survfit(Surv(time, status == 1) ~ 1, data = escohort)

# HIV+ cohort
hivpos <- patient_data %>% filter(hiv_positive)
fit_hiv <- survfit(Surv(time, status == 1) ~ 1, data = hivpos)

p1 <- ggsurvplot(fit_es, data = escohort, risk.table = TRUE, conf.int = TRUE,
                 title = "Esophageal Cancer: Overall Survival", ggtheme = theme_minimal())
p2 <- ggsurvplot(fit_hiv, data = hivpos, risk.table = TRUE, conf.int = TRUE,
                 title = "HIV+ Patients: Overall Survival", ggtheme = theme_minimal())

subplot(ggplotly(p1$plot), ggplotly(p2$plot), nrows = 1, shareY = TRUE, titleX = TRUE)
```

## ⚠️ 3. Lost-to-Follow-Up Analysis

### Cumulative Incidence Function (Competing Risks)

```{r cif-analysis}
# Convert to competing risks format
patient_data_cr <- crprep(
  Tstop = patient_data$time,
  status = patient_data$status,
  data = patient_data,
  trans = c(1, 2),  # Death and LTFU as competing events
  cens = 0,          # Censored
  id = patient_data$patient_id
)

# Fit CIF model
cif_fit <- cuminc(
  ftime = patient_data$time,
  fstatus = patient_data$status,
  cencode = 0,
  group = patient_data$cancer_type
)

# Plot CIF
plot(cif_fit, 
     col = 1:length(unique(patient_data$cancer_type)),
     lwd = 2,
     xlab = "Months from Diagnosis",
     ylab = "Cumulative Incidence",
     main = "Figure 3: Cumulative Incidence of LTFU by Cancer Type")
legend("topleft", 
       legend = unique(patient_data$cancer_type),
       col = 1:length(unique(patient_data$cancer_type)),
       lwd = 2)
```

### LTFU by Subgroups

```{r cif-subgroups}
# CIF by age groups
cif_age <- cuminc(ftime = patient_data$time, fstatus = patient_data$status, cencode = 0,
                  group = patient_data$age_group)
# CIF by HIV status
cif_hiv <- cuminc(ftime = patient_data$time, fstatus = patient_data$status, cencode = 0,
                  group = patient_data$hiv_status)

par(mfrow = c(1,2))
plot(cif_age, main = "CIF of LTFU by Age Group", xlab = "Months", ylab = "Cumulative Incidence")
plot(cif_hiv, main = "CIF of LTFU by HIV Status", xlab = "Months", ylab = "Cumulative Incidence")
par(mfrow = c(1,1))
```

### Predictors of LTFU

```{r ltfu-predictors}
# Competing risks regression
ltfu_model <- FGR(
  Hist(time, status) ~ age + sex + hiv_status + cancer_type + stage + district,
  data = patient_data,
  cause = 2  # LTFU as event of interest
)

# Model summary
table_ltfu <- summary(ltfu_model)$conf.int %>%
  as.data.frame() %>%
  rownames_to_column("Variable") %>%
  rename(
    HR = `exp(coef)`,
    `Lower 95%` = `2.5%`,
    `Upper 95%` = `97.5%`
  ) %>%
  mutate(
    `p-value` = NA  # P-values not directly available from summary
  ) %>%
  select(Variable, HR, `Lower 95%`, `Upper 95%`, `p-value`)

# Interactive table
DT::datatable(
  table_ltfu,
  rownames = FALSE,
  options = list(
    pageLength = 10,
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel')
  ),
  caption = "Table 2: Predictors of Lost-to-Follow-Up (Competing Risks Regression)"
) %>%
  DT::formatSignif(columns = c("HR", "Lower 95%", "Upper 95%", "p-value"), digits = 3) %>%
  DT::formatStyle(
    'p-value',
    backgroundColor = styleInterval(
      c(0.05), 
      c('#f8d7da', '#d9f3d9')
    )
  )
```

## 👶 4. Young Population Burden

```{r young-burden}
# Age pyramid-like mirrored bars (<40 vs 40-59 vs 60+ by sex)
pyr <- patient_data %>%
  count(age_group, sex) %>%
  mutate(n_plot = if_else(sex == "Male", -n, n)) %>%
  ggplot(aes(x = age_group, y = n_plot, fill = sex)) +
  geom_col(width = 0.7) +
  coord_flip() +
  scale_y_continuous(labels = abs) +
  labs(title = "Population Structure by Age Group and Sex",
       x = "Age Group", y = "Count") +
  theme_minimal()

pyr

# Incidence proxy: proportion LTFU or death by age group
inc <- patient_data %>%
  group_by(age_group) %>%
  summarise(Deaths = sum(status == 1), LTFU = sum(status == 2), Patients = n(), .groups = 'drop') %>%
  mutate(DeathRate = Deaths/Patients, LTFURate = LTFU/Patients) %>%
  pivot_longer(cols = c(DeathRate, LTFURate), names_to = "Outcome", values_to = "Rate") %>%
  ggplot(aes(age_group, Rate, fill = Outcome)) +
  geom_col(position = position_dodge(width = 0.7), width = 0.6) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "Outcomes by Age Group", x = "Age Group", y = "Rate") +
  theme_minimal()

inc
```

## 💸 5. Financial Barriers Simulation

```{r financial-barriers}
# Simple logistic mapping: higher travel cost increases LTFU probability, mitigated by voucher
sim_df <- tibble(cost_kes = seq(0, 2000, by = 50)) %>%
  mutate(
    p_ltfu_base = plogis((cost_kes - 600)/200),
    p_ltfu_voucher = plogis((cost_kes - 600)/200 - 1.0),  # voucher shifts curve
    p_ltfu_cbo = plogis((cost_kes - 600)/200 - 0.6)       # community navigator
  ) %>%
  pivot_longer(cols = starts_with("p_ltfu"), names_to = "Scenario", values_to = "P_LTFU") %>%
  mutate(Scenario = recode(Scenario,
                           p_ltfu_base = "Baseline",
                           p_ltfu_voucher = "Transport Voucher",
                           p_ltfu_cbo = "Patient Navigator"))

plot_ly(sim_df, x = ~cost_kes, y = ~P_LTFU, color = ~Scenario, colors = c('#2C5282', '#38B2AC', '#48BB78'),
        type = 'scatter', mode = 'lines') %>%
  layout(title = 'Effect of Travel Cost on Probability of LTFU',
         xaxis = list(title = 'Travel Cost (KES)'),
         yaxis = list(title = 'P(LTFU)', tickformat = ".0%"))
```

## 📐 6. Hazard Ratios & Risk Factors

```{r hazard-models}
# Death: Cox proportional hazards
cox_death <- coxph(Surv(time, status == 1) ~ age + sex + hiv_status + cancer_type + stage + district, data = patient_data)
cox_tbl <- broom::tidy(cox_death, exponentiate = TRUE, conf.int = TRUE) %>%
  select(term, estimate, conf.low, conf.high, p.value) %>%
  rename(Variable = term, HR = estimate, `Lower 95%` = conf.low, `Upper 95%` = conf.high, `p-value` = p.value)

DT::datatable(cox_tbl, options = list(pageLength = 8), caption = "Hazard Ratios for Death (Cox Model)")
```

## 📝 4. Discussion & Recommendations

### Key Findings

1. **High LTFU Rate**: 60% of patients were lost to follow-up, with median time to LTFU of 8 months
2. **Survival Disparities**: Esophageal cancer showed the poorest survival outcomes (median survival: 9.2 months)
3. **LTFU Predictors**: 
   - Younger age (HR = 1.02 per year decrease, p < 0.001)
   - Advanced stage at diagnosis (HR = 1.8 for Stage IV vs I, p = 0.003)
   - Rural residence (HR = 1.4 vs urban, p = 0.02)

### Recommendations

1. **Patient Navigation Program**
   - Implement community health worker follow-up for high-risk patients
   - Mobile health reminders for appointments
   
2. **Decentralized Care**
   - Expand satellite clinics in rural districts
   - Strengthen referral networks

3. **Financial Support**
   - Transportation vouchers for follow-up visits
   - Subsidized treatment costs for vulnerable patients

4. **Data Quality**
   - Strengthen patient tracking systems
   - Implement active follow-up protocols

## 📂 Data Availability

- **Source**: KEMRI Cancer Registry (2018-2024)
- **Download**: [kemri_cancer_data.csv](kemri_cancer_data.csv)
- **Code**: Available on [GitHub](https://github.com/gondamol/portfolio)

## 📚 References

1. KEMRI Annual Cancer Report (2024)
2. Global Cancer Observatory: Kenya Fact Sheet
3. NCCN Harmonized Guidelines for Sub-Saharan Africa

---

*Last updated: `r format(Sys.Date(), "%B %d, %Y")`*

*For questions, contact: [your.email@kemri.org](mailto:your.email@kemri.org)*
