Simulate a Phase II Oncology Trial in R

Introduction

30 Second Summary

Every cancer drug you have ever heard of went through a clinical trial before it reached patients. Behind those trials are biostatisticians who analyze survival data, quantify treatment effects, and produce the reports that regulators use to approve life-saving therapies.

In this project, you will build a complete Phase II oncology trial analysis pipeline in R. You will generate synthetic patient data, perform survival analysis, and produce publication-ready outputs that mirror real pharmaceutical deliverables.

What You'll Build

A self-contained R script that simulates a 200-patient randomized controlled trial of an experimental drug vs. placebo in advanced lung cancer, then performs the full statistical analysis and reporting workflow.

By the end of this project, you'll have:

• A reproducible synthetic dataset with 200 patients including demographics, treatment arms, overall survival times with censoring, and adverse event records.
• A complete primary efficacy analysis featuring a Kaplan-Meier survival curve with risk table, a Cox Proportional Hazards model extracting the treatment's Hazard Ratio, 95% confidence interval, and p-value.
• A formatted Adverse Event summary table grouped by System Organ Class and severity grade, plus a programmatic executive summary following Clinical Study Report conventions.
• Secret Mission: a Forest Plot showing subgroup Hazard Ratios to identify whether specific patient populations drive the treatment effect.

Are there any prerequisites?

You should have basic R familiarity (data frames, pipes, installing packages). No prior biostatistics or clinical trials experience is required.

Before We Start

Before We Start...

Set Up Your R Environment

Before you can simulate a clinical trial, you need the right tools on your machine. R is the programming language used by biostatisticians across the pharmaceutical industry for survival analysis and regulatory reporting.

In this step you will install R, set up RStudio Desktop as your development environment, and load the four packages that power the entire analysis pipeline.

In this step, get ready to:

• Install R and RStudio Desktop.
• Install the tidyverse, survival, survminer, and gtsummary packages.
• Create your analysis script file with the library loading block.

Install R and RStudio Desktop

R is a free, open-source language for statistical computing. This project requires version 4.2 or higher because the gtsummary package depends on it.

• Open your terminal. Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
• Type Terminal (macOS) or Command Prompt (Windows) and press Enter to open it.
• Check whether R is already installed by running this command:

R --version

Look at the first line of output. It shows the R version number in the format R version X.Y.Z.

✔️ I see version 4.2 or higher

Your R installation is ready. Move on to installing RStudio Desktop below.

ⓧ I see an older version

You have R installed, but you need version 4.2 or higher.

• Visit the CRAN download page at https://cran.r-project.org/.
• Click the download link for your operating system (macOS or Windows).
• Download and run the latest installer. It will replace your older version.
• Close and reopen your terminal, then run R --version again to confirm the new version.

macOS users with Homebrew

If you installed R via Homebrew, you can upgrade by running brew upgrade r in your terminal instead.

ⓧ Command not found

R is not installed on your machine yet.

• Visit the CRAN download page at https://cran.r-project.org/.
• Click the download link for your operating system (macOS or Windows).
• Download the latest installer and run it. Accept the default settings when prompted.
• Close and reopen your terminal, then run R --version to confirm R is installed.

What is CRAN?

CRAN stands for the Comprehensive R Archive Network. It is the official repository where R and thousands of R packages are hosted.

Now install RStudio Desktop, the IDE that gives you a code editor, console, and plot viewer all in one window.

• Visit the Posit download page at https://posit.co/download/rstudio-desktop/.
• Download the free RStudio Desktop installer for your operating system.
• Run the installer and accept the default settings.
• Open RStudio. Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
• Type RStudio and press Enter to launch it.

What is RStudio?

RStudio is an IDE (Integrated Development Environment) for R. It gives you a code editor with syntax highlighting, a console to run commands interactively, and built-in panels for viewing plots, data, and help documentation.

RStudio says it cannot find R?

RStudio auto-detects your R installation on startup. If it shows an error, close RStudio, confirm R is on your path by running R --version in your terminal, then reopen RStudio.

help me troubleshoot RStudio not finding my R installation

Install the required R packages

This project uses four packages that are not included in base R. You will install them from CRAN using the R console inside RStudio.

• In RStudio, look at the bottom-left panel labeled Console. This is where you type R commands.
• Install all four packages by pasting and running this command in the Console:

install.packages(c("tidyverse", "survival", "survminer", "gtsummary"))

What are these packages?

tidyverse gives you tools for data manipulation and visualization (dplyr, ggplot2, purrr, and more).

survival provides the core functions for survival analysis, including Surv(), survfit(), and coxph().

survminer creates publication-ready Kaplan-Meier survival curves with risk tables.

gtsummary formats regression tables and summary statistics into clean, presentation-ready output.

R may ask you to select a CRAN mirror. Choose any mirror close to your location and press Enter. Installation takes a few minutes because tidyverse bundles many sub-packages.

• Verify all packages load correctly by running this command in the Console:

library(tidyverse); library(survival); library(survminer); library(gtsummary)

You should see some startup messages from tidyverse listing the attached packages. No error messages should appear.

Seeing an installation error?

If a package fails to install, check that your internet connection is active and try running the install.packages() command again. On macOS, you may be prompted to install from source if a binary is not available. Type yes if prompted.

help me fix an R package installation error

Create your analysis script

All of your analysis code will live in a single R script file. Creating it now gives you a place to add code as you work through each step.

• In RStudio, click File in the top menu bar.
• Click New File, then R Script.
• A blank editor tab opens in the top-left panel. Add the library loading block by pasting this code:

# ============================================================================
# Phase II Oncology Clinical Trial Simulation and Analysis
# Drug: ONC-4587 (Experimental) vs. Placebo
# Indication: Advanced Non-Small Cell Lung Cancer (NSCLC)
# ============================================================================

# --- Load Libraries ---
library(tidyverse)
library(survival)
library(survminer)
library(gtsummary)

What does this code do?

The comment block at the top identifies the study (drug name, indication, study phase). This mirrors how real clinical trial scripts are labeled so collaborators know which protocol the code belongs to.

The four library() calls load the packages you just installed, making their functions available for the rest of the script.

• Save the file by pressing Cmd+S (macOS) or Ctrl+S (Windows).
• In the save dialog, name the file oncology_trial_analysis.R and save it to a location you can find easily (your Desktop works well).
• Verify the file is saved by checking that the tab title in RStudio shows oncology_trial_analysis.R without an asterisk.

Tab still shows an asterisk?

An asterisk next to the filename means unsaved changes. Press Cmd+S (macOS) or Ctrl+S (Windows) again. If the save dialog reappears, confirm the filename and click Save.

help me save my R script file in RStudio

✔️ Awesome, I've got everything!

Great. Double check you have saved oncology_trial_analysis.R and that it has no asterisk in the tab title.

ⓧ I'd like to double check the full code

Your oncology_trial_analysis.R file should contain exactly this:

# ============================================================================
# Phase II Oncology Clinical Trial Simulation and Analysis
# Drug: ONC-4587 (Experimental) vs. Placebo
# Indication: Advanced Non-Small Cell Lung Cancer (NSCLC)
# ============================================================================

# --- Load Libraries ---
library(tidyverse)
library(survival)
library(survminer)
library(gtsummary)

Your environment is locked in. Next up, you will generate a full 200-patient synthetic clinical trial dataset with demographics, survival endpoints, and adverse event records.

Generate the Synthetic Trial Dataset

Your R environment is ready with all packages loaded. Now you need data to analyze. Real clinical trial data is protected by strict regulations, so biostatisticians routinely generate synthetic datasets that mirror the statistical properties of real trials.

In this step, you will build a complete 200-patient Phase II oncology trial dataset from scratch. The dataset includes patient demographics, survival times with realistic censoring, and adverse event records. Every variable mirrors what you would find in a real Clinical Study Report.

In this step, get ready to:

• Create a demographics data frame with 200 patients randomized to ONC-4587 or Placebo.
• Simulate overall survival times with a built-in treatment effect and realistic censoring.
• Generate an adverse events dataset with System Organ Class, preferred terms, and severity grades.

Create the patient demographics data frame

Every clinical trial starts with a patient population. You need to define who is in the study, which arm they were randomized to, and their baseline characteristics.

• In oncology_trial_analysis.R, add the following code below your library loading block:

# --- Set Seed for Reproducibility ---
set.seed(42)

# ============================================================================
# SECTION 1: SYNTHETIC DATA GENERATION
# ============================================================================

n_patients <- 200

# --- Patient Demographics ---
trial_data <- tibble(
  patient_id = sprintf("PT-%04d", 1:n_patients),
  treatment_arm = rep(c("ONC-4587", "Placebo"), each = n_patients / 2),
  age = round(rnorm(n_patients, mean = 62, sd = 9)),
  sex = sample(c("Male", "Female"), n_patients, replace = TRUE, prob = c(0.58, 0.42)),
  ecog_status = sample(0:2, n_patients, replace = TRUE, prob = c(0.3, 0.5, 0.2)),
  tumor_stage = sample(c("IIIB", "IV"), n_patients, replace = TRUE, prob = c(0.35, 0.65))
)

What does this code do?

• set.seed(42) locks the random number generator so every run produces identical data. This is essential for reproducible research.
• n_patients <- 200 sets the trial size. 100 patients per arm is realistic for a Phase II oncology study.
• treatment_arm uses rep(..., each = 100) for a balanced 1:1 randomization between ONC-4587 (experimental) and Placebo.
• ecog_status represents the ECOG Performance Status (0 = fully active, 2 = ambulatory but unable to work). This is a standard baseline variable in oncology trials.

• Select all the new code you just added.
• Press Cmd+Enter (macOS) or Ctrl+Enter (Windows) to run the selection.
• Check the Environment pane in the top-right of RStudio. You should see trial_data listed as a tibble with 200 observations and 6 variables.

Don't see trial_data in the Environment pane?

• Make sure you selected ALL the new code, including set.seed(42) at the top and the closing ) of the tibble at the bottom.
• If you see an error about tibble not being found, your libraries did not load. Run library(tidyverse) in the console first.
• Check that every opening parenthesis has a matching close. A missing ) after the tibble will cause R to wait for more input (you will see a + in the console instead of >).

