calculate the change over a time variable, relative to one of the time points

relChange.R

#' Calculate relative change across time
#'
#' @param value a vector with the values to calculate the change over.
#' @param time a vector with the time points.
#' @param ref_time the reference time point to scale the change by.
relChange <- function(value, time, ref_time = min(time)){

  # check if reference time is contained in the given vector
  if(any(is.na(value))) stop("Missing values in 'value' not supported. Remove them first.")
  if(any(is.na(time))) stop("Missing values in 'time' not supported. Remove them first.")
  if(length(ref_time) != 1) stop("Provide a single value for 'ref_time'.")
  if(!(ref_time %in% time)) stop("'ref_time' is not present in time vector.")
  if(length(value) != length(time)) stop("'value' and 'time' should be of same length.")

  # # check if time is in order:
  # if(any(diff(time) < 0)) stop("Please make sure time is in ascending order")

  # Calculate the average per time point
  mean_by_time <- tapply(value, time, mean, simplify = TRUE)

  # get the value for the reference time point
  ref_value <- mean_by_time[names(mean_by_time) == ref_time]

  # calculate the scaled relative change
  #rel_change <- diff(mean_by_time)/ref_value
  rel_change <- (mean_by_time-ref_value)/ref_value

  # return as a data.frame
  data.frame(time = as.numeric(names(rel_change)), change = rel_change, stringsAsFactors = FALSE)
}

#' Calculate relative change with bootstrapped intervals
#'
#' @param value a vector with the values to calculate the change over.
#' @param time a vector with the time points.
#' @param ref_time the reference time point to scale the change by.
#' @param nboot number of bootstrap samples.
#' @param summary_quantiles the quantiles to summarise the bootstrap values.
relChangeBoot <- function(value, time, ref_time = min(time),
                          nboot = 1000,
                          summary_quantiles = c(0.005, 0.025, 0.975, 0.995)){
  # make sure value and time are sorted
  value <- value[order(time)]
  time <- time[order(time)]

  # resampling function
  stratified_sampling <- function(value, time){
    boot_sample <- tapply(value, time, sample, replace = TRUE, simplify = FALSE)
    boot_sample <- unlist(boot_sample)
    names(boot_sample) <- NULL
    return(boot_sample)
  }

  # resample the data
  boot_samples <- replicate(nboot, stratified_sampling(value, time))

  # calculate relative change for each of them
  boot_changes <- apply(boot_samples, 2, relChange, time = time, ref_time = ref_time)
  boot_changes <- do.call("rbind", boot_changes)

  # summarise bootstraps
  boot_summary <- tapply(boot_changes$change, boot_changes$time, quantile,
                         probs = summary_quantiles)

  # tidy into a data.frame
  boot_df <- as.data.frame(do.call("rbind", boot_summary))
  names(boot_df) <- paste0("q", summary_quantiles)

  # return result with mean change as well
  cbind(relChange(value, time, ref_time), boot_df)

}