Neel Mishra
← All Posts
Deep Learning · Data Preparation and Preprocessing · Handling Missing Data

Iterative (Model-Based) Imputation

September 2026 · 11 min read

Iterative imputation predicts an incomplete feature from the other features, then repeats that process across columns. Its central idea is simple: an updated fill in column B may make column B a more useful predictor when we next fill column A.

The method needs a starting guess, a record of which entries were originally missing, and a prediction model for each column. This chapter follows those three pieces through one small table. Multiple imputation is a separate next step when we also need to represent uncertainty about the fills.

One column becomes the target; the others become predictors

Suppose two numerical features, A and B, tend to move together. For the example below, start every gap at its column’s observed mean: A gets $(1+2+4)/3\approx2.333$, and B gets $(2+6+8)/3\approx5.333$.

To update B, fit a regression using rows whose original B value was observed. A supplies the predictor, using its current fill wherever needed. Replace only B’s missing entries with predictions. Then update A using the newly updated B. A complete B-then-A sweep is one round.

Update one column at a time
RowAB
112
225.3333
32.33336
448
Initialization: dashed cells contain column means. All six observed values remain fixed.
Dashed cells are estimates throughout, even after they stop changing. Step through four rounds of ordinary least-squares updates. These are two feature columns; B is not a downstream prediction label.

On the first B update, rows 1, 3, and 4 train the regression. Row 3 participates because its B is known, even though its A is currently an estimate. The fitted line is approximately B = 0.5246 + 1.9672 × A, so row 2’s B becomes 4.4590. The next A update moves row 3’s A to 2.9289.

This is how information passes between columns. We never promote a previously imputed target into an observed training target. Keep the original missingness mask even after every table cell contains a number.

The same update in compact notation

For row $r$ and target column $c$, let $x_{r,-c}$ denote its other feature values, including their current fills. Fit a model $f_c$ using only rows where column $c$ was observed. Its prediction supplies an updated missing entry:

$$\hat x_{rc}=f_c(x_{r,-c}).$$

The hat marks an estimate. The notation hides the main implementation detail: each later column uses the fills currently available, including updates earlier in the same round.

Stable fills are not proof of correct fills

After four rounds in this illustration, the two fills are approximately 4.000281 and 2.999960. The table was deliberately constructed around the relationship B = 2 × A. Real measurements may have noise, nonlinear relationships, or missingness patterns that leave little evidence about a gap.

A stopping rule asks whether the working table changed by less than a chosen tolerance, or whether a round limit was reached. A wrong model can settle on stable but wrong predictions. More rounds do not repair a model that excludes an important predictor or extrapolates beyond the observed data.

With stochastic imputation, successive fills are intended to vary. Do not apply “the numbers stopped changing” as the standard for such a chain. The multiple-imputation chapter explains the different purpose of repeated plausible completions.

Transform new rows by replaying the fitted sequence

01InitializeLearn the starting fills from training entries.
02Fit each updateStore the model fitted at each column update, in order.
03Transform new rowsUse the stored initial fills, then replay those models without refitting.
Fit the chain once per training fold. Scikit-learn retains the fitted sequence across rounds. It does not simply keep one final model per column and refit it on a test batch.

Here is the two-column example using IterativeImputer. Ordinary linear regression and the right-to-left order are explicit so the numerical example matches the figure. The library’s default estimator is Bayesian ridge regression, which need not produce these exact intermediate values.

Python: fit, preserve observed entries, and transform new rows 
import numpy as np
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
from sklearn.linear_model import LinearRegression

train = np.array([[1, 2], [2, np.nan], [np.nan, 6], [4, 8.]])
observed = ~np.isnan(train)
imputer = IterativeImputer(
    estimator=LinearRegression(),  # Matches the two-column illustration.
    initial_strategy="mean", imputation_order="arabic",  # B, then A.
    max_iter=20, tol=1e-6, random_state=0,
)
filled = imputer.fit_transform(train)
assert np.array_equal(filled[observed], train[observed])
assert np.allclose(filled, [[1, 2], [2, 4], [3, 6], [4, 8]], atol=1e-5)

new_rows = np.array([[3, np.nan], [np.nan, 10.]])
new_filled = imputer.transform(new_rows)
assert np.allclose(new_filled, [[3, 6], [5, 10]], atol=1e-5)
print(np.round(new_filled, 3))

[[ 3.  6.]
 [ 5. 10.]]

The second new row requires a prediction beyond the observed A range of 1–4. The code can produce 5; whether that extrapolation is defensible depends on the data. Check both output plausibility and downstream validation performance.

Choose models that respect what each feature means

