Missing Data

In 1976, the statistician Donald B. Rubin published a short paper in Biometrika that changed how scientists think about incomplete data [1]. Before Rubin, the usual practice was pragmatic: if a measurement was missing, you either deleted the sample or plugged in the column mean and moved on. Rubin's contribution was to show that what you can do about missing data depends on why it is missing. His framework introduced the idea of data "missing at random" (MAR); together with the later labels Missing Completely At Random (MCAR) and Missing Not At Random (MNAR), it separates harmless gaps from dangerous ones. It makes it possible to reason about when a simple fix is safe and when it quietly biases every result that follows.

Rubin went on to develop multiple imputation: filling in each missing value several times, with plausible values drawn at random, analysing every completed dataset, and combining the results [2]. Together with Roderick J. A. Little, he wrote Statistical Analysis with Missing Data (first edition 1987, third edition 2019), which became the standard reference on the subject [3]. Its ideas now shape practice in fields from clinical trials to economics.

In chemometrics, missing data has always been a practical reality rather than a theoretical curiosity. Detectors saturate. Lamps fail in the middle of a scan. Concentrations fall below the limit of detection. Datasets from different instruments cover different wavelength ranges and leave rectangular holes when they are merged. The field did have an early advantage: the algorithm now known as NIPALS (nonlinear iterative partial least squares), described by Herman Wold in 1966 for computing principal components [4], can simply skip missing entries while it iterates. PCA-based workflows in chemometrics could therefore tolerate gaps in the data long before general-purpose solutions became routine elsewhere [5],[6],[10].

Missing values in analytical chemistry are not an abstract statistical concern. They come from concrete physical and experimental situations, and knowing the source is the first step towards handling them correctly.

Detector saturation. Every detector has a dynamic range. When a signal goes beyond its upper limit (for example, an absorbance above about 3 in UV-Vis spectroscopy), the reading is clipped or flagged as invalid. The true value exists, but we cannot see it. This is common for strong absorption bands in concentrated samples.

Failed measurements. Instruments misbehave. A lamp flickers, a cuvette leaks, a fibre-optic probe loses contact with the sample. The result is a spectrum with corrupted channels, or a missing time point in a series. In process analytical technology, where instruments run continuously for days or weeks, intermittent failures are routine.

Below the detection limit. When a concentration is too low to distinguish from noise, the result is often recorded as "< LOD" or left blank. This is censored data: we know the value lies somewhere between zero and the detection limit, but not where. Replacing it with zero, with the LOD, or with LOD/2 each introduces a different bias. A simple practical improvement over plain imputation is to impute and then cap the filled values at the LOD, since we know they cannot be larger. For environmental and clinical data, Helsel's book is the standard reference on handling such values [14].

Merged datasets. When spectra from two instruments are combined (say, one covering 1100 to 2500 nm and another covering 800 to 2200 nm), each instrument contributes blanks in the range it did not measure. The missing values form large, predictable blocks.

Lost samples. In designed experiments a vial breaks, a reaction fails or a sample degrades before measurement. A whole row of the data matrix is missing, which is a different situation from a few missing wavelengths inside an otherwise complete spectrum.

Rubin's framework classifies missingness by how it relates to the data [1],[3]. The class determines which methods are safe.

The probability that a value is missing has nothing to do with the data, observed or not. It is as if a coin were flipped for every measurement.

Chemistry example: a technician accidentally skips a few vials in a routine batch because of an interruption, not because of anything about those vials.

Consequence: MCAR is the most benign case. Even deleting incomplete samples gives unbiased results (you only lose precision), because the remaining samples are still a random subset of the whole.

How to check: compare the observed variables between samples with and without missing values. If they look alike, MCAR is plausible. Little's test formalises this comparison [7]. No test can tell MAR from MNAR, though, because the evidence would be in the values you did not observe. That distinction must come from knowing the instrument and the chemistry.

The probability of being missing depends on observed information, but not on the missing value itself. The gaps are not random, but you can explain them with data you already have.

Chemistry example: the expensive reference analysis (say, HPLC of the active ingredient) is skipped when the NIR prediction, which is always available, is comfortably high, and run only when that prediction comes close to the lower specification limit. Whether the reference value $y$ is missing depends on the NIR spectrum (observed), not on the reference value itself. The demo below uses the same mechanism: $y$ goes missing more often when $x$ is high.

