# Daily clinical scores, vehicle vs drug, 10 mice each, days 7-28. R 4.3.3 with survival.
suppressMessages(library(survival))
d <- read.csv("scores.csv")
mice <- unique(d[, c("group", "mouse")])
per <- do.call(rbind, lapply(seq_len(nrow(mice)), function(i) {
  x <- d[d$mouse == mice$mouse[i] & d$group == mice$group[i], ]; on <- if (any(x$score >= 1)) min(x$day[x$score >= 1]) else NA
  data.frame(group = mice$group[i], mouse = mice$mouse[i], cumulative = sum(x$score), maximum = max(x$score),
             onset = on, diseased = as.integer(!is.na(on)), onset_time = ifelse(is.na(on), max(d$day), on)) }))
v <- per[per$group == "vehicle", ]; g <- per[per$group == "drug", ]
mw <- function(a, b) suppressWarnings(wilcox.test(b, a, exact = FALSE, correct = TRUE))$p.value
cat(sprintf("Incidence (score >= 1 on any day): vehicle %d/10, drug %d/10, Fisher p = %.4f\n", sum(v$diseased), sum(g$diseased),
            fisher.test(matrix(c(sum(v$diseased), 10 - sum(v$diseased), sum(g$diseased), 10 - sum(g$diseased)), 2, byrow = TRUE))$p.value))
cat(sprintf("Cumulative score (sum of daily scores): median vehicle %.1f, drug %.1f; Mann-Whitney p = %.4f\n",
            median(v$cumulative), median(g$cumulative), mw(v$cumulative, g$cumulative)))
q <- function(x) sprintf("%.1f to %.1f", quantile(x, 0.25), quantile(x, 0.75))
cat(sprintf("Cumulative score, interquartile range: vehicle %s, drug %s\n", q(v$cumulative), q(g$cumulative)))
cat(sprintf("Cumulative score, Hodges-Lehmann shift (drug - vehicle, median of all pairwise differences): %.1f\n", median(outer(g$cumulative, v$cumulative, "-"))))
cat(sprintf("Maximum score: median vehicle %.2f, drug %.2f; Mann-Whitney p = %.4f\n", median(v$maximum), median(g$maximum), mw(v$maximum, g$maximum)))
per$drug <- as.integer(per$group == "drug")
lr <- survdiff(Surv(onset_time, diseased) ~ drug, data = per)
cat(sprintf("Day of onset: median vehicle %.1f, drug %.1f (diseased mice); log-rank chi-square = %.2f, p = %.4f\n",
            median(v$onset, na.rm = TRUE), median(g$onset, na.rm = TRUE), lr$chisq, pchisq(lr$chisq, 1, lower.tail = FALSE)))
days <- sort(unique(d$day))
ps <- sapply(days, function(dd) { x <- d[d$day == dd, ]; a <- x$score[x$group == "vehicle"]; b <- x$score[x$group == "drug"]
  if (length(unique(c(a, b))) == 1) 1 else mw(a, b) })
cat(sprintf("Separate Mann-Whitney test on each of %d days: p < 0.05 on %d days (days %s)\n", length(days), sum(ps < 0.05),
            paste(days[ps < 0.05], collapse = ", ")))