R/as-ped.R

In pammtools: Piece-Wise Exponential Additive Mixed Modeling Tools for Survival Analysis

#' Transform data to piece-wise exponential data (PED)
#'
#' This is the general data transformation function provided by the
#' \code{pammtools} package. The following main applications must be distinguished:
#' \enumerate{
#'  \item Transformation of standard time-to-event data.
#'  \item Transformation of left-truncated time-to-event data.
#'  \item Transformation of time-to-event data with time-dependent covariates (TDC).
#'  \item Transformation of competing risks data (single or stacked data sets).
#'  \item Transformation of recurrent events and multi-state data.
#' }
#' For TDC data, the type of effect one wants to estimate is also
#' important for the data transformation step. In case of TDCs, the
#' right-hand-side of the formula can contain formula specials
#' \code{concurrent} and \code{cumulative}.
#'
#' For competing risks data, \code{as_ped} can return either:
#' \itemize{
#'  \item A list of cause-specific data sets (\code{combine = FALSE}), where
#'  each element corresponds to one event type and uses cause-specific
#'  interval split points. This is suitable for cause-specific hazards
#'  models without shared effects.
#'  \item A single stacked data set (\code{combine = TRUE}, the default),
#'  where all cause-specific data sets are combined with a \code{cause}
#'  column as covariate. Common split points are derived from all event
#'  times. This is required for models with shared covariate effects across
#'  causes, estimated via interaction terms (e.g.,
#'  \code{s(tend, by = cause)}).
#' }
#'
#' For multi-state data, \code{as_ped} extends the standard PED transformation
#' to each transition type. The follow-up of each subject is split at all
#' observed transition times across the entire dataset, and a row is added for
#' every interval-transition combination the subject is at risk for. Two key
#' differences arise compared to the single-event case:
#' \itemize{
#'  \item Delayed entry into the risk set is handled automatically, since
#'  subjects are only at risk for transitions out of a state after they have
#'  entered it.
#'  \item Competing events are treated as censoring for all other transitions
#'  within the same interval.
#' }
#'
#' In any case, the data transformation is specified by a two-sided formula.
#' See the
#' \href{https://adibender.github.io/pammtools/articles/data-transformation.html}{data-transformation},
#' \href{https://adibender.github.io/pammtools/articles/competing-risks.html}{competing-risks},
#' and
#' \href{https://adibender.github.io/pammtools/articles/recurrent-events.html}{recurrent-events}
#' vignettes for details.
#'
#' @rdname as_ped
#' @param data Either an object inheriting from data frame or in case of
#' time-dependent covariates a list of data frames (of length 2), where the
#' first data frame contains the time-to-event information and static covariates
#' while the second (and potentially further data frames) contain information on
#' time-dependent covariates and the times at which they have been observed.
#' @param formula A two sided formula with a \code{\link[survival]{Surv}} object
#' on the left-hand-side and covariate specification on the right-hand-side (RHS).
#' The RHS can be an extended formula, which specifies how TDCs should be
#' transformed using specials \code{concurrent} and \code{cumulative}. The
#' left-hand-side can be in start-stop notation. This, however, is only used
#' to create left-truncated data and does not support the full functionality.
#' @param cut Split points, used to partition the follow-up into intervals.
#' If unspecified, all unique event times will be used. For competing risks,
#' when \code{combine = TRUE} split points are derived from all event types
#' combined.
#' @param max_time If \code{cut} is unspecified, this will be the last
#' possible event time. All event times after \code{max_time} will be
#' administratively censored at \code{max_time}.
#' @param tdc_specials A character vector of names of potential specials in
#' \code{formula} for concurrent and/or cumulative effects.
#' @param censor_code Specifies the value of the status variable that indicates
#' censoring. Often this will be \code{0}, which is the default.
#' @param transition Character string. Name of the column in \code{data} that
#' identifies the transition type in multi-state models. When supplied,
#' \code{as_ped} performs the multi-state PED transformation, stacking
#' interval-transition rows for each subject.
#' @param timescale Character string, either \code{"gap"} (time since
#' last transition) or \code{"calendar"} (time since study entry, not reset
#' after each transition).
#' @param ... Further arguments passed to the \code{data.frame} method and
#' eventually to \code{\link[survival]{survSplit}}. Notably, \code{id}
#' (character string) sets the name of the subject identifier variable in
#' \code{data}, and, for competing risks data, \code{combine} (logical, default
#' \code{TRUE}) controls whether cause-specific data sets are stacked into a
#' single data frame with an additional \code{cause} column (\code{TRUE}) or
#' returned as a list of cause-specific data sets (\code{FALSE}); see the
#' \href{https://adibender.github.io/pammtools/articles/competing-risks.html}{competing-risks}
#' vignette for details.
#' @importFrom Formula Formula
#' @examples
#' # Standard single-event transformation
#' tumor[1:3, ]
#' tumor[1:3, ] %>% as_ped(Surv(days, status) ~ age + sex, cut = c(0, 500, 1000))
#' tumor[1:3, ] %>% as_ped(Surv(days, status) ~ age + sex)
#'
#' # Competing risks: stacked data set (combine = TRUE, default)
#' # Suitable for cause-specific hazards models with shared effects,
#' # estimated via interaction terms e.g. s(tend, by = cause)
#' \dontrun{
#' data("fourD", package = "etm")
#' ped_stacked <- fourD %>%
#'   as_ped(Surv(time, status) ~ ., id = "id")
#' head(ped_stacked)
#'
#' # Competing risks: list output (combine = FALSE)
#' # Suitable for cause-specific hazards models without shared effects
#' ped_list <- fourD %>%
#'   as_ped(Surv(time, status) ~ ., id = "id", combine = FALSE)
#' # ped_list[[1]]: data for cause 1 (cardiovascular death)
#' # ped_list[[2]]: data for cause 2 (death from other causes)
#' head(ped_list[[1]])
#' head(ped_list[[2]])
#'
#' # Multi-state: illness-death model on calendar timescale
#' # Uses the prothr data (liver cirrhosis patients, n = 488) from mstate.
#' # Patients can transition between normal (1) and abnormal (2) prothrombin
#' # levels and death (3): transitions 1->2, 1->3, 2->1, 2->3.
#' # Calendar timescale is used because hazards depend on overall disease
#' # duration, not time since last transition.
#' data("prothr", package = "mstate")
#' ped_msm <- prothr %>%
#'   filter(Tstart != Tstop) %>%
#'   as_ped(
#'     formula    = Surv(Tstart, Tstop, status) ~ .,
#'     transition = "trans",
#'     id         = "id",
#'     timescale  = "calendar",
#' )
#' head(ped_msm)
#' }
#' @return For standard and left-truncated data, a data frame of class
#' \code{ped} in piece-wise exponential data format. For competing risks data,
#' either a stacked data frame of class \code{ped_cr} (when
#' \code{combine = TRUE}) or a list of cause-specific \code{ped} data frames
#' of class \code{ped_cr_list} (when \code{combine = FALSE}). For multistate data,
#' the result is a stacked long-format dataset with one row per subject,
#' interval, and transition, which can be passed directly to a Poisson
#' regression model.
#' @export