Consequence: simply deleting incomplete samples now biases summaries such as the mean of $y$, because the low-potency samples are over-represented among the complete ones. (A regression of $y$ on the NIR data, fitted to the complete samples, is still fine, because the selection depends only on what that regression already uses.) Methods that use the observed information, such as regression-based imputation, multiple imputation or maximum likelihood, can correct for it. MAR is the assumption behind most modern imputation methods.

The probability of being missing depends on the missing value itself. This is the dangerous case: the value is missing because of what it would have been.

Chemistry example: concentrations below the detection limit are left blank. Whether a value is missing depends directly on it being low. Detector saturation is the mirror image: values go missing because they are high.

Consequence: no standard imputation method can fully correct MNAR without extra assumptions about how values go missing. You need to model that process, for example with a censored (Tobit-type) model: left-censored for values below the LOD (we know only that the true value is smaller), right-censored for saturated signals (we know only that it is larger). Treating MNAR data as if they were MAR introduces bias.

The demo below has two variables: an absorbance $x$ that is always measured, and a reference value $y$ that is sometimes missing. Because the data are simulated, we know the hidden values and can check what each imputation does to the mean, the standard deviation and the correlation.

Things to try:

1. Under MCAR with mean imputation, the mean of $y$ is about right, but the standard deviation shrinks and the correlation with $x$ drops: all the imputed points sit on one horizontal line.
2. Switch to MAR. Now the missing values are mostly at high $x$, so the observed mean of $y$ is too low, and mean imputation inherits that bias. Regression on x fixes the mean, because it uses the observed information that explains the gaps.
3. Look closely at Regression on x: every filled value sits exactly on a straight line, which no real measurement does. The spread of $y$ is still a bit too small. Regression + noise adds a random residual to each filled value, so the imputed points scatter like real data. Press New draw of the noise a few times: the numbers wobble from one draw to another. Averaged over many draws (the last row of numbers), the standard deviation comes back closer to the truth than with plain regression, at the price of a slightly weaker correlation, because the added noise knows nothing about $x$. Multiple imputation repeats the draw in the same way and combines the results.
4. Switch to MNAR. Every method now overestimates the mean, because the low values are exactly the ones that are missing. Regression on $x$ does much better than the plain mean, since $x$ still carries information about $y$, but it cannot remove the bias completely.

When the amount of missing data is small and the mechanism is MCAR or close to it, simple methods can be enough. They are also useful baselines against which to compare better methods.

Remove every sample (row) that has at least one missing value.

When it works: MCAR data with few incomplete samples. If 5 % of 200 spectra miss a few wavelengths, deleting those 10 spectra may be perfectly acceptable.

When it fails: if the missingness is MAR or MNAR, the remaining samples are no longer representative. And if many samples have at least one gap (typical after merging instruments), deletion can throw away most of the data.

Instead of discarding whole samples, use every pair of variables that is available when computing a statistic. The covariance between variables $j$ and $k$, for example, uses all samples in which both are observed.

When it works: computing correlation or covariance matrices when different variables are missing in different samples.

When it fails: each entry of the matrix is computed from a different subset of samples, so the result may not be a valid covariance matrix (it can have negative eigenvalues). That breaks PCA and other methods that decompose it. Use with caution.

Replace each missing value with the mean (or median) of the observed values of the same variable:

where $\text{obs}(j)$ is the set of samples in which variable $j$ was observed and $n_{\text{obs},j}$ is how many there are.

When it works: a few values scattered randomly through the matrix. If only a small percentage is missing, the distortion is usually minor.

When it fails: every imputed value is the same number, so the variance of the variable is artificially reduced and its correlations with other variables are weakened, exactly as you saw in the demo. In PCA, this shrinks scores towards the centre and distorts loadings.

When a substantial fraction is missing, or the mechanism is MAR, simple methods introduce unacceptable bias. The following methods use the structure of the data to make better guesses.

For a sample with a missing value, find the $k$ most similar samples that do have that value, and use a (weighted) average of theirs [8].

1. Measure similarity with a distance (usually Euclidean), computed only over the variables observed in both samples.
2. Find the $k$ nearest neighbours among the samples in which the missing variable was observed.
3. Average their values, often giving closer neighbours more weight:

In practice a small constant is added to the distance to avoid division by zero when a neighbour is identical. Plain, unweighted averaging (all $w_l=1$) is also common: it is the default in scikit-learn's KNNImputer (weights="uniform"), and the Python example below switches to the distance weights shown here.

