Time to Dermatologic Event by Treatment Group

id: f-derm-time-to-event
title: "Time to Dermatologic Event by Treatment Group"
regulatory_id: AEF01
type: figure
# ADTTE is the study's own time-to-event dataset, and only the CDISC pilot's own
# ADaM package ships it. `sources: phuse` takes the whole study from that
# packaging so the denominators, the arms and the event times come from one
# coherent submission rather than two — the pilot's ADSL puts 86 / 84 / 84
# subjects in the three arms, which is what the reference report reports.
sources: phuse
dataset: adtte
analysis_set: safety
denominator: adsl
group:
  - TRT01A
total: false
analyses:
  - name: km
    method: figure
    custom: km_time_to_event
    label: "Kaplan-Meier estimates of time to first dermatologic event"
    # define.xml, Figure_14-1: ParamCD "TTDE", AnalysisVariables AVAL and CNSR,
    # `time aval*cnsr(1)` — CNSR = 1 is the censoring value.
    param: TTDE
    time_var: AVAL
    censor_var: CNSR
    censor_value: 1
    # The grid for the numbers-at-risk strip. Declared here so the counts the
    # figure prints are committed statistics in the ARD, checkable against it,
    # rather than arithmetic the renderer does on its way to the page.
    risk_times: [0, 50, 100, 150, 200]

id: f-derm-time-to-event
title: "Time to Dermatologic Event by Treatment Group"
study: CDISCPILOT01
population_label: "Safety Analysis Set"
cutoff: "2014-07-01"
source: "Source: adtte, adsl (CDISC pilot ADaM package, vendored from phuse-org/phuse-scripts). Data cut-off: 2014-07-01."
footnotes:
  - "Dermatologic events are the adverse events the study flagged as its first customised query, CQ01NAM = 'DERMATOLOGIC EVENTS'; ADTTE carries one record per subject for the first such treatment-emergent event."
  - "Subjects without a dermatologic event are censored at their study completion date (CNSR = 1); the event time is the day of the first event (CNSR = 0)."
  - "Survival probability is the Kaplan-Meier product-limit estimate of remaining free of a dermatologic event. Tick marks on a curve are censored subjects, and the counts beneath it are the subjects still at risk at each of those days."
  - "The median and its 95% confidence limits use the linear transformation, reproducing the limits in the reference report; NE means the curve never reached 0.5 and the median is not estimable."
  - "The test annotated on the figure is the log-rank test of equality across the three treatment groups. Its chi-square and degrees of freedom are stated alongside the p-value so the p-value can be checked rather than taken on trust; the reference report states only that the difference was significant. A p-value smaller than the precision this display declares is reported at that boundary — '<0.0001' — rather than rounded to '0.0000', which would assert a probability of zero. The unrounded value is retained in the analysis results dataset."
  - "The curve, the tick marks, the numbers at risk and the annotated test are drawn from the same committed analysis results dataset as the table beneath them, so no number on this page is a second calculation of another."
columns:
  order: [Placebo, "Xanomeline Low Dose", "Xanomeline High Dose"]
format:
  digits:
    p: 0
    n_censor: 0
    p_censor: 0
    median: 0
    median_lcl: 0
    median_ucl: 0
    chisq: 2
    df: 0
    p_value: 4
  # `p_censor` is a proportion this display's custom.R computes, so it is
  # declared here rather than pre-scaled in the ARD (contract §3).
  proportions: [p, prop, p_censor]
  median_ci: "{median} ({median_lcl}, {median_ucl})"
rows:
  - { analysis: km, pattern: n_pct, label: "Subjects with a dermatologic event, n (%)" }
  - { analysis: km, pattern: "{n_censor} ({p_censor}%)", label: "Censored, n (%)" }
  - { analysis: km, pattern: median_ci, label: "Median time to event, days (95% CI)",
      na_text: "NE" }
figure:
  analysis: km
  plot_title: "Product-Limit Survival Function Estimates"
  width: 900
  height: 460
  x_axis:
    label: "Time to Dermatologic Event or End of Study (days)"
    min: 0
    max: 200
    ticks: [0, 50, 100, 150, 200]
  y_axis:
    label: "Survival Probability"
    min: 0
    max: 1
    ticks: [0, 0.2, 0.4, 0.6, 0.8, 1.0]
    digits: 1
  stats:
    time: time
    value: surv
    censor_time: censor_time
    censor_value: censor_surv
  # Series are named and ordered here; their colours and dash patterns come from
  # the engine's palette, which is colour-vision-safe and has a dark-scheme
  # counterpart. A hex value written here would be presentation masquerading as
  # part of the display's definition.
  series:
    - level: Placebo
    - level: "Xanomeline Low Dose"
    - level: "Xanomeline High Dose"
  at_risk:
    times: risk_time
    counts: risk_n
    label: "Number at risk"
  annotation:
    bindings:
      chisq: "km:chisq"
      df: "km:df"
      p: "km:p_value"
    template: "Log-rank chi-square {chisq} on {df} df, p {p}"
