esci2: Effect Sizes and Confidence Intervals

Documented in CI_mdiff_contrast_bs

#' Estimate the mean difference for an independent groups contrast.
#'
#' @description 
#' \loadmathjax
#' `CI_mdiff_contrast_bs` returns the point estimate and confidence interval
#' for the mean difference in a linear contrast:
#' \mjdeqn{ \psi = \sum_{i=1}^{a}c_iM_i }{psi = sum(contrasts*means)}
#' Where there are *a* groups, and *M* is each group mean and *c* is each group
#' weight; see Kline, equation 7.1
#'
#' @param means A vector of 2 or more means
#' @param sds A vector of standard deviations, same length as means
#' @param ns A vector of sample sizes, same length as means
#' @param contrast A vector of group weights, same length as means. 
#' @param conf_level The confidence level for the confidence interval in
#'   decimal form.  Defaults to 0.95.
#' @param assume_equal_variance Defaults to FALSE
#'
#' @return Returns a list with these named elements: 
#' * effect_size - the point estimate from the sample 
#' * lower - lower bound of the CI
#' * upper - upper bound of the CI 
#' * df - degrees of freedom for the contrast
#' * se - standard error
#' * moe - margin of error
#' * variability_component - if equal variance is assumed, returns sd_pooled;
#'     otherwise returns sd_avg
#' 
#' This function only takes in summary data and provides only the 
#'   minimal output.  For a friendlier version of this function, 
#'   see \code{\link{estimate_mdiff_contrast_bs}}
#'   
#' 
#' @section Details:
#' ## With equal variance assumed: ##
#' If equal variance is assumed, the standard error is calculated as given in
#' Kline, equation 7.8:
#' \mjdeqn{ 
#' s_{ \psi } = \sqrt{sd_{pooled} * (\sum_{i=1}^{a}c_i^2/n_i) } 
#' }{s_psi = sd_pooled * sqrt(sum(contrast^2/ns))}
#' where n is each group size and sd_pooled is the pooled standard deviation
#' as given in Kline, equation 3.11:
#' \mjdeqn{
#' sd_{pooled} = { \frac{ \sum_{i=1}^{a} (n_i -1) s_i^2 } 
#' { \sum_{i=1}^{a} (n_i-1) } }
#' }{sqrt(sum(variances*dfs) / sum(dfs))}
#' 
#' where *s* is the standard deviation for each group.
#' 
#' Degrees of freedom for the analysis are calculated as given in Kline,
#' equation 3.9:
#' \mjdeqn{
#' df = \sum_{i=1}^{a}{n_i -1}
#' }{sum(n-1)}
#' 
#'   
#'
#' ## Without assuming equal variance: ##
#' If equal variance is not assumed, the standard error is calculated as given
#' in Baguley, equation 15.11:
#' \mjdeqn{
#' s_{ \psi } = \sqrt{ \sum_{i=1}^{a} \frac{ c_i^2 s_i^2 }{ n_i }  }
#' }{s_psi = sqrt(sum(contrast^2 * variances / ns))}
#' 
#' The degrees of freedom for the analysis are calculated as in Baguley,
#' equation 15.12:
#' \mjdeqn{
#' df = \frac { s_{ \psi }^4 }
#' { \sum_{i=1}^{a} \frac{c_i^4 s_i^4} {n_i^2 (n_i -1)} }
#' }{se_psi^4 / sum( (contrast^4) * (sds^4) / (ns^2*dfs))}

