MCMC Settings

In JointAI, models are estimated in the Bayesian framework, using MCMC (Markov Chain Monte Carlo) sampling. The sampling is done by the software JAGS (“Just Another Gibbs Sampler”), which performs Gibbs sampling. JointAI pre-processes the data to get it into a form that can be passed to JAGS and writes the JAGS model. The R package rjags is used as an interface to JAGS.

Here, we describe how to specify the arguments in the main functions that control MCMC related settings. To learn more about how to specify the other arguments in JointAI models or the theoretical background of the statistical approach, check out the vignettes Model Specification and TheoreticalBackground.

In this vignette, we use the NHANES data for examples in cross-sectional data and the dataset simLong for examples in longitudinal data. For more info on these datasets, check out the vignette Visualizing Incomplete Data, in which the distributions of variables and missing values in both sets is explored.

Note:
In several examples we use progress.bar = 'none' which prevents printing of the progress of the MCMC sampling, since this would result in lengthy output in the vignette.

MCMC related parameters in JointAI

The main functions *_imp() have a number of arguments that specify settings for the MCMC sampling:

n.chains: number of MCMC chains
n.adapt: number of iterations in the adaptive phase
n.iter: number of iterations in the sampling phase
thin: thinning degree
monitor_params: parameters/nodes to be monitored
seed: optional seed value for reproducibility
inits: initial values
parallel: whether the MCMC chains should be sampled in parallel
n.cores: how many cores are used for parallel sampling

And some other parameters passed to functions from rjags:

quiet: should printing of information be suppressed?
progress.bar: type of progress bar

Chains and Iterations

In MCMC sampling, values are drawn from a probability distribution. The distribution of the current value is drawn from depends on the previously drawn value (but not on values before that). Values, thus, form a chain. Once the chain has converged, its elements can be seen as a sample from the target posterior distribution.

Number of chains

To evaluate the convergence of MCMC chains it is helpful to create multiple chains that have different starting values. The argument n.chains selects the number of chains (by default, n.chains = 3).

For calculating the model summary, multiple chains are merged.

Adaptive phase

JAGS has an adaptive mode, in which samplers are optimized (for example the step size is adjusted). Samples obtained during the adaptive mode do not form a Markov chain and are discarded. The argument n.adapt controls the length of this adaptive phase.

The default value for n.adapt is 100, which works well in many of the examples considered here. Complex models may require longer adaptive phases. If the adaptive phase is not sufficient for JAGS to optimize the samplers, a warning message will be printed (see example below).

Sampling iterations

n.iter specifies the number of iterations in the sampling phase, i.e., the length of the MCMC chain. How many samples are required to reach convergence and to have sufficient precision depends on the complexity of data and model, and may range from as few as 100 to several million.

Side note: How to evaluate convergence?

Convergence can, for instance, be evaluated visually using a traceplot() or using the Gelman-Rubin diagnostic criterion¹ (implemented in GR_crit(), but also returned with the model summary). The latter compares within and between chain variability and requires the JointAI object to have at least two chains.

Side note: How to check precision?

The precision of the MCMC sample can be checked with the function MC_error(). It calculates 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.²

Thinning

In settings with high autocorrelation, i.e., there are no large jumps in the chain but sampled values are always close to the previous value, it may take many iterations before a sample is created that sufficiently represents the whole range of the posterior distribution.

Processing of such long chains can be slow and take a lot of memory. The parameter thin allows the user to specify if and how much the MCMC chains should be thinned out before storing them. By default thin = 1 is used, which corresponds to keeping all values. A value thin = 10 would result in keeping every 10th value and discarding all other values.

Examples

Default settings

n.adapt = 100 and thin = 1 with 100 sampling iterations

mod1 <- lm_imp(SBP ~ alc, data = NHANES, n.iter = 100, progress.bar = 'none')

The relevant part of the model summary (obtained with summary()) shows that the first 100 iterations (adaptive phase) were discarded, and the 100 iterations that follow form the posterior sample. Thinning was set to 1 and there are 3 chains.

#> [...]
#> MCMC settings:
#> Iterations = 101:200
#> Sample size per chain = 100 
#> Thinning interval = 1 
#> Number of chains = 3

Insufficient adaptation phase

mod2 <- lm_imp(SBP ~ alc, data = NHANES, n.adapt = 10, n.iter = 100, progress.bar = 'none')
#> Warning in rjags::jags.model(file = modelfile, data = data_list, inits = inits, : Adaptation
#> incomplete
#> NOTE: Stopping adaptation

Specifying n.adapt = 10 results in a warning message. The relevant part of the model summary from the resulting model is:

#> [...]
#> MCMC settings:
#> Iterations = 11:110
#> Sample size per chain = 100 
#> Thinning interval = 1 
#> Number of chains = 3

Thinning

