---
title: "Missing data that depends on the missing value"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Missing data that depends on the missing value}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- Render time: ~7 s under rmarkdown::render() with ggplot2 installed;
     macOS arm64 (Apple silicon), R 4.6.1, one core. The sweep over five
     slopes takes most of it. The comparison with mice, Amelia,
     sampleSelection, AER and survival is read from
     results/missing_data_mnar.rds, built by
     data-raw/vignette_results/missing_data_mnar.R. -->

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  dpi = 150,
  out.width = "100%"
)
```

```{r library}
library(proxymix)
```

```{r engines}
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
```

```{r 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_mnar.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_mnar.rds was built under proxymix ",
       res$proxymix_version, ", but this is proxymix ",
       packageVersion("proxymix"), ". Rerun the simulation and ",
       "data-raw/vignette_results/missing_data_mnar.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)
}
```

## The problem

*Imputing missing data with a mixture* shows how to fill the holes in a
dataset under one assumption, called missing at random. The chance that a value is
missing may depend on the other values in its row, but not on the missing
value itself.

Often the missing values are mostly the large ones, or mostly the small
ones. An assay saturates above some level. A yield is never recorded on the
paddocks that did badly. The patients who are most unwell miss their
follow-up visit. The values that remain are then unrepresentative, and
their mean is biased. Imputation under missing at random removes only part
of this bias, because its model is fitted to the same unrepresentative
values.

Two such cases are common. The first is censoring: a value is missing
because it fell below (or above) a known limit, such as the detection limit
of an assay. The missing value is then known to lie on the far side of the
limit. The second is missing not at random: the chance that a value is
missing rises or falls with the value itself, by an unknown amount. The
observed data cannot settle how strong the link is. The usual remedy
is a sensitivity analysis (Little, 1993): the analysis is repeated under
several assumed strengths, and all the results are reported.

This vignette works through both cases on simulated data. Because the
deleted values are kept, every estimate can be checked against the truth.
The vignette then compares proxymix with established packages for each
case.

## Package capabilities

- `gmm_impute()` fits a Gaussian mixture, a sum of a few normal
  distributions, to data with holes. It returns `m` completed datasets, as
  in *Imputing missing data with a mixture*. Its `mechanism` argument
  specifies how the values came to be missing.
- `mar()` is missing at random, the default.
- `censored()` is for a missing value known to lie in a range. For
  example, `censored("y", upper = 0.3)` specifies that each missing `y`
  lies below 0.3.
- `mnar()` is for a value of `y` whose chance of being missing depends on
  `y` itself. You supply the strength of the link.
- `proxy_mnar_sensitivity()` repeats the imputation for a range of
  strengths and returns the pooled mean of `y` for each.
- `proxy_pool()` pools the mean of a column over the completed datasets.

Each missing value is drawn from its conditional distribution: the
distribution of the missing entry given the observed entries in the same
row. Under `censored()`, this distribution is cut off at the limit, and
every imputed value falls on the censored side of it. The cut-off distribution has an exact
formula. Under `mnar()`, the conditional distribution is weighted by the
chance of being missing at each value. Values that were more likely to go
missing are then drawn more often.

`mnar()` models the chance that `y` is missing as `plogis(alpha + beta * y)`,
the logistic curve used in logistic regression. This is the selection model
of Diggle and Kenward (1994). The slope `beta` is the strength of the link,
and `beta = 0` is missing at random. You supply `beta`. The package sets the
intercept `alpha` so that the model gives the observed share of missing
values.

```{r mechanism-table, echo = FALSE}
knitr::kable(
  data.frame(
    mechanism = c("missing at random", "censored", "missing not at random"),
    known = c("nothing beyond the rest of its row",
              "it lies beyond a known limit",
              "its chance of going missing depended on its size"),
    call = c("`mar()`", "`censored()`", "`mnar()`"),
    settled = c("yes, once the mixture is assumed",
                "yes, the limit is known",
                "only through the assumed shape of `y`"),
    stringsAsFactors = FALSE
  ),
  col.names = c("Mechanism", "What is known about a missing value",
                "Function", "Can the observed data fix the imputation model?"),
  caption = "The three mechanisms that `gmm_impute()` accepts."
)
```

## 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.

```{r dgp}
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"] <- NA
```

### Impute 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.

```{r recover}
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")$estimate
```

For 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.

```{r recover-table, echo = FALSE}
rec_tbl <- data.frame(
  data_used = c("complete data, before deletion", "rows not deleted",
                "imputed, missing at random",
                "imputed, missing not at random (slope 0.7)"),
  estimate = c(truth, mean(dat[!miss, "y"]), mar_est, mnar_est),
  stringsAsFactors = FALSE
)
rec_tbl$diff <- rec_tbl$estimate - truth
knitr::kable(
  rec_tbl, digits = 3L,
  col.names = c("Data used", "Mean of y", "Difference from complete data"),
  caption = paste0(
    "The mean of y from the complete data, from the rows not deleted, and ",
    "pooled over each imputation. ", sum(miss), " of ", n, " values of y ",
    "were deleted."
  )
)
```

### 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.

```{r sweep}
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]
```

```{r sweep-table, echo = FALSE}
sweep_tbl <- as.data.frame(sweep)[, c("beta", "estimate", "conf.low",
                                      "conf.high", "loglik", "converged")]
