---
title: "Ad hoc methods and `mice`"
author: "Gerko Vink and Stef van Buuren"
date: "Multiple Imputation in Practice"
output: 
  html_document:
    toc: true
    toc_depth: 5
    toc_float: true
    number_sections: false
---
<style type="text/css">

body{ /* Normal  */
      font-size: 12px;
  }
td {  /* Table  */
  font-size: 12px;
}
h1.title {
  font-size: 18px;
  color: DarkBlue;
}
h1 { /* Header 1 */
  font-size: 18px;
}
h2 { /* Header 2 */
    font-size: 18px;
}
h3 { /* Header 3 */
  font-size: 18px;
}
code.r{ /* Code block */
    font-size: 12px;
}
pre { /* Code block - determines code spacing between lines */
    font-size: 14px;
}
</style>

---

**The required `mice` version for this practical exercise is `mice > 3.1.0`. We'll install the latest version in this exercise.**

---

This is the first exercise in the series. It will give you an introduction to the `R`-package `mice`, an open-source tool for flexible imputation of incomplete data. Over the last decade, `mice` has become an important piece of imputation software, offering a very flexible environment for dealing with incomplete data. Moreover, the ability to integrate `mice` with other packages in `R`, and vice versa, offers many options for applied researchers.

The aim of this introduction is to enhance your understanding of multiple imputation, in general. You will learn how to multiply impute simple datasets and how to obtain the imputed data for further analysis. The main objective is to increase your knowledge and understanding on applications of multiple imputation. 

All the best, 