mod3 <- lm_imp(SBP ~ alc, data = NHANES, n.iter = 500, thin = 10, progress.bar = 'none')

Here, iterations 110 until 600 are used in the output, but due to thinning interval of ten, the resulting MCMC chains contain only 50 samples instead of 500, that is, the samples from iteration 110, 120, 130, …

#> [...]
#> MCMC settings:
#> Iterations = 110:600
#> Sample size per chain = 50 
#> Thinning interval = 10 
#> Number of chains = 3

Parameters to follow

JAGS only saves the values of MCMC chains for those nodes for which the user has specified that they should be monitored. This is also the case in JointAI.

What are nodes?

Nodes are variables in the Bayesian framework, i.e., everything observed or unobserved. This includes the data and parameters that are estimated, but also missing values in the data or parts of the data that are generally unobserved, such as random effects or latent class indicators.

Specifying which nodes should be monitored

For this purpose, the main analysis functions *_imp have an argument monitor_params.

For details, explanation and examples see the vignette Parameter Selection.

Initial values

Initial values are the starting point for the MCMC sampler. Setting good initial values, i.e., initial values that are likely under the posterior distribution, can speed up convergence. By default, the argument inits = NULL, which means that initial values for the nodes are generated by JAGS. JointAI only sets the initial values for the random number generators to allow reproducibility of the results.

The argument seed allows the specification of a seed value with which the starting values of the random number generator can be reproduced.

User-specified initial values

It is possible to supply initial values directly as

a list or
a function.

Initial values can be specified for every unobserved node, that is, parameters and missing values, but it is possible to only specify initial values for a subset of nodes.

Unless the user-specified initial values contain initial values for the random number generator (named .RNG.name and .RNG.seed), JointAI will add these to the initial values. This is necessary for reproducibility of the results but also when the MCMC chains are sampled in parallel (i.e, when the argument parallel = TRUE).

Initial values in a list of lists

A list containing initial values should have the same length as the number of chains, where each element is a named list of initial values. Initial values should differ between the chains.

For example, to create initial values for the parameter vector beta and the precision parameter tau_SBP for three chains:

init_list <- lapply(1:3, function(i) {
  list(beta = rnorm(4), 
       tau_SBP = rgamma(1, 1, 1))
})

init_list
#> [[1]]
#> [[1]]$beta
#> [1]  0.4194091 -0.2709566 -0.6318248 -0.2284119
#> 
#> [[1]]$tau_SBP
#> [1] 1.680774
#> 
#> 
#> [[2]]
#> [[2]]$beta
#> [1] -0.2662727  0.5281408 -1.7686592 -0.4897908
#> 
#> [[2]]$tau_SBP
#> [1] 1.994782
#> 
#> 
#> [[3]]
#> [[3]]$beta
#> [1] -1.0752636  0.2923947 -0.2066645 -0.3594423
#> 
#> [[3]]$tau_SBP
#> [1] 1.677375

The list is then passed to the argument inits. The amended version of the user provided lists of initial values are stored in the JointAI object:

mod4a <- lm_imp(SBP ~ gender + age + WC, data = NHANES, progress.bar = 'none',
                inits = init_list)

mod4a$mcmc_settings$inits
#> [[1]]
#> [[1]]$beta
#> [1]  0.4194091 -0.2709566 -0.6318248 -0.2284119
#> 
#> [[1]]$tau_SBP
#> [1] 1.680774
#> 
#> [[1]]$.RNG.name
#> [1] "base::Mersenne-Twister"
#> 
#> [[1]]$.RNG.seed
#> [1] 15419
#> 
#> 
#> [[2]]
#> [[2]]$beta
#> [1] -0.2662727  0.5281408 -1.7686592 -0.4897908
#> 
#> [[2]]$tau_SBP
#> [1] 1.994782
#> 
#> [[2]]$.RNG.name
#> [1] "base::Wichmann-Hill"
#> 
#> [[2]]$.RNG.seed
#> [1] 79827
#> 
#> 
#> [[3]]
#> [[3]]$beta
#> [1] -1.0752636  0.2923947 -0.2066645 -0.3594423
#> 
#> [[3]]$tau_SBP
#> [1] 1.677375
#> 
#> [[3]]$.RNG.name
#> [1] "base::Wichmann-Hill"
#> 
#> [[3]]$.RNG.seed
#> [1] 83115

Initial values as a function

Initial values can be specified by a function. The function should either take no arguments or a single argument called chain, and return a named list that supplies values for one chain.

For example, to create initial values for the parameter vectors beta and alpha:

inits_fun <- function() {
  list(beta = rnorm(4),
       alpha = rnorm(3))
}

inits_fun()
#> $beta
#> [1] -0.5874427 -0.5274437  0.5960212 -1.2713124
#> 
#> $alpha
#> [1] -0.2012545  0.4556339 -0.9348269


