bayesics: Core Statistical Methods via Bayesian Inference in R

JSS Journal of Statistical Software MMMMMM YYYY, V olume VV, Issue II. doi: 10.18637/jss.v000.i00 ba y esics: Core Statistical Metho ds via Ba y esian Inference in R Daniel K. Sew ell Univ ersit y of Iow a Alan Arakkal JANE LLC Abstract Ba yesian statistics is an in tegral part of contemporary applied science. ba yesics pro- vides a single framew ork, unified in syntax and output, for p erforming the most commonly used statistical pro cedures, ranging from one- and tw o-sample inference to general me- diation analysis. bay esics leans hard a w a y from the requirement that users be familiar with sampling algorithms b y using closed-form solutions whenev er possible, and auto- matically selecting the num b er of p osterior samples required for accurate inference when suc h solutions are not possible. ba y esics fo cuses on providing key inferential quan tities: p oin t estimates, credible in terv als, probability of direction, region of practical equiv alance (R OPE), and, when applicable, Bay es factors. While algorithmic assessment is not re- quired in bay esics , mo del assessmen t is still critical; tow ards that, ba y esics pro vides di- agnostic plots for parametric inference, including Bay esian p-v alues. Finally , ba y esics pro vides extensions to mo dels implemented in alternativ e R pac kages and, in the case of mediation analysis, correction to existing implemen tations. K eywor ds : Ba yesian inference, one p opulation inference, t w o p opulation inference, regression, Ba y esian non-parametrics, R . 1. In tro duction As is done collo quially , Bay esian statistics equates pr ob ability with unc ertainty . Subsequently , a Ba yesian definition of probabilit y in tends to quantify the amoun t of uncertain t y , or lac k of kno wledge, about a particular truth or ev en t. Bay esian analyses attempt to directly answer scien tific queries, such as, “Is my hypothesis true?” or “Ho w certain are we that there is a p ositive relationship b etw een x and y ?”, leading to directly interpretable quantities. In con trast, frequen tist inference is often difficult to interpret or requires additional steps to mak e scien tific conclusions. F or example, confidence interv als only allo w for in terpretation ab out the pro cess and not ab out an y sp ecific realized confidence in terv al. The other frequen- 2 ba y esics: Core Statistical Methods via Bay esian Inference in R Figure 1: Prop ortion of p eer-reviewed articles catalogued by Eur op e PMC that con tain either “credible interv al” or “confidence interv al” from 1981 - 2025. tist mainstay is p-v alues, which are “to o often misundersto o d and misused in the broader researc h comm unit y” ( W asserstein and Lazar 2016 ) and ha v e con tributed to significan t re- pro ducibilit y concerns ( Colquhoun 2017 ). Indeed, Greenland, Senn, Rothman, Carlin, Poole, Go o dman, and Altman ( 2016 ) lists a large num ber of wa ys frequen tist estimates and pro- cedures are commonly misin terpreted, man y of which corresp ond to correct interpretations from a Bay esian framew ork. Thanks to prominent figures (de Finetti, DeGro ot, Jeffreys, Sa v age, etc.) the Bay esian paradigm is well founded philosophically . It has a strong grounding in b oth theory and compu- tation, is em bedded in con temp orary applied science and machine learning ( Bon, Bretherton, Buc hhorn, Cramb, Dro v andi, Hassan, Jenner, Mayfield, McGree, Mengersen, Price, Salomone, San tos-F ernandez, V ercelloni, and W ang 2023 ), and is becoming more widely accepted in the life sciences by both industry and go v ernmen tal regulatory b o dies ( Rosner, Laud, and John- son 2021 ). Figure 1 sho ws the proportion of all p eer-review ed articles catalogued b y Eur op e PMC that include the phrase “credible interv al” as a proxy (and at least a low er b ound) for how often Ba y esian inference is being conducted in biomedical research. Starting around 2003, Bay esian inference has b ecome drastically more widespread. Ho w ev er, Figure 1 also sho ws the trend for the search term “confidence interv al,” illustrating that despite the dra- matic up tic k in the use of Ba y esian statistics, frequentist statistics still dominate the life sciences. While the reasons for this are w ell b eyond the scop e of this pap er, the authors con- jecture that a primary driv er ma y b e the circularity of science and teaching described in the op ening of W asserstein and Lazar ( 2016 )- frequen tist statistics is taugh t b ecause that is what regulators and journal editors exp ect, and regulators and editors expect frequentist statistics b ecause that is what they w ere taugh t. On the other hand, the increased usage of Ba y esian statistics has b een asserted to be in large part the product of increased computational to ols (e.g., Hagan and W est 2010 ) and statistical softw are implemen ting Bay esian methodology (e.g., Štrumbelj, Bouchard-Côté, Corander, Gelman, Rue, Murray , Pesonen, Plummer, and V ehtari 2024 ). There are a large num ber of R pac kages a v ailable for obtaining Bay esian inference for very Journal of Statistical Softw are 3 sp ecific tasks, man y of which corresp ond to adv anced statistical metho ds (causal analysis, net w ork analysis, spatial statistics, etc.); for example, bama ( Rix, Kleinsasser, and Song 2025 ) has the sole purp ose of p erforming Ba y esian inference for high-dimensional linear me- diation mo dels. On the other extreme lie adv anced R packages aimed at building mo dels from scratch and then obtaining Bay esian inference using them, most n otably greta ( Golding 2019 ), LaplacesDemon ( Statisticat and LLC. 2021 ), NIMBLE ( de V alpine, T urek, P aciorek, Anderson-Bergman, T emple Lang, and Bodik 2017 ), and rstan ( Stan Dev elopmen t T eam 2025 ). More closely related to the aims of the newly dev eloped R package ba y esics , the subject of this man uscript, are those R packages which aim to do a set of commonly implemented analyses using Ba y esian inferential tec hniques. Ba y esF actor ( Morey and Rouder 2024 ) can p erform Ba y esian h ypothesis testing for ANO V A, linear regression, and correlation (for t wo normally distributed samples), but, as the name suggests, provides only Bay es factors, and do esn’t provide point or interv al estimation. F or common one- and tw o-sample inferen tial pro cedures, DFBA ( Barch and Chec hile 2023 ) pro vides non-parametric Ba y esian approac hes, but no regression metho ds. In contrast, arm ( Gelman and Su 2024 ) pro vides only regression metho ds. While arm relies primarily on the large sample normal approximation for p osterior inference, more recen t R packages that are fo cused on regression rely on Stan , a probabilistic programming language for specifying statistical mo dels used primarily to implemen t No- U-T urn Hamiltonian Monte Carlo sampling ( Carp en ter, Gelman, Hoffman, Lee, Go o drich, Betancourt, Brubak er, Guo, Li, and Riddell 2017 ). T wo of the most notable suc h R packages are rstanarm whic h has a large array of precompiled statistical regression mo dels ( Go o drich, Gabry , Ali, and Brilleman 2025 ), and brms whic h pro vides a highly flexible approach to building complex regression and hierarc hical mo dels ( Bürkner 2017 ). These t w o packages are v ery p o w erful and useful for regression metho ds, yet it is all to o easy to use these pac kages without obtaining reliable inference. As discussed in Section 2.2 , kno wing th e Mark o v chains ha v e conv erged is insufficient for determining if one has reliable inference, and the Monte Carlo standard error (MCSE) is quite misleading in this regard as it only p ertains to the p oint, rather than in terv al estimate. F urthermore, these packages do not automatically pro vide the critical inferen tial quan tities needed to p erform a Ba y esian analysis (see Section 2.1 ). It should b e noted that this has largely b een rectified b y the well-dev elop ed R package ba y estestR ( Mak o wski, Ben-Shac har, and Lüdec k e 2019 ), although some issues remain with the default R OPE b ounds, namely adjusting according to the scale of the co v ariates. T o address these issues, w e hav e dev elop ed the ba y esics package. As its name suggests, its o v erall goal is to pro vide a single framew ork, unified in syn tax and output, for p erforming the most commonly used statistical analyses using Ba y esian inference. The analyses provided in ba y esics range from simple one- and t w o-sample analyses to somewhat more adv anced regres- sion tec hniques such as model a v eraging and general mediation analysis. T able 1 provides a full list of the functionality of bay esics . Bey ond this, one of the most important motiv ations for bay esics was to emphasize inference ov er algorithms, and to provide wa ys to b etter un- derstand/in terpret the output, including visualization to ols. It seems non-contro versial that to p erform go o d science, it is sufficient to be able to use and correctly interpret the output of, e.g., stats::glm() without in-depth kno wledge of ho w stats::glm() obtains its results via iterativ ely rew eigh ted least squares. Thus, the authors feel strongly that Bay esian infer- en tial procedures ought to b e accessible in the same w a y . T o wards that, ba yesics implements metho ds that either require no algorithmic tuning (i.e., closed-form solutions) or do so auto- 4 ba y esics: Core Statistical Methods via Bay esian Inference in R matically with minimal input from the user (namely , accuracy tolerance). Imp ortan tly , w e fo cus on the in terv al estimates rather than the p oint estimates, as these are often, rightly or wrongly , what are used to determine a “statistically significant” result (e.g., if the 95% credible interv al for a regression co efficient do es not include zero) and require considerably more posterior samples to obtain sufficien t precision. bay esics also pro vides diagnostic plots for model assessmen t (none are needed for algorithmic assessmen t) including the implemen- tation of Ba y esian p-v alues for generalized linear models (GLMs), as w ell as non-parametric alternativ es to regression mo dels should mo del assumptions fail. Finally , bay esics provides extensions to mo dels implemented in alternativ e R pac kages and, in the case of mediation analysis, corrections to existing methods. In the remainder of the paper, w e go into more detail on the emphasis on inference and understanding o v er algorithms (Section 2 ), mo del diagnostic plots and remediations (Section 3 ), and finally metho d extensions and corrections (Section 4 ). While Sections 2 - 4 include a few examples of v arious bay esics functions, Section 5 provides a fuller demonstration. W e end with a discussion in Section 6 . T able 1: F unctions, analyses, and a v ailable generics. Dep en- dencies are noted in the footnotes. F unction Analysis Generics aov_b 1-w a y analysis of v ariance coef , credint 1 , *IC 2 , plot , predict , print , summary , vcov bma_inference 3 Ba y esian mo del a v eraging for linear regression mo dels coef , credint , plot , predict , print , summary case_control_b Analysis of 2x2 con tingency tables from a case-control study chisq_test_b Indep endence analysis for 2- w a y contingency tables cor_test_b 4 Analysis of Kendall’s tau cor- relation co efficient find_beta_parms Find shap e parameters of the Beta distribution whic h tries to matc h a user-sp ecified mean and quantile Journal of Statistical Softw are 5 F unction Analysis Generics find_invgamma_parms Find shap e and rate param- eters of the Inv erse Gamma distribution that either tries to match tw o user-sp ecified quan tiles, or else finds an In v erse Gamma distribution that matc hes a priori confi- dence in the co efficien t of de- termination ( R 2 ) glm_b Generalized linear mo dels bayes_factors , coef , credint , *IC , plot , predict , print , summary , vcov heteroscedasticity_test T est whether tw o or more p opulations ha v e equal v ari- ances via Bay es factors lm_b Linear mo dels bayes_factors , coef , credint , get_posterior_draws , *IC , plot , predict , print , summary , vcov mediate_b Mediation analysis using the framew ork of Imai, Keele, and Tingley ( 2010 ) plot , print , summary np_glm_b Generalized Ba y esian infer- ence via the loss-likelihoo d b o otstrap coef , credint , plot , predict , print , summary , vcov poisson_test_b Mak e inference on one or t w o p opulations using Poisson dis- tributed count data prop_test_b Mak e inference on a sin- gle p opulation prop ortion or compare tw o p opulation pro- p ortions sign_test_b Sign test for paired data surv_fit_b Semi-parametric estimation of surviv al curves for right- censored data for one or more p opulations bayes_factors , plot , print t_test_b Mak e inference on one or t w o p opulation means using nor- mally distributed data 6 ba y esics: Core Statistical Methods via Bay esian Inference in R F unction Analysis Generics wilcoxon_test_b 5 Wilco xon Rank Sum (aka Mann-Whitney U) or the Wilco xon signed rank analy- ses 1 credint() is the Bay esian equiv alent of stats::confint() 2 *IC represents the follo wing information criteria: AIC, BIC, DIC, and W AIC 3 Dep ends on BMS::bms 4 Dep ends on DFBA::dfba_bivariate_concordance 5 Dep ends on DFBA::dfba_mann_whitney , DFBA::dfba_wilcoxon 2. Inference, not Algorithms 2.1. What can Bay esians do for you? W e assume that the reader has at least a reasonable familiarit y with Bay esian statistics. F or those who need a primer, there exist coun tless in tro ductory texts on the sub ject (e.g., Rosner et al. 2021 ). F or those w ell v ersed in Ba yesian analysis, w e recommend skipping to Section 2.2 . Before pro ceeding, we will first cov er at a very high level the fundamental inferential quan tities used b y the practicing Bay esian statistician, as these are referred to throughout the remainder of the manuscript. The first component is, of course, the point estimate. Typically used is the p osterior mo de, posterior median, or the p osterior mean. Eac h of these has a strong theoretical justification ro oted in decision theory , but it is sufficient here to sa y that our p osterior mean, for example, is our b est guess at the truth. This has very limited practical implications, how ever, as our b est guess is undoubtedly wrong 1 . Our second key inferen tial quantit y is credible interv als (CIs), whic h provide a b etter ap- proac h to understanding what the truth migh t be. CIs are precisely what practitioners wish confidence interv als were, and unfortunately , the latter are often misinterpreted as the former (this is, in fact, the very first misin terpretation of confidence interv als listed by Greenland et al. 2016 ). If a Ba y esian analysis pro vides a 95% CI for an estimand, we sa y correctly that w e are 95% confident that the truth lies in our in terv al. Third, when w e are p erforming a causal or asso ciation study , we wish to quantify our certaint y that w e understand the direction of the relationship. In regression analysis, this corresp onds to statemen ts of confidence ab out whether the regression co efficien t is, e.g., p ositiv e. This quan tit y is called, appropriately enough, the probabilit y of direction, or simply PDir. While the PDir is imp ortant for c haracterizing our knowledge of the direction of a relationship, it is wholly insufficien t for telling us whether that relationship is practically meaningful. T o address this, practitioners ought to define a region of practical equiv alence (R OPE). With this, a Ba yesian analysis can determine ho w certain we are that the estimand lies within this R OPE. Kruschk e ( 2018 ) provides a well thought out discussion on default ROPE bounds, whic h bay esics adheres to for default settings. 1 Unless our estimand takes one of a discrete set of v alues. Journal of Statistical Softw are 7 Finally , w e sometimes wish to compare the lik elihoo d of one statemen t’s fidelit y against that of a con tradictory statemen t. F or example, we ma y wish to compare head-to-head the statemen t, “exp osure to c hemical A is asso ciated with an increased c hance of disease B,” with the statement “exposure to c hemical A is not asso ciated with disease B. ” Bay es factors tell us how our prior b eliefs ab out the veracit y of these t w o comp eting statemen ts change based on the data, regardless of what those prior beliefs are. More sp ecifically , the Bay es factor will scale the prior odds that statement 1 is true vs. statemen t 2 is true to obtain these same o dds after having seen the data. These are the most well used and imp ortan t pieces of a Bay esian analysis, and it is the goal of ba y esics to provide and facilitate in terpretation of these pieces while allo wing the user to neglect the metho d of ho w these quan tities were obtained algorithmically . 2.2. Monte Carlo accuracy Ba y esian inference is often p erformed in practice through p osterior sampling- usually MCMC. Th us it is imp ortant to understand the lev el of precision obtained b y a specified num ber of p osterior samples. Most implemented algorithms use some measure of conv ergence diagnostic, e.g., the Gewek e diagnostic ( Gew ek e 1991 ) or ˆ R ( Gelman, Carlin, Stern, Dunson, V ehtari, and Rubin 2004 ), and the Monte Carlo standard error (MCSE) for the posterior mean. The former is clearly critical to kno w whether the samples are v alid, and the latter is intended to indicate whether or not sufficient n um b ers of v alid samples ha v e b een obtained. Y et the authors susp ect that unless the user is well v ersed in p osterior sampling and MCMC, the former measures of chain conv ergence are used to determine whether the rep orted inference is accurate. How ev er, ev en w ere the MCSE to be used to determine whether sufficient samples w ere obtained, this measure can b e highly misleading for anything other than p oint estimation. The imp ortant work in Doss, Flegal, Jones, and Neath ( 2014 ) provides a central limit theorem for order statistics for MCMC-derived samples. F or explication purp oses, we will use the ideal setting of iid samples, although the problem b ecomes further exacerbated b y high auto corre- lation in the c hain. F or a parameter of interest θ , let ˆ θ L,p denote the p th empirical quan tile of an iid sample of size L (note that L represents the n umber of p osterior samples, not the sample size of the original data), let θ p denote the corresp onding true posterior quan tile, and let π ( ·| y ) denote the posterior probability densit y function of θ | y . Then √ L ( ˆ θ L,p − θ p ) d → N 0 , p (1 − p ) π 2 ( θ p | y ) ! . (1) Hence, to estimate, sa y , the low er b ound of a (1 − α )100 % CI, i.e., the ( α / 2) th p osterior quan tile of θ | y , up to ± some error ϵ with (high) probability s , i.e., s = Pr( θ α 2 − ϵ < ˆ θ L, α 2 < θ α 2 + ϵ ) , w e set L to be L = α 2  1 − α 2  Z 1 − s 2 ϵπ ( θ α 2 | y ) ! 2 (2) The MCSE comes in to pla y when determining the num ber of samples required to obtain accurate estimation of the p osterior mean. Let ¯ θ denote the true p osterior mean and ˆ ¯ θ M 8 ba y esics: Core Statistical Methods via Bay esian Inference in R denote the estimate of the p osterior mean from M iid posterior samples. Then to estimate ¯ θ up to ± ϵ with probabilit y s , the n um b er of samples required is M = V ar  θ ( m )  Z 1 − s 2 ϵ ! 2 . (3) where θ ( m ) is the m th iid p osterior dra w of θ . Putting these together, w e quic kly see that for whatever probability s and Mon te Carlo (MC) margin of error ϵ , the ratio of samples required to obtain such precision for the low er credible interv al bound to that required for the p osterior mean is giv en b y α 2  1 − α 2  V ar  θ ( m )  π 2  θ α 2 | y  (4) Figure 2 illustrates the egregious nature of this problem when considering ho w man y more samples are required for tail quan tile estimation than for p osterior mean estimation. F rom this, w e can see there is the p oten tial for several orders of magnitude more samples to b e required for equiv alent precision in the in terv al estimates as for the p oin t estimates. T o further illustrate this, we simulated data according to a simple linear regression mo del with in tercept = 0 , slope = 0 . 25 , and residual standard deviation = 1 , and a sample size of 25. W e then emplo y ed rstanarm , using their default settings, to estimate using 500 simulations the standard error for the low er b ound of the 95% credible interv al for the slop e and compared that with the MCSE reported for the p osterior mean. A cross 500 sim ulations, the rep orted MCSE w as consisten t, ranging from 0.0034 to 0.0047. Ho w ev er, the empirical standard error of the lo w er b ound of the 95% credible interv al w as 0.0133, roughly 3-4 times the error reported by the MCSEs. 1 10 100 0.900 0.925 0.950 0.975 CI le vel Sample size ratio Distribution Beta(5,2) Gamma(10,1) N(0,1) t(5) Figure 2: Illustration of the ratio of posterior samples required for the lo w er credible in terv al b ound vs. the p osterior mean to be within the same margin of error. CI lev els range from 90% to 99%. Journal of Statistical Softw are 9 2.3. bay esics automated approac h When p ossible, ba yesics av oids the requiremen t for Monte Carlo n umerical approximation by using conjugacy . A simple example of this is lm_b , whic h uses a normal-in v erse gamma prior on the regression coefficients and the residual v ariance for linear regression. That is, for the familiar mo del given b y y = X β + ϵ , we use the priors β | σ 2 ∼ N ( µ , σ 2 V − 1 ) , σ 2 ∼ Γ − 1 ( a/ 2 , b/ 2) , (5) whic h yields a multiv ariate t distribution for the marginal prior and marginal p osterior for β . The default setting in bay esics is a Zellner’s g prior with g equalling the sample size. Other priors are a v ailable, notably that which corresp onds to the statemen t, “we are 95% sure a priori that a standard deviation increase in an y co v ariate will not change the mean of y b y more than a factor of 5 standard deviations of y ,” (obtained when prior="conjugate" ). The default hyperparameters for σ 2 are motiv ated by placing a prior on the co efficien t of determination R 2 . Sp ecifically , since the residual v ariance is (1 − R 2 ) times the observed v ariance, w e set 50% prior probabilit y that σ 2 is b et w een (1 − 0 . 9 2 ) s 2 y and (1 − 0 . 1 2 ) s 2 y , where s 2 y is the sample v ariance of the resp onse v ariable, i.e., we are 50% confident a priori that the correlation b etw een the observed and the fitted v alues will be b etw een 0.1 and 0.9 in magnitude. In the case of GLMs, no conjugate prior exists. Instead, bay esics implemen ts through glm_b a fixed form v ariational Bay es (VB) (as w ell as imp ortance sampling and large sample ap- pro ximations, although VB is the default), whic h finds the multiv ariate normal distribution with unconstrained posterior co v ariance matrix whic h is as close as p ossible to the p osterior with resp ect to the Kullbac k-Leibler div ergence ( Salimans and Knowles 2013 ). W e note that this approac h is also implemented in, e.g., rstanarm::stan_glm if one sets the argumen t algorithm="fullrank" . In other cases where there is no conjugate prior a v ailable, obtaining iid p osterior samples is trivially accomplished, and in such cases bay esics automatically determines the num ber of p osterior dra ws. F or example, supp ose w e wish to study the rate ratio from tw o indep endent samples which are counts, i.e., w e hav e y j ∼ P oisson ( λ j ω j ) , j = 1 , 2 , where λ j is the pop- ulation rate and ω j is the offset term. Using a conjugate gamma prior for each individual p opulation rate allo ws large n um b ers of p osterior dra ws of λ j | y j , and th us p osterior dra ws of λ 1 /λ 2 | y 1 , y 2 , to b e obtained trivially . The poisson_test_b draws a preliminary sample of size 500 whic h is then used to estimate the densit y at the desired quan tile corresponding to the lo w er and upp er bounds of the credible in terv al. This is the k ey quan tit y in determining the n um b er of p osterior samples as seen in Equation 2 , and poisson_test_b will automatically dra w the remaining samples necessary for the desired level of precision. Sometimes a mix of these tw o approac hes- closed-form solutions and sampling- is implemen ted. In aov_b , whic h implements one-wa y ANOV A, some estimands suc h as individual factor level means ha v e closed-form posterior distributions, while others such as exceedance in pairs rate (EPR) 2 do not. F or these latter cases w e can still easily obtain iid posterior samples and, as with poisson_test_b , the n um b er of suc h draws is determined in an automated fashion. In this w a y , the user of bay esics should v ery rarely need to b e concerned with an y 2 EPR is a quantit y describ ed in Rosner et al. ( 2021 ) p.136 measuring the probability that a random v alue from one p opulation is greater than a random v alue from another p opulation, i.e., how well separated the tw o p opulation distributions are. 10 ba y esics: Core Statistical Methods via Bay esian Inference in R asp ect of the algorithms used, and can trust that the output will b e correct and accurate. 3. Diagnostics and Remediations 3.1. Diagnostic plots While diagnostic plots relating to algorithms, e.g., MCMC trace plots, are not required by ba y esics algorithms, it is still crucial to assess a mo del’s go o dness of fit. T o w ards that, for one-w a y ANOV A ( aov_b ) and linear regression ( lm_b ) mo dels, the plot function provides the standard residual vs. fitted v alues plot as w ell as a qq-norm plot (among other options). F or GLMs, it is just as imp ortan t to assess mo del fit. bay esics addresses this need through the use of Ba y esian p-v alues. Briefly , Bay esian p-v alues, also called p osterior predictive p- v alues, assess the compatibilit y of a certain feature of a mo del to the observed data. Letting y obs ∈ Y denote the observed data, y pred ∈ Y denote data coming from the p osterior predictive distribution, θ ∈ Θ b e the prop osed mo del’s parameters, and T : Y × Θ 7→ ℜ b e the test statistic, the Bay esian p-v alue is defined to be Pr( T ( y pred , θ ) < T ( y obs , θ ) | y obs ) . (6) Note that the test statistic can b e a function of both the observ ables and the model param- eters. If the mo del is correct, or at least correct in terms of what the test statistic T ( · , · ) is assessing, then the mo del ough t to predict v alues similar to the observed data. In such a case, w e would hop e that the Bay esian p-v alue would near 0.5. If, how ev er, we obtain a p-v alue near 0 or 1, say , < 0 . 05 or > 0 . 95 , then we must conclude that the asp ect of the mo del captured b y our test statistic is inconsistent with the observ ed data. F or more details, see, e.g., Gelman et al. ( 2004 ). bay esics pro vides scatterplots of p osterior dra ws of ( T ( y pred , θ ) , T ( y obs , θ )) pairs, and provides the p-v alue estimates. As a simple example, consider regression data generated according to the negative binomial distribution. R> set.seed(2026) R> N = 500 R> nb_data <- + tibble(x1 = rnorm(N), + x2 = rnorm(N), + x3 = rep(letters[1:5],each = N/5), + time = rexp(N)) |> + mutate(outcome = + rnbinom(N, + mu = + exp(-2 + x1 + 2 * (x3 %in% c("d","e"))) * + time, + size = 0.7)) If w e fit a GLM b ased on the Poisson distribution, the ov erdispersion ought to be reflected Journal of Statistical Softw are 11 in a p-v alue near 0 or 1, whereas when w e fit a GLM based on the correct distribution 3 , the p-v alue ought to b e near 0.5. Indeed, this is precisely what we see in Figure 3 . R> poisson_fit <- + glm_b(outcome ~ ., + data = nb_data, + family = poisson()) R> nb_fit <- + glm_b(outcome ~ ., + data = nb_data, + family = negbinom()) R> p1 <- + plot(poisson_fit, + type = "dx") + + labs(subtitle = "Poisson assumed") R> p2 <- + plot(nb_fit, + type = "dx") + + labs(subtitle = "Neg. binomial assumed") R> patchwork::wrap_plots(p1,p2) 1160 1170 1180 700 800 900 T ( y pred , β ) T ( y obs , β ) P oisson assumed Bay esian p−v alue based on deviance = 1 810 820 830 600 700 800 900 1000 T ( y pred , β ) T ( y obs , β ) Neg. binomial assumed Bay esian p−v alue based on deviance = 0.443 Figure 3: Bay esian p-v alues for regression data generated according to the negative binomial Finally , Bay esian p-v alues are also implemented for Bay esian mo del a v eraging (more in Section 4.2 ), where the test statistics used are quan tiles of the data. 3.2. Non-parametric metho ds If no candidate mo dels are consistent with the observed data in terms of the Ba y esian p- v alues or other diagnostic chec ks, it is b etter to opt for a non-parametric approac h. F or simple scenarios suc h as paired analyses, we can use pro cedures such as the Ba y esian sign 3 The negbinom() family is provided in bay esics . 12 ba y esics: Core Statistical Methods via Bay esian Inference in R test ( sign_test_b ). The DFBA pac kage pro vides go o d implemen tations of many suc h one and tw o sample tec hniques, and rather than rein ven t the wheel, we ha ve incorp orated several of their functions in to the common in terface and the inferen tial output seen elsewhere in ba y esics . F or fitting surviv al curv es to one or more samples and/or testing group equiv alence, we ha v e implemen ted the semi-parametric approac h of Qing, Thall, and Y uan ( 2023 ). In this setup, the hazard function is assumed to b e a piecewise step function, or equiv alen tly , the time- to-ev en t data follow a piecewise exp onen tial distribution. This clever mod eling approach is highly flexible, yet pro vides closed-form p osteriors for the hazard rates, as well as a closed- form solution for the marginal lik eliho o d. The latter is imp ortant as it provides a ready approac h for (1) determining the num ber of knots in the hazard function, and (2) testing whether all samples’ time-to-ev en t outcomes follo w the same or differen t distributions via Ba y es factors. The co de given b elow sho ws the analysis of surviv al times from breast cancer patien ts in the GBSG2 study (see ?TH.data::GBSG2 for more information and references), testing to see if surviv al rates v ary by whether the patien t receiv ed hormonal therap y , with the estimated surviv al curves giv en in Figure 4 . Note that whenever ba y esics pro vides Ba y es factors, interpretations and direction are alw ays made clear. R> data(GBSG2, + package = "TH.data") R> GBSG2_fit_horTh <- + survfit_b(Surv(time,cens) ~ horTh, + data = GBSG2) R> GBSG2_fit_aggregated <- + survfit_b(Surv(time,cens) ~ 1, + data = GBSG2) R> bayes_factors(GBSG2_fit_aggregated, + GBSG2_fit_horTh) F or linear mo dels and GLMs, ba y esics implemen ts the loss-likelihoo d b o otstrap, a non- parametric metho d using general Bay esian up dating ( Lyddon, Holmes, and W alk er 2019 ). Here, we consider the mo del parameters to b e a functional of the true distribution of our observ ables inv olving a loss function, i.e., β ( F ) := argmin ˜ β Z Y ℓ ( ˜ β , y ) dF ( y ) , (7) where β ( · ) is our functional parameter, F is the true distribution of our observ ables y ∈ Y , and ℓ ( · , · ) is our loss function mapping the parameter space and Y to the p ositive reals. Uncertain t y in β is solely the result of uncertain t y in the distribution of y . Thus, rather than needing to sp ecify a prior on β , the loss-lik eliho o d b o otstrap uses the degenerate Dirichlet pro- cess prior on the distribution of the observ ables that corresp onds to the Ba y esian b o otstrap. Therefore the approac h of Lyddon et al. ( 2019 ) has the arguable double b enefit of not needing to specify a prior on the mo del parameters and not having to assume a correct lik elihoo d. As a default, the self-information loss is used (i.e., the loss function is the negativ e of the log likelihoo d), although an y user-specified loss functions can b e implemented. T o illustrate, consider a regression data with mis-specified error distribution and co v ariate structure. Journal of Statistical Softw are 13 R> plot(GBSG2_fit_horTh) 0.00 0.25 0.50 0.75 1.00 0 1000 2000 Time S(t) Group horTh: no horTh: yes Figure 4: Semi-parametric estimated surviv al curves for adv anced lung cancer, disaggregated b y sex. R> set.seed(2026) R> linreg_data <- + tibble(x = runif(100)) |> + mutate(y = 20 * x^2 + rgamma(100,shape = 2, rate = 0.5)) R> lm_fit <- + lm_b(y ~ x, + data = linreg_data) R> np_fit = + np_glm_b(y ~ x, + data = linreg_data, + family = gaussian()) F rom Figure 5 , we can see that the diagnostic plots demonstrate violated assumptions (as exp ected), and the credible interv als are appropriately wider for the non-parametric compared to the parametric approac h. Note that when it comes to in terpretation, the regression line estimated by np_glm_b is not a “true” regression line, but rather the regression line that w ould b e fit w ere w e to p erfectly know the population distribution. 14 ba y esics: Core Statistical Methods via Bay esian Inference in R R> plot(lm_fit, + type = "dx") + + ( + plot(lm_fit, + type = "cred band") + + labs(subtitle = "(Parametric)") + ) + + ( + plot(np_fit, + type = "cred band") + + labs(subtitle = "(Non-parametric)") + ) −5 0 5 10 15 5 10 15 20 y ^ ε ^ Fitted vs. Residuals −5 0 5 10 15 −2 −1 0 1 2 Theoretical quantiles Empirical quantiles QQ norm plot 0 10 20 30 0.00 0.25 0.50 0.75 1.00 x y (P arametric) Credible band f or x 0 10 20 30 0.00 0.25 0.50 0.75 1.00 x as.numeric(y) (Non−parametric) Credible band f or x Figure 5: A mis-sp ecified linear regression mo del fit using a parametric approach and the non-parametric loss-likelihoo d b o otstrap. Journal of Statistical Softw are 15 4. Metho d Extensions and Corrections 4.1. Information criteria Information criteria play an important role in model selection and ev aluation. ba y esics pro- vides four of the most commonly used information criteria (IC): Akaik e IC (AIC), Bay esian IC (BIC), Deviance IC (DIC), and Widely Applicable IC (W AIC). Despite the popularity of these four criteria, with the exception of W AIC they are not all readily av ailable as part of other widely used Ba y esian packages. In bay esics , all four are implemen ted for aov_b , lm_b , and glm_b . 4.2. Bay esian mo del a v eraging A p o w erful alternative to selecting and p erforming inference for a single mo del is Ba y esian mo del a v eraging. This approach yields a mixture distribution as the p osterior, where the mixture comp onents are the p osteriors obtained from eac h sp ecific regression model and the comp onen t w eigh ts are the p osterior probabilities that the corresp onding models are correct. The BMS package pro vides functionality to p erform Bay esian mo del av eraging. Its ep onymous function bms has considerable flexibility in terms of the prior on the mo del set, and the pac kage pro vides man y wa ys to describe the p osterior o v er the mo del space and parameter inclusion probabilities. Y et the capacit y for inference on the co efficients themselves is minimal, whic h is where ba yesics comes in. ba y esics obtains the p osterior model probabilities from wrapping BMS::bms , and then leverages lm_b to obtain inference. P osterior sampling is required for this, but again bay esics automatically determines the necessary n um b er of p osterior samples to obtain accurate interv al estimates. 4.3. Mediation analysis Mediation analysis is an important method to gain a more n uanced understanding betw een predictor v ariables and an outcome of in terest. The type of relationships assessed in a me- diation analysis is depicted in Figure 6 . Imai et al. ( 2010 ) describ ed a highly flexible causal mediation framew ork in whic h the relationship betw een the treatmen t and the mediator and the relationship b etw een the outcome and b oth the treatment and mediator can correspond to any sort of mo del. How ev er, their prop osed estimation approach, so-called quasi-Bay esian, is lacking a rigorous framew ork. In fact, their estimation approac h corresp onds to a Ba y esian analysis using improp er uniform priors. The accompan ying R package mediation , ho w ev er, refers to “quasi-Bay esian confidence in terv als” and p-v alues. In bay esics , w e hav e provided fully Ba yesian implementation of the causal mediation framework of Imai et al. ( 2010 ) for linear and generalized linear models. 5. Demonstration As a brief demonstration of bay esics , we consider the randomized clinical trial describ ed in Elm unzer, Sc heiman, Lehman, Chak, Mosler, Higgins, Hayw ard, Romagnuolo, Elta, Sherman, W aljee, Repaka, A tkinson, Cote, K w on, McHenry , Piraka, W amsteker, W atkins, K orsnes, Sc hmidt, T urner, Nicholson, and F ogel ( 2012 ) studying the effects of indomethacin on post- 16 ba y esics: Core Statistical Methods via Bay esian Inference in R T reatment/ Exp osure ( X ) Outcome ( Y ) Confounders ( W ) Mediator ( M ) Figure 6: Schematic of relationships assessed via a mediation analysis. endoscopic retrograde cholangio-pancreatogram pancreatitus (PEP). These data are publicly a v ailable through the R pac kage medicaldata . While this study con tains a large n um b er of v ariables, we fo cus only on a subset for p edagogical purp oses. Sp ecifically , w e consider as confounder v ariables age, sex, and previous history of PEP (found in the v ariable risk ). The treatmen t arm is contained in the v ariable rx and equals 0 if placeb o and 1 if treatmen t. The binary outcome v ariable is contained in the v ariable outcome and indicates whether or not PEP o ccurred. While it not recommended in randomized trials to p erform hypothesis tests to see if cov ariates are balanced, w e shall do so here purely for demonstration purp oses. Age, treated as a con tin uous cov ariate, can b e inv estigated via t_test_b after loading in the data. Journal of Statistical Softw are 17 R> data(indo_rct, + package = "medicaldata") R> options(pillar.width = 100) R> # Effective randomization for age R> t_test_b(age ~ rx, + data = indo_rct) --- Bayes factor in favor of the full vs. null model: 1.21e-06; =>Level of evidence: Decisive --- Summary of factor level means --- # A tibble: 4 × 5 Variable ` Post Mean ` Lower Upper ` Prob Dir ` 1 Mean : rx : 0_placebo 46.0 44.6 47.5 1 2 Mean : rx : 1_indomethacin 44.5 42.9 46.0 1 3 Var : rx : 0_placebo 172. 147. 201. NA 4 Var : rx : 1_indomethacin 183. 155. 215. NA --- Summary of pairwise differences --- # A tibble: 1 × 9 Comparison ` Post Mean ` Lower Upper ` Prob Dir ` ` ROPE (0.1) ` EPR 1 0_placebo-1_indomethacin 1.56 -0.461 3.64 0.933 0.412 0.533 ` EPR Lower ` ` EPR Upper ` 1 0.490 0.577 *Note: EPR (Exceedence in Pairs Rate) for a Comparison of g-h = Pr(Y_(gi) > Y_(hi)|parameters) −25 0 25 0_placebo 1_indomethacin rx ε ^ Residual plot by group −50 −25 0 25 50 −2 0 2 Theoretical quantiles Empirical quantiles QQ norm plot 20 40 60 80 0_placebo 1_indomethacin rx age Cred. and Pred. intervals Figure 7: Analysis of randomization of age b etw een treatmen t arms. 18 ba y esics: Core Statistical Methods via Bay esian Inference in R These results indicate effective randomization through the Ba y es factor, credible interv al of the difference in means, PDir, R OPE, and EPR. It is important to recognize that t_test_b automatically pro vides diagnostic plots, and in this case, the normality assumption is a little off due to the outcomes being righ t sk ew ed. If w e w ere sufficien tly concerned, w e could instead run wilcoxon_test_b to implemen t the Bay esian rank-sum analysis (pow ered b y DFBA::dfba_mann_whitney ). Next, we can assess whether or not sex was effectiv ely randomized. T o ac hiev e this, we can p erform a simple test comparing tw o sample proportions. prop_test_b provides the p osterior means and credible in terv als (for illustrativ e purp oses, this time we elected to use 99% CIs) for eac h p opulation. Additionally , the probabilit y that the o dds ratio is in the R OPE is pro vided, where the ROPE is, as suggested in the excellen t discus sion of Krusc hk e ( 2018 ), half of a small effect size (p er FD A guidance on rate ratios). As seen in Figure 8 , the default prior is a Jeffrey’s prior, although this can b e c hanged to a uniform prior b y setting prior="uniform" or to a general Beta prior b y using the prior_shapes argumen t. R> gender_table = + table(indo_rct$gender, + indo_rct$rx) R> prop_test_b(gender_table[1,], + gender_table[2,], + CI_level = 0.99) 0 5 10 15 0.00 0.25 0.50 0.75 1.00 Distribution P oster ior (P op1) P oster ior (P op2) Prior P opulation propor tion Figure 8: Assessing randomization of sex (in this dataset, sex is binary) b et w een treatmen t arms. Finally , we can c hec k for randomization in the previous PEP risk v ariable, and this time, b ecause risk is a discrete ordered risk score, w e will use the Wilcoxon rank-sum analysis rather than t_test_b . While the Ba y es factor provides some evidence that randomization failed, the probabilit y of direction and the probabilit y of falling in the R OPE suggests this is not a ma jor concern. Journal of Statistical Softw are 19 R> wilcoxon_test_b( + indo_rct$risk[which(indo_rct$rx == "1_indomethacin")], + indo_rct$risk[which(indo_rct$rx == "0_placebo")] + ) 0 5 10 15 0.00 0.25 0.50 0.75 1.00 Distribution P oster ior Prior Ω x Figure 9: Assessing randomization of risk from previous PEP b etw een treatmen t arms. W e now mov e to the study’s primary analysis (or at least, this p edagogical example’s primary analysis): a cov ariate-adjusted generalized linear regression model to assess the effect of the indomethacin treatment on the PEP outcome. W e fit the model using glm_b , and then first lo ok at the diagnostic plots. 20 ba y esics: Core Statistical Methods via Bay esian Inference in R R> pep_fit = + glm_b(outcome ~ age + gender + risk + rx, + data = indo_rct, + family = binomial()) R> plot(pep_fit, + type = "diagnostics") 450 455 460 465 470 300 400 500 600 T ( y pred , β ) T ( y obs , β ) Ba y esian p−v alue based on de viance = 0.506 Figure 10: Diagnostic plots for the primary PEP analysis. Journal of Statistical Softw are 21 As we hav e a Ba y esian p-v alue close 0.5, w e can mo v e on to inference. If the p-v alue w ere close to 0 or 1, w e w ould ha v e instead used np_glm_b with iden tical syn tax. R> summary(pep_fit) ---------- Values given in terms of odds ratios ---------- # A tibble: 4 × 9 Variable ` Post Mean ` Lower Upper ` Prob Dir ` ROPE ` ROPE bounds ` 1 age 0.993 0.975 1.01 0.754 0.145 (0.998,1.002) 2 gender2_male 1.10 0.611 2.00 0.629 0.288 (0.889,1.125) 3 risk 1.53 1.17 2.00 0.999 0.00172 (0.967,1.034) 4 rx1_indomethacin 0.472 0.288 0.773 0.999 0.00572 (0.889,1.125) ` BF favoring alternative ` Interpretation 1 0.157 Substantial (in favor of exluding from the model 2 0.129 Substantial (in favor of exluding from the model 3 14.3 Strong (in favor of keeping in the model) 4 10.6 Strong (in favor of keeping in the model) The summary , unless the argument interpretable_scale is set to FALSE , provides estimates in terms of the o dds ratios. F rom this, we see that there is an estimated 52.8% (95% CI: 22.7% - 71.2%) reduction in odds of PEP due to the treatment ( rx1_indomethacin ). The p osterior probabilit y that this effect is b eneficial (i.e., a reduction in PEP) is v ery high (the probabilit y of direction equals 0.999), and this is corrob orated by strong evidence from the Ba y es factor. Moreo v er, this statistic al ly significant effect is also a pr actic al ly significant result, demonstrated by the posterior probabilit y that the treatment effect falls in the R OPE is 0.006. These effects can easily b e visualized. Figure 11 sho ws the credible bands for the four co- v ariates; note that these visualizations are giv en on the 0 to 1 scale for easier interpretation. When bay esics creates suc h credible or prediction bands, other co v ariates must b e fixed. Users can fix these other v alues themselv es via the exemplar_covariates , but otherwise bay esics finds and uses the medoid observ ation from the data used to fit the model, ensuring that the co v ariate combination is a go o d representation of the data while ensuring it is not p o orly extrap olated as might be the case w ere the means of the co v ariate v alues to be used. 6. Discussion The ba y esics R package p ro vides a unified framework for p erforming and in terpreting results from Bay esian analyses. It fo cuses on the most commonly used approaches (hence the package 22 ba y esics: Core Statistical Methods via Bay esian Inference in R R> plot(pep_fit, + type = "cred band") 0.00 0.25 0.50 0.75 1.00 20 40 60 80 age as.numeric(outcome) Credible band f or age 0.0 0.1 0.2 1_female 2_male gender outcome Credible band f or gender 0.00 0.25 0.50 0.75 1.00 1 2 3 4 5 risk score as.numeric(outcome) Credible band f or risk 0.00 0.05 0.10 0.15 0.20 0_placebo 1_indomethacin rx outcome Credible band f or rx Figure 11: Credible bands for analyzing the PEP outcome using a GLM based on the Bernoulli distribution. Journal of Statistical Softw are 23 name’s homophone with “basics”), and aims to provide the most imp ortant Ba y esian inferen- tial quantities, namely p oint estimates, interv al estimates, probability of direction, probability that an estimand falls in the R OPE, and when applicable Bay es factors. ba y esics aims to eliminate the need for algorithmic kno wledge so that the user can focus on inference. This is ac hiev ed by using closed-form solutions when p ossible, and when p oste- rior sampling is re quired automatically selecting the num b er of samples required for accurate Mon te Carlo estimates. This is a concern that the authors feel is widely ignored or un- derestimated b y practitioners, as algorithmic quantities measuring MCMC conv ergence and accuracy of p oin t estimates provided in alternative softw are ma y mislead practitioners in to a false sense of accuracy in in terv al estimates. While algorithmic assessment is not applicable to bay esics , mo del assessmen t is still essential. ba y esics attempts to provide easy access to diagnostic plots, and should model assumptions not hold, bay esics pro vides sev eral non-parametric or semi-parametric alternativ es. It is the authors’ hope that bay esics will serv e as a useful and accurate soft w are implemen- tation of highly used Ba y esian analyses. Ho w ev er, w e an ticipate that there will be coun tless bugs to correct and b eneficial features not y et implemented, and we hum bly ask the statistical comm unit y to rep ort suc h issues at https://github.com/dksewell/bayesics/issues . References Barc h DH, Chechile RA (2023). DFBA: Distribution-F r e e Bayesian A nalysis . doi: 10.32614/CRAN.package.DFBA . R package v ersion 0.1.0, URL https://CRAN.R- project. org/package=DFBA . Bon JJ, Bretherton A, Buc hhorn K, Cram b S, Drov andi C, Hassan C, Jenner AL, Mayfield HJ, McGree JM, Mengersen K, Price A, Salomone R, Santos- F ernandez E, V ercelloni J, W ang X (2023). “Being Ba y esian in the 2020s: opp ortunities and challenges in the practice of mo dern applied Ba y esian statis- tics. ” Philosophic al T r ansactions of the R oyal So ciety A: Mathematic al, Physi- c al and Engine ering Scienc es , 381 (2247), 20220156. ISSN 1364-503X. doi:10. 1098/rsta.2022.0156 . https://royalsocietypublishing.org/rsta/article- pdf/ doi/10.1098/rsta.2022.0156/1327203/rsta.2022.0156.pdf , URL https://doi.org/ 10.1098/rsta.2022.0156 . Bürkner PC (2017). “brms: An R P ac kage for Ba y esian Multilev el Mo dels Using Stan. ” Journal of Statistic al Softwar e , 80 (1), 1–28. doi:10.18637/jss.v080.i01 . Carp en ter B, Gelman A, Hoffman MD, Lee D, Go o drich B, Betancourt M, Brubak er M, Guo J, Li P , Riddell A (2017). “Stan: A Probabilistic Programming Language. ” Journal of Statistic al Softwar e , 76 (1), 1–32. doi:10.18637/jss.v076.i01 . URL https://www. jstatsoft.org/index.php/jss/article/view/v076i01 . Colquhoun D (2017). “The reproducibility of research and the misin terpretation of p- v alues. ” R oyal So ciety Op en Scienc e , 4 (12), 171085. ISSN 2054-5703. doi:10.1098/rsos. 171085 . https://royalsocietypublishing.org/rsos/article- pdf/doi/10.1098/ rsos.171085/953838/rsos.171085.pdf , URL https://doi.org/10.1098/rsos.171085 . 24 ba y esics: Core Statistical Methods via Bay esian Inference in R de V alpine P , T urek D, Paciorek C, Anderson-Bergman C, T emple Lang D, Bodik R (2017). “Programming with models: writing statistical algorithms for general mo del structures with NIMBLE. ” Journal of Computational and Gr aphic al Statistics , 26 , 403–413. doi: 10.1080/10618600.2016.1172487 . Doss CR, Flegal JM, Jones GL, Neath RC (2014). “Marko v chain Monte Carlo estimation of quan tiles. ” Ele ctr onic Journal of Statistics , 8 (2), 2448 – 2478. doi:10.1214/14- EJS957 . URL https://doi.org/10.1214/14- EJS957 . Elm unzer BJ, Sc heiman JM, Lehman GA, Chak A, Mosler P , Higgins PD, Hayw ard RA, Romagn uolo J, Elta GH, Sherman S, W aljee AK, Repaka A, Atkinson MR, Cote GA, Kw on RS, McHenry L, Piraka CR, W amsteker EJ, W atkins JL, K orsnes SJ, Schmidt SE, T urner SM, Nic holson S, F ogel EL (2012). “A Randomized T rial of Rectal Indomethacin to Prev en t P ost-ER CP P ancreatitis. ” New England Journal of Me dicine , 366 (15), 1414–1422. doi:10. 1056/NEJMoa1111103 . https://www.nejm.org/doi/pdf/10.1056/NEJMoa1111103 , URL https://www.nejm.org/doi/full/10.1056/NEJMoa1111103 . Gelman A, Carlin JB, Stern HS, Dunson DB, V ehtari A, R ubin DB (2004). Bayesian Data Analysis . Third edition. Chapman & Hall/CR C, Bo ca Raton, USA. Gelman A, Su YS (2024). arm: Data A nalysis U sing R e gr ession and Multilevel/Hier ar chic al Mo dels . doi:10.32614/CRAN.package.arm . R pac kage v ersion 1.14-4, URL https:// CRAN.R- project.org/package=arm . Gew ek e J (1991). “Ev aluating the accuracy of sampling-based approaches to the calculation of posterior momen ts. ” Staff R ep ort 148 , F ederal Reserv e Bank of Minneap olis. URL https://EconPapers.repec.org/RePEc:fip:fedmsr:148 . Golding N (2019). “greta: simple and scalable statistical mo delling in R. ” Journal of Op en Sour c e Softwar e , 4 (40), 1601. doi:10.21105/joss.01601 . Go o drich B, Gabry J, Ali I, Brilleman S (2025). “rstanarm: Bay esian applied regression mo deling via Stan. ” R pac kage v ersion 2.32.2, URL https://mc- stan.org/rstanarm/ . Greenland S, Senn SJ, Rothman KJ, Carlin JB, P o ole C, Go o dman SN, Altman DG (2016). “Statistical tests, P v alues, confidence in terv als, and p ow er: a guide to misin terpreta- tions. ” Eur op e an Journal of Epidemiolo gy , 31 (4), 337–350. ISSN 1573-7284. doi: 10.1007/s10654- 016- 0149- 3 . URL https://doi.org/10.1007/s10654- 016- 0149- 3 . Hagan A, W est M (2010). The Oxfor d Handb o ok of Applie d Bayesian A nalysis . Oxford Hand- b o oks. OUP Oxford. ISBN 9780191582820. URL https://books.google.com/books?id= pnAWEAAAQBAJ . Imai K, Keele L, Tingley D (2010). “A general approac h to causal mediation analysis. ” Psycholo gic al Metho ds , 15 (4), 309–334. Krusc hk e JK (2018). “Rejecting or Accepting P arameter V alues in Bay esian Estimation. ” A dvanc es in Metho ds and Pr actic es in Psycholo gic al Scienc e , 1 (2), 270–280. doi:10.1177/ 2515245918771304 . https://doi.org/10.1177/2515245918771304 , URL https://doi. org/10.1177/2515245918771304 . Journal of Statistical Softw are 25 Lyddon SP , Holmes CC, W alker SG (2019). “General Bay esian up dating and the loss- lik eliho o d b o otstrap. ” Biometrika , 106 (2), 465–478. ISSN 0006-3444. doi:10. 1093/biomet/asz006 . https://academic.oup.com/biomet/article- pdf/106/2/465/ 28575480/asz006.pdf , URL https://doi.org/10.1093/biomet/asz006 . Mak o wski D, Ben-Shachar MS, Lüdec k e D (2019). “bay estestR: Describing Effects and their Uncertain t y , Existence and Significance within the Bay esian F ramew ork. ” Journal of Op en Sour c e Softwar e , 4 (40), 1541. doi:10.21105/joss.01541 . URL https://joss.theoj. org/papers/10.21105/joss.01541 . Morey RD, Rouder JN (2024). BayesF actor: Computation of Bayes F actors for Common Designs . doi:10.32614/CRAN.package.BayesFactor . R pac kage v ersion 0.9.12-4.7, URL https://CRAN.R- project.org/package=BayesFactor . Qing Y, Thall PF, Y uan Y (2023). “A Bay esian piecewise exp onen tial phase I I design for monitoring a time-to-ev en t endp oint. ” Pharmac eutic al Statistics , 22 (1), 34–44. doi:https: //doi.org/10.1002/pst.2256 . https://onlinelibrary.wiley.com/doi/pdf/10.1002/ pst.2256 , URL https://onlinelibrary.wiley.com/doi/abs/10.1002/pst.2256 . Rix A, Kleinsasser M, Song Y (2025). b ama: High Dimensional Bayesian Me diation A nalysis . doi:10.32614/CRAN.package.bama . R pac kage version 1.3.1, URL https: //CRAN.R- project.org/package=bama . Rosner GL, Laud PW, Johnson WO (2021). Bayesian Thinking in Biostatistics . Chapman and Hall/CRC. Salimans T, Kno wles D A (2013). “Fixed-F orm V ariational P osterior Approximation through Sto c hastic Linear Regression. ” Bayesian A nalysis , 8 (4), 837 – 882. doi:10.1214/13- BA858 . URL https://doi.org/10.1214/13- BA858 . Stan Dev elopmen t T eam (2025). “RStan: the R interface to Stan. ” R pac kage version 2.32.7, URL https://mc- stan.org/ . Statisticat, LLC (2021). L aplac esDemon: Complete Envir onment for Bayesian Infer enc e . R pac kage v ersion 16.1.6, URL https://web.archive.org/web/20150206004624/http: //www.bayesian- inference.com/software . W asserstein RL, Lazar NA (2016). “The ASA Statemen t on p-V alues: Con text, Process, and Purpose. ” The A meric an Statistician , 70 (2), 129–133. doi:10.1080/00031305.2016. 1154108 . https://doi.org/10.1080/00031305.2016.1154108 , URL https://doi.org/ 10.1080/00031305.2016.1154108 . Štrum b elj E, Bouc hard-Côté A, Corander J, Gelman A, Rue H, Murra y L, P esonen H, Plummer M, V eh tari A (2024). “P ast, Present and F uture of Softw are for Ba y esian In- ference. ” Statistic al Scienc e , 39 (1), 46 – 61. doi:10.1214/23- STS907 . URL https: //doi.org/10.1214/23- STS907 . 26 ba y esics: Core Statistical Methods via Bay esian Inference in R Affiliation: Daniel K. Sewell Departmen t of Biostatistics Univ ersit y of Iow a 145 n. Riverside Dr., Io w a City , IA 52242, USA E-mail: daniel-sewell@uiowa.edu URL: https://sewell.lab.uiowa.edu/ Journal of Statistical Software http://www.jstatsoft.org/ published by the F oundation for Open Access Statistics http://www.foastat.org/ MMMMMM YYYY, V olume VV, Issue I I Submitte d: yyyy-mm-dd doi:10.18637/jss.v000.i00 A c c epte d: yyyy-mm-dd