CONFIDENTIAL CASE FILE

Crime Report

Motivation

Breast cancer remains one of the most prevalent cancers worldwide, accounting for a substantial proportion of cancer incidence and health burden among women. Despite significant advances in early detection and treatment modalities, disparities in survival outcomes persist across demographic and biological groups. These disparities raise critical questions about equity in healthcare delivery and the factors that drive differential prognosis.

Clinical outcomes in breast cancer are influenced by multiple interacting factors:

Tumor characteristics: Size, stage, grade, and differentiation
Biological markers: Hormone receptor status (estrogen, progesterone)
Lymph node involvement: Regional spread and extent of nodal metastasis
Patient demographics: Age, race, and socioeconomic factors

Understanding how these factors jointly influence survival is essential for:

Improving risk stratification to identify high-risk patients
Guiding treatment decisions based on individualized prognosis
Identifying structural inequities that may be addressed through policy interventions
Building fair predictive models that perform equitably across populations

This project takes a creative “prosecution” approach—treating breast cancer as a defendant on trial—to engage audiences while rigorously analyzing survival patterns and predictors. Our goal is to determine whether the evidence is strong enough to “convict” breast cancer of systematic harm, and critically, whether that harm falls equally across all communities.

Related Work

Our analysis draws on established methodological frameworks and prior research in breast cancer epidemiology:

Methodological Foundations:

Cox Proportional Hazards Regression (Cox, 1972): The gold-standard approach for modeling time-to-event data with multiple covariates, allowing estimation of hazard ratios while accounting for censoring.
Kaplan-Meier Survival Analysis (Kaplan & Meier, 1958): Non-parametric estimation of survival functions, enabling visualization of survival trajectories across groups.
AJCC Cancer Staging System (6th Edition): The standardized TNM classification system used to categorize tumor extent (T), nodal involvement (N), and metastasis status.

Epidemiological Context:

DeSantis et al. (2019) documented persistent racial disparities in breast cancer mortality, with Black women experiencing 40% higher death rates than White women despite similar incidence rates.
Obermeyer et al. (2019) demonstrated how algorithmic bias in healthcare prediction can perpetuate inequities, motivating our fairness analysis.
Newman & Kaljee (2017) explored the intersection of biology and social determinants in explaining triple-negative breast cancer disparities among African American women.

Data Source:

This project uses the SEER (Surveillance, Epidemiology, and End Results) Breast Cancer Dataset curated by Reihaneh Namdari. The dataset contains 4,024 female breast cancer patients diagnosed with infiltrating duct and lobular carcinoma between 2006-2010, with follow-up for survival outcomes.

Initial Questions

Our investigation began with the following primary questions:

What clinical and demographic factors are associated with breast cancer mortality?
- Which tumor characteristics most strongly predict death?
- Does patient age have a linear or non-linear relationship with survival?
Do survival outcomes differ significantly across racial groups?
- Is there evidence of racial disparities in breast cancer survival?
- If so, what is the magnitude and statistical significance of these differences?
Can we build accurate predictive models for mortality risk?
- How well can Cox regression and logistic regression discriminate survivors from non-survivors?
- Are these models robust under cross-validation and bootstrap resampling?

Evolution of Questions

As our analysis progressed, new questions emerged:

Are our predictive models algorithmically fair?
- Does the model perform equally well for all racial groups?
- If disparities exist, are they attributable to model bias or underlying structural factors?
Do the data satisfy modeling assumptions?
- Does the Cox model satisfy the proportional hazards assumption?
- If not, what remediation strategies are appropriate?

Data

Data Source and Description

# Load and inspect the raw data
cancer_raw <- read_csv("data/cancer_data.csv")

# Display structure
cat("Dataset Dimensions:", nrow(cancer_raw), "observations x", ncol(cancer_raw), "variables\n\n")

## Dataset Dimensions: 4024 observations x 16 variables

The dataset contains 4024 breast cancer patients from the SEER registry. Key variables include Age, Race, Marital Status, TNM Stage (T-Stage, N-Stage), Tumor Differentiation, Estrogen/Progesterone Status, Survival Months, and Vital Status.

Data Cleaning and Preprocessing

# Clean and prepare data
cancer <- cancer_raw %>%
  janitor::clean_names() %>%
  mutate(
    # Convert categorical variables to factors with meaningful ordering
    race = factor(race, levels = c("White", "Black", "Other")),
    marital_status = factor(marital_status, 
                           levels = c("Single", "Married", "Separated", "Divorced", "Widowed")),
    t_stage = factor(t_stage, levels = c("T1", "T2", "T3", "T4")),
    n_stage = factor(n_stage, levels = c("N1", "N2", "N3")),
    x6th_stage = as.factor(x6th_stage),
    a_stage = factor(a_stage, levels = c("Regional", "Distant")),
    differentiate = factor(differentiate, 
                          levels = c("Well differentiated", "Moderately differentiated", 
                                    "Poorly differentiated", "Undifferentiated")),
    estrogen_status = factor(estrogen_status, levels = c("Positive", "Negative")),
    progesterone_status = factor(progesterone_status, levels = c("Positive", "Negative")),
    status = factor(status, levels = c("Alive", "Dead")),
    
    # Create binary outcome for modeling
    status_binary = ifelse(status == "Dead", 1, 0)
  ) %>%
  # Remove grade (redundant with differentiate - nearly perfect correlation)
  select(-grade)

# Check for missing data
missing_summary <- colSums(is.na(cancer))
cat("Missing Data Check:\n")

## Missing Data Check:

if(sum(missing_summary) == 0) {
  cat("✓ No missing values detected. Dataset is complete.\n")
} else {
  print(missing_summary[missing_summary > 0])
}

## ✓ No missing values detected. Dataset is complete.

Preprocessing Rationale:

Categorical variables were converted to factors with clinically meaningful level ordering (e.g., T1 < T2 < T3 < T4) to ensure proper interpretation in regression models. The grade variable was removed due to multicollinearity with differentiate. A binary outcome variable was created for logistic regression.

Descriptive Statistics

