---
title: "Multivariate amputation using ampute"
author: "Rianne Schouten [aut, cre], Peter Lugtig [ctb], Gerko Vink [aut]"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Create simulation data}
  %\usepackage[utf8]{inputenc}
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 5, fig.height = 5, 
                      fig.align = "center", dev = "pdf")
```

When new multiple imputation techniques are tested, simulated data sets have to be made incomplete. Often, a univariate amputation procedure is followed, generating missing values one variable at a time. This procedure is repeated for every variable that should have missing values. Van Buuren (2012, pp. 63, 64) explains in detail how this can be done. However, there are some drawbacks of this approach. First, a univariate amputation procedure makes it difficult to relate the missingness on one variable to the missingness on any other variable. Second, both simulated and real data have a multivariate structure. Applying a univariate amputation procedure to multivariate data would not do justice to the complicated nature of data sets. For these reasons, function `ampute` is developed to perform multivariate amputation according the user's specific desires.

```{r, out.width = "600px", echo = FALSE, fig.cap="Figure 1: Step-by-step flowchart of R-function ampute()"}
knitr::include_graphics("Figures/Flowchart.pdf")
```

Function `ampute()` enables the generation of missing values based on different underlying missingness mechanisms. To this end, `ampute()` has multiple arguments and the most important ones are displayed in Figure 1. Figure 1 also gives an overview of the order in which decisions need to be taken if you would like to generate a missingness mechanism that cannot easily be traced back. This is especially useful for the evaluation and comparison of missing data methods. Nevertheless, the default settings create genuine multivariate missingness as well. 

The underlying method of `ampute()` is based on Brand's (1999) multivariate amputation procedure. An adaptation to this method is made to make continuous amputation possible. The possibility to create MNAR missingness is an extra feature as well. 

In this vignette, Figure 1 serves as a guideline to explain how `ampute`'s arguments can be used to generate multivariate missingness. We will discuss the possiblities by examining each step of Figure 1. Concurrently, we will explain the underlying amputation method since this might give the user the opportunity to not only generate missing values but also use the method's features to generate complicated missingness mechanisms. Our goal is to make the use of `ampute()` easy and straightforward. 

### 1. Data and quick amputation

First, we need a complete data set. In simulation settings, multivariate data can be generated by using `mvrnorm()` from the package `MASS`. Be aware that the covariance matrix should be semi definite. 

```{r}
require("mice")
set.seed(2016)
testdata <- MASS::mvrnorm(n = 500, mu = c(10, 5, 0), 
                    Sigma = matrix(data = c(1.0, 0.2, 0.2, 0.2, 1.0, 0.2, 
                                            0.2, 0.2, 1.0), nrow = 3, byrow = T))