A linear conditional model assumes an approximately linear relationship with its predictors. A count, binary variable, or bounded measurement may need a different model or constraints. Simply rounding a continuous prediction does not provide a suitable categorical imputation model, and clipping values to a range does not establish calibration.

The available features must also match the task. In a predictive pipeline, the outcome you are trying to predict is not available at serving time, so it cannot be an ordinary input to this feature imputer. For statistical inference on an existing dataset, including an observed analysis outcome in a carefully specified imputation model can be appropriate; that is a different workflow.

MAR alone does not guarantee good imputations: the relevant observed predictors must be available, the models must describe their relationships adequately, and the observed rows must cover the cases being filled. With MNAR, standard observed-data regressions can miss a systematic difference in the hidden values. Return to the mechanism discussion rather than inferring a mechanism from a model’s fit.

Practical controls worth checking
  • Cost: each round fits several models. Reduce unnecessary predictors or rounds only after assessing the effect on results. In scikit-learn, n_nearest_features selects predictor columns, not neighboring rows as KNN does.
  • New gaps: with skip_complete=True, a feature complete during fitting uses only its initial fill if it becomes missing later. The default False learns models for such features too.
  • Entirely missing columns: there are no observed targets to fit. Configure whether to retain them, and treat a constant fallback as a schema decision.
  • Reproducibility: fix the estimator, feature order, seed, and library version. Scikit-learn still marks this estimator experimental; the enabling import is intentional.

A repeated update is not automatically multiple imputation

The deterministic loop above returns one completed dataset. Fitting a Bayesian model inside that loop does not by itself make the output multiple imputation: posterior-mean predictions are still single fills.

MICE means multivariate imputation by chained equations. For multiple imputation, conditional models generate plausible draws and the process produces several completed datasets for analysis. Suitable stochastic models and pooling matter; repeating the same deterministic fit with different seed numbers is insufficient. Continue to multiple imputation and uncertainty for that workflow.

C++: reproduce the four-round illustration

This teaching version has exactly two numerical columns, retains a separate original mask, and uses ordinary least squares in both directions. If the predictor has no variation, it falls back to the observed target mean. It demonstrates the training loop on finite values; the Python version above supplies the reusable fit/transform workflow.

Compile the two-column chain
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <vector>

using Row = std::array<double, 2>;
using Mask = std::array<bool, 2>;

void updateColumn(std::vector<Row>& filled, const std::vector<Mask>& missing,
                  std::size_t target) {
    const std::size_t predictor = 1 - target;
    double meanX = 0, meanY = 0;
    std::size_t n = 0;
    for (std::size_t r = 0; r < filled.size(); ++r) {
        if (missing[r][target]) continue;  // Only observed training targets.
        meanX += filled[r][predictor];
        meanY += filled[r][target];
        ++n;
    }
    if (!n) throw std::domain_error("no observed target");
    meanX /= n;
    meanY /= n;
    double xx = 0, xy = 0;
    for (std::size_t r = 0; r < filled.size(); ++r) {
        if (missing[r][target]) continue;
        const double dx = filled[r][predictor] - meanX;
        xx += dx * dx;
        xy += dx * (filled[r][target] - meanY);
    }
    // No variation in the predictor: use an intercept-only mean prediction.
    const double slope = xx > 0 ? xy / xx : 0;
    const double intercept = meanY - slope * meanX;
    for (std::size_t r = 0; r < filled.size(); ++r)
        if (missing[r][target])
            filled[r][target] = intercept + slope * filled[r][predictor];
}

int main() {
    const double NA = std::numeric_limits<double>::quiet_NaN();
    const std::vector<Row> raw = {{1,2}, {2,NA}, {NA,6}, {4,8}};
    std::vector<Mask> missing;
    auto filled = raw;
    for (const auto& row : raw)
        missing.push_back({std::isnan(row[0]), std::isnan(row[1])});
    for (std::size_t c = 0; c < 2; ++c) {
        double sum = 0;
        std::size_t count = 0;
        for (std::size_t r = 0; r < raw.size(); ++r)
            if (!missing[r][c]) { sum += raw[r][c]; ++count; }
        if (!count) throw std::domain_error("entire column missing");
        for (std::size_t r = 0; r < raw.size(); ++r)
            if (missing[r][c]) filled[r][c] = sum / count;
    }
    for (int round = 0; round < 4; ++round) {
        updateColumn(filled, missing, 1);  // B from current A.
        updateColumn(filled, missing, 0);  // A from updated B.
    }
    std::cout << std::fixed << std::setprecision(6)
              << "row 2, B: " << filled[1][1] << "\n"
              << "row 3, A: " << filled[2][0] << "\n";
}

row 2, B: 4.000281
row 3, A: 2.999960

References and further reading