Help me debug my trial_data creation in R

Simulate overall survival times

The primary endpoint in an oncology trial is typically Overall Survival (OS). You need to simulate how long each patient lives, with the experimental drug producing longer survival times than placebo.

• Add the following code below your trial_data tibble creation:

# --- Simulate Overall Survival (Time-to-Event) ---
trial_data <- trial_data |>
  mutate(
    os_months = if_else(
      treatment_arm == "ONC-4587",
      rweibull(n(), shape = 1.2, scale = 14),
      rweibull(n(), shape = 1.2, scale = 10)
    ),
    censor_time = runif(n(), min = 6, max = 24),
    event_status = if_else(os_months <= censor_time, 1L, 0L),
    os_months = pmin(os_months, censor_time)
  ) |>
  select(-censor_time)

What does this code do?

This block simulates realistic time-to-event data using a Weibull distribution, the standard parametric model for survival data in oncology.

• rweibull(n(), shape = 1.2, scale = 14) gives the ONC-4587 arm longer expected survival (scale = 14 months) compared to Placebo (scale = 10 months). This embeds a true treatment effect.
• censor_time simulates the data cutoff date. Patients still alive at cutoff are "censored." This is how real trials work.
• event_status is 1 if the patient died before censoring (an event occurred) and 0 if they were still alive at data cutoff.
• pmin(os_months, censor_time) caps the observed time at whichever came first: death or censoring. Then select(-censor_time) removes the helper column.

What is censoring?

In a real clinical trial, not every patient has died by the time you analyze the data. Some patients are still alive (their survival time is unknown beyond the observation window). These patients are "right-censored." Your analysis must account for them without throwing away their partial information.

• Select and run this new code block with Cmd+Enter (macOS) or Ctrl+Enter (Windows).
• In the Environment pane, trial_data should now show 200 observations and 8 variables (the original 6 plus os_months and event_status).

Seeing an error about the pipe operator?

The native pipe |> requires R 4.1 or higher. If you see an unexpected symbol error, confirm your R version by typing R.version.string in the console. You need R 4.2 or above.

Help me fix the pipe operator error in my R script

Generate the adverse events dataset

Beyond efficacy, every trial must report safety data. Adverse events (AEs) are categorized by System Organ Class and graded using the CTCAE scale (Grade 1 = mild, Grade 4 = life-threatening).

• Add the following code below the survival simulation block. This defines the vocabulary of possible adverse events:

# --- Generate Adverse Events Dataset ---
soc_options <- c(
  "Gastrointestinal disorders",
  "Skin and subcutaneous tissue disorders",
  "General disorders and administration site conditions",
  "Nervous system disorders",
  "Musculoskeletal and connective tissue disorders",
  "Respiratory, thoracic and mediastinal disorders",
  "Metabolism and nutrition disorders"
)

ae_terms <- list(
  "Gastrointestinal disorders" = c("Nausea", "Vomiting", "Diarrhea", "Constipation"),
  "Skin and subcutaneous tissue disorders" = c("Rash", "Pruritus", "Dry skin"),
  "General disorders and administration site conditions" = c("Fatigue", "Pyrexia", "Injection site reaction"),
  "Nervous system disorders" = c("Headache", "Dizziness", "Peripheral neuropathy"),
  "Musculoskeletal and connective tissue disorders" = c("Arthralgia", "Myalgia", "Back pain"),
  "Respiratory, thoracic and mediastinal disorders" = c("Cough", "Dyspnea", "Pneumonitis"),
  "Metabolism and nutrition disorders" = c("Decreased appetite", "Hyponatremia", "Hyperglycemia")
)

What does this code do?

• soc_options defines seven System Organ Classes (SOCs). These are the MedDRA-standard top-level categories used in real regulatory submissions.
• ae_terms is a named list that maps each SOC to its specific adverse event terms (called "preferred terms" in clinical reporting). When generating AEs, the code picks a SOC first, then randomly selects a preferred term within that SOC.

• Now add the code that generates individual AE records for each patient:

ae_data <- map_dfr(1:n_patients, function(i) {
  n_ae <- sample(1:4, 1, prob = c(0.3, 0.35, 0.25, 0.1))
  socs <- sample(soc_options, n_ae, replace = TRUE)
  tibble(
    patient_id = trial_data$patient_id[i],
    treatment_arm = trial_data$treatment_arm[i],
    system_organ_class = socs,
    preferred_term = map_chr(socs, ~ sample(ae_terms[[.x]], 1)),
    ctcae_grade = sample(1:4, n_ae, replace = TRUE, prob = c(0.40, 0.35, 0.20, 0.05))
  )
})

What does this code do?

• map_dfr(1:n_patients, ...) loops over all 200 patients and row-binds each patient's AE tibble into one combined data frame.
• Each patient gets 1 to 4 AEs, weighted so most patients have 1-2 events (realistic for a Phase II trial).
• ctcae_grade is heavily weighted toward Grade 1-2 (75% combined), with only 5% chance of Grade 4. This mirrors real-world AE severity distributions.

• Finally, add this code to give the experimental arm a slightly higher rate of skin toxicity (a realistic drug class effect) and print a summary:

# Make experimental arm have slightly higher AE rates for skin events
ae_data <- ae_data |>
  bind_rows(
    trial_data |>
      filter(treatment_arm == "ONC-4587") |>
      slice_sample(n = 20) |>
      transmute(
        patient_id,
        treatment_arm,
        system_organ_class = "Skin and subcutaneous tissue disorders",
        preferred_term = sample(c("Rash", "Pruritus", "Dry skin"), 20, replace = TRUE),
        ctcae_grade = sample(1:3, 20, replace = TRUE, prob = c(0.4, 0.4, 0.2))
      )
  )

cat("=== Data Generation Complete ===\n")
cat(sprintf("Patients: %d | Events: %d | Censored: %d\n",
            nrow(trial_data),
            sum(trial_data$event_status == 1),
            sum(trial_data$event_status == 0)))

Why add extra skin events to the experimental arm?

In real oncology trials, experimental drugs often have a signature toxicity profile. Skin reactions are a common class effect for targeted therapies and checkpoint inhibitors. By adding 20 extra skin AE records to the ONC-4587 arm, you create a detectable safety signal that will appear in your AE summary table later.

• Save oncology_trial_analysis.R.
• Select all the code in the adverse events section (from soc_options down to the cat() lines) and run it with Cmd+Enter (macOS) or Ctrl+Enter (Windows).
• In the console, you should see output like:

=== Data Generation Complete ===
Patients: 200 | Events: 143 | Censored: 57

Your exact numbers for Events and Censored may differ slightly depending on your R version, but the total should always be 200.

Seeing an error about trial_data not found?

The ae_data code references trial_data from the previous substep. If you skipped running the demographics and survival chunks earlier, those objects do not exist in memory yet. Select everything from set.seed(42) down to the end and run it all together.

Help me debug my adverse events data generation

✔️ Awesome, I've got everything!

Great. Double check you have saved oncology_trial_analysis.R before moving on.

ⓧ I'd like to double check the full code

Here is the complete oncology_trial_analysis.R file up to this point:

# ============================================================================
# Phase II Oncology Clinical Trial Simulation and Analysis
# Drug: ONC-4587 (Experimental) vs. Placebo
# Indication: Advanced Non-Small Cell Lung Cancer (NSCLC)
# ============================================================================

# --- Load Libraries ---
library(tidyverse)
library(survival)
library(survminer)
library(gtsummary)

# --- Set Seed for Reproducibility ---
set.seed(42)

# ============================================================================
# SECTION 1: SYNTHETIC DATA GENERATION
# ============================================================================

n_patients <- 200

# --- Patient Demographics ---
trial_data <- tibble(
  patient_id = sprintf("PT-%04d", 1:n_patients),
  treatment_arm = rep(c("ONC-4587", "Placebo"), each = n_patients / 2),
  age = round(rnorm(n_patients, mean = 62, sd = 9)),
  sex = sample(c("Male", "Female"), n_patients, replace = TRUE, prob = c(0.58, 0.42)),
  ecog_status = sample(0:2, n_patients, replace = TRUE, prob = c(0.3, 0.5, 0.2)),
  tumor_stage = sample(c("IIIB", "IV"), n_patients, replace = TRUE, prob = c(0.35, 0.65))
)

# --- Simulate Overall Survival (Time-to-Event) ---
trial_data <- trial_data |>
  mutate(
    os_months = if_else(
      treatment_arm == "ONC-4587",
      rweibull(n(), shape = 1.2, scale = 14),
      rweibull(n(), shape = 1.2, scale = 10)
    ),
    censor_time = runif(n(), min = 6, max = 24),
    event_status = if_else(os_months <= censor_time, 1L, 0L),
    os_months = pmin(os_months, censor_time)
  ) |>
  select(-censor_time)

# --- Generate Adverse Events Dataset ---
soc_options <- c(
  "Gastrointestinal disorders",
  "Skin and subcutaneous tissue disorders",
  "General disorders and administration site conditions",
  "Nervous system disorders",
  "Musculoskeletal and connective tissue disorders",
  "Respiratory, thoracic and mediastinal disorders",
  "Metabolism and nutrition disorders"
)

ae_terms <- list(
  "Gastrointestinal disorders" = c("Nausea", "Vomiting", "Diarrhea", "Constipation"),
  "Skin and subcutaneous tissue disorders" = c("Rash", "Pruritus", "Dry skin"),
  "General disorders and administration site conditions" = c("Fatigue", "Pyrexia", "Injection site reaction"),
  "Nervous system disorders" = c("Headache", "Dizziness", "Peripheral neuropathy"),
  "Musculoskeletal and connective tissue disorders" = c("Arthralgia", "Myalgia", "Back pain"),
  "Respiratory, thoracic and mediastinal disorders" = c("Cough", "Dyspnea", "Pneumonitis"),
  "Metabolism and nutrition disorders" = c("Decreased appetite", "Hyponatremia", "Hyperglycemia")
)