Strengths: kNN uses local structure. A missing value in a high-absorbance sample is estimated from other high-absorbance samples, not from the global mean. It makes no linearity assumption.

Weaknesses: it needs a meaningful distance, which becomes unreliable when many variables are missing at once (the distance is computed from fewer and fewer values). It also scales poorly to very large datasets, because each incomplete sample must be compared with every candidate.

Choosing $k$: values between 3 and 10 are common. A small $k$ follows local structure but is sensitive to noise in individual neighbours; a large $k$ is smoother but may mix distinct groups of samples.

PCA describes the main patterns of co-variation in the data. If a spectrum is missing a few wavelengths, the wavelengths that are present already tell us where that spectrum sits in score space, and the PCA model can then predict what the missing wavelengths should have been. The algorithm alternates between filling the gaps and refitting the model, in the style of expectation-maximisation (EM) [6],[8]:

1. Initialise the missing entries with column means (or any reasonable guess).
2. Fit PCA with $A$ components to the centred, completed matrix: $\mathbf{X}\approx 1{\bar{\mathbf{x}}}^T+\mathbf{T}{\mathbf{P}}^T$.
3. Predict the missing entries from the model and replace only those entries. Observed values never change:

4. Check convergence. If the imputed values changed by less than a small threshold (for example $10^{-6}$ in relative terms), stop. Otherwise go back to step 2.

In practice the PCA is computed with a truncated singular value decomposition (SVD), which is numerically stable and available in every linear algebra library. That is why the same method is also called iterative SVD.

The key decision is the number of components $A$. Too few, and the model cannot reconstruct the missing values. Too many, and it starts modelling noise, which then leaks into the imputed values.

Rather than a single "best guess", multiple imputation generates $M$ completed datasets, each with different imputed values drawn from a plausible distribution [2]. You analyse each one separately and combine the results. For a quantity of interest $\theta$ (a mean, a regression coefficient, a prediction):

The uncertainty of the combined estimate has two parts: the average uncertainty within each completed dataset, $\bar{W}$, and the spread between the datasets, $B$, which reflects how much we do not know about the missing values:

Multiple imputation is the standard in clinical trials and the social sciences. It is less common in chemometrics, where iterative PCA and related methods are preferred because they exploit the strong structure of spectral data directly.

PCA-based imputation deserves special attention in chemometrics because spectroscopic data are usually low-rank: hundreds or thousands of wavelengths are driven by a handful of underlying chemical components. That redundancy is exactly what makes reconstruction possible [5],[9]. A systematic comparison of the PCA-based methods, NIPALS and iterative PCA among them, is given by Folch-Fortuny, Arteaga and Ferrer [12]; Arteaga and Ferrer discuss what each method implies when the PCA model is later used for process monitoring [13].

NIPALS computes principal components one at a time, alternating between two small regressions: scores from loadings, and loadings from scores [4]. Each of those regressions can simply leave out the missing entries.

For component $a$, NIPALS alternates:

where the sums run only over the observed entries of row $i$ or column $j$. No explicit imputation is needed: the scores and loadings are estimated from what is available, and the missing entries can afterwards be reconstructed as ${\hat{x}}_{ij}={\sum }_a{t}_{ia}{p}_{ja}$ (plus the column mean, if the data were centred). The properties and limits of this and related approaches were studied in detail by Grung and Manne, and by Nelson, Taylor and MacGregor [9],[5].

The demo below lets you compare the three families on a small set of simulated NIR spectra built from two chemical components. The heat map shows the whole data matrix with the missing values in red; click a row (or use the Sample shown slider) to see that spectrum. The three error values at the bottom are always computed, so you can compare methods directly.

Things to try:

1. With random dropouts, column-mean imputation is clearly the worst: it ignores everything we know about the sample. kNN does better, and iterative PCA with $A=2$ (the true number of chemical components) is best.
2. Set $A=1$, then $A=3$ or more. One component is not enough to describe two chemical components. With three or more, the extra components start fitting noise, and the error grows several-fold. The algorithm also converges very slowly in that case; the demo tells you when it stops after its maximum number of iterations. This is why the number of components must be chosen with care. Regularised iterative PCA, which shrinks the small components instead of keeping them fully, reduces this problem [11].
3. Switch to merged instruments. Whole blocks are missing. kNN and PCA still work as long as the other samples cover the missing region, which is why the overlap between instruments matters so much.
4. Switch to below detection limit (the New random gaps button disappears, because this pattern is not random at all). All methods now overestimate the missing values, because the smallest values have been removed systematically. Even at the default amount, a few wavelengths vanish from every spectrum (the demo counts them), and no method can recover them from the data: for those columns all three methods fall back on the same guess. This is MNAR in action: the observed data alone cannot tell the methods how low the missing values were. Now switch on Use the known detection limit. Capping the imputed values at the LOD uses what we know about the mechanism, and the error drops several-fold.