variants:
  post_text: {}

# Kaplan-Meier statistics for Figure 14-1 of the CDISCPILOT01 clinical study
# report, "Time to Dermatologic Event by Treatment Group".
#
# The analysis is not invented here. It is stated twice in the study's own
# submission package, and this file implements exactly what those two documents
# say and nothing more:
#
#   define.xml, Analysis Results Metadata, ResultDisplay "Figure_14-1"
#     Documentation ....... "Kaplan-Meier estimates and log-rank analysis of
#                            time to first dermatological adverse event, safety
#                            population"
#     AnalysisDataset ..... ADTTE where SAFFL = "Y"
#     ParameterList ....... PARAMCD = "TTDE"
#     AnalysisVariables ... ADTTE.AVAL, ADTTE.CNSR
#     ProgrammingCode ..... proc lifetest data = adtte (where=(saffl="Y")) plots=s;
#                             id usubjid;
#                             strata trtan;
#                             time aval*cnsr(1);
#                             test trtan;
#                           run;
#
#   Statistical analysis plan, section 11.2
#     "The time to the first dermatological event will be compared across the
#      treatment groups using Kaplan-Meier methods. Graphical displays of the
#      survival curves will be presented."
#
# Two decisions follow from the SAS code and are worth stating, because both are
# reproductions of a 2006 result rather than modern defaults:
#
#   `time aval*cnsr(1)` makes CNSR = 1 the censoring value, so an event is
#   CNSR = 0. ADTTE's define.xml derivation confirms it: "if ADAE.TRTEMFL = 'Y'
#   then CNSR = 0, else CNSR = 1".
#
#   The median's confidence limits use the linear ("plain") transformation,
#   which is what PROC LIFETEST produced in the SAS release that generated the
#   report. It is not survfit()'s default. Chosen because it reproduces the
#   published limits exactly (27-48 and 24-46 days); the log-log transformation
#   does not (it gives 23-46 for the high-dose arm). Recorded here so the choice
#   is legible rather than buried in an argument.
#
# `test trtan` — the one-degree-of-freedom rank test for association with dose —
# is deliberately NOT computed. See the display's report; both tests the SAS
# step produces are p < 0.0001, and the published figure annotates the
# stratified log-rank, which is what `strata trtan` yields and what survdiff()
# reproduces.

# `custom.R` is sourced into a bare environment, so the package's own helpers
# are not in scope; this display needs one of them.
`%||%` <- function(x, y) if (is.null(x)) y else x