mod4b <- lm_imp(SBP ~ gender + age + WC, data = NHANES, progress.bar = 'none',
                inits = inits_fun)

mod4b$mcmc_settings$inits
#> [[1]]
#> [[1]]$beta
#> [1] -0.6377097  0.8351830  0.5890226  0.1700938
#> 
#> [[1]]$alpha
#> [1] -0.43717159 -0.08020806 -0.92754276
#> 
#> [[1]]$.RNG.name
#> [1] "base::Wichmann-Hill"
#> 
#> [[1]]$.RNG.seed
#> [1] 31227
#> 
#> 
#> [[2]]
#> [[2]]$beta
#> [1]  1.51886715  0.08573045  2.15279476 -1.36593150
#> 
#> [[2]]$alpha
#> [1]  0.345906225  0.565180699 -0.006097837
#> 
#> [[2]]$.RNG.name
#> [1] "base::Mersenne-Twister"
#> 
#> [[2]]$.RNG.seed
#> [1] 29634
#> 
#> 
#> [[3]]
#> [[3]]$beta
#> [1] -0.4387811  0.6888911 -0.4994192 -1.8686146
#> 
#> [[3]]$alpha
#> [1]  0.74342539 -0.08917493  1.08116368
#> 
#> [[3]]$.RNG.name
#> [1] "base::Mersenne-Twister"
#> 
#> [[3]]$.RNG.seed
#> [1] 53420

When a function is supplied, the function will be evaluated once per chain and the resulting list is stored.

For which nodes can initial values be specified?

Initial values can be specified for all unobserved stochastic nodes, i.e., parameters or unobserved data for which a distribution is specified in the JAGS model.

Initial values have to be supplied in the format the parameter or unobserved value is used in the JAGS model.

To find out which nodes there are in a model, the function coef() from package rjags can be used on a JAGS model object. It returns a list with the current values of the MCMC chains, by default the first chain. Elements of the initial values should have the same structure as the elements in this list.

Example:
We are using a longitudinal model and the simLong data in this example. The output is abbreviated to show the relevant characteristics.

mod4c <- lme_imp(bmi ~ time + HEIGHT_M + hc + SMOKE, random = ~ time | ID,
                 data = simLong, no_model = 'time', progress.bar = 'none')

coef(mod4c$model)

#> $RinvD
#>           [,1]      [,2]
#> [1,] 0.7659268        NA
#> [2,]        NA 0.2308901
#> 
#> $RinvD_hc
#>           [,1]         [,2]
#> [1,] 0.0400501           NA
#> [2,]        NA 0.0001842381
#> 
#> $Xc
#>        [,1]       [,2] [,3] [,4]
#>   [1,]   NA         NA   NA   NA
#>   [2,]   NA         NA   NA   NA
#> 
#> [...]
#> 
#>  [82,]   NA         NA   NA   NA
#>  [83,]   NA  0.4091066   NA   NA
#>  [84,]   NA         NA   NA   NA
#> 
#> [...]
#> 
#> [199,]   NA         NA   NA   NA
#> [200,]   NA         NA   NA   NA
#> 
#> $Xcat
#>        [,1]
#>   [1,]   NA
#>   [2,]   NA
#> 
#> [...]
#> 
#> [139,]   NA
#> [140,]    1
#> [141,]   NA
#> [142,]    1
#> [143,]   NA
#> 
#> [...]
#> 
#> [199,]   NA
#> [200,]   NA
#> 
#> $Xl
#>                 [,1]
#>    [1,]           NA
#>    [2,]           NA
#>    [3,]  0.000000000
#>    [4,]           NA
#>    [5,]           NA
#>    [6,]           NA
#>    [7,] -0.214817175
#> 
#> [...]
#> 
#> [2399,]           NA
#> [2400,]           NA
#> 
#> $alpha
#> [1] -0.023398669 -0.011456948  0.050830346  0.002318776 -0.195963349  0.013870602  0.834277915
#> 
#> $b
#>            [,1]       [,2]
#>   [1,] 15.89715 -1.0462718
#>   [2,] 16.92166 -1.4536040
#> 
#> [...]
#> 
#> [200,] 16.98724 -1.3006194
#> 
#> $b_hc
#>                 [,1]      [,2]
#>   [1,] -4.812325e-01 0.8058138
#>   [2,] -3.945752e-02 0.8195184
#> 
#> [...]
#> 
#> [199,] -1.455908e-01 0.8114976
#> [200,] -1.024603e-01 0.8528284
#> 
#> $beta
#> [1] 16.47083724  0.05746001 -0.21654280 -0.09097774 -1.37860072  1.44401948
#> 
#> $bmi
#>    [1] NA NA  0  0 NA NA NA NA NA NA NA NA  0  0 NA  0 NA  0 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 
#> [...]
#> 
#> [2388] NA NA  0 NA  0 NA NA NA NA NA NA NA NA
#> 
#> $delta_SMOKE
#> [1] -1.185836
#> 
#> $gamma_SMOKE
#> [1] 1.035457       NA
#> 
#> $invD
#>           [,1]      [,2]
#> [1,]  1.563199 -1.320860
#> [2,] -1.320860  8.350886
#> 
#> $invD_hc
#>            [,1]      [,2]
#> [1,]   27.14434 -181.7469
#> [2,] -181.74693 2700.2794
#> 
#> $mu_b
#>        [,1]      [,2]
#>   [1,]   NA -1.378601
#>   [2,]   NA -1.378601
#> 
#> [...]
#> 
#> [199,]   NA -1.378601
#> [200,]   NA -1.378601
#> 
#> $mu_b_hc
#>        [,1]      [,2]
#>   [1,]   NA 0.8342779
#>   [2,]   NA 0.8342779
#> 
#> [...]
#> 
#> [200,]   NA 0.8342779
#> 
#> $tau_HEIGHT_M
#> [1] 1.053073
#> 
#> $tau_bmi
#> [1] 4.897738
#> 
#> $tau_hc
#> [1] 5.353787

