# =====================================================================
# CV-Ratio Screen - interactive companion app (R / Shiny)
# =====================================================================
# Companion software for:
#   Dwyer WJ. "The Coefficient of Variation as a Cutoff for Detecting
#   Non-Normal Distributions from Summary Statistics."
#
# Run with:
#   install.packages(c("shiny", "ggplot2"))
#   shiny::runApp("app.R")      # or click "Run App" in RStudio
#
# Tabs:
#   1. Reproduce the findings   - Monte-Carlo simulator
#   2. Analyze your data        - two groups -> recommended test
#   3. Summary-statistic screen - published mean & SD -> log-normality flag
# =====================================================================

library(shiny)
library(ggplot2)

# ---------------------------------------------------------------------
# CORE STATISTICS
# ---------------------------------------------------------------------
lognormal_params <- function(mu_n, cv) {
  s2 <- log(cv^2 + 1)
  c(mu_e = log(mu_n) - s2 / 2, sigma_e = sqrt(s2))
}

draw_lognormal <- function(mu_n, cv, n, reps, logshift = 0) {
  p <- lognormal_params(mu_n, cv)
  matrix(rlnorm(n * reps, meanlog = p["mu_e"] + logshift, sdlog = p["sigma_e"]),
         nrow = reps, ncol = n)
}

permutation_t_p <- function(a, b, B = 2000, use_log = TRUE) {
  if (use_log) { a <- log(a); b <- log(b) }
  obs <- abs(mean(a) - mean(b)); pool <- c(a, b); n1 <- length(a)
  cnt <- sum(replicate(B, {
    p <- sample(pool); abs(mean(p[1:n1]) - mean(p[(n1 + 1):length(pool)])) >= obs - 1e-12
  }))
  (cnt + 1) / (B + 1)
}

cv_ratio_screen <- function(mu1, sd1, mu2, sd2) {
  cv1 <- sd1 / mu1; cv2 <- sd2 / mu2
  cv_hi <- max(cv1, cv2); cv_lo <- min(cv1, cv2)
  ratio <- if (cv_lo > 0) cv_hi / cv_lo else Inf
  larger_mean_has_larger_var <- (sd1 >= sd2) == (mu1 >= mu2)
  suspect <- (cv_hi >= 0.5) || (ratio > 1.5 && larger_mean_has_larger_var)
  verdict <- if (cv_hi < 0.3) "Effectively normal - raw-scale t-test is appropriate."
    else if (cv_hi < 0.5) "Borderline - log transform begins to matter; check the raw data if possible."
    else if (cv_hi < 1.0) "Suspect non-normality - prefer a log-scale analysis."
    else "Strongly suspect log-normality - a raw-scale t-test is likely deficient; use the log scale."
  list(cv1 = cv1, cv2 = cv2, cv_ratio = ratio, suspect = suspect, verdict = verdict)
}

analyze_two_groups <- function(x, y, n_perm = 2000) {
  positive <- all(x > 0) && all(y > 0)
  out <- list(
    n1 = length(x), n2 = length(y),
    mean1 = mean(x), sd1 = sd(x), mean2 = mean(y), sd2 = sd(y),
    cv1 = sd(x) / mean(x), cv2 = sd(y) / mean(y), positive = positive,
    sw_raw1 = if (length(x) >= 3) shapiro.test(x)$p.value else NA,
    sw_raw2 = if (length(y) >= 3) shapiro.test(y)$p.value else NA,
    p_raw_t = t.test(x, y)$p.value,
    p_mw = suppressWarnings(wilcox.test(x, y, exact = FALSE)$p.value)
  )
  if (positive) {
    out$sw_log1 <- if (length(x) >= 3) shapiro.test(log(x))$p.value else NA
    out$sw_log2 <- if (length(y) >= 3) shapiro.test(log(y))$p.value else NA
    out$p_log_t <- t.test(log(x), log(y))$p.value
    out$p_perm_log <- permutation_t_p(x, y, B = n_perm, use_log = TRUE)
  }
  scr <- cv_ratio_screen(out$mean1, out$sd1, out$mean2, out$sd2)
  out$screen <- scr
  out$recommended <- if (!positive)
      "Data contain non-positive values: log transform not applicable. Use the raw t-test if approximately normal, else Mann-Whitney or a permutation test."
    else if (scr$suspect)
      "Log-scale t-test (most powerful under log-normality); report the permutation t-test on the log scale as confirmation."
    else "Raw-scale t-test is appropriate."
  out
}