There is no universal threshold, but some rough guidelines help frame the decision. These thresholds are rules of thumb, not results.

Less than about 5 % missing, scattered randomly: almost any method works, and the choice rarely changes downstream results.

About 5 to 20 % missing: simple methods start to introduce noticeable bias. Prefer kNN or iterative PCA. The pattern matters more than the percentage: 15 % missing concentrated in a few variables is more problematic than 15 % spread evenly.

About 20 to 50 % missing: use iterative PCA or multiple imputation, and check how much your conclusions change with the method.

More than about 50 % missing: the imputed values are driven more by model assumptions than by data. Consider dropping the affected variables, or reformulating the question.

Before imputing, look at where the data are missing. A heat map of the missingness indicator matrix (observed or missing), like the one in the demo above, can reveal:

• Scattered pattern: random dropouts, consistent with MCAR. Safe for most methods.
• Block pattern: entire variable ranges or entire samples missing. Common when merging datasets. Iterative PCA handles this well if the blocks overlap enough to constrain the model.
• Monotone pattern: if you sort samples by completeness, the missing entries form a staircase. This appears, for example, when samples drop out of a time series. Specialised methods exist for it, such as sequential regression imputation, in which each variable is imputed from the more complete ones.
• Systematic pattern: missing values concentrate in specific variables for specific groups of samples. Strong evidence for MAR or MNAR. Investigate before imputing.

The examples simulate low-rank NIR spectra, remove 10 % of the values at random, and compare mean imputation, kNN (Python only) and iterative PCA. With the Python version you should see an RMSE of about 0.17 for mean imputation, 0.05 for kNN and 0.011 for iterative PCA. The last one is close to the noise level of 0.01 that was put into the data, which is as good as any method can do. The MATLAB and R versions use different random numbers, so their values differ slightly.

import numpy as np
from sklearn.impute import KNNImputer
import matplotlib.pyplot as plt

def iterative_pca_imputation(X, n_components=2, max_iter=100, tol=1e-6):
    """
    Impute missing values using iterative PCA (EM-PCA).

    Parameters:
    -----------
    X : array-like, shape (n_samples, n_variables)
        Data matrix with np.nan for missing values
    n_components : int
        Number of PCA components for reconstruction
    max_iter : int
        Maximum number of iterations
    tol : float
        Convergence tolerance (relative change in imputed values)

    Returns:
    --------
    X_filled : array
        Completed data matrix with imputed values
    n_iter : int
        Number of iterations until convergence
    """
    X_filled = X.copy()
    missing = np.isnan(X)

    # Initialize missing values with column means
    col_means = np.nanmean(X, axis=0)
    for j in range(X.shape[1]):
        X_filled[missing[:, j], j] = col_means[j]

    for iteration in range(max_iter):
        old_imputed = X_filled[missing].copy()

        # Center the data
        means = X_filled.mean(axis=0)
        X_centered = X_filled - means

        # Truncated SVD (PCA)
        U, s, Vt = np.linalg.svd(X_centered, full_matrices=False)
        U = U[:, :n_components]
        s = s[:n_components]
        Vt = Vt[:n_components, :]

        # Reconstruct
        X_reconstructed = U @ np.diag(s) @ Vt + means

        # Replace only missing values
        X_filled[missing] = X_reconstructed[missing]

        # Check convergence
        new_imputed = X_filled[missing]
        change = np.sqrt(np.sum((new_imputed - old_imputed)**2))
        scale = np.sqrt(np.sum(new_imputed**2)) + 1e-10
        if change / scale < tol:
            return X_filled, iteration + 1

    return X_filled, max_iter


# Example: NIR spectra with missing values
rng = np.random.default_rng(42)
n_samples = 50
n_wavelengths = 200
wavelengths = np.linspace(1000, 2500, n_wavelengths)