testdata <- as.data.frame(testdata)
summary(testdata)
```

The function `ampute()` immediately works when the data are entered into the function. Storing the result allows you to work with the amputed data. 

```{r}
result <- ampute(testdata)
result
```

#### Class

The return object is of class `mads` (Multivariate Amputed Data Set) and contains all the useful information about the amputation. In total, the object contains the following:

```{r}
names(result)
```

#### Inspect amputed data

The amputed data is stored under `amp`. To see whether the amputation has gone according plan, a quick investigation can be done by using function `md.pattern()`.

```{r}
md.pattern(result$amp)
```

The rows of the table show the different missing data patterns with the number of cases accordingly. The first row always refers to the complete cases. The last column contains the number of variables with missings in that specific pattern. Consequently, each column total describes the number of cells with missing values for that variable. A more thorough explanation of `md.pattern()` can be found in its help file (`?md.pattern`). Note that because `md.pattern()` sorts the columns in increasing amounts of missing information, the order of the variables is different from the order in the data. 

### 2. Proportion of missingness

The proportion of missingness is specified under:

```{r}
result$prop
```

In the default setting, this means that 50% of the cases will have missing values. It is easy to change this proportion by using the argument `prop`. One might also want to specify the percentage of missing cells. For this, the argument `bycases` should be `FALSE`. 

```{r}
result <- ampute(testdata, prop = 0.2, bycases = FALSE)
md.pattern(result$amp)
``` 

An inspection of the result shows that the proportion of missing cells is approximately 20%, as requested (the data set contains 10000 * 3 = 30000 cells, in total, 5964 cells are made missing). `ampute()` automatically calculates the proportion of missing cases that belongs to this setting. 

```{r}
result$prop
```

### 3. Patterns

The basic idea of `ampute()` is the generation of missingness patterns. Each pattern is a combination of missingness on specific variables while other variables remain complete. For example, someone could have forgotten the last page of a questionnaire, resulting in missingness on a specific set of questions. Another missingness pattern could occur when someone is not willing to answer private questions. Or when a participant misses a wave in a longitudinal study. Consequently, each pattern is a specific combination of missing and complete variables. 

The default missingness patterns can by obtained by:  

```{r}
mypatterns <- result$patterns
mypatterns
```

In the `patterns` matrix, each row refers to a missing data pattern and each column to a variable. `0` is used for variables that should have missing values in a particular pattern. `1` is used otherwise. Here, three missing data patterns are specified with missing values on 1 variable only. Note that as a result of this, none of the cases will have missingness on more than one variable. A case either has missingness on V1, V2 or V3 or remains complete.  

Subsequently, the default `patterns` matrix can be changed according your desires. For example, the missingness patterns might be changed into:

```{r}
mypatterns[2, 1] <- 0
mypatterns <- rbind(mypatterns, c(0, 1, 0))
mypatterns
```

By doing this, a missingness pattern is created where cases will have missingness on V1 and V2 but not on V3 (pattern 2). Also, I have added a fourth missing data pattern to create a combinaton of missingness on V1 and V3. 

Now, I can ampute the data again, with the desired `patterns` matrix as its third argument. Note that I have changed the desired proportion of missingness to `bycases = TRUE` (default) and 30%. 

```{r}
result <- ampute(testdata, prop = 0.3, patterns = mypatterns)
md.pattern(result$amp)
```

### 4. Relative frequency

The function `ampute()` works by assigning cases to the different missing data patterns. Each case is assigned to one pattern only, and all cases are divided among the patterns. In other words, every case is \textit{candidate} for a certain missing data pattern. This does not automatically mean that each case will obtain missing values. That decision will be made later on. 

The argument `freq` specifies with which frequency the cases should be divided over the patterns. As a default, equal frequencies are used for the patterns.  

```{r}
result$freq
```

The specifications below are an example of a situation where one would like to impose missing data pattern 1 with a higher frequency than the other missing data patterns. When changing the `freq` argument, one should keep in mind that the sum of the relative frequencies should always be 1 (in order to divide all the cases over the patterns). 

```{r}
result <- ampute(testdata, prop = 0.3, patterns = mypatterns, 
                 freq = c(0.7, 0.1, 0.1, 0.1))
md.pattern(result$amp)
```

An inspection of the missingness patterns using `md.pattern()` shows that indeed pattern 1 has received 7 times as many candidates as pattern 2 to 4. In total, approximately 3000 cases have missing values (which is what was specified under `prop`). 

### 5. Mechanism

At this point, the total proportion of missingness is defined, the missing data patterns are specified as well as the relative frequency with which they should occur. All cases are candidate for a missing data pattern. 

Whether a case will be made missing eventually depends on the missingness mechanism. If each case should have an equal probability of having missing values, the argument `mech` should be changed to `"MCAR"` (Missing Completely At Random). Then, step 6 to 9 do not need to be specified. 

```{r}
result <- ampute(testdata, prop = 0.3, patterns = mypatterns, 
                 freq = c(0.7, 0.1, 0.1, 0.1), mech = "MCAR")
