README.Rmd

---
output:
  github_document:
    pandoc_args: --webtex=https://render.githubusercontent.com/render/math?math=
bibliography: ref.bib
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%")
options(digits = 3)
```

# mdgc

[![R-CMD-check](https://github.com/boennecd/mdgc/workflows/R-CMD-check/badge.svg)](https://github.com/boennecd/mdgc/actions) 
[![](https://www.r-pkg.org/badges/version/mdgc)](http://cran.rstudio.com/web/packages/mdgc/index.html) 
[![CRAN RStudio mirror
downloads](http://cranlogs.r-pkg.org/badges/mdgc)](http://cran.rstudio.com/web/packages/mdgc/index.html)

This package contains a marginal likelihood approach to estimating the 
model discussed by @hoff07, @zhao20, and @zhao20Mat. That is, a missing data approach
where one uses Gaussian copulas in the latter case. 
We have modified the Fortran
code by @Genz02 to supply an approximation of the gradient for the log 
marginal likelihood and to use an approximation of the marginal likelihood 
similar to the CDF approximation in @Genz02. We have also used the same 
Fortran code to perform the imputation conditional on a covariance matrix
and the observed data. The method is described by @Christoffersen21 which can 
be found at [arxiv.org](https://arxiv.org/abs/2102.02642).

Importantly, we also extend the model used by @zhao20 to support multinomial 
variables. Thus, our model supports both continuous, binary, ordinal, and 
multinomial variables which makes it applicable to a large number of data sets.

The package can be useful for a lot of other models. For instance, 
the methods are directly applicable to other Gaussian copula models and some 
mixed effect models. All methods are implemented in C++, support computation 
in parallel, and should easily be able to be ported to other languages. 

## Installation
The package can be installed from Github by calling:

```{r, eval = FALSE}
remotes::install_github("boennecd/mdgc")
```

or from CRAN by calling:

```{r, eval = FALSE}
install.packages("mdgc")
```

The code benefits from being build with automatic vectorization so having e.g.  
`-O3 -mtune=native` in the `CXX11FLAGS` flags in your Makevars file may be 
useful.

## The Model
We observe four types of variables for each observation: continuous, 
binary, ordinal, and multinomial variables. Let $\vec X_i$ be a K
dimensional vector for the i'th observation. The variables $X_{ij}$ are 
continuous 
if $j\in\mathcal C$, binary if $j\in\mathcal B$ with probability $p_j$ of 
being true, ordinal if $j\in\mathcal O$ with $m_j$ levels and borders 
$\alpha_{j0} = -\infty < \alpha_1<\cdots < \alpha_{m_j} = \infty$, and 
multinomial if $j\in\mathcal M$ with $m_j$ levels. $\mathcal C$, 
$\mathcal B$, $\mathcal O$, and $\mathcal M$ are mutually exclusive.

We assume that there is
a latent variable $\vec Z_i$ which is multivariate normally distributed such 
that:

<!-- $$ -->
<!-- \begin{align*} -->
<!-- \vec Z_i & \sim N\left(\vec\mu, -->
<!--   \Sigma\right) \nonumber\\ -->
<!-- X_{ij} &= f_j(Z_{ih(j)}) & j &\in \mathcal C \\ -->
<!-- X_{ij} &= \begin{cases} -->
<!--   1 & Z_{ij} > \underbrace{-\Phi^{-1}(p_{j})}_{\mu_{h(j)}} \\ -->
<!--   0 & \text{otherwise}   -->
<!-- \end{cases} & j &\in \mathcal B \\ -->
<!-- X_{ij} &= k\Leftrightarrow \alpha_{jk} < Z_{ih(j)} \leq \alpha_{j,k + 1} -->
<!--   & j &\in \mathcal O\wedge k = 0,\dots m_j -1 \\ -->
<!-- X_{ij} &= k \Leftrightarrow Z_{i,h(j) + k} \geq -->
<!--   \max(Z_{ih(j)},\cdots,Z_{i,h(j) + m_j - 1}) -->
<!--   & j&\in \mathcal M \wedge k = 0,\dots m_j -1 -->
<!-- \end{align*} -->
<!-- $$ -->

$$\begin{align*}  \vec Z_i & \sim N\left(\vec\mu,    \Sigma\right) \nonumber\\  X_{ij} &= f_j(Z_{ih(j)}) & j &\in \mathcal C \\  X_{ij} &= 1_{\{Z_{ih(j)} > \underbrace{-\Phi^{-1}(p_{j})}_{\mu_{h(j)}}\}} & j &\in \mathcal B \\  X_{ij} &= k\Leftrightarrow \alpha_{jk} < Z_{ih(j)} \leq \alpha_{j,k + 1}    & j &\in \mathcal O\wedge k = 0,\dots m_j -1 \\  X_{ij} &= k \Leftrightarrow Z_{i,h(j) + k} \geq    \max(Z_{ih(j)},\cdots,Z_{i,h(j) + m_j - 1})    & j&\in \mathcal M \wedge k = 0,\dots m_j -1  \end{align*}$$

where $1_{\{\cdot\}}$ is one if the condition in the subscript is true and zero 
otherwise, $h(j)$ is a map to the index of the first latent variable associated with 
the j'th variable in $\vec X_i$ and $f_j$ is a bijective function. We only 
estimate some of the means, the $\vec\mu$, and some of the covariance 
parameters. Furthermore, we set $Z_{ih(j)} = 0$ if $j\in\mathcal M$ and assume
that the variable is uncorrelated with all the other $\vec Z_i$'s.

In principle, we could use other distributions than a multivariate normal 
distribution for $\vec Z_i$. However, the multivariate normal distribution
has the advantage that it is very easy to marginalize which is convenient when 
we have to estimate the model with missing entries and it is also has some 
computational advantages for approximating the log marginal likelihood as 
similar intractable problem have been thoroughly studied.

## Examples
Below, we provide an example similar to @zhao20 [Section 7.1]. The authors 
use a data set with a random correlation matrix, 5 continuous variables, 
5 binary variables, and 5 ordinal variables with 5 levels. There is a total 
of 2000 observations and 30% of the variables are missing completely at 
random.

To summarize @zhao20 results, they show that their approximate EM algorithm 
converges in what seems to be 20-25 seconds (this is with a pure R 
implementation to be fair) while it takes more than 150
seconds for the MCMC algorithm used by @hoff07. These figures
should be kept in mind when looking at the results below. Importantly, 
@zhao20 use an approximation in the E-step of an EM algorithm which is 
fast but might be crude in some settings. Using a potentially arbitrarily 
precise approximation of the log  marginal likelihood is useful if this can 
be done quickly enough.

We will provide a [quick example](#quick-example) and 
[an even shorter example](#an-even-shorter-example)
where we show how to use the methods in the package to estimate the 
correlation matrix and to perform the imputation. We then show a 
[simulation study](#simulation-study) where we compare with the method 
suggested by @zhao20.

The last section called [adding multinomial variables](#adding-multinomial-variables)
covers data sets which also have multinomial variables.

## Quick Example

We first simulate a data set and provide an example which shows how to use 
the package. The [an even shorter example](#an-even-shorter-example) section 
shows a shorter example then what is shown here. You may want to see this 
first if you just want to perform some quick imputation.

```{r load_pkgs}
# load the packages we need
library(bench)
library(mdgc)
library(missForest, quietly = TRUE)
# remotes::install_github("udellgroup/mixedgcImp", ref = "5ad6d523d")
library(mixedgcImp)
library(doParallel)
```

```{r sim_dat, cache = 1, fig.height = 3, fig.width = 7}
# simulates a data set and mask some of the data.
#
# Args: 
#   n: number of observations. 
#   p: number of variables. 
#   n_lvls: number of levels for the ordinal variables. 
# 
# Returns: 
#   Simulated masked data, the true data, and true covariance matrix. 
sim_dat <- function(n, p = 3L, n_lvls = 5L){
  # get the covariance matrix
  Sig <- cov2cor(drop(rWishart(1L, p, diag(p))))
    
  # draw the observations
  truth <- matrix(rnorm(n * p), n) %*% chol(Sig)
  
  # determine the type
  n_rep <- floor((p + 3 - 1) / 3)
  type <- rep(1:3, each = n_rep)[1:p]
  is_con <- type == 1L
  is_bin <- type == 2L
  is_ord <- type == 3L
  col_nam <- c(outer(1:n_rep, c("C", "B", "O"), 
                     function(x, y) paste0(y, x)))[1:p]
  
  # sample which are masked data 
  is_mask <- matrix(runif(n * p) < .3, n)
  
  # make sure we have no rows with all missing data
  while(any(all_nans <- rowSums(is_mask) == NCOL(is_mask)))
    is_mask[all_nans, ] <- runif(sum(all_nans) * p) < .3
  
  # create observed data
  truth_obs <- data.frame(truth)
  colnames(truth_obs) <- col_nam
  truth_obs[, is_con] <- qexp(pnorm(as.matrix(truth_obs[, is_con])))
  
  bs_border <- 0
  truth_obs[, is_bin] <- 
    truth_obs[, is_bin] > rep(bs_border, each = NROW(truth_obs))
  
  bs_ord <- qnorm(seq(0, 1, length.out = n_lvls + 1L))
  truth_obs[, is_ord] <- as.integer(cut(truth[, is_ord], breaks = bs_ord))
  for(i in which(is_ord)){
    truth_obs[, i] <- ordered(truth_obs[, i])
    levels(truth_obs[, i]) <- 
      LETTERS[seq_len(length(unique(truth_obs[, i])))]
  }

  # mask the data
  seen_obs <- truth_obs
  seen_obs[is_mask] <- NA
  
  list(truth = truth, truth_obs = truth_obs, seen_obs = seen_obs, 
       Sigma = Sig)
}

