# Cosinor analysis (24-h period) by least squares, per genotype and compared between genotypes. R 4.3.3, base R.
d <- read.csv("rhythm.csv"); w <- 2 * pi / 24
d$c <- cos(w * d$zt); d$s <- sin(w * d$zt); d$ko <- as.numeric(d$genotype == "KO")
describe <- function(b1, b2, V) { A <- sqrt(b1^2 + b2^2); phi <- atan2(b2, b1)
  gA <- c(b1 / A, b2 / A); gP <- c(-b2 / A^2, b1 / A^2)
  c(A, sqrt(drop(t(gA) %*% V %*% gA)), (phi / w) %% 24, sqrt(drop(t(gP) %*% V %*% gP)) / w) }
for (g in c("WT", "KO")) {
  x <- d[d$genotype == g, ]; f1 <- lm(expr ~ c + s, data = x); f0 <- lm(expr ~ 1, data = x)
  b <- coef(f1); r <- describe(b[["c"]], b[["s"]], vcov(f1)[2:3, 2:3]); dfe <- df.residual(f1); tq <- qt(0.975, dfe)
  a <- anova(f0, f1)
  cat(sprintf("%s: mesor %.2f; amplitude %.3f (95%% CI %.3f to %.3f); peak at ZT %.1f (95%% CI %.1f to %.1f); zero-amplitude F(2, %d) = %.1f, p = %s\n",
              g, b[[1]], r[1], r[1] - tq * r[2], r[1] + tq * r[2], r[3], (r[3] - tq * r[4]) %% 24, (r[3] + tq * r[4]) %% 24, dfe, a$F[2], formatC(signif(a$`Pr(>F)`[2], 2), format = "g")))
}
ff <- lm(expr ~ c + s + ko + ko:c + ko:s, data = d); fr <- lm(expr ~ c + s + ko, data = d); a <- anova(fr, ff)
cat(sprintf("Rhythm differs between genotypes (KO x cos and KO x sin terms): F(2, %d) = %.2f, p = %s\n", df.residual(ff), a$F[2], formatC(signif(a$`Pr(>F)`[2], 2), format = "g")))
bf <- coef(ff); Vf <- vcov(ff); tq <- qt(0.975, df.residual(ff))
bw <- bf[c("c", "s")]; bk <- bw + bf[c("c:ko", "s:ko")]; Aw <- sqrt(sum(bw^2)); Ak <- sqrt(sum(bk^2))
nm <- names(bf); g <- setNames(rep(0, length(bf)), nm)
g[c("c", "s")] <- bk / Ak - bw / Aw; g[c("c:ko", "s:ko")] <- bk / Ak
sd <- sqrt(drop(t(g) %*% Vf %*% g))
cat(sprintf("Amplitude, KO - WT: %.3f (95%% CI %.3f to %.3f)\n", Ak - Aw, Ak - Aw - tq * sd, Ak - Aw + tq * sd))
dp <- (atan2(bk[2], bk[1]) - atan2(bw[2], bw[1]) + pi) %% (2 * pi) - pi
gp <- setNames(rep(0, length(bf)), nm)
gp[c("c", "s")] <- c(-bk[2], bk[1]) / Ak^2 - c(-bw[2], bw[1]) / Aw^2; gp[c("c:ko", "s:ko")] <- c(-bk[2], bk[1]) / Ak^2
sp <- sqrt(drop(t(gp) %*% Vf %*% gp)) / w
cat(sprintf("Peak time, KO - WT: %.1f h (95%% CI %.1f to %.1f)\n", dp / w, dp / w - tq * sp, dp / w + tq * sp))
pt <- sort(unique(d$zt))
ps <- sapply(pt, function(t) t.test(d$expr[d$zt == t & d$genotype == "KO"], d$expr[d$zt == t & d$genotype == "WT"], var.equal = TRUE)$p.value)
cat(sprintf("For comparison, separate t-tests at each of %d time points: p < 0.05 at ZT %s\n", length(pt), paste(pt[ps < 0.05], collapse = ", ")))
bon <- pt[ps < 0.05 / length(pt)]
cat(sprintf("With a Bonferroni correction for %d tests (p < 0.05/%d): ZT %s\n", length(pt), length(pt), if (length(bon)) paste(bon, collapse = ", ") else "none"))