#' 
#' And the analysis also returns sd_avg, the square root of the average of 
#' the group variances, as given in Bonett, explanation of equation 6:
#' \mjdeqn{
#' sd_{avg} = \sqrt{ \frac{ \sum_{i=1}^{a} s_i^2 }{ a } }
#' }{sqrt(mean(variances))}
#' 
#' 
#' @references 
#' * Baguley, T. (2012). Serious Stats. In Serious stats: A guide to advanced 
#' statistics for the behavioral sciences. Macmillan Education UK.
#'  https://doi.org/10.1007/978-0-230-36355-7
#' * Bonett, D. G. (2018).  R code posted to personal website.
#' https://people.ucsc.edu/~dgbonett/psyc204.html
#' * Bonett, D. G. (2008). Confidence Intervals for Standardized Linear 
#' Contrasts of Means. *Psychological Methods*, 13(2), 99–109. 
#' https://doi.org/10.1037/1082-989X.13.2.99
#' * Kline, R. B. (2013). *Beyond significance testing: Statistics reform in 
#' the behavioral sciences* (2nd ed.). American Psychological Association. 
#' https://doi.org/10.1037/14136-000
#'
#'
#' @seealso
#' * \code{\link{estimate_mdiff_contrast_bs}} for friendly version of
#'   this function 
#' * \code{\link{CI_smd_contrast_bs}} to obtain a standardized
#'   mean difference
#'
#'  
#' @examples
#' # Example from Kline, 2013
#' #  Data in Table 3.4
#' #  Worked out in Chapter 7
#' #  See note to table 7.2
#' # With equal variance assumed and no correction, should give:
#' #   psi = -2 95% CI [-5.23, 1.23]  (see note to Table 7.2)
#'
#' CI_mdiff_contrast_bs(
#'   means = c(13, 11, 15),
#'   sds = c(2.738613, 2.236068, 2.000000),
#'   ns = c(5, 5, 5),
#'   contrast = contrast <- c(1, 0, -1),
#'   conf_level = 0.95,
#'   assume_equal_variance = TRUE
#' )
#' 
#' # Example from Kline, 2013
#' #  Data in Table 3.4
#' #  Worked out in Chapter 7
#' #  See note to table 7.2
#' # With equal variance assumed and no correction, should give:
#' #   psi = 3 95% CI [0.20, 5.80]  (see note to Table 7.2)
#'
#' CI_mdiff_contrast_bs(
#'   means = c(13, 11, 15),
#'   sds = c(2.738613, 2.236068, 2.000000),
#'   ns = c(5, 5, 5),
#'   contrast = contrast <- c(1/2, -1, 1/2),
#'   conf_level = 0.95,
#'   assume_equal_variance = TRUE
#' )
#' 
#' # Example from Bonett, 2018, ci.lc.mean.bs, 
#' #  https://people.ucsc.edu/~dgbonett/psyc204.html
#' # Should give:
#' #                              Estimate       SE       df        LL        UL
#' # Equal Variances Assumed:        -5.35 1.300136 36.00000 -7.986797 -2.713203
#' # Equal Variances Not Assumed:    -5.35 1.300136 33.52169 -7.993583 -2.706417
#'
#' CI_mdiff_contrast_bs(
#'   means = c(33.5, 37.9, 38.0, 44.1),
#'   sds = c(3.84, 3.84, 3.65, 4.98),
#'   ns = c(10,10,10,10),
#'   contrast = contrast <- c(.5, .5, -.5, -.5),
#'   conf_level = 0.95,
#'   assume_equal_variance = FALSE
#' )
#'
#' @export
CI_mdiff_contrast_bs <- function(
  means,
  sds,
  ns,
  contrast,
  conf_level = 0.95,
  assume_equal_variance = FALSE
) {

  # Global calculations ------------------------------------------------------
  # n_groups - number of different groups passed
  n_groups <- length(means)
  
  # degrees of freedom for each mean
  dfs <- ns-1

  # pooled and average standard deviations
  # sd_pooled from Kline, 2013, equation 3.11
  # sd_avg from Bonett, 2008, explanation of equation 6
  variances <- sds^2
  sd_pooled <- sqrt(sum(variances*dfs) / sum(dfs))
  sd_avg <- sqrt(mean(variances))
  
  
  # Do the analysis -----------------------------------------------------------
  # Effect size - Kline, 2013, equation 7.2
  effect_size <- sum(means*contrast)
  
  # Standard error; save se for contrast
  se <- if(assume_equal_variance)
    # Kline, 2013, equation 7.8
    sd_pooled * sqrt(sum(contrast^2/ns))
  else
    # Baguley, 2012, equation 15.11
    sqrt(sum(contrast^2 * variances / ns))
  
  # Degrees of freedom
  df <- if(assume_equal_variance)
    # Kline, 2013, explanation of equation 7.8
    sum(dfs)
  else
    # Baguley, 2012, equation 15.12
    se^4 / sum( (contrast^4) * (sds^4) / (ns^2*dfs))
  
  # t multiplier
  t_multiplier <- stats::qt(1-((1 - conf_level)/2), df)
  
  # MoE and confidence interval
  moe <- se * t_multiplier
  lower <- effect_size - moe
  upper <- effect_size + moe
  
  
  # Prepare results -----------------------------------------------------------
  res <- list(
    effect_size = effect_size,
    lower = lower,
    upper = upper,
    df = df,
    se = se,
    moe = moe,
    variability_component = if(assume_equal_variance) sd_pooled else sd_avg
  )
  
  return(res)
}