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. 
 A: 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


A: 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.
