R/compute_thresholds_v2.R

In CJangelo/COA34: Psychometric Analyses Aligned with FDA COA Guidances 3 and 4

#' Compute Thresholds
#'
#' Note that this computes quantiles using the `quantile()` function in R.
#' In order to align with SAS, the type 2 definition is used
#' Refer to https://blogs.sas.com/content/iml/2017/05/24/definitions-sample-quantiles.html
#'
#' @param dat pass the dataframe
#' @param anchor.group variable of the anchor group you want used here
#' @param time.var variable of the Time in the dataframe - PLEASE CHECK THAT THIS IS CORRECTLY ORDERED, default is Time
#' @param timepoint default here is to just use the final timepoint, e.g., "Time_4". Please be sure to have
#' ordered your time variable correctly.
#' @param change.score the PRO change score you want the threshold computed on
#' @return returns a dataframe with two additional variables, PRO score at baseline and PRO change score
#' @export


compute_thresholds_v2 <- function(dat = NULL,
                                  anchor.group = NULL,
                                  time.var = NULL,
                                  timepoint = NULL,
                                  change.score = NULL,
                                  subject.id = 'USUBJID'
){


  if (is.null(dat)) stop('Please specify dataframe in `compute_thresholds()` ')
  if (is.null(anchor.group)) stop('Please specify anchor.group in `compute_thresholds()` ')
  if (is.null(time.var)) stop('Please specify time.var in `compute_thresholds()` ')
  if (is.null(change.score)) stop('Please specify change.score in `compute_thresholds()` ')



  # Select Final Timepoint:
  if (is.null(timepoint)) {
    final.timepoint <- sort(unique(dat[, time.var, drop = T]), decreasing = T)[1]
  } else {
    final.timepoint <- timepoint
  }


  dat <- dat[which(dat[ , time.var] == final.timepoint), ]
  #N <- length(unique(dat$USUBJID))
  N <- length(unique(dat[ , subject.id, drop = T]))

  # Create table:
  tf <- as.formula(paste0(change.score, '~', anchor.group))
  out <- aggregate(tf,
                   function(x) c('mean' = mean(x, na.rm = T),
                                 'median' = quantile(x, probs = c(0.25, 0.5, 0.75), type = 2, na.rm = T),
                                 'n' = sum(!is.na(x)),
                                 'percent' = 100*sum(!is.na(x))/N),
                   data = dat,
                   na.action = na.pass)
  out <- do.call(data.frame, out)
  out <- cbind(change.score, out)
  colnames(out) <- c('PRO Change Score',
                     'Anchor Group',
                     'Mean',
                     '25th percentile',
                     '50th percentile - Median',
                     '75th percentile',
                     'N',
                     'Percent')

  return(out)

}