# Simulate a fictitious exposure-response data set for a pharmacometrics
# analysis, used to illustrate problems with dplyr::ntile() when there are
# ties at the median of a variable recorded with limited (1 decimal place)
# precision -- body weight.
#
# Three studies:
#   STUDY001 - Phase 1 SAD dose-escalation, all male
#   STUDY002 - Phase 1 DDI with an oral contraceptive, all female, single dose
#   STUDY003 - Phase 2 dose-finding, mixed sex
#
# Exposure (AUC, Cmax) comes from a simple one-compartment oral PK model with
# allometric weight scaling and log-normal inter-individual variability. A
# continuous response variable is generated with a simple linear relationship
# to AUC plus noise.

library(tidyverse)

set.seed(8342)

# ---- study designs -------------------------------------------------------

studies <- tribble(
  ~study_id,  ~phase, ~description,
  "S01",        "1",    "SAD dose-escalation",
  "S02",        "1",    "DDI oral contraceptive",
  "S03",        "2",    "Dose-finding"
)

study001 <- tibble(
  study_id = "S01",
  dose_mg = rep(c(10, 30, 100, 300), each = 8),
  sex = "M"
)

study002 <- tibble(
  study_id = "S02",
  dose_mg = 100,
  sex = "F"
) |> slice(rep(1, 24))

study003 <- tibble(
  study_id = "S03",
  dose_mg = rep(c(50, 100, 200), each = 40)
) |> mutate(sex = rep(c("M", "F"), length.out = n()))

subjects <- bind_rows(study001, study002, study003) |>
  left_join(studies, by = "study_id") |>
  # drop one subject so the total N is odd and the median weight is a single
  # observed value rather than an average of two
  slice(-1) |>
  mutate(subject_id = sprintf("%s-%03d", study_id, row_number()), .after = study_id)

# ---- body weight ----------------------------------------------------------

# Study/sex-specific true (continuous) weight distributions. Means differ
# across studies but ranges overlap substantially in the middle, so that
# after rounding, ties at the median can plausibly come from more than one
# study.
weight_params <- tribble(
  ~study_id,  ~sex, ~mean_wt, ~sd_wt,
  "S01",      "M",  83,       9,
  "S02",      "F",  67,       8,
  "S03",      "M",  82,       10,
  "S03",      "F",  68,       9
)

subjects <- subjects |>
  left_join(weight_params, by = c("study_id", "sex")) |>
  mutate(
    wt_kg_true = rnorm(n(), mean_wt, sd_wt),
    wt_kg = round(wt_kg_true, 1)
  ) |>
  select(-mean_wt, -sd_wt)

# ---- engineer ties at the median weight ------------------------------------

# Take the current median of the rounded weight, then nudge a handful of
# near-median subjects -- deliberately drawn from more than one study -- so
# their true weight falls in the rounding bin that lands exactly on that
# median value. This keeps wt_kg_true continuous/distinct while forcing
# a cluster of ties in wt_kg, with more than one study represented among
# the ties.
median_target <- median(subjects$wt_kg)
n_tie_target <- 14 # aim for roughly this many ties at the median, spread across studies

candidates <- subjects |>
  mutate(dist = abs(wt_kg_true - median_target)) |>
  arrange(dist) |>
  # take the nearest candidates from each study in turn (round-robin) so the
  # tie is not dominated by a single study
  mutate(rank_within_study = row_number(), .by = study_id) |>
  arrange(rank_within_study, dist) |>
  slice(1:n_tie_target)

subjects <- subjects |>
  mutate(
    nudge = subject_id %in% candidates$subject_id,
    # jitter within +/- 0.04 so all of these still round to median_target
    wt_kg_true = if_else(
      nudge,
      median_target + runif(n(), -0.04, 0.04),
      wt_kg_true
    ),
    wt_kg = round(wt_kg_true, 1)
  ) |>
  select(-nudge)

# ---- exposure (AUC, Cmax) from a simple 1-compartment oral PK model -------

theta_cl <- 5  # typical CL/F, L/h
theta_v <- 50  # typical V/F, L
ka <- 1        # absorption rate constant, 1/h

subjects <- subjects |>
  mutate(
    eta_cl = rnorm(n(), 0, 0.3),
    eta_v = rnorm(n(), 0, 0.3),
    cl = theta_cl * (wt_kg_true / 70)^0.75 * exp(eta_cl),
    v = theta_v * (wt_kg_true / 70) * exp(eta_v),
    ke = cl / v,
    tmax = log(ka / ke) / (ka - ke),
    auc_pred = dose_mg / cl,
    cmax_pred = (dose_mg / v) * (ka / (ka - ke)) * (exp(-ke * tmax) - exp(-ka * tmax)),
    # residual (assay/formulation) variability on top of the model prediction
    auc = auc_pred * exp(rnorm(n(), 0, 0.2)),
    cmax = cmax_pred * exp(rnorm(n(), 0, 0.2))
  ) |>
  select(-eta_cl, -eta_v, -cl, -v, -ke, -tmax, -auc_pred, -cmax_pred)

# ---- continuous response, linear in AUC ------------------------------------

beta0 <- 20
beta1 <- 1.2
sigma_response <- 8

subjects <- subjects |>
  mutate(response = beta0 + beta1 * auc + rnorm(n(), 0, sigma_response))

# ---- tidy up and save -------------------------------------------------------

exposure_response_data <- subjects |>
  arrange(study_id, dose_mg, subject_id) |>
  mutate(row_id = row_number()) |> 
  select(
    row_id, study_id, phase, description, subject_id, sex, 
    dose_mg, wt_kg_true, wt_kg, auc, cmax, response
  )

write_csv(
  exposure_response_data,
  here::here("posts", "2026-09-20_the-trouble-with-ntile", "exposure_response_data.csv")
)

exposure_response_data