---
title: "Minimal Example"
date: "2019-08-29"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Minimal Example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

knitr::opts_knit$set(progress = FALSE)
library(JointAI)
```


This vignette gives a first and very brief overview of how the package **JointAI**
can be used. The different settings and options are explained in more depth in
the [help pages](https://nerler.github.io/JointAI/reference/index.html) and 
other [vignettes](https://nerler.github.io/JointAI/articles/index.html).

Here, we use the NHANES data that are part of the **JointAI** package.
For more info on this data, check the
[help file](https://nerler.github.io/JointAI/reference/NHANES.html) for the NHANES data,
go to the [web page of the National Health and Nutrition Examination Survey (NHANES)](https://www.cdc.gov/nchs/nhanes/index.htm) and check out the vignette 
[*Visualizing Incomplete Data*](https://nerler.github.io/JointAI/articles/VisualizingIncompleteData.html),
in which the NHANES data is explored.

##  Fitting a linear regression model
Fitting a linear regression model with **JointAI** is straightforward with the
function `lm_imp()`:
```{r fit_lm1, message = F}
lm1 <- lm_imp(SBP ~ gender + age + race + WC + alc + educ + albu + bili,
              data = NHANES, n.iter = 500, progress.bar = 'none')
```

The specification of `lm_imp()` is similar to the specification of a
linear regression model for complete data using `lm()`.
In this minimal example the only difference is that for `lm_imp()` the number
of iterations `n.iter` has to be specified.
Of course there are many more parameters that can or should
be specified.
In the vignette [*Model Specification*](https://nerler.github.io/JointAI/articles/ModelSpecification.html) 
many of these parameters are explained in detail.

`n.iter` specifies the length of the [Markov Chain](https://en.wikipedia.org/wiki/Markov_chain),
i.e., the number of draws from the posterior distribution of the parameter
or unobserved value.
How many iterations are necessary depends on the data and complexity of the
model and can vary from as few as 100 up to thousands or millions. 

One important criterion is that the Markov chains need to have converged.
This can be evaluated visually with a traceplot.

## Traceplot
```{r results_lm1, fig.width = 8, fig.height = 5, out.width = '100%'}
traceplot(lm1)
```

The function [`traceplot()`](https://nerler.github.io/JointAI/reference/traceplot.html)
produces a plot of the sampled values across 
iterations per parameter. By default, three^[The number of chains can be changed with the argument `n.chain`.] Markov chains are produced for
each parameter and represented by different colors.

When the sampler has converged the chains show a 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 the function
[`add_samples()`](https://nerler.github.io/JointAI/reference/add_samples.html)).

When convergence has been achieved, we can obtain the result of the model
from the model summary.

## Model Summary
Results from a model fitted with **JointAI** can be printed using `summary()`:
```{r, echo = FALSE}
options(width = 80)
```

```{r}
summary(lm1)
```

The output gives the posterior summary, i.e., the summary of the 
MCMC ([Markov Chain Monte Carlo](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo))
sample (which consists of all chains combined).

By default, [`summary()`](https://nerler.github.io/JointAI/reference/summary.JointAI.html)
will only print the posterior summary for the main
model parameters of the analysis model. How to select which parameters are shown
is described in the vignette 
[*Selecting Parameters*](https://nerler.github.io/JointAI/articles/SelectingParameters.html).

The summary consists of the posterior mean, the standard deviation and the 
2.5% and 97.5% quantiles of the MCMC sample, the tail probability and the 
Gelman-Rubin criterion for convergence.


The tail probability is a measure of how likely the value 0 is under the 
estimated posterior distribution, and is calculated as 
$$2\times\min\left\{Pr(\theta > 0), Pr(\theta < 0)\right\}$$
(where $\theta$ is the parameter of interest).

In the following graphics, the shaded areas represent the minimum of 
$Pr(\theta > 0)$ and $Pr(\theta < 0)$:
```{r, echo = F, fig.width = 8, fig.height = 2}
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 = "",
       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 = "#18bc9c", border = NA)
  } else {
    polygon(x = c(x[x > 0], min(x[x > 0])),
            y = c(y[x > 0], min(y)), col = "#18bc9c", border = NA)
  }
  lines(x,y)
  axis(side = 1, at = 0)
  abline(v = 0, lty = 2)
}
par(op)
```

The Gelman-Rubin^[Gelman, A and Rubin, DB (1992) Inference from iterative simulation using multiple sequences, Statistical Science, 7, 457-511.<br>
Brooks, SP. and Gelman, A. (1998) General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7, 434-455.] criterion,
also available via the function [`GR_crit()`](https://nerler.github.io/JointAI/reference/GR_crit.html),
compares the within and between chain variation. When it is close enough to
1^[for example < 1.1; but this is not a generally accepted cut-off], the chains
can be assumed to have converged.

In the model summary, additionally, some important characteristics of the MCMC
samples on which the summary is based, are given.
This includes the range and number of iterations
(= `Sample size per chain`), thinning interval and number of chains.

Furthermore, the number of observations (the sample size of the data) is given.

With the arguments `start`, `end` and `thin` it is possible to
select which iterations from the MCMC sample are included in the summary.

For example:<br>
When the traceplot shows that the chains only converged after
1500 iterations, `start = 1500` should be specified in
[`summary()`](https://nerler.github.io/JointAI/reference/summary.JointAI.html).

## Plot of the posterior distributions
The posterior distributions can be visualized using the function `densplot()`:
```{r, fig.width = 8, fig.height = 5}
densplot(lm1)
```

By default, `densplot()` plots the empirical distribution of each of the chains
separately. When `joined = TRUE` the distributions of the combined chains
are plotted.


## Other types of models
Besides linear regression models, it is also possible to fit

* **generalized linear models**: `glm_imp()` (follows the specification of
[`glm()`](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/glm.html))
* **cumulative logit models** for ordinal outcomes: `clm_imp()` (follows the
specification of `clm()` from package [**ordinal**](https://CRAN.R-project.org/package=ordinal))
* **linear mixed models**: `lme_imp()` (follows the specification of [`lme()`](https://stat.ethz.ch/R-manual/R-devel/library/nlme/html/lme.html)
from package [**nlme**](https://CRAN.R-project.org/package=nlme))
* **generalized linear mixed models**: `glme_imp()` (analogue the specification used
in `lme_imp()` and `glm_imp()`)
* **cumulative logit mixed models** for longitudinal ordinal outcomes: `clmm_imp()` (follows the
specification of `clmm2()` from package [**ordinal**](https://CRAN.R-project.org/package=ordinal))
* **parametric (Weibull) survival models**: `survreg_imp()` (follows the specification
of [`survreg()`](https://stat.ethz.ch/R-manual/R-devel/library/survival/html/survreg.html)
from package [**survival**](https://CRAN.R-project.org/package=survival))
* **Cox proportional hazards survival models**: `coxph_imp()` (follows the specification
of [`coxph()`](https://stat.ethz.ch/R-manual/R-devel/library/survival/html/coxph.html)
from package [**survival**](https://CRAN.R-project.org/package=survival))