# simulate and show the data
set.seed(1)
p <- 15L
dat <- sim_dat(2000L, p = p)

# how an observed data set could look
head(dat$seen_obs)

# assign objects needed for model estimation
mdgc_obj <- get_mdgc(dat$seen_obs)
log_ml_ptr <- get_mdgc_log_ml(mdgc_obj)
start_val <- mdgc_start_value(mdgc_obj)

# this is very fast so we can neglect this when we consider the computation 
# time
mark(`Setup time` = {
  mdgc_obj <- get_mdgc(dat$seen_obs)
  log_ml_ptr <- get_mdgc_log_ml(mdgc_obj)
  start_val <- mdgc_start_value(mdgc_obj)
}, min_iterations = 10)

# fit the model using three different methods
set.seed(60941821)
system.time(
  fit_Lagran_start <- mdgc_fit(
    ptr = log_ml_ptr, vcov = start_val, mea = mdgc_obj$means, 
    n_threads = 4L, maxit = 100L, method = "aug_Lagran", rel_eps = 1e-3, 
    maxpts = 200L))
system.time(
  fit_Lagran <- mdgc_fit(
    ptr = log_ml_ptr, vcov = fit_Lagran_start$result$vcov, 
    mea = fit_Lagran_start$result$mea, 
    n_threads = 4L, maxit = 100L, method = "aug_Lagran", rel_eps = 1e-3, 
    maxpts = 5000L, mu = fit_Lagran_start$mu, 
    lambda = fit_Lagran_start$lambda))

system.time(
  fit_adam <- mdgc_fit(
    ptr = log_ml_ptr, vcov = start_val, mea = mdgc_obj$means, 
    n_threads = 4L, lr = 1e-3, maxit = 25L, batch_size = 100L, 
    method = "adam", rel_eps = 1e-3, maxpts = 5000L))

set.seed(fit_seed <- 19570958L)
system.time(
  fit_svrg <- mdgc_fit(
    ptr = log_ml_ptr, vcov = start_val, mea = mdgc_obj$means, 
    n_threads = 4L, lr = 1e-3, maxit = 25L, batch_size = 100L, 
    method = "svrg", verbose = TRUE, rel_eps = 1e-3, maxpts = 5000L))

# compare the log marginal likelihood 
print(rbind(
  `Augmented Lagrangian` = 
    mdgc_log_ml(vcov = fit_Lagran$result$vcov, mea = fit_Lagran$result$mea, 
                ptr = log_ml_ptr, rel_eps = 1e-3),
  ADAM = 
    mdgc_log_ml(vcov = fit_adam$result$vcov  , mea = fit_adam$result$mea, 
                ptr = log_ml_ptr, rel_eps = 1e-3),
  SVRG =
    mdgc_log_ml(vcov = fit_svrg$result$vcov  , mea = fit_svrg$result$mea, 
                ptr = log_ml_ptr, rel_eps = 1e-3),
  Truth = 
    mdgc_log_ml(vcov = dat$Sigma             , mea = numeric(5), 
                ptr = log_ml_ptr, rel_eps = 1e-3)),
  digits = 10)

# we can use an approximation in the method
set.seed(fit_seed)
system.time(
  fit_svrg_aprx <- mdgc_fit(
    ptr = log_ml_ptr, vcov = start_val, mea = mdgc_obj$means, 
    n_threads = 4L, lr = 1e-3, maxit = 25L, batch_size = 100L, 
    method = "svrg", rel_eps = 1e-3, maxpts = 5000L, use_aprx = TRUE))

# essentially the same estimates
norm(fit_svrg_aprx$result$vcov - fit_svrg$result$vcov, "F") 
sd(fit_svrg_aprx$result$mea - fit_svrg$result$mea)

