suppressPackageStartupMessages(library(afex))
d <- read.csv("gtt.csv")
d$group <- factor(d$group, levels = c("chow", "HFD"))
fmt <- function(p) if (p < 0.001) "< 0.001" else sprintf("= %.3f", p)
trap <- function(t, y) sum(diff(t) * (head(y, -1) + tail(y, -1)) / 2)   # trapezoids over the actual times
per <- do.call(rbind, lapply(split(d, d$mouse), function(g) {
  g <- g[order(g$time), ]
  data.frame(group = g$group[1], mouse = g$mouse[1], fasting = g$glucose[1],
             auc = trap(g$time, g$glucose), iauc = trap(g$time, g$glucose - g$glucose[1]))
}))
labels <- c(fasting = "Fasting glucose (mmol/L)", auc = "Total AUC (mmol/L x min)", iauc = "Incremental AUC (mmol/L x min)")
for (col in names(labels)) {
  r <- t.test(per[[col]][per$group == "HFD"], per[[col]][per$group == "chow"])   # Welch by default
  cat(sprintf("%s: chow %.1f, HFD %.1f, difference %.1f (95%% CI %.1f to %.1f), Welch t = %.2f, df = %.1f, p %s\n",
              labels[[col]], mean(per[[col]][per$group == "chow"]), mean(per[[col]][per$group == "HFD"]),
              diff(rev(r$estimate)), r$conf.int[1], r$conf.int[2], r$statistic, r$parameter, fmt(r$p.value)))
}
a <- suppressMessages(aov_ez("mouse", "glucose", d, between = "group", within = "time"))
s <- summary(a)
uni <- s$univariate.tests; gg <- s$pval.adjustments
cat(sprintf("Group (curve height): F(%.0f, %.0f) = %.2f, p %s\n", uni["group", "num Df"], uni["group", "den Df"], uni["group", "F value"], fmt(uni["group", "Pr(>F)"])))
cat(sprintf("Group x time: F(%.0f, %.0f) = %.2f, Greenhouse-Geisser epsilon = %.2f, corrected p = %.3f (uncorrected %.3f)\n",
            uni["group:time", "num Df"], uni["group:time", "den Df"], uni["group:time", "F value"],
            gg["group:time", "GG eps"], gg["group:time", "Pr(>F[GG])"], uni["group:time", "Pr(>F)"]))

# What a t-test at every time point would show (shown for comparison; not recommended)
for (t in sort(unique(d$time))) {
  h <- d$glucose[d$time == t & d$group == "HFD"]; c <- d$glucose[d$time == t & d$group == "chow"]
  cat(sprintf("%3d min: HFD - chow = %.2f mmol/L, Welch p %s\n", t, mean(h) - mean(c), fmt(t.test(h, c)$p.value)))
}