# Fear extinction, knockout vs wild-type: percent freezing per tone analysed five ways. R 4.3.3, base R.
d <- read.csv("freezing.csv")                        # one row per mouse and tone: genotype, mouse, tone, seconds_freezing, seconds_scored
f2 <- function(p) formatC(signif(p, 2), format = "g")
d$pct <- 100 * d$seconds_freezing / d$seconds_scored
Y <- tapply(d$pct, list(d$mouse, d$tone), mean)   # mouse IDs must be unique across genotypes
g <- tapply(d$genotype, d$mouse, `[`, 1)[rownames(Y)]
ko <- Y[g == "knockout", ]; wt <- Y[g == "wild-type", ]; k <- ncol(Y)
cat(sprintf("Freezing, wild-type: tone 1 %.1f%%, tone %d %.1f%%; knockout: tone 1 %.1f%%, tone %d %.1f%%\n",
            mean(wt[, 1]), k, mean(wt[, k]), mean(ko[, 1]), k, mean(ko[, k])))
p_tone <- sapply(1:k, function(j) t.test(ko[, j], wt[, j], var.equal = TRUE)$p.value)
hits <- if (any(p_tone < 0.05)) paste(which(p_tone < 0.05), collapse = ", ") else "none"
cat(sprintf("A. A t-test at each of the %d tones: %d with p < 0.05 (tones: %s); smallest p = %s\n", k, sum(p_tone < 0.05), hits, f2(min(p_tone))))
t <- t.test(as.vector(ko), as.vector(wt), var.equal = TRUE)
cat(sprintf("B. All tones pooled, tones as n (%d vs %d), t-test: knockout - wild-type %.1f points, p = %s\n", length(ko), length(wt), mean(ko) - mean(wt), f2(t$p.value)))
mk <- rowMeans(ko); mw <- rowMeans(wt); t <- t.test(mk, mw)
cat(sprintf("C. Mean freezing per mouse (%d vs %d), Welch t-test: knockout - wild-type %.1f points (95%% CI %.1f to %.1f), p = %s\n",
            length(mk), length(mw), mean(mk) - mean(mw), t$conf.int[1], t$conf.int[2], f2(t$p.value)))
x <- 1:k; slope <- function(m) apply(m, 1, function(y) coef(lm(y ~ x))[[2]])
sk <- slope(ko); sw <- slope(wt); t <- t.test(sk, sw)
cat(sprintf("D. Extinction slope per mouse (points per tone), Welch t-test: wild-type %.2f, knockout %.2f, difference %.2f (95%% CI %.2f to %.2f), p = %s\n",
            mean(sw), mean(sk), mean(sk) - mean(sw), t$conf.int[1], t$conf.int[2], f2(t$p.value)))
# Repeated-measures ANOVA, genotype x tone interaction, with and without the Greenhouse-Geisser correction
n <- length(mk); N <- 2 * n; Yall <- rbind(wt, ko); gm <- mean(Yall)
mg <- rbind(matrix(colMeans(wt), n, k, byrow = TRUE), matrix(colMeans(ko), n, k, byrow = TRUE))
rg <- rep(c(mean(wt), mean(ko)), each = n); rm_ <- rowMeans(Yall); ct <- matrix(colMeans(Yall), N, k, byrow = TRUE)
ss_int <- sum((mg - rg - ct + gm)^2); ss_err <- sum((Yall - rm_ - mg + rg)^2)
df_int <- k - 1; df_err <- (N - 2) * (k - 1); F <- (ss_int / df_int) / (ss_err / df_err)
S <- (cov(wt) + cov(ko)) / 2; H <- diag(k) - 1 / k; Sd <- H %*% S %*% H
eps <- sum(diag(Sd))^2 / ((k - 1) * sum(Sd * Sd))
cat(sprintf("E. Repeated-measures ANOVA, genotype x tone: F(%d, %d) = %.2f, p = %s; Greenhouse-Geisser epsilon %.2f, corrected p = %s\n",
            df_int, df_err, F, f2(pf(F, df_int, df_err, lower.tail = FALSE)), eps, f2(pf(F, eps * df_int, eps * df_err, lower.tail = FALSE))))