# compare the estimated correlation matrix with the true value
do_plot <- function(est, truth, main){
  par_old <- par(mfcol = c(1, 3), mar  = c(1, 1, 4, 1))
  on.exit(par(par_old))
  sc <- colorRampPalette(c("Red", "White", "Blue"))(201)
  
  f <- function(x, main)
    image(x[, NCOL(x):1], main = main, col = sc, zlim = c(-1, 1), 
          xaxt = "n", yaxt = "n", bty = "n")
  f(est, main)
  f(truth, "Truth")
  f(est - truth, "Difference")
}

do_plot(fit_Lagran$result$vcov, dat$Sigma, "Estimates (Aug. Lagrangian)")
do_plot(fit_adam  $result$vcov, dat$Sigma, "Estimates (ADAM)")
do_plot(fit_svrg  $result$vcov, dat$Sigma, "Estimates (SVRG)")

norm(fit_Lagran$result$vcov - dat$Sigma, "F")
norm(fit_adam  $result$vcov - dat$Sigma, "F")
norm(fit_svrg  $result$vcov - dat$Sigma, "F")

# perform the imputation
system.time(imp_res <- mdgc_impute(
  mdgc_obj, fit_svrg$result$vcov, mea = fit_svrg$result$mea, rel_eps = 1e-3, 
  maxit = 10000L, n_threads = 4L))

# look at the result for one of the observations
imp_res[2L]

# compare with the observed and true data
rbind(truth = dat$truth_obs[2L, ], observed = dat$seen_obs[2L, ])

# we can threshold the data like this
threshold <- function(org_data, imputed){
  # checks
  stopifnot(NROW(org_data) == length(imputed), 
            is.list(imputed), is.data.frame(org_data))
  
  # threshold
  is_cont <- which(sapply(org_data, is.numeric))
  is_bin  <- which(sapply(org_data, is.logical)) 
  is_ord  <- which(sapply(org_data, is.ordered))
  stopifnot(
    length(is_cont) + length(is_bin) + length(is_ord) == NCOL(org_data))
  is_cat <- c(is_bin, is_ord)
  
  trans_to_df <- function(x){
    if(is.matrix(x))
      as.data.frame(t(x))
    else
      as.data.frame(  x )
  }
  
  out_cont <- trans_to_df(sapply(imputed, function(x) unlist(x[is_cont])))
  out_cat <- trans_to_df(sapply(imputed, function(x) 
    sapply(x[is_cat], which.max)))
  out <- cbind(out_cont, out_cat)
  
  # set factor levels etc. 
  out <- out[, order(c(is_cont, is_bin, is_ord))]
  if(length(is_bin) > 0)
    out[, is_bin] <- out[, is_bin] > 1L
  if(length(is_ord) > 0)
    for(i in is_ord)
      out[[i]] <- ordered(
        unlist(out[[i]]), labels = levels(org_data[, i]))
  
  colnames(out) <- colnames(org_data)
  out
}
thresh_dat <- threshold(dat$seen_obs, imp_res)

# compare thresholded data with observed and true data
head(thresh_dat)
head(dat$seen_obs)  # observed data
head(dat$truth_obs) # true data

# compare correct categories
get_classif_error <- function(impu_dat, truth = dat$truth_obs, 
                              observed = dat$seen_obs){
  is_cat <- sapply(truth, function(x)
    is.logical(x) || is.ordered(x))
  is_match <- impu_dat[, is_cat] == truth[, is_cat]
  is_match[!is.na(observed[, is_cat])] <- NA_integer_
  1 - colMeans(is_match, na.rm = TRUE)
}
get_classif_error(thresh_dat)

# compute RMSE
get_rmse <- function(impu_dat, truth = dat$truth_obs,
                     observed = dat$seen_obs){
  is_con <- sapply(truth, is.numeric)
  err <- as.matrix(impu_dat[, is_con] - truth[, is_con])
  err[!is.na(observed[, is_con])] <- NA_real_
  sqrt(colMeans(err^2, na.rm = TRUE))
}
get_rmse(thresh_dat)

# we can compare this with missForest
miss_forest_arg <- dat$seen_obs
is_log <- sapply(miss_forest_arg, is.logical)
miss_forest_arg[, is_log] <- lapply(miss_forest_arg[, is_log], as.factor)
set.seed(1)
system.time(miss_res <- missForest(miss_forest_arg))

# turn binary variables back to logicals
miss_res$ximp[, is_log] <- lapply(
  miss_res$ximp[, is_log], function(x) as.integer(x) > 1L)

rbind(mdgc       = get_classif_error(thresh_dat),
      missForest = get_classif_error(miss_res$ximp))
rbind(mdgc       = get_rmse(thresh_dat),
      missForest = get_rmse(miss_res$ximp))
```

## An Even Shorter Example

Here is an example where we use the `mdgc` function to do the model 
estimation and the imputation:

```{r very_quick_example, cache = 1, eval = TRUE}
# have a data set with missing continuous, binary, and ordinal variables
head(dat$seen_obs)
  
# perform the estimation and imputation
set.seed(1)
system.time(res <- mdgc(dat$seen_obs, verbose = TRUE, maxpts = 5000L, 
                        n_threads = 4L, maxit = 25L, use_aprx = TRUE))

# compare the estimated correlation matrix with the truth
norm(dat$Sigma - res$vcov, "F") / norm(dat$Sigma, "F")

# compute the classifcation error and RMSE
get_classif_error(res$ximp)
get_rmse(res$ximp)
```

We can compare this with the `mixedgcImp` which uses the method described 
in @zhao20:

```{r use_zhao19, eval = TRUE, cache = 1}
# turn the data to a format that can be based
dat_pass <- dat$seen_obs
is_cat <- sapply(dat_pass, function(x) is.logical(x) | is.ordered(x))
dat_pass[, is_cat] <- lapply(dat_pass[, is_cat], as.integer)

system.time(imp_apr_em <- impute_mixedgc(dat_pass, eps = 1e-4))

# compare the estimated correlation matrix with the truth
get_rel_err <- function(est, keep = seq_len(NROW(truth)), truth = dat$Sigma)
  norm(truth[keep, keep] - est[keep, keep], "F") / 
  norm(truth, "F")

c(mdgc                     = get_rel_err(res$vcov), 
  mixedgcImp               = get_rel_err(imp_apr_em$R), 
  `mdgc bin/ordered`       = get_rel_err(res$vcov    , is_cat),
  `mixedgcImp bin/ordered` = get_rel_err(imp_apr_em$R, is_cat),
  `mdgc continuous`        = get_rel_err(res$vcov    , !is_cat),
  `mixedgcImp continuous`  = get_rel_err(imp_apr_em$R, !is_cat))

# compare the classifcation error and RMSE
imp_apr_res <- as.data.frame(imp_apr_em$Ximp)
is_bin <- sapply(dat$seen_obs, is.logical)
imp_apr_res[, is_bin] <- lapply(imp_apr_res[, is_bin], `>`, e2 = 0)
is_ord <- sapply(dat$seen_obs, is.ordered)
imp_apr_res[, is_ord] <- mapply(function(x, idx)
  ordered(x, labels = levels(dat$seen_obs[[idx]])), 
  x = imp_apr_res[, is_ord], i = which(is_ord), SIMPLIFY = FALSE)

