# Recreate the simulated participant-level dataset used in Part 10 of the
# NajibMY Biostatistics Learning Hub. The script uses base R only.

make_survival_data <- function(n = 240) {
  id <- seq_len(n)
  group <- as.integer(id %% 2 == 0)
  group_label <- ifelse(group == 0, "Usual care", "Intervention")
  age <- 35 + ((7 * id) %% 41)
  stage3 <- as.integer(id %% 3 == 0)

  deterministic_uniform <- function(multiplier) {
    (((id * multiplier) %% 997) + 0.5) / 998
  }

  u_recurrence <- deterministic_uniform(37)
  u_nonrecurrence_death <- deterministic_uniform(73)
  u_postrecurrence_death <- deterministic_uniform(109)
  u_dropout <- deterministic_uniform(149)
  u_dropout_time <- deterministic_uniform(181)

  recurrence_rate <- 0.035 * exp(
    log(0.55) * group + 0.025 * (age - 55) + 0.55 * stage3
  )
  nonrecurrence_death_rate <- 0.010 * exp(
    log(0.85) * group + 0.030 * (age - 55) + 0.40 * stage3
  )
  postrecurrence_death_rate <- 0.060 * exp(
    log(0.85) * group + 0.020 * (age - 55) + 0.25 * stage3
  )

  recurrence_time_latent <- -log(u_recurrence) / recurrence_rate
  nonrecurrence_death_time_latent <-
    -log(u_nonrecurrence_death) / nonrecurrence_death_rate
  postrecurrence_death_delay <-
    -log(u_postrecurrence_death) / postrecurrence_death_rate

  overall_death_time <- pmin(
    nonrecurrence_death_time_latent,
    recurrence_time_latent + postrecurrence_death_delay
  )
  dropout_time <- ifelse(
    u_dropout < 0.22,
    12 + 24 * u_dropout_time,
    Inf
  )
  censor_time <- pmin(36, dropout_time)
  overall_time <- pmin(overall_death_time, censor_time)
  overall_event <- as.integer(overall_death_time <= censor_time)
  competing_time <- pmin(
    recurrence_time_latent,
    nonrecurrence_death_time_latent,
    censor_time
  )
  competing_status <- ifelse(
    recurrence_time_latent <= nonrecurrence_death_time_latent &
      recurrence_time_latent <= censor_time,
    1L,
    ifelse(
      nonrecurrence_death_time_latent < recurrence_time_latent &
        nonrecurrence_death_time_latent <= censor_time,
      2L,
      0L
    )
  )

  data.frame(
    id,
    group,
    group_label,
    age,
    stage3,
    overall_time,
    overall_event,
    competing_time,
    competing_status,
    recurrence_time_latent,
    nonrecurrence_death_time_latent,
    censor_time,
    check.names = FALSE
  )
}

arguments <- commandArgs(trailingOnly = TRUE)
output_file <- if (length(arguments) >= 1 && nzchar(arguments[[1]])) {
  arguments[[1]]
} else {
  "part-10-simulated-survival-data.csv"
}

output_directory <- dirname(output_file)
if (!dir.exists(output_directory)) {
  dir.create(output_directory, recursive = TRUE)
}

write.csv(
  make_survival_data(),
  file = output_file,
  row.names = FALSE,
  na = ""
)

message("Wrote 240 simulated records to: ", normalizePath(output_file))