# Simulate low-rank spectral data (3 chemical components)
components = np.array([
    np.exp(-((wavelengths - 1400)**2) / 20000),
    np.exp(-((wavelengths - 1800)**2) / 30000),
    np.exp(-((wavelengths - 2100)**2) / 15000)
])
concentrations = rng.random((n_samples, 3))
X_true = concentrations @ components + rng.normal(0, 0.01, (n_samples, n_wavelengths))

# Introduce 10% missing values (MCAR)
X_missing = X_true.copy()
mask = rng.random(X_true.shape) < 0.10
X_missing[mask] = np.nan

print(f"Missing: {mask.sum()} values ({100 * mask.mean():.1f}%)")

# Method 1: Mean imputation
X_mean = X_missing.copy()
col_means = np.nanmean(X_missing, axis=0)
for j in range(n_wavelengths):
    X_mean[np.isnan(X_mean[:, j]), j] = col_means[j]

# Method 2: kNN imputation, weighted by 1/distance (the default is "uniform")
knn = KNNImputer(n_neighbors=5, weights="distance")
X_knn = knn.fit_transform(X_missing)

# Method 3: Iterative PCA
X_pca, n_iter = iterative_pca_imputation(X_missing, n_components=3)
print(f"PCA converged in {n_iter} iterations")

# Compare reconstruction errors (only at missing positions)
rmse_mean = np.sqrt(np.mean((X_mean[mask] - X_true[mask])**2))
rmse_knn = np.sqrt(np.mean((X_knn[mask] - X_true[mask])**2))
rmse_pca = np.sqrt(np.mean((X_pca[mask] - X_true[mask])**2))

print(f"RMSE (mean imputation): {rmse_mean:.6f}")
print(f"RMSE (kNN, k=5):        {rmse_knn:.6f}")
print(f"RMSE (iterative PCA):   {rmse_pca:.6f}")

# Visualize one spectrum
idx = 0  # First sample
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

axes[0].plot(wavelengths, X_true[idx], 'k-', linewidth=2, label='True')
axes[0].plot(wavelengths[mask[idx]], X_true[idx, mask[idx]], 'rx',
             markersize=8, label='Missing positions')
axes[0].set_title('Original spectrum with missing value locations')
axes[0].set_xlabel('Wavelength (nm)')
axes[0].set_ylabel('Absorbance')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(wavelengths, X_true[idx], 'k-', linewidth=2, alpha=0.3, label='True')
axes[1].plot(wavelengths, X_mean[idx], '--', label=f'Mean (RMSE={rmse_mean:.4f})')
axes[1].plot(wavelengths, X_knn[idx], '--', label=f'kNN (RMSE={rmse_knn:.4f})')
axes[1].plot(wavelengths, X_pca[idx], '--', label=f'PCA (RMSE={rmse_pca:.4f})')
axes[1].set_title('Imputation comparison')
axes[1].set_xlabel('Wavelength (nm)')
axes[1].set_ylabel('Absorbance')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

% Example: NIR spectra with missing values
rng(42);
n_samples = 50;
n_wavelengths = 200;
wavelengths = linspace(1000, 2500, n_wavelengths);

% Simulate low-rank data (3 components)
comp1 = exp(-((wavelengths - 1400).^2) / 20000);
comp2 = exp(-((wavelengths - 1800).^2) / 30000);
comp3 = exp(-((wavelengths - 2100).^2) / 15000);
components = [comp1; comp2; comp3];

concentrations = rand(n_samples, 3);
X_true = concentrations * components + randn(n_samples, n_wavelengths) * 0.01;

% Introduce 10% MCAR missing values
mask = rand(size(X_true)) < 0.10;
X_missing = X_true;
X_missing(mask) = NaN;

fprintf('Missing: %d values (%.1f%%)\n', sum(mask(:)), 100 * mean(mask(:)));

% Mean imputation
X_mean = X_missing;
col_means = mean(X_missing, 1, 'omitnan');
for j = 1:n_wavelengths
    X_mean(isnan(X_mean(:, j)), j) = col_means(j);
end

% Iterative PCA
[X_pca, n_iter] = iterative_pca_imputation(X_missing, 3);
fprintf('PCA converged in %d iterations\n', n_iter);

% Compare errors
rmse_mean = sqrt(mean((X_mean(mask) - X_true(mask)).^2));
rmse_pca = sqrt(mean((X_pca(mask) - X_true(mask)).^2));