rbind(mdgc       = get_classif_error(res$ximp),
      mixedgcImp = get_classif_error(imp_apr_res))
rbind(mdgc       = get_rmse(res$ximp),
      mixedgcImp = get_rmse(imp_apr_res))
```

## Simulation Study

```{r before_sim_clean, echo = FALSE}
rm(list = setdiff(ls(), c(
  "dat", "get_rmse", "get_rel_err", "get_classif_error", "sim_dat", 
  "threshold", "log_ml_ptr", "mdgc_obj", "p", "do_plot")))
```

We will perform a simulation study in this section to compare different 
methods in terms of their computation time and performance. We first 
perform the simulation.

```{r sim_study, message = FALSE}
# the seeds we will use
seeds <- c(293498804L, 311878062L, 370718465L, 577520465L, 336825034L, 661670751L, 750947065L, 255824398L, 281823005L, 721186455L, 251974931L, 643568686L, 273097364L, 328663824L, 490259480L, 517126385L, 651705963L, 43381670L, 503505882L, 609251792L, 643001919L, 244401725L, 983414550L, 850590318L, 714971434L, 469416928L, 237089923L, 131313381L, 689679752L, 344399119L, 330840537L, 6287534L, 735760574L, 477353355L, 579527946L, 83409653L, 710142087L, 830103443L, 94094987L, 422058348L, 117889526L, 259750108L, 180244429L, 762680168L, 112163383L, 10802048L, 440434442L, 747282444L, 736836365L, 837303896L, 50697895L, 231661028L, 872653438L, 297024405L, 719108161L, 201103881L, 485890767L, 852715172L, 542126886L, 155221223L, 18987375L, 203133067L, 460377933L, 949381283L, 589083178L, 820719063L, 543339683L, 154667703L, 480316186L, 310795921L, 287317945L, 30587393L, 381290126L, 178269809L, 939854883L, 660119506L, 825302990L, 764135140L, 433746745L, 173637986L, 100446967L, 333304121L, 225525537L, 443031789L, 587486506L, 245392609L, 469144801L, 44073812L, 462948652L, 226692940L, 165285895L, 546908869L, 550076645L, 872290900L, 452044364L, 620131127L, 600097817L, 787537854L, 15915195L, 64220696L)

# gather or compute the results (you may skip this)
res <- lapply(seeds, function(s){
  file_name <- file.path("sim-res", sprintf("seed-%d.RDS", s))
  
  if(file.exists(file_name)){
    message(sprintf("Reading '%s'", file_name))
    out <- readRDS(file_name)
  } else { 
    message(sprintf("Running '%s'", file_name))
    
    # simulate the data
    set.seed(s)
    dat <- sim_dat(2000L, p = 15L)
    
    # fit models and impute
    mdgc_time <- system.time(
      mdgc_res <- mdgc(dat$seen_obs, verbose = FALSE, maxpts = 5000L, 
                        n_threads = 4L, maxit = 25L, use_aprx = TRUE))
    
    dat_pass <- dat$seen_obs
    is_cat <- sapply(dat_pass, function(x) is.logical(x) | is.ordered(x))
    dat_pass[, is_cat] <- lapply(dat_pass[, is_cat], as.integer)
    mixedgc_time <- 
      system.time(mixedgc_res <- impute_mixedgc(dat_pass, eps = 1e-4))
    
    miss_forest_arg <- dat$seen_obs
    is_log <- sapply(miss_forest_arg, is.logical)
    miss_forest_arg[, is_log] <- lapply(
      miss_forest_arg[, is_log], as.factor)
    sink(tempfile())
    miss_time <- system.time(
      miss_res <- missForest(miss_forest_arg, verbose = FALSE))
    sink()
    
    miss_res$ximp[, is_log] <- lapply(
      miss_res$ximp[, is_log], function(x) as.integer(x) > 1L)
    
    # impute using the other estimate
    mdgc_obj <- get_mdgc(dat$seen_obs)
    impu_mixedgc_est <- mdgc_impute(mdgc_obj, mixedgc_res$R, mdgc_obj$means)
    impu_mixedgc_est <- threshold(dat$seen_obs, impu_mixedgc_est)
    
    # gather output for the correlation matrix estimates
    vcov_res <- list(truth = dat$Sigma, mdgc = mdgc_res$vcov, 
                     mixedgc = mixedgc_res$R)
    get_rel_err <- function(est, truth, keep = seq_len(NROW(truth)))
      norm(truth[keep, keep] - est[keep, keep], "F") / norm(truth, "F")
    
    vcov_res <- within(vcov_res, {
      mdgc_rel_err    = get_rel_err(mdgc   , truth)
      mixedgc_rel_err = get_rel_err(mixedgc, truth)
    })
    
    # gather the estimated means
    mea_ests <- list(marginal = mdgc_obj$means, 
                     joint    = mdgc_res$mea)
    
    # gather output for the imputation 
    mixedgc_imp_res <- as.data.frame(mixedgc_res$Ximp)
    is_bin <- sapply(dat$seen_obs, is.logical)
    mixedgc_imp_res[, is_bin] <- 
      lapply(mixedgc_imp_res[, is_bin, drop = FALSE], `>`, e2 = 0)
    is_ord <- sapply(dat$seen_obs, is.ordered)
    mixedgc_imp_res[, is_ord] <- mapply(function(x, idx)
      ordered(x, labels = levels(dat$seen_obs[[idx]])), 
      x = mixedgc_imp_res[, is_ord, drop = FALSE], 
      i = which(is_ord), SIMPLIFY = FALSE)
    
    get_bin_err <- function(x){
      . <- function(z) z[, is_bin, drop = FALSE]
      get_classif_error(
        .(x), truth = .(dat$truth_obs), observed = .(dat$seen_obs))
    }
    get_ord_err <- function(x){
      . <- function(z) z[, is_ord, drop = FALSE]
      get_classif_error(
        .(x), truth = .(dat$truth_obs), observed = .(dat$seen_obs))
    }
          
    err <- list(
      mdgc_bin = get_bin_err(mdgc_res$ximp), 
      mixedgc_bin = get_bin_err(mixedgc_imp_res), 
      mixed_bin = get_bin_err(impu_mixedgc_est),
      missForest_bin = get_bin_err(miss_res$ximp),
      
      mdgc_class = get_ord_err(mdgc_res$ximp), 
      mixedgc_class = get_ord_err(mixedgc_imp_res), 
      mixed_class = get_ord_err(impu_mixedgc_est),
      missForest_class = get_ord_err(miss_res$ximp),
      
      mdgc_rmse = get_rmse(
        mdgc_res$ximp, truth = dat$truth_obs, observed = dat$seen_obs),
      mixedgc_rmse = get_rmse(
        mixedgc_imp_res, truth = dat$truth_obs, observed = dat$seen_obs),
      mixed_rmse = get_rmse(
        impu_mixedgc_est, truth = dat$truth_obs, observed = dat$seen_obs), 
      missForest_rmse = get_rmse(
        miss_res$ximp, truth = dat$truth_obs, observed = dat$seen_obs))
    
    # gather the times
    times <- list(mdgc = mdgc_time, mixedgc = mixedgc_time, 
                  missForest = miss_time)
    
    # save stats to check convergence
    conv_stats <- list(mdgc = mdgc_res$logLik, 
                       mixedgc = mixedgc_res$loglik)
    
    # save output 
    out <- list(vcov_res = vcov_res, err = err, times = times, 
                conv_stats = conv_stats, mea_ests = mea_ests)
    saveRDS(out, file_name)
  }
  
  # print summary stat to the console while knitting
  out <- readRDS(file_name)
  . <- function(x)
    message(paste(sprintf("%8.3f", x), collapse = " "))
  with(out, {
    message(paste(
      "mdgc    logLik", 
      paste(sprintf("%.2f", conv_stats$mdgc), collapse = " ")))
    message(paste(
      "mixedgc logLik", 
      paste(sprintf("%.2f", conv_stats$mixedgc), collapse = " ")))
    message(sprintf(
      "Relative correlation matrix estimate errors are %.4f %.4f", 
      vcov_res$mdgc_rel_err, vcov_res$mixedgc_rel_err))
    message(sprintf(
      "Times are %.2f %.2f %.2f", 
      times$mdgc["elapsed"], times$mixedgc["elapsed"], 
      times$missForest["elapsed"]))
    
    message(sprintf(
      "Binary classifcation errors are %.2f %.2f %.2f (%.2f)", 
      mean(err$mdgc_bin), mean(err$mixedgc_bin), 
      mean(err$missForest_bin), mean(err$mixed_bin)))
    
    message(sprintf(
      "Ordinal classifcation errors are %.2f %.2f %.2f (%.2f)", 
      mean(err$mdgc_class), mean(err$mixedgc_class), 
      mean(err$missForest_class), mean(err$mixed_class)))
    
    message(sprintf(
      "Mean RMSEs are %.2f %.2f %.2f (%.2f)",
      mean(err$mdgc_rmse), mean(err$mixedgc_rmse), 
      mean(err$missForest_rmse), mean(err$mixed_rmse)))
    message("")
  })
  
  out  
})
```

The difference in computation time is given below:

```{r time_diff_est}
# assign function to show the summary stats
show_sim_stats <- function(v1, v2, v3, what, sub_ele = NULL){
  vals <- sapply(res, function(x) 
    do.call(rbind, x[[what]][c(v1, v2, v3)]), 
    simplify = "array")
  if(!is.null(sub_ele))
    vals <- vals[, sub_ele, , drop = FALSE]
    
  cat("Means and standard errors:\n")
  mea_se <- function(x)
    c(mean = mean(x), SE = sd(x) / sqrt(length(x)))
  print(t(apply(vals, 1L, mea_se)))
  
  cat("\nDifference:\n")
  print(t(apply(
    c(vals[v1, , ]) - 
      aperm(vals[c(v2, v3), , , drop = FALSE], c(3L, 2L, 1L)), 
    3L, mea_se)))
}

