Missing data can be a not so trivial problem when analysing a dataset and accounting for it is usually not so straightforward either. If the amount of missing data is very small relatively to the size of the dataset, then leaving out the few samples with missing features may be the best strategy in order not to bias the analysis, however leaving out available datapoints deprives the data of some amount of information and depending on the situation you face, you may want to look for other fixes before wiping out potentially useful datapoints from your dataset.

Update

A simplified approach to impute missing data with MICE package can be found there: Handling missing data with MICE package; a simple approach.

While some quick fixes such as mean-substitution may be fine in some cases, such simple approaches usually introduce bias into the data, for instance, applying mean substitution leaves the mean unchanged (which is desirable) but decreases variance, which may be undesirable.

The mice package in R, helps you imputing missing values with plausible data values. These plausible values are drawn from a distribution specifically designed for each missing datapoint.

In this post we are going to impute missing values using a the airquality dataset (available in R).
For the purpose of the article I am going to remove some datapoints from the dataset.

data <- airquality
data[4:10,3] <- rep(NA,7)
data[1:5,4] <- NA

As far as categorical variables are concerned, replacing categorical variables is usually not advisable. Some common practice include replacing missing categorical variables with the mode of the observed ones, however, it is questionable whether it is a good choice. Even though in this case no datapoints are missing from the categorical variables, we remove them from our dataset (we can add them back later if needed) and take a look at the data using summary().

data <- data[-c(5,6)]
summary(data)

     Ozone           Solar.R           Wind             Temp      
 Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :57.00  
 1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:73.00  
 Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
 Mean   : 42.13   Mean   :185.9   Mean   : 9.806   Mean   :78.28  
 3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
 Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
 NA's   :37       NA's   :7       NA's   :7        NA's   :5  

Apparently Ozone is the variable with the most missing datapoints. Below we are going to dig deeper into the missing data patterns.

Quick classification of missing data

There are two types of missing data:

Assuming data is MCAR, too much missing data can be a problem too. Usually a safe maximum threshold is 5% of the total for large datasets. If missing data for a certain feature or sample is more than 5% then you probably should leave that feature or sample out. We therefore check for features (columns) and samples (rows) where more than 5% of the data is missing using a simple function

pMiss <- function(x){sum(is.na(x))/length(x)*100}
apply(data,2,pMiss)
apply(data,1,pMiss)

    Ozone   Solar.R      Wind      Temp 
24.183007  4.575163  4.575163  3.267974 

  [1]  25  25  25  50 100  50  25  25  25  50  25   0   0   0   0   0   0   0   0   0   0
 [22]   0   0   0  25  25  50   0   0   0   0  25  25  25  25  25  25   0  25   0   0  25
 [43]  25   0  25  25   0   0   0   0   0  25  25  25  25  25  25  25  25  25  25   0   0
 [64]   0  25   0   0   0   0   0   0  25   0   0  25   0   0   0   0   0   0   0  25  25
 [85]   0   0   0   0   0   0   0   0   0   0   0  25  25  25   0   0   0  25  25   0   0
[106]   0  25   0   0   0   0   0   0   0  25   0   0   0  25   0   0   0   0   0   0   0
[127]   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
[148]   0   0  25   0   0   0

We see that Ozone is missing almost 25% of the datapoints, therefore we might consider either dropping it from the analysis or gather more measurements. The other variables are below the 5% threshold so we can keep them. As far as the samples are concerned, missing just one feature leads to a 25% missing data per sample. Samples that are missing 2 or more features (>50%), should be dropped if possible.

Using mice for looking at missing data pattern

The mice package provides a nice function md.pattern() to get a better understanding of the pattern of missing data

library(mice)
md.pattern(data)

    Temp Solar.R Wind Ozone   
104    1       1    1     1  0
 34    1       1    1     0  1
  4    1       0    1     1  1
  3    1       1    0     1  1
  3    0       1    1     1  1
  1    1       0    1     0  2
  1    1       1    0     0  2
  1    1       0    0     1  2
  1    0       1    0     1  2
  1    0       0    0     0  4
       5       7    7    37 56

The output tells us that 104 samples are complete, 34 samples miss only the Ozone measurement, 4 samples miss only the Solar.R value and so on.

A perhaps more helpful visual representation can be obtained using the VIM package as follows

library(VIM)
aggr_plot <- aggr(data, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(data), cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))

Rplot
The plot helps us understanding that almost 70% of the samples are not missing any information, 22% are missing the Ozone value, and the remaining ones show other missing patterns. Through this approach the situation looks a bit clearer in my opinion.

Another (hopefully) helpful visual approach is a special box plot

marginplot(data[c(1,2)])

Rplot02
Obviously here we are constrained at plotting 2 variables at a time only, but nevertheless we can gather some interesting insights.
The red box plot on the left shows the distribution of Solar.R with Ozone missing while the blue box plot shows the distribution of the remaining datapoints. Likewhise for the Ozone box plots at the bottom of the graph.
If our assumption of MCAR data is correct, then we expect the red and blue box plots to be very similar.

Imputing the missing data

The mice() function takes care of the imputing process

tempData <- mice(data,m=5,maxit=50,meth='pmm',seed=500)
summary(tempData)

Multiply imputed data set
Call:
mice(data = data, m = 5, method = "pmm", maxit = 50, seed = 500)
Number of multiple imputations:  5
Missing cells per column:
  Ozone Solar.R    Wind    Temp 
     37       7       7       5 