result$mech
```

Two other opties are `"MAR"` (Missing At Random), where the probability to be missing depends on the values of the variables that will not be made missing (the 'observed' variables) or `"MNAR"` (Missing Not At Random), where these probabilities depend on the values of the variables that will be made missing. For a more thorough explanation of these terms, Van Buuren (2012, pp. 6, 7, 31, 32) is useful. 

### 6. Weights

In case of `MAR` (default) or `MNAR` missingness, a weighted sum score will be calculated for each case. This is an important part in the process because the probability that a case will be made missing depends on this score. The way in which the sum scores are used, will be explained in step 8 and step 9. 

The weighted sum scores are built from the variable values and certain pre-specified weights. For each case, the value on a certain variable is multiplied with a weight. This is done for all variables and the resulting values are summed: a weighted sum score. The formula that describes the calculation is: 

\begin{equation}
s_{ik} = \sum\limits_{j=1}^J w_{jk}*c_{ij}
\end{equation}

where $s_{ik}$ is the weighted sum score of a case $i$ in a certain pattern $k$, $w_{jk}$ is the pre-specified weight of a certain variable $j$ in a certain pattern $k$ and $c_{ij}$ is the value of case $i$ on variable $j$. In the example, $j\in\{1, 2, 3\}$ and $k\in\{1, 2, 3, 4\}$ because there are three variables and four missing data patterns. 

The `weights` matrix stores the $w_{jk}$ and is of size #patterns by #variables. The default `weights` matrix for `MAR` missingness is as follows. 

```{r}
result <- ampute(testdata, prop = 0.3, patterns = mypatterns, 
                 freq = c(0.7, 0.1, 0.1, 0.1))
myweights <- result$weights
myweights
```

At first sight, this matrix might be complicated, but if you compare the matrix with the `patterns` matrix that we used, the contents make more sense. 

```{r}
mypatterns
```

As you see, the `weights` and `patterns` matrices are exactly similar. The reason for this is that in case of `MAR` missingness, the variables that will be made missing should not be weighted (`0` in the `weights` matrix). Those are the variables with `0` in the `patterns` matrix. 

The other variables receive a weight of `1`, as a default. Thus, the variables have equal weights. 

In case of `MNAR` missingness, the default `weights` matrix is as follows:

```{r}
result <- ampute(testdata, prop = 0.3, patterns = mypatterns, 
                 freq = c(0.7, 0.1, 0.1, 0.1), mech = "MNAR")
result$weights
```

If this matrix is compared with the `patterns` matrix that was used, it is easy to see that the values are reversed. Since a `MNAR` missingness mechanism means that the the missingness will depend on the variables that will be amputed, those variables are weighted in the `weight` matrix. Those are the variables with a value `0` in the `patterns` matrix. Again, the weighted variables receive equal weights of `1`. The variables that will remain complete are not weighted and receive a `0`.  

Naturally, the idea of the `weights` matrix is to weight variables differently from each other. From now on, we will impose a `MAR` mechanism to the data. We therefore focus at the default `weights` matrix which we stored under `myweights`. 
For instance, we could desire to give the values on variable V2 a higher weight than the values on variable V3. Therefore, for pattern 1, we could change the `weights` matrix into something as: 

```{r}
myweights[1, ] <- c(0, 0.8, 0.4)
```

By choosing the values 0.8 and 0.4, variable V2 is weighted twice as heavy as variable V3. The weight values are relative values, meaning that choosing the values 8 and 4 would have the same effect on the amputation process. In order to clearly see the result of the `weights` setting, one can specify weights for just a few variables and with a relative big distance from each other. 

For pattern 3, variable V1 will be weighted three times as heavy as variable V2. 

```{r}
myweights[3, ] <- c(3, 1, 0)
myweights
```

We will now apply these settings and inspect the results in two ways: boxplots and scatterplots.   

```{r}
result <- ampute(testdata, prop = 0.3, patterns = mypatterns, 
                 freq = c(0.7, 0.1, 0.1, 0.1), weights = myweights)
