demo/andova_replicates.R

In MRS: Multi-Resolution Scanning for Cross-Sample and Cross-Group Differences

library(MRS)

local({
  # This consolidates the former andova_example.R and nested-model examples.
  # It compares pooled and replicate-aware fits at three levels of replication,
  # using the same four-view diagnostic sequence for every scenario.
  set.seed(12345)

  blue <- "#0072B2"
  orange <- "#D55E00"
  green <- "#009E73"
  grey <- "#4D4D4D"

  n_groups <- 2L
  replicate_counts <- c(Small = 4L, Medium = 10L, Large = 20L)
  max_replicates <- max(replicate_counts)
  n_per_replicate <- 40L
  X <- numeric()
  G <- integer()
  H <- integer()

  # Generate one 20-replicate study. The small and medium scenarios use its
  # first 4 and 10 replicates, so changes in the results reflect added
  # independent replication rather than three unrelated simulated data sets.
  for (h in seq_len(max_replicates)) {
    for (g in seq_len(n_groups)) {
      replicate_shift <- stats::rnorm(1L, sd = 0.22)
      is_tail <- stats::rbinom(n_per_replicate, size = 1L, prob = 0.22) == 1L
      values <- stats::rnorm(n_per_replicate, mean = replicate_shift, sd = 0.82)
      tail_location <- if (g == 1L) 2.35 else 2.90
      values[is_tail] <- stats::rnorm(
        sum(is_tail), mean = tail_location + replicate_shift, sd = 0.30
      )
      X <- c(X, values)
      G <- c(G, rep(g, n_per_replicate))
      H <- c(H, rep(h, n_per_replicate))
    }
  }
  X <- matrix(X, ncol = 1L)

  fit_scenario <- function(n_replicates) {
    keep <- H <= n_replicates
    scenario_X <- X[keep, , drop = FALSE]
    scenario_G <- G[keep]
    scenario_H <- H[keep]

    list(
      newton = andova(
        scenario_X, scenario_G, scenario_H, K = 5L,
        nu_vec = c(0.2, 1, 5), baseline = 1L,
        method = "newton", n_grid_theta = 4L
      ),
      riemann = andova(
        scenario_X, scenario_G, scenario_H, K = 5L,
        nu_vec = c(0.2, 1, 5), baseline = 1L,
        method = "riemann", n_grid_theta = 32L
      ),
      pooled = mrs(scenario_X, scenario_G, K = 5L),
      X = scenario_X,
      G = scenario_G,
      H = scenario_H
    )
  }
  scenario_fits <- lapply(unname(replicate_counts), fit_scenario)
  names(scenario_fits) <- names(replicate_counts)

  method_probabilities <- cbind(
    `Pooled MRS` = vapply(
      scenario_fits, function(x) 1 - x$pooled$PostGlobNull, numeric(1L)
    ),
    `ANDOVA: Newton` = vapply(
      scenario_fits, function(x) 1 - x$newton$PostGlobNull, numeric(1L)
    ),
    `ANDOVA: Riemann` = vapply(
      scenario_fits, function(x) 1 - x$riemann$PostGlobNull, numeric(1L)
    )
  )
  group_colors <- c(blue, orange)

  plot_scenario <- function(scenario_name, scenario_data, n_replicates) {
    grid <- seq(
      min(scenario_data$X), max(scenario_data$X), length.out = 300L
    )
    density_values <- array(
      NA_real_, dim = c(length(grid), n_replicates, n_groups)
    )
    for (g in seq_len(n_groups)) {
      for (h in seq_len(n_replicates)) {
        values <- scenario_data$X[
          scenario_data$G == g & scenario_data$H == h, 1
        ]
        density_values[, h, g] <- stats::density(
          values, from = min(grid), to = max(grid), n = length(grid)
        )$y
      }
    }

    density_par <- graphics::par(no.readonly = TRUE)
    graphics::par(mar = c(4.2, 4.3, 3.2, 1), las = 1)
    graphics::matplot(
      grid, density_values[, , 1], type = "n",
      ylim = c(0, max(density_values) * 1.08),
      xlab = "Observed value", ylab = "Density",
      main = sprintf(
        "%s study: replicate and group-average structure", scenario_name
      )
    )
    for (g in seq_len(n_groups)) {
      graphics::matlines(
        grid, density_values[, , g], lty = 1, lwd = 1.2,
        col = grDevices::adjustcolor(group_colors[g], alpha.f = 0.30)
      )
      graphics::lines(
        grid, rowMeans(density_values[, , g]),
        col = group_colors[g], lwd = 3
      )
    }
    graphics::legend(
      "topleft",
      c("Reference mean", "Shifted-tail mean", "Individual replicates"),
      col = c(blue, orange, grey), lwd = c(3, 3, 1.2),
      lty = 1, bty = "n"
    )
    graphics::mtext(
      sprintf(
        "%d thin curves per group show the variation retained by ANDOVA",
        n_replicates
      ),
      side = 3, line = 0.35, adj = 1, col = grey, cex = 0.80
    )
    graphics::par(density_par)

    plot1D(
      scenario_data$newton, type = "prob", legend = TRUE,
      main = sprintf("%s study: ANDOVA probability map", scenario_name)
    )

    plot1D(
      scenario_data$newton, type = "eff", group = 2L,
      abs = FALSE, legend = TRUE, eff_scale = "odds-ratio",
      main = sprintf(
        "%s study: shifted-tail odds ratio versus reference", scenario_name
      )
    )

    plotTree(
      scenario_data$newton, type = "prob", legend = TRUE,
      main = sprintf(
        "%s study: ANDOVA multiresolution tree", scenario_name
      )
    )
  }

  for (scenario_name in names(scenario_fits)) {
    plot_scenario(
      scenario_name,
      scenario_fits[[scenario_name]],
      unname(replicate_counts[scenario_name])
    )
  }

  cat("\nReplicate-count comparison\n")
  cat("  Scenario   Replicates   Pooled MRS   ANDOVA/Newton   ANDOVA/Riemann\n")
  for (scenario in seq_along(replicate_counts)) {
    cat(sprintf(
      "  %-8s %10d %12.4f %16.4f %17.4f\n",
      names(replicate_counts)[scenario], replicate_counts[scenario],
      method_probabilities[scenario, "Pooled MRS"],
      method_probabilities[scenario, "ANDOVA: Newton"],
      method_probabilities[scenario, "ANDOVA: Riemann"]
    ))
  }
  cat("  Interpretation: pooled MRS is overconfident for the small and medium\n")
  cat("  studies; ANDOVA uncertainty decreases with independent replication.\n")

  # Pooled MRS discards H and treats every observation as independent evidence.
  # With only 4--10 genuine replicates this is pseudoreplication: the pooled
  # posterior is already near one and therefore overconfident. ANDOVA treats
  # replicates as the sampling units. It retains substantial uncertainty with
  # four replicates, then becomes progressively more certain as the number of
  # independent replicates increases to 10 and 20.
  method_colors <- c("#999999", green, blue)
  method_pch <- c(15, 19, 17)
  method_lty <- c(2, 1, 3)
  comparison_par <- graphics::par(no.readonly = TRUE)
  graphics::par(mar = c(5.2, 4.5, 3.6, 1), las = 1)
  graphics::plot(
    NA, xlim = c(3, 21), ylim = c(0, 1.04), xaxt = "n",
    xlab = "Independent replicates per group",
    ylab = "Posterior P(any cross-group difference)",
    main = "Replication reveals pooled-MRS overconfidence"
  )
  graphics::axis(
    1, at = unname(replicate_counts),
    labels = sprintf("%s\n(%d)", names(replicate_counts), replicate_counts)
  )
  graphics::abline(h = c(0.25, 0.5, 0.75, 1), col = "#E5E5E5", lwd = 1)
  for (method in seq_len(ncol(method_probabilities))) {
    graphics::lines(
      unname(replicate_counts), method_probabilities[, method],
      type = "b", lwd = 2.4, lty = method_lty[method],
      pch = method_pch[method], cex = 1.1, col = method_colors[method]
    )
  }
  graphics::legend(
    "bottomright", colnames(method_probabilities),
    col = method_colors, lty = method_lty, pch = method_pch,
    lwd = 2.4, pt.cex = 1.1, bty = "n"
  )
  graphics::text(
    7.1, 0.91,
    "Pooled MRS: premature confidence\n(replicate structure is ignored)",
    cex = 0.78, col = grey
  )
  graphics::arrows(
    5.25, 0.95, 4.25, method_probabilities[1L, "Pooled MRS"] - 0.01,
    length = 0.07, col = grey, lwd = 1.1
  )
  graphics::text(
    12.8, 0.34,
    "ANDOVA uncertainty decreases\nas independent replicates accumulate",
    cex = 0.78, col = green
  )
  graphics::arrows(
    15.0, 0.43, 19.1, 0.90,
    length = 0.07, col = green, lwd = 1.1
  )
  graphics::mtext(
    "The 4-replicate data are nested within the 10- and 20-replicate studies",
    side = 3, line = 0.30, adj = 1, cex = 0.76, col = grey
  )
  graphics::par(comparison_par)
})