# Overall summary statistics
desc_table <- cancer %>%
  select(
    age, race, marital_status, t_stage, n_stage, differentiate,
    estrogen_status, progesterone_status, tumor_size,
    regional_node_examined, regional_node_positive,
    survival_months, status
  ) %>%
  tbl_summary(
    by = status,
    missing = "no",
    label = list(
      age ~ "Age (years)",
      race ~ "Race",
      marital_status ~ "Marital Status",
      t_stage ~ "Tumor Stage (T)",
      n_stage ~ "Nodal Stage (N)",
      differentiate ~ "Tumor Differentiation",
      estrogen_status ~ "Estrogen Receptor",
      progesterone_status ~ "Progesterone Receptor",
      tumor_size ~ "Tumor Size (mm)",
      regional_node_examined ~ "Nodes Examined",
      regional_node_positive ~ "Nodes Positive",
      survival_months ~ "Survival Time (months)"
    ),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    )
  ) %>%
  add_overall() %>%
  add_p() %>%
  bold_labels()

desc_table %>%
  as_gt() %>%
  tab_options(
    table.background.color = "#EFEDE2",
    heading.background.color = "#EFEDE2",
    column_labels.background.color = "#EFEDE2",
    table.border.top.color = "#d2c9a4",
    table.border.bottom.color = "#d2c9a4",
    column_labels.border.top.color = "#d2c9a4",
    column_labels.border.bottom.color = "#d2c9a4",
    data_row.padding = px(6),
    table.font.names = "Courier New"
  ) %>%
  cols_align(
    align = "left",
    columns = everything()
  ) %>%
  tab_style(
    style = cell_borders(color = "#d2c9a4", sides = "all"),
    locations = cells_body()
  ) %>%
  tab_style(
    style = cell_borders(color = "#d2c9a4", sides = "all"),
    locations = cells_column_labels()
  )

Characteristic	Overall N = 4,024¹	Alive N = 3,408¹	Dead N = 616¹	p-value²
Age (years)	54 (9)	54 (9)	55 (10)	<0.001
Race				<0.001
White	3,413 (85%)	2,903 (85%)	510 (83%)
Black	291 (7.2%)	218 (6.4%)	73 (12%)
Other	320 (8.0%)	287 (8.4%)	33 (5.4%)
Marital Status				<0.001
Single	615 (15%)	511 (15%)	104 (17%)
Married	2,643 (66%)	2,285 (67%)	358 (58%)
Separated	45 (1.1%)	30 (0.9%)	15 (2.4%)
Divorced	486 (12%)	396 (12%)	90 (15%)
Widowed	235 (5.8%)	186 (5.5%)	49 (8.0%)
Tumor Stage (T)				<0.001
T1	1,603 (40%)	1,446 (42%)	157 (25%)
T2	1,786 (44%)	1,483 (44%)	303 (49%)
T3	533 (13%)	417 (12%)	116 (19%)
T4	102 (2.5%)	62 (1.8%)	40 (6.5%)
Nodal Stage (N)				<0.001
N1	2,732 (68%)	2,462 (72%)	270 (44%)
N2	820 (20%)	655 (19%)	165 (27%)
N3	472 (12%)	291 (8.5%)	181 (29%)
Tumor Differentiation				<0.001
Well differentiated	543 (13%)	504 (15%)	39 (6.3%)
Moderately differentiated	2,351 (58%)	2,046 (60%)	305 (50%)
Poorly differentiated	1,111 (28%)	848 (25%)	263 (43%)
Undifferentiated	19 (0.5%)	10 (0.3%)	9 (1.5%)
Estrogen Receptor				<0.001
Positive	3,755 (93%)	3,247 (95%)	508 (82%)
Negative	269 (6.7%)	161 (4.7%)	108 (18%)
Progesterone Receptor				<0.001
Positive	3,326 (83%)	2,914 (86%)	412 (67%)
Negative	698 (17%)	494 (14%)	204 (33%)
Tumor Size (mm)	30 (21)	29 (20)	37 (24)	<0.001
Nodes Examined	14 (8)	14 (8)	15 (8)	0.067
Nodes Positive	4.2 (5.1)	3.6 (4.4)	7.2 (7.3)	<0.001
Survival Time (months)	71 (23)	76 (19)	46 (24)	<0.001
¹ Mean (SD); n (%)
² Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

Interpretation of Descriptive Statistics

The descriptive analysis reveals significant baseline differences between survivors and non-survivors. The cohort is predominantly White (approximately 81%), with Black patients comprising about 16% and Other races about 3%.

Key differences between groups (p < 0.001):

Tumor Characteristics: Non-survivors presented with more advanced disease at diagnosis. A higher proportion of deceased patients had T3-T4 tumors compared to survivors. Similarly, deceased patients showed higher rates of N2-N3 nodal involvement compared to survivors.
Biological Markers: Estrogen receptor (ER) positivity was significantly higher in survivors compared to non-survivors, reflecting the availability of targeted hormonal therapies for ER+ disease.
Tumor Size: Mean tumor size was markedly larger in the deceased group compared to survivors, validating tumor size as a crude but effective proxy for disease progression.

Exploratory Analysis

Cohort Overview

# Key statistics
total_n <- nrow(cancer)
deaths <- sum(cancer$status == "Dead")
survivors <- sum(cancer$status == "Alive")
fatality_rate <- round(deaths / total_n * 100, 1)
median_survival <- median(cancer$survival_months)
mean_age <- round(mean(cancer$age), 1)
age_iqr <- quantile(cancer$age, c(0.25, 0.75))

cohort_summary <- tibble(
  Metric = c("Total Patients", "Deaths", "Survivors", "Fatality Rate", 
             "Median Survival", "Mean Age"),
  Value = c(format(total_n, big.mark = ","), 
            deaths, 
            survivors, 
            paste0(fatality_rate, "%"),
            paste0(median_survival, " months"),
            paste0(mean_age, " years"))
)

kable(cohort_summary, caption = "Table 1: Cohort Summary Statistics")

Table 1: Cohort Summary Statistics
Metric	Value
Total Patients	4,024
Deaths	616
Survivors	3408
Fatality Rate	15.3%
Median Survival	73 months
Mean Age	54 years

Interpretation

