For statistical inference, extra variability introduced by imputation has to be accounted for. This is usually done by multiple imputation.

One of the standard approaches is to impute the dataset multiple
times, generating, e.g., 10 or 100 versions of the complete data. Then,
the intended analysis (t-test, linear model etc.) is performed with each
of the datasets. Their results are then pooled, usually by Rubin’s rule
(Rubin 1987):
Parameter *estimates* are averaged. Their *variances* are
avaraged as well, and corrected upwards by adding the variance of the
paramrter estimates across imputations.

The package {mice} (Buuren and Groothuis-Oudshoorn 2011) takes
care of this pooling step. The creation of multiple complete data sets
can be done by {mice} or also by {missRanger}. In the latter case, in
order to keep the variance of imputed values at a realistic level, we
suggest to use predictive mean matching with relatively large
`pmm.k`

on top of the random forest imputation.

```
library(missRanger)
library(mice)
set.seed(19)
irisWithNA <- generateNA(iris, p = c(0, 0.1, 0.1, 0.1, 0.1))
# Generate 20 complete data sets with relatively large pmm.k
filled <- replicate(
20,
missRanger(irisWithNA, verbose = 0, num.trees = 50, pmm.k = 10),
simplify = FALSE
)
# Run a linear model for each of the completed data sets
models <- lapply(filled, function(x) lm(Sepal.Length ~ ., x))
# Pool the results by mice
summary(pooled_fit <- pool(models))
# term estimate std.error statistic df p.value
# (Intercept) 2.4927868 0.34949255 7.1325891 80.75311 3.785713e-10
# Sepal.Width 0.4337064 0.10952987 3.9597088 78.74820 1.635231e-04
# Petal.Length 0.7362439 0.09201064 8.0017258 50.94966 1.453609e-10
# Petal.Width -0.1618100 0.18820061 -0.8597739 66.35077 3.930096e-01
# Speciesversicolor -0.7050250 0.30099170 -2.3423403 69.11545 2.204973e-02
# Speciesvirginica -0.9125055 0.41666494 -2.1900224 66.90964 3.201281e-02
# Compare with model on original data
summary(lm(Sepal.Length ~ ., data = iris))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 2.17127 0.27979 7.760 1.43e-12 ***
# Sepal.Width 0.49589 0.08607 5.761 4.87e-08 ***
# Petal.Length 0.82924 0.06853 12.101 < 2e-16 ***
# Petal.Width -0.31516 0.15120 -2.084 0.03889 *
# Speciesversicolor -0.72356 0.24017 -3.013 0.00306 **
# Speciesvirginica -1.02350 0.33373 -3.067 0.00258 **
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
# Residual standard error: 0.3068 on 144 degrees of freedom
# Multiple R-squared: 0.8673, Adjusted R-squared: 0.8627
# F-statistic: 188.3 on 5 and 144 DF, p-value: < 2.2e-16
```

As expected, standard errors and p values of the multiple imputation are larger than of the original data set.

Buuren, Stef van, and Karin Groothuis-Oudshoorn. 2011. “Mice:
Multivariate Imputation by Chained Equations in r.” *Journal
of Statistical Software, Articles* 45 (3): 1–67. https://doi.org/10.18637/jss.v045.i03.

Rubin, D. B. 1987. *Multiple Imputation for Nonresponse in
Surveys*. Wiley Series in Probability and Statistics. Wiley.