suppressPackageStartupMessages(library(multcomp))
d <- read.csv("doses.csv")
d$level_f <- factor(d$level)
means <- tapply(d$response, d$level, mean); n <- tapply(d$response, d$level, length)
cat(sprintf("Group means (vehicle, 1, 3, 10 mg/kg): %s\n", paste(sprintf("%.2f", means), collapse = ", ")))
set.seed(1)
fit <- aov(response ~ level_f, data = d)
dn <- summary(glht(fit, linfct = mcp(level_f = "Dunnett")))
cat(sprintf("Dunnett, each dose vs vehicle: p = %s\n", paste(sprintf("%.2f", dn$test$pvalues), collapse = ", ")))
cc <- c(-3, -1, 1, 3)
df <- sum(n) - length(n); mse <- sum(tapply(d$response, d$level, function(x) sum((x - mean(x))^2))) / df
est <- sum(cc * means); se <- sqrt(mse * sum(cc^2 / n)); t <- est / se; p <- 2 * pt(-abs(t), df)
cat(sprintf("Linear trend contrast: estimate = %.2f, t = %.2f, df = %d, p = %.4f\n", est, t, df, p))
k <- sum(cc * 0:3); tq <- qt(0.975, df)
cat(sprintf("Change per dose step (scores 0-3): %.2f (95%% CI %.2f to %.2f)\n", est / k, (est - tq * se) / k, (est + tq * se) / k))