Calculate and plot the Monte Carlo error

Calculate, print and plot the Monte Carlo error of the samples from a 'JointAI' model, combining the samples from all MCMC chains.

MC_error(x, subset = NULL, exclude_chains = NULL, start = NULL,
  end = NULL, thin = NULL, digits = 2, warn = TRUE, mess = TRUE, ...)

# S3 method for MCElist
plot(x, data_scale = TRUE, plotpars = NULL,
  ablinepars = list(v = 0.05), minlength = 20, ...)

Arguments

x	object inheriting from class 'JointAI'
subset	subset of parameters/variables/nodes (columns in the MCMC sample). Follows the same principle as the argument `monitor_params` in `*_imp`.
exclude_chains	optional vector of the index numbers of chains that should be excluded
start	the first iteration of interest (see `window.mcmc`)
end	the last iteration of interest (see `window.mcmc`)
thin	thinning interval (integer; see `window.mcmc`). For example, `thin = 1` (default) will keep the MCMC samples from all iterations; `thin = 5` would only keep every 5th iteration.
digits	number of digits for the printed output
warn	logical; should warnings be given? Default is `TRUE`.
mess	logical; should messages be given? Default is `TRUE`.
...	Arguments passed on to `mcmcse::mcse.mat` `size` represents the batch size in ``bm'' and the truncation point in ``bartlett'' and ``tukey''. Default is `NULL` which implies that an optimal batch size is calculated using the `batchSize()` function. Can take character values of ``sqroot'' and ``cuberoot'' or any numeric value between 1 and n/2. ``sqroot'' means size is floor(n^(1/2)) and ``cuberoot'' means size is floor(n^(1/3)). `g` a function such that \(E(g(x))\) is the quantity of interest. The default is `NULL`, which causes the identity function to be used. `method` any of ``bm'', ``obm'', ``bartlett'', ``tukey''. ``bm'' represents batch means estimator, ``obm'' represents overlapping batch means estimator with, ``bartlett'' and ``tukey'' represents the modified-Bartlett window and the Tukey-Hanning windows for spectral variance estimators. `r` the lugsail parameter that converts a lag window into its lugsail equivalent. Larger values of ``r'' will typically imply less underestimation of ``cov'', but higher variability of the estimator. Default is ``r = 3'' and ``r = 1,2'' are good choices. ``r > 5'' is not recommended. Non-integer values are ok.
data_scale	logical; show the Monte Carlo error of the sample transformed back to the scale of the data (`TRUE`) or on the sampling scale (this requires the argument `keep_scaled_mcmc = TRUE` to be set when fitting the model)
plotpars	optional; list of parameters passed to `plot()`
ablinepars	optional; list of parameters passed to `abline()`
minlength	number of characters the variable names are abbreviated to

Value

An object of class MCElist with elements unscaled, scaled and digits. The first two are matrices with columns est (posterior mean), MCSE (Monte Carlo error), SD (posterior standard deviation) and MCSE/SD (Monte Carlo error divided by post. standard deviation.)

Methods (by generic)

plot: plot Monte Carlo error

Note

Lesaffre & Lawson (2012; p. 195) suggest the Monte Carlo error of a parameter should not be more than 5% of the posterior standard deviation of this parameter (i.e., \(MCSE/SD \le 0.05\)).

Long variable names:
The default plot margins may not be wide enough when variable names are longer than a few characters. The plot margin can be adjusted (globally) using the argument "mar" in par.

References

Lesaffre, E., & Lawson, A. B. (2012). Bayesian Biostatistics. John Wiley & Sons.

Examples

mod <- lm_imp(y ~ C1 + C2 + M2, data = wideDF, n.iter = 100)

MC_error(mod)
#>                 est   MCSE    SD MCSE/SD
#> (Intercept)  38.705 1.4511 25.13   0.058
#> C1          -28.963 1.0131 17.55   0.058
#> C2            0.949 0.0372  0.64   0.058
#> M22          -0.421 0.0650  0.72   0.090
#> M23           0.042 0.0391  0.68   0.057
#> M24          -1.173 0.0385  0.67   0.058
#> sigma_y       2.236 0.0091  0.16   0.058

plot(MC_error(mod), ablinepars = list(lty = 2),
     plotpars = list(pch = 19, col = 'blue'))