as_ped <- function(data, ...) {
  UseMethod("as_ped", data)
}

#' @rdname as_ped
#' @export
as_ped.data.frame <- function(
  data,
  formula,
  cut = NULL,
  max_time = NULL,
  tdc_specials = c("concurrent", "cumulative"),
  censor_code = 0L,
  transition = character(),
  timescale = c("gap", "calendar"),
  min_events = 1L,
  ...
) {
  # interval-censored data are detected before status_error(), which assumes a
  # 0/1 (or competing) status column and would choke on survival's interval
  # status codes. Right-censored and left-truncated (counting-process) responses
  # both return "none" here and fall through to the standard pipeline unchanged.
  if (detect_ic(formula, data) != "none") {
    ped <- as_ped_ic(
      data = data,
      formula = formula,
      cut = cut,
      max_time = max_time,
      ...
    )
    return(ped)
  }

  status_error(data, formula, censor_code)
  assert_subset(tdc_specials, c("concurrent", "cumulative"))

  if (test_character(transition, min.chars = 1L, min.len = 1L)) {
    ped <- as_ped_multistate(
      data = data,
      formula = formula,
      cut = cut,
      max_time = max_time,
      tdc_specials = tdc_specials,
      censor_code = censor_code,
      transition = transition,
      timescale = timescale,
      min_events = min_events,
      ...
    )
    return(ped)
  }

  event_types <- get_event_types(data, formula, censor_code)
  if (length(event_types) > 1) {
    ped <- as_ped_cr(
      data = data,
      formula = formula,
      cut = cut,
      max_time = max_time,
      tdc_specials = tdc_specials,
      censor_code = censor_code,
      ...
    )
  } else {
    dots <- list(...)
    dots$data <- data
    dots$formula <- get_ped_form(
      formula,
      data = data,
      tdc_specials = tdc_specials
    )
    dots$cut <- cut
    dots$max_time <- max_time

    ped <- do.call(split_data, dots)
    attr(ped, "time_var") <- get_lhs_vars(dots$formula)[1]
    attr(ped, "status_var") <- get_lhs_vars(dots$formula)[2]
  }

  ped
}