The overall fatality rate in this cohort is 15.3%. The median survival time of 73 months (approximately 6.1 years) indicates a relatively long follow-up period, suitable for survival analysis. This fatality rate is consistent with contemporary breast cancer statistics, reflecting improvements in treatment while acknowledging that breast cancer remains a significant cause of mortality.

Age Distribution by Outcome

# Age distribution
age_plot <- ggplot(cancer, aes(x = age, fill = status)) +
  geom_histogram(binwidth = 5, alpha = 0.7, color = "white", position = "identity") +
  scale_fill_manual(values = c("Alive" = "#2E9FDF", "Dead" = "#E74C3C"),
                   labels = c("Survived", "Deceased")) +
  geom_vline(xintercept = quantile(cancer$age, c(0.25, 0.75)), 
             linetype = "dashed", color = "gray40") +
  labs(
    title = "Figure 1: Age Distribution of Breast Cancer Patients by Vital Status",
    subtitle = paste0("Dashed lines indicate IQR: ", age_iqr[1], " - ", age_iqr[2], " years"),
    x = "Age at Diagnosis (years)",
    y = "Number of Patients",
    fill = "Outcome"
  ) +
  theme_minimal(paper ="#EFEDE2") +
  theme(legend.position = "bottom")

age_plot

Interpretation

Figure 1 approximates a normal distribution centered around 50-55 years. Notably, the distribution of deceased patients (red) appears slightly right-shifted compared to survivors (blue), suggesting that advanced age is associated with higher mortality. The overlap between distributions, however, indicates that age alone is not a definitive predictor, motivating the need for multivariate modeling. The interquartile range of 47-61 years reflects the typical age range for breast cancer diagnosis.

Racial Disparities in Outcomes

# Racial breakdown
racial_summary <- cancer %>%
  group_by(race) %>%
  summarise(
    N = n(),
    `% of Cohort` = round(n() / nrow(cancer) * 100, 1),
    Deaths = sum(status == "Dead"),
    `Fatality Rate (%)` = round(Deaths / n() * 100, 1),
    .groups = 'drop'
  )

kable(racial_summary, caption = "Table 2: Patient Characteristics by Race")

Table 2: Patient Characteristics by Race
race	N	% of Cohort	Deaths	Fatality Rate (%)
White	3413	84.8	510	14.9
Black	291	7.2	73	25.1
Other	320	8.0	33	10.3

ggplot(racial_summary, aes(x = race, y = `Fatality Rate (%)`, fill = race)) +
  geom_col(alpha = 0.8, color = "black", width = 0.7) +
  geom_text(aes(label = paste0(`Fatality Rate (%)`, "%")), vjust = -0.5, fontface = "bold", size = 5) +
  scale_fill_manual(values = c("White" = "#E7B800", "Black" = "#2E9FDF", "Other" = "#8E44AD")) +
  labs(
    title = "Figure 2: Mortality Rate by Race",
    subtitle = "Black patients show disproportionately higher mortality",
    x = "Race",
    y = "Fatality Rate (%)"
  ) +
  theme_minimal(paper ="#EFEDE2") +
  theme(legend.position = "none") +
  ylim(0, 30)

Interpretation

Figure 2 and Table 2 highlight a striking racial disparity. Black patients have a fatality rate of **

%, which is approximately 1.7 times higher** than that of White patients (

%).

This disparity aligns with national SEER data trends and suggests systemic issues affecting outcomes for Black women, potentially including:

Later stage at diagnosis due to screening access barriers
Tumor biology differences (e.g., higher rates of triple-negative breast cancer)
Disparities in treatment access and quality
Socioeconomic factors not captured in this dataset

Tumor Stage Distribution

# Calculate fatality rates by stage for annotation
t_stage_rates <- cancer %>%
  group_by(t_stage) %>%
  summarise(
    fatality = round(sum(status == "Dead") / n() * 100, 1),
    .groups = 'drop'
  )

n_stage_rates <- cancer %>%
  group_by(n_stage) %>%
  summarise(
    fatality = round(sum(status == "Dead") / n() * 100, 1),
    .groups = 'drop'
  )

p1 <- cancer %>%
  count(t_stage, status) %>%
  ggplot(aes(x = t_stage, y = n, fill = status)) +
  geom_col(position = "dodge", color = "black", alpha = 0.8) +
  scale_fill_manual(values = c("Alive" = "#2E9FDF", "Dead" = "#E74C3C"),
                   labels = c("Survived", "Deceased")) +
  labs(title = "Figure 3A: Distribution by Tumor Stage (T)", 
       x = "T Stage", y = "Count", fill = "Status") +
  theme_minimal(paper ="#EFEDE2") + 
  theme(legend.position = "none")

p2 <- cancer %>%
  count(n_stage, status) %>%
  ggplot(aes(x = n_stage, y = n, fill = status)) +
  geom_col(position = "dodge", color = "black", alpha = 0.8) +
  scale_fill_manual(values = c("Alive" = "#2E9FDF", "Dead" = "#E74C3C"), 
                    labels = c("Survived", "Deceased")) +
  labs(title = "Figure 3B: Distribution by Nodal Stage (N)", 
       x = "N Stage", y = "Count", fill = "Status") +
  theme_minimal(paper ="#EFEDE2") + 
  theme(legend.position = "bottom")

grid.arrange(p1, p2, ncol = 2)

# Combined fatality by stage
stage_fatality <- bind_rows(
  t_stage_rates %>% mutate(Type = "Tumor (T)", Stage = as.character(t_stage)) %>% 
    select(Type, Stage, `Fatality Rate (%)` = fatality),
  n_stage_rates %>% mutate(Type = "Nodal (N)", Stage = as.character(n_stage)) %>% 
    select(Type, Stage, `Fatality Rate (%)` = fatality)
)

kable(stage_fatality, caption = "Table 3: Fatality Rate by Cancer Stage")

Table 3: Fatality Rate by Cancer Stage
Type	Stage	Fatality Rate (%)
Tumor (T)	T1	9.8
Tumor (T)	T2	17.0
Tumor (T)	T3	21.8
Tumor (T)	T4	39.2
Nodal (N)	N1	9.9
Nodal (N)	N2	20.1
Nodal (N)	N3	38.3

Interpretation

Figures 3A-B and Table 3 demonstrate a clear “dose-response” relationship between cancer stage and mortality:

