#' 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")
#
#
#
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.