sweep_tbl$converged <- ifelse(sweep_tbl$converged, "yes", "no")
knitr::kable(
  sweep_tbl,
  digits = c(1L, 3L, 3L, 3L, 1L, 0L),
  col.names = c("Assumed slope", "Pooled mean of y", "CI lower",
                "CI upper", "Log-likelihood", "Converged"),
  caption = paste0(
    "The pooled mean of y at each assumed slope. The mean of y in the ",
    "complete data is ", round(truth, 3), "."
  )
)
```

```{r fig-sweep, eval = has_ggplot2, echo = has_ggplot2, fig.cap = "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.", fig.alt = "Pooled mean of y against the assumed slope, rising from left to right, with a shaded confidence band, a point at each assumed slope, a dashed horizontal line at the complete-data mean, and a dotted vertical line at the generating slope of 0.7."}
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)
```

```{r fig-sweep-skip, eval = !has_ggplot2, echo = FALSE, results = "asis"}
cat("ggplot2 is not installed on this build, so the sensitivity figure",
    "is skipped. The table above gives the same values.\n")
```

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 `r fixed(sweep$estimate[sweep$beta == 0], 3)`
therefore differs from the earlier `r fixed(mar_est, 3)` by
`r fixed(abs(sweep$estimate[sweep$beta == 0] - mar_est), 3)`.

### 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.

```{r censor}
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)
```

```{r censor-table, echo = FALSE}
cens_tbl <- data.frame(
  data_used = c("complete data, before censoring", "rows not censored",
                "hidden values set to half the limit",
                "imputed, censored below the limit"),
  estimate = c(truth, mean(cdat[!cmiss, "y"]), half_est, cens_est),
  stringsAsFactors = FALSE
)
cens_tbl$diff <- cens_tbl$estimate - truth
knitr::kable(
  cens_tbl, digits = 3L,
  col.names = c("Data used", "Mean of y", "Difference from complete data"),
  caption = paste0(
    "The mean of y when values below ", thr, " are hidden. ", sum(cmiss),
    " of ", n, " values were hidden."
  )
)
```

### Comparison with mice, Amelia, sampleSelection and the Tobit model

```{r 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]
}
mv <- function(method, what) sim_value("mnar", "mean", method, what)
cv <- function(method, what) sim_value("censored", "slope", method, what)
mar_methods <- c("proxymix, slope 0", "mice, shift 0", "Amelia")
mar_bias <- vapply(mar_methods, mv, numeric(1L), what = "bias")
mar_cov <- vapply(mar_methods, mv, numeric(1L), what = "coverage")
and_list <- function(v) {
  paste(paste(v[-length(v)], collapse = ", "), "and", v[length(v)])
}
```

The examples above use one dataset each. To check proxymix against
established tools, a simulation repeated both analyses on `r res$n_rep`
datasets of `r res$n` 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
`r and_list(res$beta_grid)`. `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 `r and_list(res$delta_grid)`, 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 `r res$m` 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.

