Neural-network imputation with deepImp

Why neural-network imputation?

Most real data sets contain missing values, and discarding incomplete rows wastes information and can bias an analysis. Imputation fills the gaps with plausible values so the data can be analysed as a whole. Classical imputation methods (mean imputation, $k$-nearest neighbours, regression/PMM) assume a fixed, usually linear, relationship between a variable and its predictors.

deepImp instead learns that relationship with a neural network, which can capture non-linear and interaction structure without the user specifying it. It imputes one variable at a time, using all the others as predictors, and cycles through the variables repeatedly until the imputations stabilise (an iterative, chained scheme). The package handles two settings with the same engine:

mixed-type data – numeric, count, semi-continuous, binary and nominal variables – via impNNet();
compositional data with rounded zeros – values that fall below a detection limit – via impNNetCoDa().

This vignette walks through a complete analysis: variable types, the network and its training settings, point imputation and its diagnostics, multiple imputation with pooled inference, imputing new data with stored models, and compositional data.

Networks run on a native torch backend by default (no Python required), or optionally on keras3.

library(deepImp)

Step 1: the data and the variable types

We use the sleep data from the VIM package – mammal sleep measurements that already contain genuine missing values.

data(sleep, package = "VIM")
str(sleep)
#> 'data.frame':    62 obs. of  10 variables:
#>  $ BodyWgt : num  6654 1 3.38 0.92 2547 ...
#>  $ BrainWgt: num  5712 6.6 44.5 5.7 4603 ...
#>  $ NonD    : num  NA 6.3 NA NA 2.1 9.1 15.8 5.2 10.9 8.3 ...
#>  $ Dream   : num  NA 2 NA NA 1.8 0.7 3.9 1 3.6 1.4 ...
#>  $ Sleep   : num  3.3 8.3 12.5 16.5 3.9 9.8 19.7 6.2 14.5 9.7 ...
#>  $ Span    : num  38.6 4.5 14 NA 69 27 19 30.4 28 50 ...
#>  $ Gest    : num  645 42 60 25 624 180 35 392 63 230 ...
#>  $ Pred    : int  3 3 1 5 3 4 1 4 1 1 ...
#>  $ Exp     : int  5 1 1 2 5 4 1 5 2 1 ...
#>  $ Danger  : int  3 3 1 3 4 4 1 4 1 1 ...
colSums(is.na(sleep))
#>  BodyWgt BrainWgt     NonD    Dream    Sleep     Span     Gest     Pred 
#>        0        0       14       12        4        4        4        0 
#>      Exp   Danger 
#>        0        0

impNNet() first infers the measurement scale of each column with guessType(). The type decides the head of the network: numeric, count and semi-continuous (“mixed”) columns get a regression head, binary and nominal columns a classification head. Override the guess with the vartypes argument when it is wrong for your data.

guessType(sleep)$type
#>  [1] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
#>  [8] "count"   "count"   "count"

Step 2: the network and its training

The network is described by a deepimp_arch() object: the depth and width (hidden), the regularisation (dropout, batchnorm, weight_decay), the activation, and the optimizer with its learning_rate.

arch <- deepimp_arch(hidden = c(128, 64, 32), dropout = 0.1, activation = "relu")
arch
#> <deepimp_arch>
#>   hidden layers : 128 -> 64 -> 32 
#>   dropout       : 0.1 
#>   activation    : relu 
#>   batchnorm     : TRUE 
#>   optimizer     : adam (lr = 0.001 )
#>   GLM offset    : yes (ridge = 0.001 )

Every head is a generalised linear model plus a network correction (glm_offset = TRUE, the default): a linear model under the regression heads and a multinomial logit under the classification heads. The correction starts at zero, and training keeps it only where it improves the fit on held-out rows, so a head is never worse there than its GLM. The network then only has to learn what the GLM misses, such as non-linear effects and interactions. For small data sets deepimp_arch_small() is a lighter preset:

deepimp_arch_small()
#> <deepimp_arch>
#>   hidden layers : 64 -> 32 
#>   dropout       : 0.1 
#>   activation    : relu 
#>   batchnorm     : TRUE 
#>   optimizer     : adam (lr = 0.001 )
#>   GLM offset    : yes (ridge = 0.001 )

The training settings live in a deepimp_control() object: the number of epochs, early stopping (patience), the minibatch size, the learning-rate schedule, the device and thread count, warm starts across sweeps, and parallel replicates.

ctl <- deepimp_control(epochs = 50)
ctl
#> <deepimp_control>
#>   epochs / patience : 50 / 10 
#>   batch size        : auto 
#>   validation split  : 0.2 
#>   restore best      : TRUE 
#>   lr schedule       : none 
#>   device            : cpu 
#>   threads           : unchanged 
#>   warm start        : point mode only (default) 
#>   parallel          : FALSE 
#>   min. observed     : 20

By default, early stopping waits 10 epochs for an improvement, the minibatch size grows with the data (about 32 minibatches per epoch), and in point mode each column model continues from its weights of the previous sweep instead of starting over (warm_start), which saves most of the training time after the first sweep. For quick runs, impNNet() also accepts epochs, patience, batch_size and validation_split directly.

