---
title: "Multilevel Imputation"
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: 8px;
}
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 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) # Imputation
library(lme4) # Multilevel modeling
library(mitml) # Multilevel imputation
library(lattice) # Plotting device
library(broom.mixed) # Helper functions for multilevel estimation
set.seed(123)
```

All the best, 

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

---

# The `popularity` data set

We are going to work with the popularity data from Joop Hox (2010). The variables in this data set are described as follows:

<table>
<tbody>
<tr class="odd">
<td align="left"><strong>pupil</strong></td>
<td align="left"> $\quad$ Pupil number within class</td>
</tr>
<tr class="even">
<td align="left"><strong>class</strong></td>
<td align="left"> $\quad$ Class number</td>
</tr>
<tr class="odd">
<td align="left"><strong>extrav</strong></td>
<td align="left"> $\quad$ Pupil extraversion</td>
</tr>
<tr class="even">
<td align="left"><strong>sex</strong></td>
<td align="left"> $\quad$ Pupil gender</td>
</tr>
<tr class="odd">
<td align="left"><strong>texp</strong></td>
<td align="left"> $\quad$ Teacher experience (years)</td>
</tr>
<tr class="even">
<td align="left"><strong>popular</strong></td>
<td align="left"> $\quad$ Popularity as rated by the pupil's peers</td>
</tr>
<tr class="odd">
<td align="left"><strong>popteach</strong></td>
<td align="left"> $\quad$ Popularity as rated by the teacher</td>
</tr>
</tbody>
</table>

---

## Obtaining the popularity data
A workspace with complete and incomplete versions of the popularity data [can be obtained here](https://www.gerkovink.com/mimp/popular.RData) or can be loaded into the Global Environment by running:
```{r eval = "FALSE"}
con <- url("https://www.gerkovink.com/mimp/popular.RData")
load(con)
```

This workspace contains several datasets and functions that, when loaded, are available to you in R. If you’d like to see what is in your Gloval Environment after importing the workspace: run the following code
```{r cache = FALSE}
ls()
```

The `popular` data set is the complete version. The other `popXXX` sets are incomplete versions of this data set with varying levels of complexity. The `icc()` function is a generic function from package `multilevel` (function `ICC1()`). 

Having the complete version of the data (the `popular` set) is a luxury that we do not have in practice. However, for educational purposes we can use this simulated set and its incomplete versions to inspect, evaluate and verify our imputation methodology and approach. 

---

The dataset `popNCR` is a variation on the Hox (2010) data, where the missingness in the variables is either missing at random (MAR) or missing not at random (MNAR). We will be working with this `popNCR` set in this exercise.

---

## Inspecting the `popNCR` data

---

**1. Check with the functions `head()`, `dim()` - alternatively one could use `nrow()` and `ncol()` instead of `dim()` - and `summary()` how large the dataset is, of which variables the data frame consists and if there are missing values in a variable.**

```{r cache = FALSE}
head(popNCR)
dim(popNCR)
nrow(popNCR)
ncol(popNCR)
summary(popNCR)
```
The data set has 2000 rows and 7 columns (variables). The variables `extrav`, `sex`, `texp`, `popular` and `popteach` contain missings. About a quarter of these variables is missing, except for `texp` where 50 % is missing.

---

**2. As we have seen before, the function `md.pattern()` is used to display all different missing data patterns. How many different missing data patterns are present in the `popNCR` dataframe and which patterns occur most frequently in the data?**

```{r cache = FALSE, echo=TRUE}
md.pattern(popNCR, plot = FALSE)
```
```{r eval = FALSE}
md.pattern(popNCR)
```
![](fig/mdpattern.png){ width=40% }

There are 32 unique patterns. The pattern where everything is observed and the pattern where only `texp` is missing occur most frequently.

If we omit `texp`, then half the patterns disappear:

```{r cache = FALSE, echo=TRUE}
md.pattern(popNCR[, -5], plot = FALSE)
```
```{r eval = FALSE}
md.pattern(popNCR[, -5])
```
![](fig/mdpattern2.png){ width=40% }

Without the fifth column `[, -5]` which is `texp`, there are only 16 patterns.

---

# Where does the missingness come from?

---

## Missingness inspection of `popular`

---

**3. Let’s focus more precisely on the missing data patterns. Does the missing data of `popular` depend on `popteach`? One could for example check this by making a histogram of `popteach` separately for the pupils with known popularity and missing popularity.**

In R the missingness indicator
```{r cache = FALSE, eval = FALSE}
is.na(popNCR$popular)
```
is a dummy variable of the same length as `popular` with value 0 (`FALSE`) for observed pupil popularity and 1 (`TRUE`) for missing pupil popularity. A histogram can be made with the function `histogram()`. The code for a conditional histogram of `popteach` given the missingness indicator for popular is
```{r cache = FALSE}
histogram(~ popteach | is.na(popular), data=popNCR)
```

We do see that the histogram for the missing `popular` (`TRUE`) is further to the right than the histogram for observed `popular` (`FALSE`). This would indicate a right-tailed MAR missingness. In fact this is exactly what happens, because I (Gerko) created the missingness in these data. However, we can make it observable by examining the relations between the missingness in `popular` and the observed data in `popteach`. This is a clear indication that the missingness in `popular` depends on  the observed values for `popteach`. This finding makes the `MCAR` assumption improbable. 

---

## Missingness inspection of the other variables. 

---

**4. Does the missingness of the other incomplete variables depend on `popteach`? If yes, what is the direction of the relation?**

```{r cache = FALSE}
histogram(~ popteach | is.na(sex), data = popNCR)  
```

There seems to be a left-tailed relation between `popteach` and the missingness in `sex`.
```{r cache = FALSE}
histogram(~ popteach | is.na(extrav), data = popNCR)
```

There also seems to be a left-tailed relation between `popteach` and the missingness in `extrav`.
```{r cache = FALSE}
histogram(~ popteach | is.na(texp), data = popNCR)
```

There seems to be no observable relation between `popteach` and the missingness in `texp`. It might be random or even MNAR.

---

**5. Find out if the missingness in teacher popularity depends on pupil popularity.**

```{r cache = FALSE}
histogram(~ popular | is.na(popteach), data = popNCR)
```

Yes: there is a dependency. The relation seems to be right-tailed.

---

We now know that some variables (columns) seem to be related to the missingness in other variables. This is a clear indication that for the relevant variables, the predictor relations during multiple imputation should be set. 

---

**6. Have a look at the intraclass correlation (ICC) for the incomplete variables `popular`, `popteach` and `texp`.**

```{r cache = FALSE}
icc(aov(popular ~ as.factor(class), data = popNCR))
icc(aov(popteach ~ class, data = popNCR))
icc(aov(texp ~ class, data = popNCR))
```
Please note that the function `icc()` comes from the package `multilevel` (function `ICC1()`), but is included in the workspace `popular.RData`. Make a note of the ICC’s, you’ll need them later.

---

**7. Do you think it is necessary to take the multilevel structure into account?**

YES! There is a strong cluster structure going on. If we ignore the clustering in our imputation model, we may run into invalid inference. To stay as close to the true data model, we must take the cluster structure into account during imputation.

---

# Multiple Imputation of the `popNCR` data set

---

**8. Impute the `popNCR` dataset with `mice` using imputation method `norm` for `popular`, `popteach`, `texp` and `extrav`. Exclude `class` as a predictor ***for all variables***. Name the multiply imputed data set (`mids`) resulting from the `mice` call `imp1`. Use $m=10$ imputations and $maxit = 15$ iterations.**
```{r cache = FALSE}
#create adjusted imputation method vector
meth <- make.method(popNCR)
meth[c("popular", "popteach", "texp", "extrav")] <- "norm"
meth
#create adjusted predictorMatrix
pred <- make.predictorMatrix(popNCR)
pred[, "class"] <- 0
pred[, "pupil"] <- 0
pred
#generate imputations
imp1 <- mice(popNCR, m = 10, maxit = 15, 
             method = meth, 
             predictorMatrix = pred, 
             print = FALSE)
