Reading the numbers without the algebra

In what order to look at things

There is a sequence, and skipping to the middle of it is how results get overinterpreted.

Step 1 – is the subspace itself trustworthy? Run mgcca_stability() first. If the axes move substantially when you resample your own participants, nothing that follows should be interpreted: you would be reading a structure that a slightly different sample would not have produced. Section @ref(is-the-subspace-reproducible) shows one trustworthy case and one that is not.

Step 2 – does the omic signal know about your variable at all? Read T against the noise control. If it is not clearly above it, stop: there is nothing to explain.

Step 3 – what is holding it up? Now R_between and the per-block profile.

The order matters: each step only makes sense if the one above it passed.

Does the omic signal know anything about my variable?

A query is any participant-level vector or matrix. It must carry individual identifiers: participants are matched by name, never by position. That is not fussiness. A misaligned query returns a perfectly plausible number rather than an error, and nothing in the output would let you catch it.

We ask about two clinical variables, and add pure noise as a control – a reference that should show almost no alignment, so the others can be read against something.

# 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 <- mgcca_sensitivity(fit, Q)
sens

## mgcca sensitivity
##   queries   : 7
##   blocks    : 3
##   grouping  : none -- total sensitivity only
##   individuals: 646 per query (by_individual, exploratory)
##   backend   : memory ; L = 2 ; ridge = guarded rule
##   balancing : top-L block weights, methylation 0.781, clinical 1.04, other 1.18
## 
##    query n_used         T  S_total
##      LDL    633 0.0086246 0.019333
##  glucose    633 0.0032796 0.007906
##   noise1    633 0.0039000 0.007755
##   noise2    633 0.0005551 0.001287
##   noise3    633 0.0043890 0.009611
##   noise4    633 0.0048719 0.011738
##   noise5    633 0.0048777 0.009840

T is the share of the query that lies in the retained subspace, bounded in \([0,1]\). S_total is its first-order sensitivity: how much that alignment would move under a structured perturbation of the blocks.

A note on the paragraphs that follow a printed table in this vignette, here and below: they are generated from the object printed above them, not typed in, so they cannot come to disagree with the numbers they describe – including the verdicts, which are derived from thresholds rather than asserted. Their code is folded away only because five screens of sprintf() between a table and its reading would make the vignette unreadable; it is in the vignette source, and every chunk that computes something a reader needs is shown.

Reading this table. The shared axes account for 0.86% of the variation in LDL and 0.33% of glucose.

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 0.06% and 0.49%. 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.

(Five columns rather than one, because a single random column is itself a random draw: it might land at 0.06% or at 0.49%. Reading a real variable against one lucky draw is exactly how a null result turns into a finding.)

LDL is above the ceiling, but only just – 0.86% against 0.49%. 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.

Glucose does not clear the ceiling at all (0.33% against 0.49%). 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.

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.

S_between and R_between are NA here, and deliberately so: no grouping was supplied, so the group-dependent quantities are absent rather than filled with a zero that would look like a measurement. cardiovascular carries no batch, centre or run variable, which is exactly the common situation in which only the total is available. Section @ref(a-real-grouping) uses a data set that does have one.

plot(sens, type = "overall", noise = paste0("noise", 1:5))

Only a query whose point sits outside the grey band is saying anything the machinery would not have produced from random numbers.

How to read this figure. One row per variable, and the grey band is the level the control columns reached – what chance alone gives in this analysis. A point inside the band is indistinguishable from noise; a point clear of it is not. It is the same comparison the table above makes, but you can take it in at a glance, and it scales: with two hundred exposures the table becomes unreadable and this does not.

Which data type is the finding coming from?

prof <- mgcca_block_profile(sens)
prof