Step 3: point imputation

By default (draw = "none") every missing cell gets the network’s best guess: the conditional mean for a regression head, the most probable class for a classification head. The chain cycles through the incomplete variables until the imputations stop changing (or iterations sweeps are done). The seed makes the run reproducible, exactly so with dropout = 0.

imp <- impNNet(sleep, arch = deepimp_arch_small(dropout = 0), control = ctl,
               iterations = 5, seed = 1)
imp
#> <deepimp>
#>   observations : 62
#>   variables    : 10
#>   imputations  : 1
#>   backend      : torch
#>   architecture :
#> <deepimp_arch>
#>   hidden layers : 64 -> 32 
#>   dropout       : 0 
#>   activation    : relu 
#>   batchnorm     : TRUE 
#>   optimizer     : adam (lr = 0.001 )
#>   GLM offset    : yes (ridge = 0.001 )

getImputed() (or as.data.frame()) returns the completed data. The non-dreaming sleep variable NonD had missing entries; here are the values that were filled in:

completed <- getImputed(imp)
missing_rows <- which(is.na(sleep$NonD))
round(completed$NonD[missing_rows], 2)
#>  [1]  2.52 10.10 14.09  2.41  3.49 10.49 10.59  8.50  8.89  2.09 10.69  1.90
#> [13]  9.20 13.50

Point imputations suit prediction-type uses. They understate the uncertainty of the imputed values, so standard errors computed from one completed data set are too small; Step 5 shows how to get valid inference.

Step 4: diagnostics

summary() reports the mode, the sweeps run, whether the chain converged, the epochs used per fit, and any variable that a safeguard caught (for example, too few observed values for a network fit):

summary(imp)
#> <deepimp> summary
#>   observations : 62
#>   variables    : 10
#>   imputations  : 1 (point mode)
#>   backend      : torch
#>   missing per variable:
#>  NonD Dream Sleep  Span  Gest 
#>    14    12     4     4     4 
#>   sweeps run   : 5 of 5
#>   converged    : 0 of 1 replicate(s); last change 0.0428
#>   epochs/fit   : median 4 (range 2-50)

plot() shows the convergence of the chain, the mean standardised change of the imputed cells from one sweep to the next:

plot(imp)

plot of chunk impnnet-plot

Here the chain stopped at the sweep cap: the change levels off at a few hundredths instead of falling below eps = 0.01. The trace of the imputed cells shows where the movement comes from. The means of the imputed NonD and Dream values are stable, while those of Span and Gest (life span and gestation time, skewed, four missing values each) move from sweep to sweep, because their models are refitted on 58 rows in every sweep:

plot(imp, type = "trace", vars = c("NonD", "Span", "Gest"))

plot of chunk impnnet-trace

The training history of every fit is kept in imp$info$history:

head(imp$info$history)
#>   imputation sweep variable epochs best_epoch   best_valid status
#> 1          1     1     NonD     10          0 4.310555e-09   <NA>
#> 2          1     1    Dream     19          9 1.817723e-06   <NA>
#> 3          1     1    Sleep     10          0 5.589238e-09   <NA>
#> 4          1     1     Span     50         49 1.147252e+00   <NA>
#> 5          1     1     Gest     11          1 3.586687e-01   <NA>
#> 6          1     2     NonD      2          0 3.724258e-09   <NA>

Step 5: multiple imputation and pooled inference

For inference, use draw mode. With draw = "pmm" (predictive mean matching) the regression heads copy the observed value of a donor whose prediction is close to the prediction for the missing cell, chosen at random among the pmm_k nearest. With draw = "normal" they add normal noise with the residual standard deviation of the held-out rows. The classification heads sample from their predicted class probabilities. By default every fit also uses a bootstrap resample of the observed rows, with the left-out rows for early stopping, so that the m completed data sets reflect the uncertainty of the model fit as well. Draw mode runs a fixed number of sweeps.

mi <- impNNet(sleep, arch = deepimp_arch_small(dropout = 0), m = 5,
              draw = "pmm", iterations = 5,
              control = deepimp_control(epochs = 50), seed = 1)
summary(mi)
#> <deepimp> summary
#>   observations : 62
#>   variables    : 10
#>   imputations  : 5 (pmm mode, multiple imputation)
#>   backend      : torch
#>   missing per variable:
#>  NonD Dream Sleep  Span  Gest 
#>    14    12     4     4     4 
#>   sweeps run   : 5 of 5
#>   convergence  : fixed sweeps (draw mode); check plot(x, type = "trace")
#>   epochs/fit   : median 14 (range 10-50)

In draw mode, plot(type = "trace") shows the mean and standard deviation of the imputed cells per sweep, one line per imputation. Lines that mix without a trend indicate that the chains have settled.

plot(mi, type = "trace", vars = c("NonD", "Dream", "Span"))

plot of chunk mi-trace

as_mids() hands the completed data sets to the mice package, which fits the analysis model to each of them and pools the results with Rubin’s rules:

fit <- with(as_mids(mi), lm(Sleep ~ log(BodyWgt) + Danger))
summary(mice::pool(fit))
#>           term   estimate std.error statistic       df      p.value
#> 1  (Intercept) 15.3451389 0.8757023 17.523236 55.46903 2.823592e-24
#> 2 log(BodyWgt) -0.5853983 0.1415588 -4.135371 56.89067 1.180026e-04
#> 3       Danger -1.5899571 0.3092113 -5.141977 55.71009 3.631827e-06

getImputed(mi, m = "long") returns the same data in the long format of mice for other tools.

Step 6: imputing new data

With keep_models = TRUE the fitted models are stored in the result (as plain R arrays, so the object can be saved with saveRDS()). predict() then imputes new observations without retraining, in the mode of the fit:

train <- sleep[1:50, ]
newdata <- sleep[51:62, ]
fit_train <- impNNet(train, arch = deepimp_arch_small(dropout = 0), control = ctl,
                     iterations = 5, seed = 1, keep_models = TRUE)
new_imp <- predict(fit_train, newdata)
getImputed(new_imp)[is.na(newdata$NonD), c("NonD", "Dream", "Sleep")]
#>         NonD     Dream    Sleep
#> 53  1.902203 0.6962252  2.60000
#> 55  8.701982 2.2962568 11.00000
#> 62 13.508102 3.8966081 17.41293

Step 7: compositional data and rounded zeros

Compositional data carry only relative information: what matters is the ratios between parts (for example, the proportions of chemical elements in a soil sample), not their absolute magnitudes. A rounded zero is a part recorded as zero only because its true value fell below a measurement detection limit – it is not a true absence. Ordinary imputation is inappropriate here for two reasons: it ignores the relative (log-ratio) geometry of the data, and it may return values above the detection limit, which contradicts the very reason the value was recorded as zero.

impNNetCoDa() (which needs the robCompositions package) addresses both: it works in pivot log-ratio coordinates and censors each imputed value at the detection limit. We illustrate with a small three-part composition in which some values of part a fall below a limit of 1 and are therefore recorded as zero:

set.seed(1)
n <- 60
x <- data.frame(
  a = runif(n, 5, 10),
  b = runif(n, 5, 10),
  c = runif(n, 5, 10)
)
dl <- c(1, 1, 1)          # detection limit per part
x$a[1:6] <- 0             # rounded zeros: recorded 0, truly below the limit
head(x)
#>   a        b        c
#> 1 0 9.564380 9.959193
#> 2 0 6.468017 7.477968
#> 3 0 7.295329 7.421748
#> 4 0 6.661973 5.867212
#> 5 0 8.254352 8.774105
#> 6 0 6.290084 7.269477

We pass the detection limits in dl and tell the function which value marks a rounded zero with label = 0:

impc <- impNNetCoDa(x, dl = dl, label = 0, correction = "truncate",
                    arch = deepimp_arch_small(), epochs = 50,
                    iterations = 2, seed = 1)
filled <- getImputed(impc)
head(filled)
#>   a        b        c
#> 1 1 9.564380 9.959193
#> 2 1 6.468017 7.477968
#> 3 1 7.295329 7.421748
#> 4 1 6.661973 5.867212
#> 5 1 8.254352 8.774105
#> 6 1 6.290084 7.269477

The rounded zeros in part a have been replaced by strictly positive values that respect the detection limit:

filled$a[1:6]
#> [1] 1 1 1 1 1 1
all(filled$a[1:6] > 0 & filled$a[1:6] <= dl[1])
#> [1] FALSE

The correction argument controls how the prediction is censored:

"truncate" (default) sets any imputed value above the limit to the limit, as in Templ (2021);
"expectation" uses the mean of a normal distribution truncated at the limit, pulling values strictly below it;
"none" applies no censoring (useful for diagnostics or ablation).

Setting coda = FALSE bypasses the log-ratio transform and imputes on the raw scale, which is occasionally useful for comparison.

Choosing a backend

Both functions accept backend = "torch" (the default) or backend = "keras". torch is the primary backend: it runs on a native libtorch library without Python, and every feature of the package is available with it.

install.packages("torch")
torch::install_torch()          # downloads libtorch once

keras3 is an optional, secondary backend for users who already work with Keras/TensorFlow. It supports point and draw mode, but not keep_models, parallel replicates or the device setting. It needs a Python environment with Keras and TensorFlow, which keras3 sets up once:

install.packages("keras3")
keras3::install_keras()         # creates a Python virtual environment
imp_keras <- impNNet(sleep, backend = "keras",
                     arch = deepimp_arch_small(), epochs = 50, seed = 1)

The two backends implement the same architecture, split the data the same way and apply the same early-stopping and restore-best rules, but they do not produce numerically identical results: their weight initialisation, random-number streams and optimiser internals differ.

Reproducibility

A seed makes a run repeatable within a backend. Reproducibility is exact when dropout = 0; with dropout > 0 the backend’s training-time mask sampling is not fully pinned by the seed in the current torch build, so repeated runs are close but may not be bit-for-bit identical. Parallel replicates (deepimp_control(parallel = TRUE)) give the same results as a sequential run with the same seed and the same number of torch threads (num_threads).