ae_data <- map_dfr(1:n_patients, function(i) {
  n_ae <- sample(1:4, 1, prob = c(0.3, 0.35, 0.25, 0.1))
  socs <- sample(soc_options, n_ae, replace = TRUE)
  tibble(
    patient_id = trial_data$patient_id[i],
    treatment_arm = trial_data$treatment_arm[i],
    system_organ_class = socs,
    preferred_term = map_chr(socs, ~ sample(ae_terms[[.x]], 1)),
    ctcae_grade = sample(1:4, n_ae, replace = TRUE, prob = c(0.40, 0.35, 0.20, 0.05))
  )
})

# Make experimental arm have slightly higher AE rates for skin events
ae_data <- ae_data |>
  bind_rows(
    trial_data |>
      filter(treatment_arm == "ONC-4587") |>
      slice_sample(n = 20) |>
      transmute(
        patient_id,
        treatment_arm,
        system_organ_class = "Skin and subcutaneous tissue disorders",
        preferred_term = sample(c("Rash", "Pruritus", "Dry skin"), 20, replace = TRUE),
        ctcae_grade = sample(1:3, 20, replace = TRUE, prob = c(0.4, 0.4, 0.2))
      )
  )

cat("=== Data Generation Complete ===\n")
cat(sprintf("Patients: %d | Events: %d | Censored: %d\n",
            nrow(trial_data),
            sum(trial_data$event_status == 1),
            sum(trial_data$event_status == 0)))

You now have a complete synthetic trial dataset sitting in memory. Next up, you will feed this data into a Kaplan-Meier survival analysis and produce your first publication-ready survival curve.

Perform the Kaplan-Meier Survival Analysis

In the previous step, you built a complete synthetic trial dataset with 200 patients, survival times, and adverse event records. Now it is time to answer the central question of any oncology trial: does the experimental drug actually help patients live longer?

The Kaplan-Meier method estimates the probability of survival over time for each treatment arm. Combined with a log-rank test and a publication-quality plot, this is the primary efficacy analysis you would find in a real Clinical Study Report.

In this step, get ready to:

• Fit Kaplan-Meier survival curves and extract median survival times.
• Perform a log-rank test to formally compare the two treatment arms.
• Generate a publication-ready survival plot with a number-at-risk table and save it as a PNG file.

Fit the Kaplan-Meier curves stratified by treatment arm

A Surv object pairs each patient's survival time with their event indicator (did the event occur, or were they censored?). Passing this into survfit() stratified by treatment_arm gives you a separate survival curve for each group.

• In oncology_trial_analysis.R, scroll to the bottom of your file (after the data generation complete message).
• Add the section header and KM fitting code by pasting the following:

# ============================================================================
# SECTION 2: KAPLAN-MEIER SURVIVAL ANALYSIS
# ============================================================================