```

---

## Inspecting convergence

---

**9. Inspect convergence of the `mice` algorithm by studying the traceplots**
```{r}
plot(imp1)
```

Convergence has been reached for some variables, but is not convincing for all. However, we've not done any justice to the structure of the data - we've imputed the data as a flat file without paying any attention to the clustering in the data. The imputation model can be greatly improved, which we have seen before in Exercise 3.

---

## Pooling the analysis

---

**9. Run the following model on the imputed data, the incomplete data and compare the model to the 'true' data**

- For the true data:
```{r, message=FALSE}
library(lme4)
fit.true <- with(popular, lmer(popular ~ 1 + (1|class)))
summary(fit.true)
```

- For the incomplete data:
```{r}
fit.true <- with(popNCR, lmer(popular ~ 1 + (1|class)))
summary(fit.true)
```

- For the imputed data:
```{r}
fit <- with(imp1, lmer(popular ~ 1 + (1|class)))
pool <- pool(fit)
summary(pool)
```

We may obtain the variance components by the `testEstimates()` function from package `mitml`:
```{r}
library(mitml)
testEstimates(as.mitml.result(fit), var.comp = TRUE)$var.comp
```

There are quite some differences between the models. With multiple imputation, we are able to correct the inference quite to some extend (e.g. the fixed effects), but we miss out on correcting the inference on the level of the random effects. Then again, we completely ignored the clustering structure in our data. 

---

# Investigating the imputations

---

## Simple comparison of column means

---

**10. Compare the means of the variables in the first imputed dataset and in the incomplete dataset.**
```{r cache = FALSE}
summary(complete(imp1))
summary(popNCR)
```

---

**11. The missingness in `texp` is MNAR: higher values for `texp` have a larger probability to be missing. Can you see this in the imputed data? Do you think this is a problem?**

Yes, we can see this in the imputed data: teacher experience increases slightly after imputation. However, `texp` is the same for all pupils in a class. But not all pupils have this information recorded (as if some pupils did not remember, or were not present during data collection). This is not a problem, because as long as at least one pupil in each class has teacher experience recorded, we can deductively impute the correct (i.e. true) value for every pupil in the class.

---

## Comparing the intraclass correlations

**12. Compare the ICC’s of the variables in the first imputed dataset with those in the incomplete dataset (use `popular`, `popteach` and `texp`). Make a notation of the ICC’s after imputation.**
```{r cache = FALSE}
data.frame(vars = names(popNCR[c(6, 7, 5)]), 
           observed = c(icc(aov(popular ~ class, popNCR)), 
                        icc(aov(popteach ~ class, popNCR)), 
                        icc(aov(texp ~ class, popNCR))), 
           norm     = c(icc(aov(popular ~ class, complete(imp1))), 
                        icc(aov(popteach ~ class, complete(imp1))), 
                        icc(aov(texp ~ class, complete(imp1)))))