fprintf('RMSE (mean imputation): %.6f\n', rmse_mean);
fprintf('RMSE (iterative PCA):   %.6f\n', rmse_pca);

% Plot
figure('Position', [100 100 1000 600]);
subplot(2, 1, 1);
plot(wavelengths, X_true(1, :), 'k-', 'LineWidth', 2);
hold on;
idx_miss = find(mask(1, :));
plot(wavelengths(idx_miss), X_true(1, idx_miss), 'rx', 'MarkerSize', 8);
title('Original spectrum with missing positions');
xlabel('Wavelength (nm)');
ylabel('Absorbance');
legend('True', 'Missing');
grid on;

subplot(2, 1, 2);
plot(wavelengths, X_true(1, :), '-', 'LineWidth', 2, 'Color', [0.7 0.7 0.7]);
hold on;
plot(wavelengths, X_mean(1, :), '--', 'LineWidth', 1.5);
plot(wavelengths, X_pca(1, :), '--', 'LineWidth', 1.5);
title('Imputation comparison');
xlabel('Wavelength (nm)');
ylabel('Absorbance');
legend('True', 'Mean', 'PCA');
grid on;

function [X_filled, n_iter] = iterative_pca_imputation(X, n_components, max_iter, tol)
    % Impute missing values using iterative PCA (EM-PCA)
    %
    % Parameters:
    %   X            - Data matrix with NaN for missing values
    %   n_components - Number of PCA components
    %   max_iter     - Maximum iterations (default: 100)
    %   tol          - Convergence tolerance (default: 1e-6)
    %
    % Returns:
    %   X_filled - Completed matrix
    %   n_iter   - Iterations until convergence

    if nargin < 3, max_iter = 100; end
    if nargin < 4, tol = 1e-6; end

    X_filled = X;
    missing = isnan(X);

    % Initialize with column means
    col_means = mean(X, 1, 'omitnan');
    for j = 1:size(X, 2)
        X_filled(missing(:, j), j) = col_means(j);
    end

    for iter = 1:max_iter
        old_imputed = X_filled(missing);

        % Center
        means = mean(X_filled, 1);
        X_centered = X_filled - means;

        % Truncated SVD
        [U, S, V] = svds(X_centered, n_components);

        % Reconstruct
        X_reconstructed = U * S * V' + means;

        % Replace only missing values
        X_filled(missing) = X_reconstructed(missing);

        % Check convergence
        new_imputed = X_filled(missing);
        change = norm(new_imputed - old_imputed);
        scale = norm(new_imputed) + 1e-10;
        if change / scale < tol
            n_iter = iter;
            return;
        end
    end
    n_iter = max_iter;
end

iterative_pca_imputation <- function(X, n_components = 2,
                                      max_iter = 100, tol = 1e-6) {
  # Impute missing values using iterative PCA (EM-PCA)
  #
  # Parameters:
  #   X            - Data matrix with NA for missing values
  #   n_components - Number of PCA components
  #   max_iter     - Maximum iterations
  #   tol          - Convergence tolerance
  #
  # Returns:
  #   list with X_filled (completed matrix) and n_iter

  missing <- is.na(X)
  X_filled <- X

  # Initialize with column means
  col_means <- colMeans(X, na.rm = TRUE)
  for (j in seq_len(ncol(X))) {
    X_filled[missing[, j], j] <- col_means[j]
  }

  for (iter in seq_len(max_iter)) {
    old_imputed <- X_filled[missing]

    # Center
    means <- colMeans(X_filled)
    X_centered <- sweep(X_filled, 2, means)

    # Truncated SVD
    sv <- svd(X_centered, nu = n_components, nv = n_components)
    U <- sv$u[, 1:n_components, drop = FALSE]
    d <- sv$d[1:n_components]
    V <- sv$v[, 1:n_components, drop = FALSE]

    # Reconstruct
    X_reconstructed <- U %*% diag(d, nrow = n_components) %*% t(V)
    X_reconstructed <- sweep(X_reconstructed, 2, means, "+")

    # Replace only missing values
    X_filled[missing] <- X_reconstructed[missing]

    # Convergence check
    new_imputed <- X_filled[missing]
    change <- sqrt(sum((new_imputed - old_imputed)^2))
    scale <- sqrt(sum(new_imputed^2)) + 1e-10
    if (change / scale < tol) {
      return(list(X_filled = X_filled, n_iter = iter))
    }
  }

  return(list(X_filled = X_filled, n_iter = max_iter))
}