#' @rdname as_ped
#' @export
as_ped.nested_fdf <- function(
  data,
  formula,
  ...
) {
  dots <- list(...)
  # update interval break points (if necessary)
  cut <- dots$cut
  if (is.null(cut)) {
    cut <- attr(data, "breaks")
  }
  ccr_breaks <- attr(data, "ccr_breaks")
  cut <- union(cut, ccr_breaks[ccr_breaks <= max(cut)]) %>% sort()

  # ped <- data %>%
  #   select_if(is.atomic) %>%
  #   as.data.frame() %>%
  #   as_ped(
  #     formula  = formula,
  #     id       = dots$id,
  #     cut      = cut,
  #     max_time = dots$max_time,
  #     ...)
  dots$formula <- formula
  dots$data <- as.data.frame(select_if(data, is.atomic))
  dots$cut <- cut
  ped <- do.call(as_ped, dots)

  # replace updated attributes
  attr(data, "breaks") <- attr(ped, "breaks")
  attr(data, "id_n") <- ped %>%
    group_by(!!sym(attr(data, "id_var"))) %>%
    summarize(id_n = n()) %>%
    pull("id_n") %>%
    as_vector()
  attr(data, "id_tseq") <- ped %>%
    group_by(!!sym(attr(data, "id_var"))) %>%
    transmute(id_tseq = row_number()) %>%
    pull("id_tseq") %>%
    as_vector()
  attr(data, "id_tz_seq") <- rep(
    seq_len(nrow(data)),
    times = attr(data, "id_n")
  )

  if (has_special(formula, "concurrent")) {
    ped <- ped %>% add_concurrent(data = data, id_var = dots$id)
  }

  if (has_special(formula, "cumulative")) {
    ped <- add_cumulative(ped, data = data, formula = formula)
    attr(ped, "ll_weights") <- imap(
      attr(ped, "tz"),
      ~ bind_cols(!!.y := .x, ll_weight = c(mean(abs(diff(.x))), abs(diff(.x))))
    )
    class(ped) <- c("fped", class(ped))
  }
  attr(ped, "time_var") <- get_lhs_vars(formula)[1]
  attr(ped, "func_mat_names") <- make_mat_names(
    attr(ped, "func"),
    attr(ped, "time_var")
  )

  ped
}

#' @rdname as_ped
#' @export
as_ped.list <- function(
  data,
  formula,
  tdc_specials = c("concurrent", "cumulative"),
  censor_code = 0L,
  ...
) {
  assert_class(data, "list")
  assert_class(formula, "formula")

  status_error(data[[1]], formula, censor_code)

  nl <- length(data)
  # form  <- Formula(formula)
  has_tdc <- has_tdc_form(formula, tdc_specials = tdc_specials)

  if (nl == 1 && !has_tdc) {
    ped <- data[[1]] %>%
      as_ped(formula = formula, tdc_specials = tdc_specials, ...)
  } else {
    if (nl == 2 && !has_tdc) {
      stop(
        "Two data sets provided in 'data' but no specification of
      time-dependent covariate effects in 'formula'"
      )
    } else {
      nested_fdf <- nest_tdc(data, formula, ...)
      ped <- as_ped(nested_fdf, formula, ...)
    }
  }
  lhs_vars <- get_lhs_vars(formula)
  attr(ped, "time_var") <- lhs_vars[1]
  attr(ped, "trafo_args")$formula <- formula

  ped
}

#' @rdname as_ped
#' @param x any R object.
#' @export
is.ped <- function(x) inherits(x, "ped")


#' @rdname as_ped
#' @param newdata A new data set (\code{data.frame}) that contains the same
#' variables that were used to create the PED object (\code{data}).
#' @export
as_ped.ped <- function(data, newdata, ...) {
  if (is.ped(newdata)) {
    stop("newdata already in ped format.")
  }

  trafo_args <- attr(data, "trafo_args")
  trafo_args[["data"]] <- newdata
  do.call(as_ped, trafo_args)
}


#' @rdname as_ped
#' @export
as_ped.pamm <- function(data, newdata, ...) {
  if (is.ped(newdata)) {
    stop("newdata already in ped format.")
  }
  trafo_args <- data[["trafo_args"]]
  trafo_args$data <- newdata
  do.call(split_data, trafo_args)
}

## Competing risks