## mgcca block profile
##   queries : 7 ; blocks : 3 ; grouping : none
## 
##    query       block   S_total share_total
##      LDL methylation 0.0057460      0.2972
##      LDL    clinical 0.0041883      0.2166
##      LDL       other 0.0093988      0.4862
##  glucose methylation 0.0024817      0.3139
##  glucose    clinical 0.0019507      0.2467
##  glucose       other 0.0034738      0.4394
##   noise1 methylation 0.0028779      0.3711
##   noise1    clinical 0.0009187      0.1185
##   noise1       other 0.0039584      0.5104
##   noise2 methylation 0.0004082      0.3172
##   noise2    clinical 0.0003093      0.2404
##   noise2       other 0.0005692      0.4424
##   noise3 methylation 0.0027124      0.2822
##   noise3    clinical 0.0023281      0.2422
##   noise3       other 0.0045701      0.4755
##   noise4 methylation 0.0043291      0.3688
##   noise4    clinical 0.0024251      0.2066
##   noise4       other 0.0049835      0.4246
##   noise5 methylation 0.0028775      0.2924
##   noise5    clinical 0.0019665      0.1998
##   noise5       other 0.0049963      0.5077
## 
##   No grouping was supplied to mgcca_sensitivity(), so the between-group
##   columns are not shown: they are undefined here, not missing.

mgcca_block_profile() extracts; it does not recompute. The block-level quantities come back from the same pass that produced the totals, so there is one route to them. Two independent routes would eventually disagree, and the disagreement would be silent.

Read the shares next to the absolute values. A block carrying most of a negligible sensitivity is still carrying most of it, and the share alone would not tell you which situation you are in.

plot(sens)

How to read this figure. One panel per variable, one bar per table. The bar height is how much that table contributes to the variable’s sensitivity, so the tallest bar names the data type the finding actually rests on. In a paper this is the difference between “the integrated signal relates to LDL” and “the relation is essentially methylation, with the clinical table adding little” – and only the second is a statement a reviewer can act on.

Here the largest contributor to LDL is other (49% of its total sensitivity).

Is my result leaning on those groups?

Now the grouped quantities are defined. We ask about tumour stage as an ordered score, against the same noise control – five columns of it, for the reason spelt out below.

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 <- mgcca_sensitivity(acc, Q, group = stage)
sens_acc

## mgcca sensitivity
##   queries   : 6
##   blocks    : 3
##   grouping  : 4 levels (stage i, stage ii, stage iii, stage iv)
##   individuals: 78 per query (by_individual, exploratory)
##   backend   : memory ; L = 2 ; ridge = guarded rule
##   balancing : top-L block weights, mRNA 0.916, protein 1.27, miRNA 0.817
## 
##        query n_used        T  S_total S_between R_between
##  stage_score     78 0.018475 0.026627 0.0041339   0.15525
##       noise1     80 0.007661 0.022952 0.0007974   0.03474
##       noise2     80 0.039144 0.106382 0.0070629   0.06639
##       noise3     80 0.001145 0.002458 0.0001064   0.04330
##       noise4     80 0.076076 0.173660 0.0049349   0.02842
##       noise5     80 0.007696 0.021166 0.0006524   0.03082

Reading this table. The three assays together account for 1.85% of the variation in tumour stage, against 0.11–7.61% for the five random columns. Of everything that could shift that reading, 15.5% follows stage itself as a grouping (R_between = 0.155) – a low value.

Caution: note where the alignment sits relative to the noise columns: inside their range. 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.

What to do with that number, whatever it had turned out to be. This is the part worth having in advance, because the temptation is always to decide what a value means after seeing it.

The case in front of us. A high R_between means the alignment and the group structure move together. That is what you would expect from a genuine biological gradient – patients at the same stage really do resemble each other – and also exactly what a batch confounded with stage would produce. The analysis cannot separate those two. You can, because you know how the samples were collected.

If you cannot rule the technical explanation out, the honest options are: report the alignment as conditional on the grouping, or re-run with the grouping residualised out and report what is left. Both are defensible. Reporting the raw alignment without mentioning the grouping is not.

R_between is the share of the sensitivity carried by the grouping. Read it with its two components, never alone. A high ratio can sit on a negligible absolute sensitivity, and a low ratio on a large one; the ratio says how the sensitivity is distributed, not how much of it there is. That is why S_total and S_between are printed beside it rather than being folded away.

plot(sens_acc)

How to read this figure. Each bar is one assay. The height says which assay the finding rests on; the inner bar and its percentage say how much of that assay’s contribution follows the grouping. An assay with a tall bar and a small percentage is carrying signal that does not depend on the groups – the most comfortable position to be in. A tall bar with a large percentage is the one to look at twice.

