Clinical Trials

---
name: clinical-trials
description: Clinical trial design and analysis methods in R, including randomization, estimands, multiplicity, and reporting.
---

# Clinical Trials Statistical Methods

## Overview

Comprehensive clinical trial design and analysis methods in R covering sample size calculation, randomization, interim analyses, multiplicity adjustment, and regulatory-compliant statistical methods.

## Sample Size Calculation

### Two-Group Comparisons

```r
library(pwr)

# Two-sample t-test
pwr.t.test(
  d = 0.5,           # Cohen's d effect size
  sig.level = 0.05,
  power = 0.80,
  type = "two.sample",
  alternative = "two.sided"
)

# Proportions (chi-square)
pwr.2p.test(
  h = ES.h(p1 = 0.6, p2 = 0.4),  # Cohen's h
  sig.level = 0.05,
  power = 0.80
)

# Two proportions (unequal groups)
pwr.2p2n.test(
  h = ES.h(p1 = 0.6, p2 = 0.4),
  n1 = 100,
  sig.level = 0.05
)
```

### Survival Endpoints

```r
library(gsDesign)

# Log-rank test sample size
nSurv(
  lambda1 = log(2)/12,  # Control median = 12 months
  lambda2 = log(2)/18,  # Treatment median = 18 months (HR = 0.67)
  Ts = 24,              # Study duration
  Tr = 12,              # Accrual duration
  alpha = 0.025,        # One-sided

  beta = 0.20,          # 80% power
  ratio = 1             # 1:1 randomization
)

# Using rpact
library(rpact)
getSampleSizeSurvival(
  hazardRatio = 0.67,
  lambda1 = log(2)/12,
  accrualTime = 12,
  followUpTime = 12,
  alpha = 0.025,
  beta = 0.20,
  allocationRatioPlanned = 1
)
```

### Non-Inferiority Trials

```r
library(TrialSize)

# Non-inferiority for proportions
TwoSampleProportion.NIS(
  p = 0.80,           # Expected proportion in both groups
  delta = 0.10,       # Non-inferiority margin
  alpha = 0.025,      # One-sided
  power = 0.80
)

# Non-inferiority for means
TwoSampleMean.NIS(
  sigma = 10,         # SD
  delta = 5,          # Non-inferiority margin
  alpha = 0.025,
  power = 0.80
)
```

## Randomization

### Simple Randomization

```r
# Base R simple randomization
set.seed(123)
n <- 100
treatment <- sample(c("A", "B"), n, replace = TRUE)

# blockrand package
library(blockrand)
randomization <- blockrand(
  n = 100,
  num.levels = 2,
  levels = c("Treatment", "Control"),
  id.prefix = "PAT",
  block.prefix = "BLK"
)
```

### Stratified Block Randomization

```r
library(blockrand)

# Generate lists for each stratum
strata <- expand.grid(
  sex = c("Male", "Female"),
  age_group = c("<65", ">=65")
)

rand_lists <- lapply(1:nrow(strata), function(i) {
  blockrand(
    n = 50,
    num.levels = 2,
    levels = c("Treatment", "Control"),
    block.sizes = c(2, 4, 6),  # Variable block sizes
    stratum = paste(strata[i, ], collapse = "_")
  )
})

full_list <- do.call(rbind, rand_lists)
```

### Minimization

```r
library(Minirand)

# Minimization randomization
minirand(
  covariates = data.frame(
    sex = c("M", "F", "M"),
    age = c("young", "old", "young"),
    center = c("A", "A", "B")
  ),
  treatment = c("A", "B"),
  ratio = c(1, 1),
  p = 0.85  # Probability of assigning to minimizing treatment
)
```

## Baseline Characteristics

For randomized trials, baseline tables should primarily describe the randomized groups and assess clinically meaningful imbalance. Routine baseline hypothesis tests are usually not appropriate because any baseline differences arise after randomization and p-values mostly reflect sample size.