#' Kaplan-Meier estimates and the stratified log-rank test
#'
#' @param data   ADTTE, already restricted to the analysis set by build_ard().
#' @param spec   The analysis entry; `spec$group` names the stratification
#'   variable, `spec$param` the PARAMCD to analyse, and `spec$risk_times` the
#'   grid at which subjects at risk are counted for the strip under the curve.
#' @param denom  The subject-level denominator dataset (unused: ADTTE is already
#'   one record per subject per parameter, so the risk set is the denominator).
#' @return A cards-shaped data frame.
km_time_to_event <- function(data, spec, denom) {
  if (!requireNamespace("survival", quietly = TRUE)) {
    stop("The `survival` package is required by km_time_to_event().", call. = FALSE)
  }
  group <- spec$group
  if (length(group) != 1) {
    stop("km_time_to_event() needs exactly one grouping variable.", call. = FALSE)
  }
  param <- spec$param %||% "TTDE"
  if (!is.null(data$PARAMCD)) {
    data <- data[!is.na(data$PARAMCD) & data$PARAMCD == param, , drop = FALSE]
  }
  if (!nrow(data)) {
    stop("km_time_to_event(): no records for PARAMCD '", param, "'.", call. = FALSE)
  }
  dup <- anyDuplicated(data$USUBJID)
  if (dup) {
    stop(
      "km_time_to_event(): ADTTE has more than one '", param, "' record for ",
      data$USUBJID[dup], "; the risk set would count a subject twice.",
      call. = FALSE
    )
  }

  # `time aval*cnsr(1)`: censored when CNSR is the censoring value, event
  # otherwise. Read from the spec so the mapping is declared, not assumed.
  censor_value <- as.numeric(spec$censor_value %||% 1)
  time <- as.numeric(data[[spec$time_var %||% "AVAL"]])
  cnsr <- as.numeric(data[[spec$censor_var %||% "CNSR"]])
  if (anyNA(time) || anyNA(cnsr)) {
    stop("km_time_to_event(): AVAL or CNSR is missing for at least one subject.", call. = FALSE)
  }
  event <- as.integer(cnsr != censor_value)

  risk_times <- as.numeric(spec$risk_times %||% numeric(0))

  strata <- data[[group]]
  levels <- if (is.factor(strata)) levels(droplevels(strata)) else sort(unique(as.character(strata)))
  strata <- factor(as.character(strata), levels = levels)

  df <- data.frame(.time = time, .event = event, .strata = strata, stringsAsFactors = FALSE)
  surv <- survival::Surv(df$.time, df$.event)

  # conf.type = "plain": see the header note on reproducing the published limits.
  fit <- survival::survfit(surv ~ .strata, data = df, conf.type = "plain")
  tbl <- summary(fit)$table
  if (is.null(dim(tbl))) tbl <- t(as.matrix(tbl))
  rownames(tbl) <- levels

  pieces <- list()
  for (lv in levels) {
    keep <- df$.strata == lv
    n_sub <- sum(keep)
    n_evt <- sum(df$.event[keep])
    n_cns <- n_sub - n_evt
    part <- if (length(levels) == 1) fit else fit[match(lv, levels)]

    # The step function, as coordinates: (0, 1) then one point per distinct
    # time survfit reports, which includes censoring-only times so the curve
    # runs to the end of follow-up rather than stopping at the last event.
    step_time <- c(0, as.numeric(part$time))
    step_surv <- c(1, as.numeric(part$surv))
    cens <- part$n.censor > 0
    row <- function(stat_name, stat_label, value, context) {
      out <- data.frame(
        group1 = group, group1_level = lv,
        group2 = NA_character_, group2_level = NA_character_,
        variable = spec$time_var %||% "AVAL", variable_level = param,
        context = context, stat_name = stat_name, stat_label = stat_label,
        warning = NA_character_, error = NA_character_,
        stringsAsFactors = FALSE
      )
      out$stat <- list(value)
      out
    }
    # The numbers-at-risk strip prints committed statistics, not a second
    # calculation: the count at each declared time is the size of the risk set,
    # the subjects whose observed time has not yet passed. It travels in the ARD
    # so a reader can check the strip against `ard.json` rather than against the
    # renderer's arithmetic.
    risk_n <- vapply(risk_times, function(t) sum(df$.time[keep] >= t), numeric(1))

    pieces <- c(pieces, list(
      row("time", "Time (days)", step_time, "survival"),
      row("surv", "Survival probability", step_surv, "survival"),
      row("n_risk", "Number at risk", c(n_sub, as.numeric(part$n.risk)), "survival"),
      row("n_event_at", "Events at time", c(0, as.numeric(part$n.event)), "survival"),
      row("censor_time", "Censoring time", as.numeric(part$time[cens]), "survival"),
      row("censor_surv", "Survival at censoring", as.numeric(part$surv[cens]), "survival"),
      row("risk_time", "Times at which subjects at risk are counted (days)", risk_times, "survival"),
      row("risk_n", "Subjects at risk", risk_n, "survival"),
      row("N", "Number of subjects", n_sub, "survival_summary"),
      row("n", "Subjects with an event", n_evt, "survival_summary"),
      row("p", "Proportion with an event", n_evt / n_sub, "survival_summary"),
      row("n_censor", "Subjects censored", n_cns, "survival_summary"),
      row("p_censor", "Proportion censored", n_cns / n_sub, "survival_summary"),
      row("median", "Median time to event", unname(tbl[lv, "median"]), "survival_summary"),
      row("median_lcl", "Median 95% LCL", unname(tbl[lv, "0.95LCL"]), "survival_summary"),
      row("median_ucl", "Median 95% UCL", unname(tbl[lv, "0.95UCL"]), "survival_summary")
    ))
  }

  # `strata trtan` in PROC LIFETEST: the test of equality over strata. survdiff()
  # with rho = 0 is the log-rank form, on (number of strata - 1) df.
  if (length(levels) > 1) {
    sd <- survival::survdiff(surv ~ .strata, data = df, rho = 0)
    dfree <- length(sd$n) - 1
    stat_row <- function(stat_name, stat_label, value) {
      out <- data.frame(
        group1 = NA_character_, group1_level = NA_character_,
        group2 = NA_character_, group2_level = NA_character_,
        variable = spec$time_var %||% "AVAL", variable_level = param,
        context = "survival_test", stat_name = stat_name, stat_label = stat_label,
        warning = NA_character_, error = NA_character_,
        stringsAsFactors = FALSE
      )
      out$stat <- list(value)
      out
    }
    pieces <- c(pieces, list(
      stat_row("chisq", "Log-rank chi-square", unname(sd$chisq)),
      stat_row("df", "Degrees of freedom", dfree),
      # Named `p_value`, never `p`: the renderer scales any statistic called `p`
      # to a percentage, and a p-value multiplied by 100 would publish silently.
      stat_row("p_value", "Log-rank p-value", stats::pchisq(sd$chisq, dfree, lower.tail = FALSE))
    ))
  }

  out <- do.call(rbind, pieces)
  rownames(out) <- NULL
  out
}