The bars are the total sensitivity per block; the inner bar and the label are the part carried by the grouping. They are drawn overlaid rather than side by side for a plain reason: the between-group part is routinely two orders of magnitude shorter, so side by side it renders as an invisible sliver and the figure fails to show the one comparison it exists for.

Seeing it instead of inferring it

R_between says how much of the reading follows the grouping. The map view says in which direction, and that is often the more persuasive of the two.

plot(sens_acc, type = "map")

Individuals on the reference axes, group centroids as diamonds, and each query as an arrow.

Each point is a patient placed on the axes the assays agree on; the diamonds are the group centroids; each arrow is the direction a query represents.

What to look for. The dashed line is the direction the groups separate in. If a query’s arrow runs along it, the query and the grouping are describing the same movement – which is exactly what a high R_between measures. If it points elsewhere, they are not. Arrow lengths mean nothing here: they are all scaled alike, and only the angle carries information.

The angles quoted below are measured, not read off the picture: two people looking at the same arrows will not agree to within ten degrees, and this comparison lives or dies on about that much. Everything the calculation needs is on the result itself – sens_acc$reference carries the axes, the grouping and the query directions the figure was drawn from.

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

## stage_score      noise1      noise2      noise3      noise4      noise5 
##        11.3        83.2        35.6        79.2        89.3        69.3

Measured from the figure’s own coordinates: tumour stage sits 11° from the dashed line, while the five noise columns land at 83°, 36°, 79°, 89° and 69°. 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.

And something the numbers do not tell you. Look at how much the ellipses overlap. Stage has a direction in this space, but the groups are not separated along it – individual patients of different stages are thoroughly mixed. A clinician takes that in immediately, and it tempers the reading in a way no single coefficient can: a real but weak gradient, not a clean separation.

With more than two components, call plot(..., comps = c(1, 3)) and so on. Three dimensions have to be rotated to be read and do not survive a printed page; beyond three there is nothing to rotate.

Which participants is it coming from?

The per-block profile says which data type carries the sensitivity. There is a second cut of exactly the same quantity: S_total is a sum over individuals as much as over blocks, and by_individual reports the per-participant terms of that sum. Same arithmetic, different axis.

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

##         query           id block      S_total    S_between  share_total
## 1 stage_score TCGA-OR-A5J1  mRNA 3.334040e-06 1.398997e-05 3.322697e-04
## 2 stage_score TCGA-OR-A5J2  mRNA 1.283954e-05 7.349934e-05 1.279586e-03
## 3 stage_score TCGA-OR-A5J3  mRNA 6.877110e-07 8.754045e-06 6.853714e-05
## 4 stage_score TCGA-OR-A5J4  mRNA 2.895208e-08 7.349934e-05 2.885358e-06
##   share_between
## 1   0.006743561
## 2   0.035428769
## 3   0.004219698
## 4   0.035428769

One row per participant, per block, per query. S_total is that participant’s contribution to that block’s sensitivity; share_total is the same thing as a fraction of the block total, so within one query and block the shares add to one.

Caution: key this table by query, id and block, never by row position. Row order is not part of the contract. It is the same discipline the query itself is held to, and for the same reason: a positional read produces a plausible number rather than an error.

Caution: S_between is a property of the group, not of the person. The between-group part of the sensitivity is built from group means, so every member of a group carries the identical value – here twelve distinct numbers for four stages times three assays, repeated across the participants in each:

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

## [1] 12

So read S_total per participant, and S_between as the group’s contribution carried in each of its members’ rows. Ranking individuals by S_between ranks their groups, which is a different statement.

Caution: the rows are the participants of the grouped set, not of the fit. Anyone the grouping does not cover is not in the calculation at all, so the numbers are conditional on that set:

length(unique(bi$id))          # participants in the decomposition

## [1] 78

nrow(acc$Y)                    # participants in the fit

## [1] 80

sum(is.na(stage))              # not covered by the grouping

## [1] 2

That is a change of question, not a defect. A grouping covering fewer people gives different contributions because it is a different conditioning set.

summary() prints the leading participants for every query, next to the per-block table they decompose:

summary(sens_acc)

