demo/mrs_2d.R

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

library(MRS)

local({
  # This consolidates the former examples/example1.R workflow: fit a known
  # group contrast, compare it with shuffled labels, and inspect the fitted
  # multiresolution structure.
  set.seed(12345)

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

  n <- 220L
  x1 <- cbind(stats::rnorm(n), stats::rnorm(n))

  localized <- stats::rbinom(n, size = 1L, prob = 0.24) == 1L
  x2 <- cbind(stats::rnorm(n), stats::rnorm(n))
  x2[localized, 1] <- stats::rnorm(sum(localized), mean = 2.15, sd = 0.34)
  x2[localized, 2] <- stats::rnorm(sum(localized), mean = -1.75, sd = 0.34)

  X <- rbind(x1, x2)
  G <- rep(1:2, each = n)
  fit <- mrs(X, G, K = 5L)
  shuffled_fit <- mrs(X, sample(G), K = 5L)

  op <- graphics::par(mar = c(4.2, 4.3, 3.2, 1), las = 1)
  graphics::plot(
    X[, 1], X[, 2], type = "n", asp = 1,
    xlab = expression(X[1]), ylab = expression(X[2]),
    main = "A localized two-dimensional difference"
  )
  graphics::points(
    x1, pch = 1, cex = 0.85,
    col = grDevices::adjustcolor(blue, alpha.f = 0.58)
  )
  graphics::points(
    x2, pch = 16, cex = 0.65,
    col = grDevices::adjustcolor(orange, alpha.f = 0.48)
  )
  graphics::legend(
    "topleft", c("Reference", "Localized cluster"),
    pch = c(1, 16), col = c(blue, orange), pt.cex = c(1, 0.8), bty = "n"
  )
  graphics::text(
    2.15, -1.75, "localized excess", pos = 4, offset = 0.55,
    col = grey, cex = 0.82
  )
  graphics::arrows(
    2.02, -1.70, 1.72, -1.56, length = 0.08,
    col = grey, lwd = 1.2
  )
  graphics::mtext(
    sprintf("Posterior P(any difference) = %.3f", 1 - fit$PostGlobNull),
    side = 3, line = 0.35, adj = 1, col = grey, cex = 0.82
  )

  comparison <- c(
    `Observed labels` = 1 - fit$PostGlobNull,
    `Shuffled labels` = 1 - shuffled_fit$PostGlobNull
  )
  y <- rev(seq_along(comparison))
  graphics::par(mar = c(4.2, 8.2, 3.2, 1), las = 1)
  graphics::plot(
    NA, xlim = c(0, 1), ylim = c(0.5, length(comparison) + 0.5),
    xlab = "Posterior P(any difference)", ylab = "", yaxt = "n",
    main = "Observed signal versus a label-shuffled control"
  )
  graphics::abline(v = c(0.25, 0.5, 0.75), col = "#E5E5E5", lwd = 1)
  graphics::segments(0, y, comparison, y, col = c(orange, "#999999"), lwd = 4)
  graphics::points(comparison, y, pch = 21, cex = 1.8,
                   bg = c(orange, "#999999"), col = "white", lwd = 1.5)
  graphics::axis(2, at = y, labels = names(comparison), tick = FALSE, las = 1)
  graphics::text(
    comparison, y, labels = sprintf("%.3f", comparison),
    pos = ifelse(comparison > 0.90, 2, 4), offset = 0.65,
    xpd = NA, col = grey
  )
  graphics::par(op)

  plot2D(
    fit, type = "prob", levels = c(0L, 2L, 4L),
    data.points = "differential", background = "smeared", legend = TRUE,
    main = "Posterior probability of a difference"
  )

  plot2D(
    fit, type = "eff", group = 2L, abs = FALSE,
    levels = c(1L, 3L, 5L), data.points = "none", legend = TRUE,
    main = "Signed effect for the localized-cluster group"
  )

  plotTree(
    fit, type = "prob", main = "Multiresolution tree: evidence by node",
    legend = TRUE
  )

  cat("\nTwo-dimensional MRS demo\n")
  cat(sprintf("  Observations: %d (%d per group)\n", nrow(X), n))
  cat(sprintf("  Posterior P(any difference), observed labels: %.4f\n",
              1 - fit$PostGlobNull))
  cat(sprintf("  Posterior P(any difference), shuffled labels: %.4f\n",
              1 - shuffled_fit$PostGlobNull))
  cat("  Probability maps show how the signal resolves across tree levels.\n")

  # For global evidence without a tree, use return_tree = FALSE. Conversely,
  # return_global_null = FALSE omits the global-null calculation.
})