```

---

# Incorporating the fixed effects

---

**13. Now impute the `popNCR` dataset again with `mice` using imputation method `norm` for `popular`, `popteach`, `texp` and `extrav`, but now include `class` as a predictor ***for all variables***. Call the `mids`-object `imp2`.**
```{r cache = FALSE}
pred <- make.predictorMatrix(popNCR)
pred[, "pupil"] <- 0
imp2 <- mice(popNCR, meth = meth, pred = pred, print = FALSE)
```
We exclude `pupil` here to avoid overfitting our data with the pupil identifier. Including pupil would result in a model with zero residual variance. 

---

## Comparing the ICC's

---

**14. Compare the ICC’s of the variables in the first imputed dataset from `imp2` with those of `imp1` and the incomplete dataset (use `popular`, `popteach` and `texp`). Make a notation of the ICC’s after imputation.**
```{r cache = FALSE}
data.frame(vars = names(popNCR[c(6, 7, 5)]), 
           observed  = c(icc(aov(popular ~ class, popNCR)), 
                         icc(aov(popteach ~ class, popNCR)), 
                         icc(aov(texp ~ class, popNCR))), 
           norm      = c(icc(aov(popular ~ class, complete(imp1))), 
                         icc(aov(popteach ~ class, complete(imp1))), 
                         icc(aov(texp ~ class, complete(imp1)))), 
           normclass = c(icc(aov(popular ~ class, complete(imp2))), 
                         icc(aov(popteach ~ class, complete(imp2))), 
                         icc(aov(texp ~ class, complete(imp2)))))