# compare estimation time
show_sim_stats(1L, 2L, 3L, "times", "elapsed")
```

The summary stats for the relative Frobenius norm between the estimated and 
true correlation matrix is given below:

```{r F_norm_diff_est}
# relative norms
show_sim_stats("mixedgc_rel_err", "mdgc_rel_err", NULL, "vcov_res")
```

Finally, here are the results for the classification error for the binary 
and ordinal outcomes and the root mean square error:

```{r impu_diff_est}
# the binary variables
show_sim_stats("mdgc_bin", "mixedgc_bin", "missForest_bin", "err")

# the ordinal variables
show_sim_stats("mdgc_class", "mixedgc_class", "missForest_class", "err")

# the continuous variables
show_sim_stats("mdgc_rmse", "mixedgc_rmse", "missForest_rmse", "err")
```

It is important to emphasize that missForest is not estimating the true 
model.

## Adding Multinomial Variables
We extend the model suggested by @zhao20 in this section. The example is 
very similar to the previous one but with multinomial variables.

```{r mult_sim, cache = 1, fig.height = 3, fig.width = 7}
# simulates a data set and mask some of the data.
#
# Args: 
#   n: number of observations. 
#   p: number of variables. 
#   n_lvls: number of levels for the ordinal and multinomial variables. 
#   verbose: print status during the simulation.
# 
# Returns: 
#   Simulated masked data, the true data, and true covariance matrix. 
sim_dat <- function(n, p = 4L, n_lvls = 5L, verbose = FALSE){
  # determine the type
  n_rep <- floor((p + 4 - 1) / 4)
  type <- rep(1:4, n_rep)[1:p]
  is_con  <- type == 1L
  is_bin  <- type == 2L
  is_ord  <- type == 3L
  is_mult <- type == 4L
  
  col_nam <- c(outer(c("C", "B", "O", "M"), 1:n_rep, paste0))[1:p]
  idx <- head(cumsum(c(1L, ifelse(type == 4, n_lvls, 1L))), -1L)
  
  # get the covariance matrix
  n_latent <- p + (n_lvls - 1L) * (p %/% 4)
  Sig <- drop(rWishart(1L, 2 * n_latent, diag(1 / n_latent / 2, n_latent)))
  
  # essentially set the reference level to zero
  for(i in idx[is_mult]){
    Sig[i,  ] <- 0
    Sig[ , i] <- 0
  }
  
  # rescale some rows and columns
  sds <- sqrt(diag(Sig))
  for(i in idx[is_mult]){
    sds[i] <- 1
    sds[i + 3:n_lvls - 1] <- 1
  }
  Sig <- diag(1/sds) %*% Sig %*% diag(1/sds)
      
  # draw the observations
  truth <- mvtnorm::rmvnorm(n, sigma = Sig)
  truth[, idx[is_mult]] <- 0
  
  # sample which are masked data 
  is_mask <- matrix(runif(n * p) < .3, n)
  
  # make sure we have no rows with all missing data
  while(any(all_nans <- rowSums(is_mask) == NCOL(is_mask)))
    is_mask[all_nans, ] <- runif(sum(all_nans) * p) < .3
  
  # create the observed data
  truth_obs <- lapply(type, function(i) if(i == 1L) numeric(n) else integer(n))
  truth_obs <- data.frame(truth_obs)
  colnames(truth_obs) <- col_nam
  
  bs_ord <- qnorm(seq(0, 1, length.out = n_lvls + 1L))
  for(i in 1:p){
    idx_i <- idx[i]
    switch(
      type[i],
      # continuous
      truth_obs[, i] <- qexp(pnorm(truth[, idx_i])),
      # binary
      truth_obs[, i] <- truth[, idx_i] > 0,
      # ordinal
      {
        truth_obs[, i] <- 
          ordered(as.integer(cut(truth[, idx_i], breaks = bs_ord)))
        levels(truth_obs[, i]) <- 
          LETTERS[seq_len(length(unique(truth_obs[, i])))]
      },
      # multinomial
      {
        truth_obs[, i] <- apply(
          truth[, idx_i + 1:n_lvls - 1L], 1L, which.max)
        truth_obs[, i] <- factor(truth_obs[, i], 
                                 labels = paste0("T", 1:n_lvls))
      }, 
      stop("Type is not implemented"))
  }

  # mask the data
  seen_obs <- truth_obs
  seen_obs[is_mask] <- NA
  
  list(truth = truth, truth_obs = truth_obs, seen_obs = seen_obs, 
       Sigma = Sig)
}

