## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
                      fig.width = 7, fig.height = 3.6)

## ----flow, echo=FALSE, fig.width=7.6, fig.height=4.8, fig.cap="The order matters: each step only makes sense if the one above it passed."----
# Drawn with ggplot2 rather than a diagram package: ggplot2 is already a hard
# dependency, so the figure cannot quietly vanish on a machine where an optional
# one is missing.
pal <- c("#0072B2", "#D55E00", "#009E73")
BOX <- function(x, y, w, label, fill, size = 3.0, h = 1.25)
  list(ggplot2::annotate("tile", x = x, y = y, width = w, height = h, fill = fill),
       ggplot2::annotate("text", x = x, y = y, label = label, colour = "white",
                         size = size, lineheight = 1.0))
AR <- function(x0, y0, x1, y1, lab = NULL)
  c(list(ggplot2::annotate("segment", x = x0, y = y0, xend = x1, yend = y1,
           arrow = grid::arrow(length = grid::unit(0.15, "cm"), type = "closed"),
           colour = "grey45", linewidth = 0.4)),
    if (!is.null(lab)) list(ggplot2::annotate("text", x = (x0 + x1)/2 - 0.30,
           y = (y0 + y1)/2, label = lab, size = 2.7, colour = "grey30", hjust = 1)))
sp <- 2.0
ggplot2::ggplot() +
  BOX(sp, 9.0, 5.4, "1.  Is the subspace stable?\nmgccaStability()", pal[1]) +
  AR(sp - 1.4, 8.35, -4.0, 7.75, "no") + AR(sp, 8.35, sp, 7.75, "yes") +
  BOX(-4.0, 7.1, 4.0, "STOP\nnothing below\nis interpretable", "grey55", 2.8) +
  BOX(sp, 7.1, 5.4, "2.  Is T above the\nnoise floor?", pal[1]) +
  AR(sp - 1.4, 6.45, -4.0, 5.85, "no") + AR(sp, 6.45, sp, 5.85, "yes") +
  BOX(-4.0, 5.2, 4.0, "STOP\nthe omics have\nnothing to say", "grey55", 2.8) +
  BOX(sp, 5.2, 5.4, "3.  Is R_between\nhigh?", pal[1]) +
  AR(sp - 1.4, 4.55, -1.0, 3.95, "no") + AR(sp + 1.4, 4.55, 6.4, 3.95, "yes") +
  BOX(-1.0, 3.3, 5.0, "STRONGEST\na real share, and it does\nnot lean on your groups", pal[3], 2.8) +
  BOX(6.4, 3.3, 5.0, "CHECK\nwhat IS your grouping?", pal[2], 2.8) +
  AR(5.2, 2.65, 4.0, 2.05) + AR(7.6, 2.65, 8.8, 2.05) +
  BOX(4.0, 1.4, 4.4, "biological\n(stage, tissue)\n= a finding", pal[3], 2.7) +
  BOX(8.8, 1.4, 4.4, "technical\n(plate, batch)\n= suspect an artefact", pal[2], 2.7) +
  ggplot2::coord_cartesian(xlim = c(-6.3, 11.3), ylim = c(0.5, 9.9)) +
  ggplot2::theme_void() +
  ggplot2::labs(title = "In what order to read the numbers",
                subtitle = "each step only makes sense if the one above it passed") +
  ggplot2::theme(plot.title = ggplot2::element_text(face = "bold", size = 11),
                 plot.subtitle = ggplot2::element_text(size = 9, colour = "grey35"))

## ----install, eval=FALSE------------------------------------------------------
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("mgcca")

## ----load_library-------------------------------------------------------------
library(mgcca)

## ----cardio_fit---------------------------------------------------------------
data(cardiovascular)
X <- list(methylation = as.matrix(X1),
          clinical    = as.matrix(X2),
          other       = as.matrix(X3))
X <- lapply(X, function(m) { storage.mode(m) <- "double"; m })