```

By simply forcing the algorithm to use the class variable during estimation we adopt a *fixed effects approach*. This conforms to formulating seperate regression models for each `class` and imputing within classes from these models.

---

## Checking Convergence of the imputations

**15. Inspect the trace lines for the variables `popular`, `texp` and `extrav`.**
```{r cache = FALSE}
plot(imp2, c("popular", "texp", "popteach"))
```

It seems not quite OK. For example, there are clear trends for `popteach` and some of the streams do not even intermingle between iterations 4 and 5 for `texp`. Adding another 20 iterations gives a more convincing graph: 

```{r cache = FALSE}
imp2b <- mice.mids(imp2, maxit = 20, print = FALSE)
plot(imp2b, c("popular", "texp", "popteach"))
```

We can use function `mice.mids()` to continue from the state where the object `imp2` stopped. This is fortunate, as we do not have to re-run the first `maxit = 5` default iterations.

---

## Visually inspecting the completed data

**16. Plot the densities of the observed and imputed data (use `imp2`) with the function `densityplot()`.**

To obtain all densities of the different imputed datasets use
```{r cache = FALSE}
densityplot(imp2)
```

---

# Intermezzo

---

<div style="background-color:lightblue">

It is always wise to ask yourself the following questions when *attacking* a multilevel missing data problem - [see Table 7.1 in Van Buuren (2018)](https://stefvanbuuren.name/fimd/sec-missmult.html):

---

Q1.	**Will the complete-data model include random slopes?:** *Thus far we have not formulated any complicated analysis model, but if our analysis model would include random slopes than we should definitely account for those terms in our imputation model. For now, let's assume that our analysis model remains:*
`lme4::lmer(popular ~ (1 | class))`

---

Q2.	**Will the data contain systematically missing values?** *With systematically missing data there are no observed values in the cluster*

We can simply verify this by asking the length of the complete cases for each cluster
```{r}
aggregate(. ~ class, data = popNCR, FUN=function(x) length(is.na(x)))
```
None of the clusters `class` are completely unobserved.

---

Q3.	**Will the distribution of the residuals be non-normal?**
```{r}
fit <- with(popNCR, lmer(popular ~ (1 | class)))
res <- residuals(fit)
hist(res, main = "Residuals")
```
The residuals seem to be quite normally distributed. 

---

Q4.	**Will the error variance differ over clusters?**
```{r}
dotplot(ranef(fit, condVar = TRUE))
```
The error variance seems quite constant from a quick visual inspection

---

Q5.	**Will there be small clusters?**
```{r}
aggregate(pupil ~ class, data = popNCR, length)
```

We can see that the average cluster size is approximately 20 pupils in every cluster with no cluster having fewer than 16 pupils.

---

Q6.	**Will there be a small number of clusters?** *NO, there are 100 clusters*

---

Q7.	**Will the complete-data model have cross-level interactions?** *No*

---

Q8.	**Will the dataset be very large?** *No, only 2000 cases over 100 clusters*

---

** End of intermezzo**

</div>

---

# Replacing `norm` with `pmm`

---

**17. Impute the `popNCR` data once more where you use predictive mean matching and exclude only `pupil` as a predictor. Name the object `imp4`.**
```{r cache = FALSE}
imp4 <- mice(popNCR[, -1])
```

---

## Visually inspecting the completed data

**18. Plot again the densities of the observed and imputed data with the function `densityplot()`, but now use `imp4`. Is there a difference between the imputations obtained with `pmm` and `norm` and can you explain this?**
```{r cache = FALSE}
densityplot(imp4)
```

Yes, `pmm` samples from the observed values and this clearly shows: imputations follow the shape of the observed data.

---

## Comparing the ICC's

**19. Compare the ICC’s of the variables in the first imputed dataset from `imp4` with those of `imp1`, `imp2` and the incomplete dataset (use `popular`, `popteach` and `texp`).**

See **Exercise 12** for the solution. 

---

**20. Finally, compare the ICC’s of the imputations to the ICC’s in the original data. The original data can be found in dataset `popular`. What do you conclude?**
```{r cache = FALSE}
data.frame(vars      = names(popNCR[c(6, 7, 5)]), 
           observed  = c(icc(aov(popular ~ class, popNCR)), 
                         icc(aov(popteach ~ class, popNCR)), 
                         icc(aov(texp ~ class, popNCR))), 
           norm      = c(icc(aov(popular ~ class, complete(imp1))), 
                         icc(aov(popteach ~ class, complete(imp1))), 
                         icc(aov(texp ~ class, complete(imp1)))), 
           normclass = c(icc(aov(popular ~ class, complete(imp2))), 
                         icc(aov(popteach ~ class, complete(imp2))), 
                         icc(aov(texp ~ class, complete(imp2)))), 
           pmm       = c(icc(aov(popular ~ class, complete(imp4))), 
                         icc(aov(popteach ~ class, complete(imp4))), 
                         icc(aov(texp ~ class, complete(imp4)))), 
           orig      = c(icc(aov(popular ~ as.factor(class), popular)), 
                         icc(aov(popteach ~ as.factor(class), popular)), 
                         icc(aov(texp ~ as.factor(class), popular))))
