R code examples

Simple randomisation with different treatments by domain

Simple randomisation function with 50% probability (1:1 ratio) to randomise to A or B in domain 1 and 80% vs. 20% probability (4:1 ratio) to randomise to A or C in domain 2. The domain comes from the randomisation form and can be accessed in the R code as randomisation$domain. The allocation sequence uses a seed to make it reproducible. Metadata is used to stored the seed state between function calls.

if (is.null(metadata$seed)) {
  # Set seed on first run
  set.seed(57897231)
} else {
  # Restore SEED state
  .Random.seed <- metadata$seed
}

if (randomisation$domain == "Domain 1") {
  # Domain 1
  group <- sample(c("A", "B"), size = 1, prob = c(0.5, 0.5))
} else {
  # Domain 2
  group <- sample(c("A", "C"), size = 1, prob = c(0.8, 0.2))
}

list(
  group = group,
  metadata = list(
    seed = .Random.seed
  )
)

Stratified random permuted block lists

This method uses the rpbrPar method from the randomizeR package to create a blocked list for each stratum and keep track of the location in the list to return the next available treatment.

library("randomizeR")
library("digest")

# List length per stratum - recommended to make this same as trial sample size
# so that we don't run out of allocations during the trial. Also note that if
# you change this mid-trial you will start a different blocked list!
N_PER_STRATUM = 200

# Block sizes
BLOCK_SIZES = c(4, 6)

# Treatment groups
GROUPS = c("Intervention", "Control")

# Treatment group ratios
RATIO = c(1, 1)

# Stratification factors
STRATA = c(randomisation$Sex, randomisation$AgeGroup) 

# Initialise metadata
if (!is.list(metadata$counts)) {
  message("initialise metadata")
  metadata = list("counts" = list(), "seeds" = list())
}

# Keep track of position in randomisation list for each stratum
get_list_position <- function(stratum, metadata) {
  if (!(stratum %in% names(metadata$counts))) {
    # First randomisation in this stratum
    message("stratum not in metadata")
    metadata$counts[stratum] = 1
  }
  return(metadata$counts[[stratum]])
}

# Initialise if new stratum
initialise_seed <- function(stratum, metadata) {
  if (!(stratum %in% names(metadata$seeds))) {
    # First randomisation in this stratum
    message("seed not in metadata")
    metadata$seeds[stratum] = sample(1:99999999, 1)
  }
  return(metadata)
}

increment_list_position <- function(stratum, metadata) {
  metadata$counts[stratum] <- get_list_position(stratum, metadata) + 1
  return(metadata)
}

randomise <- function(stratum, metadata) {
  n <- get_list_position(stratum, metadata)
  message("List position ", n)
  r <- rpbrPar(N_PER_STRATUM, BLOCK_SIZES, K = 2, ratio = RATIO, filledBlock = FALSE)
  metadata <- increment_list_position(stratum, metadata)
  metadata <- initialise_seed(stratum, metadata)
  seed <- metadata$seeds[[stratum]]
  message("Using seed ", seed)
  s <- genSeq(r, 1, seed)$M[,n] + 1
  return(list(group = GROUPS[s], metadata = metadata))
}

stratum <- digest(STRATA, algo="sha256")
message(paste("Randomising in strata:", toString(STRATA), "=", stratum))
randomise(stratum, metadata)

Response adaptive randomisation

Response adaptive randomisation function to group A or B, with c = n/2N. Based on Practical Bayesian Adaptive Randomization in Clinical Trials (Peter F Thall, J Kyle Wathen).

Outcome data comes from a follow-up form (a Yes/No field called Success) and is passed into the R-code as auxiliary_data$outcomes. All previous randomisations (containing the group  field) are passed in as auxiliary_data$randomisations .

The message function is used to log comments, so they can be viewed on the View R randomisation log page.

rar <- function(outcome_data, randomisation, randomisations, metadata, sample_size) {
  # --- Calculate Success Rates (Response-Adaptive Component) ---
  
  # Handle the initial state where there is no outcome data yet.
  # We assume a neutral success rate of 0.5 for each group to start.
  if (nrow(outcome_data) == 0) {
    sr_a <- 0.5
    sr_b <- 0.5
  } else {
    # Subset outcomes for each group
    outcomes_a <- outcome_data[outcome_data$group == "A", ]
    outcomes_b <- outcome_data[outcome_data$group == "B", ]
    
    # Calculate success rate for Group A, using 0.5 if no data exists
    if (nrow(outcomes_a) == 0) {
      sr_a <- 0.5
    } else {
      sr_a <- sum(outcomes_a$Success) / nrow(outcomes_a)
    }
    
    # Calculate success rate for Group B, using 0.5 if no data exists
    if (nrow(outcomes_b) == 0) {
      sr_b <- 0.5
    } else {
      sr_b <- sum(outcomes_b$Success) / nrow(outcomes_b)
    }
  }
  
  # Normalize success rates into probabilities. Handle case where both are 0.
  total_sr <- sr_a + sr_b
  if (total_sr == 0) {
    prob_a_adaptive <- 0.5
  } else {
    prob_a_adaptive <- sr_a / total_sr
  }
  prob_b_adaptive <- 1 - prob_a_adaptive
  
  # --- Determine Final Probabilities ---
      
  # Use the adaptive probabilities with c=n/2N
  c = nrow(randomisations) / (2 * sample_size)
  final_prob_a <- prob_a_adaptive^c / (prob_a_adaptive^c + prob_b_adaptive^c)
  final_prob_b <- 1 - final_prob_a
  message("ID: ", randomisation$id)
  message("Using adaptive probabilities with c = ", c)
  message(
    "Success Rates -> A: ", round(sr_a, 2), ", B: ", round(sr_b, 2),
    " | Final Probs -> A: ", final_prob_a, ", B: ", final_prob_b 
  )
 
  # --- Randomize and Return Assignment ---
  
  # Use the final probabilities to randomly sample the assigned group
  group <- sample(c("A", "B"), size = 1, prob = c(final_prob_a, final_prob_b))
  
  return(list(group = group, metadata = metadata))
}

# Convert Yes/No to booleans
auxiliary_data$outcomes$Success = auxiliary_data$outcomes$Success == "Yes"

# Convert auxiliary data into format required by the function
if (is.data.frame(auxiliary_data$randomisations)) {
  randomisations <- auxiliary_data$randomisations
} else {
  randomisations <- data.frame(
    subjectId = integer(),
    group = character()
  )
}

if (is.data.frame(auxiliary_data$outcomes) && is.data.frame(auxiliary_data$randomisations)) {
  historical_data <- merge(
    auxiliary_data$outcomes, 
    auxiliary_data$randomisations, 
    by = "subjectId"
  )
  
  outcome_data <- data.frame(
    subjectId = historical_data$subjectId,
    group = historical_data$group,
    Success = historical_data$Success
  )
} else {
  outcome_data <- data.frame(
    subjectId = integer(),
    group = character(),
    Success = logical()
  )
}

rar(
  outcome_data,
  randomisation,
  randomisations,
  metadata,
  100 
)