---
title: "`mice`: Algorithmic convergence and inference pooling"
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>

---

This is the second exercise in the series.

The aim of this exercise is to enhance your understanding of multiple imputation, in general. You will learn how to pool the results of analyses performed on multiply-imputed data, how to approach different types of data and how to avoid the pitfalls researchers may fall into. The main objective is to increase your knowledge and understanding on applications of multiple imputation. 

Again, we start by loading (with `library()`) the necessary packages and fixing the random seed to allow for our outcomes to be replicable. 
```{r, message=FALSE, warning=FALSE}
library(mice)
set.seed(123)
```

All the best, 

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

---

# Multiple imputation with `mice`

---

## Varying $m$

---

**1. Vary the number of imputations. **

The number of imputed data sets can be specified by the `m = ...` argument. For example, to create just three imputed data sets, specify
```{r, cache = FALSE}
imp <- mice(nhanes, m = 3, print=F)
```

---

## Changing `predictorMatrix`

---

**2. Change the predictor matrix**

The predictor matrix is a square matrix that specifies the variables that are used to impute each incomplete variable. Let us have a look at the predictor matrix that was used
```{r, cache = FALSE}
imp$pred
```
Each variable in the data has a row and a column in the predictor matrix. A value `1` indicates that the column variable was used to impute the row variable. For example, the `1` at entry `[bmi, age]` indicates that variable `age` was used to impute the incomplete variable `bmi`. Note that the diagonal is zero because a variable is not allowed to impute itself. The row of `age` is redundant, because there were no missing values in `age`. Even though predictor relations are specified for `age`, `mice` will not use these relations because it will never overwrite the observed values with imputations. `mice` gives you complete control over the predictor matrix, enabling you to choose your own predictor relations. This can be very useful, for example, when you have many variables or when you have clear ideas or prior knowledge about relations in the data at hand. 

There are two ways in which you can create a predictor matrix in `mice`:

- A. You can grab the `predictorMatrix` from every object returned by `mice()`:

For example, we can use `mice()` to give you us initial predictor matrix, and change it afterwards, without running the algorithm. This can be done by typing
```{r, cache = FALSE}
ini <- mice(nhanes, maxit=0, print=F)
pred <- ini$pred
pred
```
The object `pred` contains the predictor matrix from an initial run of `mice` with zero iterations, specified by `maxit = 0`. 

- B. We can use `make.predictorMatrix()` to generate a predictor matrix from any incomplete data set.

For example, 
```{r}
pred <- make.predictorMatrix(nhanes)
pred
```

Altering the predictor matrix and returning it to the mice algorithm is very simple. For example, the following code removes the variable `hyp` from the set of predictors, but still leaves it to be predicted by the other variables.
```{r, cache = FALSE}
pred[, "hyp"] <- 0
pred
```

Use your new predictor matrix in `mice()` as follows
```{r, cache = FALSE}
imp <- mice(nhanes, predictorMatrix =  pred, print = F)
```