```

Note: these display only the first imputed data set.

---

# True multilevel imputation

---

Mice includes several imputation methods for imputing multilevel data:

- **2l.norm**: Imputes univariate missing data using a two-level normal model with heterogeneous within group variances
- **2l.pan**: Imputes univariate missing data using a two-level normal model with homogeneous within group variances
- **2l.lmer**: Imputes univariate systematically and sporadically missing data using a two-level normal model using `lme4::lmer()`
- **2l.bin**: Imputes univariate systematically and sporadically missing data using a two-level logistic model using `lme4::glmer()`
- **2lonly.mean**: Imputes the mean of the class within the class
- **2lonly.norm**: Imputes univariate missing data at level 2 using Bayesian linear regression analysis
- **2lonly.pmm**: Imputes univariate missing data at level 2 using predictive mean matching

The latter two methods aggregate level 1 variables at level 2, but in combination with `mice.impute.2l.pan`, allow switching regression imputation between level 1 and level 2 as described in Yucel (2008) or Gelman and Hill (2006, p. 541). For more information on these imputation methods see the help.

There are also many useful imputation routines available from packages `miceadds` and `micemd`. See [Tables 7.2, 7.3 and 7.4 in Van Buuren (2018)](https://stefvanbuuren.name/fimd/sec-multioutcome.html) for a detailed overview of these functions. 

---

## Imputing the data

**21. Impute the variable `popular` by means of `2l.norm`. Use dataset `popNCR2`.**
```{r cache = FALSE}
ini <- mice(popNCR2, maxit = 0)
pred <- ini$pred
pred["popular", ] <- c(0, -2, 2, 2, 2, 0, 2)
```

In the predictor matrix, `-2` denotes the class variable, a value `1` indicates a fixed effect and a value `2` indicates a random effect. However, the currently implemented algorithm does not handle predictors that are specified as fixed effects (type = `1`). When using `mice.impute.2l.norm()`, the current advice is to specify all predictors as random effects (type = `2`).
```{r cache = FALSE}
meth <- ini$meth
meth <- c("", "", "", "", "", "2l.norm", "")
imp5 <- mice(popNCR2, pred = pred, meth=meth, print = FALSE)
```

---

## Inspecting the imputations

**22. Inspect the imputations. Did the algorithm converge?**
```{r cache = FALSE}
densityplot(imp5, ~popular, ylim = c(0, 0.35), xlim = c(-1.5, 10))
densityplot(imp4, ~popular, ylim = c(0, 0.35), xlim = c(-1.5, 10))
```

The imputations generated with `2l.norm` are very similar to the ones obtained by `pmm` with `class` as a fixed effect. If we plot the first imputed dataset from `imp4` and `imp5` against the original (true) data:
```{r cache = FALSE}
plot(density(popular$popular))  #true data 
lines(density(complete(imp5)$popular), col = "red", lwd = 2)  #2l.norm
lines(density(complete(imp4)$popular), col = "green", lwd = 2)  #PMM
```

We can see that the imputations are very similar. When studying the convergence
```{r cache = FALSE}
plot(imp5)
```

we conclude that it may be wise to run additional iterations. Convergence is not apparent from this plot.
```{r cache = FALSE}
imp5.b <- mice.mids(imp5, maxit = 10, print = FALSE)
plot(imp5.b)
```

After running another 10 iterations, convergence is more convincing.

---

## Homogeneous group variances

**23. In the original data, the group variances for `popular` are homogeneous. Use `2l.pan` to impute the variable `popular` in dataset `popNCR2`. Inspect the imputations. Did the algorithm converge?**
```{r cache = FALSE}
ini <- mice(popNCR2, maxit = 0)
pred <- ini$pred
pred["popular", ] <- c(0, -2, 2, 2, 1, 0, 2)
meth <- ini$meth
meth <- c("", "", "", "", "", "2l.pan", "")
imp6 <- mice(popNCR2, pred = pred, meth = meth, print = FALSE)
```

Let us create the densityplot for `imp6`
```{r cache = FALSE}
densityplot(imp6, ~popular, ylim = c(0, 0.35), xlim = c(-1.5, 10))
```

and compare it to the one for `imp4`
```{r cache = FALSE}
densityplot(imp4, ~popular, ylim = c(0, 0.35), xlim = c(-1.5, 10))
```

If we plot the first imputed dataset from both objects against the original (true) density, we obtain the following plot:
```{r cache = FALSE}
plot(density(popular$popular), main = "black = truth | green = PMM | red = 2l.pan")  # 
lines(density(complete(imp6)$popular), col = "red", lwd = 2)  #2l.pan
lines(density(complete(imp4)$popular), col = "green", lwd = 2)  #PMM
```

We can see that the imputations are very similar. When studying the convergence
```{r cache = FALSE}
plot(imp6)
```

we conclude that it may be wise to run additional iterations. Convergence is not apparent from this plot.

```{r cache = FALSE}
imp6.b <- mice.mids(imp5, maxit = 10, print = FALSE)
plot(imp6.b)
```

Again, after running another 10 iterations, convergence is more convincing.

---

# Multiple imputaiton of `popNCR3`

---

## True multilevel imputation
**24. Now inspect dataset `popNCR3` and impute the incomplete variables according to the following imputation methods:**

Variable |	Method
:--------|:--------
extrav |	2l.norm
texp |	2lonly.mean
sex |	2l.bin
popular	| 2l.pan
popteach |	2l.lmer

```{r cache = FALSE}
ini <- mice(popNCR3, maxit = 0)
pred <- ini$pred
pred["extrav", ] <- c(0, -2, 0, 2, 2, 2, 2)  #2l.norm
pred["sex", ] <- c(0, -2, 1, 0, 1, 1, 1)  #2l.bin
pred["texp", ] <- c(0, -2, 1, 1, 0, 1, 1)  #2lonly.mean
pred["popular", ] <- c(0, -2, 2, 2, 1, 0, 2)  #2l.pan
pred["popteach", ] <- c(0, -2, 2, 2, 1, 2, 0)  #2l.lmer
meth <- ini$meth
meth <- c("", "", "2l.norm", "2l.bin", "2lonly.mean", "2l.pan", "2l.lmer")
imp7 <- mice(popNCR3, pred = pred, meth = meth, print = FALSE)
```

```{r}
fit <- with(imp7, lmer(popular ~ 1 + (1|class)))
testEstimates(as.mitml.result(fit), var.comp = TRUE)
```

---

### Evaluating the imputations

**25. Evaluate the imputations by means of convergence, distributions and plausibility.**
```{r cache = FALSE}
densityplot(imp7)
```

```{r cache = FALSE}
stripplot(imp7)
```

Given what we know about the missingness, the imputed densities look very reasonable.
```{r cache = FALSE}
plot(imp7)
```

Convergence has not yet been reached. more iterations are advisable.

---

## `pmm` imputation with `class` as a predictor

**26 . Repeat the same imputations as in the previous step, but now use `pmm` for all non-binary incomplete variables. Include `class` as a predictor. **
```{r cache = FALSE}
pmmdata <- popNCR3
pmmdata$class <- as.factor(popNCR3$class)
imp8 <- mice(pmmdata, m = 5, print = FALSE)
```

A warning is printed that there are 90 logged events. When we inspect 
```{r cache = FALSE}
imp8$loggedEvents
```
We find that `texp` has been excluded as a predictor in 90 instances. This is to be expected because when `class` is added as a factor (categorical variable) to the model, a seperate model will be fitted for each `class`. In each of these models, observed `texp` is a constant and, hence, will automatically be removed by `mice` to avoid estimation problems because of a redundant parameter.

With `pmm`, the imputations are very similar and conform to the shape of the observed data.
```{r cache = FALSE}
densityplot(imp8)
```

```{r cache = FALSE}
stripplot(imp8)
```

When looking at the convergence of `pmm`, more iterations are advisable:
```{r cache = FALSE}
plot(imp8)
```

```{r}
fit <- with(imp8, lmer(popular ~ 1 + (1|class)))
testEstimates(as.mitml.result(fit), var.comp = TRUE)
```

If the primary interest is on the fixed effects, using `pmm` in a fixed effects setup where the class is added as a cluster dummy variable, may be an easily implementable alternative to running a full multilevel model: see e.g. [Vink, Lazendic & Van Buuren, 2015](https://dspace.library.uu.nl/bitstream/handle/1874/325909/07_Vink.pdf?sequence=1&isAllowed=y). In our example, when compared to the true data set, parameters for the model `popular ~ 1 + (1|class)` are unbiased. However, some authors have reported that the bias in random slopes and variance components can be substantial. If the focus lies on the random components in the model, a parametric full multilevel imputation would be more efficient if it is properly fitted. See [Drechsler (2015)](http://journals.sagepub.com/doi/abs/10.3102/1076998614563393?journalCode=jebb), [Lüdtke, Robitzsch, and Grund (2017)](http://psycnet.apa.org/doiLanding?doi=10.1037%2Fmet0000096) and [Speidel, Drechsler, and Sakshaug (2017)](https://www.research.manchester.ac.uk/portal/en/publications/biases-in-multilevel-analyses-caused-by-cluster-specific-fixed-effects-imputation(0f103b85-f428-4e15-a9a5-0da23f09c93b).html) for more detail on using cluster dummies in multilevel imputation and estimation.

---

# Conclusions

There are ways to ensure that imputations are not just ‘guesses of unobserved values’. Imputations can be checked by using a standard of reasonability - [see Aboyami, Gelman, and Levy (2010)](https://pdfs.semanticscholar.org/29b6/701a9794daf08fc7566b13e03f0381b25115.pdf) for a wonderful introduction to diagnostics for multiple imputations. We are able to check the differences between observed and imputed values, the differences between their distributions as well as the distribution of the completed data as a whole. If we do this, we can see whether imputations make sense in the context of the problem being studied.

---

# References

Abayomi, K. , Gelman, A. and Levy, M. (2008), Diagnostics for multivariate imputations. *Journal of the Royal Statistical Society: Series C (Applied Statistics)*, 57: 273-291. [Article](https://pdfs.semanticscholar.org/29b6/701a9794daf08fc7566b13e03f0381b25115.pdf)

Drechsler, J. (2015). Multiple Imputation of Multilevel Missing Data: Rigor Versus Simplicity. *Journal of Educational and Behavioral Statistics* 40 (1): 69–95. [Article](http://journals.sagepub.com/doi/abs/10.3102/1076998614563393?journalCode=jebb)

Gelman, A., & Hill, J. (2006). *Data analysis using regression and multilevel/hierarchical models*. [Cambridge University Press](http://www.cambridge.org/nl/academic/subjects/statistics-probability/statistical-theory-and-methods/data-analysis-using-regression-and-multilevelhierarchical-models?format=HB&isbn=9780521867061).

Hox, J. J., Moerbeek, M., & van de Schoot, R. (2010). *Multilevel analysis: Techniques and applications*. [Routledge](https://www.routledge.com/products/9781848728462).

Lüdtke, O., Robitzsch, A., & Grund, S. (2017). Multiple imputation of missing data in multilevel designs: A comparison of different strategies. *Psychological Methods*, 22(1), 141-165. [Article](http://psycnet.apa.org/doiLanding?doi=10.1037%2Fmet0000096)

Speidel, M., J. Drechsler, and J. W. Sakshaug. (2017). Biases in Multilevel Analyses Caused by Cluster-Specific Fixed-Effects Imputation. *Behavior Research Methods*, 1–17. [Article](https://www.research.manchester.ac.uk/portal/en/publications/biases-in-multilevel-analyses-caused-by-cluster-specific-fixed-effects-imputation(0f103b85-f428-4e15-a9a5-0da23f09c93b).html)

Van Buuren, S. *Flexible imputation of missing data*. Chapman & Hall/CRC. 2018. [online book](https://stefvanbuuren.name/fimd/)

Vink, G., Lazendic, G., Van Buuren, S. (2015). Psychological Test and Assessment Modeling, volume 57, issue 4, pp. 577 - 594. [Article](https://dspace.library.uu.nl/bitstream/handle/1874/325909/07_Vink.pdf?sequence=1&isAllowed=y)

Yucel, R. M. (2008). Multiple imputation inference for multivariate multilevel continuous data with ignorable non-response. *Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 366(1874), 2389-2403. [Article](http://rsta.royalsocietypublishing.org/content/366/1874/2389)

---

**- End of exercise**

---
