Simulate MAR (Missing at Random) data

I am trying to generate MCAR, MAR and MNAR data. MCAR and MNAR are relatively easy. However I am struggling with MAR data.

I generate 500 observations with 2 variables (Y and X) out of a multivariate normal distribution with correlation 0.7.

library(mvtnorm)
set.seed(1994)
nobs = 500
nvar = 2
corr = 0.7
miss.prop = 0.4
mu <- rep(0,nvar)
Sigma <- matrix(corr, nrow=nvar, ncol=nvar) + diag(nvar)*(1-corr)  #draw from a multivariate normal
y <- rmvnorm(n=nobs, mean=mu, sigma=Sigma)


Then I use the package CoImp such that

 y.miss <- MAR(y, perc.miss = 0.4, setseed = 1994)@db.missing


Then I make a plot of the whole data. Black dots are the complete observations and the orange-ish triangles are the the missing data (so in theory we wouldn't observe them, but because it is simulated data we can see them :P ).

However, I am not sure how this algorithm is working because I do not really see any pattern on the data. As long as I understand missing values on Y should be dependent on the value of X..

I would really appreciate some clarification on this and maybe some short code on how to manually simulate MAR data.

• The plot you posted shows a clear pattern to me, namely a correlation and hence dependence between X and Y. That's exactly what You expected. So where is the problem? – matus Apr 18 '17 at 16:08
• @matus I am not looking for a pattern in the data itself as I already know that the correlation is 0.7 because I constructed myself as you can see in the code posted. I am looking for a pattern in the underlying missing data generating process, e.g. Do higher values of X lead to more NA's? – adrian1121 Apr 18 '17 at 16:10

With your data, MAR could be generated as follows:

# Create a normal data frame (not necessary, but makes the following easier)
data <- data.frame(y)
colnames(data) <- c("y", "x")

# Generate missing dummy y.miss
x_noise <- data$x + rnorm(length(data$x), 0, 0.5)
y.miss <- rep(0, length(data$x)) y.miss[x_noise < mean(data$x)] <- 1 # 1 corresponds to observed values in y; 0 corresponds to missing values in y

# Plot missings
plot(data$x[y.miss == 1], data$y[y.miss == 1], xlab = "x", ylab = "y", xlim = c(-4, 4), ylim = c(-4, 4))
points(data$x[y.miss == 0], data$y[y.miss == 0], col = 2)


UPDATE: Multivariate Case

# Correlated data
N <- 1000
y <- rnorm(N)
x1 <- y + rnorm(N, 7, 3)
x2 <- y + rnorm(N, 0, 5)
x3 <- y + rnorm(N, 10, 2)
x4 <- y + rnorm(N, -100, 1)

data <- data.frame(y, x1, x2, x3, x4)

# Calculate response propensity
mod <- - 1 * x1 - 0.08 * x2 + 1 * x3 - 0.0001 * x4 # Response model
rp <- exp(mod) / (exp(mod) + 1) # Suppress values between 0 and 1 via inverse-logit

# rp can be seen as probability to respond in y.
# See literature about reponse propensity for more details

# Create missings based on rp
y.miss <- rbinom(N, 1, rp)

# Plot missings y and x1
plot(data$x1[y.miss == 1], data$y[y.miss == 1], xlab = "x1", ylab = "y", xlim = c(-4, 20), ylim = c(-5, 5))
points(data$x1[y.miss == 0], data$y[y.miss == 0], col = 2)

# If you want to have missings with a specific response propensity, you could add a constant to your model
# e.g.
mod2 <- - 5 - 1 * x1 - 0.08 * x2 + 1 * x3 - 0.0001 * x4 # Response model
rp2 <- exp(mod2) / (exp(mod2) + 1) # Suppress values between 0 and 1 via inverse-logit

mean(rp) # Response rate without constant
mean(rp2) # Response rate with constant


• Thanks for your response, upvote! And how would you do the Multivariate case, for instance with 4 variables? Because i am interested in generating multivariate missing data as well. – adrian1121 Apr 20 '17 at 8:34
• I updated my answer with another example for multivariate data. In case of multivariate data I would calculate a response propensity first and insert the missings based on that. Just google for response propensity, you will find a lot of literature about the topic. – JSP Apr 20 '17 at 9:03
• Hi thank you very much, your response is really useful, although I may have not expressed correctly. I meant if it is possible to generate missing data not only Y but on X1 X2 X3 and X4. – adrian1121 Apr 20 '17 at 10:39
• Nice to hear that it helped! You can do the same procedure for every variable, for which you want to insert missings. If you calculate the response propensity several times (with different models if you want), you will get different response indicators for the different variables. – JSP Apr 20 '17 at 10:52

If you supply a single parameter value (0.4) to a function design to generate missingness-at-random, then this model is underspecified. According to the package documentation, the ... in MAR receives extra arguments to parametrize a copula used to generate a probability distribution for missingness over the space of observed values. You should consult the parameters argument and run examples from the fitCopula function in the copula package. The only advantage of using a copula is that it's easy to set the percentage missing values.

Generating MAR data is as easy as specifying a logistic model for missingness where other, non-missing features are supplied as covariates in that model. You can randomly generate missingness indicators from your RNG and recode those values to missing. To control the overall proportion of missing values, you can recalibrate the intercept of the logistic model using the uniroot command.