R/experiments.R

In mellux-project/melluxdrc: melluxdrc

#' Generates samples of measured melatonin values at two lux levels
#'
#' The samples are generated by sampling (without replacement) individuals from a
#' supplied tibble representing a virtual experiment. These samples can either be
#' from different individuals (for a between-type study) or from the same individuals
#' measured twice (for a within-type study).
#'
#' @param is_between a Boolean indicating whether experiment is within or between type
#' @param lux_1 first measured lux value
#' @param lux_2 second measured lux value
#' @param n number of individuals in sample
#' @param population_df a tibble representing a virtual experiment
#'
#' @return a list of two sets of measured melatonin values
generate_two_samples <- function(is_between, lux_1, lux_2, n, population_df) {
  if(!lux_1 %in% population_df$lux)
    stop("lux_1 must be in set of luxes in population_df")
  if(!lux_2 %in% population_df$lux)
    stop("lux_2 must be in set of luxes in population_df")

  n_id <- dplyr::n_distinct(population_df$id)
  if(2 * n > n_id)
    stop("number of individuals in study must be less than half population size.")

  if(is_between) {
    sample_ids <- sample(1:n_id, 2 * n)
    sample_data <- population_df %>%
      dplyr::filter(.data$id %in% sample_ids)
    df_1 <- sample_data %>%
      dplyr::filter(.data$lux == lux_1)
    df_2 <- sample_data %>%
      dplyr::filter(.data$lux == lux_2)

    # pick 1st / 2nd half of ids as population 1 / 2
    vals_1 <- df_1$y[1:n]
    vals_2 <- df_2$y[(n + 1):(2 * n)]
  } else{ # within
    sample_ids <- sample(1:n_id, n)
    sample_data <- population_df %>%
      dplyr::filter(.data$id %in% sample_ids)
    df_1 <- sample_data %>%
      dplyr::filter(.data$lux == lux_1)
    df_2 <- sample_data %>%
      dplyr::filter(.data$lux == lux_2)

    # pick 1st / 2nd half of ids as population 1 / 2
    vals_1 <- df_1$y
    vals_2 <- df_2$y
  }
  list(vals_1=vals_1, vals_2=vals_2)
}

#' Determines if a t-test has found a significant difference of the correct sign
#'
#' @inheritParams generate_two_samples
#' @param vals_1 set of melatonin values at first lux
#' @param vals_2 set of melatonin values at second lux
#' @param fit a result of running t.test
#'
#' @return a named list comprising 'result': a binary value indicating test success (if=1) or failure (if=0) to detect difference of correct sign;
#' and 'p_value': the p-value from the t test
is_comparison_successful <- function(vals_1, vals_2, lux_1, lux_2, fit) {
  res <- 0
  if(lux_2 > lux_1) {
    if(mean(vals_2) > mean(vals_1))
        res <- 1
  } else if (lux_2 < lux_1){
    if(mean(vals_2) < mean(vals_1))
        res <- 1
  } else { # correct detection means no difference detected
      res <- 1
  }
  list(result=res, p_value=fit$p.value)
}

#' Performs t test between two sets of values
#'
#' Ensures that tests are conducted assuming paired values if required
#'
#' @inheritParams generate_two_samples
#' @param vals a named list of 'vals_1' and 'vals_2', each of the same length
#'
#' @return a list with class "htest" resultant from performing a t.test
t_test <- function(is_between, vals) {
  vals_1 <- vals$vals_1
  vals_2 <- vals$vals_2
  if(is_between)
    is_paired <- FALSE
  else
    is_paired <- TRUE

  fit <- stats::t.test(vals_1, vals_2, paired=is_paired)
  fit
}

#' Performs a single between- or within-individual experiments comparing melatonin suppression at two
#' lux levels using t tests
#'
#' For a between-individual experiment, a t-test is used to compare the melatonin suppression level at two lux values for two different
#' samples of individuals. This function requires based an (ideally large) simulated population of individual dose-response data.
#'
#' For a within-individual, a paired t-test is used to compare the melatonin suppression level at two lux values for a sample of individuals
#' measured at each of those values. This function requires based an (ideally large) simulated population of individual dose-response data.
#'
#' @inheritParams generate_two_samples
#' @inheritParams is_comparison_successful
#'
#' @return a named list comprising 'result': a binary value indicating test success (if=1) or failure (if=0) to detect difference of correct sign;
#' and 'p_value': the p-value from the t test
comparison_test_single <- function(is_between, lux_1, lux_2, n, population_df) {

  vals <- generate_two_samples(is_between, lux_1, lux_2, n, population_df)

  fit <- t_test(is_between, vals)

  is_success_df <- is_comparison_successful(vals$vals_1, vals$vals_2, lux_1, lux_2, fit)
  is_success_df
}

