d <- read.csv("pressure.csv")
d$group <- factor(d$group, levels = c("vehicle", "drug"))
d$change <- d$week4 - d$baseline
d$pct <- 100 * d$change / d$baseline
fmt <- function(p) if (p < 0.001) "< 0.001" else sprintf("= %.3f", p)
cat(sprintf("Baseline means: vehicle %.1f, drug %.1f; baseline-week 4 correlation r = %.2f\n",
            mean(d$baseline[d$group == "vehicle"]), mean(d$baseline[d$group == "drug"]), cor(d$baseline, d$week4)))
labels <- c(week4 = "Week 4 only", change = "Change from baseline", pct = "Percent change")
for (col in names(labels)) {
  a <- d[[col]][d$group == "drug"]; b <- d[[col]][d$group == "vehicle"]
  r <- t.test(a, b)
  cat(sprintf("%s: drug - vehicle = %.1f (95%% CI %.1f to %.1f), Welch t = %.2f, df = %.1f, p %s\n",
              labels[[col]], mean(a) - mean(b), r$conf.int[1], r$conf.int[2], r$statistic, r$parameter, fmt(r$p.value)))
}
m <- lm(week4 ~ group + baseline, data = d)
ci <- confint(m)["groupdrug", ]; cs <- coef(summary(m))["groupdrug", ]
cat(sprintf("ANCOVA (week 4 adjusted for baseline): drug - vehicle = %.1f (95%% CI %.1f to %.1f), t = %.2f, df = %d, p %s\n",
            cs[1], ci[1], ci[2], cs[3], m$df.residual, fmt(cs[4])))