R/adherence.R

In jpasquier/ad_test: What the Package Does (One Line, Title Case)

#' Compute the adherence
#'
#' Compute the adherence which is defined as ...
#'
#' @param time individual (and time dependent) times
#' @param comp individual and time dependent compliances (0/1)
#' @param id ids in the longitudinal dataset
#' @param cens censoring indicator in the longitudinal dataset
#'
#' @return estimated compliance and adherence as a list of two dataframes
#' \itemize{
#'   \item ad: a dataframe containing three estimates of adherence
#'     \itemize{
#'       \item time = all observed time points
#        \item imp  = empirical implementation
#        \item pers = estimatated persistence
#        \item w0   = weights for "non-compliant" observation
#        \item w1   = weights for "compliant" observation
#'       \item adh1 = empirical estimate
#'       \item adh2 = indirect estimate
#'       \item adh3 = weighted estimate
#'     }
#'   \item DatLong: a dataframe containing individual and time dependent
#'     compliance data for adherence
#'     \itemize{
#'       \item id   = patients ids
#'       \item time = observation time
#'       \item comp = observed compliance
#'       \item w    = associated weight (to correct adherence estimations)
#'     }
#' }
#' 
#' @examples
#' data(onco)
#' onco$compliance <- onco$observed >= onco$expected
#' adherence(onco$drel, onco$compliance, onco$id, onco$cens)
#'
#' @importFrom survival Surv survfit
#'
#' @export

adherence <- function(time, comp, id, cens) {

  # Compute the observed implementation with the implementation function
  imp <- implementation(time, comp)

  # Estimate the persistence with the Kaplan-Meier estimator
  ## 1) Generate the survival data
  sdat <- SurvDat(id, time, cens)
  ## 2) Estimate the persistence with the Kaplan-Meier estimator
  m <- survival::survfit(survival::Surv(sdat$surv, sdat$cens) ~ 1)
  pers <- data.frame(time = m$time, pers = m$surv)
  ## 3) Complete persistence data
  pers <- do.call(rbind, lapply(imp$time, function(u) {
    if (any(pers$time <= u)) {
      r <- data.frame(
        time = u,
        pers = pers[pers$time == max(pers$time[pers$time <= u]), "pers"]
      )
    } else {
      r <- data.frame(time = u, pers = 1)
    }
  }))
  ## 4) Remove temporary objects
  rm(m)

  # Missing values and discontinuations
  # In order to estimate the adherence, we 
  #   i)  add missing values;
  #   ii) add "compliance=0" in the data when a discontinuation is observed.
  #
  ## 1) Put compliances in a dataframe with times and ids
  cdat <- data.frame(id = id, time = time, comp = comp)
  ## 2) Add missing compliances (merge cdat with all ids and observed times
  ##    combinations using the option all = TRUE)
  cdat <- merge(cdat, expand.grid(id = unique(id), time = unique(time)),
                by = c("id", "time"), all = TRUE)
  ## 3) Replace missing values of comp by 0 in case of discontinuation
  for (i in sdat[sdat$cens == 1 & sdat$surv <= max(time), "id"]) {
    cdat[cdat$id == i & cdat$time >= sdat[sdat$id == i, "surv"], "comp"] <- 0
  }
  rm(i)

  # Probability to simultaneously be non-/observed and compliant
  # ... will be used to compute the weights
  c0obs <- aggregate(comp ~ time, cdat, function(x) mean(!is.na(x) & x == 0))
  names(c0obs) <- c("time", "c0obs")
  c1obs <- aggregate(comp ~ time, cdat, function(x) mean(!is.na(x) & x == 1))
  names(c1obs) <- c("time", "c1obs")

  # Empirical (biased) adherence
  emp_adh <- aggregate(comp ~ time, cdat, mean, na.rm = TRUE)
  names(emp_adh) <- c("time", "adh1")

  # Put all together
  adh <- merge(imp, pers, by = "time")
  adh <- merge(adh, c0obs, by = "time")
  adh <- merge(adh, c1obs, by = "time")
  adh <- merge(adh, emp_adh, by = "time")
  
  # Compute the weights
  adh$w0 <- with(adh, c1obs / (imp * pers))
  adh$w1 <- with(adh, ifelse(imp * pers < 1, c0obs / (1 - imp * pers), 1))

  # Compute the corrected adherence
  ## Indirect adherence
  adh$adh2 <- with(adh, imp * pers)
  ## Weighted adherence
  adh$adh3 <- with(adh, w1 * adh1 / (w1 * adh1 + w0 * (1 - adh1)))

  # Format adh
  adh <- adh[order(adh$time),
             c("time", "imp", "pers", "w0", "w1", "adh1", "adh2", "adh3")]

  # Add weights to cdat
  cdat <- merge(cdat, adh[c("time", "w0", "w1")], by = "time", all.x = TRUE)
  cdat$w <- with(cdat, (1 - comp) * w0 + comp * w1)
  cdat <- cdat[order(cdat$id, cdat$time), c("id", "time", "comp", "w")]

  # Return the results
  return(list(ad = adh, DatLong = cdat))

}