# Example: NIR spectra with missing values
set.seed(42)
n_samples <- 50
n_wavelengths <- 200
wavelengths <- seq(1000, 2500, length.out = n_wavelengths)

# Simulate low-rank data (3 components)
comp1 <- exp(-((wavelengths - 1400)^2) / 20000)
comp2 <- exp(-((wavelengths - 1800)^2) / 30000)
comp3 <- exp(-((wavelengths - 2100)^2) / 15000)
components <- rbind(comp1, comp2, comp3)

concentrations <- matrix(runif(n_samples * 3), ncol = 3)
X_true <- concentrations %*% components +
          matrix(rnorm(n_samples * n_wavelengths, sd = 0.01),
                 nrow = n_samples)

# Introduce 10% MCAR missing
mask <- matrix(runif(n_samples * n_wavelengths) < 0.10,
               nrow = n_samples)
X_missing <- X_true
X_missing[mask] <- NA

cat(sprintf("Missing: %d values (%.1f%%)\n",
            sum(mask), 100 * mean(mask)))

# Mean imputation
X_mean <- X_missing
col_means <- colMeans(X_missing, na.rm = TRUE)
for (j in seq_len(n_wavelengths)) {
  X_mean[is.na(X_mean[, j]), j] <- col_means[j]
}

# Iterative PCA
result <- iterative_pca_imputation(X_missing, n_components = 3)
X_pca <- result$X_filled
cat(sprintf("PCA converged in %d iterations\n", result$n_iter))

# Compare errors
rmse_mean <- sqrt(mean((X_mean[mask] - X_true[mask])^2))
rmse_pca <- sqrt(mean((X_pca[mask] - X_true[mask])^2))

cat(sprintf("RMSE (mean imputation): %.6f\n", rmse_mean))
cat(sprintf("RMSE (iterative PCA):   %.6f\n", rmse_pca))

# Plot
par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))

plot(wavelengths, X_true[1, ], type = "l", lwd = 2,
     xlab = "Wavelength (nm)", ylab = "Absorbance",
     main = "Original spectrum with missing positions")
miss_idx <- which(mask[1, ])
points(wavelengths[miss_idx], X_true[1, miss_idx],
       col = "red", pch = 4, cex = 1.2)
legend("topright", legend = c("True", "Missing"),
       lty = c(1, NA), pch = c(NA, 4),
       col = c("black", "red"), lwd = c(2, NA))
grid()

plot(wavelengths, X_true[1, ], type = "l", lwd = 2,
     col = rgb(0, 0, 0, 0.3),
     xlab = "Wavelength (nm)", ylab = "Absorbance",
     main = "Imputation comparison")
lines(wavelengths, X_mean[1, ], col = "#2563eb", lwd = 1.5, lty = 2)
lines(wavelengths, X_pca[1, ], col = "#059669", lwd = 1.5, lty = 2)
legend("topright",
       legend = c("True", "Mean", "PCA"),
       col = c(rgb(0, 0, 0, 0.3), "#2563eb", "#059669"),
       lwd = c(2, 1.5, 1.5), lty = c(1, 2, 2))
grid()

Always investigate the mechanism. Before choosing a method, find out why data are missing. Plot the missingness pattern and check whether it relates to observed variables. MCAR, MAR or MNAR decides which methods are valid.

Do not just delete. Listwise deletion is tempting but wasteful: it throws away entire spectra because of a few missing wavelengths, and it introduces bias whenever the data are not MCAR.

Use the structure of spectra. Spectroscopic data are low-rank. PCA-based imputation exploits that and usually beats generic methods such as mean imputation on spectral data.

Validate the imputation. Hide known values, impute them, and compare with the truth. Report the error. If downstream results (PCA scores, PLS predictions) change a lot when you change the imputation method, your conclusions depend on assumptions rather than data.

Missing values are never just "gaps to fill". Rubin's three mechanisms tell you how dangerous they are: MCAR loses information but not accuracy, MAR can be corrected with the data you have, and MNAR needs a model of why values went missing. MAR and MNAR cannot be told apart from the data alone: that judgement comes from knowing the instrument and the chemistry.

Mean imputation is easy but shrinks variance and weakens correlations. Methods that use the structure of the data, such as kNN and especially iterative PCA for low-rank spectra, do much better. Whatever you choose, look at the pattern first and validate the result on values you deliberately hid.