Imputation methods:
  Ozone Solar.R    Wind    Temp 
  "pmm"   "pmm"   "pmm"   "pmm" 
VisitSequence:
  Ozone Solar.R    Wind    Temp 
      1       2       3       4 
PredictorMatrix:
        Ozone Solar.R Wind Temp
Ozone       0       1    1    1
Solar.R     1       0    1    1
Wind        1       1    0    1
Temp        1       1    1    0
Random generator seed value:  500 

A couple of notes on the parameters:

If you would like to check the imputed data, for instance for the variable Ozone, you need to enter the following line of code

tempData$imp$Ozone

      1  2   3   4   5
5    13 20  28  12   9
10    7 16  28  14  20
25    8 14  14   1   8
26    9 19  32   8  37
...

The output shows the imputed data for each observation (first column left) within each imputed dataset (first row at the top).
If you need to check the imputation method used for each variable, mice makes it very easy to do

tempData$meth

  Ozone Solar.R    Wind    Temp 
  "pmm"   "pmm"   "pmm"   "pmm"

Now we can get back the completed dataset using the complete() function. It is almost plain English:

completedData <- complete(tempData,1)

The missing values have been replaced with the imputed values in the first of the five datasets. If you wish to use another one, just change the second parameter in the complete() function.

Inspecting the distribution of original and imputed data

Let’s compare the distributions of original and imputed data using a some useful plots.
First of all we can use a scatterplot and plot Ozone against all the other variables

xyplot(tempData,Ozone ~ Wind+Temp+Solar.R,pch=18,cex=1)

Here it is
Rplot03
What we would like to see is that the shape of the magenta points (imputed) matches the shape of the blue ones (observed). The matching shape tells us that the imputed values are indeed “plausible values”.
Another helpful plot is the density plot:

densityplot(tempData)

Rplot04
The density of the imputed data for each imputed dataset is showed in magenta while the density of the observed data is showed in blue. Again, under our previous assumptions we expect the distributions to be similar.
Another useful visual take on the distributions can be obtained using the stripplot() function that shows the distributions of the variables as individual points

stripplot(tempData, pch = 20, cex = 1.2)

Rplot05

Pooling

Suppose that the next step in our analysis is to fit a linear model to the data. You may ask what imputed dataset to choose. The mice package makes it again very easy to fit a a model to each of the imputed dataset and then pool the results together

modelFit1 <- with(tempData,lm(Temp~ Ozone+Solar.R+Wind))
summary(pool(modelFit1))

                     est         se         t       df     Pr(>|t|)
(Intercept) 72.812078768 2.95380500 24.650266 84.18464 0.000000e+00
Ozone        0.163094287 0.02607674  6.254397 57.78569 5.236295e-08
Solar.R      0.009679676 0.00789576  1.225933 37.48960 2.278691e-01
Wind        -0.352582008 0.21639828 -1.629320 92.89136 1.066321e-01
                   lo 95       hi 95 nmis       fmi    lambda
(Intercept) 66.938301817 78.68585572   NA 0.1477818 0.1277731
Ozone        0.110891894  0.21529668   37 0.2155848 0.1888975
Solar.R     -0.006311604  0.02567095    7 0.3004189 0.2640672
Wind        -0.782312735  0.07714872    0 0.1300747 0.1115442

The variable modelFit1 containts the results of the fitting performed over the imputed datasets, while the pool() function pools them all together. Apparently, only the Ozone variable is statistically significant.
Note that there are other columns aside from those typical of the lm() model: fmi contains the fraction of missing information while lambda is the proportion of total variance that is attributable to the missing data. For more information I suggest to check out the paper cited at the bottom of the page.

Remember that we initialized the mice function with a specific seed, therefore the results are somewhat dependent on our initial choice. To reduce this effect, we can impute a higher number of dataset, by changing the default m=5 parameter in the mice() function as follows

tempData2 <- mice(data,m=50,seed=245435)
modelFit2 <- with(tempData2,lm(Temp~ Ozone+Solar.R+Wind))
summary(pool(modelFit2))

                     est          se         t       df     Pr(>|t|)
(Intercept) 73.156084276 2.803010282 26.099114 129.3154 0.000000e+00
Ozone        0.166242781 0.024926976  6.669192 118.4408 8.645631e-10
Solar.R      0.009046835 0.007374103  1.226839 114.5471 2.223989e-01
Wind        -0.382700790 0.202976584 -1.885443 136.6735 6.149264e-02
                   lo 95       hi 95 nmis        fmi    lambda
(Intercept) 67.610387851 78.70178070   NA 0.11141367 0.0977762
Ozone        0.116882484  0.21560308   37 0.16290744 0.1488906
Solar.R     -0.005560458  0.02365413    7 0.18096774 0.1667911
Wind        -0.784081566  0.01867999    0 0.07425875 0.0608104

After having taken into account the random seed initialization, we obtain (in this case) more or less the same results as before with only Ozone showing statistical significance.

A gist with the full code for this post can be found here.

Thank you for reading this post, leave a comment below if you have any question.

Note: I learnt this technique in a paper entitled mice: Multivariate Imputation by Chained Equations in R by Stef van Buuren. It is a great paper and I highly recommend to read it if you are interested in multiple imputation!