## mgcca sensitivity
##   queries   : 6
##   blocks    : 3
##   grouping  : 4 levels (stage i, stage ii, stage iii, stage iv)
##   individuals: 78 per query (by_individual, exploratory)
##   backend   : memory ; L = 2 ; ridge = guarded rule
##   balancing : top-L block weights, mRNA 0.916, protein 1.27, miRNA 0.817
## 
##        query n_used        T  S_total S_between R_between
##  stage_score     78 0.018475 0.026627 0.0041339   0.15525
##       noise1     80 0.007661 0.022952 0.0007974   0.03474
##       noise2     80 0.039144 0.106382 0.0070629   0.06639
##       noise3     80 0.001145 0.002458 0.0001064   0.04330
##       noise4     80 0.076076 0.173660 0.0049349   0.02842
##       noise5     80 0.007696 0.021166 0.0006524   0.03082
## 
## Per-block decomposition
##        query   block   S_total S_between share_total share_between
##  stage_score    mRNA 0.0100341 2.075e-03      0.3768       0.50185
##  stage_score protein 0.0088889 7.669e-04      0.3338       0.18551
##  stage_score   miRNA 0.0077041 1.292e-03      0.2893       0.31264
##       noise1    mRNA 0.0081930 3.504e-04      0.3570       0.43946
##       noise1 protein 0.0073492 8.561e-05      0.3202       0.10736
##       noise1   miRNA 0.0074098 3.614e-04      0.3228       0.45318
##       noise2    mRNA 0.0332003 3.360e-03      0.3121       0.47576
##       noise2 protein 0.0380220 1.290e-03      0.3574       0.18258
##       noise2   miRNA 0.0351596 2.413e-03      0.3305       0.34166
##       noise3    mRNA 0.0009452 5.028e-05      0.3846       0.47248
##       noise3 protein 0.0008498 1.937e-05      0.3458       0.18201
##       noise3   miRNA 0.0006628 3.677e-05      0.2697       0.34552
##       noise4    mRNA 0.0579982 2.302e-03      0.3340       0.46637
##       noise4 protein 0.0719257 8.950e-04      0.4142       0.18136
##       noise4   miRNA 0.0437366 1.738e-03      0.2519       0.35227
##       noise5    mRNA 0.0080281 3.945e-04      0.3793       0.60475
##       noise5 protein 0.0073872 3.920e-05      0.3490       0.06008
##       noise5   miRNA 0.0057510 2.187e-04      0.2717       0.33517
## 
## Per-individual sensitivity contributions (exploratory)
##   query 'stage_score' -- top 5 of 78 participants by S_total
##            id         S share_of_query
##  TCGA-OR-A5KZ 0.0044005        0.16526
##  TCGA-P6-A5OG 0.0015148        0.05689
##  TCGA-OR-A5J8 0.0013237        0.04971
##  TCGA-PA-A5YG 0.0011402        0.04282
##  TCGA-OR-A5LD 0.0009808        0.03684
##   query 'noise1' -- top 5 of 78 participants by S_total
##            id         S share_of_query
##  TCGA-P6-A5OG 0.0071935        0.31341
##  TCGA-OR-A5J8 0.0040138        0.17488
##  TCGA-OR-A5JO 0.0013022        0.05674
##  TCGA-OR-A5JB 0.0012947        0.05641
##  TCGA-OR-A5J2 0.0007724        0.03365
##   query 'noise2' -- top 5 of 78 participants by S_total
##            id        S share_of_query
##  TCGA-OR-A5KZ 0.009491        0.08922
##  TCGA-P6-A5OG 0.006853        0.06442
##  TCGA-OR-A5JO 0.004406        0.04141
##  TCGA-OR-A5KY 0.004021        0.03780
##  TCGA-OR-A5J2 0.003945        0.03708
##   query 'noise3' -- top 5 of 78 participants by S_total
##            id         S share_of_query
##  TCGA-P6-A5OG 8.478e-04        0.34493
##  TCGA-OR-A5J8 4.997e-04        0.20331
##  TCGA-OR-A5JB 1.440e-04        0.05860
##  TCGA-OR-A5JO 1.294e-04        0.05266
##  TCGA-PK-A5H9 7.803e-05        0.03175
##   query 'noise4' -- top 5 of 78 participants by S_total
##            id        S share_of_query
##  TCGA-P6-A5OG 0.054812        0.31563
##  TCGA-OR-A5J8 0.030404        0.17508
##  TCGA-OR-A5JO 0.009048        0.05210
##  TCGA-OR-A5JB 0.007863        0.04528
##  TCGA-PK-A5H9 0.005841        0.03363
##   query 'noise5' -- top 5 of 78 participants by S_total
##            id         S share_of_query
##  TCGA-P6-A5OG 0.0049069        0.23183
##  TCGA-OR-A5J8 0.0025041        0.11831
##  TCGA-OR-A5JO 0.0012698        0.05999
##  TCGA-PK-A5H9 0.0007358        0.03476
##  TCGA-OR-A5JB 0.0007227        0.03414
## 
##   S is summed over blocks; share_of_query is its share of the query's
##   total sensitivity (the by_individual table's own share_total is a
##   share WITHIN one block, and sums to one there).
##   EXPLORATORY: the block-level parent has MIXTO calibration evidence,
##   and no individual-level calibration has been established. These are
##   sensitivity contributions -- not variance, uncertainty, causal
##   influence, or reliability guarantees -- conditioned on the grouped
##   participant set.
## 
##   T is the share of the query lying in the retained subspace.
##   S_total is its first-order sensitivity; S_between the part carried by the
##   grouping. Read R_between WITH its components: a high ratio can accompany a
##   negligible absolute sensitivity, and a low one a large signal.

