R/pca_sample_distance.R
pca_sample_distance.RdComputes a robust multivariate distance for each sample from a center of
the PCA score cloud, across the first n_pcs principal components.
Residuals from the center are converted to robust z-scores per PC; the
distance is the Euclidean norm of those z-scores, i.e. the Mahalanobis
distance in PC space with robust location/scale. The value is descriptive
guidance for sample review — it does not classify or exclude samples.
pca_sample_distance(
pca,
n_pcs = 5,
reference_quantile = 0.975,
groups = NULL,
min_group_size = 3
)A prcomp or summary.prcomp result (any list with a
score matrix in $x; samples in rows, PCs in columns; no NA values).
Missing row/column names fall back to row numbers / PC1, PC2, ....
Number of leading PCs to use. Clamped to the number of PCs available. Default 5.
Chi-squared quantile used for the unusual
reference value. Default 0.975.
Optional vector (one value per row of pca$x) assigning
each sample to a group. If named, it is aligned to rownames(pca$x)
by name (an error if the names do not match); if unnamed, it is taken in
row order. Numeric vectors are treated as continuous covariates: a warning
is issued and the global center is used (pass as.character() for
numeric category codes). Default NULL (global center for all samples).
Minimum group size for a group to get its own center. Default 3.
A list:
Data frame with sampleName,
distance, p_value (upper-tail probability of
distance^2 under the same chi-squared reference used for
typical/unusual), p_adjusted (p_value
after Benjamini-Hochberg correction across every sample in this call —
corrects for testing all plotted samples at once, unlike the fixed
unusual threshold), largest_pc (the used PC with the
largest absolute robust z-score, i.e. where the sample deviates most),
and center ("group" when the sample was measured from its own
group's center, "global" otherwise).
Character vector of the PCs used.
Number of PCs used (chi-squared degrees of freedom).
Distance of a typical sample: sqrt(qchisq(0.5, n_pcs)).
Reference distance beyond which a sample is unusual:
sqrt(qchisq(reference_quantile, n_pcs)).
When groups is supplied, each sample is measured from the center of
its own group (per-group median per PC) instead of the global center, so
members of distinct-but-healthy clusters (runs, cohorts) are not scored as
far just for belonging to a small group. The within-group scale is pooled
across all qualifying groups using pairwise differences (see Details), so
small groups do not bias it. Samples whose group is too small (or has an NA
label) are scored exactly as in the ungrouped run — global center, global
scale — and marked "global" in the center column.
The global scale per PC is the MAD of the scores. The grouped
scale per PC is the median of all within-group pairwise absolute
differences across qualifying groups, divided by
sqrt(2) * qnorm(0.75) for normal consistency — unlike a MAD of
median-centered residuals, this has no structural zeros and stays
calibrated even for groups of 3 (simulated false-"unusual" rates stay at
the nominal level for group sizes 3-10).
PCs with (near-)zero scale on either estimator are skipped — both exactly
degenerate score distributions and the float-noise components a
rank-deficient PCA leaves behind (relative tolerance 1e-8 of the
largest scale); the chi-squared references use the number of PCs actually
used. Under approximate normality the squared distance follows a
chi-squared distribution with n_pcs degrees of freedom, which
provides the typical (median) and unusual
(reference_quantile) reference values. For very small groups the
group-median center makes distances slightly conservative (a group-of-3
midpoint scores 0 on each PC).
mat <- matrix(rnorm(200), nrow = 20)
rownames(mat) <- paste0("S", 1:20)
pca <- prcomp(mat, center = TRUE, scale. = TRUE)
pca_sample_distance(pca)
#> $sample_distances
#> sampleName distance p_value p_adjusted largest_pc center
#> 1 S1 2.4165591 0.3221248105 0.585681474 PC2 global
#> 2 S2 2.4342795 0.3135142810 0.585681474 PC4 global
#> 3 S3 3.7355694 0.0159011742 0.079505871 PC4 global
#> 4 S4 0.8234224 0.9841529413 0.984152941 PC1 global
#> 5 S5 2.4646835 0.2990213939 0.585681474 PC1 global
#> 6 S6 2.5379120 0.2656402524 0.585681474 PC1 global
#> 7 S7 3.5188982 0.0299046500 0.119618600 PC3 global
#> 8 S8 1.3265028 0.8813060160 0.927690543 PC3 global
#> 9 S9 2.5740973 0.2499754533 0.585681474 PC4 global
#> 10 S10 4.3778062 0.0017905920 0.017905920 PC4 global
#> 11 S11 1.7340258 0.6989303872 0.877293203 PC1 global
#> 12 S12 1.4824007 0.8211952816 0.927690543 PC2 global
#> 13 S13 4.8175185 0.0003079502 0.006159003 PC4 global
#> 14 S14 4.0661748 0.0054743798 0.036495865 PC4 global
#> 15 S15 1.3965488 0.8559717481 0.927690543 PC2 global
#> 16 S16 2.0211601 0.5372315739 0.767473677 PC3 global
#> 17 S17 2.6125147 0.2339658107 0.585681474 PC3 global
#> 18 S18 2.1007601 0.4915812442 0.756278837 PC2 global
#> 19 S19 1.7285859 0.7018345624 0.877293203 PC4 global
#> 20 S20 2.1176545 0.4819652423 0.756278837 PC5 global
#>
#> $pcs_used
#> [1] "PC1" "PC2" "PC3" "PC4" "PC5"
#>
#> $n_pcs
#> [1] 5
#>
#> $typical
#> [1] 2.086015
#>
#> $unusual
#> [1] 3.582248
#>
pca_sample_distance(pca, groups = rep(c("Run1", "Run2"), each = 10))
#> $sample_distances
#> sampleName distance p_value p_adjusted largest_pc center
#> 1 S1 2.1496973 0.46383303 0.6678954 PC2 group
#> 2 S2 1.6012008 0.76685083 0.9444060 PC4 group
#> 3 S3 2.2368233 0.41546800 0.6678954 PC1 group
#> 4 S4 0.8712030 0.97957473 0.9853670 PC1 group
#> 5 S5 2.1431477 0.46752678 0.6678954 PC1 group
#> 6 S6 2.3293639 0.36613562 0.6678954 PC1 group
#> 7 S7 3.0253765 0.10311860 0.6678954 PC3 group
#> 8 S8 1.5244460 0.80274508 0.9444060 PC2 group
#> 9 S9 2.3141412 0.37407941 0.6678954 PC5 group
#> 10 S10 2.8125670 0.16123574 0.6678954 PC2 group
#> 11 S11 2.2019708 0.43462597 0.6678954 PC5 group
#> 12 S12 1.8141776 0.65518302 0.8735774 PC2 group
#> 13 S13 2.4026577 0.32896167 0.6678954 PC4 group
#> 14 S14 3.2837473 0.05585639 0.6678954 PC5 group
#> 15 S15 0.8090755 0.98536696 0.9853670 PC2 group
#> 16 S16 2.1843027 0.44443932 0.6678954 PC3 group
#> 17 S17 2.2438585 0.41163570 0.6678954 PC3 group
#> 18 S18 2.2680312 0.39856447 0.6678954 PC2 group
#> 19 S19 1.3379105 0.87737130 0.9748570 PC4 group
#> 20 S20 2.5428297 0.26347848 0.6678954 PC5 group
#>
#> $pcs_used
#> [1] "PC1" "PC2" "PC3" "PC4" "PC5"
#>
#> $n_pcs
#> [1] 5
#>
#> $typical
#> [1] 2.086015
#>
#> $unusual
#> [1] 3.582248
#>