Tumor Stage (T): The fatality rate increases progressively from T1 (9.8%) to T4 (39.2%), representing a 4-fold increase.
Nodal Stage (N): Nodal involvement is a critical prognostic indicator. The fatality rate increases from N1 (9.9%) to N3 (38.3%), confirming that early detection before nodal spread is paramount for survival.

This validates the biological basis of the TNM staging system and emphasizes the importance of screening programs for early detection.

Correlation Analysis

# Correlation matrix for continuous variables
cor_matrix <- cancer %>%
  select(age, tumor_size, regional_node_examined, regional_node_positive, survival_months, status_binary) %>%
  cor()
par(bg = "#EFEDE2")

corrplot(
  cor_matrix,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  number.cex = 0.7,
  tl.col = "black",
  tl.srt = 45,
  title = "Figure 4: Correlation Matrix of Continuous Variables",
  mar = c(0, 0, 2, 0),
  col = colorRampPalette(c("#E74C3C", "white", "#2E9FDF"))(100)
)

Interpretation

The correlation matrix (Figure 4) reveals important relationships:

Nodes & Mortality: regional_node_positive has a moderate positive correlation with status_binary (r = 0.26), confirming that nodal metastasis is a strong driver of mortality.
Size & Nodes: There is a moderate positive correlation between tumor_size and regional_node_positive (r = 0.24), consistent with the biological understanding that larger tumors are more likely to metastasize.
Survival & Mortality: As expected, survival_months is negatively correlated with status_binary (r = -0.48), as death events terminate the survival interval.
Multicollinearity: No predictor correlations exceed |0.7|, suggesting multicollinearity will not substantially affect our regression models.

Additional Analysis

Kaplan-Meier Survival Analysis

Kaplan-Meier analysis provides non-parametric estimates of survival probability over time, accounting for censored observations.

Overall Survival

# Overall KM curve
survfit_overall <- survfit(Surv(survival_months, status_binary) ~ 1, data = cancer)

# Calculate 5-year survival
surv_5yr <- summary(survfit_overall, times = 60)$surv

ggsurvplot(
  survfit_overall,
  data = cancer,
  conf.int = TRUE,
  risk.table = TRUE,
  surv.median.line = "hv",
  title = "Figure 5: Overall Kaplan-Meier Survival Curve",
  subtitle = paste0("5-year survival: ", round(surv_5yr * 100, 1), "%"),
  xlab = "Time Since Diagnosis (months)",
  ylab = "Survival Probability",
  palette = "#2C3E50",
  ggtheme = theme_minimal(paper ="#EFEDE2")
)

Interpretation

Figure 5 shows the overall survival trajectory for the cohort. The 5-year survival rate is approximately 88.1%, reflecting the generally favorable prognosis for breast cancer with modern treatment. The survival curve shows a gradual decline over time, with the steepest drops occurring in the first 60 months after diagnosis.

Survival by Race

# KM by race
survfit_race <- survfit(Surv(survival_months, status_binary) ~ race, data = cancer)
logrank_race <- survdiff(Surv(survival_months, status_binary) ~ race, data = cancer)
logrank_p <- 1 - pchisq(logrank_race$chisq, df = 2)

# 5-year survival by race
surv_5yr_race <- summary(survfit_race, times = 60)

ggsurvplot(
  survfit_race,
  data = cancer,
  conf.int = TRUE,
  risk.table = TRUE,
  pval = TRUE,
  title = "Figure 6: Kaplan-Meier Survival Curves by Race",
  subtitle = "Significant racial disparities in survival (log-rank p < 0.0001)",
  xlab = "Time Since Diagnosis (months)",
  ylab = "Survival Probability",
  palette = c("#E7B800", "#2E9FDF", "#8E44AD"),
  legend.labs = c("White", "Black", "Other"),
  ggtheme = theme_minimal(paper ="#EFEDE2")
)

Interpretation

The Kaplan-Meier curves (Figure 6) provide statistical validation of the racial disparities observed in the exploratory analysis:

Statistical Significance: The log-rank test p-value (< 0.0001) confirms that the survival distributions are significantly different across racial groups.
Persistent Disparity: The curve for Black patients (blue) diverges early and remains consistently lower than White patients (yellow) throughout follow-up. This indicates a persistent survival disadvantage that is not limited to a specific timepoint in the disease trajectory.
5-Year Survival Estimates:
- White: 88.4%
- Black: 79.5%
- Other: 93%
Clinical Implications: The 8.9 percentage point gap between White and Black patients represents a substantial disparity demanding intervention.

Survival by Tumor Stage

# KM by T stage
survfit_tstage <- survfit(Surv(survival_months, status_binary) ~ t_stage, data = cancer)

ggsurvplot(
  survfit_tstage,
  data = cancer,
  conf.int = FALSE,
  risk.table = TRUE,
  pval = TRUE,
  title = "Figure 7: Kaplan-Meier Survival Curves by Tumor Stage (T)",
  subtitle = "Clear dose-response relationship: survival decreases with advancing stage",
  xlab = "Time Since Diagnosis (months)",
  ylab = "Survival Probability",
  palette = "Set1",
  legend.labs = c("T1", "T2", "T3", "T4"),
  ggtheme = theme_minimal(paper ="#EFEDE2")
)

Interpretation

Figure 7 illustrates distinct risk stratification by tumor stage with several key findings:

Perfect Ordinal Separation: The survival curves show perfect separation by stage (T1 > T2 > T3 > T4), with no crossing. This validates the biological basis of the TNM staging system.
Magnitude of Effect: The 5-year survival probability drops dramatically from T1 tumors (>90%) to T4 tumors (<70%), representing a clinically meaningful difference of over 20 percentage points.
Early Detection Imperative: These results provide strong evidence-based justification for screening mammography—detecting tumors at T1 rather than T4 could substantially improve survival outcomes.

Survival by Nodal Stage

# KM by N stage
survfit_nstage <- survfit(Surv(survival_months, status_binary) ~ n_stage, data = cancer)