```

#### Boxplots

First, there is the function `bwplot()` that can be imposed on the `mads` object immediately. `bwplot()` plots the distributions of the amputed and non-amputed data for several variables. This is useful since these boxplots show the relation between the missingness and the variables values. In other words, by examining the boxplots one can see for which values of a certain variable the data will be amputed. Note that not necessarily the values of the variables themselves will be amputed. The boxplots merely show the relation between the amputations and the variables. 

In the function `bwplot()`, the argument `which.pat` can be used to define the patterns you are interested in (default: all patterns). The argument `yvar` should contain the variables names (default: all variables). Besides, the function returns the mean, variance and n of the amputed and non-amputed data for each variable and each pattern requested. In the column `Amp`, a `1` refers to the amputed data and `0` to the non-amputed data. If the descriptives are not required, the argument `descriptives` can be set to `FALSE`. 

```{r}
lattice::bwplot(result, which.pat = c(1, 3), descriptives = TRUE)
```

The medians and boundaries of the boxes show that in pattern 1, the amputed data is shifted to the right with respect to the non-amputed data. For variable V2, this effect is the largest, due to the weight value that was specified. For V1, there is a very small difference between the boxplots of the amputed and non-amputed data. This makes sense, because variable V1 was amputed in the first pattern and therefore set to `0` in the `weights` matrix. The small difference that is visible is due to the positive correlation between V1 on the one side and V2 and V3 on the other side. These correlations were created during the simulation of the data. 

If desired, one could use the function `tsum.test()` from package `BSDA` to perform a t-test on the amputed and non-amputed data. The data returned in the descriptives table can be used for that. For example, to know whether the mean difference between the amputed and non-amputed data for variable V2 in pattern 1 is significant, one could run:

```{r, include = FALSE}
require(BSDA)
```

```{r}
BSDA::tsum.test(mean.x = 0.52879, mean.y = -0.22978,
                s.x = sqrt(0.85752), s.y = sqrt(0.89401),
                n.x = 2099, n.y = 4896)
```

As is visible, there is a significant difference between the amputed and non-amputed data of variable V2 in pattern 1. 

For pattern 3, the difference between the distributions of the amputed and non-amputed data is largest for variable V1, as can be expected due to the weight values in pattern 3. 

#### Scatterplots

Scatterplots might also help to investigate the effect of the specifications. We can directly impose the function `xyplot()` on the `mads` object. The function contains arguments comparable to `bwplot()`. 

For example, we can investigate the weighted sum scores of pattern 1 as follows:

```{r}
xyplot(result, which.pat = 1, colors = mdc(1:2))
```

The scatterplots show that there is a very small relation between V1 and the weighted sum scores. Furthermore, the relation between V2 and the weighted sum scores is very strong, meaning that a case's value on V2 is very important in the generation of the weighted sum score. Actually, this is what causes the differences between the amputed and non-amputed data in the boxplots above. For V3 and the weighted sum scores, the relation is a bit weaker than for V2 but more present than for V1. 

### 7. Continuous

As a default, `ampute()` creates continuous missingness. This means that logit probability functions are used to define a candidate's probability of having missing values. The type of continuous amputation can be decided under `type`. A more thorough explanation of how the logit functions work will be given in part 8 of this vignette.

Instead of using continuous formulas to specify the missingness probabilities, these probabilities can also be defined by hand. Thereto, the argument `cont` should be set to `FALSE`. The `odds` argument can be used to define the probabilities (see part 9). 

### 8. Type

```{r, fig.cap = "Figure 2: Adaptations of continuous logit functions", echo = FALSE}
logistic <- function(x){
  exp(x)/(1+exp(x))
} 
x1 <- seq(-3, 3, 0.01)
y1 <- logistic(-mean(x1)+x1[])
data2 <- as.data.frame(matrix(c(x1, y1), 601, 2))
names(data2) <- c("X1", "Y")
data2[602:1202, "X1"] <- x1
data2[602:1202, "Y"] <- logistic(-abs(x1[]-mean(x1))+0.75)
data2[1203:1803, "X1"] <- x1
data2[1203:1803, "Y"] <- logistic(abs(x1[]-mean(x1))-0.75)
data2[1804:2404, "X1"] <- x1
data2[1804:2404, "Y"] <- logistic(mean(x1)-x1[])
data2["Type"] <- c(rep("RIGHT", 601), rep("MID", 601), rep("TAIL", 601), rep("LEFT", 601))
plot1 <- xyplot(Y ~ X1, data2, group = Type, t = 'l', 
                lty=c(1,2,3,4), 
                cex = 0.5, col = rep("black", 4),
                xlab = "Standardized weighted sum scores", 
                ylab = "Probability",
                key=list(space="top", columns=4, lines = list(lty = c(1, 2, 3, 4)), text = list(c("LEFT", "MID", "RIGHT", "TAIL"), cex = 0.5),
                colors = mdc(1:2)))
