|
1 | 1 | #' TOST function for a one-sample t-test (Cohen's d) |
| 2 | +#' @description |
| 3 | +#' `r lifecycle::badge('superseded')` |
| 4 | +#' |
| 5 | +#' Development on this function is complete, and for new code we recommend |
| 6 | +#' switching to [tsum_TOST], which is easier to use, more featureful, |
| 7 | +#' and still under active development. |
| 8 | +#' |
2 | 9 | #' @param m mean |
3 | 10 | #' @param mu value to compare against |
4 | 11 | #' @param sd standard deviation |
5 | 12 | #' @param n sample size |
6 | 13 | #' @param low_eqbound_d lower equivalence bounds (e.g., -0.5) expressed in standardized mean difference (Cohen's d) |
7 | 14 | #' @param high_eqbound_d upper equivalence bounds (e.g., 0.5) expressed in standardized mean difference (Cohen's d) |
| 15 | +#' @param low_eqbound lower equivalence bounds (e.g., -0.5) expressed in raw units |
| 16 | +#' @param high_eqbound upper equivalence bounds (e.g., 0.5) expressed in raw units |
8 | 17 | #' @param alpha alpha level (default = 0.05) |
9 | 18 | #' @param plot set whether results should be plotted (plot = TRUE) or not (plot = FALSE) - defaults to TRUE |
10 | 19 | #' @param verbose logical variable indicating whether text output should be generated (verbose = TRUE) or not (verbose = FALSE) - default to TRUE |
|
18 | 27 | #' @export |
19 | 28 | #' |
20 | 29 |
|
21 | | -TOSTone<-function(m,mu,sd,n,low_eqbound_d, high_eqbound_d, alpha, plot = TRUE, verbose = TRUE){ |
22 | | - message("Note: this function is defunct. Please use tsum_TOST instead") |
| 30 | +TOSTone <- |
| 31 | + function(m, |
| 32 | + mu, |
| 33 | + sd, |
| 34 | + n, |
| 35 | + low_eqbound_d, |
| 36 | + high_eqbound_d, |
| 37 | + alpha, |
| 38 | + plot = TRUE, |
| 39 | + verbose = TRUE) { |
| 40 | + lifecycle::deprecate_soft("0.4.0", "TOSTone()", "tsum_TOST()") |
23 | 41 | if(missing(alpha)) { |
24 | 42 | alpha<-0.05 |
25 | 43 | } |
@@ -113,3 +131,107 @@ TOSTone<-function(m,mu,sd,n,low_eqbound_d, high_eqbound_d, alpha, plot = TRUE, v |
113 | 131 | # Return results in list() |
114 | 132 | 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)) |
115 | 133 | } |
| 134 | + |
| 135 | +#' @rdname TOSTone |
| 136 | +#' @export |
| 137 | + |
| 138 | +TOSTone.raw <- |
| 139 | + function(m, |
| 140 | + mu, |
| 141 | + sd, |
| 142 | + n, |
| 143 | + low_eqbound, |
| 144 | + high_eqbound, |
| 145 | + alpha, |
| 146 | + plot = TRUE, |
| 147 | + verbose = TRUE) { |
| 148 | + lifecycle::deprecate_soft("0.4.0", "TOSTone.raw()", "tsum_TOST()") |
| 149 | + if(missing(alpha)) { |
| 150 | + alpha<-0.05 |
| 151 | + } |
| 152 | + 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.") |
| 153 | + if(n < 2) stop("The sample size should be larger than 1.") |
| 154 | + if(1 <= alpha | alpha <= 0) stop("The alpha level should be a positive value between 0 and 1.") |
| 155 | + if(sd <= 0) stop("The standard deviation should be a positive value.") |
| 156 | + # Calculate TOST, t-test, 90% CIs and 95% CIs |
| 157 | + degree_f<-n-1 |
| 158 | + t1<-(m-mu-low_eqbound)/(sd/sqrt(n))# t-test |
| 159 | + p1<-pt(t1, degree_f, lower.tail=FALSE) |
| 160 | + t2<-(m-mu-high_eqbound)/(sd/sqrt(n)) #t-test |
| 161 | + p2<-pt(t2, degree_f, lower.tail=TRUE) |
| 162 | + t<-(m-mu)/(sd/sqrt(n)) |
| 163 | + pttest<-2*pt(-abs(t), df=degree_f) |
| 164 | + LL90<-(m-mu)-qt(1-alpha, degree_f)*(sd/sqrt(n)) |
| 165 | + UL90<-(m-mu)+qt(1-alpha, degree_f)*(sd/sqrt(n)) |
| 166 | + LL95<-(m-mu)-qt(1-(alpha/2), degree_f)*(sd/sqrt(n)) |
| 167 | + UL95<-(m-mu)+qt(1-(alpha/2), degree_f)*(sd/sqrt(n)) |
| 168 | + ptost<-max(p1,p2) #Get highest p-value for summary TOST result |
| 169 | + ttost<-ifelse(abs(t1) < abs(t2), t1, t2) #Get lowest t-value for summary TOST result |
| 170 | + results<-data.frame(t1,p1,t2,p2,degree_f,LL90,UL90) |
| 171 | + colnames(results) <- c("t-value 1","p-value 1","t-value 2","p-value 2","df", paste("Lower Limit ",100*(1-alpha*2),"% CI",sep=""),paste("Upper Limit ",100*(1-alpha*2),"% CI",sep="")) |
| 172 | + dif<-(m-mu) |
| 173 | + testoutcome<-ifelse(pttest<alpha,"significant","non-significant") |
| 174 | + TOSToutcome<-ifelse(ptost<alpha,"significant","non-significant") |
| 175 | + |
| 176 | + # Plot results |
| 177 | + if (plot == TRUE) { |
| 178 | + plot(NA, ylim=c(0,1), xlim=c(min(LL95,low_eqbound,dif)-max(UL95-LL95, high_eqbound-low_eqbound,dif)/10, max(UL95,high_eqbound,dif)+max(UL95-LL95, high_eqbound-low_eqbound, dif)/10), bty="l", yaxt="n", ylab="",xlab="Mean Difference") |
| 179 | + points(x=dif, y=0.5, pch=15, cex=2) |
| 180 | + abline(v=high_eqbound, lty=2) |
| 181 | + abline(v=low_eqbound, lty=2) |
| 182 | + abline(v=0, lty=2, col="grey") |
| 183 | + segments(LL90,0.5,UL90,0.5, lwd=3) |
| 184 | + segments(LL95,0.5,UL95,0.5, lwd=1) |
| 185 | + 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) |
| 186 | + } |
| 187 | + |
| 188 | + if(missing(verbose)) { |
| 189 | + verbose <- TRUE |
| 190 | + } |
| 191 | + if(verbose == TRUE){ |
| 192 | + cat("TOST results:\n") |
| 193 | + 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)) |
| 194 | + cat("\n") |
| 195 | + 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)) |
| 196 | + cat("\n") |
| 197 | + cat("degrees of freedom :",round(degree_f, digits = 2)) |
| 198 | + cat("\n\n") |
| 199 | + cat("Equivalence bounds (raw scores):") |
| 200 | + cat("\n") |
| 201 | + cat("low eqbound:", paste0(round(low_eqbound, digits = 4)),"\nhigh eqbound:",paste0(round(high_eqbound, digits = 4))) |
| 202 | + cat("\n\n") |
| 203 | + cat("TOST confidence interval:") |
| 204 | + cat("\n") |
| 205 | + 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 = "") |
| 206 | + cat("\n\n") |
| 207 | + cat("NHST confidence interval:") |
| 208 | + cat("\n") |
| 209 | + cat("lower bound ",100*(1-alpha),"% CI: ", paste0(round(LL95, digits = 3)),"\nupper bound ",100*(1-alpha),"% CI: ",paste0(round(UL95,digits = 3)), sep = "") |
| 210 | + cat("\n\n") |
| 211 | + cat("Equivalence Test Result:\n") |
| 212 | + 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="")) |
| 213 | + cat("\n") |
| 214 | + cat("Null Hypothesis Test Result:\n") |
| 215 | + 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="")) |
| 216 | + if(pttest <= alpha && ptost <= alpha){ |
| 217 | + combined_outcome <- paste0("NHST: reject null significance hypothesis that the effect is equal to ", 0," \n", |
| 218 | + "TOST: reject null equivalence hypothesis") |
| 219 | + } |
| 220 | + if(pttest < alpha && ptost > alpha){ |
| 221 | + combined_outcome <- paste0("NHST: reject null significance hypothesis that the effect is equal to ",0," \n", |
| 222 | + "TOST: don't reject null equivalence hypothesis") |
| 223 | + } |
| 224 | + if(pttest > alpha && ptost <= alpha){ |
| 225 | + combined_outcome <- paste0("NHST: don't reject null significance hypothesis that the effect is equal to ",0," \n", |
| 226 | + "TOST: reject null equivalence hypothesis") |
| 227 | + } |
| 228 | + if(pttest > alpha && ptost > alpha){ |
| 229 | + combined_outcome <- paste0("NHST: don't reject null significance hypothesis that the effect is equal to ",0," \n", |
| 230 | + "TOST: don't reject null equivalence hypothesis") |
| 231 | + } |
| 232 | + cat("\n") |
| 233 | + message(combined_outcome) |
| 234 | + } |
| 235 | + # Return results in list(),NHST_t = t, NHST_p = pttest |
| 236 | + 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)) |
| 237 | + } |
0 commit comments