#' TOST function for an independent t-test (Cohen's d)
#' @description
#' `r lifecycle::badge('superseded')`
#'
#' Development on `TOSTtwo` is complete, and for new code we recommend
#' switching to [tsum_TOST], which is easier to use, more featureful,
#' and still under active development.
#'
#' @param m1 mean of group 1
#' @param m2 mean of group 2
#' @param sd1 standard deviation of group 1
#' @param sd2 standard deviation of group 2
#' @param n1 sample size in group 1
#' @param n2 sample size in group 2
#' @param low_eqbound lower equivalence bounds (e.g., -0.5) expressed in raw scale units (e.g., scalepoints)
#' @param high_eqbound upper equivalence bounds (e.g., 0.5) expressed in raw scale units (e.g., scalepoints)
#' @param low_eqbound_d lower equivalence bounds (e.g., -0.5) expressed in standardized mean difference (Cohen's d)
#' @param high_eqbound_d upper equivalence bounds (e.g., 0.5) expressed in standardized mean difference (Cohen's d)
#' @param alpha alpha level (default = 0.05)
#' @param var.equal logical variable indicating whether equal variances assumption is assumed to be TRUE or FALSE. Defaults to FALSE.
#' @param plot set whether results should be plotted (plot = TRUE) or not (plot = FALSE) - defaults to TRUE
#' @param verbose logical variable indicating whether text output should be generated (verbose = TRUE) or not (verbose = FALSE) - default to TRUE
#' @return Returns TOST t-value 1, TOST p-value 1, TOST t-value 2, TOST p-value 2, degrees of freedom, low equivalence bound, high equivalence bound, low equivalence bound in Cohen's d, high equivalence bound in Cohen's d, Lower limit confidence interval TOST, Upper limit confidence interval TOST
#' @importFrom stats pnorm pt qnorm qt
#' @importFrom graphics abline plot points segments title
#' @examples
#' ## Eskine (2013) showed that participants who had been exposed to organic
#' ## food were substantially harsher in their moral judgments relative to
#' ## those exposed to control (d = 0.81, 95% CI: [0.19, 1.45]). A
#' ## replication by Moery & Calin-Jageman (2016, Study 2) did not observe
#' ## a significant effect (Control: n = 95, M = 5.25, SD = 0.95, Organic
#' ## Food: n = 89, M = 5.22, SD = 0.83). Following Simonsohn's (2015)
#' ## recommendation the equivalence bound was set to the effect size the
#' ## original study had 33% power to detect (with n = 21 in each condition,
#' ## this means the equivalence bound is d = 0.48, which equals a
#' ## difference of 0.384 on a 7-point scale given the sample sizes and a
#' ## pooled standard deviation of 0.894). Using a TOST equivalence test
#' ## with default alpha = 0.05, not assuming equal variances, and equivalence
#' ## bounds of d = -0.43 and d = 0.43 is significant, t(182) = -2.69,
#' ## p = 0.004. We can reject effects larger than d = 0.43.
#'
#' TOSTtwo(m1=5.25,m2=5.22,sd1=0.95,sd2=0.83,n1=95,n2=89,low_eqbound_d=-0.43,high_eqbound_d=0.43)
#' @section References:
#' Berger, R. L., & Hsu, J. C. (1996). Bioequivalence Trials, Intersection-Union Tests and Equivalence Confidence Sets. Statistical Science, 11(4), 283-302.
#'
#' Gruman, J. A., Cribbie, R. A., & Arpin-Cribbie, C. A. (2007). The effects of heteroscedasticity on tests of equivalence. Journal of Modern Applied Statistical Methods, 6(1), 133-140, formula for Welch's t-test on page 135
#' @export
#'
TOSTtwo <-
function(m1,
m2,
sd1,
sd2,
n1,
n2,
low_eqbound_d,
high_eqbound_d,
alpha,
var.equal,
plot = TRUE,
verbose = TRUE) {
lifecycle::deprecate_soft("0.4.0", "TOSTtwo()", "tsum_TOST()")
if(missing(alpha)) {
alpha<-0.05
}
if(missing(var.equal)) {
var.equal<-FALSE
}
if(low_eqbound_d >= high_eqbound_d) warning("The lower bound is equal to or larger than the upper bound. Check the plot and output to see if the bounds are specified as you intended.")
if(n1 < 2 | n2 < 2) stop("The sample size should be larger than 1.")
if(1 <= alpha | alpha <= 0) stop("The alpha level should be a positive value between 0 and 1.")
if(sd1 <= 0 | sd2 <= 0) stop("The standard deviation should be a positive value.")
# Calculate TOST, t-test, 90% CIs and 95% CIs
if(var.equal==TRUE) {
sdpooled<-sqrt((((n1 - 1)*(sd1^2)) + (n2 - 1)*(sd2^2))/((n1+n2)-2)) #calculate sd pooled
low_eqbound<-low_eqbound_d*sdpooled
high_eqbound<-high_eqbound_d*sdpooled
degree_f<-n1+n2-2
t1<-((m1-m2)-low_eqbound)/(sdpooled*sqrt(1/n1 + 1/n2)) #students t-test lower bound
p1<-pt(t1, degree_f, lower.tail=FALSE)
t2<-((m1-m2)-high_eqbound)/(sdpooled*sqrt(1/n1 + 1/n2)) #students t-test upper bound
p2<-pt(t2, degree_f, lower.tail=TRUE)
t<-(m1-m2)/(sdpooled*sqrt(1/n1 + 1/n2))
pttest<-2*pt(-abs(t), df=degree_f)
LL90<-(m1-m2)-qt(1-alpha, n1+n2-2)*(sdpooled*sqrt(1/n1 + 1/n2))
UL90<-(m1-m2)+qt(1-alpha, n1+n2-2)*(sdpooled*sqrt(1/n1 + 1/n2))
LL95<-(m1-m2)-qt(1-(alpha/2), n1+n2-2)*(sdpooled*sqrt(1/n1 + 1/n2))
UL95<-(m1-m2)+qt(1-(alpha/2), n1+n2-2)*(sdpooled*sqrt(1/n1 + 1/n2))
} else {
sdpooled<-sqrt((sd1^2 + sd2^2)/2) #calculate sd root mean squared for Welch's t-test
low_eqbound<-low_eqbound_d*sdpooled
high_eqbound<-high_eqbound_d*sdpooled
degree_f<-(sd1^2/n1+sd2^2/n2)^2/(((sd1^2/n1)^2/(n1-1))+((sd2^2/n2)^2/(n2-1))) #degrees of freedom for Welch's t-test
t1<-((m1-m2)-low_eqbound)/sqrt(sd1^2/n1 + sd2^2/n2) #welch's t-test upper bound
p1<-pt(t1, degree_f, lower.tail=FALSE) #p-value for Welch's TOST t-test
t2<-((m1-m2)-high_eqbound)/sqrt(sd1^2/n1 + sd2^2/n2) #welch's t-test lower bound
p2<-pt(t2, degree_f, lower.tail=TRUE) #p-value for Welch's TOST t-test
t<-(m1-m2)/sqrt(sd1^2/n1 + sd2^2/n2) #welch's t-test NHST
pttest<-2*pt(-abs(t), df=degree_f) #p-value for Welch's t-test
LL90<-(m1-m2)-qt(1-alpha, degree_f)*sqrt(sd1^2/n1 + sd2^2/n2) #Lower limit for CI Welch's t-test
UL90<-(m1-m2)+qt(1-alpha, degree_f)*sqrt(sd1^2/n1 + sd2^2/n2) #Upper limit for CI Welch's t-test
LL95<-(m1-m2)-qt(1-(alpha/2), degree_f)*sqrt(sd1^2/n1 + sd2^2/n2) #Lower limit for CI Welch's t-test
UL95<-(m1-m2)+qt(1-(alpha/2), degree_f)*sqrt(sd1^2/n1 + sd2^2/n2) #Upper limit for CI Welch's t-test
}
ptost<-max(p1,p2) #Get highest p-value for summary TOST result
ttost<-ifelse(abs(t1) < abs(t2), t1, t2) #Get lowest t-value for summary TOST result
dif<-(m1-m2)
testoutcome<-ifelse(pttest<alpha,"significant","non-significant")
TOSToutcome<-ifelse(ptost<alpha,"significant","non-significant")
# Plot results
if (plot == TRUE) {
plot(NA, ylim=c(0,1), xlim=c(min(LL90,low_eqbound)-max(UL90-LL90, high_eqbound-low_eqbound)/10, max(UL90,high_eqbound)+max(UL90-LL90, high_eqbound-low_eqbound)/10), bty="l", yaxt="n", ylab="",xlab="Mean Difference")
points(x=dif, y=0.5, pch=15, cex=2)
abline(v=high_eqbound, lty=2)
abline(v=low_eqbound, lty=2)
abline(v=0, lty=2, col="grey")
segments(LL90,0.5,UL90,0.5, lwd=3)
segments(LL95,0.5,UL95,0.5, lwd=1)
title(main=paste("Equivalence bounds ",round(low_eqbound,digits=3)," and ",round(high_eqbound,digits=3),"\nMean difference = ",round(dif,digits=3)," \n TOST: ", 100*(1-alpha*2),"% CI [",round(LL90,digits=3),";",round(UL90,digits=3),"] ", TOSToutcome," \n NHST: ", 100*(1-alpha),"% CI [",round(LL95,digits=3),";",round(UL95,digits=3),"] ", testoutcome,sep=""), cex.main=1)
}
if(missing(verbose)) {
verbose <- TRUE
}
if(verbose == TRUE){
cat("TOST results:\n")
cat("t-value lower bound:",format(t1, digits = 3, nsmall = 2, scientific = FALSE),"\tp-value lower bound:",format(p1, digits = 1, nsmall = 3, scientific = FALSE))
cat("\n")
cat("t-value upper bound:",format(t2, digits = 3, nsmall = 2, scientific = FALSE),"\tp-value upper bound:",format(p2, digits = 1, nsmall = 3, scientific = FALSE))
cat("\n")
cat("degrees of freedom :",round(degree_f, digits = 2))
cat("\n\n")
cat("Equivalence bounds (Cohen's d):")
cat("\n")
cat("low eqbound:", paste0(round(low_eqbound_d, digits = 4)),"\nhigh eqbound:",paste0(round(high_eqbound_d, digits = 4)))
cat("\n\n")
cat("Equivalence bounds (raw scores):")
cat("\n")
cat("low eqbound:", paste0(round(low_eqbound, digits = 4)),"\nhigh eqbound:",paste0(round(high_eqbound, digits = 4)))
cat("\n\n")
cat("TOST confidence interval:")
cat("\n")
cat("lower bound ",100*(1-alpha*2),"% CI: ", paste0(round(LL90, digits = 3)),"\nupper bound ",100*(1-alpha*2),"% CI: ",paste0(round(UL90,digits = 3)), sep = "")
cat("\n\n")
cat("NHST confidence interval:")
cat("\n")
cat("lower bound ",100*(1-alpha),"% CI: ", paste0(round(LL95, digits = 3)),"\nupper bound ",100*(1-alpha),"% CI: ",paste0(round(UL95,digits = 3)), sep = "")
cat("\n\n")
cat("Equivalence Test Result:\n")
message(cat("The equivalence test was ",TOSToutcome,", t(",round(degree_f, digits=2),") = ",format(ttost, digits = 3, nsmall = 3, scientific = FALSE),", p = ",format(ptost, digits = 3, nsmall = 3, scientific = FALSE),", given equivalence bounds of ",format(low_eqbound, digits = 3, nsmall = 3, scientific = FALSE)," and ",format(high_eqbound, digits = 3, nsmall = 3, scientific = FALSE)," (on a raw scale) and an alpha of ",alpha,".",sep=""))
cat("\n")
cat("Null Hypothesis Test Result:\n")
message(cat("The null hypothesis test was ",testoutcome,", t(",round(degree_f, digits=2),") = ",format(t, digits = 3, nsmall = 3, scientific = FALSE),", p = ",format(pttest, digits = 3, nsmall = 3, scientific = FALSE),", given an alpha of ",alpha,".",sep=""))
if(pttest <= alpha && ptost <= alpha){
combined_outcome <- paste0("NHST: reject null significance hypothesis that the effect is equal to ", 0," \n",
"TOST: reject null equivalence hypothesis")
}
if(pttest < alpha && ptost > alpha){
combined_outcome <- paste0("NHST: reject null significance hypothesis that the effect is equal to ",0," \n",
"TOST: don't reject null equivalence hypothesis")
}
if(pttest > alpha && ptost <= alpha){
combined_outcome <- paste0("NHST: don't reject null significance hypothesis that the effect is equal to ",0," \n",
"TOST: reject null equivalence hypothesis")
}
if(pttest > alpha && ptost > alpha){
combined_outcome <- paste0("NHST: don't reject null significance hypothesis that the effect is equal to ",0," \n",
"TOST: don't reject null equivalence hypothesis")
}
cat("\n")
message(combined_outcome)
}
# Return results in list()
invisible(list(diff=dif,TOST_t1=t1,TOST_p1=p1,TOST_t2=t2,TOST_p2=p2, TOST_df=degree_f,alpha=alpha,low_eqbound=low_eqbound,high_eqbound=high_eqbound,low_eqbound_d=low_eqbound_d,high_eqbound_d=high_eqbound_d, LL_CI_TOST=LL90,UL_CI_TOST=UL90,LL_CI_TTEST=LL95, UL_CI_TTEST=UL95,NHST_t = t, NHST_p = pttest))
}
#' @rdname TOSTtwo
#' @export
TOSTtwo.raw <-
function(m1,
m2,
sd1,
sd2,
n1,
n2,
low_eqbound,
high_eqbound,
alpha,
var.equal,
plot = TRUE,
verbose = TRUE) {
lifecycle::deprecate_soft("0.4.0", "TOSTtwo.raw()", "tsum_TOST()")
if (missing(alpha)) {
alpha <- 0.05
}
if (missing(var.equal)) {
var.equal <- FALSE
}
if (low_eqbound >= high_eqbound)
warning(
"The lower bound is equal to or larger than the upper bound. Check the plot and output to see if the bounds are specified as you intended."
)
if (n1 < 2 |
n2 < 2)
stop("The sample size should be larger than 1.")
if (1 <= alpha |
alpha <= 0)
stop("The alpha level should be a positive value between 0 and 1.")
if (sd1 <= 0 |
sd2 <= 0)
stop("The standard deviation should be a positive value.")
# Calculate TOST, t-test, 90% CIs and 95% CIs
if (var.equal == TRUE) {
sdpooled <-
sqrt((((n1 - 1) * (sd1 ^ 2)) + (n2 - 1) * (sd2 ^ 2)) / ((n1 + n2) - 2)) #calculate sd pooled
t1 <- ((m1 - m2) - low_eqbound) / (sdpooled * sqrt(1 / n1 + 1 / n2))
degree_f <- n1 + n2 - 2
p1 <- pt(t1, degree_f, lower.tail = FALSE)
t2 <- ((m1 - m2) - high_eqbound) / (sdpooled * sqrt(1 / n1 + 1 / n2))
p2 <- pt(t2, degree_f, lower.tail = TRUE)
LL90 <- (m1 - m2) - qt(1 - alpha, degree_f) * (sdpooled * sqrt(1 / n1 + 1 /
n2))
UL90 <- (m1 - m2) + qt(1 - alpha, degree_f) * (sdpooled * sqrt(1 / n1 + 1 /
n2))
LL95 <-
(m1 - m2) - qt(1 - (alpha / 2), degree_f) * (sdpooled * sqrt(1 / n1 + 1 /
n2))
UL95 <-
(m1 - m2) + qt(1 - (alpha / 2), degree_f) * (sdpooled * sqrt(1 / n1 + 1 /
n2))
t <- (m1 - m2) / (sdpooled * sqrt(1 / n1 + 1 / n2))
pttest <- 2 * pt(-abs(t), df = degree_f)
} else {
sdpooled <-
sqrt((sd1 ^ 2 + sd2 ^ 2) / 2) #calculate sd root mean squared for Welch's t-test
t1 <-
((m1 - m2) - low_eqbound) / sqrt(sd1 ^ 2 / n1 + sd2 ^ 2 / n2) #welch's t-test lower bound
degree_f <-
(sd1 ^ 2 / n1 + sd2 ^ 2 / n2) ^ 2 / (((sd1 ^ 2 / n1) ^ 2 / (n1 - 1)) + ((sd2 ^
2 / n2) ^ 2 / (n2 - 1))) #degrees of freedom for Welch's t-test
p1 <-
pt(t1, degree_f, lower.tail = FALSE) #p-value for Welch's TOST t-test
t2 <-
((m1 - m2) - high_eqbound) / sqrt(sd1 ^ 2 / n1 + sd2 ^ 2 / n2) #welch's t-test upper bound
p2 <-
pt(t2, degree_f, lower.tail = TRUE) #p-value for Welch's TOST t-test
t <- (m1 - m2) / sqrt(sd1 ^ 2 / n1 + sd2 ^ 2 / n2) #welch's t-test NHST
pttest <- 2 * pt(-abs(t), df = degree_f) #p-value for Welch's t-test
LL90 <-
(m1 - m2) - qt(1 - alpha, degree_f) * sqrt(sd1 ^ 2 / n1 + sd2 ^ 2 / n2) #Lower limit for CI Welch's t-test
UL90 <-
(m1 - m2) + qt(1 - alpha, degree_f) * sqrt(sd1 ^ 2 / n1 + sd2 ^ 2 / n2) #Upper limit for CI Welch's t-test
LL95 <-
(m1 - m2) - qt(1 - (alpha / 2), degree_f) * sqrt(sd1 ^ 2 / n1 + sd2 ^ 2 /
n2) #Lower limit for CI Welch's t-test
UL95 <-
(m1 - m2) + qt(1 - (alpha / 2), degree_f) * sqrt(sd1 ^ 2 / n1 + sd2 ^ 2 /
n2) #Upper limit for CI Welch's t-test
}
ptost <- max(p1, p2) #Get highest p-value for summary TOST result
ttost <-
ifelse(abs(t1) < abs(t2), t1, t2) #Get lowest t-value for summary TOST result
dif <- (m1 - m2)
testoutcome <- ifelse(pttest < alpha, "significant", "non-significant")
TOSToutcome <- ifelse(ptost < alpha, "significant", "non-significant")
# Plot results
if (plot == TRUE) {
plot(
NA,
ylim = c(0, 1),
xlim = c(
min(LL90, low_eqbound) - max(UL90 - LL90, high_eqbound - low_eqbound) /
10,
max(UL90, high_eqbound) + max(UL90 - LL90, high_eqbound - low_eqbound) /
10
),
bty = "l",
yaxt = "n",
ylab = "",
xlab = "Mean Difference"
)
points(
x = dif,
y = 0.5,
pch = 15,
cex = 2
)
abline(v = high_eqbound, lty = 2)
abline(v = low_eqbound, lty = 2)
abline(v = 0,
lty = 2,
col = "grey")
segments(LL90, 0.5, UL90, 0.5, lwd = 3)
segments(LL95, 0.5, UL95, 0.5, lwd = 1)
title(
main = paste(
"Equivalence bounds ",
round(low_eqbound, digits = 3),
" and ",
round(high_eqbound, digits = 3),
"\nMean difference = ",
round(dif, digits = 3),
" \n TOST: ",
100 * (1 - alpha * 2),
"% CI [",
round(LL90, digits = 3),
";",
round(UL90, digits = 3),
"] ",
TOSToutcome,
" \n NHST: ",
100 * (1 - alpha),
"% CI [",
round(LL95, digits = 3),
";",
round(UL95, digits = 3),
"] ",
testoutcome,
sep = ""
),
cex.main = 1
)
}
if (missing(verbose)) {
verbose <- TRUE
}
if (verbose == TRUE) {
cat("TOST results:\n")
cat(
"t-value lower bound:",
format(
t1,
digits = 3,
nsmall = 2,
scientific = FALSE
),
"\tp-value lower bound:",
format(
p1,
digits = 1,
nsmall = 3,
scientific = FALSE
)
)
cat("\n")
cat(
"t-value upper bound:",
format(
t2,
digits = 3,
nsmall = 2,
scientific = FALSE
),
"\tp-value upper bound:",
format(
p2,
digits = 1,
nsmall = 3,
scientific = FALSE
)
)
cat("\n")
cat("degrees of freedom :", round(degree_f, digits = 2))
cat("\n\n")
cat("Equivalence bounds (raw scores):")
cat("\n")
cat("low eqbound:",
paste0(round(low_eqbound, digits = 4)),
"\nhigh eqbound:",
paste0(round(high_eqbound, digits = 4)))
cat("\n\n")
cat("TOST confidence interval:")
cat("\n")
cat(
"lower bound ",
100 * (1 - alpha * 2),
"% CI: ",
paste0(round(LL90, digits = 3)),
"\nupper bound ",
100 * (1 - alpha * 2),
"% CI: ",
paste0(round(UL90, digits = 3)),
sep = ""
)
cat("\n\n")
cat("NHST confidence interval:")
cat("\n")
cat(
"lower bound ",
100 * (1 - alpha),
"% CI: ",
paste0(round(LL95, digits = 3)),
"\nupper bound ",
100 * (1 - alpha),
"% CI: ",
paste0(round(UL95, digits = 3)),
sep = ""
)
cat("\n\n")
cat("Equivalence Test Result:\n")
message(
cat(
"The equivalence test was ",
TOSToutcome,
", t(",
round(degree_f, digits = 2),
") = ",
format(
ttost,
digits = 3,
nsmall = 3,
scientific = FALSE
),
", p = ",
format(
ptost,
digits = 3,
nsmall = 3,
scientific = FALSE
),
", given equivalence bounds of ",
format(
low_eqbound,
digits = 3,
nsmall = 3,
scientific = FALSE
),
" and ",
format(
high_eqbound,
digits = 3,
nsmall = 3,
scientific = FALSE
),
" (on a raw scale) and an alpha of ",
alpha,
".",
sep = ""
)
)
cat("\n")
cat("Null Hypothesis Test Result:\n")
message(
cat(
"The null hypothesis test was ",
testoutcome,
", t(",
round(degree_f, digits = 2),
") = ",
format(
t,
digits = 3,
nsmall = 3,
scientific = FALSE
),
", p = ",
format(
pttest,
digits = 3,
nsmall = 3,
scientific = FALSE
),
", given an alpha of ",
alpha,
".",
sep = ""
)
)
if (pttest <= alpha && ptost <= alpha) {
combined_outcome <-
paste0(
"NHST: reject null significance hypothesis that the effect is equal to ",
0,
" \n",
"TOST: reject null equivalence hypothesis"
)
}
if (pttest < alpha && ptost > alpha) {
combined_outcome <-
paste0(
"NHST: reject null significance hypothesis that the effect is equal to ",
0,
" \n",
"TOST: don't reject null equivalence hypothesis"
)
}
if (pttest > alpha && ptost <= alpha) {
combined_outcome <-
paste0(
"NHST: don't reject null significance hypothesis that the effect is equal to ",
0,
" \n",
"TOST: reject null equivalence hypothesis"
)
}
if (pttest > alpha && ptost > alpha) {
combined_outcome <-
paste0(
"NHST: don't reject null significance hypothesis that the effect is equal to ",
0,
" \n",
"TOST: don't reject null equivalence hypothesis"
)
}
cat("\n")
message(combined_outcome)
}
# Return results in list()
invisible(
list(
diff = dif,
TOST_t1 = t1,
TOST_p1 = p1,
TOST_t2 = t2,
TOST_p2 = p2,
TOST_df = degree_f,
alpha = alpha,
low_eqbound = low_eqbound,
high_eqbound = high_eqbound,
LL_CI_TOST = LL90,
UL_CI_TOST = UL90,
LL_CI_TTEST = LL95,
UL_CI_TTEST = UL95,
NHST_t = t,
NHST_p = pttest
)
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.