ggsurvplot(
  survfit_nstage,
  data = cancer,
  conf.int = FALSE,
  risk.table = TRUE,
  pval = TRUE,
  title = "Figure 8: Kaplan-Meier Survival Curves by Nodal Stage (N)",
  subtitle = "Greater nodal involvement associated with worse prognosis",
  xlab = "Time Since Diagnosis (months)",
  ylab = "Survival Probability",
  palette = "Dark2",
  legend.labs = c("N1", "N2", "N3"),
  ggtheme = theme_minimal(paper ="#EFEDE2")
)

Interpretation

Figure 8 demonstrates the critical prognostic importance of lymph node involvement:

Progressive Decline: Survival probability decreases progressively with increasing nodal involvement (N1 > N2 > N3).
Biological Rationale: Lymph node metastasis indicates that cancer has begun to spread beyond the primary tumor site. Extensive nodal involvement (N3) suggests widespread regional dissemination and higher likelihood of distant metastasis.
Treatment Implications: Nodal status guides treatment decisions—patients with N2-N3 disease typically receive more aggressive adjuvant chemotherapy and radiation to address the higher recurrence risk.

Cox Proportional Hazards Regression

Cox regression allows estimation of the independent effect of each predictor on survival, controlling for confounding.

Proportional Hazards Assumption Test

# Define formula for initial model
cox_formula <- as.formula(
  Surv(survival_months, status_binary) ~ 
    age + race + marital_status + t_stage + n_stage + 
    differentiate + tumor_size + estrogen_status + progesterone_status + 
    regional_node_examined + regional_node_positive
)

# Fit initial model
cox_model <- coxph(cox_formula, data = cancer)

# Test PH assumption
ph_test <- cox.zph(cox_model)
print(ph_test)

##                          chisq df       p
## age                     0.1051  1 0.74584
## race                    0.9508  2 0.62165
## marital_status          2.9355  4 0.56867
## t_stage                 0.2729  3 0.96505
## n_stage                 1.6048  2 0.44826
## differentiate           1.8576  3 0.60248
## tumor_size              0.9924  1 0.31915
## estrogen_status        29.3770  1 6.0e-08
## progesterone_status    32.3417  1 1.3e-08
## regional_node_examined  0.0104  1 0.91881
## regional_node_positive  0.0164  1 0.89818
## GLOBAL                 51.1608 20 0.00015

# Identify violations
violations <- rownames(ph_test$table)[ph_test$table[, "p"] < 0.05]
if(length(violations) > 0) {
  cat("\n⚠ Variables with PH violations (p < 0.05):", paste(violations, collapse = ", "), "\n")
  cat("Remedy: Stratification by violating variables.\n")
}

## 
## ⚠ Variables with PH violations (p < 0.05): estrogen_status, progesterone_status, GLOBAL 
## Remedy: Stratification by violating variables.

Interpretation

The proportional hazards test reveals violations for estrogen_status and progesterone_status (p < 0.05). This means the effect of hormone receptor status on survival changes over time—a biologically plausible finding since hormonal therapies have different efficacy profiles at different disease stages. We address this by stratifying the model by hormone receptor status.

Stratified Cox Model

# Stratified model to address PH violations
cox_stratified <- coxph(
  Surv(survival_months, status_binary) ~ 
    age + race + marital_status + t_stage + n_stage + 
    differentiate + tumor_size + regional_node_examined + regional_node_positive +
    strata(estrogen_status, progesterone_status),
  data = cancer
)

# Verify PH assumption
ph_test_strat <- cox.zph(cox_stratified)
global_p <- ph_test_strat$table["GLOBAL", "p"]

cat("Stratified Model - Global PH Test p-value:", round(global_p, 4), "\n")

## Stratified Model - Global PH Test p-value: 0.8413

if(global_p > 0.05) {
  cat("✓ Proportional hazards assumption now satisfied.\n")
}

## ✓ Proportional hazards assumption now satisfied.

# Model summary
summary(cox_stratified)

## Call:
## coxph(formula = Surv(survival_months, status_binary) ~ age + 
##     race + marital_status + t_stage + n_stage + differentiate + 
##     tumor_size + regional_node_examined + regional_node_positive + 
##     strata(estrogen_status, progesterone_status), data = cancer)
## 
##   n= 4024, number of events= 616 
## 
##                                             coef exp(coef)  se(coef)      z
## age                                     0.019810  1.020007  0.004862  4.074
## raceBlack                               0.373448  1.452735  0.129523  2.883
## raceOther                              -0.366107  0.693429  0.180670 -2.026
## marital_statusMarried                  -0.160807  0.851457  0.113294 -1.419
## marital_statusSeparated                 0.322121  1.380052  0.289343  1.113
## marital_statusDivorced                  0.029555  1.029996  0.146015  0.202
## marital_statusWidowed                  -0.011422  0.988643  0.179705 -0.064
## t_stageT2                               0.351044  1.420550  0.113835  3.084
## t_stageT3                               0.419737  1.521561  0.219809  1.910
## t_stageT4                               0.828155  2.289092  0.233050  3.554
## n_stageN2                               0.453907  1.574452  0.109757  4.136
## n_stageN3                               0.622174  1.862973  0.179063  3.475
## differentiateModerately differentiated  0.450487  1.569077  0.171091  2.633
## differentiatePoorly differentiated      0.795512  2.215575  0.176940  4.496
## differentiateUndifferentiated           1.502353  4.492249  0.378095  3.973
## tumor_size                              0.001439  1.001440  0.003157  0.456
## regional_node_examined                 -0.032959  0.967578  0.006494 -5.075
## regional_node_positive                  0.057181  1.058848  0.011375  5.027
##                                        Pr(>|z|)    
## age                                    4.62e-05 ***
## raceBlack                              0.003936 ** 
## raceOther                              0.042726 *  
## marital_statusMarried                  0.155790    
## marital_statusSeparated                0.265586    
## marital_statusDivorced                 0.839596    
## marital_statusWidowed                  0.949320    
## t_stageT2                              0.002044 ** 
## t_stageT3                              0.056191 .  
## t_stageT4                              0.000380 ***
## n_stageN2                              3.54e-05 ***
## n_stageN3                              0.000512 ***
## differentiateModerately differentiated 0.008463 ** 
## differentiatePoorly differentiated     6.93e-06 ***
## differentiateUndifferentiated          7.08e-05 ***
## tumor_size                             0.648617    
## regional_node_examined                 3.87e-07 ***
## regional_node_positive                 4.98e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                        exp(coef) exp(-coef) lower .95 upper .95
## age                                       1.0200     0.9804    1.0103    1.0298
## raceBlack                                 1.4527     0.6884    1.1270    1.8726
## raceOther                                 0.6934     1.4421    0.4866    0.9881
## marital_statusMarried                     0.8515     1.1745    0.6819    1.0632
## marital_statusSeparated                   1.3801     0.7246    0.7827    2.4332
## marital_statusDivorced                    1.0300     0.9709    0.7737    1.3713
## marital_statusWidowed                     0.9886     1.0115    0.6951    1.4061
## t_stageT2                                 1.4205     0.7040    1.1365    1.7756
## t_stageT3                                 1.5216     0.6572    0.9890    2.3409
## t_stageT4                                 2.2891     0.4369    1.4497    3.6144
## n_stageN2                                 1.5745     0.6351    1.2697    1.9523
## n_stageN3                                 1.8630     0.5368    1.3116    2.6462
## differentiateModerately differentiated    1.5691     0.6373    1.1220    2.1942
## differentiatePoorly differentiated        2.2156     0.4514    1.5663    3.1340
## differentiateUndifferentiated             4.4922     0.2226    2.1411    9.4254
## tumor_size                                1.0014     0.9986    0.9953    1.0077
## regional_node_examined                    0.9676     1.0335    0.9553    0.9800
## regional_node_positive                    1.0588     0.9444    1.0355    1.0827
## 
## Concordance= 0.709  (se = 0.013 )
## Likelihood ratio test= 340.9  on 18 df,   p=<2e-16
## Wald test            = 363.7  on 18 df,   p=<2e-16
## Score (logrank) test = 405.1  on 18 df,   p=<2e-16

