MinimalExample.Rmd
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 and other vignettes.
Here, we use the NHANES data that are part of the JointAI package. For more info on this data, check the help file for the NHANES data, go to the web page of the National Health and Nutrition Examination Survey (NHANES) and check out the vignette Visualizing Incomplete Data, in which the NHANES data is explored.
Fitting a linear regression model with JointAI is straightforward with the function lm_imp()
:
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 many of these parameters are explained in detail.
n.iter
specifies the length of the 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(lm1)
The function traceplot()
produces a plot of the sampled values across iterations per parameter. By default, three^{1} 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()
).
When convergence has been achieved, we can obtain the result of the model from the model summary.
Results from a model fitted with JointAI can be printed using summary()
:
summary(lm1)
#>
#> Linear model fitted with JointAI
#>
#> Call:
#> lm_imp(formula = SBP ~ gender + age + race + WC + alc + educ +
#> albu + bili, data = NHANES, n.iter = 500, progress.bar = "none")
#>
#> Posterior summary:
#> Mean SD 2.5% 97.5% tail-prob. GR-crit
#> (Intercept) 60.637 22.8464 16.936 105.739 0.00800 0.999
#> genderfemale -3.149 2.2634 -7.561 1.466 0.16133 1.000
#> age 0.364 0.0751 0.211 0.519 0.00000 1.017
#> raceOther Hispanic 0.688 5.1184 -9.485 11.165 0.89200 1.002
#> raceNon-Hispanic White -1.512 3.1741 -7.757 4.487 0.64667 1.004
#> raceNon-Hispanic Black 8.906 3.7582 1.562 16.408 0.01867 1.000
#> raceother 3.756 3.5452 -3.315 10.889 0.29600 1.002
#> WC 0.242 0.0808 0.092 0.398 0.00267 1.005
#> alc>=1 7.267 2.3031 2.636 11.726 0.00267 1.001
#> educhigh -3.403 2.1868 -7.864 0.640 0.10400 1.003
#> albu 5.254 4.0617 -2.869 13.053 0.19467 1.000
#> bili -5.589 5.0370 -16.047 3.904 0.28400 1.003
#>
#> Posterior summary of residual std. deviation:
#> Mean SD 2.5% 97.5% GR-crit
#> sigma_SBP 13.2 0.692 12 14.6 1
#>
#>
#> MCMC settings:
#> Iterations = 101:600
#> Sample size per chain = 500
#> Thinning interval = 1
#> Number of chains = 3
#>
#> Number of observations: 186
The output gives the posterior summary, i.e., the summary of the MCMC (Markov Chain Monte Carlo) sample (which consists of all chains combined).
By default, summary()
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.
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)\):
The Gelman-Rubin^{2} criterion, also available via the function GR_crit()
, compares the within and between chain variation. When it is close enough to 1^{3}, 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:
When the traceplot shows that the chains only converged after 1500 iterations, start = 1500
should be specified in summary()
.
The posterior distributions can be visualized using the function densplot()
:
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.
Besides linear regression models, it is also possible to fit
glm_imp()
(follows the specification of glm()
)clm_imp()
(follows the specification of clm()
from package ordinal)lme_imp()
(follows the specification of lme()
from package nlme)glme_imp()
(analogue the specification used in lme_imp()
and glm_imp()
)clmm_imp()
(follows the specification of clmm2()
from package ordinal)survreg_imp()
(follows the specification of survreg()
from package survival)coxph_imp()
(follows the specification of coxph()
from package survival)The number of chains can be changed with the argument n.chain
.↩︎
Gelman, A and Rubin, DB (1992) Inference from iterative simulation using multiple sequences, Statistical Science, 7, 457-511.
Brooks, SP. and Gelman, A. (1998) General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7, 434-455.↩︎
for example < 1.1; but this is not a generally accepted cut-off↩︎