R/modify_standata.R

In bentrueman/bgamcar1: Fit Bayesian GAMs With CAR(1) Errors To Censored Data

#' Modify Stan data to include left-censoring limits on predictor variables.
#'
#' @param sdata A list generated by `brms::make_standata()`
#' @param data A dataframe containing the model variables.
#' @param lcl A numeric vector of censoring limits for left-censored variables. Alternatively, a list of numeric vectors with censoring limits for
#' each variable and observation; the order should correspond to that of `lcl` and `cens_ind`.
#' @param var_xcens A character vector containing the names of predictor variables with left-censoring; the order should
#' correspond to that of `lcl` and `cens_ind`.
#' @param cens_ind A character vector containing the names of left-censoring indicators variables corresponding
#' to the contents of `var_xcens`.
#'
#' @return A list that can be passed to `rstan::stan()` or `cmdstanr::cmdstan_model`.
#' @export
#'
#' @examples
#' library("brms")
#' data <- read.csv(paste0(system.file("extdata", package = "bgamcar1"), "/data_gam.csv"))
#' data$f_cens <- ifelse(data$f < 5, "left", "none")
#' sdata <- make_standata(bf(y ~ f), data)
#' sdata_modified <- modify_standata(sdata, data, lcl = 5, var_xcens = "f", cens_ind = "f_cens")
#' str(sdata_modified)
modify_standata <- function(sdata, data, lcl, var_xcens, cens_ind) {

  lcl <- as.list(lcl)

  stopifnot(
    "lengths of 'var_xcens' and 'lcl' must be equal" = length(lcl) == length(var_xcens),
    "lengths of 'var_xcens' and 'cens_var' must be equal" = length(var_xcens) == length(cens_ind)
  )

  varstan <- make.names(str_remove_all(var_xcens, "_"))

  for(i in seq_along(var_xcens)) {
    # make logical censoring indicator:
    cens_logical <- data[[cens_ind[i]]] == "left"

    if (sum(cens_logical) > 0) { # append to standata list only if there are left-censored values
      sdata[[paste0("Ncens_", varstan[i])]] <- sum(cens_logical) # number of left-censored
      # positions of left-censored:
      sdata[[paste0("Jcens_", varstan[i])]] <- as.array(seq_len(nrow(data))[cens_logical])
      sdata[[paste0("U_", varstan[i])]] <- lcl[[i]] # left-censoring limit
    } else message(glue("No left-censored {var} values."))
  }

  sdata
}