# simulate and show the data
set.seed(1)
p <- 8L
dat <- sim_dat(2000L, p = p, verbose = TRUE, n_lvls = 4)

# show the first rows of the observed data
head(dat$seen_obs)

# assign object to perform the estimation and the imputation
obj <- get_mdgc(dat$seen_obs)
ptr <- get_mdgc_log_ml(obj)

# get starting values
start_vals <- mdgc_start_value(obj)

# plot the starting values and the true values
do_plot <- function(est, truth, main){
  par_old <- par(mfcol = c(1, 3), mar  = c(1, 1, 4, 1))
  on.exit(par(par_old))
  sc <- colorRampPalette(c("Red", "White", "Blue"))(201)

  ma <- max(abs(est), max(abs(truth)))  
  f <- function(x, main)
    image(x[, NCOL(x):1], main = main, col = sc, zlim = c(-ma, ma), 
          xaxt = "n", yaxt = "n", bty = "n")
  f(est, main)
  f(truth, "Truth")
  f(est - truth, "Difference")
}

do_plot(start_vals, dat$Sigma, "Starting values")
# check the log marginal likelihood at the starting values and compare with
# the true values at the starting values
mdgc_log_ml(ptr, start_vals, mea = obj$means, n_threads = 1L)
# and at the true values
mdgc_log_ml(ptr, dat$Sigma , mea = numeric(length(obj$means)), 
            n_threads = 1L)

# much better than using a diagonal matrix!
mdgc_log_ml(ptr, diag(NROW(dat$Sigma)), mea = obj$means, n_threads = 1L)

# estimate the model
system.time(
  ests <- mdgc_fit(ptr, vcov = start_vals, mea = obj$means, 
                   method = "aug_Lagran",
                   n_threads = 4L, rel_eps = 1e-2, maxpts = 1000L, 
                   minvls = 200L, use_aprx = TRUE, conv_crit = 1e-8))

# refine the estimates
system.time(
  ests <- mdgc_fit(ptr, vcov = ests$result$vcov, 
                   mea = ests$result$mea, 
                   method = "aug_Lagran",
                   n_threads = 4L, rel_eps = 1e-3, maxpts = 10000L, 
                   minvls = 1000L, mu = ests$mu, lambda = ests$lambda, 
                   use_aprx = TRUE, conv_crit = 1e-8))

# use ADAM
system.time(
  fit_adam <- mdgc_fit(
    ptr, vcov = start_vals, mea = obj$means, minvls = 200L,
    n_threads = 4L, lr = 1e-3, maxit = 25L, batch_size = 100L, 
    method = "adam", rel_eps = 1e-3, maxpts = 5000L, 
    use_aprx = TRUE))

# use SVRG
system.time(
  fit_svrg <- mdgc_fit(
    ptr, vcov = start_vals, mea = obj$means,  minvls = 200L,
    n_threads = 4L, lr = 1e-3, maxit = 25L, batch_size = 100L, 
    method = "svrg", verbose = TRUE, rel_eps = 1e-3, maxpts = 5000L, 
    use_aprx = TRUE, conv_crit = 1e-8))

# compare log marginal likelihood
print(rbind(
  `Augmented Lagrangian` = 
    mdgc_log_ml(ptr, ests$result$vcov    , mea = ests$result$mea, 
                n_threads = 1L),
  ADAM = 
    mdgc_log_ml(ptr, fit_adam$result$vcov, mea = fit_adam$result$mea, 
                n_threads = 1L),
  SVRG = 
    mdgc_log_ml(ptr, fit_svrg$result$vcov, mea = fit_svrg$result$mea, 
                n_threads = 1L),
  Truth = 
    mdgc_log_ml(ptr, dat$Sigma           , mea = numeric(length(obj$means)), 
                n_threads = 1L)), digits = 10)

# compare the estimated and the true values (should not match because of
# overparameterization? See https://stats.stackexchange.com/q/504682/81865)
do_plot(ests$result$vcov    , dat$Sigma, "Estimates (Aug. Lagrangian)")
do_plot(fit_adam$result$vcov, dat$Sigma, "Estimates (ADAM)")
do_plot(fit_svrg$result$vcov, dat$Sigma, "Estimates (SVRG)")
# after rescaling
do_plot_rescale <- function(x, lab){
  trans <- function(z){
    scal <- diag(NCOL(z))
    m <- obj$multinomial[[1L]]
    for(i in seq_len(NCOL(m))){
      idx <- m[3, i] + 1 + seq_len(m[2, i] - 1)
      scal[idx, idx] <- solve(t(chol(z[idx, idx])))
    }
    tcrossprod(scal %*% z, scal)
  }
  
  do_plot(trans(x), trans(dat$Sigma), lab)  
}
do_plot_rescale(ests$result$vcov    , "Estimates (Aug. Lagrangian)")
do_plot_rescale(fit_adam$result$vcov, "Estimates (ADAM)")
do_plot_rescale(fit_svrg$result$vcov, "Estimates (SVRG)")
# perform the imputation
system.time(
  imp_res <- mdgc_impute(obj, ests$result$vcov, mea = ests$result$mea, 
                         rel_eps = 1e-3, maxit = 10000L, n_threads = 4L))

# look at the result for one of the observations
imp_res[1L]

# compare with the observed and true data
rbind(truth = dat$truth_obs[1L, ], observed = dat$seen_obs[1L, ])

