MCAR

Data which is Missing Completely At Random has nothing systematic about which observations are missing values. There is no relationship between missingness and either observed or unobserved covariates.

MAR

Missing At Random is weaker than MCAR. The missingness is still random, but due entirely to observed variables. For example, those from a lower socioeconomic status may be less willing to provide salary information (but we know their SES status). The key is that the missingness is not due to the values which are not observed. MCAR implies MAR but not vice-versa.

MNAR

If the data are Missing Not At Random, then the missingness depends on the values of the missing data. Censored data falls into this category. For example, individuals who are heavier are less likely to report their weight. Another example, the device measuring some response can only measure values above .5. Anything below that is missing.

Determining which type of missingness

Unfortunately, there is no statistical way to determine which category your data will fall. MCAR and MAR are assumptions that need to be made based upon your knowledge of the data and its collection mechanisms.

There is nothing in the data that can inform MNAR - If the missingness is due to the unobserved values, there is nothing in the data that can help us. MNAR should however be the default assumption.

To gain some insight, you can fit a logistic model predicting missingness based upon the covariates. If you find some relationships, that implies that MAR is more realistic than MCAR. However, it does not guarantee that we don’t have MNAR instead! Similarly, if we find no relationships, it may be because

• The missingness is MCAR.
• The missingness is MNAR.
• The missingness is really MAR, but you haven’t modeled it appropriately.

Methods to handle missingness

If the data is MCAR:

• Complete case analysis is valid.
• Multiple imputation or any other imputation method is valid.

If the data is MAR:

• Some complete cases analyses are valid under weaker assumptions than MCAR.
  • E.g. linear regression is unbiased if missingness is independent of the response, conditional on the predictors.
• Multiple imputation is valid (it is biased, but the bias is negligible).

If the data is MNAR:

• You must model the missingness explicitly; jointly modeling the response and missingness.
• In some specific cases (e.g. survival analysis), MNAR data (e.g. censored data) is handled appropriately.
• Generally, we assume MAR whenever possible just to avoid this situation.

mice

The mice function performs the imputation via chained equations. Chained equations are a variation of a Gibbs Sampler (an MCMC approach) that iterates between drawing estimates of missing values and estimates of parameters for distribution of the variable (both conditional on the other variables). Chained equations have fast convergence compared to most MCMC due to conditional independence of \(Y_t\) and \(Y_{t-1}\).

The default number of imputations, controlled by the m argument, is 5. This comes from the original work back in the 80’s and 90’s from a closed form estimate of the “efficiency” of the estimate - with 25% missing data and 5 imputations, efficiency is over 95%.

Update 2025: However, modern work has shown those results are only for the point estimates. If you similarly want your standard errors to be replicable, you need a much larger sample size. von Hippel has released an approach to estimate the required number of imputations; the number is often 100-200. There is an R package to carry out the calculations. For these notes I will stick with 5 imputations for ease, but in practice, in absence of other information or demands, starting with 100 imputations is more realistic.

The other major argument to be concerned with, method, controls how each imputed value is computed. While all use a chained equation approach, the details differ greatly. Without diving into the details too much, it suffices to say there are quite a few, which are documented well in the original mice paper. For our purposes, we’ll stick with the default, pmm (predictive mean matching), which has the benefits that it is semi-parametric in a sense (non-linear relationships are preserved even if the imputation model is mis-specified) and imputed values are restricted to observed values, so we cannot impute “impossible” values.

library(mice)
data(nhanes)
nhanes.imp <- mice(nhanes, m = 5, method = 'pmm')

## 
##  iter imp variable
##   1   1  bmi  hyp  chl
##   1   2  bmi  hyp  chl
##   1   3  bmi  hyp  chl
##   1   4  bmi  hyp  chl
##   1   5  bmi  hyp  chl
##   2   1  bmi  hyp  chl
##   2   2  bmi  hyp  chl
##   2   3  bmi  hyp  chl
##   2   4  bmi  hyp  chl
##   2   5  bmi  hyp  chl
##   3   1  bmi  hyp  chl
##   3   2  bmi  hyp  chl
##   3   3  bmi  hyp  chl
##   3   4  bmi  hyp  chl
##   3   5  bmi  hyp  chl
##   4   1  bmi  hyp  chl
##   4   2  bmi  hyp  chl
##   4   3  bmi  hyp  chl
##   4   4  bmi  hyp  chl
##   4   5  bmi  hyp  chl
##   5   1  bmi  hyp  chl
##   5   2  bmi  hyp  chl
##   5   3  bmi  hyp  chl
##   5   4  bmi  hyp  chl
##   5   5  bmi  hyp  chl

