---
title: "After Fitting"
subtitle: "Summary, visualization and evaluation of the results"
date: "2019-08-29"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 4
vignette: >
  %\VignetteIndexEntry{After fitting}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.align = 'center'
)
library(JointAI)
options(width = 100)
```


This vignette describes how the results from a `JointAI` model can be visualized,
summarized and evaluated.
We use the [NHANES](https://nerler.github.io/JointAI/reference/NHANES.html)
data for examples in cross-sectional data and the 
dataset [simLong](https://nerler.github.io/JointAI/reference/simLong.html) for
examples in longitudinal data.
For more info on these datasets, check out the vignette 
[*Visualizing Incomplete Data*](https://nerler.github.io/JointAI/articles/VisualizingIncompleteData.html),
in which the distributions of variables and missing values in both sets is
explored.


The functions described in this section use, by default, the full MCMC sample
and show only the parameters of the analysis model.
A number of arguments are available to select a subset of the MCMC
samples that is used to calculate the summary.
The argument `subset` allows controlling
which part of the MCMC sample is returned and follows the same logic as
the argument [`monitor_params`](https://nerler.github.io/JointAI/articles/SelectingParameters.html) in [`*_imp()`](https://nerler.github.io/JointAI/reference/model_imp.html).
The use of these arguments is further explained [below](#sec:subset).

## Visualizing the posterior sample
The posterior sample can be visualized by two commonly used plots: a
traceplot, showing samples across iterations, or a plot of the empirical
density of the posterior sample.


### Traceplot
A traceplot shows the sampled values per chain and node throughout
iterations. It allows us to visually evaluate convergence and mixing of the chains,
and can be obtained with the function 
[`traceplot()`](https://nerler.github.io/JointAI/reference/traceplot.html):
```{r, message = FALSE, fig.width = 7, fig.height = 4, out.width = '100%'}
mod13a <- lm_imp(SBP ~ gender + WC + alc + creat, data = NHANES, n.iter = 500, seed = 2020)

traceplot(mod13a)
```


When the sampler has converged the chains show one horizontal band, as in the above figure.
Consequently, when traces show a trend convergence has not been
reached and more iterations are necessary (e.g., using [`add_samples()`](https://nerler.github.io/JointAI/reference/add_samples.html)).

Graphical aspects of the traceplot can be controlled by specifying 
standard graphical arguments via the dots argument `"..."`, which are passed to [`matplot()`](https://stat.ethz.ch/R-manual/R-devel/library/graphics/html/matplot.html).
This allows changing color, linetype and -width, limits, etc.
Arguments `nrow` and/or `ncol` can be supplied to set specific numbers of 
rows and columns for the layout of the grid of plots.

With the argument `use_ggplot` it is possible to get a [**ggplot2**](https://CRAN.R-project.org/package=ggplot2)
version of the traceplot. It can be extended using standard **ggplot2** syntax.
```{r ggtrace15a, fig.width = 7, fig.height = 3.5, out.width = '100%'}
library(ggplot2)
traceplot(mod13a, ncol = 3, use_ggplot = TRUE) +
  theme(legend.position = 'bottom') +
  scale_color_brewer(palette = 'Dark2')
```


### Density plot
The posterior distributions can also be visualized using the function 
[`densplot()`](https://nerler.github.io/JointAI/reference/densplot.html),
which plots the empirical density per node per chain, or combining chains
(when `joined = TRUE`).
```{r, echo = F, eval = F, fig.width = 7, fig.height = 3.5, out.width = '100%'}
densplot(mod13a, joined = TRUE, ncol = 3)
```

The argument `vlines` takes a list of lists, containing specifications passed
to `abline`, and allows us to add (vertical) lines to the plot, e.g., marking zero:
```{r, fig.width = 7, fig.height = 3.5, out.width = '100%'}
densplot(mod13a, ncol = 3, col = c("darkred", "darkblue", "darkgreen"),
         vlines = list(list(v = c(rep(0, nrow(summary(mod13a)$stats) - 1), NA),
                            col = grey(0.8))))

