R/CI_mdiff_one.R

In rcalinjageman/esci2: Effect Sizes and Confidence Intervals

#' Estimate the mean difference from a sample mean to a population mean.
#'
#'
#' @description 
#' \loadmathjax
#' `CI_mdiff_one` returns the point estimate and confidence interval
#' for the mean difference between a sample mean and a population/reference mean
#' \mjdeqn{ M_{diff} = X - \mu}{Mdiff = X - mu}
#' 
#' 
#' @param comparison_m Mean from the sample
#' @param comparison_s Standard deviation from the sample
#' @param comparison_n Sample size
#' @param population_m Optional value for population mean; defaults to 0
#' @param population_s Optional value for population sd; defaults to NULL.  If NULL
#'   function returns a CI based on the sample sd and a t-distribution.  If
#'   a numeric value is passed, returns a CI based on the population_sd and
#'   a z-distribution
#' @param conf_level The confidence level in
#'   decimal form.  Defaults to 0.95.
#' 
#' 
#' @return Returns a list with these named elements: 
#' * effect_size - Calculated as comparison_m - population_m
#' * lower - lower bound of the CI
#' * upper - upper bound of the CI 
#' * df - degrees of freedom; NA if population_s is passed; otherwise 
#'    comparison_n - 1
#' * se - standard error
#' * moe - margin of error
#' * variability_component - population_s if passed, otherwise comparison_s
#' 
#' 
#' @section Details:
#' If population_s is passed:
#' * se = population_s / sqrt(comparison_n)
#' * MoE = se * z_critical
#' 
#' If population_s is NULL:
#' * se = comparison_s / sqrt(comparison_n)
#' * MoE = se * t_critical
#' 
#' 
#' @examples 
#'   CI_mdiff_one(
#'     comparison_m = 10, 
#'     comparison_s = 3, 
#'     comparison_n = 25
#'   )   
#' 
#' @export
CI_mdiff_one <- function(
  comparison_m,
  comparison_s,
  comparison_n,
  population_m = 0,
  population_s = NULL,
  conf_level = 0.95
) {
  
  effect_size <- comparison_m - population_m
  
  if(is.null(population_s)) {
    variability_component <- comparison_s
    df = comparison_n - 1
    multiplier <- stats::qt(1-((1 - conf_level)/2), df)
  } else {
    variability_component <- population_s
    df <- NA
    multiplier <- stats::qnorm(1-((1 - conf_level)/2))
  }
  
  se <- variability_component  / sqrt(comparison_n)
  moe <- se * multiplier
  lower <- effect_size - moe
  upper <- effect_size + moe
  
  res <- list(
    effect_size = effect_size,
    lower = lower,
    upper = upper,
    df = df,
    se = se,
    moe = moe,
    variability_component = variability_component 
  )
  
  return(res)
}


CI_mdiff_one_tests <- function() {
  myd <- c(3,7,11,0,7,0,4,5,6,2)
  
  
  #Basic call
  myres <- CI_mdiff_one(
    comparison_m = mean(myd), 
    comparison_s = sd(myd), 
    comparison_n = length(myd)
  )
  tes <- t.test(x = myd)
  all.equal(tes$conf.int[[1]], myres$lower)
  all.equal(tes$conf.int[[2]], myres$upper)
  all.equal(tes$parameter, myres$df, check.attributes = FALSE)
  all.equal(tes$estimate, myres$effect_size, check.attributes = FALSE)

  # 99% CI
  myres <- CI_mdiff_one(
    comparison_m = mean(myd), 
    comparison_s = sd(myd), 
    comparison_n = length(myd), 
    conf_level = .99
  )
  tes <- t.test(x = myd, conf.level = 0.99)
  all.equal(tes$conf.int[[1]], myres$lower)
  all.equal(tes$conf.int[[2]], myres$upper)
  all.equal(tes$parameter, myres$df, check.attributes = FALSE)
  all.equal(tes$estimate, myres$effect_size, check.attributes = FALSE)
  
  
  # 99% CI and reference mean not 0
  myres <- CI_mdiff_one(
    comparison_m = mean(myd), 
    comparison_s = sd(myd), 
    comparison_n = length(myd), 
    population_m = 2.5,
    conf_level = .99
  )
  tes <- t.test(x = myd - 2.5, conf.level = 0.99)
  all.equal(tes$conf.int[[1]], myres$lower)
  all.equal(tes$conf.int[[2]], myres$upper)
  all.equal(tes$parameter, myres$df, check.attributes = FALSE)
  all.equal(tes$estimate, myres$effect_size, check.attributes = FALSE)
  
  
  # 95% CI and z-test with sigma = 2
  myres <- CI_mdiff_one(
    comparison_m = 10, 
    comparison_s = 3, 
    comparison_n = 25, 
    population_s = 2,
    conf_level = .95
  )  
  all.equal(myres$lower, 9.216, tolerance = 0.01)
  all.equal(myres$upper, 10.784, tolerance = 0.01)
  all.equal(myres$effect_size, 10)
  
  # 99% and z-test
  myres <- CI_mdiff_one(
    comparison_m = 10, 
    comparison_s = 3, 
    comparison_n = 25, 
    population_s = 2,
    conf_level = .99
  )  
  all.equal(myres$lower, 8.97, tolerance = 0.01)
  all.equal(myres$upper, 11.03, tolerance = 0.01)
  all.equal(myres$effect_size, 10)
}