h5file <- tempfile(fileext = ".h5")
fit <- mgcca(X, filename = h5file, nfac = 2, method = "solve", scores = TRUE)
fit

## ----cardio_query-------------------------------------------------------------
# The clinical table covers 633 of the 646 individuals in the fit. A query does
# NOT have to cover everyone: those it misses simply contribute nothing, and
# `n_used` reports how many actually did.
cids <- intersect(rownames(fit$Y), rownames(X2))

# The noise floor is ESTIMATED, not taken from one draw. A single random column is
# itself a random number: it can land high or low, and reading a real variable
# against one lucky draw is how a null result becomes a finding. Five columns give
# the floor a range instead of a point.
set.seed(2026)
noise <- replicate(5, rnorm(length(cids)))
colnames(noise) <- paste0("noise", 1:5)

Q <- cbind(LDL     = X2[cids, "LDL"],
           glucose = X2[cids, "Gluc"],
           noise)
rownames(Q) <- cids

sens <- mgccaSensitivity(fit, Q)
sens

## ----cardio_read, echo=FALSE, results='asis'----------------------------------
o     <- sens$overall
tl    <- o[["T"]][o$query == "LDL"]
tg    <- o[["T"]][o$query == "glucose"]
floor <- o[["T"]][grepl("^noise", o$query)]
hi    <- max(floor)
say <- function(x, nm) {
  if (x > 3 * hi)
    sprintf(paste0("**%s clears the ceiling comfortably.** The omic signal accounts ",
                   "for %.2f%% of it, against a ceiling of %.2f%%, so what it picks ",
                   "up is not something random numbers would have produced."),
            nm, 100 * x, 100 * hi)
  else if (x > hi)
    sprintf(paste0("**%s is above the ceiling, but only just** -- %.2f%% against ",
                   "%.2f%%. A gap that narrow is worth noting and not worth ",
                   "building on: treat it as suggestive, and see whether it holds ",
                   "in a second cohort before it goes in an abstract."),
            nm, 100 * x, 100 * hi)
  else
    sprintf(paste0("**%s does not clear the ceiling at all** (%.2f%% against ",
                   "%.2f%%). The integrated signal accounts for no more of it than ",
                   "it would of random numbers, so there is nothing here to explain ",
                   "and no follow-up analysis will change that."),
            nm, 100 * x, 100 * hi)
}
cat(sprintf(paste0(
  "**Reading this table.** The shared axes account for %.2f%% of the variation in ",
  "LDL and %.2f%% of glucose.\n\n",
  "Now, is that a lot or a little? On its own, neither -- the number has no ",
  "conventional threshold. That is what the five `noise` rows are for. They are ",
  "columns of random numbers put through exactly the same machinery, and they come ",
  "back between %.2f%% and %.2f%%. **A variable with no relationship whatsoever ",
  "still scores in that band**, simply because projecting anything onto a couple ",
  "of axes catches a little by chance. So that band is the ceiling that ",
  "meaningless variables reach in this analysis, and it is the yardstick.\n\n",
  "(Five columns rather than one, because a single random column is itself a ",
  "random draw: it might land at %.2f%% or at %.2f%%. Reading a real variable ",
  "against one lucky draw is exactly how a null result turns into a finding.)\n\n",
  "%s\n\n%s\n\n",
  "What this does **not** say: nothing about methylation values being higher or ",
  "lower, nothing about direction, and nothing about causing or being caused. It ",
  "says how much of a variable's variation the integrated signal can reproduce, ",
  "and whether that is more than chance alone would give.\n"),
  100 * tl, 100 * tg, 100 * min(floor), 100 * hi,
  100 * min(floor), 100 * hi,
  say(tl, "LDL"), say(tg, "Glucose")))

## ----cardio_floor, fig.height=3.6, fig.cap="Only a query whose point sits outside the grey band is saying anything the machinery would not have produced from random numbers."----
plot(sens, type = "overall", noise = paste0("noise", 1:5))