plot1
```

The logit functions are more thoroughly explained by Van Buuren (2012, pp. 63, 64), but a quick overview will be given here. Four missingness functions are known: RIGHT, MID, TAIL and LEFT missingness. Figure 2 shows the course of the probabilities for standardized values. 

In `ampute()`, the logit functions will be applied to the weighted sum scores. Consequently, in the situation of RIGHT missingness, cases with high weighted sum scores will have a higher probability to have missing values, compared to cases with low weighted sum scores. For MID missingness, the higher probabilities are given to the cases with weighted sum scores around the average.

For each pattern, a different missingness type can be chosen. In our example, we have four patterns, so four type specifications are required. It is advised to inspect the result with `bwplot()` (below this is done for pattern 2), although the scatterplots give insight as well (as an example, a plot for pattern 4 is shown). 

```{r}
result <- ampute(testdata, prop = 0.3, patterns = mypatterns, 
                 freq = c(0.7, 0.1, 0.1, 0.1), weights = myweights,
                 type = c("RIGHT", "TAIL", "MID", "LEFT"))
```

```{r echo = FALSE}
lattice::bwplot(result, which.pat = 2, descriptives = FALSE)
```

From the boxplots of pattern 2, it becomes visible that the interquartile range (IQR) is much larger for the amputed V3 values compared to the non-amputed data. This is due to the fact that in pattern 2 only V3 defines the missingness. Besides, we requested a TAIL missingness type, which means that all the cases with values at the tails of the distribution of the weighted sum scores (based on merely V3), will be made missing.

```{r echo = FALSE}
xyplot(result, which.pat = 4, colors = mdc(1:2))
```

First, notice that there are much fewer dots in these scatterplots compared to the scatterplots we saw earlier. This is due to the `freq` setting: we specified that only 10 percent of the cases with missing values should have missingness pattern 4. Second, the scatterplots show that all the amputed data are at the left hand side of the weighted sum scores due to the `"LEFT"` setting in the `type` argument. Third, these figures show that there is a perfect relation between variable V2 and the weighted sum scores. Clearly, pattern 4 depends on variable V2 only, which we can remember from the `weights` matrix we used.   

```{r}
result$weights
```

### 9. Odds

If the missingness probabilities should not depend on a continuous probability function, you will have to define the probabilities to be missing yourself. This method was first described by Brand (1999) and consists of two steps. First, one has to decide in how many groups the weighted sum scores should be divided. Second, for each group, an odds value defines the relative probability of having missing values. In other words, each candidate case will be assigned to a group based on his weighted sum score. 

Let us first have a look at the working of the odds values. The default `odds` matrix is as follows:  

```{r}
myodds <- result$odds
myodds
```

This matrix means that for each pattern, the candidates will be divided into four groups. This division occurs based on quantiles, in order to obtain equally sized groups. The values `c(1, 2, 3, 4)` mean that a case with a weighted sum score in the highest quantile will have a probability of having missing values that is four times higher than a candidate with a weighted sum score in the lowest quantile. In Figure 3 the different probabilities that belong to this setting are shown for 100 candidates of pattern 1. 

```{r, include = FALSE}
len <- length(result$scores[[1]])
R <- sample(x = c(1, 0), size = len, 
            prob = c(100 / len, (len - 100) / len), replace = TRUE)
data3 <- matrix(NA, length(R[R == 1]), 3)
data3[, 1] <- result$scores[[1]][R == 1]
ng <- length(result$odds[1, ][!is.na(result$odds[1, ])])
quantiles <- quantile(data3[, 1], probs = seq(0, 1, by = 1 / ng))
Q <- rep(NA, length(data3[, 1]))
for (k in 1:ng) {
  Q <- replace(Q, data3[, 1] >= quantiles[k] 
               & data3[, 1] <= quantiles[k + 1], k)
}
data3[, 2] <- Q
for (l in 1:ng) {
data3[data3[, 2] == l, 3] <- (ng * result$prop * result$odds[1, l]) / 
  sum(result$odds[1, ], na.rm = TRUE)
}
data3 <- as.data.frame(data3)
names(data3) <- c("scores", "group", "prob")
plot2 <- xyplot(prob ~ scores, data = data3, groups = group,
                ylab = "Probability", xlab = "Standardized weighted sum scores",
                pch=c(1,2,3,4), 
                cex = 0.5, col = rep("black", 4),
                key=list(space="top", columns=4, title="", 
                         cex = 1, points = list(pch = c(1, 2, 3, 4)), 
                         text = list(c("Group1", "Group2", "Group3", "Group4"), 
                                     cex = 0.5)))