And the whole distribution is easier seen than tabulated:

plot(sens_acc, type = "individuals", top_n = 10)

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.

How to read this figure. One bar per participant, stacked by assay, so you see who carries the sensitivity and which assay it reaches them through at the same time. Compare the real query against its noise controls: a query whose sensitivity rests on a handful of people looks very different from one spread evenly across the cohort, and that difference is invisible in S_total alone.

The pale “all others” bar is the point of the figure, not decoration. A top-N chart without the remainder silently overstates concentration – ten tall bars look like the whole story whether they carry most of the total or a tenth of it. Drawn together, the comparison is immediate:

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)

##      top10 all_others 
##      0.497      0.503

Caution: participants absent from a block contribute a numerical zero there, not an exact one. miniACC is full of them – most participants are missing from at least one assay – and their contribution to that assay comes back at the scale of rounding error rather than at 0:

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

## [1] 5.82801e-38 5.82801e-38

They are deliberately not rounded down: the value is defined, it simply happens to be negligible, and thresholding it would hide the difference between “absent” and “present and unimportant”.

What these per-participant numbers are, and what they are not. They are an exploratory decomposition of the first-order sensitivity diagnostic. The block-level quantity they decompose has empirical calibration evidence, with a MIXTO verdict; no individual-level calibration has been established, and the block-level evidence does not transfer to the individual level on its own. They are sensitivity contributions – not variance, not uncertainty, not causal influence, and not a guarantee of correctness or reliability for any participant. A large bar says the fitted alignment is more sensitive to that participant’s measurements under this configuration; it says nothing about that person, and it does not rank people in any way that survives a change of fit, ridge, query, block or grouping.

Read this view as a diagnostic pointer: if one or two participants carry a large share of a finding, that is worth looking at – an outlying sample, a processing anomaly – before the finding is reported as a property of the cohort.

A result you can lean on

The cardiovascular fit has 646 individuals. Holding out 10% of them at a time:

stab_ok <- mgcca_stability(fit, B = 30, fraction = 0.90, seed = 1)
stab_ok

## mgcca subspace stability
##   method   : subsample ; L = 2 ; backend = memory
##   subsample: B = 30, fraction = 0.900
## 
##     method n_resamples n_failed overlap_min overlap_median rank_margin_min
##  subsample          30        0      0.9973         0.9987          0.9598
## 
##   overlap = agreement with the reference subspace on the SHARED participants.
##   rank_margin = ratio of the smallest to the largest singular value of the
##   reference frame restricted to them (its reciprocal condition number): it
##   says whether that comparison is well posed, NOT whether the components are
##   separated. A random share of an orthonormal frame's rows keeps this ratio
##   near 1, so a value well below 1 signals genuine ill-conditioning.