We can examine the imputations using the following:

# Show imputation 3
head(complete(nhanes.imp, 3))
##   age  bmi hyp chl
## 1   1 27.2   1 187
## 2   2 22.7   1 187
## 3   1 35.3   1 187
## 4   3 21.7   1 187
## 5   1 20.4   1 113
## 6   3 21.7   1 184

# Show just the imputed values. Columns are imputations, rows are observations
nhanes.imp$imp
## $age
## [1] 1 2 3 4 5
## <0 rows> (or 0-length row.names)
## 
## $bmi
##       1    2    3    4    5
## 1  27.2 27.5 27.2 25.5 27.4
## 3  33.2 27.4 35.3 27.2 27.2
## 4  27.4 27.4 21.7 24.9 22.5
## 6  24.9 25.5 21.7 20.4 27.5
## 10 25.5 21.7 22.5 20.4 27.4
## 11 33.2 26.3 29.6 20.4 26.3
## 12 27.5 28.7 20.4 25.5 20.4
## 16 30.1 22.5 29.6 28.7 28.7
## 21 28.7 30.1 27.4 30.1 30.1
## 
## $hyp
##    1 2 3 4 5
## 1  1 1 1 1 1
## 4  1 2 1 1 2
## 6  1 1 1 1 1
## 10 1 1 1 1 1
## 11 1 1 1 1 1
## 12 1 1 1 2 1
## 16 1 1 2 1 1
## 21 1 1 1 1 1
## 
## $chl
##      1   2   3   4   5
## 1  118 131 187 131 238
## 4  199 204 187 186 184
## 10 187 131 186 238 199
## 11 199 187 187 131 187
## 12 186 184 238 218 206
## 15 187 187 131 131 238
## 16 238 187 187 187 187
## 20 184 218 206 218 206
## 21 187 238 187 238 187
## 24 284 218 184 206 204

# Extract the "tall" matrix which stacks the imputations
nhanes.comp <- complete(nhanes.imp, "long", include = TRUE)
# Now, the .imp variable identifies which imputation each belongs to.
# 0 (included because of `include = TRUE` above) is the un-imputed data.
table(nhanes.comp$.imp)
## 
##  0  1  2  3  4  5 
## 25 25 25 25 25 25

# Let's visualize the distribution of imputed and observed values.
# `cci` returns logical whether its input is complete at each observation.
nhanes.comp$bmi.NA <- cci(nhanes$bmi)
# Note that the above uses the recycling properties of matrices/data.frame:
#  The `cci` call returns length 25; but because values are recycled to the total
#  number of rows in nhanes.comp, it replicates 6 times.
head(nhanes.comp[, c("bmi", "bmi.NA")])
##    bmi bmi.NA
## 1   NA  FALSE
## 2 22.7   TRUE
## 3   NA  FALSE
## 4   NA  FALSE
## 5 20.4   TRUE
## 6   NA  FALSE
library(ggplot2)
ggplot(nhanes.comp,
       aes(x = .imp, y = bmi, color = bmi.NA)) +
  geom_jitter(show.legend = FALSE,
              width = .1)
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

Because we used the default method (pmm), we’re guaranteed not to have any extreme values. That isn’t the case for other methods, so this is a nice check.

with

with is a general use function in R that replicates using attach, running code, then running detach. However, when we pass a mids object (which is the output of a mice() call) to with, it becomes a special mice version which properly handles the imputations.

with can take any valid R command as its second argument, and applies it to each imputed data set. Some examples:

with(nhanes.imp, mean(bmi))

## call :
## with.mids(data = nhanes.imp, expr = mean(bmi))
## 
## call1 :
## mice(data = nhanes, m = 5, method = "pmm")
## 
## nmis :
## [1]  0  9  8 10
## 
## analyses :
## [[1]]
## [1] 27.308
## 
## [[2]]
## [1] 26.484
## 
## [[3]]
## [1] 26.416
## 
## [[4]]
## [1] 25.924
## 
## [[5]]
## [1] 26.5

with(nhanes.imp, t.test(bmi ~ hyp))