# --- Fit KM Curves ---
km_fit <- survfit(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
cat("\n=== Kaplan-Meier Median Survival ===\n")
print(km_fit)

What does this code do?

Surv(os_months, event_status) creates a survival object where os_months is the observed time and event_status indicates whether the patient died (1) or was censored (0).

The formula ~ treatment_arm tells survfit() to compute separate curves for ONC-4587 and Placebo. Printing the result shows the median overall survival and 95% confidence interval for each arm.

• Save oncology_trial_analysis.R.
• Select the new code you just added and click **Run** (or press **Ctrl+Enter** on Windows, **Cmd+Enter** on macOS) to execute it.
• In the console, you should see a table showing the number of patients, events, and median survival for each treatment arm.

Seeing an error about Surv not found?

Make sure library(survival) is loaded at the top of your script and has been executed in your current R session. If you restarted RStudio, re-run the full script from the top to reload all libraries and regenerate the data.

help me debug a Surv not found error in my R survival analysis script

Perform the log-rank test

The Kaplan-Meier curves show you the shape of survival, but you need a formal statistical test to determine whether the difference between arms is significant. The log-rank test compares the entire survival distributions under the null hypothesis that both arms are identical.

• Add the log-rank test code below your KM fitting block:

# --- Log-Rank Test ---
logrank_test <- survdiff(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
logrank_p <- 1 - pchisq(logrank_test$chisq, df = 1)
cat(sprintf("\nLog-Rank Test p-value: %.4f\n", logrank_p))

What does this code do?

survdiff() performs the log-rank test (also called the Mantel-Haenszel test) by comparing observed vs. expected events in each group. It returns a chi-squared statistic stored in logrank_test$chisq.

The p-value is calculated from the chi-squared distribution with 1 degree of freedom. A p-value below 0.05 indicates a statistically significant difference in survival between the ONC-4587 and Placebo arms.

• Save oncology_trial_analysis.R and run the new lines.
• In the console, you should see a line reading Log-Rank Test p-value: followed by a four-decimal number (likely below 0.05 given the built-in treatment effect).

Getting a p-value of exactly 0 or NaN?

Confirm that your event_status column contains both 0s and 1s. If all patients are censored (all 0s), the test cannot compute a valid statistic. Re-run the data generation section from the top of the script to regenerate trial_data.

help me troubleshoot a log-rank test returning NaN in my R script

Generate the publication-ready KM plot with a risk table

A survival curve without a number-at-risk table is incomplete for publication. The ggsurvplot() function from survminer combines both into a single figure with confidence intervals, p-value annotation, and a clean theme.

• Add the plotting code below your log-rank test block:

# --- Publication-Ready KM Plot with Risk Table ---
km_plot <- ggsurvplot(
  km_fit,
  data = trial_data,
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  risk.table.col = "strata",
  palette = c("#E7B800", "#2E9FDF"),
  xlab = "Time (Months)",
  ylab = "Overall Survival Probability",
  title = "Kaplan-Meier Estimate of Overall Survival\nONC-4587 vs. Placebo in Advanced NSCLC",
  legend.title = "Treatment Arm",
  legend.labs = c("ONC-4587", "Placebo"),
  ggtheme = theme_minimal(base_size = 12),
  risk.table.height = 0.25,
  surv.median.line = "hv"
)

print(km_plot)
ggsave("km_survival_plot.png", plot = print(km_plot), width = 10, height = 7, dpi = 300)
cat("\nKM plot saved as 'km_survival_plot.png'\n")

What does this code do?

• pval = TRUE overlays the log-rank p-value directly on the plot.
• conf.int = TRUE draws shaded 95% confidence bands around each curve.
• risk.table = TRUE adds a number-at-risk table below the figure, which regulatory agencies expect to see.
• surv.median.line = "hv" draws horizontal and vertical reference lines at the median survival point for each arm.
• ggsave() exports the composite figure as a 300 DPI PNG file suitable for publication or a CSR appendix.

• Save oncology_trial_analysis.R and run the new lines.
• In the RStudio **Plots** pane, you should see a Kaplan-Meier survival curve with two lines (gold for ONC-4587, blue for Placebo), confidence bands, a p-value annotation, and a risk table below.
• Check your working directory for the file km_survival_plot.png by running:

file.exists("km_survival_plot.png")

You should see TRUE printed in the console.

Plot not appearing or getting a data argument error?

The data = trial_data argument in ggsurvplot() is required. Without it, survminer cannot extract the original data frame from the fit object and you will see an error. Make sure trial_data exists in your environment (check the **Environment** pane in RStudio).

help me fix a ggsurvplot data argument error in my KM analysis

✔️ Awesome, I've got everything!

Great. Double check you have saved oncology_trial_analysis.R and that km_survival_plot.png exists in your working directory.

ⓧ I'd like to double check the full code

Here is the complete oncology_trial_analysis.R file up to this point:

# ============================================================================
# Phase II Oncology Clinical Trial Simulation and Analysis
# Drug: ONC-4587 (Experimental) vs. Placebo
# Indication: Advanced Non-Small Cell Lung Cancer (NSCLC)
# ============================================================================

# --- Load Libraries ---
library(tidyverse)
library(survival)
library(survminer)
library(gtsummary)

# --- Set Seed for Reproducibility ---
set.seed(42)

# ============================================================================
# SECTION 1: SYNTHETIC DATA GENERATION
# ============================================================================

n_patients <- 200

trial_data <- tibble(
  patient_id = sprintf("PT-%04d", 1:n_patients),
  treatment_arm = rep(c("ONC-4587", "Placebo"), each = n_patients / 2),
  age = round(rnorm(n_patients, mean = 62, sd = 9)),
  sex = sample(c("Male", "Female"), n_patients, replace = TRUE, prob = c(0.58, 0.42)),
  ecog_status = sample(0:2, n_patients, replace = TRUE, prob = c(0.3, 0.5, 0.2)),
  tumor_stage = sample(c("IIIB", "IV"), n_patients, replace = TRUE, prob = c(0.35, 0.65))
)

trial_data <- trial_data |>
  mutate(
    os_months = if_else(
      treatment_arm == "ONC-4587",
      rweibull(n(), shape = 1.2, scale = 14),
      rweibull(n(), shape = 1.2, scale = 10)
    ),
    censor_time = runif(n(), min = 6, max = 24),
    event_status = if_else(os_months <= censor_time, 1L, 0L),
    os_months = pmin(os_months, censor_time)
  ) |>
  select(-censor_time)

soc_options <- c(
  "Gastrointestinal disorders",
  "Skin and subcutaneous tissue disorders",
  "General disorders and administration site conditions",
  "Nervous system disorders",
  "Musculoskeletal and connective tissue disorders",
  "Respiratory, thoracic and mediastinal disorders",
  "Metabolism and nutrition disorders"
)

ae_terms <- list(
  "Gastrointestinal disorders" = c("Nausea", "Vomiting", "Diarrhea", "Constipation"),
  "Skin and subcutaneous tissue disorders" = c("Rash", "Pruritus", "Dry skin"),
  "General disorders and administration site conditions" = c("Fatigue", "Pyrexia", "Injection site reaction"),
  "Nervous system disorders" = c("Headache", "Dizziness", "Peripheral neuropathy"),
  "Musculoskeletal and connective tissue disorders" = c("Arthralgia", "Myalgia", "Back pain"),
  "Respiratory, thoracic and mediastinal disorders" = c("Cough", "Dyspnea", "Pneumonitis"),
  "Metabolism and nutrition disorders" = c("Decreased appetite", "Hyponatremia", "Hyperglycemia")
)

ae_data <- map_dfr(1:n_patients, function(i) {
  n_ae <- sample(1:4, 1, prob = c(0.3, 0.35, 0.25, 0.1))
  socs <- sample(soc_options, n_ae, replace = TRUE)
  tibble(
    patient_id = trial_data$patient_id[i],
    treatment_arm = trial_data$treatment_arm[i],
    system_organ_class = socs,
    preferred_term = map_chr(socs, ~ sample(ae_terms[[.x]], 1)),
    ctcae_grade = sample(1:4, n_ae, replace = TRUE, prob = c(0.40, 0.35, 0.20, 0.05))
  )
})

ae_data <- ae_data |>
  bind_rows(
    trial_data |>
      filter(treatment_arm == "ONC-4587") |>
      slice_sample(n = 20) |>
      transmute(
        patient_id,
        treatment_arm,
        system_organ_class = "Skin and subcutaneous tissue disorders",
        preferred_term = sample(c("Rash", "Pruritus", "Dry skin"), 20, replace = TRUE),
        ctcae_grade = sample(1:3, 20, replace = TRUE, prob = c(0.4, 0.4, 0.2))
      )
  )

cat("=== Data Generation Complete ===\n")
cat(sprintf("Patients: %d | Events: %d | Censored: %d\n",
            nrow(trial_data),
            sum(trial_data$event_status == 1),
            sum(trial_data$event_status == 0)))

# ============================================================================
# SECTION 2: KAPLAN-MEIER SURVIVAL ANALYSIS
# ============================================================================

km_fit <- survfit(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
cat("\n=== Kaplan-Meier Median Survival ===\n")
print(km_fit)

logrank_test <- survdiff(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
logrank_p <- 1 - pchisq(logrank_test$chisq, df = 1)
cat(sprintf("\nLog-Rank Test p-value: %.4f\n", logrank_p))

km_plot <- ggsurvplot(
  km_fit,
  data = trial_data,
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  risk.table.col = "strata",
  palette = c("#E7B800", "#2E9FDF"),
  xlab = "Time (Months)",
  ylab = "Overall Survival Probability",
  title = "Kaplan-Meier Estimate of Overall Survival\nONC-4587 vs. Placebo in Advanced NSCLC",
  legend.title = "Treatment Arm",
  legend.labs = c("ONC-4587", "Placebo"),
  ggtheme = theme_minimal(base_size = 12),
  risk.table.height = 0.25,
  surv.median.line = "hv"
)

print(km_plot)
ggsave("km_survival_plot.png", plot = print(km_plot), width = 10, height = 7, dpi = 300)
cat("\nKM plot saved as 'km_survival_plot.png'\n")

You now have a complete Kaplan-Meier analysis with a formal statistical test and a publication-quality figure. Next up, you will fit a Cox Proportional Hazards model to quantify the treatment effect as a Hazard Ratio while adjusting for patient characteristics.

Fit the Cox Proportional Hazards Model

Your Kaplan-Meier analysis showed whether the survival curves differ between treatment arms. Now you need to quantify how much the treatment reduces the risk of death, while controlling for patient characteristics that could confound the result.

The Cox Proportional Hazards model is the standard regression method for time-to-event data in clinical trials. It produces a Hazard Ratio that tells you how the risk of death in one group compares to another, adjusted for covariates like age and performance status.

In this step, get ready to:

• Fit a multivariable Cox PH model and extract the Hazard Ratio with 95% CI.
• Validate the proportional hazards assumption using Schoenfeld residuals.
• Create a publication-ready regression summary table.

Fit the Cox PH model and extract the treatment effect

The KM curve gives you a visual comparison, but regulators need a single number with a confidence interval. The Cox model provides exactly that by estimating how each covariate shifts the hazard (instantaneous risk of death).

• In oncology_trial_analysis.R, scroll to the bottom of your file (after the KM plot section).
• Add the Cox model section header and model fitting code by pasting this below your existing code:

# ============================================================================
# SECTION 3: COX PROPORTIONAL HAZARDS MODEL
# ============================================================================

# --- Fit Multivariable Cox Model ---
cox_model <- coxph(
  Surv(os_months, event_status) ~ treatment_arm + age + sex + ecog_status,
  data = trial_data
)

cat("\n=== Cox Proportional Hazards Model Summary ===\n")
cox_summary <- summary(cox_model)
print(cox_summary)

What does this code do?

coxph() fits a Cox Proportional Hazards regression model. The formula Surv(os_months, event_status) ~ treatment_arm + age + sex + ecog_status tells R to model the hazard of death as a function of treatment arm, age, sex, and ECOG status simultaneously.

By including covariates beyond just treatment arm, you get an adjusted estimate of the treatment effect. This is the standard approach in clinical trial primary analyses.

• Source your script by pressing Cmd+Shift+S (macOS) or Ctrl+Shift+S (Windows) in RStudio.
• You should see a full model summary in the console, including coefficient estimates, standard errors, and p-values for each covariate.

Seeing an error about the formula?

Double-check that your trial_data object exists in the environment. If you restarted R without sourcing the full script, the data generation code needs to run first.

Also verify the column names match exactly: os_months, event_status, treatment_arm, age, sex, ecog_status.

Help me debug my coxph() formula error

Now you need to extract the specific treatment effect statistic. Regulators and reviewers want a single Hazard Ratio with its confidence interval and p-value.

• Add the treatment effect extraction code directly below the print(cox_summary) line:

# --- Extract Treatment Effect ---
hr <- exp(coef(cox_model)["treatment_armPlacebo"])
hr_ci <- exp(confint(cox_model)["treatment_armPlacebo", ])
hr_p <- cox_summary$coefficients["treatment_armPlacebo", "Pr(>|z|)"]

cat(sprintf("\n--- Primary Endpoint Result ---\n"))
cat(sprintf("Hazard Ratio (Placebo vs. ONC-4587): %.3f\n", hr))
cat(sprintf("95%% CI: [%.3f, %.3f]\n", hr_ci[1], hr_ci[2]))
cat(sprintf("p-value: %.4f\n", hr_p))

Why treatment_armPlacebo?

R creates coefficient names by combining the variable name with the non-reference level. Since "ONC-4587" sorts alphabetically before "Placebo", R treats ONC-4587 as the reference. The coefficient treatment_armPlacebo represents the effect of being in the Placebo arm relative to ONC-4587.

exp(coef(cox_model)) converts the log-hazard coefficient to a Hazard Ratio. An HR greater than 1.0 means the Placebo group has a higher hazard (worse survival) compared to the experimental drug.

• Save oncology_trial_analysis.R and source the script again.
• You should see a Primary Endpoint Result block in the console showing the Hazard Ratio, 95% CI, and p-value.

Getting a subscript out of bounds error?

This usually means the coefficient name doesn't match. Run names(coef(cox_model)) in the console to see the exact coefficient names R created. Make sure you're using treatment_armPlacebo (capital P, no spaces).

Help me fix my coefficient extraction error

Test the proportional hazards assumption

The Cox model assumes that the hazard ratio between groups stays constant over time. If this assumption is violated, your HR estimate may be misleading.

• Add the proportional hazards assumption test below your extraction code:

# --- Test Proportional Hazards Assumption ---
ph_test <- cox.zph(cox_model)
cat("\n=== Proportional Hazards Assumption Test ===\n")
print(ph_test)

How do I read the cox.zph() output?

cox.zph() performs the Schoenfeld residuals test for each covariate in your model. It produces a table with columns for the chi-squared statistic, degrees of freedom, and p-value.

The interpretation is straightforward: if the p-value for a covariate is greater than 0.05, the proportional hazards assumption is met for that variable. If any p-value is below 0.05, the hazard ratio for that variable changes over time and you may need a time-varying coefficient.

• Save oncology_trial_analysis.R and source the script.
• You should see a table with one row per covariate (treatment_arm, age, sex, ecog_status) plus a GLOBAL row, each showing a p-value.

Seeing a p-value below 0.05?

Because our data is synthetically generated, you might occasionally see a borderline p-value depending on random seed behavior across R versions. For this project, any result is valid learning material. In real trials, a violation would trigger discussions with the statistical team about stratified models or time-varying coefficients.

Help me interpret my cox.zph() results

Create a publication-ready regression table

Raw R output is fine for analysis, but clinical study reports require formatted tables. The gtsummary package can convert your Cox model into a presentation-ready table with one function call.

• Add the formatted regression table code below your PH assumption test:

# --- Formatted Regression Table ---
tbl_cox <- cox_model |>
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      treatment_arm ~ "Treatment Arm",
      age ~ "Age (years)",
      sex ~ "Sex",
      ecog_status ~ "ECOG Performance Status"
    )
  ) |>
  bold_labels() |>
  bold_p(t = 0.05)

print(tbl_cox)

What does exponentiate = TRUE do?

Cox models estimate coefficients on the log-hazard scale by default. Setting exponentiate = TRUE in tbl_regression() automatically exponentiates the coefficients to display Hazard Ratios instead of log-hazard ratios.

bold_labels() highlights the variable group names, and bold_p(t = 0.05) bolds any p-value that reaches statistical significance. These formatting touches are expected in regulatory submissions.

• Save oncology_trial_analysis.R and source the script.
• You should see a neatly formatted table in the console (or RStudio Viewer pane) showing each covariate's Hazard Ratio, 95% CI, and p-value with clean labels.

Table not rendering cleanly?

If the table appears as raw text instead of a formatted display, try running print(tbl_cox) directly in the R console. The Viewer pane rendering depends on your RStudio version. The console output still shows the correct values regardless of formatting.

Help me get my gtsummary table to display correctly

✔️ Awesome, I've got everything!

Great. Double check you've saved your file before moving on.

ⓧ I'd like to double check the full code

Here is the complete oncology_trial_analysis.R file up to this point:

# ============================================================================
# Phase II Oncology Clinical Trial Simulation and Analysis
# Drug: ONC-4587 (Experimental) vs. Placebo
# Indication: Advanced Non-Small Cell Lung Cancer (NSCLC)
# ============================================================================

# --- Load Libraries ---
library(tidyverse)
library(survival)
library(survminer)
library(gtsummary)

# --- Set Seed for Reproducibility ---
set.seed(42)

# ============================================================================
# SECTION 1: SYNTHETIC DATA GENERATION
# ============================================================================

n_patients <- 200

trial_data <- tibble(
  patient_id = sprintf("PT-%04d", 1:n_patients),
  treatment_arm = rep(c("ONC-4587", "Placebo"), each = n_patients / 2),
  age = round(rnorm(n_patients, mean = 62, sd = 9)),
  sex = sample(c("Male", "Female"), n_patients, replace = TRUE, prob = c(0.58, 0.42)),
  ecog_status = sample(0:2, n_patients, replace = TRUE, prob = c(0.3, 0.5, 0.2)),
  tumor_stage = sample(c("IIIB", "IV"), n_patients, replace = TRUE, prob = c(0.35, 0.65))
)

trial_data <- trial_data |>
  mutate(
    os_months = if_else(
      treatment_arm == "ONC-4587",
      rweibull(n(), shape = 1.2, scale = 14),
      rweibull(n(), shape = 1.2, scale = 10)
    ),
    censor_time = runif(n(), min = 6, max = 24),
    event_status = if_else(os_months <= censor_time, 1L, 0L),
    os_months = pmin(os_months, censor_time)
  ) |>
  select(-censor_time)

soc_options <- c(
  "Gastrointestinal disorders",
  "Skin and subcutaneous tissue disorders",
  "General disorders and administration site conditions",
  "Nervous system disorders",
  "Musculoskeletal and connective tissue disorders",
  "Respiratory, thoracic and mediastinal disorders",
  "Metabolism and nutrition disorders"
)

ae_terms <- list(
  "Gastrointestinal disorders" = c("Nausea", "Vomiting", "Diarrhea", "Constipation"),
  "Skin and subcutaneous tissue disorders" = c("Rash", "Pruritus", "Dry skin"),
  "General disorders and administration site conditions" = c("Fatigue", "Pyrexia", "Injection site reaction"),
  "Nervous system disorders" = c("Headache", "Dizziness", "Peripheral neuropathy"),
  "Musculoskeletal and connective tissue disorders" = c("Arthralgia", "Myalgia", "Back pain"),
  "Respiratory, thoracic and mediastinal disorders" = c("Cough", "Dyspnea", "Pneumonitis"),
  "Metabolism and nutrition disorders" = c("Decreased appetite", "Hyponatremia", "Hyperglycemia")
)

ae_data <- map_dfr(1:n_patients, function(i) {
  n_ae <- sample(1:4, 1, prob = c(0.3, 0.35, 0.25, 0.1))
  socs <- sample(soc_options, n_ae, replace = TRUE)
  tibble(
    patient_id = trial_data$patient_id[i],
    treatment_arm = trial_data$treatment_arm[i],
    system_organ_class = socs,
    preferred_term = map_chr(socs, ~ sample(ae_terms[[.x]], 1)),
    ctcae_grade = sample(1:4, n_ae, replace = TRUE, prob = c(0.40, 0.35, 0.20, 0.05))
  )
})

ae_data <- ae_data |>
  bind_rows(
    trial_data |>
      filter(treatment_arm == "ONC-4587") |>
      slice_sample(n = 20) |>
      transmute(
        patient_id,
        treatment_arm,
        system_organ_class = "Skin and subcutaneous tissue disorders",
        preferred_term = sample(c("Rash", "Pruritus", "Dry skin"), 20, replace = TRUE),
        ctcae_grade = sample(1:3, 20, replace = TRUE, prob = c(0.4, 0.4, 0.2))
      )
  )

cat("=== Data Generation Complete ===\n")
cat(sprintf("Patients: %d | Events: %d | Censored: %d\n",
            nrow(trial_data),
            sum(trial_data$event_status == 1),
            sum(trial_data$event_status == 0)))

# ============================================================================
# SECTION 2: KAPLAN-MEIER SURVIVAL ANALYSIS
# ============================================================================

km_fit <- survfit(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
cat("\n=== Kaplan-Meier Median Survival ===\n")
print(km_fit)

logrank_test <- survdiff(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
logrank_p <- 1 - pchisq(logrank_test$chisq, df = 1)
cat(sprintf("\nLog-Rank Test p-value: %.4f\n", logrank_p))

km_plot <- ggsurvplot(
  km_fit,
  data = trial_data,
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  risk.table.col = "strata",
  palette = c("#E7B800", "#2E9FDF"),
  xlab = "Time (Months)",
  ylab = "Overall Survival Probability",
  title = "Kaplan-Meier Estimate of Overall Survival\nONC-4587 vs. Placebo in Advanced NSCLC",
  legend.title = "Treatment Arm",
  legend.labs = c("ONC-4587", "Placebo"),
  ggtheme = theme_minimal(base_size = 12),
  risk.table.height = 0.25,
  surv.median.line = "hv"
)

print(km_plot)
ggsave("km_survival_plot.png", plot = print(km_plot), width = 10, height = 7, dpi = 300)
cat("\nKM plot saved as 'km_survival_plot.png'\n")

# ============================================================================
# SECTION 3: COX PROPORTIONAL HAZARDS MODEL
# ============================================================================

cox_model <- coxph(
  Surv(os_months, event_status) ~ treatment_arm + age + sex + ecog_status,
  data = trial_data
)

cat("\n=== Cox Proportional Hazards Model Summary ===\n")
cox_summary <- summary(cox_model)
print(cox_summary)

hr <- exp(coef(cox_model)["treatment_armPlacebo"])
hr_ci <- exp(confint(cox_model)["treatment_armPlacebo", ])
hr_p <- cox_summary$coefficients["treatment_armPlacebo", "Pr(>|z|)"]

cat(sprintf("\n--- Primary Endpoint Result ---\n"))
cat(sprintf("Hazard Ratio (Placebo vs. ONC-4587): %.3f\n", hr))
cat(sprintf("95%% CI: [%.3f, %.3f]\n", hr_ci[1], hr_ci[2]))
cat(sprintf("p-value: %.4f\n", hr_p))

ph_test <- cox.zph(cox_model)
cat("\n=== Proportional Hazards Assumption Test ===\n")
print(ph_test)

tbl_cox <- cox_model |>
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      treatment_arm ~ "Treatment Arm",
      age ~ "Age (years)",
      sex ~ "Sex",
      ecog_status ~ "ECOG Performance Status"
    )
  ) |>
  bold_labels() |>
  bold_p(t = 0.05)

print(tbl_cox)

You now have a complete Cox regression analysis with a validated assumption test and a table ready for publication. Next up, you'll summarize the safety profile by building an Adverse Event table grouped by severity.

Create the Adverse Event Summary Table

Your Cox model quantified the survival benefit of ONC-4587. Efficacy is only half the story for regulatory reviewers, though.

Every Clinical Study Report must include a safety analysis showing which adverse events occurred, how severe they were, and whether the experimental arm caused more harm than placebo. In this step, you will summarize adverse events by System Organ Class and CTCAE severity grade, then format the results into a publication-ready table.

In this step, get ready to:

• Summarize adverse events by System Organ Class and CTCAE Grade per treatment arm.
• Format the AE summary into a presentation-ready table using gtsummary.

Summarize AEs by System Organ Class and grade

A single patient can have the same adverse event recorded multiple times at different grades. Regulatory tables report each patient only once per System Organ Class, counting their worst (maximum) grade.

• In oncology_trial_analysis.R, scroll to the bottom of your script (below the print(tbl_cox) line).
• Add the Section 4 header and the patient-level AE rollup by copying this code:

# ============================================================================
# SECTION 4: ADVERSE EVENT SUMMARY TABLE
# ============================================================================

cat("\n=== Adverse Event Summary ===\n")

# Summarize AEs by SOC and treatment arm (one AE per patient per SOC, max grade)
ae_summary <- ae_data |>
  group_by(patient_id, treatment_arm, system_organ_class) |>
  summarise(max_grade = max(ctcae_grade), .groups = "drop")

What does this code do?

The ae_data frame has multiple rows per patient (one per AE event). This pipeline collapses them down to one row per patient per System Organ Class.

• group_by(patient_id, treatment_arm, system_organ_class) groups so that each patient-SOC combination becomes one unit.
• summarise(max_grade = max(ctcae_grade)) keeps only the worst grade the patient experienced in that organ class. This is the standard "maximum grade per patient per SOC" convention used in regulatory AE tables.

Now you need to turn these max grades into a count table that shows how many patients hit each severity level in each organ class.

• Add the following code directly below the ae_summary block you just pasted:

ae_table <- ae_summary |>
  mutate(
    grade_category = case_when(
      max_grade == 1 ~ "Grade 1",
      max_grade == 2 ~ "Grade 2",
      max_grade >= 3 ~ "Grade 3-4"
    )
  ) |>
  count(treatment_arm, system_organ_class, grade_category) |>
  pivot_wider(names_from = grade_category, values_from = n, values_fill = 0) |>
  arrange(treatment_arm, desc(`Grade 3-4`))

What does this code do?

• case_when() bins the numeric grades into clinical categories. Grades 3 and 4 are combined because both represent "severe" events that regulators scrutinize.
• count() tallies patients per arm, SOC, and grade category.
• pivot_wider() reshapes the long data into a wide table with one column per grade category, matching the layout regulators expect.
• arrange(treatment_arm, desc(`Grade 3-4`)) sorts so the most concerning organ classes (highest Grade 3-4 counts) appear first within each arm.

• Save oncology_trial_analysis.R.
• Select all the new Section 4 code you just added and press Cmd+Enter (macOS) or Ctrl+Enter (Windows) to run it.
• In the R console, type ae_table and press Enter. You should see a wide tibble with columns for treatment_arm, system_organ_class, Grade 1, Grade 2, and Grade 3-4.

Seeing an error about ae_data not found?

The ae_summary pipeline depends on ae_data from Section 1. Make sure you have run all prior sections first so that ae_data exists in your environment. You can check by typing ls() in the console to list all objects.

If objects are missing, select all code from the top of the file and re-run it with Cmd+Enter (macOS) or Ctrl+Enter (Windows).

Help me debug my AE summary code in R

Format the AE table with gtsummary

The wide ae_table is useful for internal review, but regulatory submissions need polished formatting. The gtsummary package produces presentation-ready tables that match CSR standards with minimal code.

• Add the following code below your ae_table block to create the formatted gtsummary table:

# gtsummary AE table
ae_for_tbl <- ae_summary |>
  mutate(
    any_ae = TRUE,
    grade_3_4 = max_grade >= 3
  )

ae_gtsummary <- ae_for_tbl |>
  select(treatment_arm, system_organ_class, grade_3_4) |>
  tbl_summary(
    by = treatment_arm,
    label = list(
      system_organ_class ~ "System Organ Class",
      grade_3_4 ~ "Grade >= 3 Event"
    )
  ) |>
  bold_labels() |>
  modify_header(label = "**Adverse Event Category**")

print(ae_gtsummary)

What does this code do?

• ae_for_tbl adds two boolean columns: any_ae (always TRUE, used later in the executive summary) and grade_3_4 (TRUE when the max grade is 3 or higher).
• tbl_summary(by = treatment_arm) creates a summary table stratified by treatment arm, showing counts and percentages for each SOC and grade category.
• bold_labels() and modify_header() apply formatting that matches regulatory table conventions.

• Save oncology_trial_analysis.R.
• Select the new gtsummary code and run it with Cmd+Enter (macOS) or Ctrl+Enter (Windows).
• You should see a formatted table in the console showing adverse event counts and percentages split by ONC-4587 and Placebo columns. Skin and subcutaneous tissue disorders should appear more frequently in the ONC-4587 arm.

Table not displaying?

If you see an error about tbl_summary not being found, make sure library(gtsummary) is loaded at the top of your script. Re-run that line if your R session restarted.

If the table displays in the RStudio Viewer pane instead of the console, that is normal behavior. The content is the same.

Help me fix my gtsummary AE table

Why collapse Grades 3 and 4 together?

CTCAE Grade 3 means "severe" and Grade 4 means "life-threatening." Regulatory reviewers focus on the combined Grade 3-4 rate because both grades trigger dose modifications or treatment discontinuation. Reporting them together gives a single "serious toxicity" signal that drives clinical decision-making.

✔️ Awesome, I've got everything!

Great. Double check you have saved your file before moving on.

ⓧ I'd like to double check the full code

Here is the complete oncology_trial_analysis.R file including all sections up to Section 4:

# ============================================================================
# Phase II Oncology Clinical Trial Simulation and Analysis
# Drug: ONC-4587 (Experimental) vs. Placebo
# Indication: Advanced Non-Small Cell Lung Cancer (NSCLC)
# ============================================================================

# --- Load Libraries ---
library(tidyverse)
library(survival)
library(survminer)
library(gtsummary)

# --- Set Seed for Reproducibility ---
set.seed(42)

# ============================================================================
# SECTION 1: SYNTHETIC DATA GENERATION
# ============================================================================

n_patients <- 200

trial_data <- tibble(
  patient_id = sprintf("PT-%04d", 1:n_patients),
  treatment_arm = rep(c("ONC-4587", "Placebo"), each = n_patients / 2),
  age = round(rnorm(n_patients, mean = 62, sd = 9)),
  sex = sample(c("Male", "Female"), n_patients, replace = TRUE, prob = c(0.58, 0.42)),
  ecog_status = sample(0:2, n_patients, replace = TRUE, prob = c(0.3, 0.5, 0.2)),
  tumor_stage = sample(c("IIIB", "IV"), n_patients, replace = TRUE, prob = c(0.35, 0.65))
)

trial_data <- trial_data |>
  mutate(
    os_months = if_else(
      treatment_arm == "ONC-4587",
      rweibull(n(), shape = 1.2, scale = 14),
      rweibull(n(), shape = 1.2, scale = 10)
    ),
    censor_time = runif(n(), min = 6, max = 24),
    event_status = if_else(os_months <= censor_time, 1L, 0L),
    os_months = pmin(os_months, censor_time)
  ) |>
  select(-censor_time)

soc_options <- c(
  "Gastrointestinal disorders",
  "Skin and subcutaneous tissue disorders",
  "General disorders and administration site conditions",
  "Nervous system disorders",
  "Musculoskeletal and connective tissue disorders",
  "Respiratory, thoracic and mediastinal disorders",
  "Metabolism and nutrition disorders"
)

ae_terms <- list(
  "Gastrointestinal disorders" = c("Nausea", "Vomiting", "Diarrhea", "Constipation"),
  "Skin and subcutaneous tissue disorders" = c("Rash", "Pruritus", "Dry skin"),
  "General disorders and administration site conditions" = c("Fatigue", "Pyrexia", "Injection site reaction"),
  "Nervous system disorders" = c("Headache", "Dizziness", "Peripheral neuropathy"),
  "Musculoskeletal and connective tissue disorders" = c("Arthralgia", "Myalgia", "Back pain"),
  "Respiratory, thoracic and mediastinal disorders" = c("Cough", "Dyspnea", "Pneumonitis"),
  "Metabolism and nutrition disorders" = c("Decreased appetite", "Hyponatremia", "Hyperglycemia")
)

ae_data <- map_dfr(1:n_patients, function(i) {
  n_ae <- sample(1:4, 1, prob = c(0.3, 0.35, 0.25, 0.1))
  socs <- sample(soc_options, n_ae, replace = TRUE)
  tibble(
    patient_id = trial_data$patient_id[i],
    treatment_arm = trial_data$treatment_arm[i],
    system_organ_class = socs,
    preferred_term = map_chr(socs, ~ sample(ae_terms[[.x]], 1)),
    ctcae_grade = sample(1:4, n_ae, replace = TRUE, prob = c(0.40, 0.35, 0.20, 0.05))
  )
})

ae_data <- ae_data |>
  bind_rows(
    trial_data |>
      filter(treatment_arm == "ONC-4587") |>
      slice_sample(n = 20) |>
      transmute(
        patient_id,
        treatment_arm,
        system_organ_class = "Skin and subcutaneous tissue disorders",
        preferred_term = sample(c("Rash", "Pruritus", "Dry skin"), 20, replace = TRUE),
        ctcae_grade = sample(1:3, 20, replace = TRUE, prob = c(0.4, 0.4, 0.2))
      )
  )

cat("=== Data Generation Complete ===\n")
cat(sprintf("Patients: %d | Events: %d | Censored: %d\n",
            nrow(trial_data),
            sum(trial_data$event_status == 1),
            sum(trial_data$event_status == 0)))

# ============================================================================
# SECTION 2: KAPLAN-MEIER SURVIVAL ANALYSIS
# ============================================================================

km_fit <- survfit(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
cat("\n=== Kaplan-Meier Median Survival ===\n")
print(km_fit)

logrank_test <- survdiff(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
logrank_p <- 1 - pchisq(logrank_test$chisq, df = 1)
cat(sprintf("\nLog-Rank Test p-value: %.4f\n", logrank_p))

km_plot <- ggsurvplot(
  km_fit,
  data = trial_data,
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  risk.table.col = "strata",
  palette = c("#E7B800", "#2E9FDF"),
  xlab = "Time (Months)",
  ylab = "Overall Survival Probability",
  title = "Kaplan-Meier Estimate of Overall Survival\nONC-4587 vs. Placebo in Advanced NSCLC",
  legend.title = "Treatment Arm",
  legend.labs = c("ONC-4587", "Placebo"),
  ggtheme = theme_minimal(base_size = 12),
  risk.table.height = 0.25,
  surv.median.line = "hv"
)

print(km_plot)
ggsave("km_survival_plot.png", plot = print(km_plot), width = 10, height = 7, dpi = 300)
cat("\nKM plot saved as 'km_survival_plot.png'\n")

# ============================================================================
# SECTION 3: COX PROPORTIONAL HAZARDS MODEL
# ============================================================================

cox_model <- coxph(
  Surv(os_months, event_status) ~ treatment_arm + age + sex + ecog_status,
  data = trial_data
)

cat("\n=== Cox Proportional Hazards Model Summary ===\n")
cox_summary <- summary(cox_model)
print(cox_summary)

hr <- exp(coef(cox_model)["treatment_armPlacebo"])
hr_ci <- exp(confint(cox_model)["treatment_armPlacebo", ])
hr_p <- cox_summary$coefficients["treatment_armPlacebo", "Pr(>|z|)"]

cat(sprintf("\n--- Primary Endpoint Result ---\n"))
cat(sprintf("Hazard Ratio (Placebo vs. ONC-4587): %.3f\n", hr))
cat(sprintf("95%% CI: [%.3f, %.3f]\n", hr_ci[1], hr_ci[2]))
cat(sprintf("p-value: %.4f\n", hr_p))

ph_test <- cox.zph(cox_model)
cat("\n=== Proportional Hazards Assumption Test ===\n")
print(ph_test)

tbl_cox <- cox_model |>
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      treatment_arm ~ "Treatment Arm",
      age ~ "Age (years)",
      sex ~ "Sex",
      ecog_status ~ "ECOG Performance Status"
    )
  ) |>
  bold_labels() |>
  bold_p(t = 0.05)

print(tbl_cox)

# ============================================================================
# SECTION 4: ADVERSE EVENT SUMMARY TABLE
# ============================================================================

cat("\n=== Adverse Event Summary ===\n")

# Summarize AEs by SOC and treatment arm (one AE per patient per SOC, max grade)
ae_summary <- ae_data |>
  group_by(patient_id, treatment_arm, system_organ_class) |>
  summarise(max_grade = max(ctcae_grade), .groups = "drop")

ae_table <- ae_summary |>
  mutate(
    grade_category = case_when(
      max_grade == 1 ~ "Grade 1",
      max_grade == 2 ~ "Grade 2",
      max_grade >= 3 ~ "Grade 3-4"
    )
  ) |>
  count(treatment_arm, system_organ_class, grade_category) |>
  pivot_wider(names_from = grade_category, values_from = n, values_fill = 0) |>
  arrange(treatment_arm, desc(`Grade 3-4`))

# gtsummary AE table
ae_for_tbl <- ae_summary |>
  mutate(
    any_ae = TRUE,
    grade_3_4 = max_grade >= 3
  )

ae_gtsummary <- ae_for_tbl |>
  select(treatment_arm, system_organ_class, grade_3_4) |>
  tbl_summary(
    by = treatment_arm,
    label = list(
      system_organ_class ~ "System Organ Class",
      grade_3_4 ~ "Grade >= 3 Event"
    )
  ) |>
  bold_labels() |>
  modify_header(label = "**Adverse Event Category**")

print(ae_gtsummary)

Your safety analysis is complete. Next up, you will tie the efficacy and safety results together into a programmatic executive summary that mirrors what biostatisticians deliver in a Clinical Study Report.

Write the Executive Summary

Your analysis is complete. You have a survival curve, a Cox model with Hazard Ratio, and a formatted Adverse Event table.

Now it is time to bring everything together into a single narrative. In this step, you will programmatically generate a Clinical Study Report (CSR) executive summary that pulls key statistics directly from your analysis objects. No copy-pasting numbers by hand.

In this step, get ready to:

• Build a programmatic executive summary using glue::glue() that references analysis objects directly.
• Print the summary to the console following CSR conventions (Background, Methods, Results, Conclusion).

Build the executive summary with glue::glue()

Real CSR executive summaries are written by medical writers who pull statistics from tables and insert them into a structured narrative. You will automate this by using glue::glue() to inject R expressions directly into a template string.

• In oncology_trial_analysis.R, scroll to the bottom of your script (below the SECTION 4 AE code).
• Add the following code to extract median survival times and build the executive summary string:

# ============================================================================
# SECTION 5: EXECUTIVE SUMMARY
# ============================================================================

median_os <- surv_median(km_fit)

exec_summary <- glue::glue("
============================================================
EXECUTIVE SUMMARY - CLINICAL STUDY REPORT
Study ONC-4587-PII-NSCLC
============================================================

BACKGROUND:
This Phase II, randomized, double-blind, placebo-controlled study
evaluated the efficacy and safety of ONC-4587 in patients with
advanced Non-Small Cell Lung Cancer (Stage IIIB/IV) who had
progressed on first-line platinum-based chemotherapy.

METHODS:
A total of {nrow(trial_data)} patients were randomized 1:1 to
ONC-4587 (n={sum(trial_data$treatment_arm == 'ONC-4587')}) or
Placebo (n={sum(trial_data$treatment_arm == 'Placebo')}).
The primary endpoint was Overall Survival (OS).
A Cox Proportional Hazards model was used for the primary analysis,
adjusted for age, sex, and ECOG performance status.

RESULTS:
Median OS was {round(median_os$median[1], 1)} months in the ONC-4587
arm vs. {round(median_os$median[2], 1)} months in the Placebo arm.
The Hazard Ratio for death (Placebo vs. ONC-4587) was {round(hr, 2)}
(95% CI: {round(hr_ci[1], 2)}-{round(hr_ci[2], 2)}; p={format.pval(hr_p, digits=3)}).
The log-rank test confirmed a statistically significant difference
between survival curves (p={format.pval(logrank_p, digits=3)}).

The proportional hazards assumption was {if_else(all(ph_test$table[, 'p'] > 0.05), 'met (all p > 0.05)', 'violated for at least one covariate')}.

SAFETY:
The most commonly reported adverse events in the ONC-4587 arm were
Gastrointestinal disorders and Skin and subcutaneous tissue disorders.
Rates of Grade 3-4 adverse events were {round(mean(ae_for_tbl$grade_3_4[ae_for_tbl$treatment_arm == 'ONC-4587']) * 100, 1)}%
in the experimental arm vs. {round(mean(ae_for_tbl$grade_3_4[ae_for_tbl$treatment_arm == 'Placebo']) * 100, 1)}%
in the placebo arm.

CONCLUSION:
ONC-4587 demonstrated a clinically meaningful and statistically
significant improvement in overall survival compared to placebo in
patients with advanced NSCLC. The safety profile was manageable,
with skin toxicity being the most notable class effect.
============================================================
")

What does this code do?

surv_median(km_fit) extracts the median overall survival time per treatment arm from the Kaplan-Meier fit. It returns a data frame with columns for strata, median, and confidence limits.

glue::glue() lets you embed R expressions inside curly braces {} within a string. At runtime, each expression evaluates and its result is pasted into the text. This means the summary always reflects the actual analysis output, not manually typed numbers.

• {nrow(trial_data)} inserts the total patient count.
• {round(hr, 2)} inserts the Hazard Ratio you extracted in the previous step.
• {format.pval(hr_p, digits=3)} formats the p-value following standard reporting conventions (e.g. displaying very small values as <2e-16).
• The if_else() block checks the PH assumption test results and inserts the appropriate interpretation.

Why automate the summary?

In pharmaceutical biostatistics, manually transcribing numbers from output into a report is a known source of errors. Regulatory agencies have rejected submissions over transcription mistakes. Programmatic report generation ensures every statistic traces directly back to the analysis code.

Print the summary to the console

The summary string is stored in exec_summary but has not been displayed yet. You need to explicitly output it using cat().

• Add the following print calls at the bottom of your script:

cat(exec_summary)

cat("\n\n=== Analysis Complete ===\n")

What does this code do?

cat(exec_summary) prints the fully interpolated executive summary to the R console. Unlike print(), cat() outputs raw text without surrounding quotes or escape characters. This gives you clean, readable output that matches what would appear in a real Clinical Study Report.

• Save oncology_trial_analysis.R.
• Select all the code in Section 5 (from the section header comment down to the final cat() call) and press Cmd+Enter (macOS) or Ctrl+Enter (Windows) to run it.
• In the console, you should see a formatted executive summary with four sections: BACKGROUND, METHODS, RESULTS, and CONCLUSION. The RESULTS section will display your actual median OS values, Hazard Ratio, confidence interval, and p-values pulled from the analysis objects.

Seeing an error about object not found?

The glue::glue() template references objects from earlier sections (hr, hr_ci, logrank_p, ph_test, ae_for_tbl). If any are missing, run the full script from the top by pressing Cmd+Shift+Enter (macOS) or Ctrl+Shift+Enter (Windows) to source the entire file.

help me debug an object not found error in my glue executive summary

✔️ Awesome, I've got everything!

Great. Double check you have saved oncology_trial_analysis.R and confirmed the executive summary prints correctly.

ⓧ I'd like to double check the full code

Here is the complete oncology_trial_analysis.R file including all sections from data generation through the executive summary:

# ============================================================================
# Phase II Oncology Clinical Trial Simulation and Analysis
# Drug: ONC-4587 (Experimental) vs. Placebo
# Indication: Advanced Non-Small Cell Lung Cancer (NSCLC)
# ============================================================================

# --- Load Libraries ---
library(tidyverse)
library(survival)
library(survminer)
library(gtsummary)

# --- Set Seed for Reproducibility ---
set.seed(42)

# ============================================================================
# SECTION 1: SYNTHETIC DATA GENERATION
# ============================================================================

n_patients <- 200

trial_data <- tibble(
  patient_id = sprintf("PT-%04d", 1:n_patients),
  treatment_arm = rep(c("ONC-4587", "Placebo"), each = n_patients / 2),
  age = round(rnorm(n_patients, mean = 62, sd = 9)),
  sex = sample(c("Male", "Female"), n_patients, replace = TRUE, prob = c(0.58, 0.42)),
  ecog_status = sample(0:2, n_patients, replace = TRUE, prob = c(0.3, 0.5, 0.2)),
  tumor_stage = sample(c("IIIB", "IV"), n_patients, replace = TRUE, prob = c(0.35, 0.65))
)

trial_data <- trial_data |>
  mutate(
    os_months = if_else(
      treatment_arm == "ONC-4587",
      rweibull(n(), shape = 1.2, scale = 14),
      rweibull(n(), shape = 1.2, scale = 10)
    ),
    censor_time = runif(n(), min = 6, max = 24),
    event_status = if_else(os_months <= censor_time, 1L, 0L),
    os_months = pmin(os_months, censor_time)
  ) |>
  select(-censor_time)

soc_options <- c(
  "Gastrointestinal disorders",
  "Skin and subcutaneous tissue disorders",
  "General disorders and administration site conditions",
  "Nervous system disorders",
  "Musculoskeletal and connective tissue disorders",
  "Respiratory, thoracic and mediastinal disorders",
  "Metabolism and nutrition disorders"
)

ae_terms <- list(
  "Gastrointestinal disorders" = c("Nausea", "Vomiting", "Diarrhea", "Constipation"),
  "Skin and subcutaneous tissue disorders" = c("Rash", "Pruritus", "Dry skin"),
  "General disorders and administration site conditions" = c("Fatigue", "Pyrexia", "Injection site reaction"),
  "Nervous system disorders" = c("Headache", "Dizziness", "Peripheral neuropathy"),
  "Musculoskeletal and connective tissue disorders" = c("Arthralgia", "Myalgia", "Back pain"),
  "Respiratory, thoracic and mediastinal disorders" = c("Cough", "Dyspnea", "Pneumonitis"),
  "Metabolism and nutrition disorders" = c("Decreased appetite", "Hyponatremia", "Hyperglycemia")
)

ae_data <- map_dfr(1:n_patients, function(i) {
  n_ae <- sample(1:4, 1, prob = c(0.3, 0.35, 0.25, 0.1))
  socs <- sample(soc_options, n_ae, replace = TRUE)
  tibble(
    patient_id = trial_data$patient_id[i],
    treatment_arm = trial_data$treatment_arm[i],
    system_organ_class = socs,
    preferred_term = map_chr(socs, ~ sample(ae_terms[[.x]], 1)),
    ctcae_grade = sample(1:4, n_ae, replace = TRUE, prob = c(0.40, 0.35, 0.20, 0.05))
  )
})

ae_data <- ae_data |>
  bind_rows(
    trial_data |>
      filter(treatment_arm == "ONC-4587") |>
      slice_sample(n = 20) |>
      transmute(
        patient_id,
        treatment_arm,
        system_organ_class = "Skin and subcutaneous tissue disorders",
        preferred_term = sample(c("Rash", "Pruritus", "Dry skin"), 20, replace = TRUE),
        ctcae_grade = sample(1:3, 20, replace = TRUE, prob = c(0.4, 0.4, 0.2))
      )
  )

cat("=== Data Generation Complete ===\n")
cat(sprintf("Patients: %d | Events: %d | Censored: %d\n",
            nrow(trial_data),
            sum(trial_data$event_status == 1),
            sum(trial_data$event_status == 0)))

# ============================================================================
# SECTION 2: KAPLAN-MEIER SURVIVAL ANALYSIS
# ============================================================================

km_fit <- survfit(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
cat("\n=== Kaplan-Meier Median Survival ===\n")
print(km_fit)

logrank_test <- survdiff(Surv(os_months, event_status) ~ treatment_arm, data = trial_data)
logrank_p <- 1 - pchisq(logrank_test$chisq, df = 1)
cat(sprintf("\nLog-Rank Test p-value: %.4f\n", logrank_p))

km_plot <- ggsurvplot(
  km_fit,
  data = trial_data,
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  risk.table.col = "strata",
  palette = c("#E7B800", "#2E9FDF"),
  xlab = "Time (Months)",
  ylab = "Overall Survival Probability",
  title = "Kaplan-Meier Estimate of Overall Survival\nONC-4587 vs. Placebo in Advanced NSCLC",
  legend.title = "Treatment Arm",
  legend.labs = c("ONC-4587", "Placebo"),
  ggtheme = theme_minimal(base_size = 12),
  risk.table.height = 0.25,
  surv.median.line = "hv"
)

print(km_plot)
ggsave("km_survival_plot.png", plot = print(km_plot), width = 10, height = 7, dpi = 300)
cat("\nKM plot saved as 'km_survival_plot.png'\n")

# ============================================================================
# SECTION 3: COX PROPORTIONAL HAZARDS MODEL
# ============================================================================

cox_model <- coxph(
  Surv(os_months, event_status) ~ treatment_arm + age + sex + ecog_status,
  data = trial_data
)

cat("\n=== Cox Proportional Hazards Model Summary ===\n")
cox_summary <- summary(cox_model)
print(cox_summary)

hr <- exp(coef(cox_model)["treatment_armPlacebo"])
hr_ci <- exp(confint(cox_model)["treatment_armPlacebo", ])
hr_p <- cox_summary$coefficients["treatment_armPlacebo", "Pr(>|z|)"]

cat(sprintf("\n--- Primary Endpoint Result ---\n"))
cat(sprintf("Hazard Ratio (Placebo vs. ONC-4587): %.3f\n", hr))
cat(sprintf("95%%%% CI: [%.3f, %.3f]\n", hr_ci[1], hr_ci[2]))
cat(sprintf("p-value: %.4f\n", hr_p))

ph_test <- cox.zph(cox_model)
cat("\n=== Proportional Hazards Assumption Test ===\n")
print(ph_test)

tbl_cox <- cox_model |>
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      treatment_arm ~ "Treatment Arm",
      age ~ "Age (years)",
      sex ~ "Sex",
      ecog_status ~ "ECOG Performance Status"
    )
  ) |>
  bold_labels() |>
  bold_p(t = 0.05)

print(tbl_cox)

# ============================================================================
# SECTION 4: ADVERSE EVENT SUMMARY TABLE
# ============================================================================

cat("\n=== Adverse Event Summary ===\n")

ae_summary <- ae_data |>
  group_by(patient_id, treatment_arm, system_organ_class) |>
  summarise(max_grade = max(ctcae_grade), .groups = "drop")

ae_table <- ae_summary |>
  mutate(
    grade_category = case_when(
      max_grade == 1 ~ "Grade 1",
      max_grade == 2 ~ "Grade 2",
      max_grade >= 3 ~ "Grade 3-4"
    )
  ) |>
  count(treatment_arm, system_organ_class, grade_category) |>
  pivot_wider(names_from = grade_category, values_from = n, values_fill = 0) |>
  arrange(treatment_arm, desc(`Grade 3-4`))

ae_for_tbl <- ae_summary |>
  mutate(
    any_ae = TRUE,
    grade_3_4 = max_grade >= 3
  )

ae_gtsummary <- ae_for_tbl |>
  select(treatment_arm, system_organ_class, grade_3_4) |>
  tbl_summary(
    by = treatment_arm,
    label = list(
      system_organ_class ~ "System Organ Class",
      grade_3_4 ~ "Grade >= 3 Event"
    )
  ) |>
  bold_labels() |>
  modify_header(label = "**Adverse Event Category**")

print(ae_gtsummary)

# ============================================================================
# SECTION 5: EXECUTIVE SUMMARY
# ============================================================================

median_os <- surv_median(km_fit)

exec_summary <- glue::glue("
============================================================
EXECUTIVE SUMMARY - CLINICAL STUDY REPORT
Study ONC-4587-PII-NSCLC
============================================================

BACKGROUND:
This Phase II, randomized, double-blind, placebo-controlled study
evaluated the efficacy and safety of ONC-4587 in patients with
advanced Non-Small Cell Lung Cancer (Stage IIIB/IV) who had
progressed on first-line platinum-based chemotherapy.

METHODS:
A total of {nrow(trial_data)} patients were randomized 1:1 to
ONC-4587 (n={sum(trial_data$treatment_arm == 'ONC-4587')}) or
Placebo (n={sum(trial_data$treatment_arm == 'Placebo')}).
The primary endpoint was Overall Survival (OS).
A Cox Proportional Hazards model was used for the primary analysis,
adjusted for age, sex, and ECOG performance status.

RESULTS:
Median OS was {round(median_os$median[1], 1)} months in the ONC-4587
arm vs. {round(median_os$median[2], 1)} months in the Placebo arm.
The Hazard Ratio for death (Placebo vs. ONC-4587) was {round(hr, 2)}
(95% CI: {round(hr_ci[1], 2)}-{round(hr_ci[2], 2)}; p={format.pval(hr_p, digits=3)}).
The log-rank test confirmed a statistically significant difference
between survival curves (p={format.pval(logrank_p, digits=3)}).

The proportional hazards assumption was {if_else(all(ph_test$table[, 'p'] > 0.05), 'met (all p > 0.05)', 'violated for at least one covariate')}.

SAFETY:
The most commonly reported adverse events in the ONC-4587 arm were
Gastrointestinal disorders and Skin and subcutaneous tissue disorders.
Rates of Grade 3-4 adverse events were {round(mean(ae_for_tbl$grade_3_4[ae_for_tbl$treatment_arm == 'ONC-4587']) * 100, 1)}%
in the experimental arm vs. {round(mean(ae_for_tbl$grade_3_4[ae_for_tbl$treatment_arm == 'Placebo']) * 100, 1)}%
in the placebo arm.

CONCLUSION:
ONC-4587 demonstrated a clinically meaningful and statistically
significant improvement in overall survival compared to placebo in
patients with advanced NSCLC. The safety profile was manageable,
with skin toxicity being the most notable class effect.
============================================================
")

cat(exec_summary)

cat("\n\n=== Analysis Complete ===\n")

Your analysis script is now complete. You have a full pipeline that generates data, performs survival analysis, fits a Cox model, summarizes adverse events, and produces a regulatory-style executive summary. Next, take on the Secret Mission to add a Forest Plot showing subgroup treatment effects.

Secret mission

Forest Plot of Subgroup Hazard Ratios

Your overall Hazard Ratio tells you ONC-4587 works on average. But does it work equally well for everyone? In regulatory submissions, agencies expect a subgroup forest plot showing whether age, sex, or performance status modifies the treatment effect. Build one from scratch using ggplot2.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs or cloud resources to worry about. Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

• oncology_trial_analysis.R script file.
• km_survival_plot.png and forest_plot_subgroups.png output files.
• R packages (tidyverse, survival, survminer, gtsummary).
• R and RStudio Desktop installations.

Keep everything running

No action needed. Choose this if you plan to keep working with R and clinical trial analysis.

• Your oncology_trial_analysis.R script and PNG outputs remain in your working directory, ready to re-run or modify.
• R, RStudio, and the installed packages stay on your machine for future projects.

Pause - I'll come back to this later

Close RStudio to free up memory. Your files and installations remain intact for when you return.

• Close RStudio Desktop.
• Your script file, PNG outputs, and all installed packages will still be there when you reopen RStudio.

Delete - I don't want to use this again

Remove all project resources and optionally uninstall R and RStudio.

• Delete the project files. In your terminal, navigate to the folder containing your project and run:

macOS

rm oncology_trial_analysis.R km_survival_plot.png forest_plot_subgroups.png

If you also want to remove the R packages, open R and run:

remove.packages(c("tidyverse", "survival", "survminer", "gtsummary"))

To uninstall R completely, delete the R framework and app:

sudo rm -rf /Library/Frameworks/R.framework
sudo rm -rf /Applications/R.app

To uninstall RStudio, drag it from your Applications folder to the Trash.

Windows

Remove-Item oncology_trial_analysis.R, km_survival_plot.png, forest_plot_subgroups.png

If you also want to remove the R packages, open R and run:

remove.packages(c("tidyverse", "survival", "survminer", "gtsummary"))

To uninstall R and RStudio:

• Open Settings.
• Select Apps, then Installed apps.
• Search for R for Windows and click Uninstall.
• Search for RStudio and click Uninstall.

Nice Work!

Nice Work!

You just built a clinical trial analysis pipeline that mirrors what biostatisticians deliver to regulatory agencies. From synthetic patient data to publication-ready outputs, you now have a complete Phase II oncology analysis in a single R script.

You've learned how to:

• Generated a realistic 200-patient Phase II oncology trial dataset with demographics, Weibull-distributed survival endpoints, right-censoring, and adverse event records.
• Performed Kaplan-Meier survival analysis with a log-rank test and fit a Cox Proportional Hazards model to extract a Hazard Ratio, 95% CI, and p-value for the treatment effect.
• Created regulatory-style deliverables including a publication-ready survival curve with risk table, a formatted Adverse Event summary table by System Organ Class and severity grade, and a programmatic Clinical Study Report executive summary.
• Built a Forest Plot showing subgroup Hazard Ratios by age, sex, and ECOG status to identify potential treatment effect modifiers.

Ready to quiz yourself?