R/esci_testing.R

Defines functions plot_esci_test test_mdiff_contrast_bs

Documented in plot_esci_test test_mdiff_contrast_bs

#' conduct a hypothesis test on a linear contrast of means
#' 
#' @param estimate - must be an esci_estimate object generated by
#'   estimate_mdiff_contrast_bs
#' @param rope_lower - numeric lower limit of the region of practical 
#'   equivalence (rope)
#' @param rope_upper - numeric upper limit of the rope
#' @param rope_units - conduct test on original units (raw) or sds (cohen's d);
#'  defaults to raw
#' @param alpha - alpha level to use for evaluation
#' 
#' @return Returns an esci_test object
#' 
#' @export
test_mdiff_contrast_bs <- function(estimate, 
                      rope_lower = -0.1,
                      rope_upper = 0.1, 
                      rope_units = c("raw", "sd"),
                      alpha = 0.05
                      ) {
  

  # Input Checks ---------------------
  # This function expects:
  #   estimate should be of class estimate
  #   rope_lower <= 0 
  #   rope_upper >= 0 
  #   If rope_lower and rope_upper = 0, only traditional nil test returned
  #   0 < alpha < 1
  #   rope_units = "sd" | "raw"
  esci_assert_type(estimate, "is.estimate")
  esci_assert_range(
    var = rope_lower, 
    upper = 0, 
    upper_inclusive = TRUE
  )
  esci_assert_range(
    var = rope_upper, 
    lower = 0, 
    lower_inclusive = TRUE
    )
  esci_assert_range(
    var = alpha, 
    lower = 0, 
    upper = 1, 
    lower_inclusive = FALSE, 
    upper_inclusive = FALSE
    )
  rope_units <- match.arg(rope_units)
  
  
  # Prep ------------------------------------------
  res <- list()
  
  res$properties <- list(
    rope_lower = rope_lower,
    rope_upper = rope_upper,
    rope_units = rope_units,
    alpha = alpha
  )
  
  # Reduce down to differences
  effect_sizes <- estimate$effect_sizes[
    estimate$effect_sizes$type == "Difference", 
  ]
  
  # Loop through all effect sizes
  for (my_row in 1:nrow(effect_sizes)) {
    
    # Get single effect size and contrast to work with
    es <- as.list(effect_sizes[my_row, ])
    
    if (is.null(es$outcome_variable_name)) {
      outcome_variable_name <- estimate$properties$outcome_variable_name
    } else {
      outcome_variable_name <- es$outcome_variable_name
    }
    
    #### Set rope limits
    if (rope_units == "sd") {
      this_rope_upper <- rope_upper * es$variability_component
      this_rope_lower <- rope_lower * es$variability_component
    } else {
      this_rope_upper <- rope_upper
      this_rope_lower <- rope_lower
    }
  
    # Nil hypothesis test ------------------------------
    df <- es$df
    t_nil <- es$effect_size / es$se
    if (is.na(df))
      p_nil <- 2*pnorm(-abs(t_nil))
    else
      p_nil <- 2*pt(-abs(t_nil), df=df)
    
    if (p_nil < alpha) {
      if (es$effect_size > 0) {
       my_operator <- ">" 
      } else {
       my_operator <- "<"
      }
      conclusion_nil <- glue::glue(
        "Mu_diff is {my_operator} 0"
      )  
    } else {
      conclusion_nil <- glue::glue(
        "Sign of Mu_diff is ambiguous"
      )
    }
  
    nil_result <- list(
      outcome_variable_name = outcome_variable_name, 
      effect_size_label = es$effect,
      null_hypothesis = glue::glue("Mu_diff is exactly 0"),
      t = t_nil,
      df = df,
      p = p_nil,
      alpha = alpha,
      conclusion = conclusion_nil
    )
    
    res$hypothesis_evaluations <- rbind(
      res$hypothesis_evaluations,
      as.data.frame(nil_result)
    )
    
    if (!(rope_lower == 0 & rope_upper == 0)) {
      t_upper <- (es$effect_size - this_rope_upper) / es$se
      t_lower <- (es$effect_size - this_rope_lower) / es$se
      if (is.na(df)) {
        eq_p_upper <- pnorm(t_upper, lower.tail=TRUE)
        eq_p_lower <- pnorm(t_lower, lower.tail=FALSE)
      } else {
        eq_p_upper <- pt(t_upper, df, lower.tail=TRUE)
        eq_p_lower <- pt(t_lower, df, lower.tail=FALSE)
        
      }
      
      if (abs(t_upper) > abs(t_lower)) {
        t_report <- t_upper
      } else {
        t_report <- t_lower
      }
      
      if (eq_p_upper < alpha & eq_p_lower < alpha) {
        conclusion_eq <- "All compatible values of Mu_diff are negligble"
      } else {
        conclusion_eq <- "Cannot rule out substantive values of Mu_diff"
      }
      
      eq_result <- list(
        outcome_variable_name = outcome_variable_name,
        effect_size_label = es$effect,
        null_hypothesis = glue::glue("
Mu_diff is substantive, 
outside of range {round(this_rope_lower, 2)} to {round(this_rope_upper, 2)}"
        ),
        t = t_report,
        df = df,
        p = max(eq_p_upper, eq_p_lower),
        alpha = alpha,
        conclusion = conclusion_eq
      )
      
      if (is.na(df)) {
        me_p_upper <- pnorm(t_upper, lower.tail=FALSE)
        me_p_lower <- pnorm(t_lower, lower.tail=TRUE)
      } else {
        me_p_upper <- pt(t_upper, df, lower.tail=FALSE)
        me_p_lower <- pt(t_lower, df, lower.tail=TRUE)
      }

      
      if (me_p_upper < alpha/2 | me_p_lower < alpha/2) {
        conclusion_me <- "All compatible values of Mu_diff are substantive"
      } else {
        conclusion_me <- "Cannot rule out neglible values of Mu_diff"
      }
      
      me_result <- list(
        outcome_variable_name = outcome_variable_name,
        effect_size_label = es$effect,
        null_hypothesis = glue::glue(
"Mu_diff is neglible, 
inside the range {round(this_rope_lower, 2)} to {round(this_rope_upper, 2)}"),
        t = t_report,
        df = df,
        p = min(me_p_upper, me_p_lower),
        alpha = alpha/2,
        conclusion = conclusion_me
      )
      
      res$hypothesis_evaluations <- rbind(
        res$hypothesis_evaluations,
        as.data.frame(eq_result),
        as.data.frame(me_result)
      )
    } # end of equivalence and minimal effect tests  
    
    
  } # end of looping through effect sizes
    
  class(res) <- "esci_test"
  return(res)
}

#' visual display of a  hypothesis test on a linear contrast of means
#' 
#' @param estimate - must be an esci_estimate object generated by
#'   estimate_mdiff_contrast_bs
#' @param rope_lower - numeric lower limit of the region of practical 
#'   equivalence (rope)
#' @param rope_upper - numeric upper limit of the rope
#' @param rope_units - conduct test on original units (raw) or sds (cohen's d);
#'  defaults to raw
#' @param alpha - alpha level to use for evaluation
#' @param error_layout - optional, tells shape to use for expected error
#'   distribution: halfeye (default), eye, gradient, or none
#' @param error_scale - optional number giving width allocated to expected error
#'   distribution
#' @param error_normalize - optional; the way error distribution height is
#'   normalized in the graph: groups (default), all, panels
#' @param ggtheme - an optional ggtheme object to apply to the graph.  Defaults
#'   to theme_classic
#'   
#'   
#' @return Returns an ggplot2 object
#' 
#' @export
plot_esci_test <- function(
  estimate, 
  rope_lower = -0.1,
  rope_upper = 0.1, 
  rope_units = c("sd", "raw"),
  alpha = 0.05,
  error_layout = c("halfeye", "eye", "gradient", "none"),
  error_scale = 0.3,
  error_normalize = c("groups", "all", "panels"),
  ggtheme = ggplot2::theme_classic()
) {
  
  # Input checks ------------------------------------------------
  rope_units <- match.arg(rope_units)
  error_layout <- match.arg(error_layout)
  error_normalize <- match.arg(error_normalize)
  if (is.null(ggtheme)) {ggtheme <- ggplot2::theme_classic()}
  
  
  # Data processing ---------------------------------------------
  # Reduce down to differences
  es_data <- if(rope_units == "raw")
      estimate$effect_sizes[estimate$effect_sizes$type == "Difference", ]
  else
      estimate$standardized_effect_sizes
  

  # Set CIs to mark
  my_widths <- if(rope_lower == rope_upper) 
    c((1-alpha))
  else
    c((1-alpha), 1-(2*alpha))
    
  # Set ropes
  es_data$rope_upper <- rope_upper
  es_data$rope_lower <- rope_lower

  # xlabel
  xlab <- if(rope_units == "raw")
    estimate$properties$effect_size_name_html
  else
    estimate$standardized_effect_size_properties$d_name_html
  xlab <- esci_effect_size_expression(xlab)
  
  ylab <- es_data[1, 'effect']
  
  # Outcomve variable name
  if (is.null(es_data$outcome_variable_name)) 
    es_data$outcome_variable_name <- estimate$properties$outcome_variable_name
  
  
  # Build the plot -----------------------------------------------  
  # Basic plot
  myplot <- ggplot2::ggplot(
    data = es_data,
    ggplot2::aes(
      x = effect_size,
      y = outcome_variable_name
      )
    ) + ggtheme
  
  
  # if (!is.null(es_data$outcome_variable_name)) {
  #   # Facet by outcome_variable_name
  #   myplot <- myplot + ggplot2::facet_grid(
  #     cols = ggplot2::vars(outcome_variable_name)
  #   )
  # }
  
  
  # # Effect size locations and sampling error 
  #   # fill = stat(x > rope_lower & x < rope_upper),
  # es_data$rl <- es_data$rope_lower
  # es_data$ru <- es_data$rope_upper
  
  error_glue <- 
    "myplot <- myplot + {error_call}(
    data = es_data,
    orientation = 'horizontal',
    ggplot2::aes(
      x = effect_size, 
      y = outcome_variable_name,
      dist = distributional::dist_student_t(
        df = df, 
        mu = effect_size, 
        sigma = se,
      ),
      l = {rope_lower},
      u = {rope_upper},
      fill = stat(x > l & x < u)
    ),
    scale={error_scale},
    .width = {my_widths},
    normalize = '{error_normalize}',
    alpha = 0.5
    )"
  error_call <- esci_plot_error_layouts(error_layout)
  error_expression <- parse(text = glue::glue(error_glue))
  myplot <- try(eval(error_expression))
  
  # Mark 0 point of no effect
  myplot <- myplot + ggplot2::geom_vline(
    xintercept = 0,
    linetype = "dashed"
  )
  
  # Mark rope
  myplot <- myplot + ggplot2::geom_vline(
    xintercept = rope_lower,
    linetype = "dotted"
  )
  myplot <- myplot + ggplot2::geom_vline(
    xintercept = rope_upper,
    linetype = "dotted"
  )
  
  myplot <- myplot + ggplot2::xlab(xlab)
  myplot <- myplot + ggplot2::ylab(ylab)

  return(myplot)
  
}





