How can I simulate NMAR missing data with bias?

Question

Question summary

I am trying to simulate several random variables (using R), then introduce missingness into one of them, and show that a complete-case regression model yields biased results when the missingness is not at random (NMAR), as opposed to MCAR (missing completely at random). Several textbooks argue that NMAR leads to bias, but I am failing to see this in my simulations. What am I missing?

Step 1: Simulating random variables

I simulate two independent variables (x1 and x2), a dependent variable (y) and an unmodeled variable z by sampling 1,000 draws from a multivariate normal distribution. The correlation structure is such that x1 and x2 are moderately correlated with each other (at 0.3), highly with y (at 0.7 and 0.6), and not at all with z (0.0), which is in turn correlated with y at 0.5.

library("MASS")
VCOV <- matrix(c(1.0, 0.3, 0.7, 0.0,
                 0.3, 1.0, 0.6, 0.0,
                 0.7, 0.6, 1.0, 0.5,
                 0.0, 0.0, 0.5, 1.0), nrow = 4, byrow = TRUE)
mat <- mvrnorm(n = 10000, mu = c(0, 0, 0, 0), Sigma = VCOV)
colnames(mat) <- c("x1", "x2", "y", "z")
head(mat)

This yields:

             x1          x2          y           z
[1,] -0.2296673 -0.77817017 -0.5212378 -0.13099032
[2,] -0.5500676 -1.05473085 -0.4393144  0.07021164
[3,]  0.5009282  0.21303023  0.2580359 -1.30247661
[4,]  0.1491819 -0.36735643  0.6037848  1.29704084
[5,]  0.1457850  0.06076667 -1.1044882 -0.94487507
[6,]  1.1976893  1.46720672  1.9255790  0.37833377

As expected, cor(mat) returns something like this:

            x1           x2         y            z
x1 1.000000000  0.301028137 0.7093628  0.006123339
x2 0.301028137  1.000000000 0.6011744 -0.004656056
y  0.709362832  0.601174422 1.0000000  0.491996893
z  0.006123339 -0.004656056 0.4919969  1.000000000

I estimate two linear models without intercept, with and without z, to show that z is not a confounder and can be left out:

library("texreg")
model1 <- lm(y ~ x1 + x2 - 1, data = dat)
model2 <- lm(y ~ x1 + x2 + z - 1, data = dat)
screenreg(list(model1, model2), single.row = TRUE)

The result:

=================================================
           Model 1             Model 2           
-------------------------------------------------
x1            0.57 (0.02) ***     0.58 (0.01) ***
x2            0.38 (0.02) ***     0.41 (0.01) ***
z                                 0.50 (0.01) ***
-------------------------------------------------
R^2           0.65                0.89           
Adj. R^2      0.65                0.89           
Num. obs.  1000                1000              
=================================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Step 2: Introducing NMAR missing data

NMAR missing data are characterized by observations on the dependent variable that are missing as a function of another variable that is unobserved, according to the textbooks I have read. z is supposed to be this unobserved variable. It is correlated with y but not with x1 and x2 (otherwise the model would have omitted variable bias and would be pointless to begin with), so it can be used to create a non-random missingness pattern.

s <- dat$y <= median(dat$y)
dat$y_nmar <- dat$y
dat$y_nmar[s] <- NA
model3 <- lm(y_nmar ~ x1 + x2 - 1, data = dat)

Alternatively, I can use fitted values of regressing y on z to determine which values should be missing:

s <- predict(lm(dat$y ~ dat$z))
s <- s <= median(s)
dat$y_nmar <- dat$y
dat$y_nmar[s] <- NA
model4 <- lm(y_nmar ~ x1 + x2 - 1, data = dat)
screenreg(list(model1, `Model 3` = model3, `Model 4` = model4))

But the results do not differ much although the missing data are missing not at random:

==============================================
           Model 1      Model 3     Model 4   
----------------------------------------------
x1            0.57 ***    0.54 ***    0.57 ***
             (0.02)      (0.03)      (0.03)   
x2            0.38 ***    0.39 ***    0.39 ***
             (0.02)      (0.03)      (0.03)   
----------------------------------------------
R^2           0.65        0.63        0.65    
Adj. R^2      0.65        0.63        0.65    
Num. obs.  1000         500         500       
==============================================
*** p < 0.001; ** p < 0.01; * p < 0.05

One could of course argue that this is because there is no correlation between the missingness process and the observed covariates. But then one could also argue that a reasonable scientist should not be estimating a model with significant omitted variable bias to begin with. So either there is no bias due to NMAR missingness, or the missingness does not matter in practice because the model has grave problems anyway. What am I missing here? Is there a coding problem? Or a problem with my argument? Or is the textbook definition I found incomplete/imprecise?

Answer 1

You've correctly found that the data being NMAR is not a sufficient condition for there to be bias in the estimated parameters. To make the example more obvious, consider a simple linear regression with only 1 covariate and no correlation between $x$ and $y$.

set.seed(1)
n <- 10000
x <- rnorm(n)
y <- rnorm(n)
mod1 <- lm(y ~ 0 + x)
summary(mod1)
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# x 0.004768   0.009787   0.487    0.626

As we expect, the estimated slope parameter is 0. Now let's try making the data NMAR in some extreme way. We can remove all data with $y$ values above 0. This is NMAR because the missingness depends directly on the data that is missing.

y_nmar <- y
y_nmar[y > 0] <- NA
mod2 <- lm(y_nmar ~ 0 + x)
summary(mod2)
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# x 0.008405   0.013559    0.62    0.535

There is still no relation and the estimate is unbiased, despite being NMAR.

But, if we also remove all data that has $x$ greater than 0,

x_nmar <- x
x_nmar[x  > 0] <- NA
mod3 <- lm(y_nmar ~ 0 + x_nmar)
summary(mod3)
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# x_nmar  0.60184    0.01433   42.01   <2e-16 ***

Then the estimate is extremely biased. Note that in both my example and yours, if you include an intercept in the models then the estimated intercept shows bias immediately upon making the data NMAR. So whether bias exists depends on both the missingness mechanism and the particular model that is being fit. But NMAR missingness can certainly cause bias.