NS_GRID <- c(10, 15, 20, 30, 50, 100, 200)

# ---------------------------------------------------------------------
# UI
# ---------------------------------------------------------------------
ui <- fluidPage(
  titlePanel("CV-Ratio Screen"),
  helpText("Detecting non-normal (log-normal) distributions from summary statistics, and recovering parametric inference."),
  tabsetPanel(
    # ---- Tab 1: simulate ----
    tabPanel("1. Reproduce the findings",
      sidebarLayout(
        sidebarPanel(
          selectInput("finding", "Finding",
            c("Detectability: normality-test power vs n",
              "Inference cost: power raw vs log vs Mann-Whitney",
              "The cutoff: CV where log starts to matter")),
          sliderInput("reps", "Monte-Carlo replicates", 500, 6000, 1500, step = 500),
          numericInput("seed", "Random seed", 2026),
          conditionalPanel("input.finding.indexOf('Detectability') == 0 || input.finding.indexOf('Inference') == 0",
            sliderInput("cv", "Coefficient of variation", 0.2, 3.0, 1.0, step = 0.1)),
          conditionalPanel("input.finding.indexOf('Inference') == 0",
            sliderInput("ratio", "True median ratio (effect)", 1.0, 3.0, 1.5, step = 0.1)),
          conditionalPanel("input.finding.indexOf('cutoff') > -1",
            sliderInput("n", "Sample size per group", 10, 100, 30, step = 5)),
          actionButton("run", "Run simulation", class = "btn-primary")
        ),
        mainPanel(plotOutput("simPlot", height = "430px"), textOutput("simNote"))
      )
    ),
    # ---- Tab 2: analyze ----
    tabPanel("2. Analyze your data",
      sidebarLayout(
        sidebarPanel(
          textAreaInput("g1", "Group 1 (control)", "0.4 0.5 0.6 0.3 1.2 0.5 0.7 0.4 0.6 0.5", rows = 4),
          textAreaInput("g2", "Group 2 (disease)", "1.1 2.3 0.9 4.5 1.4 1.0 6.2 1.3 2.0 1.5", rows = 4),
          sliderInput("nperm", "Permutation replicates", 500, 8000, 2000, step = 500),
          actionButton("analyze", "Analyze", class = "btn-primary")
        ),
        mainPanel(
          verbatimTextOutput("analyzeOut"),
          plotOutput("histPlot", height = "320px")
        )
      )
    ),
    # ---- Tab 3: summary screen ----
    tabPanel("3. Summary-statistic screen",
      sidebarLayout(
        sidebarPanel(
          numericInput("mu1", "Group 1 mean", 0.5),
          numericInput("sd1", "Group 1 SD", 0.5),
          numericInput("mu2", "Group 2 mean", 1.5),
          numericInput("sd2", "Group 2 SD", 3.4)
        ),
        mainPanel(verbatimTextOutput("screenOut"))
      )
    )
  ),
  hr(),
  helpText("Methods: CV-ratio screen (eq 17), log-scale and permutation t-tests, Mann-Whitney U.")
)