Worst case 0.9973 over 30 resamples, with a rank margin of 0.9598. 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.

A result you cannot

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

## mgcca subspace stability
##   method   : both ; L = 2 ; backend = memory
##   subsample: B = 40, fraction = 0.632
##   groups   : stage i, stage ii, stage iii, stage iv
## 
##     method n_resamples n_failed overlap_min overlap_median rank_margin_min
##       loco           4        0      0.9482         0.9718          0.7163
##  subsample          40        0      0.5440         0.9749          0.5057
## 
##   overlap = agreement with the reference subspace on the SHARED participants.
##   rank_margin = ratio of the smallest to the largest singular value of the
##   reference frame restricted to them (its reciprocal condition number): it
##   says whether that comparison is well posed, NOT whether the components are
##   separated. A random share of an orthonormal frame's rows keeps this ratio
##   near 1, so a value well below 1 signals genuine ill-conditioning.

Here the median is 0.9749 but the worst case is 0.5440. At the median, the axes move only slightly. The worst case is the number that matters, and at 0.5440, the axes are NOT stable here, and individual directions should not be interpreted.

Leaving out an entire stage group at a time gives a worst case of 0.9482: at least one group matters a good deal, and dropping it moves the axes appreciably.

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.

plot(stab)

Plotting overlap against rank margin makes the caution explicit: the points to distrust are the ones on the left, whatever their height.

What made the difference

The difference is not that one data set is good and the other bad. It is how hard each was pushed, and how much sample there was to absorb it. To show that, the cardiovascular fit is put through the same diagnostic at four resampling severities – from holding out 37% of the individuals at a time down to 5% – with everything else held fixed:

sw <- do.call(rbind, lapply(c(0.632, 0.80, 0.90, 0.95), function(f) {
  m <- mgcca_stability(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.")

Two things to take from this table.

The first is that stability is a property of the data and of the question you ask of it. Reporting an overlap without the fraction it was computed at says almost nothing.

The second is what the rank margin adds: it stays close to one at every severity, with its lowest value in the table at 0.9096. The reference frame therefore remains well conditioned on the shared participants even when 37% of them are held out, and the lower overlap at that severity is not an artefact of an ill-posed comparison: with fewer individuals, the axes themselves move more. The margin drifts down a little as the resamples shrink, because a smaller subset is a noisier restriction of the frame, but it does not shrink with the resample size by construction. Had it fallen well below one, a low overlap could not have been read as instability, because the comparison itself would have been in doubt. Reading the two columns together is what tells the two explanations apart.

Where the computation happens

backend = "auto" decides by block size, exactly as mgcca(route = "auto") does, but "memory" and "hdf5" force the choice. The backend that ran is recorded in the result, so a number is never ambiguous about how it was produced.

mgcca_sensitivity(acc, Q, group = stage, backend = "hdf5")

block_size forces the out-of-core block size. Its normal value is 0 (let BigDataStatMeth choose); a small positive value exists so the blocked path can be exercised on small data, which is how the package tests it.

Keeping your results

Persistence is opt-in. An analysis function should not write to your file unasked, particularly when the fit may be shared or read-only.

sens_acc <- mgcca_sensitivity(acc, Q, group = stage,
                              store = TRUE, store_name = "stage_sensitivity")

Results go under RELIABILITY/<name>/ in the fit’s own HDF5 file, with the same attribute manifest the package uses elsewhere. An existing name is refused rather than overwritten: two runs with different settings under one name would leave a manifest that no longer describes the data beside it, and nothing downstream could detect that.

What you can and cannot claim

These are first-order sensitivity coefficients of a fitted operator, under the scaling, weighting and availability conventions of the fit that produced them. They support statements about where a query’s alignment is sensitive, and about whether a subspace is reproducible. They do not support causal statements, and a high R_between does not by itself mean an association is an artefact – it means the alignment and the grouping move together in this fit.

When the retained and discarded parts of the spectrum are not separated, the first-order expansion is not meaningful at all. The result carries a validity flag and the print method says so rather than letting the numbers stand unqualified.