# test_mean <- function(estimate, 
#                       rope_lower = -0.1,
#                       rope_upper = 0.1, 
#                       rope_units = c("sd", "raw"),
#                       alpha = 0.05
#                       ) {
#   
#   
#   if(nrow(estimate$effect_size_estimates) > 1) {
#     test_list <- list()
#     for (myvar in estimate$effect_size_estimates$variable_name) {
#       res <- test_mean(estimate = esci_split_estimate(estimate, myvar),
#                                       rope_lower = rope_lower,
#                                       rope_upper = rope_upper,
#                                       rope_units = rope_units,
#                                       alpha = alpha)
#       test_list[["hypothesis_evaluations"]] <- rbind(
#         test_list[["hypothesis_evaluations"]], 
#         res$hypothesis_evaluations
#       )
#     }
#     class(test_list) <- "esci_test"
#     return(test_list)
#     
#   }
# 
#   
#   # Input Checks ---------------------
#   # This function expects:
#   #   estimate should be of class estimate
#   #   rope_lower <= 0 
#   #   rope_upper >= 0 
#   #   If rope_lower and rope_upper = 0, only traditional nil test returned
#   #   0 < alpha < 1
#   #   rope_units = "sd" | "raw"
#   esci_assert_type(estimate, "is.estimate")
#   esci_assert_range(var = rope_lower, upper = 0, upper_inclusive = TRUE)
#   esci_assert_range(var = rope_upper, lower = 0, lower_inclusive = TRUE)
#   esci_assert_range(
#     var = alpha, 
#     lower = 0, 
#     upper = 1, 
#     lower_inclusive = FALSE, 
#     upper_inclusive = FALSE
#   )
#   rope_units <- match.arg(rope_units)
#   
#   
#   es <- as.list(estimate$effect_size_estimates[1, ])
#   
#   
#   #### Set rope limits
#   if (rope_units == "sd") {
#     rope_upper <- rope_upper * es$s_comp
#     rope_lower <- rope_lower * es$s_comp
#   }   
# 
#   
#   
#   
#   df <- es$df
#   t_nil <- es$effect_size / es$se
#   p_nil <- 2*pt(-abs(t_nil), df=df)
#   if (p_nil < alpha) {
#     conclusion_nil <- glue::glue("mu is not exactly {es$reference}")
#   } else {
#     conclusion_nil <- "Sign of mu is ambiguous"
#   }
# 
#   nil_result <- list(
#     variable_name = es$variable_name,
#     null_hypothesis = glue::glue("mu is exactly {es$reference}"),
#     t = t_nil,
#     df = df,
#     p = p_nil,
#     alpha = alpha,
#     conclusion = conclusion_nil
#   )
#   
#   hypothesis_evaluations <- as.data.frame(nil_result)
#   
#   if (!(rope_lower == 0 & rope_upper == 0)) {
#     t_upper <- (es$effect_size - rope_upper) / es$se
#     t_lower <- (es$effect_size - rope_lower) / es$se
#     eq_p_upper <- pt(t_upper, df, lower.tail=TRUE)
#     eq_p_lower <- pt(t_lower, df, lower.tail=FALSE)
#     
#     if (abs(t_upper) > abs(t_lower)) {
#       t_report <- t_upper
#     } else {
#       t_report <- t_lower
#     }
#     
#     if (eq_p_upper < alpha & eq_p_lower < alpha) {
#       conclusion_eq <- "Data are only compatible with neglible values of mu"
#     } else {
#       conclusion_eq <- "Data are compatible with substantive values of mu"
#     }
#     
#     eq_result <- list(
#       variable_name = es$variable_name,
#       null_hypothesis = glue::glue(
# "mu is substantial, 
# outside of range {round(rope_lower+es$reference, 2)} to {round(rope_upper+es$reference, 2)}"
#       ),
#       t = t_report,
#       df = df,
#       p = max(eq_p_upper, eq_p_lower),
#       alpha = alpha,
#       conclusion = conclusion_eq
#     )
#     
#     me_p_upper <- pt(t_upper, df, lower.tail=FALSE)
#     me_p_lower <- pt(t_lower, df, lower.tail=TRUE)
#     
#     if (me_p_upper < alpha/2 | me_p_lower < alpha/2) {
#       conclusion_me <- "Data are only compatible with substantive values of mu"
#     } else {
#       conclusion_me <- "Data are compatible with neglible values of mu"
#     }
#     
#     me_result <- list(
#       variable_name = es$variable_name,
#       null_hypothesis = glue::glue(
# "mu is neglible, 
# inside the range {round(rope_lower+es$reference, 2)} to {round(rope_upper+es$reference, 2)}"
#       ),
#       t = t_report,
#       df = df,
#       p = min(me_p_upper, me_p_lower),
#       alpha = alpha/2,
#       conclusion = conclusion_me
#     )
#     
#     hypothesis_evaluations <- rbind(hypothesis_evaluations,
#                            as.data.frame(eq_result),
#                            as.data.frame(me_result)
#                            )
#   }  
#   
#   res <- list(hypothesis_evaluations = hypothesis_evaluations)
#   class(res) <- "esci_test"
#   return(res)
# }
# 
# 
# some_tests <- function() {
#   mydata <- data.frame(var1 = rnorm(n = 1000, m = 104, sd = 20),
#                        var2 = rnorm(n = 1000, m = 108, sd = 15),
#                        var3 = c(rnorm(n = 950, m = 95, sd = 15), rep(NA, 50))
#   )
#   estimate <- estimate_mean_one(
#     mydata, var1, var2, var3, reference_m = 100, conf.level = 0.95
#     )
#   plot_estimate(estimate, rope = list(reference = 100, upper = 96, lower = 104))
#   test <- test_mean(estimate, rope_lower = -4, rope_upper = 4, rope_units = "raw")
#   
#   test_mean(estimate, rope_lower = 0, rope_upper = 0, rope_units = "raw")
#   
#   # TOSTER::dataTOSTone(
#   #   data = mydata, 
#   #   vars = "var1", 
#   #   mu = 100, 
#   #   low_eqbound = -4/sd(mydata$var1), 
#   #   high_eqbound = 4/sd(mydata$var1), 
#   #   desc = TRUE, 
#   #   alpha = 0.05
#   # )
#   
#   test_mean(estimate, rope_lower = -1/5, rope_upper = 1/5, rope_units = "sd")
# 
#   # Meaningful errors
#   test_mean(estimate, rope_lower = 4, rope_upper = 4, rope_units = "raw")
#   test_mean(estimate, rope_lower = -4, rope_upper = -4, rope_units = "raw")
#   
#   
# 
# }
rcalinjageman/esci2 documentation built on Dec. 22, 2021, 1:02 p.m.