## ----cardio_profile-----------------------------------------------------------
prof <- mgccaBlockProfile(sens)
prof

## ----cardio_profile_plot, fig.height=4.4--------------------------------------
plot(sens)

## ----cardio_profile_read, echo=FALSE, results='asis'--------------------------
bb <- prof$profile[prof$profile$query == "LDL", ]
top <- bb$block[which.max(bb$share_total)]
cat(sprintf(
  "Here the largest contributor to LDL is **%s** (%.0f%% of its total sensitivity).\n",
  top, 100 * max(bb$share_total)))

## ----acc_build----------------------------------------------------------------
library(MultiAssayExperiment)
library(SummarizedExperiment)
data("miniACC", package = "MultiAssayExperiment")

# Assays are features x samples; mgcca wants individuals x variables with the
# patient ID as row name, so each assay is renamed through the sample map and
# transposed. A few hundred variables per block keeps the vignette quick.
sm  <- sampleMap(miniACC)
use <- c(mRNA = "RNASeq2GeneNorm", protein = "RPPAArray", miRNA = "miRNASeqGene")

tabs <- lapply(use, function(a) {
  m <- assay(experiments(miniACC)[[a]])
  k <- sm[sm$assay == a, ]
  colnames(m) <- as.character(k$primary[match(colnames(m), as.character(k$colname))])
  m <- m[, !duplicated(colnames(m)), drop = FALSE]
  Xa <- t(m); storage.mode(Xa) <- "double"
  keep <- apply(Xa, 2, function(v) all(is.finite(v)) && sd(v) > 0)
  Xa[, which(keep)[seq_len(min(sum(keep), 150))], drop = FALSE]
})

sapply(tabs, dim)
length(Reduce(union, lapply(tabs, rownames)))       # union of individuals
length(Reduce(intersect, lapply(tabs, rownames)))   # present in all three

## ----acc_lambda, message=TRUE-------------------------------------------------
sel <- mgccaSelectLambda(tabs, nfac = 2, K = 5, seed = 1)
sel

## ----acc_fit------------------------------------------------------------------
h5acc <- tempfile(fileext = ".h5")
acc <- mgcca(tabs, filename = h5acc, nfac = 2,
             method = "penalized", lambda = sel$lambda, scores = TRUE)

ids  <- rownames(acc$Y)
cd   <- as.data.frame(colData(miniACC))
stage <- factor(as.character(cd[ids, "pathologic_stage"]))
names(stage) <- ids
table(stage, useNA = "ifany")

## ----acc_sens-----------------------------------------------------------------
score <- setNames(as.numeric(stage), ids)
# FIVE noise columns, not one. A single random column is a DRAW from the noise
# level, not the noise level -- and on the map below, one draw lands wherever it
# lands, sometimes squarely along the direction the groups separate in. Five of
# them show the spread, which is the honest comparison.
set.seed(2026)
Q <- cbind(stage_score = score, replicate(5, rnorm(length(ids))))
colnames(Q) <- c("stage_score", paste0("noise", 1:5))
rownames(Q) <- ids

sens_acc <- mgccaSensitivity(acc, Q, group = stage)
sens_acc

## ----acc_sens_read, echo=FALSE, results='asis'--------------------------------
o  <- sens_acc$overall
sg <- o[o$query == "stage_score", ]; nzr <- range(o[["T"]][grepl("^noise", o$query)])
r  <- sg$R_between
band <- if (r < 0.05) "very low" else if (r < 0.20) "low" else if (r < 0.50) "moderate" else "high"
cat(sprintf(paste0(
  "**Reading this table.** The three assays together account for %.2f%% of the ",
  "variation in tumour stage, against %.2f--%.2f%% for the five random columns. ",
  "Of everything that could shift that reading, **%.1f%% follows stage itself as ",
  "a grouping** (`R_between` = %.3f) -- a %s value.\n\n",
  "**Caution:** note where the alignment sits relative to the noise columns: %s. ",
  "With this many patients and two retained axes, an arbitrary vector picks up a ",
  "share of the subspace by chance, and `T` alone does not distinguish the two. ",
  "What does is the grouping decomposition below, and the direction on the map.\n"),
  100 * sg[["T"]], 100 * nzr[1], 100 * nzr[2], 100 * r, r, band,
  if (sg[["T"]] > nzr[2]) "above every one of them"
  else if (sg[["T"]] < nzr[1]) "**below every one of them**"
  else "**inside their range**"))

