# Numbers for /guides/kruskal-wallis-vs-anova/ in R (cross-check of example.py).
# R 4.3.3, rstatix 0.7.2, coin 1.4.3. Run from this folder after example.py:
#   Rscript example.R > example_R_output.txt
g <- read.csv("grades.csv")
g$dose <- factor(g$dose, levels = c("vehicle", "low", "high"))
cat("== Kruskal-Wallis, histology grade ==\n")
print(kruskal.test(grade ~ dose, data = g))
cat("== Dunn (Holm), rstatix ==\n")
print(as.data.frame(rstatix::dunn_test(g, grade ~ dose, p.adjust.method = "holm")))
cat("== Approximate permutation p (coin, 100000 resamples, seed 1) ==\n")
set.seed(1)
print(coin::kruskal_test(grade ~ dose, data = g,
                         distribution = coin::approximate(nresample = 100000)))
cat("== Separate Mann-Whitney tests (each pair re-ranked), Holm ==\n")
pw <- pairwise.wilcox.test(g$grade, g$dose, p.adjust.method = "holm", exact = FALSE)
print(pw$p.value)

d <- read.csv("il6.csv")
d$dose <- factor(d$dose, levels = c("vehicle", "low", "high"))
cat("== IL-6: classical ANOVA (raw) ==\n")
print(summary(aov(il6 ~ dose, data = d)))
cat("== IL-6: Welch ANOVA (raw) ==\n")
print(oneway.test(il6 ~ dose, data = d))
cat("== IL-6: Welch ANOVA (log) ==\n")
print(oneway.test(log(il6) ~ dose, data = d))
cat("== IL-6: Kruskal-Wallis ==\n")
print(kruskal.test(il6 ~ dose, data = d))
cat("== IL-6: ratio high/vehicle, Welch on logs (unadjusted 95% CI) ==\n")
tt <- t.test(log(il6) ~ dose, data = droplevels(subset(d, dose %in% c("high", "vehicle"))))
print(c(ratio = exp(diff(tt$estimate)), lower = exp(-tt$conf.int[2]), upper = exp(-tt$conf.int[1]),
        p = tt$p.value))