demo/mrs_1d.R

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

library(MRS)

local({
  # A compact one-dimensional workflow combining data inspection,
  # multiresolution evidence, signed effects, and Bayesian FDR selection.
  set.seed(12345)

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

  n <- 220L
  make_group <- function(tail_location) {
    is_tail <- stats::rbinom(n, size = 1L, prob = 0.24) == 1L
    x <- stats::rnorm(n, mean = 0, sd = 0.90)
    x[is_tail] <- stats::rnorm(sum(is_tail), mean = tail_location, sd = 0.28)
    sample(x)
  }

  x1 <- make_group(2.55)
  x2 <- make_group(3.25)
  X <- matrix(c(x1, x2), ncol = 1L)
  G <- rep(1:2, each = n)

  fit <- mrs(X, G, K = 6L)

  densities <- lapply(list(x1, x2), stats::density, from = min(X), to = max(X))
  ymax <- max(vapply(densities, function(z) max(z$y), numeric(1)))
  op <- graphics::par(mar = c(4.2, 4.3, 3.2, 1), las = 1)
  graphics::plot(
    densities[[1]], type = "n", ylim = c(0, ymax * 1.08),
    xlab = "Observed value", ylab = "Density",
    main = "A localized distributional change"
  )
  graphics::polygon(
    c(densities[[1]]$x, rev(densities[[1]]$x)),
    c(densities[[1]]$y, rep(0, length(densities[[1]]$y))),
    col = grDevices::adjustcolor(blue, alpha.f = 0.20), border = NA
  )
  graphics::polygon(
    c(densities[[2]]$x, rev(densities[[2]]$x)),
    c(densities[[2]]$y, rep(0, length(densities[[2]]$y))),
    col = grDevices::adjustcolor(orange, alpha.f = 0.20), border = NA
  )
  graphics::lines(densities[[1]], col = blue, lwd = 2.5)
  graphics::lines(densities[[2]], col = orange, lwd = 2.5)
  graphics::rug(x1, side = 1, col = grDevices::adjustcolor(blue, alpha.f = 0.22))
  graphics::rug(x2, side = 3, col = grDevices::adjustcolor(orange, alpha.f = 0.22))
  graphics::legend(
    "topleft", c("Reference", "Shifted tail"), col = c(blue, orange),
    lwd = 3, bty = "n"
  )
  graphics::mtext(
    sprintf("Posterior P(any difference) = %.3f", 1 - fit$PostGlobNull),
    side = 3, line = 0.35, adj = 1, col = grey, cex = 0.82
  )
  graphics::par(op)

  plot1D(
    fit, type = "prob", legend = TRUE,
    main = "Posterior probability map: where do the groups differ?"
  )

  plot1D(
    fit, type = "eff", group = 2L, abs = FALSE, legend = TRUE,
    main = "Signed effect map for the shifted-tail group"
  )

  selected <- plot1DSigWindows(
    fit, fdr = 0.10, type = "prob", legend = TRUE,
    main = "Windows retained at 10% Bayesian FDR"
  )

  cat("\nOne-dimensional MRS demo\n")
  cat(sprintf("  Observations: %d (%d per group)\n", nrow(X), n))
  cat(sprintf("  Posterior P(any difference): %.4f\n", 1 - fit$PostGlobNull))
  cat(sprintf("  FDR-selected tree nodes: %d\n", length(selected$indices)))
  cat(sprintf("  Selection threshold: %.4f\n", selected$threshold))
})