## ----acc_plot, fig.height=4.6-------------------------------------------------
plot(sens_acc)

## ----acc_map, fig.height=4.8, fig.cap="Individuals on the reference axes, group centroids as diamonds, and each query as an arrow."----
plot(sens_acc, type = "map")

## ----acc_map_angles-----------------------------------------------------------
R2  <- sens_acc$reference
Y2  <- R2$Y[, 1:2, drop = FALSE]
gv  <- R2$group[rownames(Y2)]
cen <- t(vapply(split(seq_len(nrow(Y2)), gv),
                function(i) colMeans(Y2[i, , drop = FALSE]), numeric(2)))
gdir <- cen[nrow(cen), ] - cen[1L, ]      # the direction the groups separate in

# Angle to the LINE, not the ray: the dashed line is drawn through the origin both
# ways, so a query pointing "backwards" along it agrees just as strongly.
angl <- function(u) 180 / pi * acos(pmin(1, abs(sum(u * gdir)) /
                                    sqrt(sum(u^2) * sum(gdir^2))))
W2 <- R2$directions[1:2, , drop = FALSE]
aa <- vapply(colnames(W2), function(q) angl(W2[, q]), numeric(1))
round(aa, 1)                              # degrees from the group direction

## ----acc_map_angles_read, echo=FALSE, results='asis'--------------------------
an <- aa[grepl("^noise", names(aa))]
cat(sprintf(paste0(
  "Measured from the figure's own coordinates: **tumour stage sits %.0f\u00b0 from ",
  "the dashed line**, while the five noise columns land at %.0f\u00b0, %.0f\u00b0, ",
  "%.0f\u00b0, %.0f\u00b0 and %.0f\u00b0. %s\n"),
  aa[["stage_score"]], an[1], an[2], an[3], an[4], an[5],
  if (aa[["stage_score"]] < min(an))
    paste0("Stage is closer to the group direction than **every** noise draw -- ",
           "which is the comparison that matters, and it is the one the single ",
           "number `T` failed to make.")
  else
    paste0("**Caution:** at least one noise draw lies as close to the group ",
           "direction as tumour stage does. That is not a defect in the figure: ",
           "with this sample size a random vector can fall anywhere, which is ",
           "precisely why one control column is never enough.")))

## ----acc_ind_head-------------------------------------------------------------
bi <- sens_acc$by_individual
head(bi[bi$query == "stage_score", ], 4)

## ----acc_ind_between----------------------------------------------------------
length(unique(bi$S_between[bi$query == "stage_score"]))   # 4 stages x 3 assays

## ----acc_ind_rows-------------------------------------------------------------
length(unique(bi$id))          # participants in the decomposition
nrow(acc$Y)                    # participants in the fit
sum(is.na(stage))              # not covered by the grouping

## ----acc_ind_summary----------------------------------------------------------
summary(sens_acc)

## ----acc_ind_plot, fig.wide=TRUE, fig.height=9, dpi=120, fig.cap="One panel per query. Bars are stacked by assay, and the pale bar at the foot of each panel is everyone outside the top ten."----
plot(sens_acc, type = "individuals", top_n = 10)