There is a `quickpred()` function that applies a quick selection procedure of predictors, which can be handy for datasets containing many variables. See `?quickpred` for more info. Selecting predictors according to data relations with a minimum correlation of $\rho=.30$ can be done by
```{r, cache = FALSE}
ini <- mice(nhanes, pred=quickpred(nhanes, mincor=.3), print=F)
ini$pred
```
For large predictor matrices, it can be useful to export them to dedicated spreadsheet software like e.g. Microsoft Excel for easier configuration (e.g. see the [`xlsx` package](https://cran.r-project.org/web/packages/xlsx/index.html) for easy exporting and importing of Excel files). Importing data is straightforward in `RStudio` through `File` > `Import Dataset`. 

---

## Inspect convergence

---

**3. Inspect the convergence of the algorithm**

The `mice()` function implements an iterative Markov Chain Monte Carlo type of algorithm. Let us have a look at the trace lines generated by the algorithm to study convergence:
```{r, cache = FALSE}
imp <- mice(nhanes, print=F)
plot(imp)
```

The plot shows the mean (left) and standard deviation (right) of the imputed values only. In general, we would like the streams to intermingle and be free of any trends at the later iterations. We inspect trends for the imputed values alone, because the observed data does not change.

The `mice` algorithm uses random sampling, and therefore, the results will be (perhaps slightly) different if we repeat the imputations with different seeds. In order to get exactly the same result, use the `seed` argument 
```{r, cache = FALSE}
imp <- mice(nhanes, seed=123, print=F)
```
where `123` is some arbitrary number that you can choose yourself. Rerunning this command will always yields the same imputed values.

---

## Change `method` 

---

**4. Change the imputation method**

For each column, the algorithm requires a specification of the imputation method. To see which method was used by default:
```{r, cache = FALSE}
imp$meth
```
The variable `age` is complete and therefore not imputed, denoted by the `""` empty string. The other variables have method `pmm`, which stands for *predictive mean matching*, the default in
`mice` for numerical and integer data.  

In reality, the `nhanes` data are better described a as mix of numerical and categorical data. Let us take a look at the `nhanes2` data frame:
```{r, cache = FALSE}
summary(nhanes2)
```
and the structure of the data frame
```{r, cache = FALSE}
str(nhanes2)
```
Variable `age` consists of 3 age categories, while variable `hyp` is binary. The `mice()` function takes these properties automatically into account. Impute the `nhanes2` dataset
```{r, cache = FALSE}
imp <- mice(nhanes2, print=F)
imp$meth
```
Notice that `mice` has set the imputation method for
variable `hyp` to `logreg`, which implements multiple imputation by *logistic
regression*. 

An up-to-date overview of the methods in mice can be found by
```{r, warning=FALSE}
methods(mice)
```

Let us change the imputation method for `bmi` to Bayesian normal linear regression imputation
```{r, cache = FALSE}
meth <- make.method(nhanes2)
meth
meth["bmi"] <- "norm"
meth
```
and run the imputations again.
```{r, cache = FALSE}
imp <- mice(nhanes2, meth = meth, print=F)
```
We may now again plot trace lines to study convergence
```{r, cache = FALSE}
plot(imp)
```

---

## Add more iterations

---

**5. Extend the number of iterations**

Though using just five iterations (the default) often works well in practice, we need to extend the number of iterations of the `mice` algorithm to confirm that there is no trend and that the trace lines intermingle well. We can increase the number of iterations to 40 by running 35 additional iterations using the `mice.mids()` function. 
```{r, cache = FALSE}
imp40 <- mice.mids(imp, maxit=35, print=F)
plot(imp40)
```

---

## Diagnostic evaluation

---

**6. Further diagnostic checking. Use function `stripplot()`.**

Generally, one would prefer for the imputed data to be plausible values, i.e. values that could have
been observed if they had not been missing. In order to form an idea about plausibility, one may check the imputations and compare them against the observed values. If we are willing to assume that
the data are missing completely at random (MCAR), then the imputations should have the same distribution as the observed data. In general, distributions may be different because the missing
data are MAR (or even MNAR). However, very large discrepancies need to be screened. Let us plot the observed and imputed data of `chl` by
```{r, cache = FALSE}
stripplot(imp, chl~.imp, pch=20, cex=2)
```

The convention is to plot observed data in blue and the imputed data in red. The figure graphs the data values of `chl` before and after imputation. Since the PMM method draws imputations from the observed data, imputed values have the same gaps as in the observed data, and are always within the range of the observed data. The figure indicates that the distributions of the imputed and the observed values are similar. The observed data have a particular feature that, for some reason, thedata cluster around the value of 187. The imputations reflect this feature, and are close to the data. Under MCAR, univariate distributions of the observed and imputed data are expected to be identical. Under MAR, they can be different, both in location and spread, but their multivariate distribution is assumed to be identical. There are many other ways to look at the imputed data.

The following command creates a simpler version of the graph from the previous step and adds the plot for `bmi`. 
```{r, cache = FALSE}
stripplot(imp)
```

Remember that `bmi` was imputed by Bayesian linear regression and (the range of) imputed values may therefore be different than observed values. 

---

# Repeated analysis in mice

---

**7. Perform the following regression analysis on the multiply imputed data and assign the result to object `fit`. **

\[
\text{bmi} = \beta_0 + \beta_1 \text{chl} + \epsilon
\]
```{r, cache = FALSE}
fit <- with(imp, lm(bmi ~ chl))
fit
```
The `fit` object contains the regression summaries for each data set. 
The new object `fit` is actually of class `mira` (*multiply imputed repeated analyses*).
```{r, cache = FALSE}
class(fit)
```
Use the `ls()` function to what out what is in the object. 
```{r, cache = FALSE}
ls(fit)
```
Suppose we want to find the regression model fitted to the
second imputed data set. It can be found as
```{r, cache = FALSE}
summary(fit$analyses[[2]])
```

---

**8. Pool the analyses from object `fit`. **

Pooling the repeated regression analyses can be done simply by typing
```{r, cache = FALSE}
pool.fit <- pool(fit)
summary(pool.fit)
pool.fit
```
which gives the relevant pooled regression coefficients and
parameters, as well as the fraction of information about the
coefficients missing due to nonresponse (`fmi`) and the proportion of the variation attributable to the missing data (`lambda`). The pooled fit object is of class `mipo`, which stands for *multiply imputed pooled object*. 

`mice` can to pool many analyses from a variety of packages for you. For flexibility and in order to run custom pooling functions, mice also incorporates a function `pool.scalar()` which pools univariate estimates of $m$ repeated complete data analysis conform Rubin's pooling rules (Rubin, 1987, paragraph 3.1) 

---

# 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)

---

**- End of exercise**