---
title: "BatchVaria Quickstart: Evaluating Batch Correction Performance"
author:
- name: Nicholas Moir
  affiliation: University of Edinburgh
  email: nicholas.moir@ed.ac.uk
package: BatchVaria
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
vignette: >
  %\VignetteIndexEntry{BatchVaria Quickstart: Evaluating Batch Correction Performance}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, message=FALSE, warning=FALSE}
library(BatchVaria)
library(SummarizedExperiment)
library(BiocStyle)
```

# Overview

High-throughput transcriptomic datasets frequently contain unwanted technical variation arising from batch effects. Numerous batch correction methods exist, but selecting an appropriate method is challenging because correction can remove biological signal or distort sample relationships. BatchVaria provides a framework for quantifying variance components, tracking variance changes after correction, and evaluating correction performance using multiple metrics, enabling evidence-based selection of batch correction strategies.
<br>
<br>
While `r Biocpkg("sva")`, `r Biocpkg("limma")` and `r Biocpkg("BatchQC")` provide batch correction and diagnostic visualisations, BatchVaria focuses specifically on quantitative evaluation of correction performance across multiple assays and correction strategies within a unified framework. BatchVaria unifies variance decomposition, tracking, and comparison into a single coherent workflow, built on `r Biocpkg("SummarizedExperiment")` and using `r Biocpkg("variancePartition")` as one of its two variance engines.

If you use BatchVaria in published work, please cite Moir N, Sherwood K,
Simpson TI (2026), *BatchVaria: a variance-aware framework for evaluating
batch correction in high-dimensional omics data*, bioRxiv,
[doi:10.64898/2026.05.07.721996](https://doi.org/10.64898/2026.05.07.721996).
`citation("BatchVaria")` prints the same entry.

---

# Installation

Install the release version with `r CRANpkg("BiocManager")`:

```{r installation-release, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("BatchVaria")
```

For the development version:

```{r installation-devel, eval=FALSE}
BiocManager::install(version = "devel")
BiocManager::install("BatchVaria")
```

The development source is at `r Githubpkg("cmpmolonc/BatchVaria")`.

# Example Dataset

We begin with a synthetic dataset with *partially confounded* batch and biology. This creates a situation where effective batch correction can be challenging. Within this dataset, most control samples sit in one batch and most treated samples in the other, so effects are not cleanly separable.

```{r runCorrection}
## exampleBatchVaria() draws random data; seed for a reproducible vignette
set.seed(1)
bv <- exampleBatchVaria(confounding = 0.6)

class(bv)

table(
    batch = bv$batch,
    group = bv$group
)
```

BatchVaria leverages `r Biocpkg("SummarizedExperiment")` for storage. Expression matrices are assays, sample
annotation is `colData()`, and everything the package records goes in
`metadata()`. See `?"BatchVaria-requirements"` for further detail.

Correction methods are registered rather than hard-coded, so the set available
is extensible in the same way the variance engines are:

```{r methods}
availableCorrectionMethods()
```

This vignette uses `combat` throughout, so that the comparison below isolates
the effect of the design rather than differences between methods.

ComBat is applied twice. The first call corrects for batch alone. The second
additionally names `group` as a variable to preserve, so the batch effect is
estimated adjusted for the named preservation variable.

```{r corrections}
bv <- runCorrection(
    bv,
    method = "combat", batch = "batch",
    newAssayName = "naive"
)

bv <- runCorrection(
    bv,
    method = "combat", batch = "batch", preserve = "group",
    newAssayName = "protected"
)

assayNames(bv)
```

`raw_noise` is a deliberately distorted assay supplied with the example data,
kept here as a reference point for what a genuinely destructive transformation
looks like. `raw_center` and `raw_scale` are no-op controls: centring and
scaling change no variance fraction, so a diagnostic that reports a change on
them is reporting an artefact of the diagnostic rather than of the data.

---

# Variance Profiling 

We explicitly quantify the contribution of known covariates.

```{r profile}
assaysUsed <- c("raw", "naive", "protected", "raw_noise")

bv <- profileVariance(
    bv,
    formula = ~ batch + group,
    assays = assaysUsed,
    methods = c("anova", "variancePartition")
)
```
This step fits the specified model to each assay and decomposes variance into
batch (unwanted), group (biological) and residual (unexplained). Two engines are
run: `anova`, which fits linear models per feature and partitions variance with
Type II sums of squares, and `variancePartition`, which fits mixed models. Every
result is stored in an internal variance ledger keyed on assay, method, formula
and term.

Because the ledger now holds two decompositions per assay, summaries must say
which one they mean. Asking for a table without choosing returns an error rather than
a silent mixture, since the engines' terms are not interchangeable:

```{r ambiguous, error=TRUE}
varianceTable(bv, assays = assaysUsed)
```

## What the object is carrying

`provenance()` reports what has been done, reading both ledgers:

```{r provenance}
p <- provenance(bv)
```

Two things there are worth pointing at. The correction table's `preserve`
column is what distinguishes the two ComBat runs - without it the rows are
identical but for the assay names, and the comparison this vignette is built
on would not be reconstructable from the object. The `no_op` column flags a
correction that returned an assay indistinguishable from its input; a method
can succeed and change nothing, and recording it on the object means the
condition outlives the console session that produced it.

It prints that report and returns the same content as a data.frame, so the
history can be filtered or joined rather than only read:

```{r provenance-return}
p$corrections[, c("method", "output", "preserve", "no_op")]
```

---

# Variance Summary Tables

A core output of BatchVaria is a pair of summary tables.

```{r table}
vt <- varianceTable(bv, assays = assaysUsed, method = "anova")
```
---

# Variance Composition Table

This table shows the percentage of variance explained by each component for each assay.

```{r percent}
vt$percent
```
---

# Variance Change Relative to Raw

This table shows how each correction changed variance relative to the raw data.
Every entry is a change in a *fraction*, which constrains what it can be read
to mean; the next section sets out what that constraint is.

```{r delta}
varianceTable(
    bv,
    assays = assaysUsed,
    method = "anova",
    formatDelta = TRUE
)$delta
```
---

# Total Variance: Reading Fractions Correctly

Every number above is a *fraction*. Within an assay, the fractions sum to
one. This has a consequence that must be accounted for before any of these
tables can be interpreted: **removing variance associated with one term
mechanically inflates the fraction attributed to every remaining term**,
even when the absolute variance of those terms has not increased at all.

`assayVariance()` supplies the missing denominator.

```{r totalvar}
av <- assayVariance(bv, assays = assaysUsed)
av
```

Multiplying a fraction by its assay's total variance converts it to an
absolute quantity, which is directly comparable across assays.

```{r absolute}
vr <- varianceResults(bv, method = "anova")

absVar <- merge(
    vr[, c("assay", "term", "variance_fraction")],
    av[, c("assay", "total_variance")],
    by = "assay"
)
absVar$absolute <- round(absVar$variance_fraction * absVar$total_variance, 1)

## as.data.frame(), because pivot_wider() returns a tibble and tibble
## prints three significant figure.
as.data.frame(tidyr::pivot_wider(
    absVar[, c("assay", "term", "absolute")],
    names_from = "term",
    values_from = "absolute"
))
```

Read down the `group` column: 144.6 in `raw`, 126.0 after the naive
correction and 127.0 after the protected one. In absolute terms group variance
barely moves. Its *fraction* meanwhile climbs from 11.3 to 15.2 per cent,
because the denominator fell from 1281 to 829. The fraction rose because
something else was removed, not because biological signal grew.

---

# Two Engines, One Answer - Until the Design Stops Cooperating

BatchVaria harmonises independent variance-profiling engines onto a common
schema: the same columns, and the same names for the two terms that mean the
same thing in any engine - `residual` and `shared`. Both runs above reach
their decompositions by different routes: `anova` fits a linear model per
feature and partitions with Type II sums of squares, while
`variancePartition` fits mixed models via `lme4`.

```{r concordance}
cmp <- varianceResults(bv, assayName = "raw")[
    , c("method", "term", "variance_fraction")
]

as.data.frame(
    tidyr::pivot_wider(
        cmp,
        names_from = "method",
        values_from = "variance_fraction"
    )
)
```

Here they disagree, and the disagreement is the point. On an orthogonal design
they do not:

```{r orthogonal}
orth <- profileVariance(
    exampleBatchVaria(confounding = 0),
    ~ batch + group,
    assays = "raw",
    methods = c("anova", "variancePartition")
)

## Printed at full precision deliberately: at three significant figures
## both columns read 0.163 and the agreement below cannot be checked.
as.data.frame(
    tidyr::pivot_wider(
        varianceResults(orth, assayName = "raw")[
            , c("method", "term", "variance_fraction")
        ],
        names_from = "method",
        values_from = "variance_fraction"
    )
)
```

Every term the two engines share agrees to floating-point precision. `shared`
is exactly zero, because an orthogonal design leaves nothing unattributable,
and `variancePartition` reports `NA` for it because a mixed model has no such
parameter to report.

The divergence above is therefore a property of the design, not of either
implementation. The two engines estimate different quantities:
`variancePartition` reports how much variance is *associated with* a term,
`anova` with Type II sums of squares reports how much is *uniquely attributable
to* it. When a design is orthogonal these coincide, which is why the two tables
differ as they do. When batch and group are partially aliased they cannot coincide,
because roughly 11 per cent of the variance here is consistent with either term
and identifiable to neither.

`anova` reports that share separately, as `shared`. A per-gene mixed model has
no corresponding parameter - and could not, since the effects are not
separately identifiable in this design - so the ambiguous variance is
reflected in both terms, each rising from about 11 per cent to about 17.
Neither number is wrong for the quantity it estimates; they simply stop being
interchangeable once the design cannot separate the effects.

The gap between them is therefore a measurement of the design rather than a
disagreement between tools: it says how much of the variance the experiment
permits you to attribute at all. Running both, and seeing them coincide on a
clean design and separate on a compromised one, is more informative than
trusting either alone.

---

# Interpretation

The naive correction looks like a success. Read the fraction table on its own:
batch variance halves, from 11.0 per cent to 5.6, and group variance *rises*,
from 11.3 per cent to 15.2. Unwanted variation down, biological signal up. By
the standard before-and-after comparison - decompose variance, correct,
decompose again - this is a correction working exactly as intended.

It is not. Three things the fraction table cannot show say otherwise.

**Total variance fell by 35 per cent**, from 1281 to 829. The group *fraction*
rose because the denominator collapsed, not because biology was preserved. In
absolute terms group variance fell from 145 to 126. This is the compositional
artefact: within an assay the fractions must sum to one, so removing variance
from any term mechanically inflates the share of every other. A rising fraction
is not evidence of preserved signal, and reading it as such is the easiest way
to certify a bad correction as a good one.

**The batch effect was not removed.** 5.6 per cent remains, half of what was
there.

**`shared` went negative**, from +11.4 per cent to -3.5. The correction did not
merely fail to resolve the confounding; it distorted the design's
attributability, leaving the model terms explaining more jointly than they do
separately.

Protecting the covariate addresses the second of these. Telling ComBat that
`group` matters lets it estimate the batch effect adjusted for biology, and
batch variance drops to 1.1 per cent - essentially removed - while total
variance falls less far, to 927. It does not rescue biological signal on this
data: absolute group variance is 126 for the naive correction and 127 for the
protected one. The case for protecting the covariate here is that the correction
works, not that the biology survives better.

That balance is specific to this level of confounding. This vignette uses a
single value, `confounding = 0.6`; the trade-off moves with it, because as
batch and group approach aliasing, protecting the covariate increasingly means
protecting the batch effect alongside it. Raising or lowering the parameter
will change which correction looks better on absolute group variance, so treat
that comparison as a property of the design rather than a general ranking of
the two calls.

The sign of `shared` is the more dependable signal. Wherever there is any
confounding to resolve it is consistently negative after naive correction and
positive after protected correction, so a correction that drives it below zero
has distorted the
design's attributability regardless of what the other fractions report. Being a
sign rather than a magnitude, that reading needs no threshold and no reference
distribution - but note it is demonstrated here on simulated data with one
correction method and one design shape.

For contrast, `raw_noise` shows what a destructive transformation looks like
under the same instruments. Its batch and group fractions both fall, which reads
as damage. But its total variance *more than doubles*, to 3077, so in absolute
terms every component grew: nothing was removed, the batch effect was diluted.
`basisRetention()` is unambiguous about the difference.

```{r retention}
basisRetention(bv, assays = assaysUsed)
```

It retains 0.37 of the original variance structure against 0.99 for both ComBat
runs, because added noise occupies directions the raw data did not. Note what
this measure does and does not separate: it distinguishes a correction from a
distortion, not a good correction from a bad one. Both ComBat runs stay on the
original axes because ComBat adjusts each feature's location and scale within a
batch, which rescales the existing structure rather than introducing new
directions.

The general rule: a variance fraction is a statement about *composition*, not
about *quantity*. Use fractions to compare terms within an assay, and totals or
absolute variances to compare the same term across assays.

This is also the answer to what BatchVaria offers over calling the underlying
tools twice. Decomposing variance before and after correction is exactly what
certified the naive correction as successful. Recording each decomposition
alongside the total variance it is relative to, the share the design cannot
attribute, and the structure the correction preserved is what catches it.

---

# Summary Figure

`evaluateCorrections()` runs the diagnostics above in one call and returns them
in a form the plotting layer reads directly.

```{r summaryfig, fig.width=10, fig.height=6, out.width="100%"}
ev <- evaluateCorrections(
    bv,
    baseline = "raw",
    assays = assaysUsed,
    diagnostics = "variance"
)

plotVarianceDelta(ev)
```

One group of bars per model term, one bar per assay, and height is the signed
change in that term's fraction relative to `raw` - as a fraction here, where the
delta table above reported the same quantities as percentages. Terms are drawn
in alphabetical order: `batch`, `group`, `residual`, `shared`.

Both corrections cut `batch`, the protected run roughly twice as hard, and both
raise `residual`. The bar that separates them is `shared`. Both push it down,
and by different amounts - 0.149 for the naive run against 0.044 for the
protected one - but the amount is not the point. What distinguishes them is
where they land: the naive run takes `shared` from +0.114 to **-0.035**, while
the protected run takes it from +0.114 to **+0.070**. As above, a negative
`shared` means the terms explain more jointly than they do separately, so it is
the side of zero that carries the finding, not the size of the drop.
`raw_noise` cuts `batch` too, by dilution rather than correction, which is why
no single bar here is sufficient on its own.

# Diagnostic Plots

The argument above is made with tables, deliberately: every claim in it is a
number that can be checked. The plotting layer draws those same quantities, and is
collected here rather than interleaved so that the two are not confused. Each
plot below answers one question.

`evaluateCorrections()` computes all three diagnostic families in one call and
returns them in the form the plotting functions read:

```{r evalall}
ev <- evaluateCorrections(
    bv,
    baseline = "raw",
    assays = assaysUsed,
    varianceMethod = "anova"
)

names(ev)
```

## Where the variance sits

`plotVarianceComposition()` draws each assay's decomposition as a
composition: one stacked bar per assay, summing to 100 per cent. It answers
"what is this assay made of", and needs no baseline - an assay profiled on its
own can be drawn.

```{r composition, fig.width=8, fig.height=5, out.width="100%"}
plotVarianceComposition(bv, assays = assaysUsed, method = "anova")
```

Read it against the totals from `assayVariance()` above, not on its own: the
bars are shares, so `raw_noise` shows a smaller `batch` share than `raw` while
its absolute batch variance has *risen*, from 141 to 242. The share fell only
because the total more than doubled underneath it.

`plotVarianceRadar()` shows the same deltas as the summary figure on polar
axes, one spoke per term. It is the quicker read when several corrections are
being compared at once and the question is the *shape* of the change rather
than its size - two corrections that move the same terms in the same
proportions trace the same outline.

```{r radar, fig.width=7, fig.height=6, out.width="80%"}
plotVarianceRadar(bv, assays = assaysUsed, method = "anova")
```

## Where the samples sit

Every assay is projected onto a single basis fitted on the reference assay,
rather than each being fitted separately. This is what makes the panels
comparable: a per-assay fit gives `PC1` a different meaning in each, so
samples appear to move when only the axes have.

The coordinates are available directly, not only through a plot:

```{r basisprojection}
coords <- basisProjection(
    bv,
    assays = c("raw", "naive"),
    reference = "raw"
)

head(coords)
```

`plotBasisProjection()` draws exactly what that returns, one panel per assay:

```{r basisplot, fig.width=6, fig.height=4.5, out.width="70%", fig.show="hold"}
for (p in plotBasisProjection(
    bv,
    assays = c("raw", "naive", "protected"),
    reference = "raw"
)) {
    print(p)
}
```

The `naive` panel repays a careful look, because it does not show what a
reader expects of a corrected assay: the batches still appear separated. They
have in fact been brought together - the distance between batch means along
`PC1` falls from 1.73 to 0.54 standard deviations. What remains on the panel is
`group`, whose separation *rises* over the same step, from 1.75 to 1.81. At
`confounding = 0.6` the batch and group labels agree for sixteen of the twenty
samples, so a panel coloured by batch draws surviving group structure in batch
colours, and the four samples whose labels disagree are the only ones that give
it away.

Colouring by the other variable shows what is actually there:

```{r basisplot-group, fig.width=6, fig.height=4.5, out.width="70%"}
plotBasisProjection(
    bv,
    assays = "naive",
    reference = "raw",
    colourBy = "group"
)[[1]]
```

This is a limit of the panel rather than of the correction, and it is not
peculiar to this function: wherever batch and biology are confounded, a
projection coloured by batch cannot distinguish a batch effect that survived
from biological signal that was preserved. Here it ranks the two corrections
backwards: set the two corrected panels side by side and `protected` looks
*more* batch-separated than `naive` - 1.31 standard deviations against 0.54 -
while carrying a fifth of the batch variance, 1.1 per cent against 5.6. The
variance table settles what the panel cannot.

There is a second thing these panels cannot show. `raw_noise`, drawn on this
basis, is nearly indistinguishable from `raw`: the noise it added lies almost
entirely *off* the reference basis, so it never enters these coordinates,
despite the assay carrying more than twice the variance. Projecting onto a
basis shows what happened within the retained structure and is silent about
structure that left it, which is why `basisRetention()` is reported alongside
it - scoring `raw_noise` at 0.37 against 0.99 for both ComBat runs.

`plotPCATrajectory()` makes the same comparison as a single panel, drawing an
arrow from each sample's position before correction to its position after. Use
it when the question is which samples moved rather than whether the groups
separate.

```{r trajectory, fig.width=7, fig.height=5, out.width="80%"}
plotPCATrajectory(bv, assayBefore = "raw", assayAfter = "naive")
```

It inherits the caveat above: an arrow shows that a sample moved, not which of
the two confounded variables moved it.

## Whether the sample relationships survived

A correction can leave every variance fraction looking reasonable and still
rearrange which samples resemble which. Two diagnostics check that directly.

`plotSampleDistance()` draws the sample-sample distance matrix, so blocks of
mutually similar samples are visible as blocks:

```{r distance, fig.width=6.5, fig.height=5.5, out.width="70%"}
plotSampleDistance(bv, assayName = "raw")
```

`plotCorrelationChange()` reduces that to one number per assay: the mean
absolute change in sample-sample correlation relative to the baseline. Large
values mean the correction rewrote the relationships between samples, which is
a different failure from removing too much variance and is not visible in any
of the variance tables.

```{r corrchange, fig.width=7, fig.height=4.5, out.width="80%"}
plotCorrelationChange(ev)
```

Note which assay scores best here: `raw_noise`, at 0.100 against 0.187 for the
naive correction. The assay this vignette has used throughout as its example of
a destructive transformation is the one that disturbs sample relationships
least. That is not a defect in the diagnostic, it is what the diagnostic
measures - adding independent noise to every feature shrinks all sample
correlations towards zero fairly evenly, which is a small *change in* the
correlation structure, whereas shifting feature means per batch rearranges it
more. The two diagnostics agree that the protected run is gentler than the naive
one; where they part company is `raw_noise`, which this measure ranks first and
`basisRetention()` ranks last. That single disagreement is the clearest
illustration in this vignette of why no one diagnostic settles the question.

## Whether the batches actually mixed

`plotBatchEntropy()` measures local batch mixing: for each sample it looks at
its `k` nearest neighbours and computes the entropy of their batch labels.
High entropy means a sample's neighbourhood is batch-mixed; low entropy means
it still sits among its own batch. This is the neighbourhood-level question
the global variance fractions cannot answer.

It draws one assay at a time, so the comparison is made by drawing two:

```{r entropy, fig.width=7, fig.height=4, out.width="80%", fig.show="hold"}
plotBatchEntropy(bv, assayName = "raw", batchVar = "batch", k = 10)
plotBatchEntropy(bv, assayName = "naive", batchVar = "batch", k = 10)
```

The distribution moves right after correction, which is the batch effect
leaving the local structure.

Under partial confounding this should be read with care, and this design shows
why. With two batches the entropy is bounded above by `log(2)`, about 0.693;
after the naive correction the mean neighbourhood entropy is 0.65, close to
that ceiling. But batch and group overlap here, so neighbourhoods that are
perfectly batch-mixed are neighbourhoods in which the biological groups have
been mixed too. Entropy approaching its maximum is evidence of overcorrection
in this design rather than of a job well done - which is the same verdict the
`shared` term reached, arrived at independently.

---

# Why This Matters

Batch correction is frequently applied as a routine preprocessing step, yet it is often performed without systematic or quantitative validation. This can lead to suboptimal outcomes, where overcorrection removes meaningful biological signal or undercorrection leaves residual technical artefacts that confound downstream analyses. BatchVaria addresses this gap by providing a transparent and reproducible framework for evaluating correction strategies. By enabling direct comparison of variance components across multiple corrected assays, the package supports quantitative validation of correction outcomes and facilitates more robust, interpretable analytical pipelines.

---

# BatchVaria Workflow Summary

The typical BatchVaria workflow is:

1. Load expression data into a SummarizedExperiment
2. Apply batch correction methods with runCorrection(), naming any
   variables to preserve
3. Profile variance using profileVariance()
4. Summarise variance using varianceTable(method = ...), naming the engine
5. Convert fractions to absolute variance using assayVariance()
6. Compare corrections using delta variance and basisRetention()
7. Visualise variance composition

---

# Summary

BatchVaria provides a coherent and extensible framework for evaluating batch correction in high-dimensional expression data. By combining variance decomposition, structured tracking of results across multiple assays, and intuitive summary outputs, the package enables users to quantitatively assess how different correction strategies impact both unwanted technical variation and biological signal. The integration of variance composition tables, relative change metrics, and visual summaries supports transparent comparison of correction methods within a single workflow. Together, these features facilitate reproducible and interpretable evaluation of batch correction, addressing a key gap in many transcriptomic analysis pipelines.

---
# Session Information
```{r session}
sessionInfo()
```