```
or marking the posterior mean and 2.5% and 97.5% quantiles:
```{r densplot15a, fig.width = 7, fig.height = 3.5, out.width = '100%'}
densplot(mod13a, ncol = 3,
         vlines = list(list(v = summary(mod13a)$stats[, "Mean"], lty = 1,
                            lwd = 2),
                       list(v = summary(mod13a)$stats[, "2.5%"], lty = 2),
                       list(v = summary(mod13a)$stats[, "97.5%"], lty = 2)
         )
)
```


Like with `traceplot()` it is possible to use the **ggplot2** version of `densplot()`
when setting `use_ggplot = TRUE`. Here, vertical lines can be added as additional
layers:
```{r ggdens15a, fig.width = 7, fig.height = 3.5, out.width = '100%'}
# fit the complete-case version of the model
mod13a_cc <- lm(formula(mod13a), data = NHANES)


# make a dataset containing the quantiles of the posterior sample and
# confidence intervals from the complete case analysis:
quantDF <- rbind(data.frame(variable = rownames(summary(mod13a)$stat),
                            type = '2.5%',
                            model = 'JointAI',
                            value = summary(mod13a)$stat[, c('2.5%')]
                            ),
                 data.frame(variable = rownames(summary(mod13a)$stat),
                            type = '97.5%',
                            model = 'JointAI',
                            value = summary(mod13a)$stat[, c('97.5%')]
                 ),
                 data.frame(variable = names(coef(mod13a_cc)),
                            type = '2.5%',
                            model = 'cc',
                            value = confint(mod13a_cc)[, '2.5 %']
                 ),
                 data.frame(variable = names(coef(mod13a_cc)),
                            type = '97.5%',
                            model = 'cc',
                            value = confint(mod13a_cc)[, '97.5 %']
                 )
)


# ggplot version:
p13a <- densplot(mod13a, ncol = 3, use_ggplot = TRUE, joined = TRUE) +
  theme(legend.position = 'bottom')


# add vertical lines for the:
# - confidence intervals from the complete case analysis
# - quantiles of the posterior distribution
p13a +
  geom_vline(data = quantDF, aes(xintercept = value, color = model),
             lty = 2) +
  scale_color_manual(name = 'CI from model: ', 
                     limits = c('JointAI', 'cc'),
                     values = c('blue', 'red'),
                     labels = c('JointAI', 'compl.case'))
```


## Model Summary
A summary of the posterior distribution estimated in a `JointAI` model can be 
obtained using the function `summary()`.

The posterior summary consists of the mean, standard deviation and 
quantiles (by default the 2.5% and 97.5% quantiles) of the MCMC samples from all
chains combined, as well as the tail probability (see below) and Gelman-Rubin
criterion (see [section below](#sec:grcrit)).

Additionally, some important characteristics of the MCMC samples on which the
summary is based is given. This includes the range and number of iterations
("Sample size per chain"), thinning interval and number of chains.
Furthermore, the number of observations (number of rows in the data) is given.

```{r, message = FALSE}
summary(mod13a)
```

For mixed models, `summary()` also returns the posterior summary of the
random effects covariance matrix `D` and the number of groups:
```{r, message = FALSE}
library(splines)
mod13b <- lme_imp(bmi ~ GESTBIR + ETHN + HEIGHT_M + ns(age, df = 3),
                  random = ~ ns(age, df = 1) | ID,
                  data = subset(simLong, !is.na(bmi)),
                  n.iter = 500, no_model = 'age', seed = 2019)

summary(mod13b)
```


### Tail probability
The tail probability, calculated as
$2\times\min\left\{Pr(\theta > 0), Pr(\theta < 0)\right\},$
where $\theta$ is the parameter of interest,
is a measure of how likely the value 0 is under the 
estimated posterior distribution.
The figure visualizes three examples of posterior distributions
and the corresponding minimum of $Pr(\theta > 0)$ and $Pr(\theta < 0)$
(shaded area):
```{r tailprob, echo = F, fig.width = 7, fig.height = 1.5, out.width = '100%'}
op <- par(mfrow = c(1, 3), mgp = c(1, 0.6, 0), mar = c(2.5, 1, 2, 1))
mus <- c(1, -1.5, -2.5)