# we can threshold the data like this
threshold <- function(org_data, imputed){
  # checks
  stopifnot(NROW(org_data) == length(imputed), 
            is.list(imputed), is.data.frame(org_data))
  
  # threshold
  is_cont <- which(sapply(org_data, is.numeric))
  is_bin  <- which(sapply(org_data, is.logical)) 
  is_ord  <- which(sapply(org_data, is.ordered))
  is_mult <- which(sapply(org_data, is.factor))
  is_mult <- setdiff(is_mult, is_ord)
  stopifnot(
    length(is_cont) + length(is_bin) + length(is_ord) + length(is_mult) == 
      NCOL(org_data))
  is_cat <- c(is_bin, is_ord, is_mult)
  
  trans_to_df <- function(x){
    if(is.matrix(x))
      as.data.frame(t(x))
    else
      as.data.frame(  x )
  }
  
  out_cont <- trans_to_df(sapply(imputed, function(x) unlist(x[is_cont])))
  out_cat <- trans_to_df(sapply(imputed, function(x) 
    sapply(x[is_cat], which.max)))
  out <- cbind(out_cont, out_cat)
  
  # set factor levels etc. 
  out <- out[, order(c(is_cont, is_bin, is_ord, is_mult))]
  if(length(is_bin) > 0)
    out[, is_bin] <- out[, is_bin] > 1L
  if(length(is_ord) > 0)
    for(i in is_ord)
      out[[i]] <- ordered(
        unlist(out[[i]]), labels = levels(org_data[, i]))
  if(length(is_mult) > 0)
    for(i in is_mult)
      out[[i]] <- factor(
        unlist(out[[i]]), labels = levels(org_data[, i]))
  
  colnames(out) <- colnames(org_data)
  out
}
thresh_dat <- threshold(dat$seen_obs, imp_res)

# compare thresholded data with observed and true data
head(thresh_dat)
head(dat$seen_obs)  # observed data
head(dat$truth_obs) # true data

# compare correct categories
get_classif_error <- function(impu_dat, truth = dat$truth_obs, 
                              observed = dat$seen_obs){
  is_cat <- sapply(truth, function(x)
    is.logical(x) || is.factor(x))
  is_match <- impu_dat[, is_cat] == truth[, is_cat]
  is_match <- matrix(is_match, ncol = sum(is_cat))
  is_match[!is.na(observed[, is_cat])] <- NA_integer_
  setNames(1 - colMeans(is_match, na.rm = TRUE), 
           colnames(truth)[is_cat])
}
get_classif_error(thresh_dat)

# compute RMSE
get_rmse <- function(impu_dat, truth = dat$truth_obs,
                     observed = dat$seen_obs){
  is_con <- sapply(truth, is.numeric)
  err <- as.matrix(impu_dat[, is_con] - truth[, is_con])
  err[!is.na(observed[, is_con])] <- NA_real_
  sqrt(colMeans(err^2, na.rm = TRUE))
}
get_rmse(thresh_dat)

# we can compare this with missForest
miss_forest_arg <- dat$seen_obs
is_log <- sapply(miss_forest_arg, is.logical)
miss_forest_arg[, is_log] <- lapply(miss_forest_arg[, is_log], as.factor)
set.seed(1)
system.time(miss_res <- missForest(miss_forest_arg))

# turn binary variables back to logical variables
miss_res$ximp[, is_log] <- lapply(
  miss_res$ximp[, is_log], function(x) as.integer(x) > 1L)

# compare errors
rbind(mdgc       = get_classif_error(thresh_dat),
      missForest = get_classif_error(miss_res$ximp))
rbind(mdgc       = get_rmse(thresh_dat),
      missForest = get_rmse(miss_res$ximp))
```

### Edgar Anderson's Iris Data
We make a small example below were we take the iris data set and randomly 
mask it. Then we compare the imputation method in this package with 
missForest.

```{r comp_w_iris, cache=1}
# re-scales continuous variables to have scale 1.
#
# Args:
#   dat: data to rescale.
re_scale <- function(dat){
  is_num <- sapply(dat, is.numeric)
  if(!any(is_num))
    return(dat)
  dat[is_num] <- lapply(dat[is_num], scale)
  dat[is_num] <- lapply(dat[is_num], c)
  dat
}

# load the iris data set
data(iris)
iris <- re_scale(iris)

# assign function to produce iris data set with NAs.
# 
# Args:
#   dat: data set to mask.
#   p_na: chance of missing a value.
mask_dat <- function(dat, p_na = .3){
  is_miss <- matrix(p_na > runif(NROW(dat) * NCOL(dat)), 
                    NROW(dat), NCOL(dat))
  while(any(all_missing <- apply(is_miss, 1, all)))
    # avoid rows with all missing variables
    is_miss[all_missing, ] <- p_na > runif(sum(all_missing) * NCOL(dat))
  
  # create data set with missing values
  out <- dat
  out[is_miss] <- NA 
  out
}

# get a data set with all missing values
set.seed(68129371)
dat <- mask_dat(iris)

# use the mdgc method
system.time(
  mdgc_res <- mdgc(dat, maxpts = 10000L, minvls = 500L, n_threads = 4L, 
                   maxit = 50L, use_aprx = TRUE, conv_crit = 1e-8, 
                   method = "svrg", rel_eps = 1e-2, batch_size = 100L, 
                   iminvls = 2000L, imaxit = 20000L, irel_eps = 1e-3, 
                   lr = 1e-3))

# some of the imputed values
head(mdgc_res$ximp)

# compare with missForest
system.time(miss_res <- missForest(dat))

# the errors
rbind(
  mdgc = get_classif_error(
    impu_dat = mdgc_res$ximp, truth = iris, observed = dat), 
  missForest = get_classif_error(
    impu_dat = miss_res$ximp, truth = iris, observed = dat))

rbind(
  mdgc = get_rmse(
    impu_dat = mdgc_res$ximp, truth = iris, observed = dat), 
  missForest = 
    get_rmse(impu_dat = miss_res$ximp, truth = iris, observed = dat))
```

We repeat this a few times to get Monte Carlo estimates of the errors:

```{r mult_iris, message=FALSE, cache=1}
# function to get Monte Carlo estimates of the errors. 
# 
# Args:
#   dat: data set to use. 
#   seeds: seeds to use.
get_err_n_time <- function(dat, seeds){
  cl <- makeCluster(4L)
  registerDoParallel(cl)
  on.exit(stopCluster(cl))
    
  sapply(seeds, function(s){
    # mask data
    set.seed(s)
    dat_mask <- mask_dat(dat)
    
    # fit models
    mdgc_time <- system.time(
      mdgc_res <- mdgc(dat_mask, maxpts = 10000L, minvls = 500L, n_threads = 4L, 
                       maxit = 50L, use_aprx = TRUE, conv_crit = 1e-8, 
                       method = "svrg", rel_eps = 1e-2, batch_size = 100L, 
                       iminvls = 2000L, imaxit = 20000L, irel_eps = 1e-3, 
                       lr = 1e-3))
    
    # compare with missForest
    miss_forest_arg <- dat_mask
    is_log <- sapply(miss_forest_arg, is.logical)
    miss_forest_arg[, is_log] <- lapply(miss_forest_arg[, is_log], as.factor)
    sink(tempfile())
    miss_time <- system.time(miss_res <- missForest(
      miss_forest_arg, parallelize = "forests"))
    sink()
    
    # turn binary variables back to logicals
    miss_res$ximp[, is_log] <- lapply(
      miss_res$ximp[, is_log], function(x) as.integer(x) > 1L)
    
    # get the errors
    er_int <- rbind(
      mdgc = get_classif_error(
        impu_dat = mdgc_res$ximp, truth = dat, observed = dat_mask), 
      missForest = get_classif_error(
        impu_dat = miss_res$ximp, truth = dat, observed = dat_mask))
    er_con <- rbind(
      mdgc = get_rmse(
        impu_dat = mdgc_res$ximp, truth = dat, observed = dat_mask), 
      missForest = 
        get_rmse(impu_dat = miss_res$ximp, truth = dat, observed = dat_mask))
      
    # gather the output and return 
    er <- cbind(er_int, er_con)
    er <- er[, match(colnames(er), colnames(dat))]
    ti <- rbind(mdgc_time, miss_time)[, 1:3]
    out <- cbind(er, ti)
    message(sprintf("\nResult with seed %d is:", s))
    message(paste0(capture.output(print(out)), collapse = "\n"))
    out

  }, simplify = "array")
}