Forest Plot of Hazard Ratios

# Tidy results
cox_tidy <- broom::tidy(cox_stratified, exponentiate = TRUE, conf.int = TRUE) %>%
  mutate(
    term_clean = case_when(
      str_detect(term, "^age") ~ "Age (per year)",
      str_detect(term, "raceBlack") ~ "Race: Black vs. White",
      str_detect(term, "raceOther") ~ "Race: Other vs. White",
      str_detect(term, "marital_status.*Divorced") ~ "Marital: Divorced vs. Single",
      str_detect(term, "marital_status.*Married") ~ "Marital: Married vs. Single",
      str_detect(term, "marital_status.*Separated") ~ "Marital: Separated vs. Single",
      str_detect(term, "marital_status.*Widowed") ~ "Marital: Widowed vs. Single",
      str_detect(term, "t_stageT2") ~ "Tumor: T2 vs. T1",
      str_detect(term, "t_stageT3") ~ "Tumor: T3 vs. T1",
      str_detect(term, "t_stageT4") ~ "Tumor: T4 vs. T1",
      str_detect(term, "n_stageN2") ~ "Nodal: N2 vs. N1",
      str_detect(term, "n_stageN3") ~ "Nodal: N3 vs. N1",
      str_detect(term, "differentiate.*Moderate") ~ "Diff: Moderate vs. Well",
      str_detect(term, "differentiate.*Poor") ~ "Diff: Poor vs. Well",
      str_detect(term, "differentiate.*Undiff") ~ "Diff: Undifferentiated vs. Well",
      str_detect(term, "tumor_size") ~ "Tumor Size (per mm)",
      str_detect(term, "regional_node_examined") ~ "Nodes Examined (per node)",
      str_detect(term, "regional_node_positive") ~ "Nodes Positive (per node)",
      TRUE ~ term
    ),
    significant = ifelse(p.value < 0.05, "p < 0.05", "p ≥ 0.05")
  )

# Forest plot
ggplot(cox_tidy, aes(x = estimate, y = reorder(term_clean, estimate))) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "red", size = 1) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high, color = significant), size = 0.7) +
  scale_color_manual(values = c("p < 0.05" = "#E74C3C", "p ≥ 0.05" = "gray50")) +
  scale_x_log10(breaks = c(0.5, 1, 1.5, 2, 2.5, 3)) +
  labs(
    title = "Figure 9: Forest Plot of Adjusted Hazard Ratios",
    subtitle = "Stratified Cox Model (stratified by ER/PR status)",
    x = "Hazard Ratio (log scale)",
    y = "",
    color = "Significance"
  ) +
  theme_minimal(paper ="#EFEDE2") +
  theme(legend.position = "bottom")

# Significant results table
sig_results <- cox_tidy %>%
  filter(p.value < 0.05) %>%
  arrange(desc(estimate)) %>%
  mutate(
    `HR (95% CI)` = paste0(round(estimate, 2), " (", round(conf.low, 2), "-", round(conf.high, 2), ")"),
    `p-value` = ifelse(p.value < 0.001, "<0.001", round(p.value, 3))
  ) %>%
  select(Variable = term_clean, `HR (95% CI)`, `p-value`)

kable(sig_results, caption = "Table 4: Significant Predictors from Cox Model")

Table 4: Significant Predictors from Cox Model
Variable	HR (95% CI)	p-value
Diff: Undifferentiated vs. Well	4.49 (2.14-9.43)	<0.001
Tumor: T4 vs. T1	2.29 (1.45-3.61)	<0.001
Diff: Poor vs. Well	2.22 (1.57-3.13)	<0.001
Nodal: N3 vs. N1	1.86 (1.31-2.65)	<0.001
Nodal: N2 vs. N1	1.57 (1.27-1.95)	<0.001
Diff: Moderate vs. Well	1.57 (1.12-2.19)	0.008
Race: Black vs. White	1.45 (1.13-1.87)	0.004
Tumor: T2 vs. T1	1.42 (1.14-1.78)	0.002
Nodes Positive (per node)	1.06 (1.04-1.08)	<0.001
Age (per year)	1.02 (1.01-1.03)	<0.001
Nodes Examined (per node)	0.97 (0.96-0.98)	<0.001
Race: Other vs. White	0.69 (0.49-0.99)	0.043

Interpretation

