# How to generate MAR data with a fixed proportion of missing values?

To generate data with a missing at random (MAR) mechanism, usually we can first generate a complete data set and then model the missing probability for the variable $Y$; i.e. $Pr(Y=\text{missing}|{\bf X})$, using a logistical model, i.e. $$Pr(Y_i=1|{\bf X_i})=\frac{\exp(\beta'{\bf X_i})}{1+\exp(\beta'{\bf X_i})}.$$

Based on this probability we can generate the binary variable indicating if $Y$ is missing:

$$I_i = \text{rbinom}(1,p_i).$$

If $I_i=1$ then we delete the corresponding value of $Y$.

If we would like to control the proportion of missing values of $Y$, how can we do that? Since by using the above method the proportion of missing data depends on the values of the $X$'s then the generation of an indicator variable is random, so it's hard to control the overall proportion of missing values.

The expected number of missing values is the mean of the probabilities

$$f_X(\alpha, \beta) = \frac{1}{n}\sum_i \frac{1}{1 + \exp(-\alpha-X_i\beta^\prime)}$$

where I have written $\alpha$ for the constant term, $\beta$ as a row vector of $d$ coefficients, and followed the convention of arranging the variable values in columns of the design matrix $X$ with $d$ columns (say) and $n$ rows so that $X_i$ is the $i^\text{th}$ row, representing one case.

This is a smooth function of the $d+1$ parameters $\alpha, \beta$ whose range is the full interval $(0,1)$. Therefore, given a target proportion $p$ there exists a $d$-dimensional manifold $f_X^{-1}(p)$ of parameter values for which

$$f_X(\alpha, \beta) =p.$$

I interpret the question as a request for a way to find $\alpha, \beta$ that satisfy this equation.

One way to find a definite solution is to pick a one-dimensional path $t\to\gamma(t)$ in the $(\alpha, \beta)$ space and compute its intersection with $f^{-1}(p)$, which amounts to finding a univariate root. When the path is smooth and intersects $f^{-1}(p)$ transversally, a standard root-finding method such as Newton-Raphson will usually work very well to solve the one-variable equation

$$f(\gamma(t)) - p = 0.$$

For instance, among the infinitely many choices available, you might fix the coefficients $\beta$ and vary $\alpha$, giving the path $t\to (t, \beta)$, or you might fix all the coefficients initially and use the path $t\to (t\alpha, t\beta)$ to scale them uniformly.

As a worked example of the first method (varying only the intercept), consider the following R code. After creating some sample data x and choosing a fixed set of coefficients beta, it sets up a function f to compute $t\to f(\gamma(t))$ and then invokes a root-finding routine uniroot to find zeros of $f(\gamma(t))-p$. Any such zero will be a usable value of $\alpha$, thereby fully specifying the missingness model (which can later be used to generate simulated datasets).

n <- 10; d <- 3               # Set the dimensions of the data
set.seed(17)                  # Create a reproducible starting point
x <- matrix(rnorm(n*d), n, d) # Create independent variables
#
# Determine a model with a given expected proportion of responses
# missing at random.
#
beta <- 1:d                   # Fix the coefficients at any desired values
f <- function(t) {            # Define a path through parameter space
sapply(t, function(y) mean(1 / (1 + exp(-y -x %*% beta))))
# (sapply makes this function vectorizable)
}
#
# Find parameters (alpha, beta) yielding any specified proportions p.
#
p <- seq(.1, .9, .2)                        # Specify the target proportions
results <- sapply(p, function(p) {
alpha <- uniroot(function(t) f(t) - p, c(-1e6, 1e6), tol = .Machine$double.eps^0.5)$root
c(alpha, f(alpha))})
dimnames(results) <- list(c("alpha", "f(alpha)"), p=p)
print(results)


The output tabulates some solutions for various values of $p$ along with a check of their validity:

          p
0.1      0.3       0.5        0.7     0.9
alpha    -5.102146 -3.45776 -2.251828 -0.7976692 1.27534
f(alpha)  0.100000  0.30000  0.500000  0.7000000 0.90000


Modify p as desired. Change f and/or beta to search for other sets of parameters.

• So you meant that the expected proportion of missing values is the MEAN of the probabilities? $f_X(α,β)=\frac{1}{n}∑_i^n\frac{1}{1+exp(−α−X_iβ′)}$. Thank you for your answer, this makes sense. – askming Jul 29 '14 at 19:40
• Yes, thank you. Fortunately the code uses the mean as intended. I will change the opening description to match that. – whuber Jul 29 '14 at 19:53
• One more thing. In the R code where you define the function f, the sapply wrapper seems not necessary since later on when you use the f function, you only feed it with scalar rather than a vector. Correct me if I'm wrong. Thanks. – askming Jul 29 '14 at 20:15
• The vectorization is useful for graphing the function, for using it in an outer product, and all sorts of other things. – whuber Jul 29 '14 at 21:26