#' Competing risks trafo
#'
#' @inherit as_ped
#' @param combine Logical. If \code{TRUE} (the default), cause-specific data
#'   sets are stacked into a single data frame with an additional \code{cause}
#'   column, using split points common to all event types. If \code{FALSE}, a
#'   list of cause-specific data sets is returned.
#' @importFrom rlang .env
#'
#' @keywords internal
as_ped_cr <- function(
  data,
  formula,
  cut = NULL,
  max_time = NULL,
  tdc_specials = c("concurrent", "cumulative"),
  censor_code = 0L,
  combine = TRUE,
  ...
) {
  lhs_vars <- get_lhs_vars(formula)
  n_lhs <- length(lhs_vars)
  event_types <- get_event_types(data, formula, censor_code)

  cut <- map2(
    event_types,
    if (is.list(cut)) cut else list(cut),
    function(.event, .cut) {
      get_cut(
        data,
        formula = formula,
        cut = .cut,
        max_time = NULL,
        event = .event
      )
    }
  )
  if (length(cut) > 1 && combine) {
    cut <- list(reduce(cut, union))
  }

  ped <- map2(
    event_types,
    cut,
    function(.event, .cut) {
      ped_i <- data %>%
        mutate(
          !!lhs_vars[n_lhs] := 1L *
            (as.character(.data[[lhs_vars[n_lhs]]]) ==
              as.character(.env[[".event"]]))
        ) %>%
        as_ped(
          formula = formula,
          cut = .cut,
          max_time = max_time,
          tdc_specials = tdc_specials,
          ...
        )
      ped_i$cause <- as.factor(.event)
      ped_i
    }
  )

  if (combine) {
    ped <- do.call(rbind, ped)
    class(ped) <- c("ped_cr_union", "ped_cr", class(ped))
    attr(ped, "intvars") <- c(attr(ped, "intvars"), "cause")
    attr(ped, "breaks") <- if (length(cut) == 1) unlist(cut) else cut
  } else {
    class(ped) <- c("ped_cr_list", "ped_cr", "ped", class(ped))
    names(ped) <- paste0("cause = ", event_types)
    attributes(ped)$trafo_args$id <- attributes(ped[[1]])$trafo_args$id
    attributes(ped)$trafo_args$formula <- formula
  }

  attr(ped, "trafo_args")[["cut"]] <- if (length(cut) == 1) unlist(cut) else cut
  attr(ped, "trafo_args")[["combine"]] <- combine
  attr(ped, "trafo_args")[["censor_code"]] <- censor_code
  attr(ped, "risks") <- event_types

  ped
}

#' Extract event types
#'
#' Given a formula that specifies the status variable of the outcome, this function
#' extracts the different event types (except for censoring, specified by
#' \code{censor_code}).
#'
#' @inheritParams as_ped
#'
#' @keywords internal
get_event_types <- function(data, formula, censor_code) {
  lhs_vars <- get_lhs_vars(formula)
  status_values <- unique(data[[lhs_vars[length(lhs_vars)]]]) %>% sort()
  status_values <- status_values[status_values != censor_code]
  if (is.factor(status_values)) {
    droplevels(status_values[status_values != censor_code])
  } else {
    status_values
  }
}


#' Recurrent events trafo
#'
#' @examples
#' \dontrun{
#' data("cgd", package = "frailtyHL")
#' cgd2 <- cgd %>%
#'  select(id, tstart, tstop, enum, status, age) %>%
#'  filter(enum %in% c(1:2))
#' ped_re <- as_ped_multistate(
#'   formula = Surv(tstart, tstop, status) ~ age + enum,
#'   data = cgd2,
#'  transition = "enum",
#'  timescale = "calendar")
#' }
#' @rdname as_ped
#' @export
#' @keywords internal
as_ped_multistate <- function(
  data,
  formula,
  cut = NULL,
  max_time = NULL,
  tdc_specials = c("concurrent", "cumulative"),
  censor_code = 0L,
  transition = character(),
  timescale = c("gap", "calendar"),
  min_events = 1L,
  ...
) {
  assert_character(
    transition,
    min.chars = 1L,
    min.len = 1L,
    any.missing = FALSE,
    len = 1L
  )
  assert_integer(min_events, lower = 1L, len = 1L)

  status_error(data, formula, censor_code)
  assert_subset(tdc_specials, c("concurrent", "cumulative"))

  rhs_vars <- get_rhs_vars(formula)
  if (!(transition %in% rhs_vars)) {
    formula <- add_to_rhs(formula, transition)
  }

  dots <- list(...)
  dots$data <- data
  dots$formula <- get_ped_form(
    formula,
    data = data,
    tdc_specials = tdc_specials
  )
  dots$cut <- sort(unique(cut))
  dots$max_time <- max_time
  dots$transition <- transition
  dots$min_events <- min_events
  dots$timescale <- timescale

  ped <- do.call(split_data_multistate, dots)
  attr(ped, "time_var") <- get_lhs_vars(dots$formula)[1]

  return(ped)
}