[Gerko](https://www.gerkovink.com) and [Stef](http://www.stefvanbuuren.name)

---

`mice` uses `R`'s random number generator to draw values with a probabilistic nature. Therefore, each time we use `mice` we will get slightly different results. To avoid this, we can fix the seed value of the random number generator. 
```{r}
set.seed(123)
```

With this seed you'll get the exact same results if you follow the steps in this document. If you obtain different results, you have changed the order of the steps either by adding or re-running a step. 

---

# Installing the latest version of `mice`
We can use package `devtools` (which needs to be installed) to directly *grab* the latest version from the [`mice` Github Page](https://github.com/stefvanbuuren/mice) and compile the `mice` package from source. 
```{r eval=FALSE}
install.packages("mice")
```
You can update all dependencies with more recent versions, if asked. 

---

# Working with mice

---

**1. Open `R` and load the package `mice`**

```{r, message=FALSE, warning=FALSE}
library(mice)
```
The version number for your `mice` can be found by running
```{r}
version()
```

---

**2.  Inspect the incomplete data**

The `mice` package contains several datasets. Once the package is loaded, these datasets can be used. Have a look at the `nhanes` dataset (Schafer, 1997, Table 6.14) by typing
```{r}
nhanes
```
 The `nhanes` dataset is a small data set with non-monotone
 missing values. It contains 25 observations on four variables:
 *age group*, *body mass index*, *hypertension*
 and *cholesterol (mg/dL)*.

To learn more about the data, use one of the two following help commands:
```{r, cache = FALSE}
help(nhanes)
?nhanes
```

---

**3.  Get an overview of the data by the `summary()` command:**

```{r, cache = FALSE}
summary(nhanes)
```

Using `summary()` on data sets is often informative, because the distributional information (continuous variables) or the frequency distribution (factors) for every column in your data frame is printed to the `R` console. However, if there are too many variables, a step-by-step approach may be more useful. 

---

**4.  Inspect the missing data pattern**

Check the missingness pattern for the `nhanes` dataset
```{r, cache = FALSE}
md.pattern(nhanes)
```
The missingness pattern shows that there are 27 missing values in total: 10 for `chl` , 9 for `bmi` and 8 for `hyp`. Moreover, there are thirteen completely observed rows, four rows with 1 missing, one row with 2 missings and seven rows with 3 missings. Looking at the missing data pattern is always useful (but may be difficult for datasets with many variables). It can give you an indication on how much information is missing and how the missingness is distributed.

---

# Ad Hoc imputation methods

---

**5. Form a regression model where `age` is predicted from `bmi`. **

We can use the `with()` family of functions for this. The following function call 
```{r, cache=TRUE}
fit <- with(nhanes, lm(age ~ bmi))
```

evaluates `with(data, expression)`, so it evaluates the linear model `lm(age ~ bmi)` on data set `nhanes`. The resulting object `fit` is identical to the output from `lm(age ~ bmi, data = nhanes)`. We learn the `with()` function now, because we need it later when we start evaluating analytical models on multiply imputed data sets. 

If we ask the summary of the fitted regression analysis, we obtain:
```{r}
summary(fit)
```
No significant effect for `bmi` when we model the `age` variable. 

---

**6. Impute the missing data in the `nhanes` dataset with mean imputation. **
The following function call imputes the mean (`method = "mean"`) for every incomplete column in the `nhanes` data set and returns a single (`m = 1`) imputed data set. The algorithm - i.e. the method that generates the imputations - has been given a single iteration (`maxit = 1`) to reach convergence. We'll dive into the specifics of algorithmic convergence with `mice` in the next exercise.

```{r, cache=TRUE}
imp <- mice(nhanes, method = "mean", m = 1, maxit = 1)
```
The imputations are now done. Running only a single imputation is practically efficient, as substituting each missing data multiple times with the observed data mean would not make any sense (the inference would be equal, no matter which imputed dataset we would analyze). Likewise, more iterations would be computationally inefficient as the *observed* data mean does not change based on our imputations. We named the imputed object `imp` following the convention used in `mice`, but if you wish you can name it anything you'd like. 

---

**7. Explore the imputed data with the `complete()` function. What do you think the variable means are? What happened to the regression equation after imputation?**

We use the function `complete()`, which by default returns the first completed data set. Since we only have a single imputation for every missing datum, this makes sense and we do not have to change the default behavior of `complete()`.
```{r, cache = FALSE}
complete(imp)
```
We see the repetitive numbers `26.5625` for `bmi`, `1.2352594` for `hyp`, and `191.4` for `chl`. These can be confirmed as the means of the respective variables (columns):
```{r, cache = FALSE}
colMeans(nhanes, na.rm = TRUE)
```
We've seen during the inspection of the missing data pattern that variable `age` has no missings. Therefore nothing is imputed for `age` because we would not want to alter the observed (and bonafide) values. 

To inspect the regression model with the imputed data, run:
```{r}
fit <- with(imp, lm(age ~ bmi))
summary(fit)
```
It is clear that our inference did not change, but then again this is not surprising as variable `bmi` is more-or-less normally distributed and we are just adding weight to the mean.
```{r}
densityplot(nhanes$bmi)
```

---

**8. Impute the missing data in the `nhanes` dataset with regression imputation. **

We can use the same function call as under exercise 7, with `method = "norm.predict"`, which yields predictions from the normal linear regression model. 
```{r, cache=TRUE}
imp <- mice(nhanes, method = "norm.predict", m = 1, maxit = 1)
```
The imputations are now done. This code imputes the missing values in the data set by the regression imputation method. The argument `method = "norm.predict"` first fits a regression model for each observed value, based on the corresponding values in other variables and then imputes the missing values with the fitted (predicted) values from the normal linear regression model. 

---

**9. Again, inspect the completed data and investigate the imputed data regression model. **

The completed data:
```{r, cache = FALSE}
complete(imp)
```
The repetitive numbering we saw under mean imputation is now gone when we impute the conditional mean - i.e. the expectation of `age` for every given `bmi`. We have now obtained a more natural looking set of imputations: instead of filling in the same `bmi` for all ages, we now take `age` (as well as `hyp` and `chl`) into account when imputing `bmi`. 

To inspect the regression model with the imputed data, run:
```{r}
fit <- with(imp, lm(age ~ bmi))
summary(fit)
```
It is clear that our inference has changed. In fact, we extrapolated (part of) the regression model for the observed data to missing data in `bmi`. In other words; the relation (read: information) gets stronger and we've obtained more observations that conform exactly to the relation in the observed data. From an inferential statistics viewpoint, this approach would ***only*** be valid if we have definitive proof that the unobserved values would exactly conform to the observed data. If this assumption does not hold, `norm.predict` creates too little variation in our data set and we can not trust the resulting standard errors and p-values. 

---

**10. Impute the missing data in the `nhanes` dataset with stochastic regression imputation.** 
With stochastic regression imputation, an error term is added to the predicted values, such that the imputations show variation around the regression line. The errors are normally distributed with mean 0 and variance equal to the residual variance. 
```{r, cache=TRUE}
imp <- mice(nhanes, method = "norm.nob", m = 1, maxit = 1)
```
The imputations are now done. This code imputes the missing values in the data set by the stochastic regression imputation method. The function does not incorporate the variability of the regression weights, so it is not 'proper' in the sense of Rubin (1987). For small samples, the variability of the imputed data will be underestimated. 

---

**11. Again, inspect the completed data and investigate the imputed data regression model. **

```{r, cache = FALSE}
complete(imp)
```
We have once more obtained a more natural looking set of imputations, where instead of filling in the same `bmi` for all ages, we now take `age` (as well as `hyp` and `chl`) into account when imputing `bmi`. We also add a random error to allow for our imputations to be off the regression line.  

To inspect the regression model with the imputed data, run:
```{r}
fit <- with(imp, lm(age ~ bmi))
summary(fit)
```

---

**12. Re-run the stochastic imputation model with seed `123` and verify if your results are the same as the ones below**

```{r, echo=FALSE, warning=FALSE, message=FALSE}
imp <- mice(nhanes, method = "norm.nob", m = 1, maxit = 1, seed = 123, print=F)
fit <- with(imp, lm(age ~ bmi))
summary(fit)
```
The imputation procedure uses random sampling, and therefore, the results will be (perhaps slightly) different if we repeat the imputations. In order to get exactly the same result, you can use the seed argument
```{r, eval=FALSE}
imp <- mice(nhanes, method = "norm.nob", m = 1, maxit = 1, seed = 123)
fit <- with(imp, lm(age ~ bmi))
summary(fit)
```
where 123 is some arbitrary number that you can choose yourself. Re-running this command will always yields the same imputed values. The ability to replicate one's findings exactly is considered essential in today's reproducible science. 

---

# Multiple imputation

---

**13.  Let us impute the missing data in the `nhanes` dataset**
To do multiple imputation, we can simply call `mice()` on our data set:
```{r, cache=TRUE}
imp <- mice(nhanes)
```
The imputations are now done. As you can see, the algorithm ran for 5 iterations (the default) and presented us with 5 imputations for each missing datum. For the rest of this document we will omit printing of the iteration cycle when we run `mice`. We do so by adding `print=F` to the `mice` call. 

```{r}
imp
```

The object `imp` contains a multiply imputed data set (of class `mids`). It encapsulates all information from imputing the `nhanes` dataset, such as the original data, the imputed values, the number of missing values, number of iterations, and so on.

To obtain an overview of the information stored in the object
`imp`, use the `attributes()` function:
```{r, cache = FALSE}
attributes(imp)
```
For example, the original data are stored as
```{r, cache = FALSE}
imp$data
```
and the imputations are stored as
```{r, cache = FALSE}
imp$imp
```

---

**14. Extract the completed data**

By default, `mice()` calculates five (*m* = 5) imputed data sets. In order to get the third imputed data set, use the `complete()` function
```{r, cache=TRUE}
c3 <- complete(imp, 3) 
md.pattern(c3)
```
The collection of the $m$ imputed data sets can be exported by function
`complete()` in long, broad and repeated formats. 
For example, 
```{r, cache=TRUE}
c.long <- complete(imp, "long")  
c.long
```
and
```{r, cache=TRUE}
c.broad <- complete(imp, "broad")
c.broad
```
are completed data sets in long and broad format, respectively. See `?complete` for more detail.

---

# Conclusion

We have seen that (multiple) imputation is straightforward with `mice`. However, don't let the simplicity of the software fool you into thinking that the problem itself is also straightforward. In the next exercise we will therefore explore how the mice package can flexibly provide us the tools to assess and control the imputation of missing data.

---

# References

Rubin, D. B. *Multiple imputation for nonresponse in surveys*. John Wiley & Sons, 1987. [Amazon](http://www.amazon.com/Multiple-Imputation-Nonresponse-Surveys-Donald/dp/0471655740/ref=sr_1_1?ie=UTF8&qid=1434466788&sr=8-1&keywords=Multiple+imputation+for+nonresponse+in+surveys)

Schafer, J.L. (1997). *Analysis of Incomplete Multivariate Data*. London: Chapman & Hall. Table 6.14. [Amazon](http://www.amazon.com/Incomplete-Multivariate-Monographs-Statistics-Probability/dp/0412040611/ref=sr_1_1?ie=UTF8&qid=1434466828&sr=8-1&keywords=Analysis+of+Incomplete+Multivariate+Data)

Van Buuren, S. and Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. *Journal of Statistical Software*, 45(3), 1-67. [pdf](http://www.jstatsoft.org/v45/i03/paper)

---

**- End of exercise**

---