```{r compare-table, echo = FALSE}
cmp_rows <- data.frame(
  design = c(rep("mnar", 6L), rep("censored", 4L)),
  method = c("complete data", "available cases", "Amelia",
             "proxymix, slope 0.7", "mice, shift 0.4",
             "sampleSelection heckit",
             "complete data", "zeros at face value", "AER tobit",
             "proxymix"),
  label = c("Mean of y: complete data, before deletion",
            "Mean of y: rows not deleted",
            "Mean of y: Amelia (missing at random)",
            "Mean of y: proxymix, slope 0.7",
            "Mean of y: mice, shift 0.4",
            "Mean of y: Heckman two-step, no exclusion restriction",
            "Slope: complete data, before censoring",
            "Slope: zeros taken at face value",
            "Slope: Tobit model (AER, survival)",
            "Slope: proxymix, censored imputation"),
  stringsAsFactors = FALSE
)
cmp_rows$estimand <- ifelse(cmp_rows$design == "mnar", "mean", "slope")
cmp_tbl <- data.frame(
  label = cmp_rows$label,
  bias = mapply(sim_value, cmp_rows$design, cmp_rows$estimand,
                cmp_rows$method, "bias"),
  rmse = mapply(sim_value, cmp_rows$design, cmp_rows$estimand,
                cmp_rows$method, "rmse"),
  coverage = mapply(sim_value, cmp_rows$design, cmp_rows$estimand,
                    cmp_rows$method, "coverage"),
  width = mapply(sim_value, cmp_rows$design, cmp_rows$estimand,
                 cmp_rows$method, "width"),
  stringsAsFactors = FALSE
)
knitr::kable(
  cmp_tbl, digits = 3L, row.names = FALSE,
  align = c("l", "r", "r", "r", "r"),
  col.names = c("Target and method", "Bias", "Error", "Coverage",
                "Interval width"),
  caption = paste0(
    "Results over ", res$n_rep, " simulated datasets per design. Bias is ",
    "the average difference from the true value, and error is the root ",
    "mean squared error. The mean of y is ", res$truth[["mnar_mean"]],
    " in the population, and the true slope is ",
    res$truth[["cens_slope"]], ". The proxymix and mice rows are the grid ",
    "values closest to the truth. 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), 3L), ". The Heckman coverage is ",
    "over the ", res$n_rep - res$n_undefined, " datasets in which its ",
    "estimated variance was positive."
  )
)
```

Assuming missing at random, proxymix at slope 0, `mice` at shift 0 and
`Amelia` all understated the mean by about `r fixed(-mean(mar_bias), 2)`,
with intervals that contained it in only `r fixed(min(mar_cov), 2)` to
`r fixed(max(mar_cov), 2)` of datasets. At the grid values closest to the
truth, proxymix and `mice` did equally well, with coverages of
`r fixed(mv("proxymix, slope 0.7", "coverage"), 3)` and
`r fixed(mv("mice, shift 0.4", "coverage"), 3)`, 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
`r round(mv("sampleSelection heckit", "rmse"))` and no valid interval in
`r res$n_undefined` of `r res$n_rep` datasets. In the censored design the
Tobit model did better than proxymix, with an error of
`r fixed(cv("AER tobit", "rmse"), 3)` against
`r fixed(cv("proxymix", "rmse"), 3)` and a coverage of
`r fixed(cv("AER tobit", "coverage"), 3)` against
`r fixed(cv("proxymix", "coverage"), 3)`, even though the proxymix
intervals were wider.