The Cox regression (Figure 9 and Table 4) identifies independent predictors of mortality after adjustment:

Significant Risk Factors (HR > 1):

Advanced Tumor Stage: T4 tumors have the highest hazard ratio, confirming tumor stage as the strongest risk factor. The progressive increase from T2 to T4 validates the TNM staging system.
Nodal Involvement: N2 and N3 stages show significantly elevated hazard compared to N1, with each additional positive lymph node incrementally increasing mortality risk.
Race: Even after adjusting for clinical factors like stage and tumor size, Black race remains a significant predictor of mortality (HR > 1). This suggests that the racial survival gap is not fully explained by stage at diagnosis or tumor characteristics, pointing to unmeasured confounders (e.g., treatment quality, socioeconomics).
Poor Differentiation: Poorly differentiated and undifferentiated tumors are associated with higher hazard, reflecting aggressive tumor biology.

Protective Factors (HR < 1):

Nodes Examined: More lymph nodes examined is associated with lower hazard, likely reflecting better surgical staging and more complete treatment.
Married Status: Married patients show lower hazard than single patients, consistent with literature on social support and cancer outcomes.

Model Performance

# Concordance index
concordance <- cox_stratified$concordance[1]

cat("Model Discrimination:\n")

## Model Discrimination:

cat("C-statistic (Concordance Index):", round(concordance, 4), "\n")

## C-statistic (Concordance Index): 762041

cat("Interpretation: Model correctly ranks", round(concordance * 100, 1), "% of patient pairs by risk.\n")

## Interpretation: Model correctly ranks 76204100 % of patient pairs by risk.

cat("(0.5 = random, 0.7-0.8 = acceptable, >0.8 = excellent)\n")

## (0.5 = random, 0.7-0.8 = acceptable, >0.8 = excellent)

Logistic Regression & Fairness Analysis

We utilized logistic regression to predict binary mortality status and assessed the model’s fairness across racial groups.

Model Fitting and Performance

# Train/Test Split
set.seed(123)
train_idx <- createDataPartition(cancer$status_binary, p = 0.7, list = FALSE)
train <- cancer[train_idx, ]
test <- cancer[-train_idx, ]

cat("Data Split:\n")

## Data Split:

cat("• Training:", nrow(train), "patients (70%)\n")

## • Training: 2817 patients (70%)

cat("• Testing:", nrow(test), "patients (30%)\n")

## • Testing: 1207 patients (30%)

# Model with polynomial age term
logit_model <- glm(
  status_binary ~ poly(age, 2) + race + marital_status + t_stage + 
    differentiate + estrogen_status + progesterone_status + 
    regional_node_examined + regional_node_positive,
  data = train,
  family = binomial
)

# Predictions
test$pred_prob <- predict(logit_model, newdata = test, type = "response")
test$pred_class <- ifelse(test$pred_prob > 0.5, 1, 0)

# Overall metrics
accuracy <- mean(test$pred_class == test$status_binary)
roc_obj <- roc(test$status_binary, test$pred_prob, quiet = TRUE)
auc_val <- auc(roc_obj)

cat("\nModel Performance:\n")

## 
## Model Performance:

cat("• Accuracy:", round(accuracy, 3), "\n")

## • Accuracy: 0.854

cat("• AUC-ROC:", round(as.numeric(auc_val), 3), "\n")

## • AUC-ROC: 0.777

ROC Curve

ggroc(roc_obj, color = "#E74C3C", size = 1.2) +
  geom_abline(intercept = 1, slope = 1, linetype = "dashed", color = "gray50") +
  annotate("text", x = 0.3, y = 0.3, 
           label = paste("AUC =", round(auc_val, 3)),
           size = 6, fontface = "bold") +
  labs(
    title = "Figure 10: ROC Curve - Logistic Regression Model",
    subtitle = "Model shows good discrimination between survivors and non-survivors",
    x = "Specificity",
    y = "Sensitivity"
  ) +
  theme_minimal(paper ="#EFEDE2")

Fairness Audit

# Fairness Audit by Race
fairness <- test %>%
  group_by(race) %>%
  summarise(
    N = n(),
    Deaths = sum(status_binary),
    `Mortality (%)` = round(Deaths / N * 100, 1),
    AUC = round(as.numeric(auc(roc(status_binary, pred_prob, quiet = TRUE))), 3),
    .groups = 'drop'
  )

kable(fairness, caption = "Table 5: Algorithmic Fairness Audit (Performance by Race)")

Table 5: Algorithmic Fairness Audit (Performance by Race)
race	N	Deaths	Mortality (%)	AUC
White	1026	157	15.3	0.774
Black	82	25	30.5	0.792
Other	99	9	9.1	0.623

ggplot(fairness, aes(x = race, y = AUC, fill = race)) +
  geom_col(alpha = 0.8, color = "black", width = 0.6) +
  geom_text(aes(label = AUC), vjust = -0.5, fontface = "bold", size = 5) +
  geom_hline(yintercept = as.numeric(auc_val), linetype = "dashed", size = 1) +
  scale_fill_manual(values = c("White" = "#E7B800", "Black" = "#2E9FDF", "Other" = "#8E44AD")) +
  ylim(0, 1) +
  labs(
    title = "Figure 11: Algorithmic Fairness - Model AUC by Race",
    subtitle = paste0("Dashed line = overall AUC (", round(auc_val, 3), ")"),
    x = "Race",
    y = "AUC"
  ) +
  theme_minimal(paper ="#EFEDE2") +
  theme(legend.position = "none")

Interpretation

The fairness audit (Table 5 and Figure 11) indicates that our predictive model achieves algorithmic parity:

Comparable AUC: The AUC values for White (0.774) and Black (0.792) patients are comparable, indicating similar discriminative ability across groups.
Key Implication: The model discriminates between survivors and non-survivors equally well for both racial groups. It is not “biased” in its predictions.
Crucial Distinction: While the model is mathematically fair, the outcomes remain disparate. The model correctly predicts that Black patients have higher risk, reflecting the real-world inequalities in the data rather than algorithmic error.
Clinical Significance: This finding suggests the model can be deployed equitably without exacerbating existing disparities through differential prediction quality.

Discussion

Summary of Key Findings

