Function to help fulfill elements of the Bayesian Analysis Reporting Guidelines.
Source:R/barg.R
barg.RdThe Bayesian Analysis Reporting Guidelines were put forward by Kruschke (https://www.nature.com/articles/s41562-021-01177-7) to aide in reproducibility and documentation of Bayesian statistical analyses that are sometimes unfamiliar to reviewers or scientists. The purpose of this function is to summarize goodness of fit metrics from one or more Bayesian models made by growthSS and fitGrowth. See details for explanations of those metrics and the output.
Arguments
- fit
The brmsfit object or a list of brmsfit objects in the case that you split models to run on subsets of the data for computational simplicity.
- ss
The growthSS output used to specify the model. If fit is a list then this can either be one growthSS list in which case the priors are assumed to be the same for each model or it can be a list of the same length as fit. Note that the only parts of this which are used are the
call$startwhich is expected to be a call,pcvrForm, anddflist elements, so if you have a list of brmsfit objects and no ss object you can specify a stand-in list. This can also be left NULL (the default) and posterior predictive plots and prior predictive plots will not be made.
Value
A named list containing Rhat, ESS, NEFF, and Prior/Posterior Predictive plots. See details for interpretation.
Details
General: This includes chain number, length, and total divergent transitions per model. Divergent transitions are a marker that the MCMC had something go wrong. Conceptually it may be helpful to think about rolling a marble over a 3D curve then having the marble suddenly jolt in an unexpected direction, something happened that suggests a problem/misunderstood surface. In practice you want extremely few (ideally no) divergences. If you do have divergences then consider specifying more control parameters (see brms::brm or examples for fitGrowth). If the problem persists then the model may need to be simplified. For more information on MCMC and divergence see the stan manual (https://mc-stan.org/docs/2_19/reference-manual/divergent-transitions).
ESS: ESS stands for Effective Sample Size and is a goodness of fit metric that approximates the number of independent replicates that would equate to the same amount of information as the (autocorrelated) MCMC iterations. ESS of 1000+ is often considered as a pretty stable value, but more is better. Still, 100 per chain may be plenty depending on your applications and the inference you wish to do. One of the benefits to using lots of chains and/or longer chains is that you will get more complete information and that benefit will be shown by a larger ESS. This is separated into "bulk" and "tail" to represent the middle and tails of the posterior distribution, since those can sometimes have very different sampling behavior. A summary and the total values are returned, with the summary being useful if several models are included in a list for fit argument
Rhat: Rhat is a measure of "chain mixture". It compares the between vs within chain values to assess how well the chains mixed. If chains did not mix well then Rhat will be greater than 1, with 1.05 being a broadly agreed upon cutoff to signify a problem. Running longer chains should result in lower Rhat values. The default in brms is to run 4 chains, partially to ensure that there is a good chance to check that the chains mixed well via Rhat. A summary and the total values are returned, with the summary being useful if several models are included in a list for fit argument
NEFF: NEFF is the NEFF ratio (Effective Sample Size over Total MCMC Sample Size). Values greater than 0.5 are generally considered good, but there is a consensus that lower can be fine down to about 0.1. A summary and the total values are returned, with the summary being useful if several models are included in a list for fit argument
priorPredictive: A plot of data simulated from the prior using plotPrior. This should generate data that is biologically plausible for your situation, but it will probably be much more variable than your data. That is the effect of the mildly informative thick tailed lognormal priors. If you specified non-default style priors then this currently will not work.
posteriorPredictive: A plot of each model's posterior predictive interval over time. This is the same as plots returned from growthPlot and shows 1-99 coming to a mean yellow trend line. These should encompass the overwhelming majority of your data and ideally match the variance pattern that you see in your data. If parts of the predicted interval are biologically impossible (area below 0, percentage about 100 model should be reconsidered.
See also
plotPrior for visual prior predictive checks.
Examples
# \donttest{
simdf <- growthSim("logistic",
n = 20, t = 25,
params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5))
)
ss <- growthSS(
model = "logistic", form = y ~ time | id / group, sigma = "logistic",
df = simdf, start = list(
"A" = 130, "B" = 12, "C" = 3,
"sigmaA" = 20, "sigmaB" = 10, "sigmaC" = 2
), type = "brms"
)
fit_test <- fitGrowth(ss,
iter = 600, cores = 1, chains = 1, backend = "cmdstanr",
sample_prior = "only" # only sampling from prior for speed
)
#> Start sampling
#> Init values were only set for a subset of parameters.
#> Missing init values for the following parameters:
#> Intercept_nu
#>
#> To disable this message use options(cmdstanr_warn_inits = FALSE).
#> Running MCMC with 1 chain...
#>
#> Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
#> Chain 1 Iteration: 100 / 600 [ 16%] (Warmup)
#> Chain 1 Iteration: 200 / 600 [ 33%] (Warmup)
#> Chain 1 Iteration: 300 / 600 [ 50%] (Warmup)
#> Chain 1 Iteration: 301 / 600 [ 50%] (Sampling)
#> Chain 1 Iteration: 400 / 600 [ 66%] (Sampling)
#> Chain 1 Iteration: 500 / 600 [ 83%] (Sampling)
#> Chain 1 Iteration: 600 / 600 [100%] (Sampling)
#> Chain 1 finished in 0.0 seconds.
#> Loading required package: rstan
#> Loading required package: StanHeaders
#>
#> rstan version 2.32.6 (Stan version 2.32.2)
#> For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores()).
#> To avoid recompilation of unchanged Stan programs, we recommend calling
#> rstan_options(auto_write = TRUE)
#> For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
#> change `threads_per_chain` option:
#> rstan_options(threads_per_chain = 1)
barg(fit_test, ss)
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: Coercing LHS to a list
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: Coercing LHS to a list
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> $General
#> chains iter num.divergent model
#> 1 1 600 0 1
#>
#> $Rhat
#> $Rhat$summary
#> b_nu_Intercept b_A_groupa b_A_groupb b_B_groupa b_B_groupb b_C_groupa
#> Min. 1.007256 1.003564 1.01621 1.011313 0.9982719 1.023324
#> 1st Qu. 1.007256 1.003564 1.01621 1.011313 0.9982719 1.023324
#> Median 1.007256 1.003564 1.01621 1.011313 0.9982719 1.023324
#> Mean 1.007256 1.003564 1.01621 1.011313 0.9982719 1.023324
#> 3rd Qu. 1.007256 1.003564 1.01621 1.011313 0.9982719 1.023324
#> Max. 1.007256 1.003564 1.01621 1.011313 0.9982719 1.023324
#> b_C_groupb b_sigmaA_groupa b_sigmaA_groupb b_sigmaB_groupa
#> Min. 1.014474 1.001622 1.01393 1.001079
#> 1st Qu. 1.014474 1.001622 1.01393 1.001079
#> Median 1.014474 1.001622 1.01393 1.001079
#> Mean 1.014474 1.001622 1.01393 1.001079
#> 3rd Qu. 1.014474 1.001622 1.01393 1.001079
#> Max. 1.014474 1.001622 1.01393 1.001079
#> b_sigmaB_groupb b_sigmaC_groupa b_sigmaC_groupb Intercept_nu lprior
#> Min. 1.001481 1.002892 0.9997347 1.007256 0.9984929
#> 1st Qu. 1.001481 1.002892 0.9997347 1.007256 0.9984929
#> Median 1.001481 1.002892 0.9997347 1.007256 0.9984929
#> Mean 1.001481 1.002892 0.9997347 1.007256 0.9984929
#> 3rd Qu. 1.001481 1.002892 0.9997347 1.007256 0.9984929
#> Max. 1.001481 1.002892 0.9997347 1.007256 0.9984929
#> lp__
#> Min. 0.9987048
#> 1st Qu. 0.9987048
#> Median 0.9987048
#> Mean 0.9987048
#> 3rd Qu. 0.9987048
#> Max. 0.9987048
#>
#> $Rhat$complete
#> $Rhat$complete[[1]]
#> [1] 1.007256
#>
#> $Rhat$complete[[2]]
#> [1] 1.003564
#>
#> $Rhat$complete[[3]]
#> [1] 1.01621
#>
#> $Rhat$complete[[4]]
#> [1] 1.011313
#>
#> $Rhat$complete[[5]]
#> [1] 0.9982719
#>
#> $Rhat$complete[[6]]
#> [1] 1.023324
#>
#> $Rhat$complete[[7]]
#> [1] 1.014474
#>
#> $Rhat$complete[[8]]
#> [1] 1.001622
#>
#> $Rhat$complete[[9]]
#> [1] 1.01393
#>
#> $Rhat$complete[[10]]
#> [1] 1.001079
#>
#> $Rhat$complete[[11]]
#> [1] 1.001481
#>
#> $Rhat$complete[[12]]
#> [1] 1.002892
#>
#> $Rhat$complete[[13]]
#> [1] 0.9997347
#>
#> $Rhat$complete[[14]]
#> [1] 1.007256
#>
#> $Rhat$complete[[15]]
#> [1] 0.9984929
#>
#> $Rhat$complete[[16]]
#> [1] 0.9987048
#>
#> $Rhat$complete$model
#> [1] 1
#>
#>
#>
#> $NEFF
#> $NEFF$summary
#> b_nu_Intercept b_A_groupa b_A_groupb b_B_groupa b_B_groupb b_C_groupa
#> Min. 0.8077354 0.5884866 0.7148728 0.5149354 0.6861411 0.418486
#> 1st Qu. 0.8077354 0.5884866 0.7148728 0.5149354 0.6861411 0.418486
#> Median 0.8077354 0.5884866 0.7148728 0.5149354 0.6861411 0.418486
#> Mean 0.8077354 0.5884866 0.7148728 0.5149354 0.6861411 0.418486
#> 3rd Qu. 0.8077354 0.5884866 0.7148728 0.5149354 0.6861411 0.418486
#> Max. 0.8077354 0.5884866 0.7148728 0.5149354 0.6861411 0.418486
#> b_C_groupb b_sigmaA_groupa b_sigmaA_groupb b_sigmaB_groupa
#> Min. 0.8313222 0.8249887 0.2686725 0.9180275
#> 1st Qu. 0.8313222 0.8249887 0.2686725 0.9180275
#> Median 0.8313222 0.8249887 0.2686725 0.9180275
#> Mean 0.8313222 0.8249887 0.2686725 0.9180275
#> 3rd Qu. 0.8313222 0.8249887 0.2686725 0.9180275
#> Max. 0.8313222 0.8249887 0.2686725 0.9180275
#> b_sigmaB_groupb b_sigmaC_groupa b_sigmaC_groupb Intercept_nu lprior
#> Min. 0.5739564 0.645482 0.7389399 0.8077354 0.4609321
#> 1st Qu. 0.5739564 0.645482 0.7389399 0.8077354 0.4609321
#> Median 0.5739564 0.645482 0.7389399 0.8077354 0.4609321
#> Mean 0.5739564 0.645482 0.7389399 0.8077354 0.4609321
#> 3rd Qu. 0.5739564 0.645482 0.7389399 0.8077354 0.4609321
#> Max. 0.5739564 0.645482 0.7389399 0.8077354 0.4609321
#> lp__
#> Min. 0.4192135
#> 1st Qu. 0.4192135
#> Median 0.4192135
#> Mean 0.4192135
#> 3rd Qu. 0.4192135
#> Max. 0.4192135
#>
#> $NEFF$complete
#> $NEFF$complete[[1]]
#> [1] 0.8077354
#>
#> $NEFF$complete[[2]]
#> [1] 0.5884866
#>
#> $NEFF$complete[[3]]
#> [1] 0.7148728
#>
#> $NEFF$complete[[4]]
#> [1] 0.5149354
#>
#> $NEFF$complete[[5]]
#> [1] 0.6861411
#>
#> $NEFF$complete[[6]]
#> [1] 0.418486
#>
#> $NEFF$complete[[7]]
#> [1] 0.8313222
#>
#> $NEFF$complete[[8]]
#> [1] 0.8249887
#>
#> $NEFF$complete[[9]]
#> [1] 0.2686725
#>
#> $NEFF$complete[[10]]
#> [1] 0.9180275
#>
#> $NEFF$complete[[11]]
#> [1] 0.5739564
#>
#> $NEFF$complete[[12]]
#> [1] 0.645482
#>
#> $NEFF$complete[[13]]
#> [1] 0.7389399
#>
#> $NEFF$complete[[14]]
#> [1] 0.8077354
#>
#> $NEFF$complete[[15]]
#> [1] 0.4609321
#>
#> $NEFF$complete[[16]]
#> [1] 0.4192135
#>
#> $NEFF$complete$model
#> [1] 1
#>
#>
#>
#> $ESS
#> $ESS$summary
#> A_groupa A_groupb B_groupa B_groupb C_groupa C_groupb
#> Bulk_ESS.Min. 714.7464 593.0616 621.8739 681.5482 561.5601 743.1364
#> Bulk_ESS.1st Qu. 714.7464 593.0616 621.8739 681.5482 561.5601 743.1364
#> Bulk_ESS.Median 714.7464 593.0616 621.8739 681.5482 561.5601 743.1364
#> Bulk_ESS.Mean 714.7464 593.0616 621.8739 681.5482 561.5601 743.1364
#> Bulk_ESS.3rd Qu. 714.7464 593.0616 621.8739 681.5482 561.5601 743.1364
#> Bulk_ESS.Max. 714.7464 593.0616 621.8739 681.5482 561.5601 743.1364
#> Tail_ESS.Min. 176.5460 214.4618 154.4806 205.8423 125.5458 249.3967
#> Tail_ESS.1st Qu. 176.5460 214.4618 154.4806 205.8423 125.5458 249.3967
#> Tail_ESS.Median 176.5460 214.4618 154.4806 205.8423 125.5458 249.3967
#> Tail_ESS.Mean 176.5460 214.4618 154.4806 205.8423 125.5458 249.3967
#> Tail_ESS.3rd Qu. 176.5460 214.4618 154.4806 205.8423 125.5458 249.3967
#> Tail_ESS.Max. 176.5460 214.4618 154.4806 205.8423 125.5458 249.3967
#> nu_Intercept sigmaA_groupa sigmaA_groupb sigmaB_groupa
#> Bulk_ESS.Min. 743.1364 523.1536 625.59329 743.1364
#> Bulk_ESS.1st Qu. 743.1364 523.1536 625.59329 743.1364
#> Bulk_ESS.Median 743.1364 523.1536 625.59329 743.1364
#> Bulk_ESS.Mean 743.1364 523.1536 625.59329 743.1364
#> Bulk_ESS.3rd Qu. 743.1364 523.1536 625.59329 743.1364
#> Bulk_ESS.Max. 743.1364 523.1536 625.59329 743.1364
#> Tail_ESS.Min. 242.3206 247.4966 80.60176 275.4082
#> Tail_ESS.1st Qu. 242.3206 247.4966 80.60176 275.4082
#> Tail_ESS.Median 242.3206 247.4966 80.60176 275.4082
#> Tail_ESS.Mean 242.3206 247.4966 80.60176 275.4082
#> Tail_ESS.3rd Qu. 242.3206 247.4966 80.60176 275.4082
#> Tail_ESS.Max. 242.3206 247.4966 80.60176 275.4082
#> sigmaB_groupb sigmaC_groupa sigmaC_groupb
#> Bulk_ESS.Min. 743.1364 579.3801 743.1364
#> Bulk_ESS.1st Qu. 743.1364 579.3801 743.1364
#> Bulk_ESS.Median 743.1364 579.3801 743.1364
#> Bulk_ESS.Mean 743.1364 579.3801 743.1364
#> Bulk_ESS.3rd Qu. 743.1364 579.3801 743.1364
#> Bulk_ESS.Max. 743.1364 579.3801 743.1364
#> Tail_ESS.Min. 172.1869 193.6446 221.6820
#> Tail_ESS.1st Qu. 172.1869 193.6446 221.6820
#> Tail_ESS.Median 172.1869 193.6446 221.6820
#> Tail_ESS.Mean 172.1869 193.6446 221.6820
#> Tail_ESS.3rd Qu. 172.1869 193.6446 221.6820
#> Tail_ESS.Max. 172.1869 193.6446 221.6820
#>
#> $ESS$complete
#> par Bulk_ESS Tail_ESS model
#> 1 nu_Intercept 743.1364 242.32063 1
#> 2 A_groupa 714.7464 176.54597 1
#> 3 A_groupb 593.0616 214.46185 1
#> 4 B_groupa 621.8739 154.48062 1
#> 5 B_groupb 681.5482 205.84233 1
#> 6 C_groupa 561.5601 125.54581 1
#> 7 C_groupb 743.1364 249.39665 1
#> 8 sigmaA_groupa 523.1536 247.49661 1
#> 9 sigmaA_groupb 625.5933 80.60176 1
#> 10 sigmaB_groupa 743.1364 275.40825 1
#> 11 sigmaB_groupb 743.1364 172.18693 1
#> 12 sigmaC_groupa 579.3801 193.64461 1
#> 13 sigmaC_groupb 743.1364 221.68196 1
#>
#>
#> $priorPredictive
#> $priorPredictive[[1]]
#>
#>
#> $posteriorPredictive
#> $posteriorPredictive[[1]]
#>
#>
fit_2 <- fit_test
fit_list <- list(fit_test, fit_2)
x <- barg(fit_list, list(ss, ss))
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: Coercing LHS to a list
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: Coercing LHS to a list
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
#> Warning: The ESS has been capped to avoid unstable estimates.
# }