#' Performs between- or within-individual experiments comparing melatonin suppression at two
#' lux levels using t tests
#'
#' For a between-individual experiment, a t-test is used to compare the melatonin suppression level at two lux values for two different
#' samples of individuals. This function requires based an (ideally large) simulated population of individual dose-response data.
#'
#' For a within-individual, a paired t-test is used to compare the melatonin suppression level at two lux values for a sample of individuals
#' measured at each of those values. This function requires based an (ideally large) simulated population of individual dose-response data.
#'
#' @inheritParams generate_two_samples
#' @inheritParams is_comparison_successful
#' @param nreps integer indicating number of test replicates (which defaults to 1)
#'
#' @return a dataframe where each row corresponds to a test result in a replicate; note that 'result' indicates that the difference
#' was of the correct sign
#' @export
#' @importFrom rlang .data
#'
#' @examples
#' library(melluxdrc)
#'
#' # generate virtual data for 200 individuals
#' population_df <- virtual_experiment(200)
#'
#' # carry out between-individual experiment comparing lux_1=10 with lux_2=30
#' # for a test sample size of 20
#' is_between <- TRUE
#' comparison_test(is_between, 10, 30, 20, population_df)
#'
#' # carry out within-individual experiments comparing lux_1=10 with lux_2=30
#' # for a test sample size of 20
#' is_between <- FALSE
#' comparison_test(is_between, 10, 30, 20, population_df)
comparison_test <- function(is_between, lux_1, lux_2, n, population_df, nreps=1) {

  m_results <- matrix(nrow = nreps, ncol = 3)
  for(i in 1:nreps) {
    result <- comparison_test_single(is_between, lux_1, lux_2, n, population_df)
    m_results[i, ] <- c(i, result$result, result$p_value)
  }
  colnames(m_results) <- c("replicate", "result", "p_value")
  m_results <- as.data.frame(m_results)
  m_results
}

#' Create two samples of melatonin suppression observations from an intervention type experiment
#'
#' @inheritParams comparison_test
#' @param lux the lux value at which to conduct the comparison
#'
#' @return a list of two sets of measured melatonin values
#' @importFrom rlang .data
generate_two_samples_at_one_lux <- function(is_between, lux, n, population_df) {
  if(!lux %in% population_df$lux)
    stop("lux must be in set of luxes in population_df")
  lux_1 <- lux
  single_lux_df <- population_df %>%
    dplyr::filter(.data$lux == lux_1)

  n_id <- dplyr::n_distinct(single_lux_df$id)
  if(2 * n > n_id)
    stop("number of individuals in study must be less than half population size")

  treatment_types <- single_lux_df %>%
    dplyr::summarise(untreated=sum(!.data$treated),
                     treated=sum(.data$treated))

  enough_untreated <- treatment_types$untreated >  2 * n
  enough_treated <- treatment_types$treated >  2 * n
  if(!(enough_untreated) | !(enough_treated))
    stop("insufficient individuals for comparison at that lux (must exceed twice n)")

  treated_ids <- single_lux_df %>%
    dplyr::filter(.data$treated) %>%
    dplyr::pull(.data$id) %>%
    unique() %>%
    sort()
  untreated_ids <- single_lux_df %>%
    dplyr::filter(!.data$treated) %>%
    dplyr::pull(.data$id) %>%
    unique() %>%
    sort()
  treated_ids <- sample(treated_ids, n)
  untreated_ids <- sample(untreated_ids, n)
  if(is_between) {

    # non-overlapping individuals picked
    df_untreated <- single_lux_df %>%
      dplyr::filter(.data$id %in% untreated_ids) %>%
      dplyr::filter(!.data$treated)
    df_treated <- single_lux_df %>%
      dplyr::filter(.data$id %in% treated_ids) %>%
      dplyr::filter(.data$treated)

  } else{ # within
    sample_ids <- sample(1:n_id, n)

    # pick individuals and look at treated and untreated obs
    sample_data <- single_lux_df %>%
      dplyr::filter(.data$id %in% sample_ids)

    df_untreated <- sample_data %>%
      dplyr::filter(!.data$treated)
    df_treated <- sample_data %>%
      dplyr::filter(.data$treated)
  }

  vals_u <- df_untreated$y
  vals_t <- df_treated$y

  list(untreated=vals_u, treated=vals_t)
}

