d <- read.csv("density.csv")
d$genotype <- factor(d$genotype, levels = c("WT", "KO"))
fmt <- function(p) if (p < 0.001) "< 0.001" else sprintf("= %.3f", p)
for (g in c("WT", "KO")) {
  s <- d[d$genotype == g, ]
  r <- cor.test(s$age_weeks, s$density)
  m <- lm(density ~ age_weeks, data = s); ci <- confint(m)["age_weeks", ]
  cat(sprintf("%s: r = %.2f, p %s; slope = %.2f per week (95%% CI %.2f to %.2f), slope p %s, R^2 = %.2f\n",
              g, r$estimate, fmt(r$p.value), coef(m)[2], ci[1], ci[2], fmt(coef(summary(m))[2, 4]), summary(m)$r.squared))
}
m <- lm(density ~ age_weeks * genotype, data = d)
cs <- coef(summary(m))["age_weeks:genotypeKO", ]; ci <- confint(m)["age_weeks:genotypeKO", ]
cat(sprintf("Difference in slopes (KO - WT) = %.2f per week (95%% CI %.2f to %.2f), t = %.2f, df = %d, p %s\n",
            cs[1], ci[1], ci[2], cs[3], m$df.residual, fmt(cs[4])))
c <- lm(density ~ age_weeks + genotype, data = d); ci <- confint(c)["age_weeks", ]
cat(sprintf("Common slope = %.2f per week (95%% CI %.2f to %.2f)\n", coef(c)["age_weeks"], ci[1], ci[2]))
for (g in c("WT", "KO")) cat(sprintf("%s residual SD = %.2f\n", g, summary(lm(density ~ age_weeks, data = d[d$genotype == g, ]))$sigma))
rob <- lmtest::coeftest(m, vcov = sandwich::vcovHC(m, type = "HC3"))
cat(sprintf("Robust (HC3) test of the difference in slopes: p %s\n", fmt(rob["age_weeks:genotypeKO", 4])))