for (i in seq_along(mus)) {
  x <- seq(-3.5, 3.5, length = 1000) + mus[i]
  y <- dnorm(x, mean = mus[i])
  
  plot(x,y, type  = 'l', yaxt = 'n', xaxt = 'n',
       xlab = expression(theta), ylab = "", cex.lab = 1.5,
       main = paste0('tail prob. = ', round(2*pnorm(0, abs(mus[i])), 3)))
  
  if (mus[i] > 0) {
    polygon(x = c(x[x < 0], max(x[x < 0])),
            y = c(y[x < 0], min(y)), col = "#e30f41", border = NA)
  } else {
    polygon(x = c(x[x > 0], min(x[x > 0])),
            y = c(y[x > 0], min(y)), col = "#e30f41", border = NA)
  }
  lines(x,y)
  axis(side = 1, at = 0)
  abline(v = 0, lty = 2)
}
par(op)
```


<!-- By default, `summary()` will only print the posterior summary for the main -->
<!-- model parameters of the analysis model and use the whole MCMC sample. -->
<!-- The argument `subset` allows to control -->
<!-- which part of the MCMC sample is returned and follows the same logic as -->
<!-- the argument `monitor_params` in `*_imp()`. -->
<!-- Arguments `start`, `end` and `thin` are available to select a subset of the MCMC -->
<!-- samples that is used to calculate the summary. -->
<!-- The use of these arguments is further explained [below](#sec:subset)). -->


## Evaluation criteria
The convergence of the MCMC chains and precision of the posterior sample can also
be evaluated more formally. Implemented in **JointAI** are 
the Gelman-Rubin criterion for convergence^[Gelman, A. and Rubin, D.B. (1992).
Inference from Iterative Simulation Using Multiple Sequences.
*Statistical Science* **7**(4), 457-472. doi: 10.1214/ss/1177011136.
<br><br>
Brooks, S. P. and Gelman, A. (1998).
General Methods for Monitoring Convergence of Iterative Simulations.
*Journal of Computational and Graphical Statistics*
**7**(4), 434 - 455. doi: 10.1080/10618600.1998.10474787.]
and a comparison of the Monte Carlo Error with the posterior standard deviation.

### Gelman-Rubin criterion for convergence {#sec:grcrit}
The Gelman-Rubin criterion evaluates convergence by comparing
within and between chain variability and, thus, requires at least two MCMC chains
to be calculated.
It is implemented for `JointAI` objects in the function
[`GR_crit()`](https://nerler.github.io/JointAI/reference/GR_crit.html),
which is based on the function `gelman.diag()` from the package
[**coda**](https://CRAN.R-project.org/package=coda).
The upper limit of the confidence interval should not be much larger than 1.
```{r}
GR_crit(mod13a)
```

Besides the arguments `start`, `end`, `thin`, and `subset`, which are explained
[below](#sec:subset), `GR_crit()` also takes the arguments of `gelman.diag()`.


### Monte Carlo Error {#sec:mcerror}
The precision of the MCMC sample can be checked with the function
[`MC_error()`](https://nerler.github.io/JointAI/reference/MC_error.html).
It uses the function `mcmcse.mat()` from the package [**mcmcse**](https://CRAN.R-project.org/package=mcmcse)
to calculate the Monte Carlo error (the error that is made since the sample
is finite) and compares it to the standard deviation of the posterior sample.
A rule of thumb is that the Monte Carlo error should not be more than 5% of the
standard deviation^[Lesaffre, E. M. and A. B. Lawson (2012).
*Bayesian Biostatistics*. John Wiley & Sons. doi: 10.1002/9781119942412.].
Besides the arguments explained [below](#sec:subset), `MC_error()`
takes the arguments of `mcmcse.mat()`.
```{r}
MC_error(mod13a)
```

`MC_error()` returns an object of class `MCElist`, which is a list containing 
matrices with the posterior mean,
estimated Monte Carlo error, posterior standard deviation and the ratio of the
Monte Carlo error to the posterior standard deviation, for the scaled (if they
are part of the `JointAI` object) and unscaled
(transformed back to the scale of the data) posterior samples.
The associated print method prints only the latter.

To facilitate quick evaluation of the Monte Carlo error to posterior standard
deviation ratio, plotting of an object of class `MCElist` using `plot()` shows
this ratio for each (selected) node and automatically adds a vertical line at
the desired cut-off (by default 5%).
```{r MCE15a, fig.width = 8, fig.height = 3, out.width = '100%'}
par(mar = c(3, 5, 0.5, 0.5), mgp = c(2, 0.6, 0), mfrow = c(1, 2))
plot(MC_error(mod13a))  # left panel: all iterations 101:600
plot(MC_error(mod13a, end = 250))  # right panel: iterations 101:250
```


## Subset of output {#sec:subset}
By default, the functions `traceplot()`, `densplot()`, `summary()`, 
`GR_crit()`, `MC_Error()` and `predict()`
use all iterations of the MCMC sample and consider
only the parameters of the analysis model (if they were monitored).
In this section, we describe how the
set of iterations and parameters to display can be changed using the
arguments `subset`, `start`, `end`, `thin` and `exclude_chains`.

### Subset of parameters
As long as the main parameters have been monitored in a `JointAI` object, 
only these parameters are returned in the model summary, plots and criteria shown above.
When the main parameters of the analysis model were not monitored, i.e.,
`monitor_params = c(analysis_main = FALSE)`, and the argument `subset` is not
specified, all parameters that were monitored are displayed.

To display output for nodes other than the main parameters of the analysis model
or for a subset of nodes, the argument `subset` needs to be specified.

#### Examples
To display only the parameters of the imputation models, we
set `subset = c(analysis_main = FALSE, imp_pars = TRUE)` (after re-estimating
the model with the monitor for these parameters switched on):
```{r, message = FALSE}
mod13c <- update(mod13a, monitor_params = c(imp_pars = TRUE))
summary(mod13c, subset = c(analysis_main = FALSE, imp_pars = TRUE))
```

To select only some of the parameters, they can be specified directly by
name via the `other` element of `subset`.
```{r message = FALSE, eval = FALSE}
densplot(mod13c, subset = list(other = c('creat', 'alc>=1')))
```

This also works when a subset of the imputed values should be displayed:
```{r trace15d, message = FALSE}
# re-fit the model and monitor the imputed values
mod13d <- update(mod13a, monitor_params = c(imps = TRUE))