# get the results
seeds <- c(40574428L, 13927943L, 31430660L, 38396447L, 20137114L, 59492953L, 
           93913797L, 95452931L, 77261969L, 10996196L)
res <- get_err_n_time(iris, seeds)
```

```{r show_res_mult_iris}
# compute means and Monte Carlo standard errors
show_res <- function(res, dat){
  stats <- apply(res, 1:2, function(x) 
    c(mean = mean(x), SE = sd(x) / sqrt(length(x))))
  stats <- stats[, , c(colnames(dat), "user.self", "elapsed")]
  for(i in seq_len(dim(stats)[[3]])){
    nam <- dimnames(stats)[[3]][i]
    cap <- if(nam %in% colnames(dat)){
      if(is.ordered(dat[[nam]]))
        "ordinal"
      else if(is.factor(dat[[nam]]))
        "multinomial"
      else if(is.logical(dat[[nam]]))
        "binary"
      else
        "continuous"
    } else
      "computation time"
    
    cat(sprintf("\n%s (%s):\n", nam, cap))
    print(apply(round(
      stats[, , i], 4), 2, function(x) sprintf("%.4f (%.4f)", x[1], x[2])), 
      quote = FALSE)
  }
}

show_res(res, iris)
```

### Chemotherapy for Stage B/C Colon Cancer
We do as in the [Edgar Anderson's Iris Data](#edgar-andersons-iris-data)
section here but with a different data set.

```{r colon, cache = 1}
# prepare the data
library(survival)
colon_use <- colon[, setdiff(
  colnames(colon), c("id", "study", "time", "status", "node4", "etype"))]
colon_use <- within(colon_use, {
  sex <- sex > 0
  obstruct <- obstruct > 0
  perfor <- perfor > 0
  adhere <- adhere > 0
  differ <- ordered(differ)
  extent <- ordered(extent)
  surg <- surg > 0
})
colon_use <- colon_use[complete.cases(colon_use), ]
colon_use <- re_scale(colon_use)

# stats for the data set
summary(colon_use)

# sample missing values
set.seed(68129371)
dat <- mask_dat(colon_use)

# use the mdgc method
system.time(
  mdgc_res <- mdgc(dat, maxpts = 10000L, minvls = 500L, n_threads = 4L, 
                   maxit = 50L, use_aprx = TRUE, conv_crit = 1e-8, 
                   method = "svrg", rel_eps = 1e-2, batch_size = 100L, 
                   iminvls = 2000L, imaxit = 20000L, irel_eps = 1e-3, 
                   lr = 1e-3))

# some of the imputed values
head(mdgc_res$ximp)

# compare with missForest
miss_forest_arg <- dat
is_log <- sapply(miss_forest_arg, is.logical)
miss_forest_arg[, is_log] <- lapply(miss_forest_arg[, is_log], as.factor)
system.time(miss_res <- missForest(miss_forest_arg))

# turn binary variables back to logicals
miss_res$ximp[, is_log] <- lapply(
  miss_res$ximp[, is_log], function(x) as.integer(x) > 1L)

# the errors
rbind(
  mdgc = get_classif_error(
    impu_dat = mdgc_res$ximp, truth = colon_use, observed = dat), 
  missForest = get_classif_error(
    impu_dat = miss_res$ximp, truth = colon_use, observed = dat))

rbind(
  mdgc = get_rmse(
    impu_dat = mdgc_res$ximp, truth = colon_use, observed = dat), 
  missForest = get_rmse(
    impu_dat = miss_res$ximp, truth = colon_use, observed = dat))
```

```{r mult_colon, message=FALSE, cache=1}
# get the results
res <- get_err_n_time(colon_use, seeds)
```

```{r show_res_mult_colon}
# compute means and Monte Carlo standard errors
show_res(res, colon_use)
```

### Cholesterol Data From a US Survey
We do as in the [Edgar Anderson's Iris Data](#edgar-andersons-iris-data)
section here but with a different data set.

```{r nhanes, cache = 1}
# prepare the data
data("nhanes", package = "survey")

nhanes_use <- within(nhanes, {
  HI_CHOL <- HI_CHOL > 0
  race <- factor(race)
  agecat <- ordered(agecat)
  RIAGENDR <- RIAGENDR > 1
})[, c("HI_CHOL", "race", "agecat", "RIAGENDR")]
nhanes_use <- nhanes_use[complete.cases(nhanes_use), ]
nhanes_use <- re_scale(nhanes_use)

# summary stats for the data
summary(nhanes_use)

# sample a data set
set.seed(1)
dat <- mask_dat(nhanes_use)

# use the mdgc method
system.time(
  mdgc_res <- mdgc(dat, maxpts = 10000L, minvls = 500L, n_threads = 4L, 
                   maxit = 50L, use_aprx = TRUE, conv_crit = 1e-8, 
                   method = "svrg", rel_eps = 1e-2, batch_size = 100L, 
                   iminvls = 2000L, imaxit = 20000L, irel_eps = 1e-3, 
                   lr = 1e-3))

# some of the imputed values
head(mdgc_res$ximp)

# compare with missForest
miss_forest_arg <- dat
is_log <- sapply(miss_forest_arg, is.logical)
miss_forest_arg[, is_log] <- lapply(miss_forest_arg[, is_log], as.factor)
system.time(miss_res <- missForest(miss_forest_arg))

# turn binary variables back to logicals
miss_res$ximp[, is_log] <- lapply(
  miss_res$ximp[, is_log], function(x) as.integer(x) > 1L)

# the errors
rbind(
  mdgc = get_classif_error(
    impu_dat = mdgc_res$ximp, truth = nhanes_use, observed = dat), 
  missForest = get_classif_error(
    impu_dat = miss_res$ximp, truth = nhanes_use, observed = dat))
```

```{r mult_nhanes, message=FALSE, cache=1}
# get the results
res <- get_err_n_time(nhanes_use, seeds)
```

```{r show_res_mult_nhanes}
# compute means and Monte Carlo standard errors
show_res(res, nhanes_use)
```

## References