Why Modular MICE Scripts Matter
A MICE analysis is a sequence of related operations: load the observed data, prepare it for analysis, impute missing values, analyze each completed dataset, pool the results, and produce outputs. When all of these operations are mixed together in one long script, it becomes difficult to tell which step changed a result, reuse one step safely, or test a step without rerunning everything.
A modular script assigns each stage a defined responsibility and gives it explicit inputs and outputs. The goal is not to create the largest possible number of functions. It is to make the path from observed data to reported result visible, with each function doing one coherent job. This complements the earlier material on reproducible MICE projects and random seeds: the script’s structure should make the intended analysis inspectable and repeatable, not alter its statistical specification.
The examples use an R-style description of functions and objects, but the design applies in any language. The function signatures below are pseudocode rather than package-specific code. In particular, a function named impute_data represents the project’s chosen MICE implementation; it does not prescribe a particular conditional model or imputation method.
Set the Contracts Before Connecting the Stages
A function contract states what a function receives, what it returns, and what it is responsible for. A useful division for an applied analysis is:
- load_data(path) receives a data location and returns the raw dataset as read, without silently changing its scientific meaning.
- preprocess_data(raw_data) receives the raw dataset and returns the analysis-ready observed data, applying documented, deterministic preparation steps.
- impute_data(data, m, seed, specification) receives the prepared data and the stated imputation inputs, and returns the MICE run object or the set of completed datasets.
- analyze_imputations(imputed_data, analysis_specification) applies the same substantive analysis to each completed dataset and returns estimates and their within-imputation variances.
- write_outputs(results, output_path) receives results and an output destination, and writes agreed outputs. It should not recalculate the analysis using a different specification.
These boundaries are meaningful only if the contracts match. For example, if preprocessing converts a categorical predictor to a factor, the output should include that prepared variable rather than modifying a global dataset behind the caller’s back. If the analysis function returns only point estimates, it has not returned enough information for ordinary Rubin-rule pooling: it must also return a variance estimate for each completed-data estimate.
The statistical specification must also remain stable across the software boundaries. For instance, when the target is a treatment contrast in mean systolic blood pressure, every completed-data analysis must estimate that same contrast in the same units. The MICE specification must be appropriate for the analysis and the missing-data assumptions being used; dividing code into functions does not establish those assumptions or make an unsuitable model valid.
Worked Examples of Individual Modules
Worked Example: Loading and Preparing a Hypothetical Dataset
Suppose a hypothetical study file contains 120 participant rows, an identifier, treatment group, age, and a blood-pressure outcome. Eight outcome values are missing. Two rows have the same participant identifier because of a duplicated export record. The analysis plan says to retain one verified record per participant and not to discard participants just because their outcome is missing.
Loading contract. The loading function reads the specified file and returns the 120-row table exactly as imported. It does not decide which duplicate to remove or fill in missing outcomes. That makes it possible to inspect the imported data before preparation.
Preprocessing contract. The preprocessing function checks the identifier rule, removes the one confirmed duplicate record, and standardizes the treatment labels. Its output has 119 rows. It leaves the eight missing outcomes missing so the imputation stage can address them. The output row count follows directly: \(120-1=119\).
A useful set of returned checks would state that there are 119 rows, 119 unique identifiers, and eight missing outcome values. If preprocessing instead returned 111 rows because it had dropped every row with a missing outcome, it would conflict with the stated analysis plan. The contract and these counts expose that error before MICE is run.
Preprocessing functions should have visible, narrow responsibilities. A function that both removes records and creates imputed values is harder to review: a row-count change could arise from either operation, and a later analyst may not know whether missingness was intentionally retained for imputation.
Worked Example: Returning Enough Information for Pooling
Suppose an analysis estimates the adjusted treatment contrast in systolic blood pressure, measured in mmHg, in \(m=5\) completed datasets. The estimates are \(2.0,\ 2.4,\ 1.8,\ 2.2,\) and \(2.1\) mmHg. Assume, for this illustration, that each completed-data analysis reports a variance of \(1.44\ \text{mmHg}^2\), corresponding to a standard error of \(1.2\) mmHg.
The analysis function should return five estimate-variance pairs, not just five estimates. The pooled estimate is their mean:
The between-imputation variance is the sample variance of the five estimates. Their deviations from \(2.10\) are \(-0.10,\ 0.30,\ -0.30,\ 0.10,\) and \(0\) mmHg. Therefore,
The mean within-imputation variance is \(\bar{U}=1.44\ \text{mmHg}^2\). Under the usual multiple-imputation pooling conditions, including proper imputations and a common estimand with an appropriate complete-data variance estimate, Rubin’s total variance is
Thus the pooled standard error is \(\sqrt{1.50}=1.225\) mmHg. The result passed between modules should preserve the pooled estimate, its uncertainty, the units, and the number of imputations. Returning only \(2.10\) would omit information needed for valid inference.
Justifying the Pooling Boundary
Separating analysis from pooling is useful because the two operations have different inputs. Each completed-data analysis produces an estimate \(Q_j\) and its estimated variance \(U_j\), for \(j=1,\ldots,m\). Pooling combines those paired results; it is not equivalent to averaging the completed datasets first and analyzing that single average.
The formula separates two sources of uncertainty. The mean \(\bar{U}\) summarizes uncertainty that would remain within a completed dataset. The variance \(B\) measures how estimates differ across imputations, reflecting uncertainty associated with missing values. The multiplier \(1+1/m\) includes the finite-\(m\) adjustment. If every module uses the same definition of \(Q_j\) and \(U_j\), the pooling function can operate on a compact, auditable set of results rather than needing access to hidden analysis objects.
This boundary also makes a key check possible: the number of estimates must equal the number of variance estimates and the declared number of imputations. If the analysis stage returns five estimates but four variances, pooling should stop with an informative error rather than silently pairing incompatible results.
Connect the Modules Through One Entry Point
The entry point is the function a user calls to run the planned analysis. It owns the sequence of operations and passes each output to the next stage. Its arguments should make important run choices visible, such as the input path, the number of imputations, the seed, and the output location. The upcoming material on configuration files addresses a different way to store run choices; a modular entry point can be designed without introducing that mechanism.
Call the loading function with the input path and receive the raw data.
Pass the raw data to preprocessing and receive the analysis-ready data, retaining the missing values intended for imputation.
Pass prepared data and explicit imputation arguments to the MICE function; retain its returned run object or completed datasets.
Analyze each completed dataset under the same analysis specification, then combine the estimate-variance pairs.
Send the results to the output function and return a result object to the caller.
Worked Example: A Single Entry Point
Consider a hypothetical run with the 119 analysis-ready participants above, \(m=5\) imputations, seed 731, and an output destination named “run_A”. Suppose the imputation and analysis stages succeed and return the five estimate-variance pairs from the previous example.
The entry point receives those run choices explicitly: conceptually, run_analysis(input_path, output_path, m, seed). It calls the stages in order and returns a result containing the pooled estimate \(2.10\) mmHg, total variance \(1.50\ \text{mmHg}^2\), standard error \(1.225\) mmHg, and \(m=5\). The output function writes the requested summary for “run_A”. The calculation of the standard error is \(\sqrt{1.50}=1.225\) mmHg; the entry point does not replace this result with a separately calculated value.
If the prepared dataset has 119 rows but the imputation stage reports 118 records, the entry point should not continue as if the objects were interchangeable. It should stop and identify the failed stage or violated contract. A clear failure at the boundary is preferable to a plausible-looking output based on mismatched data.
Checks, Alternatives, and Failure Modes
Modularity does not mean that every stage must be a separate file or that every short line deserves its own function. It means the workflow’s important operations and data flow are explicit. A single entry point should orchestrate the stages, while each stage remains testable using a small, controlled input.
- Look for hidden state. If a function reads a global dataset or silently changes the working directory, its declared arguments do not explain what determines its result.
- Check contracts at boundaries. Confirm row counts, variable names, the number of imputation-specific results, and the units and meaning of estimates where appropriate.
- Watch for duplicated analysis logic. If a model is separately written in several locations, one copy can drift from the others. Keep the repeated substantive analysis in one function called for each completed dataset.
- Do not collapse stages to make the script look shorter. Dropping incomplete cases during preprocessing, for example, changes which data reach imputation and can change the analysis population.
- Avoid excessive fragmentation. A function with no clear responsibility, no useful contract, and no independent test value can make a workflow harder to follow rather than easier.
A monolithic script can be convenient for a one-off exploration, but its implicit state and interleaved operations make review and reuse harder. A notebook can support exploration and explanation, yet execution order may be unclear if cells are run out of sequence. A modular script with one entry point is often preferable for a planned repeatable analysis because it makes the order and arguments explicit. None of these formats guarantees valid inference: the MICE assumptions, analysis specification, and interpretation still require statistical judgment.
Check Your Understanding
Answer each question, then compare your reasoning with the brief solutions.
- Why should a preprocessing function generally leave outcomes missing when those values are intended for MICE?
- An analysis function returns six estimates but only five variance estimates. What should the pooling stage do, and why?
- In a run with \(m=4\), the completed-data estimates are \(1,\ 2,\ 2,\) and \(3\) units. Calculate \(\bar{Q}\) and \(B\).
- Give one example of a hidden input that weakens a function’s contract.
- Does splitting a script into functions establish that the imputation model is appropriate? Explain.
Solutions. (1) Dropping those outcomes would remove cases intended for imputation and could change the analysis population. (2) Stop with an informative error: pooling requires a matched estimate and variance for each of the same \(m\) completed datasets. (3) \(\bar{Q}=(1+2+2+3)/4=2\) units; \(B=[(1-2)^2+(2-2)^2+(2-2)^2+(3-2)^2]/(4-1)=2/3=0.667\) units\(^2\). (4) A global data object read inside a function without appearing among its arguments is a hidden input. (5) No. Modularity improves organization and checking, but the imputation model and its assumptions still need statistical justification.