Multiple imputation para substituição de missing não é uma muito técnica muito nova. O que existe de mais moderno seria usar os modelos mistos generalizados para estimar seus coeficiente desconsiderando os dados faltantes.
Caso isso não seja possível, é claro cada casa é um caso. O pacote Amélia, VIM e mi tem funções para multiple imputation.
Segue um exemplo.
# "Iterative Robust Model-based Input" através da função irmi do pacote "VIM"
library(VIM)
data(sleep, package = "VIM")
sleep$NonD[5] <- NA # 2.1
sleep$NonD[27] <- NA # 10.4
irmi_sleep <- irmi(sleep)
irmi_sleep
# "Bootstrapped expectation-maximization" através da função amelia do pacote Amelia
library(Amelia)
amelia_sleep <- amelia(sleep, bounds = matrix(c(3, 0, 30), ncol = 3))
amelia_data_sleep <- amelia_sleep$imputations$imp5
amelia_data_sleep
library(mi)
mi_sleep <- mi(sleep, info = info, n.iter=10, add.noise=FALSE)
mi_data_sleep <- mi.data.frame(mi_sleep)