---
title: "Missing values: what they mean and how to fill them"
author: "Sergio Ciordia"
package: NADIA
output:
BiocStyle::html_document:
toc: true
toc_float: true
number_sections: true
vignette: >
%\VignetteIndexEntry{Missing values: what they mean and how to fill them}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", crop = NULL)
options(width = 100)
```
# Introduction
Missing values in a proteomics matrix are not interchangeable. A value may be
absent because a protein was present but not quantified successfully, or because
its abundance was too low to be detected. These mechanisms call for different
assumptions and can coexist within the same experiment.
This vignette explains how to inspect that missingness, summarises the 20
strategies available in NADIA and describes the behaviour of the hybrid and
detection-aware options. It does not repeat the benchmarking metrics used to
rank candidate methods; those are covered in `vignette("choosing-methods")`.
# Installation
```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("NADIA")
```
```{r library, message = FALSE}
library(NADIA)
data(nadia_dia)
```
# Why values are missing
Imputation replaces missing entries with estimates so that analyses requiring a
complete matrix can be performed. The estimate should reflect a plausible reason
for the value being absent. NADIA uses the following practical distinction:
- **MAR (missing at random):** the protein is present, but its intensity was not
quantified successfully—for example because of poor chromatography, precursor
interference or a missed match. Under this assumption, the gap is not itself
evidence that the unknown intensity was low. Its value can therefore be
estimated from patterns in the available measurements, such as similar
protein profiles or related samples.
- **MNAR (missing not at random):** the protein is absent or its abundance is
below the detection limit. Missingness depends on the unobserved abundance
being low and therefore carries biological or measurement information.
Replacing such a gap with a typical neighbouring value may erase that
evidence, so a low-value or detection-aware method is generally more
plausible.
```{r mar-mnar-diagram, echo = FALSE, fig.align = "center", out.width = "100%", fig.cap = "Practical distinction between MAR and MNAR missing values in proteomics. MAR values are reconstructed from observed data structure, whereas MNAR values call for a low-value or detection-aware model.", fig.alt = "Two-panel diagram comparing MAR and MNAR missing values. MAR shows a missing replicate among similar observed intensities and recommends structure-based imputation. MNAR shows low values falling below a detection limit and recommends low-value or detection-aware imputation. A note states that an individual missing value does not reveal its cause and that both mechanisms can coexist."}
knitr::include_graphics("figures/MAR_MNAR_Missingness.svg")
```
Click the figure to enlarge it; press Esc or click outside the image to close it.
An individual `NA` does not reveal why the value is missing. MAR and MNAR are
possible explanations that must be evaluated using the experimental design, the
overall distribution of missing values and their relationship with observed
intensity. Both mechanisms may occur in the same dataset.
## How NADIA represents missing intensities
By the time the data enter the analysis, an unreported intensity is represented
as `NA`. How that `NA` is created depends on the input format. In a long
Spectronaut report, a protein that was not quantified in a particular run is
usually absent from the report; when NADIA constructs the wide `protein_quant`
table, the corresponding protein–run cell becomes `NA`. Some wide reports use
zero instead to mark an unreported intensity. Before normalisation,
`normalize_proteomics()` converts these zeros to `NA`, so downstream functions
handle both situations consistently.
# Diagnosing missingness in the example data
The `nadia_dia` example dataset included with NADIA retains the missing-value
pattern of the original experiment. Approximately 8% of its intensities are
missing, and more than a quarter of its proteins contain at least one gap.
```{r pattern}
quant <- as.matrix(nadia_dia$protein_quant[
grep("^PG.Quantity_", colnames(nadia_dia$protein_quant))])
round(100 * mean(is.na(quant)), 1)
table(missing_per_protein = rowSums(is.na(quant)))
```
The distribution of `NA` values across proteins is informative. There are 70
proteins with one gap, 67 with two and 99 with three, followed by 243 proteins
with exactly four missing values. Four is also the number of replicates per
condition. Most of those 243 proteins are absent from every replicate of one
condition:
```{r condition-dropout}
condition <- rep(c("A", "B", "D"), each = 4)
whole_condition_missing <- apply(
quant, 1,
function(x) any(tapply(is.na(x), condition, all)))
c(exactly_four_missing = sum(rowSums(is.na(quant)) == 4),
whole_condition_missing = sum(rowSums(is.na(quant)) == 4 &
whole_condition_missing))
```
Because each condition contains four replicates, a protein with four missing
values may be absent from one complete condition. This pattern suggests that at
least some gaps reflect condition-wise dropout rather than the scattered
missingness usually treated as MAR in this workflow. As a complementary
diagnostic, we next examine the distribution of the intensities that remain
observed:
```{r intensity-vs-missing, fig.width = 6, fig.height = 4, fig.alt = "Mean log2 intensity against number of missing values"}
n_missing <- rowSums(is.na(quant))
mean_log2 <- rowMeans(log2(quant), na.rm = TRUE)
boxplot(mean_log2 ~ n_missing,
xlab = "missing values in the protein",
ylab = "mean log2 intensity of the values present",
main = "Missingness is concentrated at low abundance")
```
Proteins with more missing values tend to have lower observed intensities. This
pattern supports the presence of an abundance-dependent MNAR component, in
which low-abundance proteins are more likely to fall below the detection limit.
However, this overall trend cannot determine the mechanism responsible for each
individual missing value. Taken together, the condition-wise dropout pattern
and its relationship with intensity suggest that this dataset contains both
scattered gaps and missing values caused by low abundance.
## Data processing shapes the missing-value pattern
NADIA includes DIA-NN and Spectronaut reports generated from the same twelve
LC-MS injections. The instrument data, samples and raw files are therefore the
same, allowing a descriptive comparison of the missing-value patterns produced
by the two processing workflows.
```{r engines, include = FALSE}
diann <- preprocess_diann(
system.file("extdata", "nadia_diann_report.tsv.gz", package = "NADIA"),
condition_order = c("A", "B", "D"), verbose = FALSE)
qd <- as.matrix(diann$protein_quant[
grep("^PG.Quantity_", colnames(diann$protein_quant))])
shared_protein_groups <- length(intersect(
nadia_dia$protein_quant$PG.ProteinGroups,
diann$protein_quant$PG.ProteinGroups))
```
Each trimmed file contains 2,000 protein groups selected independently, and
only `r shared_protein_groups` groups are shared. The comparison can therefore
describe overall patterns, but it is not a protein-by-protein benchmark of the
two search engines.
```{r engines-summary}
summarise_missingness <- function(q) {
n_missing <- rowSums(is.na(q))
grp <- rep(c("A", "B", "D"), each = 4)
mean_log2 <- rowMeans(log2(q), na.rm = TRUE)
median_complete <- median(mean_log2[n_missing == 0])
median_incomplete <- median(mean_log2[n_missing > 0], na.rm = TRUE)
c(
missing_intensities_pct = round(100 * mean(is.na(q)), 1),
condition_dropout_pct = round(100 * mean(apply(
q, 1, function(r) any(tapply(is.na(r), grp, all)))), 1),
median_complete_log2 = round(median_complete, 1),
median_incomplete_log2 = round(median_incomplete, 1),
complete_minus_incomplete_log2 = round(
median_complete - median_incomplete, 1)
)
}
rbind(
spectronaut = summarise_missingness(quant),
diann = summarise_missingness(qd)
)
```
The distribution of missing values per protein provides additional context:
```{r engines-shape}
missing_per_protein <- function(q) {
table(factor(rowSums(is.na(q)), levels = 0:12))
}
rbind(
spectronaut = missing_per_protein(quant),
diann = missing_per_protein(qd)
)
```
The Spectronaut example contains 8.1% missing intensities, compared with 12.5%
in the DIA-NN example. In both matrices, the distribution peaks at four missing
values per protein—the number of replicates in one condition—indicating that
condition-wise dropout is common. This affects 14.4% of proteins in Spectronaut
and 20.4% in DIA-NN; DIA-NN also contains more proteins with eight or more
missing values. Incomplete proteins (`n_missing > 0`) have a lower median
observed intensity than complete proteins (`n_missing = 0`) in both matrices,
with a difference of 1.9 log~2~ units in
Spectronaut and 2.6 log~2~ units in DIA-NN. These within-matrix differences do
not compare the absolute intensity scales of the search engines. Overall, the
DIA-NN example shows more missingness, more condition-wise dropout and a
stronger association between low abundance and missingness. This does not show
that one search engine is better; it shows that imputation decisions should be
based on the missing-value pattern of the matrix that will actually be
analysed.
# Imputation strategies available in NADIA
NADIA offers 20 selectable imputation strategies. The table below groups them
according to their main purpose and the type of missingness they are generally
intended to address. This grouping provides practical guidance; it does not
imply that the mechanism behind every missing value can be determined
unambiguously.
::: {.wide-reference-table .imputation-methods-table}
| Group | imp_method | Principle |
|---|---|---|
| Structure-based, generally used for MAR | `bpca` | Bayesian principal-component imputation using `pcaMethods`. |
| | `knn` | Estimates a value from similar proteins using k-nearest neighbours. |
| | `mice` | Multiple imputation by chained equations; NADIA averages the completed matrices. |
| | `missForest` | Non-parametric random-forest imputation using relationships among variables. |
| | `Impseq` | Sequential model-based imputation from `rrcovNA`. |
| | `Impseqrob` | Robust sequential imputation that reduces sensitivity to outlying profiles. |
| | `MLE` | Maximum-likelihood imputation under a multivariate normal model. |
| | `nbavg` | Replaces a gap with the mean of the nearest observed positions in the same protein profile. |
| Low-value, generally used for MNAR | `QRILC` | Quantile-regression imputation for left-censored data. |
| | `MinDet` | Uses a low, sample-specific quantile; the default is the first percentile. |
| | `MinProb` | Draws low values around a sample-specific low quantile. |
| | `PI` | Perseus-style sampling from a narrowed, down-shifted normal distribution. |
| | `min` | Replaces every gap with the global observed minimum. |
| | `halfmin` | Uses half the global minimum on the original scale, equivalent to one unit below the minimum on the log2 scale; this convention is used by DIA-NN. |
| | `zero` | Replaces every gap with zero on the log2 scale. |
| | `with` | Replaces every gap with a user-supplied constant. |
| Hybrid | `combo` | Classifies missing cells from condition-wise presence patterns, then applies separate MAR and MNAR methods. Defaults are `Impseqrob` and `min`. |
| | `softHybrid` | Blends estimates from MAR and MNAR methods continuously using protein missingness and mean intensity. |
| Detection model | `limpa` | Estimates a detection-probability curve and returns abundance estimates with uncertainty for downstream analysis with `de_method = "limpa"`. |
| No imputation | `none` | Leaves missing values unchanged. |
:::
Structure-based methods estimate scattered gaps from relationships among the
observed proteins and samples. In contrast, low-value methods assume that
missingness is associated with low abundance and impute values from the lower
tail of the intensity distribution. `combo`, `softHybrid` and `limpa` address
this distinction more explicitly, as described below. The `none` option leaves
the missing values unchanged and therefore requires downstream methods that can
analyse an incomplete matrix.
`halfmin` operates on log~2~-transformed data. Halving the global minimum on the
original intensity scale is equivalent to subtracting one from its log~2~ value,
so the imputed value is `log2(minimum) - 1`. Both `min` and `halfmin`
depend on the single lowest observed intensity in the matrix. Distribution-based
methods such as `MinDet`, `QRILC` and `MinProb` instead characterise the
low-intensity tail using more than this single value.
# Why one method may not represent every gap
If the imputation method does not match the likely cause of missingness, it may
distort the biological signal. When a protein is not observed in any replicate
of one condition, a MAR method may borrow information from the other conditions
and impute values that are too similar to them. This can make a genuine
presence–absence difference appear smaller. Conversely, when an isolated value
is missing because of a technical measurement failure, an MNAR method may place
it unnecessarily close to the detection limit and create an artificial
difference between conditions.
Because both patterns can coexist, applying the same method to every `NA` may
not be appropriate. NADIA's hybrid strategies address this problem by applying
different imputation assumptions to different gaps. However, a hybrid method is
not automatically the best choice for every dataset: when one missingness
mechanism predominates, a single method may be sufficient. Candidate methods
can be compared using the procedure described in
`vignette("choosing-methods")`.
# The default hybrid: `combo`
The `combo` option is NADIA's default because a single dataset may contain both
isolated missing measurements and proteins that are not quantified in most or
all replicates of a condition. For each protein, NADIA examines how its observed
and missing values are distributed across the experimental conditions. Isolated
gaps are assigned to the MAR branch, whereas patterns consistent with
condition-wise dropout are assigned to the MNAR branch. By default, these
branches use `Impseqrob` and `min`, respectively.
Four arguments control the classification:
- `prop_na_mnar = 0.51` defines the minimum proportion of missing values within
a condition for its gaps to be considered MNAR candidates.
- `prop_present_mar = 0.5` requires at least half of the replicates in another
condition to contain an observed intensity for that protein. In this example,
this means at least two of the four replicates.
- `min_present_mar = 1` requires at least one observed value in that other
condition, regardless of its number of replicates. The proportional and
absolute presence requirements must both be met.
- `require_n_conditions = 1` defines how many other conditions must meet both
presence requirements. With the default value, one qualifying condition is
sufficient.
A missing cell is assigned to the MNAR branch only when two criteria are met:
its condition reaches the required proportion of missing values, and the
required number of other conditions provide sufficient observed evidence for
that protein. All other missing cells are assigned to the MAR branch. With four
replicates, `prop_na_mnar = 0.51` requires at least three missing values in a
condition; two missing values represent only 50% and do not meet the threshold.
The following examples use conditions A and B, each with four replicates.
`Observed` denotes a measured intensity. Classification is performed separately
for the missing cells in each condition, so the protein is not assigned one
fixed MAR or MNAR label.
```{r combo-classification-diagram, echo = FALSE, fig.align = "center", out.width = "100%", fig.cap = "Examples of how the four combo arguments classify missing cells. A cell is assigned to the MNAR branch only when the missingness threshold and the evidence in other conditions are both satisfied; otherwise it enters the MAR branch or the protein may be removed by the preliminary combo filter.", fig.alt = "Four coloured examples of combo classification with four replicates per condition. The first is assigned to MNAR because one condition is 75 percent missing and another is fully observed. The second is assigned to MAR because 50 percent missing is below the 51 percent threshold. The third protein is removed because no other condition has at least 50 percent observed values. The fourth uses three conditions and shows that one qualifying supporting condition is sufficient."}
knitr::include_graphics("figures/Combo_Classification_Examples.svg")
```
Click the figure to enlarge it; press Esc or click outside the image to close it.
In these examples, each condition has four replicates. Therefore,
`prop_present_mar = 0.5` means that another condition must contain at least two
observed values before it can support an MNAR classification. Once two values
are observed, the requirement imposed by `min_present_mar = 1` is already
satisfied, so this second parameter does not change the classification here.
Before filling any missing values, `combo` also removes proteins that contain
too little observed information to apply these rules. A protein is retained if
it either contains a pattern classified as an MNAR candidate or has sufficient
observed values in at least `require_n_conditions` conditions.
The following example reports how many proteins pass this filter and how many
of the original missing cells are assigned to each branch:
```{r combo, message = FALSE}
res <- process_proteomics(nadia_dia, verbose = FALSE)
imp <- impute_proteomics(res$se_proc, "cycloess", verbose = FALSE)
x_combo <- SummarizedExperiment::assay(res$se_proc, "cycloess")
c(
proteins_before = nrow(nadia_dia$protein_quant),
proteins_after_combo_rules = nrow(res$se_proc),
MAR_missing_values = sum(imp$mar_mask & is.na(x_combo)),
MNAR_missing_values = sum(imp$mnar_mask & is.na(x_combo))
)
```
The output shows that both branches are used. These assignments reflect the
selected rules and do not establish the true mechanism of any individual `NA`.
When a distribution-based MNAR method is selected, its low-intensity parameters
are estimated from the original matrix, before MAR imputation, so reconstructed
MAR values do not shift the estimated lower tail.
# `softHybrid`: a continuous combination
`softHybrid` was introduced by [Shi et al. (2026)](https://doi.org/10.64898/2026.01.13.699212)
as a hybrid imputation method for single-cell proteomics. Unlike `combo`, which
assigns each missing cell to either a MAR or an MNAR branch, `softHybrid`
calculates both estimates and combines them continuously.
For each protein, the method uses its proportion of missing values and its mean
observed intensity to calculate a weight. Extensive missingness and low
intensity increase the contribution of the MNAR estimate, whereas a
well-observed, higher-intensity protein receives more influence from the MAR
estimate. The final value for each gap is therefore a weighted average of the
two estimates rather than the result of a binary classification. NADIA uses
`Impseqrob` and `min` as the default MAR and MNAR methods, respectively.
The transition is controlled by two sigmoid functions. Parameter `a` determines
how rapidly the weight changes with the missing-value proportion, `b` controls
the transition along the standardised intensity axis and `lambda` balances the
two signals. In most cases, the defaults provide the natural starting point:
```{r softhybrid}
imp_soft <- impute_proteomics(
res$se_proc,
normalized_assay_name = "cycloess",
imp_method = "softHybrid",
mar_method = "Impseqrob",
mnar_method = "min",
verbose = FALSE)
grep("softHybrid", SummarizedExperiment::assayNames(imp_soft$se), value = TRUE)
imp_soft$imputation_summary
```
The summary reports the initial and final missing-value rates and the
distribution of the MAR weights. Whether this continuous combination improves
the analysis remains dataset-dependent and should be assessed as described in
`vignette("choosing-methods")`.
# `limpa`: detection-aware estimates with uncertainty
The `limpa` framework was described by [Li, Cobbold and Smyth
(2025)](https://doi.org/10.1101/2025.04.28.651125). It differs from conventional
methods that simply replace each missing value with a single completed-matrix
estimate. First, `limpa` estimates a detection-probability curve that relates
protein abundance to the probability of observing a signal. It then uses this
curve to estimate the abundance of missing measurements together with their
standard errors.
Here, *precision* means statistical reliability. `limpa` converts each
abundance estimate and its uncertainty into a **precision weight** for differential
analysis: an estimate based on limited evidence receives less weight than an
observed, reproducible measurement. Missing values can therefore still
contribute information, while uncertain estimates have less influence on the
fitted model.
The estimates and weights are stored together and must remain linked through the
differential abundance step:
```{r limpa, eval = requireNamespace("limpa", quietly = TRUE)}
res_limpa <- process_proteomics(nadia_dia,
imp_method = "limpa",
de_method = "limpa",
verbose = FALSE)
class(S4Vectors::metadata(res_limpa$se_proc)$limpa_elist)
table(res_limpa$DEPs_results$Change)
```
`imp_method = "limpa"` should therefore be paired with
`de_method = "limpa"`. Using its abundance estimates with an ordinary `limma`
test would discard the uncertainty model.
# Choosing among the strategies
The missing-value pattern, biological design and intended differential model
should guide the set of plausible strategies. When the true missing abundances
are unknown, no diagnostic can establish a universally correct method. NADIA
therefore provides artificial-masking metrics, ranking summaries and plots for
comparing candidate methods, while sensitivity analysis shows whether the main
biological conclusions depend on that choice.
Those evaluation procedures, their limitations and the interpretation of the
rankings are described in `vignette("choosing-methods")`. Keeping them there
separates two questions: this vignette explains **what the missing values may
mean and how NADIA can handle them**; the companion vignette explains **how to
compare the available choices on a particular dataset**.
# References and software
- Li M, Cobbold SA, Smyth GK (2025). [Quantification and differential analysis
of mass spectrometry proteomics data with probabilistic recovery of
information from missing values](https://doi.org/10.1101/2025.04.28.651125).
*bioRxiv*, 2025.04.28.651125. This work describes the statistical framework
implemented in `limpa`; the software is available from
[Bioconductor](https://bioconductor.org/packages/limpa) and
[GitHub](https://github.com/SmythLab/limpa).
- Shi Y, Davis S, Charles PD, Taylor S, Dombi E, Berridge G, Ebner D, Fischer R
(2026). [SoftHybrid: A Hybrid Imputation Algorithm Optimised for Single-Cell
Proteomics Data](https://doi.org/10.64898/2026.01.13.699212). *bioRxiv*,
2026.01.13.699212. This work describes the `softHybrid` imputation strategy.
# Session information
```{r session}
sessionInfo()
```