dat <- result$data[result$cand == 1, "V2"]
data3["V2"] <-  dat[R == 1]
dat <- result$data[result$cand == 1, "V3"]
data3["V3"] <-  dat[R == 1]
plot3 <- xyplot(V3 ~ V2, data = data3, groups = group,
                ylab = "Variable V3", xlab = "Variable V2", 
                pch=c(1,2,3,4), 
                cex = 0.8, col = rep("black", 4),
                key=list(space="top", columns=4, title="", 
                         cex = 1, points = list(pch = c(1, 2, 3, 4)), 
                         text = list(c("Group1", "Group2", "Group3", "Group4"), 
                                     cex = 0.5)))
```

```{r fig.cap = "Figure 3: Probabilities to have missing values for the four groups in pattern 1 and the relation between the groups and the weighted sum scores", echo = FALSE}
plot2
```

As can be seen, there are indeed four groups in pattern 1. The groups have an approximately equal size, with each a certain probability to obtain missing values. The probability of group 4 is indeed four times as large as the probability of group 1.

```{r fig.cap = "Figure 4: Division of odds groups over variables $V_2$ and $V_3$", echo = FALSE}
plot3
```

The relation between the groups and the variable values is shown in Figure 4. Because the relationship between variable V2 and the weighted sum scores is high (due to the `weights` setting), the groups can be distinguished very well. Besides, for higher values of V2, the weighted sum scores are higher. These are the cases that are placed in group 4. Therefore, they are at the right hand side of the V2 scale. For variable V3, the relation between the values and the group allocation is small. This again is due to the `weights` setting. Still, because of the odds values, group 4 is much more to the right of the V3 scale than group 1, 2 and 3. 

Let us now dig deeper into the contents of the `odds` matrix. The #rows of this matrix is equal to #patterns. The #columns can be defined by the user and depends on the desired amputation procedure. 

The amount of values in each row defines the number of groups that will be created from the weighted sum scores in that specific pattern. This number can be different for the different patterns. The values themselves define the odds probabilities of having missing values. Note that the values are relative values. For instance, a setting of `c(2, 6)` is similar to `c(3, 9)`.

The cells in the `odds` matrix that are not used, should be filled with NAs. Let us define the matrix as follows. 

```{r}
myodds[3, ] <- c(1, 0, 0, 1)
myodds[4, ] <- c(1, 1, 2, 2)
myodds <- cbind(myodds, matrix(c(NA, NA, NA, 1, NA, NA, NA, 1), nrow = 4, byrow = F))
myodds
```

We keep the default setting of the first two patterns. Then, for pattern 3, the weighted sum scores will be divided into four groups. The odds values mean that candidates with low weighted sum scores will have a probability to have missing values that is equal to the probability of candidates with high weighted sum scores. However, candidates with weighted sum scores around average will not be made missing. Because pattern 3 depends on variable V1 with a weight of 3 and on variable V2 with a weight of 1, the effect will be most visible for variable V1. 

The weighted sum scores of the fourth pattern will be divided into six groups. All candidates will have a probability of having missing values, but this chance is larger for candidates with weighted sum scores around average. 

```{r}
result <- ampute(testdata, prop = 0.3, patterns = mypatterns, 
                 freq = c(0.7, 0.1, 0.1, 0.1), weights = myweights,
                 cont = FALSE, odds = myodds)