# ---------------------------------------------------------------------
# SERVER
# ---------------------------------------------------------------------
server <- function(input, output, session) {

  simData <- eventReactive(input$run, {
    set.seed(input$seed)
    f <- input$finding
    if (startsWith(f, "Detectability")) {
      df <- data.frame()
      for (n in NS_GRID) {
        X <- draw_lognormal(1, input$cv, n, input$reps)
        sw <- mean(apply(X, 1, function(r) shapiro.test(r)$p.value) < .05)
        ad <- mean(apply(X, 1, function(r) {
          a <- nortest_ad(r); a }) < .05)
        df <- rbind(df,
          data.frame(n = n, test = "Shapiro-Wilk", power = sw),
          data.frame(n = n, test = "Anderson-Darling", power = ad))
      }
      list(kind = "detect", df = df)
    } else if (startsWith(f, "Inference")) {
      shift <- log(input$ratio); df <- data.frame()
      for (n in NS_GRID) {
        g1 <- draw_lognormal(1, input$cv, n, input$reps)
        g2 <- draw_lognormal(1, input$cv, n, input$reps, logshift = shift)
        raw <- mean(sapply(1:input$reps, function(i) t.test(g1[i, ], g2[i, ])$p.value) < .05)
        lg  <- mean(sapply(1:input$reps, function(i) t.test(log(g1[i, ]), log(g2[i, ]))$p.value) < .05)
        mw  <- mean(sapply(1:input$reps, function(i) suppressWarnings(wilcox.test(g1[i, ], g2[i, ], exact = FALSE)$p.value)) < .05)
        df <- rbind(df,
          data.frame(n = n, test = "log-scale t", power = lg),
          data.frame(n = n, test = "Mann-Whitney", power = mw),
          data.frame(n = n, test = "raw-scale t", power = raw))
      }
      list(kind = "power", df = df)
    } else {
      cvg <- c(0.2, 0.3, 0.5, 0.75, 1, 1.5, 2, 2.5, 3); df <- data.frame()
      for (cc in cvg) {
        go <- draw_lognormal(1, cc, input$n, input$reps)
        gd <- draw_lognormal(1.5, cc, input$n, input$reps)
        pr <- mean(sapply(1:input$reps, function(i) t.test(go[i, ], gd[i, ])$p.value) < .05)
        pl <- mean(sapply(1:input$reps, function(i) t.test(log(go[i, ]), log(gd[i, ]))$p.value) < .05)
        df <- rbind(df, data.frame(cv = cc, gain = pl - pr))
      }
      list(kind = "cutoff", df = df)
    }
  })

  output$simPlot <- renderPlot({
    d <- simData()
    if (d$kind %in% c("detect", "power")) {
      ggplot(d$df, aes(n, power, colour = test)) +
        geom_line() + geom_point() + scale_x_log10(breaks = NS_GRID) +
        geom_hline(yintercept = 0.8, linetype = 2, colour = "grey50") +
        labs(x = "sample size per group, n", y = "power", colour = NULL) +
        theme_minimal(base_size = 14)
    } else {
      ggplot(d$df, aes(cv, gain)) +
        annotate("rect", xmin = 0.5, xmax = 1.0, ymin = -Inf, ymax = Inf, alpha = .15, fill = "gold") +
        geom_line(colour = "#c44e52") + geom_point(colour = "#c44e52") +
        geom_hline(yintercept = 0.05, linetype = 2, colour = "grey50") +
        labs(x = "CV (each group)", y = "power gained by log transform") +
        theme_minimal(base_size = 14)
    }
  })
  output$simNote <- renderText({
    d <- simData()
    switch(d$kind,
      detect = "All tests lose power at small n - the blind spot is shared, not a Shapiro-Wilk artifact.",
      power  = "The log-scale t-test recovers power the raw test forfeits and exceeds the non-parametric test.",
      cutoff = "Below CV~0.3 data are effectively normal; by CV~1.0 the raw t-test is clearly deficient.")
  })

  output$analyzeOut <- renderPrint({
    input$analyze
    isolate({
      x <- as.numeric(strsplit(gsub(",", " ", input$g1), "\\s+")[[1]])
      y <- as.numeric(strsplit(gsub(",", " ", input$g2), "\\s+")[[1]])
      x <- x[!is.na(x)]; y <- y[!is.na(y)]
      if (length(x) < 3 || length(y) < 3) { cat("Need >= 3 observations per group."); return(invisible()) }
      R <- analyze_two_groups(x, y, n_perm = input$nperm)
      cat(sprintf("Group 1: n=%d  mean=%.3g  SD=%.3g  CV=%.3g\n", R$n1, R$mean1, R$sd1, R$cv1))
      cat(sprintf("Group 2: n=%d  mean=%.3g  SD=%.3g  CV=%.3g\n", R$n2, R$mean2, R$sd2, R$cv2))
      cat(sprintf("CV ratio (eq 17): %.2f\n", R$screen$cv_ratio))
      cat("Screen verdict: ", R$screen$verdict, "\n\n", sep = "")
      cat("Two-sided p-values:\n")
      cat(sprintf("  raw-scale t (Welch):  %.4f\n", R$p_raw_t))
      if (!is.null(R$p_log_t))   cat(sprintf("  log-scale t (Welch):  %.4f\n", R$p_log_t))
      if (!is.null(R$p_perm_log)) cat(sprintf("  permutation t (log):  %.4f\n", R$p_perm_log))
      cat(sprintf("  Mann-Whitney U:       %.4f\n\n", R$p_mw))
      cat(">> Recommended: ", R$recommended, "\n", sep = "")
    })
  })

  output$histPlot <- renderPlot({
    input$analyze
    isolate({
      x <- as.numeric(strsplit(gsub(",", " ", input$g1), "\\s+")[[1]])
      y <- as.numeric(strsplit(gsub(",", " ", input$g2), "\\s+")[[1]])
      x <- x[!is.na(x) & x > 0]; y <- y[!is.na(y) & y > 0]
      if (length(x) < 1 || length(y) < 1) return(NULL)
      df <- rbind(
        data.frame(value = c(x, y), grp = rep(c("G1", "G2"), c(length(x), length(y))), scale = "raw"),
        data.frame(value = log(c(x, y)), grp = rep(c("G1", "G2"), c(length(x), length(y))), scale = "log"))
      ggplot(df, aes(value, fill = grp)) + geom_histogram(alpha = .6, bins = 12, position = "identity") +
        facet_wrap(~scale, scales = "free") + theme_minimal(base_size = 13) + labs(x = NULL, fill = NULL)
    })
  })

  output$screenOut <- renderPrint({
    if (input$mu1 <= 0 || input$mu2 <= 0) { cat("Means must be positive."); return(invisible()) }
    scr <- cv_ratio_screen(input$mu1, input$sd1, input$mu2, input$sd2)
    p1 <- lognormal_params(input$mu1, input$sd1 / input$mu1)
    p2 <- lognormal_params(input$mu2, input$sd2 / input$mu2)
    cat(sprintf("CV  Group 1: %.2f\nCV  Group 2: %.2f\nCV ratio:    %.2f\n\n", scr$cv1, scr$cv2, scr$cv_ratio))
    cat(sprintf("Implied log-domain parameters:\n  G1: mu_e=%.3f  sigma_e=%.3f\n  G2: mu_e=%.3f  sigma_e=%.3f\n\n",
                p1["mu_e"], p1["sigma_e"], p2["mu_e"], p2["sigma_e"]))
    cat("Verdict: ", scr$verdict, "\n", sep = "")
  })
}

# Anderson-Darling (self-contained; avoids requiring the 'nortest' package)
nortest_ad <- function(x) {
  x <- sort(x); n <- length(x)
  p <- pnorm((x - mean(x)) / sd(x))
  p <- pmin(pmax(p, 1e-12), 1 - 1e-12)
  S <- sum((2 * seq_len(n) - 1) * (log(p) + log(1 - rev(p))))
  A2 <- -n - S / n
  A2 <- A2 * (1 + 0.75 / n + 2.25 / n^2)      # small-sample correction
  # convert to approximate p-value (Stephens); compare to ~0.752 critical at 5%
  if (A2 < 0.2) p_val <- 1 - exp(-13.436 + 101.14 * A2 - 223.73 * A2^2)
  else if (A2 < 0.34) p_val <- 1 - exp(-8.318 + 42.796 * A2 - 59.938 * A2^2)
  else if (A2 < 0.6) p_val <- exp(0.9177 - 4.279 * A2 - 1.38 * A2^2)
  else p_val <- exp(1.2937 - 5.709 * A2 + 0.0186 * A2^2)
  p_val
}

shinyApp(ui, server)