```r
library(gtsummary)

baseline_table <- adsl |>
  dplyr::filter(ITTFL == "Y") |>
  dplyr::select(TRT01P, AGE, SEX, BMI, BASELINE_SCORE) |>
  tbl_summary(
    by = TRT01P,
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    )
  ) |>
  modify_header(label = "**Characteristic**")
```

If the SAP requires standardized mean differences, report them as descriptive diagnostics rather than randomization tests.

## Group Sequential Designs

### gsDesign Package

```r
library(gsDesign)

# O'Brien-Fleming boundaries
gs_design <- gsDesign(
  k = 3,               # Number of analyses
  test.type = 2,       # Two-sided symmetric
  alpha = 0.025,       # One-sided alpha
  beta = 0.20,         # Type II error
  sfu = "OF",          # O'Brien-Fleming spending function
  timing = c(0.5, 0.75, 1)  # Information fractions
)

# Summary
gs_design

# Boundaries
gs_design$upper$bound  # Upper efficacy boundary
gs_design$lower$bound  # Lower futility boundary

# Plot
plot(gs_design)
```

### Alpha Spending Functions

```r
# Pocock
gs_pocock <- gsDesign(k = 3, sfu = "Pocock")

# Hwang-Shih-DeCani
gs_hsd <- gsDesign(k = 3, sfu = sfHSD, sfupar = -4)

# Power family (Kim-DeMets)
gs_power <- gsDesign(k = 3, sfu = sfPower, sfupar = 2)

# Custom spending
gs_custom <- gsDesign(
  k = 3,
  sfu = sfPoints,
  sfupar = c(0.01, 0.03, 0.025),  # Cumulative alpha at each look
  timing = c(0.5, 0.75, 1)
)
```

### rpact Package

```r
library(rpact)

# Design
design <- getDesignGroupSequential(
  kMax = 3,
  alpha = 0.025,
  beta = 0.20,
  sided = 1,
  typeOfDesign = "OF",  # O'Brien-Fleming
  informationRates = c(0.5, 0.75, 1)
)

# Sample size
sampleSize <- getSampleSizeMeans(
  design = design,
  alternative = 0.5,
  stDev = 1
)

# Interim analysis
getAnalysisResults(
  design,
  dataInput = getDataset(
    n = c(50, 75),
    means = c(0.3, 0.4),
    stDevs = c(1, 1)
  )
)
```

## Multiplicity Adjustment

### P-value Adjustments

```r
# Bonferroni
p.adjust(p_values, method = "bonferroni")

# Holm (step-down)
p.adjust(p_values, method = "holm")

# Hochberg (step-up)
p.adjust(p_values, method = "hochberg")

# Benjamini-Hochberg (FDR)
p.adjust(p_values, method = "BH")

# Hommel
p.adjust(p_values, method = "hommel")
```

### Graphical Approaches

```r
library(gMCP)

# Define hypothesis graph
graph <- matrix2graph(
  # Transition matrix
  m = matrix(c(
    0, 0.5, 0.5, 0,
    0.5, 0, 0, 0.5,
    0.5, 0, 0, 0.5,
    0, 0.5, 0.5, 0
  ), nrow = 4, byrow = TRUE),
  # Initial weights
  w = c(0.5, 0.5, 0, 0)
)

# Set hypothesis names
nodeNames(graph) <- c("H1_OS", "H2_OS", "H1_PFS", "H2_PFS")

# Plot graph
plot(graph)

# Perform test
gMCP(
  graph = graph,
  pvalues = c(0.01, 0.03, 0.02, 0.04),
  alpha = 0.025
)
```

### Gatekeeping Procedures

```r
library(multcomp)

# Serial gatekeeping
# Primary must be significant before testing secondary
serial_gate <- function(p_primary, p_secondary, alpha = 0.05) {
  if (p_primary < alpha) {
    # Primary significant, test secondary at full alpha
    return(c(primary = p_primary < alpha, secondary = p_secondary < alpha))
  } else {
    return(c(primary = FALSE, secondary = FALSE))
  }
}
```

## Missing Data

### Mixed Models for Repeated Measures (MMRM)