```

```{r}
bwplot(result, which.pat = c(3, 4), descriptives = FALSE)
```

In the boxplots of pattern 3, it is visible that the IQR of the amputed data is larger than the one of the non-amputed data. This is especially the case for variable V1, a bit less for variable V2 but almost not present for variable V3.

In pattern 4, the effect of the specifications is only visible for variable V2, because the other variables will be made missing. In contrast to pattern 3, the amputation is performed in the center of the weighted sum scores, resulting to a MID-like missingness pattern. 

### 10. Combining missingness mechanisms

Function `ampute()` imposes a certain missingness mechanism on all the patterns. If `mech == "MAR"`, for instance, all patterns specified under `patterns` are created with a `MAR` mechanism. Besides, all these patterns should be made missing based on either continuous (`cont == TRUE`) or discrete (`cont == FALSE`) probability distributions. 

In reality, you might want to generate multiple kinds of missingness mechanisms within the same data set. In order to do this, I advise to ampute one data set multiple times according the desired specifications. Thereafter, for each row of the original data set, a row can be sampled from the multiple amputed data sets. In case of three missingness mechanisms, 1/3 of the data set would be the result of amputation 1, 1/3 from amputation 2 and the last third from amputation 3. In total, all three amputation mechanisms will be generated in the data set. 

Below is an example of how this can be done. 

```{r}
ampdata1 <- ampute(testdata, patterns = c(0, 1, 1), prop = 0.2, mech = "MAR")$amp
ampdata2 <- ampute(testdata, patterns = c(1, 0, 1), prop = 0.5, mech = "MNAR")$amp
ampdata3 <- ampute(testdata, patterns = c(1, 1, 0), prop = 0.8, mech = "MCAR")$amp

indices <- sample(x = c(1, 2, 3), size = nrow(testdata), replace = TRUE, 
                  prob = c(1/3, 1/3, 1/3))

ampdata <- matrix(NA, nrow = nrow(testdata), ncol = ncol(testdata))
ampdata[indices == 1, ] <- as.matrix(ampdata1[indices == 1, ])
ampdata[indices == 2, ] <- as.matrix(ampdata2[indices == 2, ])
ampdata[indices == 3, ] <- as.matrix(ampdata3[indices == 3, ])

md.pattern(ampdata)
```

The function `md.pattern()` does not show the different mechanisms (`bwplot()` or `xyplot()` are more useful for this), but the proportions show that the code above works. 

20% of cases is amputed according a MAR mechanism, 50% according a MNAR mechanism and 80 percent according a MCAR mechanism. For each mechanism, merely one missingness pattern is created to clearly see the difference between the three amputation rounds. These rounds are imposed on the data with equal probabilities (1/3 each). Consequently, 1/3 of 20% of the cases has a MAR missingness mechanism, 1/3 of 50% of the cases has a MNAR missingness mechanism and 1/3 of 80% of the cases has a MCAR missingness mechanism. In total, this results in 50% of the cases having missing values. 

### 11. Other specifications

#### Argument `run`

For large data sets or slow computers, you might not want to perform the amputation right away. When the argument `run` is set to `FALSE`, all results will be stored in the `mads` object except for the amputed data set. Subsequently, the default settings for the `patterns`, `weights` or `odds` argument can be changed easily and entered into a new run thereafter (with `run == TRUE`).

```{r}
emptyresult <- ampute(testdata, run = FALSE)
emptyresult$amp
```

#### Other `mads` contents

The return object from `ampute()` is of class `mads`. `mads` contains the amputed data set, the function specifications and some extra objects that might be useful. 

The object `cand` is a vector that contains for every case the missing data pattern it was candidate for. 

```{r}
result$cand[1:30]
```

The object `scores` is a list with, for each pattern, the weighted sum scores of the candidates. 

```{r}
result$scores[[1]][1:10]
```

Furthermore, the object `data` contains the original data. 

```{r}
head(result$data)
```

## Go ahead and `ampute`!

### References

Brand, J.P.L. (1999). Development, implementation and 
evaluation of multiple imputation strategies for the statistical analysis of incomplete data sets (pp. 110-113). Dissertation. Rotterdam: Erasmus University.

Van Buuren, S. (2012). Flexible imputation of missing data.
Boca Raton, FL.: Chapman & Hall/CRC Press.

