d <- read.csv("hearts.csv")
d$group <- factor(d$group, levels = c("control", "treated"))
d$ratio <- d$heart_mg / d$body_weight_g
fmt <- function(p) if (p < 0.001) "< 0.001" else sprintf("= %.3f", p)
labels <- c(body_weight_g = "Body weight (g)", heart_mg = "Heart weight (mg)", ratio = "Heart/body weight (mg/g)")
for (col in names(labels)) {
  tr <- d[[col]][d$group == "treated"]; co <- d[[col]][d$group == "control"]
  r <- t.test(tr, co)
  cat(sprintf("%s: control %.2f, treated %.2f, Welch t = %.2f, df = %.1f, p %s\n", labels[[col]], mean(co), mean(tr), r$statistic, r$parameter, fmt(r$p.value)))
}
slopes <- lm(log(heart_mg) ~ group * log(body_weight_g), data = d)
cat(sprintf("Slopes differ between groups? interaction p %s\n", fmt(coef(summary(slopes))["grouptreated:log(body_weight_g)", 4])))
m <- lm(log(heart_mg) ~ group + log(body_weight_g), data = d)
ci <- confint(m); cs <- coef(summary(m))
cat(sprintf("ANCOVA (log scale): common slope = %.2f (95%% CI %.2f to %.2f)\n", coef(m)[3], ci[3, 1], ci[3, 2]))
cat(sprintf("Treated vs control at the same body weight: ratio of heart weights %.3f (95%% CI %.3f to %.3f), t = %.2f, df = %d, p %s\n",
            exp(coef(m)[2]), exp(ci[2, 1]), exp(ci[2, 2]), cs[2, 3], m$df.residual, fmt(cs[2, 4])))