---
title: "Analyzing Imputed Data with Multilevel Models and merTools"
author: "Jared Knowles"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Imputation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r setup, echo = FALSE, message=FALSE, warning=FALSE, results='hide'}
knitr::opts_chunk$set(
  cache=FALSE,
  comment="#>",
  collapse=TRUE, 
  echo=TRUE, 
  fig.width = 7
)
library(knitr); library(merTools)
amelia_eval <- "Amelia" %in% rownames(installed.packages())
amelia_uneval <- !amelia_eval
```

## Introduction

Multilevel models are valuable in a wide array of problem areas that involve 
non-experimental, or observational data. In many of these cases the data on 
individual observations may be incomplete. In these situations, the analyst may 
turn to one of many methods for filling in missing data depending on the specific 
problem at hand, disciplinary norms, and prior research. 

One of the most common cases is to use multiple imputation. Multiple imputation 
involves fitting a model to the data and estimating the missing values for 
observations. For details on multiple imputation, and a discussion of some of 
the main implementations in R, look at the documentation and vignettes for the 
`mice` and `Amelia` packages. 

The key difficulty multiple imputation creates for users of multilevel models 
is that the result of multiple imputation is K replicated datasets corresponding 
to different estimated values for the missing data in the original dataset. 

For the purposes of this vignette, I will describe how to use one flavor of 
multiple imputation and the function in `merTools` to obtain estimates from a 
multilevel model in the presence of missing and multiply imputed data. 


## Missing Data and its Discontents

To demonstrate this workflow, we will use the `hsb` dataset in the `merTools` 
package which includes data on the math achievement of a wide sample of students 
nested within schools. The data has no missingness, so first we will simulate 
some missing data. 

```{r}
data(hsb)

# Create a function to randomly assign NA values

add_NA <- function(x, prob){
  z <- rbinom(length(x), 1, prob = prob)
  x[z==1] <- NA
  return(x)
}

hsb$minority <- add_NA(hsb$minority, prob = 0.05)
table(is.na(hsb$minority))

hsb$female <- add_NA(hsb$female, prob = 0.05)
table(is.na(hsb$female))

hsb$ses <- add_NA(hsb$ses, prob = 0.05)
table(is.na(hsb$ses))

hsb$size <- add_NA(hsb$size, prob = 0.05)
table(is.na(hsb$size))

```

```{r impute, message=FALSE, eval = amelia_eval}
# Load imputation library
library(Amelia)
# Declare the variables to include in the imputation data
varIndex <- names(hsb)
# Declare ID variables to be excluded from imputation
IDS <- c("schid", "meanses")
# Imputate
impute.out <- amelia(hsb[, varIndex], idvars = IDS, 
                         noms = c("minority", "female"), 
                         m = 5)
summary(impute.out)
```

```{r boot, message=FALSE, eval = amelia_uneval}
# Amelia is not available so let's just boostrap resample our data
impute.out <- vector(mode = "list", 5)

for (i in 1:5) {
  impute.out[[i]] <- hsb[sample(nrow(hsb), nrow(hsb), replace = TRUE), ]
}

# Declare the variables to include in the imputation data
summary(impute.out)
```


## Fitting and Summarizing a Model List

Fitting a model is very similar 

```{r}
fmla <- "mathach ~ minority + female + ses + meanses + (1 + ses|schid)"
mod <- lmer(fmla, data = hsb)
if(amelia_eval) {
  modList <- lmerModList(fmla, data = impute.out$imputations)
} else {
  # Use bootstrapped data instead
  modList <- lmerModList(fmla, data = impute.out)
}

```


The resulting object `modList` is a list of `merMod` objects the same length as the number 
of imputation datasets. This object is assigned the class of `merModList` and 
`merTools` provides some convenience functions for reporting the results of 
this object. 

Using this, we can directly compare the model fit with missing data excluded to 
the aggregate from the imputed models:

```{r}
fixef(mod) # model with dropped missing
fixef(modList)
```

```{r}
VarCorr(mod) # model with dropped missing
VarCorr(modList) # aggregate of imputed models
```

If you want to inspect the individual models, or you do not like taking the 
mean across the imputation replications, you can take the `merModList` apart 
easily:

```{r}
lapply(modList, fixef)
```

And, you can always operate on any single element of the list:

```{r}
fixef(modList[[1]])
fixef(modList[[2]])
```

## Output of a Model List

```{r}
print(modList)
```

```{r}
summary(modList)
```

```{r}
fastdisp(modList)
```

The standard errors reported for the model list include a correction, Rubin's 
correction (see documentation), which adjusts for the within and between imputation 
set variance as well. 

## Specific Model Information Summaries

```{r}
modelRandEffStats(modList)
modelFixedEff(modList)
VarCorr(modList)
```

### Diagnostics of List Components

```{r}
modelInfo(mod)
```


Let's apply this to our model list. 

```{r}
lapply(modList, modelInfo)
```

### Model List Generics

```{r}
summary(modList)
```

```{r}
modelFixedEff(modList)
```

```{r}
ranef(modList)
```


## Cautions and Notes

Often it is desirable to include aggregate values in the level two or level three 
part of the model such as level 1 SES and level 2 mean SES for the group. In 
cases where there is missingness in either the level 1 SES values, or in the 
level 2 mean SES values, caution and careful thought need to be given to how to 
proceed with the imputation routine. 
