Addressing the problem
Data in which the large values go missing
The data have 600 rows and two columns, x1 and
y. They are drawn from a mixture of two normal
distributions of equal size. In each, x1 and y
have unit variance and a correlation of 0.6. The chance that
y is deleted is plogis(-0.5 + 0.7 * y). Larger
values of y are therefore deleted more often. The complete
data are kept as the truth.
set.seed(20260622)
n <- 600L
comp <- sample(1:2, n, replace = TRUE)
mu <- rbind(c(0, 0), c(1.5, 0.5))
chol_r <- chol(matrix(c(1, 0.6, 0.6, 1), 2L))
z_full <- matrix(rnorm(2 * n), n, 2L) %*% chol_r + mu[comp, ]
colnames(z_full) <- c("x1", "y")
truth <- mean(z_full[, 2L])
beta_true <- 0.7
miss <- runif(n) < plogis(-0.5 + beta_true * z_full[, 2L])
dat <- z_full
dat[miss, "y"] <- NAImpute under two assumptions
The first imputation assumes missing at random. The second supplies
the slope that generated the data, beta = 0.7, to
mnar(). Both make 5 completed datasets. Both allow up to
500 fitting rounds (max_iter), the default of
proxy_mnar_sensitivity() below. Near missing at random, the
fitting can take more rounds than the gmm_impute() default
of 100.
m_draws <- 5L
mar_fit <- gmm_impute(dat, N = 2L, m = m_draws, mechanism = mar(),
seed = 1L, max_iter = 500L)
mnar_fit <- gmm_impute(dat, N = 2L, m = m_draws,
mechanism = mnar("y", beta = beta_true), seed = 1L,
max_iter = 500L)
mar_est <- proxy_pool(mar_fit, "y")$estimate
mnar_est <- proxy_pool(mnar_fit, "y", method = "rubin")$estimateFor data missing at random, proxy_pool() computes the
pooled standard error exactly from the fitted mixture. That exact
formula does not hold under censored() or
mnar(). These fits are pooled by Rubin’s rules instead
(Rubin, 1987). These rules average the estimates over the completed
datasets and add their spread to the average variance. Called without
method = "rubin", proxy_pool() switches to
Rubin’s rules for such fits and prints a message.
| Data used | Mean of y | Difference from complete data |
|---|---|---|
| complete data, before deletion | 0.285 | 0.000 |
| rows not deleted | -0.044 | -0.330 |
| imputed, missing at random | 0.061 | -0.224 |
| imputed, missing not at random (slope 0.7) | 0.244 | -0.041 |
Sweep the assumed slope
In real data the slope is unknown. Under this model, the observed
data carry some information about the slope, but only through the
assumed shape of the distribution of y. A different shape
could fit the observed data equally well with a different slope. An
estimate at a single slope would rest on that shape.
proxy_mnar_sensitivity() therefore repeats the
imputation at each slope in a grid. For each slope, it returns the
pooled mean with its 95% confidence interval (CI), the log-likelihood,
and whether the fit converged. The log-likelihood measures how well the
model with that slope fits the observed data. Higher values mean a
better fit. A fit has converged when its rounds settle within the limit
of 500.
sweep <- proxy_mnar_sensitivity(dat, "y",
beta_grid = seq(0, 1.2, by = 0.3),
N = 2L, m = m_draws, seed = 1L)
covers <- sweep$conf.low <= truth & sweep$conf.high >= truth
beta_first <- sweep$beta[min(which(covers))]
beta_best <- sweep$beta[which.max(sweep$loglik)]
ll_gain <- max(sweep$loglik) - sweep$loglik[1L]| Assumed slope | Pooled mean of y | CI lower | CI upper | Log-likelihood | Converged |
|---|---|---|---|---|---|
| 0.0 | 0.057 | -0.114 | 0.228 | -1812.4 | yes |
| 0.3 | 0.129 | -0.004 | 0.262 | -1802.9 | yes |
| 0.6 | 0.225 | 0.119 | 0.332 | -1797.6 | yes |
| 0.9 | 0.279 | 0.161 | 0.396 | -1795.5 | yes |
| 1.2 | 0.368 | 0.234 | 0.502 | -1795.9 | yes |
sweep_df <- as.data.frame(sweep)
ggplot2::ggplot(sweep_df, ggplot2::aes(beta, estimate)) +
ggplot2::geom_ribbon(
ggplot2::aes(ymin = conf.low, ymax = conf.high),
fill = "#56B4E9", alpha = 0.3
) +
ggplot2::geom_line(colour = "#0072B2", linewidth = 0.9) +
ggplot2::geom_point(colour = "#0072B2", size = 2) +
ggplot2::geom_hline(yintercept = truth, linetype = "dashed",
colour = "#000000") +
ggplot2::geom_vline(xintercept = beta_true, linetype = "dotted",
colour = "#D55E00") +
ggplot2::annotate("text", x = min(sweep_df$beta), y = truth,
label = "complete-data mean", hjust = 0, vjust = -0.6,
size = 3.2) +
ggplot2::annotate("text", x = beta_true, y = min(sweep_df$conf.low),
label = "generating slope", hjust = 1.05, vjust = 0,
size = 3.2, colour = "#D55E00") +
ggplot2::labs(
x = "assumed slope",
y = "pooled mean of y",
title = "The pooled mean of y under each assumed slope"
) +
ggplot2::theme_minimal(base_size = 11)
The pooled mean of y and its 95% confidence interval at each assumed slope. The dashed line is the mean of the complete data. The dotted line marks the slope that generated the data.
The row at slope 0 makes the same missing-at-random assumption as the
earlier table, but proxy_mnar_sensitivity() refits the
mixture and draws its own completions. Its mean of 0.057 therefore
differs from the earlier 0.061 by 0.004.
Censoring at a known limit
When values are missing because they fell beyond a known limit, the
mechanism is known and no sweep is needed. Here every value of
y below 0.3 is hidden, as if 0.3 were the detection limit
of an assay. censored("y", upper = 0.3) draws each hidden
value from the mixture’s conditional distribution cut off at 0.3. A
common alternative replaces each hidden value with half the detection
limit.
thr <- 0.3
cmiss <- z_full[, 2L] < thr
cdat <- z_full
cdat[cmiss, "y"] <- NA
cfit <- gmm_impute(cdat, N = 2L, m = m_draws,
mechanism = censored("y", upper = thr), seed = 1L)
cens_est <- proxy_pool(cfit, "y", method = "rubin")$estimate
half_est <- mean(ifelse(cmiss, thr / 2, z_full[, 2L]))
below_half <- mean(z_full[cmiss, 2L] < thr / 2)| Data used | Mean of y | Difference from complete data |
|---|---|---|
| complete data, before censoring | 0.285 | 0.000 |
| rows not censored | 1.102 | 0.817 |
| hidden values set to half the limit | 0.631 | 0.345 |
| imputed, censored below the limit | 0.347 | 0.061 |
Comparison with mice, Amelia, sampleSelection and the Tobit model
The examples above use one dataset each. To check proxymix against established tools, a simulation repeated both analyses on 300 datasets of 300 rows each.
The first design is the missing-not-at-random example above, at half
the size. The target is the mean of y. proxymix swept the
slopes 0, 0.35, 0.7 and 1.05. mice (van Buuren and
Groothuis-Oudshoorn, 2011) drew each missing y from a
normal regression on x1 and then added a fixed shift to
every imputed value. This is called delta adjustment, and a shift of 0
is missing at random. mice used the shifts 0, 0.2, 0.4 and
0.6, in units of y, in place of proxymix’s slopes.
Amelia (Honaker, King and Blackwell, 2011) assumes missing
at random. The Heckman two-step estimator (Heckman, 1979) in
sampleSelection (Toomet and Henningsen, 2008) models
whether y is observed and the value of y
together. It was run with x1 in both parts and without an
exclusion restriction, that is, without a variable that affects whether
y is observed but not y itself. This design
has no such variable, and without one the Heckman estimator is known to
be unreliable.
In the second design, y is recorded as zero whenever it
falls below zero, which happens to 30% of the values. The target is the
slope of y on a predictor x. Its true value is
1. The Tobit model (Tobin, 1958) treats each zero as a value known only
to lie at or below zero. tobit() in AER
(Kleiber and Zeileis, 2008) and survreg() in
survival (Therneau and Grambsch, 2000) fit this model and
gave identical results. With proxymix, the zeros were set to missing and
imputed with censored("y", upper = 0), and the regression
was pooled by Rubin’s rules.
Each imputer made 10 completed datasets. Every interval is a 95% interval. Coverage is the share of datasets whose interval contained the true value, and it should be close to 0.95.
| Target and method | Bias | Error | Coverage | Interval width |
|---|---|---|---|---|
| Mean of y: complete data, before deletion | 0.003 | 0.059 | 0.947 | 0.234 |
| Mean of y: rows not deleted | -0.284 | 0.294 | 0.047 | 0.295 |
| Mean of y: Amelia (missing at random) | -0.187 | 0.200 | 0.250 | 0.284 |
| Mean of y: proxymix, slope 0.7 | -0.002 | 0.077 | 0.947 | 0.305 |
| Mean of y: mice, shift 0.4 | -0.014 | 0.074 | 0.950 | 0.297 |
| Mean of y: Heckman two-step, no exclusion restriction | 0.341 | 116.203 | 0.986 | 2382.095 |
| Slope: complete data, before censoring | -0.005 | 0.062 | 0.933 | 0.229 |
| Slope: zeros taken at face value | -0.302 | 0.306 | 0.000 | 0.183 |
| Slope: Tobit model (AER, survival) | -0.004 | 0.073 | 0.937 | 0.269 |
| Slope: proxymix, censored imputation | -0.033 | 0.085 | 0.907 | 0.317 |
Assuming missing at random, proxymix at slope 0, mice at
shift 0 and Amelia all understated the mean by about 0.19,
with intervals that contained it in only 0.23 to 0.25 of datasets. At
the grid values closest to the truth, proxymix and mice did
equally well, with coverages of 0.947 and 0.950, although with real data
neither the right slope nor the right shift is known. Without an
exclusion restriction, the Heckman estimator broke down, with an error
of 116 and no valid interval in 13 of 300 datasets. In the censored
design the Tobit model did better than proxymix, with an error of 0.073
against 0.085 and a coverage of 0.937 against 0.907, even though the
proxymix intervals were wider.
proxymix is the slowest method here. On one dataset, its sweep over
four slopes took 4.3 seconds against 0.29 for the four mice
shifts, and its censored imputation took 0.46 seconds against 0.002 for
the Tobit fit (median of five runs on one computer).
The code below analyses one dataset from each design with every method. It repeats the simulation code for a single dataset.
library(proxymix)
library(mice)
library(Amelia)
library(AER)
library(survival)
library(sampleSelection)
# one dataset of 300 rows in which larger values of y are more often deleted
set.seed(1L)
n <- 300L
comp <- sample(1:2, n, replace = TRUE)
mu <- rbind(c(0, 0), c(1.5, 0.5))
chol_r <- chol(matrix(c(1, 0.6, 0.6, 1), 2L))
z <- matrix(rnorm(2 * n), n, 2L) %*% chol_r + mu[comp, ]
full <- data.frame(x1 = z[, 1L], y = z[, 2L])
obs <- full
obs$y[runif(n) < plogis(-0.5 + 0.7 * full$y)] <- NA
# proxymix: the pooled mean of y at four assumed slopes
proxy_mnar_sensitivity(obs, "y", beta_grid = c(0, 0.35, 0.7, 1.05),
m = 10L, seed = 1L)
# mice: add a fixed shift to every imputed y, then pool the mean
lapply(c(0, 0.2, 0.4, 0.6), function(delta) {
post <- make.post(obs)
post["y"] <- paste0("imp[[j]][, i] <- imp[[j]][, i] + ", delta)
imp <- mice(obs, m = 10L, method = "norm", post = post, seed = 1L,
printFlag = FALSE)
summary(pool(with(imp, lm(y ~ 1))), conf.int = TRUE)
})
# Amelia: ten completed datasets, pooled by mice::pool()
fits <- lapply(amelia(obs, m = 10L, p2s = 0L)$imputations,
function(d) lm(y ~ 1, data = d))
summary(pool(fits), conf.int = TRUE)
# sampleSelection: Heckman two-step with x1 in both parts; the mean of y
# is the fitted model for y at the mean of x1, with an approximate interval
obs$seen <- !is.na(obs$y)
fit <- heckit(seen ~ x1, y ~ x1, data = obs)
x_bar <- c(1, mean(obs$x1))
est <- sum(x_bar * coef(fit)[3:4])
se <- sqrt(as.numeric(t(x_bar) %*% vcov(fit)[3:4, 3:4] %*% x_bar))
c(estimate = est, conf.low = est - qnorm(0.975) * se,
conf.high = est + qnorm(0.975) * se)
# one dataset of 300 rows with y recorded as 0 whenever it falls below 0
set.seed(1L)
a_cens <- -sqrt(2) * qnorm(0.3)
x <- rnorm(n)
y_star <- a_cens + x + rnorm(n)
cens <- data.frame(x = x, y = pmax(y_star, 0))
# proxymix: set the zeros to missing, draw them below 0, pool the regression
holes <- cens
holes$y[cens$y == 0] <- NA
imp <- gmm_impute(holes, m = 10L, mechanism = censored("y", upper = 0),
seed = 1L)
summary(pool(lapply(complete(as_mids(imp), "all"),
function(d) lm(y ~ x, data = d))), conf.int = TRUE)
# the Tobit model, fitted by AER and by survival
coef(tobit(y ~ x, data = cens))
coef(survreg(Surv(y, y > 0, type = "left") ~ x, data = cens,
dist = "gaussian"))The code needs mice, Amelia,
AER, survival and
sampleSelection, all on CRAN, and it is not run when this
vignette is built. The extended
version of this article gives the full simulation and applies both
mechanisms to two real datasets.