R/CI_smd_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_smd_one` returns the point estimate and confidence interval
#' for the standardized mean difference between a sample mean and 
#' a population/reference mean
#' \mjdeqn{ d_{1} = \frac{X - \mu}{s}}{Mdiff = (X - mu)/s}
#' 
#' 
#' @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 standard deviation
#'   and a t-distribution.  If a numeric value is passed, returns a CI 
#'   based on population_s and a normal-distribution
#' @param conf_level The confidence level in
#'   decimal form.  Defaults to 0.95.
#' 
#' 
#' @return Returns a list with these named elements: 
#' * effect_size - the point estimate; d = (comparison_m - population_m)/s
#' * lower - lower bound of the CI 
#' * upper - upper bound of the CI 
#' * numerator - the numerator for Cohen's d_biased;  
#'    (comparison_m - population_m)
#' * denominator - the denominator for Cohen's d_biased; 
#'     population_s if passed, otherwise comparison_s
#' * df - the degrees of freedom used for correction and CI calculation
#' * se - the standard error of the estimate; **warning** not totally 
#'     sure about this yet 
#' * moe - margin of error; 1/2 length of the CI 
#' * d_biased - Cohen's d without correction applied 
#' * properties - a list of properties for the result
#'
#' Properties
#' * d_name - if equal variance assumed d_s, otherwise d_avg
#' * d_name_html - html representation of d_name
#' * denominator_name - if equal variance assumed sd_pooled otherwise sd_avg
#' * denominator_name_html - html representation of denominator name
#' * bias_corrected - TRUE/FALSE if bias correction was applied
#' * message - a message explaining denominator and correction status
#' * message_html - html representation of message
#' 
#' 
#' @examples 
#'   CI_smd_one(
#'     comparison_m = 10, 
#'     comparison_s = 3, 
#'     comparison_n = 25
#'   )   
#' 
#' @export
CI_smd_one <- function(
  comparison_m,
  comparison_s,
  comparison_n,
  population_m = 0,
  population_s = NULL,
  conf_level = 0.95,
  correct_bias = TRUE
) {

  # Options --------------------------------------------
  known_s <- !is.null(population_s)
  do_correct <- correct_bias & !known_s
  
  
  # Initial calculations -------------------------------  
  # Degrees of freedom
  df <- comparison_n - 1
  
  # Correction factor
  J <- if(do_correct)
    exp ( lgamma(df/2) - (log(sqrt(df/2)) + lgamma((df-1)/2) ) )
  else
    1
  
  
  # Calculate d -----------------------------------------------------------
  numerator <- comparison_m - population_m
  
  denominator <- if(known_s)
    population_s
  else
    comparison_s
  
  d_biased <- numerator / denominator
  
  effect_size <- d_biased * J
  
  
  # Calculate CI --------------------------------------------------------    
  lambda <- effect_size * sqrt(comparison_n)
  back_from_lambda <- 1 / sqrt(comparison_n)
  
  lambda_low <- sadists::qlambdap(1/2-conf_level/2, df = df, t = lambda )
  lambda_high <- sadists::qlambdap(1/2+conf_level/2, df = df, t = lambda )

  k <- sqrt(2/df)* exp((lgamma((df+1)/2)) - (lgamma((df)/2)))
  lambda_se <- sqrt(1 + lambda^2 * (1-k^2))
  se <- lambda_se * back_from_lambda
  
  lower <- lambda_low * back_from_lambda 
  upper <- lambda_high * back_from_lambda 
  moe <- (upper-lower)/2
  
  
  properties <- list(
    d_name = "d_1",
    d_name_html = if (do_correct | known_s) 
      "<i>d</i><sub>1</sub>" 
    else 
      "<i>d</i><sub>1.biased</sub>",
    denominator_name = if(known_s) 
      "s_population"
    else
      "s_sample",
    denominator_name_html = if(known_s)
      "<i>s</i><sub>population</sub>"
    else
      "<i>s</i><sub>sample</sub>",
    bias_corrected = do_correct
  )
  
  properties$message <- glue::glue(
    "
This standardized mean difference is called {properties$d_name} because 
the standardizer used was {properties$denominator_name}.
{properties$d_name} {if (properties$bias_corrected) 'has' else 'has *not*'} been 
corrected for bias.  
Correction for bias can be important when df < 50.
    "
  )
  
  properties$message_html <- glue::glue(
    "
This standardized mean difference is called {properties$d_name_html} 
because the standardizer used was {properties$denominator_name_html}.<br>
{properties$d_name_html} {if (properties$bias_corrected) 'has' else 'has *not*'} 
been corrected for bias.  
Correction for bias can be important when <i>df</i> < 50.<br>
    "
  )
  
  # Put in table form
  res <- list(
    effect_size = effect_size,
    lower = lower,
    upper = upper,
    numerator = numerator,
    denominator = denominator,
    df = df,
    se = se,
    moe = moe,
    d_biased = d_biased,
    properties = properties
  )
  
  
  return(res)
}