# select all imputed values for 'WC' (3rd column of Xc)
sub3 <- grep('Xc\\[[[:digit:]]+,3\\]', parameters(mod13d), value = TRUE)
sub3

# pass "sub3" to "subset" via "other", for example in a traceplot:
# traceplot(mod13d, subset = list(other = sub3), ncol = 2)
```

When the number of imputed values is large or in order to check convergence
of random effects, it may not be feasible to plot and inspect all traceplots.
In that case a random subset of, for instance the random effects, can be selected
(output not shown):
```{r}
# re-fit the model monitoring the random effects
mod13e <- update(mod13b, monitor_params = c(ranef = TRUE))

# extract random intercepts and random slopes
ri <- grep('^b\\[[[:digit:]]+,1\\]$', colnames(mod13e$MCMC[[1]]), value = T)
rs <- grep('^b\\[[[:digit:]]+,2\\]$', colnames(mod13e$MCMC[[1]]), value = T)

# to plot the chains of 12 randomly selected random intercepts and slopes:
traceplot(mod13e, subset = list(other = sample(ri, size = 12)), ncol = 4)
traceplot(mod13e, subset = list(other = sample(rs, size = 12)), ncol = 4)
```


### Subset of MCMC samples
With the arguments `start`, `end` and `thin` it is possible to
select which iterations from the MCMC sample are included in the summary.
`start` and `end` specify the first and last iterations to be used,
`thin` the thinning interval. Specification of `start`, thus, allows discarding
a "burn-in", i.e., the iterations before the MCMC chain had converged.

If a particular chain does not have converged it can be excluded from the result
summary or plot using the argument `exclude_chains` which takes a numeric
vector identifying chains to be excluded, for example:

```{r}
mod14 <- lm_imp(SBP ~ gender + WC + alc + creat, data = NHANES, n.iter = 100,
              progress.bar = 'none', n.chains = 5)