```r
library(mmrm)

# MMRM model
mmrm_fit <- mmrm(
  formula = change ~ treatment * visit + baseline + us(visit | subject),
  data = long_data,
  weights = NULL,
  reml = TRUE
)

# Least squares means
library(emmeans)
emmeans(mmrm_fit, ~ treatment | visit)

# Treatment comparison at each visit
emmeans(mmrm_fit, pairwise ~ treatment | visit)
```

### Multiple Imputation

```r
library(mice)

# Create imputations
imp <- mice(
  data = df,
  m = 20,           # Number of imputations
  method = "pmm",   # Predictive mean matching
  maxit = 10
)

# Analyze each imputed dataset
analyses <- with(imp, lm(outcome ~ treatment + covariates))

# Pool results (Rubin's rules)
pooled <- pool(analyses)
summary(pooled)
```

### Tipping Point Analysis

```r
# Sensitivity analysis for MNAR
library(rbmi)

# Define imputation method with delta adjustment
draws <- draws(
  data = data,
  data_ice = ice_data,
  method = method_bayes(),
  vars = vars
)

# Impute with different delta values
impute(draws, references = c("Control" = "Control", "Treatment" = "Control"))
```

## Subgroup Analysis

### Forest Plots for Subgroups

```r
library(forestplot)

# Calculate treatment effects by subgroup
subgroup_effects <- df |>
  group_by(subgroup) |>
  summarise(
    n = n(),
    effect = mean(outcome[trt == 1]) - mean(outcome[trt == 0]),
    se = sqrt(var(outcome[trt == 1])/sum(trt == 1) +
              var(outcome[trt == 0])/sum(trt == 0)),
    lower = effect - 1.96 * se,
    upper = effect + 1.96 * se
  )

# Create forest plot
forestplot(
  labeltext = subgroup_effects$subgroup,
  mean = subgroup_effects$effect,
  lower = subgroup_effects$lower,
  upper = subgroup_effects$upper,
  zero = 0,
  xlab = "Treatment Effect (95% CI)"
)
```

### Interaction Tests

```r
# Test for treatment-by-subgroup interaction
interaction_model <- lm(outcome ~ treatment * subgroup, data = df)
anova(interaction_model)

# Quantitative interaction test
library(QI)
qi_test(outcome ~ treatment | subgroup, data = df)
```

## Regulatory Considerations

### ICH E9 Estimands Framework

```r
# Define estimand components:
# 1. Treatment condition
# 2. Population
# 3. Variable (endpoint)
# 4. Intercurrent events and strategies
# 5. Population-level summary measure

# Example: Treatment policy estimand with MMRM
mmrm_fit <- mmrm(
  change ~ treatment * visit + baseline + us(visit | subject),
  data = data_all_randomized  # Include all randomized (ITT)
)
```

Missing-data methods should follow the estimand and the plausible missingness mechanism. Prespecify primary handling and sensitivity analyses rather than defaulting to complete-case analysis.

### CONSORT Diagram

```r
library(consort)

# Create CONSORT diagram
consort_plot(
  data = trial_data,
  orders = c(
    Assessed = "Assessed for eligibility",
    Randomized = "Randomized",
    Arm_A = "Allocated to Arm A",
    Arm_B = "Allocated to Arm B",
    Lost_A = "Lost to follow-up (Arm A)",
    Lost_B = "Lost to follow-up (Arm B)",
    Analyzed_A = "Analyzed (Arm A)",
    Analyzed_B = "Analyzed (Arm B)"
  ),
  side_box = c("Excluded", "Discontinued_A", "Discontinued_B"),
  cex = 0.8
)
```

## Key Packages Summary

| Package | Purpose |
|---------|---------|
| pwr | Power analysis |
| gsDesign | Group sequential designs |
| rpact | Adaptive designs |
| blockrand | Randomization |
| gMCP | Graphical multiplicity |
| mmrm | MMRM analysis |
| mice | Multiple imputation |
| rbmi | Reference-based imputation |
| emmeans | Least squares means |
| consort | CONSORT diagrams |