## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5, dpi = 150, out.width = "100%" ) ## ----library------------------------------------------------------------------ library(proxymix) ## ----engines------------------------------------------------------------------ has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE) has_mice <- requireNamespace("mice", quietly = TRUE) ## ----stored-results, include = FALSE------------------------------------------ ## The comparison table reads stored simulation results. They must come ## from the same major.minor version of proxymix as this build. res <- readRDS("results/missing_data.rds") major_minor <- function(v) paste(unlist(package_version(v))[1:2], collapse = ".") if (major_minor(res$proxymix_version) != major_minor(as.character(packageVersion("proxymix")))) { stop("results/missing_data.rds was built under proxymix ", res$proxymix_version, ", but this is proxymix ", packageVersion("proxymix"), ". Rerun the simulation and ", "data-raw/vignette_results/missing_data.R.", call. = FALSE) } ## Small numbers are written as plain decimals rather than in the ## scientific notation that knitr's inline hook would otherwise use. fixed <- function(v, digits) { format(round(v, digits), nsmall = digits, scientific = FALSE) } ## ----data--------------------------------------------------------------------- set.seed(20260620) n <- 600L lab <- sample(c(-1, 1), n, replace = TRUE) x1 <- 2 * lab + rnorm(n, 0, 0.6) x2 <- 2 * lab + 0.5 * (x1 - 2 * lab) + rnorm(n, 0, 0.6) truth <- cbind(x1 = x1, x2 = x2) x_holes <- truth missing <- runif(n) < plogis(0.6 * x1) x_holes[missing, "x2"] <- NA frac_missing <- mean(missing) ## ----gap-truth---------------------------------------------------------------- in_gap <- function(v) mean(abs(v) < 1) gap_truth <- in_gap(truth[missing, "x2"]) ## ----impute------------------------------------------------------------------- imp <- gmm_impute(x_holes, N = 2L, m = 20L, seed = 1L) imp ## ----complete----------------------------------------------------------------- done <- gmm_complete(imp, 1L) anyNA(done) ## ----single------------------------------------------------------------------- imp1 <- gmm_impute(x_holes, N = 1L, m = 20L, seed = 1L) done1 <- gmm_complete(imp1, 1L) ## ----gap-table, echo = FALSE-------------------------------------------------- gap_tbl <- data.frame( source = c("deleted values (truth)", "mixture imputation (N = 2)", "single-normal imputation (N = 1)"), gap_share = c(gap_truth, in_gap(done[missing, "x2"]), in_gap(done1[missing, "x2"])), stringsAsFactors = FALSE ) knitr::kable( gap_tbl, digits = 3L, col.names = c("Values", "Share with $\\lvert x_2 \\rvert < 1$"), caption = paste0( "Share of values in the empty middle, over the ", sum(missing), " deleted entries. The imputation rows use the first completed ", "dataset of each imputation." ) ) ## ----fig-modes, eval = has_ggplot2, echo = has_ggplot2, fig.height = 3.8, fig.cap = "Density of the deleted values of $x_2$ and of the values each imputation put in their place, from the first completed dataset of each. The shaded band is the empty middle, $|x_2| < 1$.", fig.alt = "Three density curves over x2: the deleted values, the mixture imputation and the single-normal imputation. All three have two peaks. The single-normal curve has lower peaks and more density in the shaded band between the two groups."---- dens_df <- function(v, label) { d <- density(v) data.frame(x2 = d$x, density = d$y, source = label, stringsAsFactors = FALSE) } plot_df <- rbind( dens_df(truth[missing, "x2"], "deleted values (truth)"), dens_df(done[missing, "x2"], "mixture (N = 2)"), dens_df(done1[missing, "x2"], "single normal (N = 1)") ) plot_df$source <- factor( plot_df$source, levels = c("deleted values (truth)", "mixture (N = 2)", "single normal (N = 1)") ) ggplot2::ggplot(plot_df, ggplot2::aes(x2, density, colour = source)) + ggplot2::geom_line(linewidth = 0.9) + ggplot2::annotate("rect", xmin = -1, xmax = 1, ymin = -Inf, ymax = Inf, fill = "grey60", alpha = 0.18) + ggplot2::scale_colour_manual( name = NULL, values = c("deleted values (truth)" = "#000000", "mixture (N = 2)" = "#009E73", "single normal (N = 1)" = "#D55E00") ) + ggplot2::labs( x = expression(x[2]), y = "density", title = "Deleted values and the values imputed in their place", subtitle = expression("Shaded band: the empty middle, " * group("|", x[2], "|") < 1) ) + ggplot2::theme_minimal(base_size = 11) + ggplot2::theme(legend.position = "top") ## ----fig-modes-skip, eval = !has_ggplot2, echo = FALSE, results = "asis"------ # cat("ggplot2 is not installed on this build, so the imputation-density", # "figure is skipped.\n") ## ----pool-mean---------------------------------------------------------------- pooled <- proxy_pool(imp, "x2") pooled1 <- proxy_pool(imp1, "x2") fmi_mix <- proxy_fmi(imp, "x2") mean_truth <- mean(truth[, "x2"]) err_mix <- abs(pooled$estimate - mean_truth) err_one <- abs(pooled1$estimate - mean_truth) ## ----pool-table, echo = FALSE------------------------------------------------- pool_tbl <- data.frame( model = c("complete data, before deletion", "rows not deleted", "mixture imputation (N = 2)", "single-normal imputation (N = 1)"), estimate = c(mean(truth[, "x2"]), mean(x_holes[!missing, "x2"]), pooled$estimate, pooled1$estimate), std_error = c(NA_real_, NA_real_, pooled$std.error, pooled1$std.error), conf_low = c(NA_real_, NA_real_, pooled$conf.low, pooled1$conf.low), conf_high = c(NA_real_, NA_real_, pooled$conf.high, pooled1$conf.high), stringsAsFactors = FALSE ) pool_tbl$abs_error <- abs(pool_tbl$estimate - mean(truth[, "x2"])) old_opt <- options(knitr.kable.NA = "") pool_out <- knitr::kable( pool_tbl, digits = 4L, col.names = c("Data used", "Mean of $x_2$", "SE", "CI lower", "CI upper", "Distance from complete data"), caption = paste( "The mean of $x_2$ from the complete data, from the rows that were", "not deleted, and pooled over each imputation." ) ) options(old_opt) pool_out ## ----pool-mice, eval = has_mice, echo = has_mice------------------------------ mice_fit <- mice::pool(with(as_mids(imp), lm(x2 ~ x1))) ## ----pool-mice-table, eval = has_mice, echo = FALSE--------------------------- mice_tbl <- summary(mice_fit) mice_tbl$p.value <- ifelse(mice_tbl$p.value < 0.001, "< 0.001", format(round(mice_tbl$p.value, 3L), nsmall = 3L)) knitr::kable( mice_tbl, digits = 3L, align = c("l", "r", "r", "r", "r", "r"), caption = paste0( "The regression of $x_2$ on $x_1$, fitted in each of the ", imp@m, " completed datasets of the mixture imputation and pooled by ", "`mice::pool()`." ) ) ## ----pool-mice-skip, eval = !has_mice, echo = FALSE, results = "asis"--------- # cat("mice is not installed on this build, so the pooled regression is", # "skipped. The pooled mean above does not need mice.\n") ## ----compare-facts, include = FALSE------------------------------------------- sim_value <- function(design, estimand, method, what) { s1 <- res$sim_tab$design == design & res$sim_tab$estimand == estimand & res$sim_tab$method == method res$sim_tab[[what]][s1] } cov_b <- function(method) fixed(sim_value("B", "slope", method, "coverage"), 3) cov_a <- function(method) fixed(sim_value("A", "slope", method, "coverage"), 3) ## ----compare-table, echo = FALSE---------------------------------------------- methods <- c("complete data", "proxymix", "mice", "Amelia") cmp_tbl <- data.frame( method = c("complete data, before deletion", "proxymix", "mice", "Amelia"), cov_a = vapply(methods, function(s1) { sim_value("A", "slope", s1, "coverage") }, numeric(1L)), cov_b = vapply(methods, function(s1) { sim_value("B", "slope", s1, "coverage") }, numeric(1L)), rmse_b = vapply(methods, function(s1) { sim_value("B", "slope", s1, "rmse") }, numeric(1L)), width_b = vapply(methods, function(s1) { sim_value("B", "slope", s1, "width") }, numeric(1L)), stringsAsFactors = FALSE ) knitr::kable( cmp_tbl, digits = 3L, row.names = FALSE, align = c("l", "r", "r", "r", "r"), col.names = c("Data used", "Coverage, one cloud", "Coverage, two groups", "Error, two groups", "Interval width, two groups"), caption = paste0( "Slope of $x_2$ on $x_1$ over ", res$n_rep, " simulated datasets per ", "design. Coverage is the share of 95% intervals that contained the ", "true slope. Error is the root mean squared error of the estimate. ", "With ", res$n_rep, " datasets, a coverage near 0.95 has a simulation ", "standard error of about ", fixed(sqrt(0.95 * 0.05 / res$n_rep), 2L), "." ) ) ## ----compare-code, eval = FALSE----------------------------------------------- # library(proxymix) # library(mice) # library(Amelia) # # # one dataset of 500 rows with two groups # set.seed(1L) # n <- 500L # grp <- runif(n) < 0.5 # rho <- ifelse(grp, -0.3, 0.6) # z1 <- rnorm(n) # full <- data.frame( # x1 = ifelse(grp, 1.5, -1.5) + z1, # x2 = ifelse(grp, 2, -1.5) + rho * z1 + sqrt(1 - rho^2) * rnorm(n) # ) # # # delete x2 with a probability that rises with x1 # obs <- full # obs$x2[runif(n) < plogis(0.4 + 0.6 * full$x1)] <- NA # # # 20 completed datasets from each package # sets <- list( # proxymix = complete(as_mids(gmm_impute(obs, m = 20L, seed = 1L)), "all"), # mice = complete(mice(obs, m = 20L, seed = 1L, printFlag = FALSE), "all"), # Amelia = amelia(obs, m = 20L, p2s = 0L)$imputations # ) # # # the same regression in every completed dataset, pooled by mice::pool() # lapply(sets, function(s) { # fits <- lapply(s, function(d) lm(x2 ~ x1, data = d)) # summary(pool(fits), conf.int = TRUE) # }) ## ----gap-all, include = FALSE------------------------------------------------- ## share in the empty middle, averaged over every completed dataset gap_over <- function(im) { mean(vapply(seq_len(im@m), function(i1) { in_gap(gmm_complete(im, i1)[missing, "x2"]) }, numeric(1L))) } gap_all <- c(mixture = gap_over(imp), single = gap_over(imp1)) ## line and spread along which each imputation model draws x2 given x1, ## from the covariance matrix of each fitted component model_line <- function(fit) { covs <- fit@covariances list(slope = vapply(covs, function(v) v[2L, 1L] / v[1L, 1L], numeric(1L)), sd = vapply(covs, function(v) sqrt(v[2L, 2L] - v[2L, 1L]^2 / v[1L, 1L]), numeric(1L))) } model_mix <- model_line(imp@point_fit) model_one <- model_line(imp1@point_fit) fit_within <- lm(x2 ~ x1 + factor(lab), data = as.data.frame(truth)) slope_within <- coef(fit_within)[["x1"]] sd_within <- sd(resid(fit_within)) ## pooled estimate of the analyst's regression = mean over completions pooled_slope <- function(im) { mean(vapply(seq_len(im@m), function(i1) { coef(lm(x2 ~ x1, data = as.data.frame(gmm_complete(im, i1))))[["x1"]] }, numeric(1L))) } fit_slope <- c(mixture = pooled_slope(imp), single = pooled_slope(imp1)) slope_complete <- coef(lm(x2 ~ x1, data = as.data.frame(truth)))[["x1"]] ## ----session-info, collapse = FALSE, class.output = "session-info"------------ sessionInfo()