densplot(mod14, ncol = 3)
densplot(mod14, exclude_chains = c(2,4), ncol = 3)
```

```{r}
traceplot(mod14, thin = 5, ncol = 3)
traceplot(mod14, start = 150, ncol = 3)
traceplot(mod14, end = 120, ncol = 3)
```



## Predicted values
Often, the aim of an analysis is not only to estimate the association between outcome
and covariates but to predict future outcomes or outcomes for new subjects.

The function [`predict()`](https://nerler.github.io/JointAI/reference/predict.JointAI.html)
allows us to obtain predicted values and corresponding
credible intervals from `JointAI` objects. Note that for mixed
models, currently, only marginal prediction but not prediction conditional on
the random effects is implemented.

A dataset containing data which the prediction should be performed is specified
via the argument `newdata`. If no `newdata` is given, the original data from the
JointAI object are used.
The argument `quantiles` allows the specification of the 
quantiles of the posterior sample that are used to obtain the credible interval
(by default the 2.5% and 97.5% quantile).
Arguments `start`, `end`, `thin` and `exclude_chains` control the subset of MCMC samples that is used.

```{r, echo = FALSE}
options(width = 70)
```

```{r}
predict(mod13a, newdata = NHANES[27, ])
```

`predict()` returns a list with elements `dat`, `fit` and `quantiles`,
containing `newdata` with the predicted values and quantiles appended, 
the predicted values and quantiles that form the credible interval.

### Prediction to visualize non-linear effects
Another reason to obtain predicted values is the visualization of non-linear
effects. 
To facilitate the generation of a dataset for such a prediction,
the function [`predDF()`](https://nerler.github.io/JointAI/reference/predDF.html)
can be used. It generates a `data.frame`
that contains a sequence of values through the range of observed values for a
covariate specified by the argument `var`, and the median or reference value for
all other continuous and categorical variables.


```{r, fig.width = 7, fig.height = 4.5, out.width = '90%'}
# create dataset for prediction
newDF <- predDF(mod13b, var = "age")

# obtain predicted values
pred <- predict(mod13b, newdata = newDF)

# plot predicted values and credible interval
matplot(pred$dat$age, pred$dat[, c('fit', '2.5%', '97.5%')],
        lty = c(1,2,2), type = 'l', col = 1,
        xlab = 'age in months', ylab = 'predicted value')

```


## Export of imputed values
Imputed datasets can be extracted from a `JointAI` object (in which a monitor for
the imputed values has been set, i.e., `monitor_params = c(imps = TRUE)`),
with the function [`get_MIdat()`](https://nerler.github.io/JointAI/reference/get_MIdat.html).

A completed dataset is created by taking the imputed values from a randomly
chosen iteration of the MCMC sample, transforming them back to the original scale
if scaling had been performed before the MCMC sampling, 
and filling them into the original incomplete data.

The argument `m` specifies the number of imputed datasets to be created,
`include` controls whether the original data are included in the long format 
`data.frame` (default is `include = TRUE`),
`start` specifies the first iteration that may be used and `minspace` is
the minimum number of iterations between iterations eligible to be selected.

To make the selection of iterations to form the imputed data reproducible, a
seed value can be specified via the argument `seed`.

When `export_to_SPSS = TRUE` the imputed data is exported to SPSS, i.e., 
a `.txt` file containing the data and a `.sps` file containing SPSS syntax to 
convert the data into an SPSS data file (with ending `.sav`) are written.
Arguments `filename` and `resdir` allow the specification of the name of the `.txt` and
`.sps` file and the directory they are written to.

`get_MIdat()` returns a long-format `data.frame` containing the imputed datasets
(and by default the original data) stacked onto each other. The imputation number
is given in the variable `Imputation_`, column `.id` contains a newly created
id variable for each observation in cross-sectional data (multi-level data should
already contain an id variable).

```{r}
impDF <- get_MIdat(mod13d, m = 10, seed = 2018)
```

The function `plot_imp_distr()` allows us to visually compare the distribution
of the observed and imputed values.
```{r, fig.width = 7, fig.height = 2.5}
plot_imp_distr(impDF, nrow = 1)
```

