# Worked example for /guides/when-can-you-exclude-an-outlier/ in base R 4.3.3 (constructed data).
# Run from this folder: Rscript example.R
control <- c(40, 44, 37, 46, 42, 39)
treated <- c(52, 58, 49, 61, 55, 118)

print(t.test(treated, control))                       # Welch, all values
print(t.test(treated[-6], control))                   # Welch, without T6
print(wilcox.test(treated, control, exact = TRUE))    # Mann-Whitney, all values

# Yuen's test on 20% trimmed means (Welch-type), written out from its definition
yuen <- function(x, y, tr = 0.2) {
  part <- function(v) {
    n <- length(v); g <- floor(tr * n); h <- n - 2 * g
    s <- sort(v); w <- s
    if (g > 0) { w[1:g] <- s[g + 1]; w[(n - g + 1):n] <- s[n - g] }
    list(tm = mean(v, trim = tr), d = (n - 1) * var(w) / (h * (h - 1)), h = h)
  }
  a <- part(x); b <- part(y)
  t <- (a$tm - b$tm) / sqrt(a$d + b$d)
  df <- (a$d + b$d)^2 / (a$d^2 / (a$h - 1) + b$d^2 / (b$h - 1))
  c(trimmed_diff = a$tm - b$tm, t = t, df = df, p = 2 * pt(-abs(t), df))
}
print(yuen(treated, control))

print(t.test(log(treated), log(control)))             # Welch on the log scale
lt <- t.test(log(treated), log(control))
cat("ratio of geometric means", exp(diff(rev(lt$estimate))), "CI", exp(lt$conf.int), "\n")

lw <- t.test(log(treated[-6]), log(control))           # planned analysis without T6 (sensitivity)
print(lw)
cat("ratio without T6", exp(diff(rev(lw$estimate))), "CI", exp(lw$conf.int), "\n")

z <- (treated - mean(treated)) / sd(treated)
cat("z scores", round(z, 3), "\n")
n <- length(treated)
cat("largest possible |z| for n = 6:", (n - 1) / sqrt(n), "\n")
tc <- qt(1 - 0.05 / (2 * n), n - 2)
cat("Grubbs G", max(abs(z)), "critical", (n - 1) / sqrt(n) * sqrt(tc^2 / (n - 2 + tc^2)), "\n")
cat("MAD scores", round((treated - median(treated)) / mad(treated), 3), "\n")