Check possible p, t, and SD of difference score for within-subject means and SDs over a range of r-values

within_check.R

library(dplyr)
library(faux)
library(tidyr)
library(purrr)
library(ggplot2)
library(glue)

within_t <- function(n = 100, mu1 = 0, mu2 = 0, sd1 = 1, sd2 = 1, r = 0, 
         alternative = c("two.sided", "less", "greater")) {
  x <- faux::rnorm_multi(n = n,
                         mu = c(A = mu1, B = mu2), 
                         sd = c(sd1, sd2),
                         r = r, 
                         empirical = TRUE)
  
  test <- t.test(x$A, x$B, 
                 paired = TRUE, 
                 alternative = match.arg(alternative))
  
  # test info
  list(
    diff_mean = mu2 - mu1,
    diff_sd = sd(x$A - x$B),
    t = test$statistic[[1]],
    p = test$p.value
  )
}


check_p <- function(n = 100, mu1 = 0, mu2 = 0, sd1 = 1, sd2 = 1, 
         r_min = -.99, r_max = .99, r_digits = 2,
         alternative = c("two.sided", "less"),
         reported_p = NULL, plot_type = c("p", "t", "diff_sd"),
         save = NULL) {
  plot_type <- match.arg(plot_type)
  # make sure r-values are in real range
  r_min <- max(-.99, r_min)
  r_max <- min(.99, r_max)

  # set up params
  params <- tidyr::crossing(
    n = n,
    mu1 = mu1,
    mu2 = mu2,
    sd1 = sd1,
    sd2 = sd2,
    r = seq(r_min, r_max, 10^(-r_digits)),
    alternative = alternative
  ) %>%
    # calculate p and t-values
    dplyr::bind_cols(purrr::pmap_df(., within_t))
  
  # get plausible p-values if reported_p is set
  if (!is.null(reported_p)) {
    params <- params %>%
      dplyr::mutate(plausible = 
             abs(round(p, r_digits) - reported_p) < 
             10^(-r_digits))
    
    r_range <- r_max - r_min
    p_x <- params %>%
      filter(plausible) %>%
      summarise(max = max(r),
                min = min(r)) %>%
      mutate(x = ifelse(max > .5, r_min + r_range*.1, r_max - r_range*.1)) %>%
      pull(x)
  }
  
  multi <- list(mu1 = mu1, mu2 = mu2, sd1 = sd1, sd2 = sd2, n = n,
                alternative = alternative) %>%
    lapply(length)
  grp <- names(multi[multi>1])
  if (length(grp) == 0) grp <- "all"
  
  # plot the results
  plot_n <- paste(n, collapse = ',')
  plot_m1 <- round(mu1, 2) %>% paste(collapse = ',')
  plot_sd1 <- round(sd1, 2) %>% paste(collapse = ',')
  plot_m2 <- round(mu2, 2) %>% paste(collapse = ',')
  plot_sd2 <- round(sd2, 2) %>% paste(collapse = ',')
  
  plot <- params %>%
    tidyr::unite(groups, grp, sep = "_", remove = FALSE) %>%
    ggplot() +
    labs(
      title = glue("N = {plot_n}; M1 = {plot_m1} ({plot_sd1}); M2 = {plot_m2} ({plot_sd2})"),
      color = paste(grp, collapse = ", "),
      x = "r-value",
      y = switch(plot_type,
                 p = "p-value", 
                 t = "t-value", 
                 diff_sd = "SD of the diference score")
    )
  
  # add reference line for p-value
  if (plot_type == "p" && !is.null(reported_p)) {
    plot <- plot + 
      geom_hline(yintercept = reported_p, alpha = 0.6) +
      annotate("label", x = p_x, y = reported_p, label = glue("reported\np = {reported_p}"))
  }
  
  plot <- plot + 
    geom_line(aes(x = r, y = .data[[plot_type]], color = groups), 
              size = 1)

  if (is.null(save)) {
    print(plot)
  } else {
    ggsave(filename = save, plot = plot, width = 5, height = 5)
  }
  
  invisible(params)
}

# test

params <- check_p(n = 24, 
        mu1 = 25, mu2 = 32, 
        sd1 = 30.51, sd2 = 30.29, 
        plot_type = "p",
        reported_p = 0.03)