## ----acc_ind_conc-------------------------------------------------------------
d   <- bi[bi$query == "stage_score", ]
per <- sort(rowsum(d$S_total, d$id)[, 1], decreasing = TRUE)
# `top_n` only changes how many bars are named; the remainder bar absorbs the
# rest either way, so the panel always accounts for the full total.
round(c(top10 = sum(head(per, 10)), all_others = sum(per[-(1:10)])) / sum(per), 3)

## ----acc_ind_absent-----------------------------------------------------------
range(d$S_total[d$S_total < 1e-20])

## ----cardio_stab--------------------------------------------------------------
stab_ok <- mgccaStability(fit, B = 30, fraction = 0.90, seed = 1)
stab_ok

## ----cardio_stab_read, echo=FALSE, results='asis'-----------------------------
s1 <- stab_ok$summary_metrics
cat(sprintf(paste0(
  "Worst case %.4f over %d resamples, with a rank margin of %.4f. Read that as: ",
  "had you recruited a somewhat different set of individuals, you would have ",
  "recovered the same axes, and the comparison itself was well posed throughout. ",
  "This is what it looks like when a result does not depend on who happened to be ",
  "sampled.\n"),
  s1$overlap_min, s1$n_resamples, s1$rank_margin_min))

## ----acc_stab-----------------------------------------------------------------
stab <- mgccaStability(acc, method = "both", group = stage,
                       B = 40, fraction = 0.632, seed = 1)
stab

## ----acc_stab_read, echo=FALSE, results='asis'--------------------------------
sm <- stab$summary_metrics
ss <- sm[sm$method == "subsample", ]; lo <- sm[sm$method == "loco", ]
# The verdict is DERIVED from the number, never printed beside it. A canned
# sentence that says "essentially the same axes" whatever the minimum turns out
# to be is not a reading: it is decoration sitting next to data, and sooner or
# later it contradicts the figure above it.
verdict <- function(x) {
  if (x >= 0.99) return("the axes are effectively unchanged")
  if (x >= 0.95) return("the axes move only slightly")
  if (x >= 0.80) return("the axes move noticeably, so individual directions should be interpreted with care")
  "the axes are NOT stable here, and individual directions should not be interpreted"
}
loco_line <- if (lo$overlap_min >= 0.95) {
  "no single group is holding the result up on its own."
} else {
  "at least one group matters a good deal, and dropping it moves the axes appreciably."
}
cat(sprintf(paste0(
  "Here the median is %.4f but the worst case is %.4f. At the median, %s. ",
  "**The worst case is the number that matters**, and at %.4f, %s.\n\n",
  "Leaving out an entire stage group at a time gives a worst case of %.4f: %s\n\n",
  "This is not a defect of the data. `miniACC` has 80 patients, and a third of ",
  "them removed at once is a large perturbation for a sample that size. The point ",
  "is that the diagnostic *said so* instead of returning a comfortable number.\n"),
  ss$overlap_median, ss$overlap_min, verdict(ss$overlap_median),
  ss$overlap_min, verdict(ss$overlap_min), lo$overlap_min, loco_line))

## ----acc_stab_plot, fig.height=4.2--------------------------------------------
plot(stab)

## ----sweep--------------------------------------------------------------------
sw <- do.call(rbind, lapply(c(0.632, 0.80, 0.90, 0.95), function(f) {
  m <- mgccaStability(fit, B = 20, fraction = f, seed = 1)$summary_metrics
  data.frame(fraction = f, worst_overlap = round(m$overlap_min, 4),
             median_overlap = round(m$overlap_median, 4),
             worst_rank_margin = round(m$rank_margin_min, 4))
}))
knitr::kable(sw, caption = "Cardiovascular fit: the same data at four resampling severities.")

## ----backend, eval=FALSE------------------------------------------------------
# mgccaSensitivity(acc, Q, group = stage, backend = "hdf5")

## ----store, eval=FALSE--------------------------------------------------------
# sens_acc <- mgccaSensitivity(acc, Q, group = stage,
#                              store = TRUE, store_name = "stage_sensitivity")

## ----session------------------------------------------------------------------
sessionInfo()