## call :
## with.mids(data = nhanes.imp, expr = t.test(bmi ~ hyp))
## 
## call1 :
## mice(data = nhanes, m = 5, method = "pmm")
## 
## nmis :
## [1]  0  9  8 10
## 
## analyses :
## [[1]]
## 
##  Welch Two Sample t-test
## 
## data:  bmi by hyp
## t = 0.39664, df = 15.871, p-value = 0.6969
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -1.982589  2.894493
## sample estimates:
## mean in group 1 mean in group 2 
##        27.38095        26.92500 
## 
## 
## [[2]]
## 
##  Welch Two Sample t-test
## 
## data:  bmi by hyp
## t = -0.62971, df = 21.936, p-value = 0.5354
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -2.876939  1.536939
## sample estimates:
## mean in group 1 mean in group 2 
##           26.35           27.02 
## 
## 
## [[3]]
## 
##  Welch Two Sample t-test
## 
## data:  bmi by hyp
## t = -0.98907, df = 19.788, p-value = 0.3346
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -4.059166  1.449166
## sample estimates:
## mean in group 1 mean in group 2 
##          26.155          27.460 
## 
## 
## [[4]]
## 
##  Welch Two Sample t-test
## 
## data:  bmi by hyp
## t = -0.76475, df = 21.732, p-value = 0.4526
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -3.323811  1.533811
## sample estimates:
## mean in group 1 mean in group 2 
##          25.745          26.640 
## 
## 
## [[5]]
## 
##  Welch Two Sample t-test
## 
## data:  bmi by hyp
## t = 0.41685, df = 11.287, p-value = 0.6846
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -2.451633  3.601633
## sample estimates:
## mean in group 1 mean in group 2 
##          26.615          26.040

with(nhanes.imp, lm(bmi ~ age + chl))

## call :
## with.mids(data = nhanes.imp, expr = lm(bmi ~ age + chl))
## 
## call1 :
## mice(data = nhanes, m = 5, method = "pmm")
## 
## nmis :
## [1]  0  9  8 10
## 
## analyses :
## [[1]]
## 
## Call:
## lm(formula = bmi ~ age + chl)
## 
## Coefficients:
## (Intercept)          age          chl  
##    23.45528     -3.33037      0.05018  
## 
## 
## [[2]]
## 
## Call:
## lm(formula = bmi ~ age + chl)
## 
## Coefficients:
## (Intercept)          age          chl  
##    20.38148     -2.41047      0.05438  
## 
## 
## [[3]]
## 
## Call:
## lm(formula = bmi ~ age + chl)
## 
## Coefficients:
## (Intercept)          age          chl  
##     23.3579      -3.8341       0.0516  
## 
## 
## [[4]]
## 
## Call:
## lm(formula = bmi ~ age + chl)
## 
## Coefficients:
## (Intercept)          age          chl  
##    21.39073     -3.20281      0.05347  
## 
## 
## [[5]]
## 
## Call:
## lm(formula = bmi ~ age + chl)
## 
## Coefficients:
## (Intercept)          age          chl  
##    21.21720     -2.48164      0.04917

As an alternative, you could iterate over .imp inside nhanes.comp, performing the analysis on each. Generally I wouldn’t recommend this over with, but if necessary you could do this.

pool

Finally, for any method which produces an object having both coef() and vcov() methods (run methods(coef) or methods(vcov) to see), such as lm, aov or lmer (from lme4), the pool function can combine the individual coefficient estimates into one “pooled” estimate, based upon Rubin’s rules for pooling. As mentioned before, the point estimates are the strict average, but the standard errors are more complicated (the estimate is a combination of the average with-in imputation standard error, and the between imputation standard error).

mod1 <- with(nhanes.imp, lm(bmi ~ age + chl))
summary(pool(mod1))

##          term    estimate  std.error statistic        df      p.value
## 1 (Intercept) 21.96051494 3.50564555  6.264328 14.473290 1.784709e-05
## 2         age -3.05188010 1.08721540 -2.807061  8.933234 2.062019e-02
## 3         chl  0.05175869 0.01796149  2.881648 19.850697 9.271112e-03

Also useful is pool.r.squared for obtaining pooled \(R^2\) estimates of linear models.

pool.r.squared(mod1)

##           est      lo 95     hi 95        fmi
## R^2 0.4096055 0.09922382 0.6900042 0.06670379