This forensic analysis of 4,024 breast cancer patients reveals consistent patterns that indict both biological aggression and structural inequities:

The “Guilty” Parties (Risk Factors)

Advanced Tumor (T) and Nodal (N) stages are the most potent predictors of mortality. Specifically, the transition from N1 to N3 represents a catastrophic shift in prognosis. The Cox model C-statistic of 7.62041^{5} indicates good discriminatory ability.

The Racial Gap

Black patients face significantly higher mortality (

% vs

%) and worse survival trajectories. Crucially, the Cox model shows this risk persists even when adjusting for tumor stage and biology. This implies that “catching it early” is necessary but insufficient to close the gap; disparities in treatment access, quality of care, or specific tumor subtypes likely play additional roles.

Algorithmic Insight

Our logistic regression model demonstrates good predictive power (AUC = 0.78) and algorithmic fairness (comparable AUC across racial groups). The disparity is not an artifact of the model, but a reflection of the reality captured in the data.

Limitations

Missing Variables: The dataset lacks detailed treatment information (chemotherapy/radiation regimens), HER2 status, and socioeconomic variables (income, insurance), which are known drivers of disparities.
“Other” Category: The heterogeneity of the “Other” race category limits specific inferences for Asian, Hispanic, or Indigenous populations.
Temporal Scope: Data from 2006-2010 may not reflect current treatment patterns, including newer targeted therapies and immunotherapies.
External Validation: Models require validation on independent, contemporary cohorts before clinical deployment.

Conclusion

VERDICT

Based on this comprehensive forensic analysis, we find Breast Cancer GUILTY on the following counts:

Count 1 - Causing Mortality: GUILTY — 616 deaths (15.3% fatality rate)

Count 2 - Disparate Impact: GUILTY — Black patients experience 1.7x higher mortality

Count 3 - Algorithmic Bias: NOT GUILTY — Models perform equitably across racial groups

RECOMMENDATION: The prosecution recommends immediate intervention in the form of:

Targeted screening programs for high-risk demographic groups
Policy reviews to address non-biological drivers of the racial survival gap
Research into the structural factors underlying persistent disparities

The prosecution rests.

References

Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B, 34(2), 187-220.
DeSantis, C. E., et al. (2019). Breast cancer statistics, 2019. CA: A Cancer Journal for Clinicians, 69(6), 438-451.
Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457-481.
Newman, L. A., & Kaljee, L. M. (2017). Health disparities and triple-negative breast cancer in African American women. JAMA Surgery, 152(5), 485-493.
Obermeyer, Z., et al. (2019). Dissecting racial bias in an algorithm used to manage the health of populations. Science, 366(6464), 447-453.

Appendix: Session Information

sessionInfo()

## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3   corrplot_0.95   rsample_1.3.1   knitr_1.50     
##  [5] broom_1.0.9     gt_1.0.0        gtsummary_2.2.0 pROC_1.19.0.1  
##  [9] caret_7.0-1     lattice_0.22-7  survminer_0.5.1 ggpubr_0.6.2   
## [13] survival_3.8-3  janitor_2.2.1   lubridate_1.9.4 forcats_1.0.0  
## [17] stringr_1.5.1   dplyr_1.1.4     purrr_1.0.4     readr_2.1.5    
## [21] tidyr_1.3.1     tibble_3.3.0    ggplot2_4.0.1   tidyverse_2.0.0
## 
## loaded via a namespace (and not attached):
##  [1] rlang_1.1.6          magrittr_2.0.3       snakecase_0.11.1    
##  [4] furrr_0.3.1          compiler_4.5.1       vctrs_0.6.5         
##  [7] reshape2_1.4.5       crayon_1.5.3         pkgconfig_2.0.3     
## [10] fastmap_1.2.0        backports_1.5.0      labeling_0.4.3      
## [13] KMsurv_0.1-6         rmarkdown_2.29       markdown_2.0        
## [16] prodlim_2025.04.28   tzdb_0.5.0           bit_4.6.0           
## [19] xfun_0.52            cachem_1.1.0         litedown_0.7        
## [22] jsonlite_2.0.0       recipes_1.3.1        parallel_4.5.1      
## [25] R6_2.6.1             bslib_0.9.0          stringi_1.8.7       
## [28] RColorBrewer_1.1-3   parallelly_1.45.1    car_3.1-3           
## [31] rpart_4.1.24         jquerylib_0.1.4      Rcpp_1.0.14         
## [34] iterators_1.0.14     future.apply_1.20.0  zoo_1.8-14          
## [37] Matrix_1.7-3         splines_4.5.1        nnet_7.3-20         
## [40] timechange_0.3.0     tidyselect_1.2.1     rstudioapi_0.17.1   
## [43] abind_1.4-8          yaml_2.3.10          timeDate_4051.111   
## [46] ggtext_0.1.2         codetools_0.2-20     listenv_0.10.0      
## [49] plyr_1.8.9           withr_3.0.2          S7_0.2.1            
## [52] evaluate_1.0.5       future_1.68.0        xml2_1.4.0          
## [55] survMisc_0.5.6       pillar_1.10.2        carData_3.0-5       
## [58] foreach_1.5.2        stats4_4.5.1         generics_0.1.4      
## [61] vroom_1.6.5          hms_1.1.3            commonmark_1.9.5    
## [64] scales_1.4.0         globals_0.18.0       xtable_1.8-4        
## [67] class_7.3-23         glue_1.8.0           tools_4.5.1         
## [70] data.table_1.17.6    ModelMetrics_1.2.2.2 gower_1.0.2         
## [73] ggsignif_0.6.4       grid_4.5.1           cards_0.7.0         
## [76] ipred_0.9-15         nlme_3.1-168         cardx_0.3.0         
## [79] Formula_1.2-5        cli_3.6.5            km.ci_0.5-6         
## [82] lava_1.8.2           gtable_0.3.6         rstatix_0.7.3       
## [85] sass_0.4.10          digest_0.6.37        farver_2.1.2        
## [88] htmltools_0.5.8.1    lifecycle_1.0.4      hardhat_1.4.2       
## [91] gridtext_0.1.5       bit64_4.6.0-1        MASS_7.3-65

⟵ Previous

To Front