0
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ What are you doing with the missing values after you assign missingness? Are you simply removing those records before running the regression? Also a nitpick on your vocabulary, it is typically written MNAR, missing not at random as opposed to NMAR. $\endgroup$
    – astel
    Feb 5, 2021 at 20:06
  • $\begingroup$ The regression drops incomplete cases. See the Num. obs. row in the table. I have seen MNAR and NMAR in different textbooks. $\endgroup$ Feb 5, 2021 at 20:08

1 Answer 1

0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.