RinvD is the scale matrix in the Wishart prior for the inverse of the random effects design matrix D. In the data that is passed JAGS (which is stored in the element data_list in a JointAI object), this matrix is specified as diagonal matrix, with unknown diagonal elements:

mod4c$data_list['RinvD']
#> $RinvD
#>      [,1] [,2]
#> [1,]   NA    0
#> [2,]    0   NA

These diagonal elements are estimated in the model and have a Gamma prior. The corresponding part of the JAGS model is:

#> [...]
#>   # Priors for the covariance of the random effects
#>   for (k in 1:2){
#>     RinvD[k, k] ~ dgamma(shape_diag_RinvD, rate_diag_RinvD)
#>   }
#>   invD[1:2, 1:2] ~ dwish(RinvD[ , ], KinvD)
#>   D[1:2, 1:2] <- inverse(invD[ , ])
#> [...]

The element RinvD in the initial values has to be a 2 $\times$ 2 matrix, with positive values on the diagonal and NA as off-diagonal elements, since these are fixed in the data.

The following element in the output of coef(mod4c$model) is Xc, the fixed effects design matrix containing baseline covariates. The corresponding matrix in the data is

mod4c$data_list$Xc[81:85, ]

#>       (Intercept)  HEIGHT_M SMOKEsmoked until[...] SMOKEcontin[...]
#> 171.1           1 0.2736946                     NA               NA
#> 172.1           1 0.1148171                     NA               NA
#> 173.1           1        NA                     NA               NA
#> 174.1           1 0.5045126                     NA               NA
#> 175.1           1 1.8822249                     NA               NA

Again, the matrix Xc in the initial values has the same dimension as Xc in the data, and values where there are missing values in Xc, e.g., Xc[83, 2] and NA elsewhere. There are no initial values specified for the third and fourth column, since these columns contain the dummy variables corresponding to the categorical variable SMOKE and are calculated from a column in the matrix Xcat, which contains the original version of SMOKE:

head(mod4c$data_list$Xcat)
#>       SMOKE
#> 1.1       1
#> 10.1      1
#> 100.1     3
#> 101.1     1
#> 102.1     1
#> 103.1    NA

The relevant part of the JAGS model is:

#> [...]
#>     # ordinal model for SMOKE
#>     Xcat[i, 1] ~ dcat(p_SMOKE[i, 1:3])
#> 
#> [...]
#> 
#>     Xc[i, 3] <- ifelse(Xcat[i, 1] == 2, 1, 0)
#>     Xc[i, 4] <- ifelse(Xcat[i, 1] == 3, 1, 0)
#> [...]

Where Xcat in the data has missing values (e.g., rows 37 and 39), initial values need to be specified.

Elements that are completely unobserved, like the parameter vectors alpha and beta, the random effects b or scalar parameters delta_SMOKE and gamma_SMOKE are entirely specified in the initial values.

Gelman, A., Meng, X. L., & Stern, H. (1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 733-760.↩︎
See p. 195 of Lesaffre, E., & Lawson, A. B. (2012). Bayesian Biostatistics. John Wiley & Sons.↩︎

MCMC Settings

2019-08-29

Chains and Iterations

Number of chains

Adaptive phase

Sampling iterations

Side note: How to evaluate convergence?

Side note: How to check precision?

Thinning

Examples

Default settings

Insufficient adaptation phase

Thinning

Parameters to follow

What are nodes?

Specifying which nodes should be monitored

Initial values

User-specified initial values

Initial values in a list of lists

Initial values as a function

For which nodes can initial values be specified?

Parallel sampling

Other arguments