#' Determines if a t-test has found a significant difference of the correct sign for a single lux
#' treated vs untreated experiment
#'
#' @param vals_untreated responses for the untreated individuals
#' @param vals_treated responses for the treated individuals
#' @param is_treated_higher Boolean indicating whether treated group has a higher response
#' @inheritParams is_comparison_successful
#'
#' @return a named list comprising 'result': a binary value indicating test success (if=1) or failure (if=0) to detect difference of correct sign;
#' and 'p_value': the p-value from the t test
is_comparison_successful_one_lux <- function(vals_untreated, vals_treated, is_treated_higher,
                                             fit) {
  res <- 0
  if(is_treated_higher) {
    if(mean(vals_treated) > mean(vals_untreated))
        res <- 1
  } else {
    if(mean(vals_treated) < mean(vals_untreated))
        res <- 1
  }
  list(result=res, p_value=fit$p.value)
}

#' Performs between- or within-individual experiments comparing melatonin suppression for
#' treated and untreated subgroups at a single lux
#'
#' For a between-individual comparison, a t-test comparing the melatonin suppression level
#' at a single lux is conducted on two subgroups (untreated and treated) comprising different
#' sets of individuals.
#'
#' For a within-individual comparison, a paired t-test comparing the melatonin suppression level
#' at a single lux level for a set of individuals before and after they are treated is conducted.
#'
#' This function requires based an (ideally large) simulated population of individual dose-response data.
#'
#' @inheritParams comparison_test
#' @param population_treated_df a virtual population data frame with individual who have been treated
#' and those who have been untreated (where treatment is indicated by a Boolean column named "treated")
#' @inheritParams generate_two_samples_at_one_lux
#' @inheritParams is_comparison_successful_one_lux
#'
#' @return a binary value indicating test success (if=1) or failure (if=0) to detect difference of correct sign
comparison_test_treatment_single <- function(is_between, lux, n, population_treated_df, is_treated_higher) {

  vals <- generate_two_samples_at_one_lux(is_between, lux, n, population_treated_df)

  # todo: fix this hack
  vals_cop <- list()
  vals_cop$vals_1 <- vals$untreated
  vals_cop$vals_2 <- vals$treated
  fit <- t_test(is_between, vals_cop)

  is_success <- is_comparison_successful_one_lux(
    vals$untreated, vals$treated, is_treated_higher,
    fit)

  is_success
}

#' Performs between- or within-individual experiments comparing melatonin suppression for
#' treated and untreated subgroups at a single lux
#'
#' For a between-individual comparison, a t-test comparing the melatonin suppression level
#' at a single lux is conducted on two subgroups (untreated and treated) comprising different
#' sets of individuals.
#'
#' For a within-individual comparison, a paired t-test comparing the melatonin suppression level
#' at a single lux level for a set of individuals before and after they are treated is conducted.
#'
#' This function requires based an (ideally large) simulated population of individual dose-response data.
#'
#' @inheritParams comparison_test_treatment_single
#' @inheritParams comparison_test
#'
#' @return a dataframe where each row corresponds to a test result in a replicate; note that 'result' indicates that the difference
#' was of the correct sign
#' @export
#'
#' @examples
#' library(melluxdrc)
#' library(purrr)
#' library(dplyr)
#'
#' # create a population where half are treated; half not
#' n <- 200
#' population_df <- virtual_treatment_experiment(n, treated_ed50_multiplier=0.5)
#'
#' # carry out a within-study test at 10 lux for 30 people
#' # note that, because the ed50 multiplier for the treated group was
#' # <1, this group should have a higher response
#' is_treated_higher <- TRUE
#' comparison_test_treatment(FALSE, 10, 30, population_df, is_treated_higher)
comparison_test_treatment <- function(
  is_between, lux, n, population_treated_df, is_treated_higher,
  nreps=1) {
  m_results <- matrix(nrow = nreps, ncol = 3)
  for(i in 1:nreps) {
    result <- comparison_test_treatment_single(is_between, lux, n, population_treated_df, is_treated_higher)
    m_results[i, ] <- c(i, result$result, result$p_value)
  }
  colnames(m_results) <- c("replicate", "result", "p_value")
  m_results <- as.data.frame(m_results)
  m_results
}