##############################################
## BOOTSTRAPPED WELCH T-TEST
## BEN MEULEMAN
## 2024-09-25
##############################################

# LOAD REQUIRED LIBRARY
library(boot)

# FUNCTION
bootwelch <- function(y,g,nsim=999) {
  tempdata <- data.frame(y,g)
  form <- as.formula("y~g")
  org <- t.test(form,data=tempdata,var.equal=FALSE)
  bootfunc <- function(d,i) { output <- t.test(form,data=d[i,],var.equal=FALSE) ; c(output$statistic,output$parameter) }
  result <- boot(tempdata,bootfunc,R=nsim,stype="i")
  boott <- boot.ci(result,type="bca",t0=result$t0[1],t=result$t[,1])  
  sd.av <- sqrt(mean(tapply(y,g,var)))

  out <- list(t.original=org$statistic,BCA=boott$bca[4:5],DF=c(mean(result$t[,2]),range(result$t[,2])),
              p.boot=(1-pt(abs(org$statistic),df=min(result$t[,2])))*2,cohen.d=-1*diff(org$estimate)/sd.av)

  cat("\n-------------------------------------------\n",
      "Bootstrapped Welch t-test ","(",nsim," resamples)\n",
      "-------------------------------------------\n\n",
      "Welch t: ",round(org$statistic,3),", BCA 95%CI: [",paste(round(boott$bca[4:5],3),collapse=", "),"]\n",
      "\n",
      "Minimal bootstrapped Welch DF: ",round(min(result$t[,2]),2)," (mean: ",round(mean(result$t[,2]),2),")\n",
      "Minimal theoretical p-value: ",round((1-pt(abs(org$statistic),df=min(result$t[,2])))*2,4),"\n\n",
      "Adjusted Cohen's d: ",round(diff(org$estimate)/sd.av,2),"\n",sep="")
  flush.console()
  invisible(out)
}