proxymix is the slowest method here. On one dataset, its sweep over four
slopes took `r fixed(res$time_secs[["proxymix sweep"]], 1)` seconds
against `r fixed(res$time_secs[["mice sweep"]], 2)` for the four `mice`
shifts, and its censored imputation took
`r fixed(res$time_secs[["proxymix censored"]], 2)` seconds against
`r fixed(res$time_secs[["tobit"]], 3)` 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.

```{r compare-code, eval = FALSE}
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](https://max578.github.io/proxymix/articles/extended/missing_data_mnar.html)
gives the full simulation and applies both mechanisms to two real
datasets.

## Interpretation

Of the `r n` values of `y`, `r sum(miss)` were deleted. Because the
larger values were deleted more often, the rows not deleted give a mean of
`r fixed(mean(dat[!miss, "y"]), 3)`, against
`r fixed(truth, 3)` in the complete data. Imputing under missing at random
gives `r fixed(mar_est, 3)`, which is still
`r fixed(truth - mar_est, 3)` too low. That imputation model is fitted to
the rows that remain, which have too few large values. Supplying the true
slope to `mnar()` gives `r fixed(mnar_est, 3)`, within
`r fixed(abs(mnar_est - truth), 3)` of the complete-data mean.

In the sweep, the pooled mean
`r if (all(diff(sweep$estimate) > 0)) "rises at every step" else "changes"`
from `r fixed(sweep$estimate[1L], 3)` at slope 0 to
`r fixed(sweep$estimate[nrow(sweep)], 3)` at slope `r max(sweep$beta)`.
Its interval first contains the complete-data mean at a slope of
`r beta_first`. The data were generated with a slope of `r beta_true`.
The fit
`r if (all(sweep$converged)) "converged at every slope" else paste0("did not converge at slope ", paste(sweep$beta[!sweep$converged], collapse = ", "))`.
The log-likelihood is highest at slope `r beta_best`, where it is
about `r round(ll_gain)` above its value at slope 0. This difference
reflects the assumed two-component shape of `y`, not why the values went
missing. A report should show the whole curve and leave the choice of a
plausible slope to the reader.

In the censoring example, `r sum(cmiss)` values fell below the limit of
`r thr`.
The rows not censored give a mean of `r fixed(mean(cdat[!cmiss, "y"]), 3)`,
which is `r fixed(mean(cdat[!cmiss, "y"]) - truth, 3)` too high. Setting
each hidden value to half the limit gives `r fixed(half_est, 3)`, which is
still `r fixed(half_est - truth, 3)` too high. Half the limit is a
convention for concentrations, which cannot fall below zero. Here `y` can
be negative, and `r round(100 * below_half)` per cent of the hidden values
lie below `r thr / 2`. The censored imputation gives
`r fixed(cens_est, 3)`, within `r fixed(abs(cens_est - truth), 3)` of the
complete-data mean.

## Limitations

`censored()` and `mnar()` act on one column. A row missing that column is
assumed to have all its other columns observed. This suits a detection
limit or a single outcome, not holes spread across many columns. Estimates
other than a column mean are pooled by passing the completed datasets to
`mice` through `as_mids()`. They then rely on the large-sample
assumptions of Rubin's rules.

The sweep checks one kind of departure from missing at random. The chance
of being missing is a logistic curve in the missing value alone, with the
intercept set from the observed share of missing values. A mechanism that
also depends on a variable not in the data is outside what the sweep
covers. The sweep is not a test, and its log-likelihood cannot tell you
which slope is right. The range of results is only as wide as the grid you
choose.

The single-dataset examples use one simulated dataset of `r n` rows with
`r m_draws` completed datasets. Their numbers would change with another
seed. The settings that matter most are the number of components and the
assumed slope. A mixture with too few components cannot represent the two
groups in the data. The pooled mean moves steadily as the assumed slope
moves away from the truth, as the sweep table shows.

With data censored at a known limit, the Tobit model was more accurate than
proxymix in the simulation. The proxymix slope was biased by
`r fixed(cv("proxymix", "bias"), 3)`, and its 95% intervals contained the
true slope in only `r fixed(cv("proxymix", "coverage"), 3)` of datasets.
When the outcome follows a normal linear regression below the limit, as
in this design, the Tobit model is the better choice. The simulation covers one mixture shape, one sample size, a
logistic mechanism on one column, censoring at one known limit, and grids
that contain a value near the truth. It does not show how either sweep
behaves when its grid does not reach the truth.

## Further reading

*Imputing missing data with a mixture* covers data missing at random, the
case this vignette sets aside. It shows how a mixture and a single normal
distribution differ in the values they impute.

*The closed-form operator calculus on a mixture* explains the formulas
behind the conditional distributions used here.

## References

Diggle, P. and Kenward, M. G. (1994). *Informative drop-out in
longitudinal data analysis.* Journal of the Royal Statistical Society C
43(1), 49--93. <https://doi.org/10.2307/2986113>.

Heckman, J. J. (1979). *Sample selection bias as a specification error.*
Econometrica 47(1), 153--161. <https://doi.org/10.2307/1912352>.

Honaker, J., King, G. and Blackwell, M. (2011). *Amelia II: A program for
missing data.* Journal of Statistical Software 45(7), 1--47.
<https://doi.org/10.18637/jss.v045.i07>.

Hoek, J. van der and Elliott, R. J. (2024). *Mixtures of multivariate
Gaussians.* Stochastic Analysis and Applications.
<https://doi.org/10.1080/07362994.2024.2372605>.

Kleiber, C. and Zeileis, A. (2008). *Applied Econometrics with R.*
Springer. <https://doi.org/10.1007/978-0-387-77318-6>.

Little, R. J. A. (1993). *Pattern-mixture models for multivariate
incomplete data.* Journal of the American Statistical Association
88(421), 125--134. <https://doi.org/10.1080/01621459.1993.10594302>.

Rubin, D. B. (1987). *Multiple Imputation for Nonresponse in Surveys.*
Wiley.

Therneau, T. M. and Grambsch, P. M. (2000). *Modeling Survival Data:
Extending the Cox Model.* Springer.
<https://doi.org/10.1007/978-1-4757-3294-8>.

Tobin, J. (1958). *Estimation of relationships for limited dependent
variables.* Econometrica 26(1), 24--36.
<https://doi.org/10.2307/1907382>.

Toomet, O. and Henningsen, A. (2008). *Sample selection models in R:
Package sampleSelection.* Journal of Statistical Software 27(7), 1--23.
<https://doi.org/10.18637/jss.v027.i07>.

van Buuren, S. and Groothuis-Oudshoorn, K. (2011). *mice: Multivariate
imputation by chained equations in R.* Journal of Statistical Software
45(3), 1--67. <https://doi.org/10.18637/jss.v045.i03>.

## Reproduce

The data are generated with seed `20260622`. Every call to `gmm_impute()`
and `proxy_mnar_sensitivity()` is given `seed = 1L`, so the completed
datasets are reproducible, and the seed does not change the random-number
state outside the call. The comparison is read from stored results of a
simulation run under proxymix `r res$proxymix_version`, `mice`
`r res$versions_desc[["mice"]]`, `Amelia` `r res$versions_desc[["Amelia"]]`,
`AER` `r res$versions_desc[["AER"]]`, `survival`
`r res$versions_desc[["survival"]]` and `sampleSelection`
`r res$versions_desc[["sampleSelection"]]`. It took about
`r round(res$elapsed_secs / 60)` minutes on one core and raised
`r if (res$n_warnings == 0L) "no warnings" else paste(res$n_warnings, "warnings")`.

```{r session-info, collapse = FALSE, class.output = "session-info"}
sessionInfo()
```