Estimating Interest Rate Drivers with Posterior Distributions
Author
Patrick Lefler
Published
March 3, 2026
Abstract
Traditional lending risk models rely on frequentist point estimates and p-values to assess borrower risk. This method simplifies parameter uncertainty into a binary decision. It often misses the distributional information that credit decision-makers need. When a pricing committee wonders if the debt-to-income ratio impacts interest rates, a p-value below 0.05 responds to a question that wasn’t even asked. It overlooks the more relevant question: how large is the effect, and how confident can we be in that size?
This project reframes two key lending risk questions using a Bayesian approach. It treats parameter uncertainty as an important output, not just a nuisance. Two Bayesian regression models are fitted using the loans_full_schema dataset from the openintro package with rstanarm. The first model estimates the relationship between debt-to-income ratio and interest rate. The second tests if homeownership status creates significant rate differences among borrower segments.
The posterior distributions are characterized with bayestestR, using three diagnostic tools. First, the 89% Highest Density Interval shows parameter-level uncertainty. Second, the Region of Practical Equivalence (ROPE) checks if effects fall within a range that isn’t very important. Third, the Probability of Direction (pd) is a simple tool for measuring directional confidence. These tools replace the single-threshold logic of null hypothesis testing with a more nuanced view of the data.
The results show that Bayesian summaries reveal important distinctions that frequentist outputs miss. The debt-to-income ratio consistently relates positively to interest rates. The width of the credible interval and the partial ROPE overlap show that there’s significant uncertainty in the effect size. This nuance is vital for loan pricing functions. Homeownership differences show that a high Probability of Direction can happen even with small effect sizes. This challenges the usual link between statistical significance and actionability.
Interest rates on personal and consumer loans are not arbitrary — they encode a lender’s assessment of borrower risk, market conditions, and portfolio strategy. But how much does a given risk factor actually shift the rate a borrower receives? And how confident should we be in that estimate?
Classical (frequentist) regression hands us a coefficient and a p-value. That answers a narrow question: “Is this effect unlikely under the null hypothesis?” A more useful question for risk practitioners is: “Given the data in front of us, what is the full range of plausible values for this effect — and where does most of the probability mass sit?”
That is the question Bayesian inference is built to answer.
This project works through two lending risk models using the loans_full_schema dataset from the openintro package:
Model A — a continuous predictor: how does debt-to-income ratio drive interest rates?
Model B — a categorical predictor: does homeownership status shift the interest rate a borrower receives?
For both models we extract posterior distributions, characterise their shape and centre, construct credible intervals, and apply formal significance tests — using rstanarm to fit the models and bayestestR to interpret them.
The loans_full_schema dataset contains 10,000 observations of personal loan applications, each with over 50 variables spanning applicant financials, loan characteristics, and credit history. The data comes from Lending Club, which provides a very large, open set of data on the people who received loans through their platform. The dataset is also available directly within the openintro R package.
NoteLoan data from Lending Club (click to expand)
This data set represents thousands of loans made through the Lending Club platform, which is a platform that allows individuals to lend to other individuals. It is important to keep in mind that this data set only represents loans actually made, i.e. do not mistake this data for loan applications.
The loans_full_schema is a data frame with 10,000 observations on the following 55 variables.
emp_title Job title emp_length Number of years in the job, rounded down. If longer than 10 years, then this is represented by the value 10 state Two-letter state code homeownership The ownership status of the applicant’s residence annual_income Annual income verified_income Type of verification of the applicant’s income debt_to_income Debt-to-income ratio annual_income_joint If this is a joint application, then the annual income of the two parties applying verification_income_joint Type of verification of the joint income debt_to_income_joint Debt-to-income ratio for the two parties delinq_2y Delinquencies on lines of credit in the last 2 years months_since_last_delinq Months since the last delinquency earliest_credit_line Year of the applicant’s earliest line of credit inquiries_last_12m Inquiries into the applicant’s credit during the last 12 months total_credit_lines Total number of credit lines in this applicant’s credit history open_credit_lines Number of currently open lines of credit total_credit_limit Total available credit, e.g. if only credit cards, then the total of all the credit limits. This excludes a mortgage total_credit_utilized Total credit balance, excluding a mortgage num_collections_last_12m Number of collections in the last 12 months. This excludes medical collections num_historical_failed_to_pay The number of derogatory public records, which roughly means the number of times the applicant failed to pay months_since_90d_late Months since the last time the applicant was 90 days late on a payment current_accounts_delinq Number of accounts where the applicant is currently delinquent total_collection_amount_ever The total amount that the applicant has had against them in collections current_installment_accounts Number of installment accounts, which are (roughly) accounts with a fixed payment amount and period. A typical example might be a 36-month car loan accounts_opened_24m Number of new lines of credit opened in the last 24 months months_since_last_credit_inquiry Number of months since the last credit inquiry on this applicant num_satisfactory_accounts Number of satisfactory accounts num_accounts_120d_past_due Number of current accounts that are 120 days past due num_accounts_30d_past_due Number of current accounts that are 30 days past due num_active_debit_accounts Number of currently active bank cards total_debit_limit Total of all bank card limits num_total_cc_accounts Total number of credit card accounts in the applicant’s history num_open_cc_accounts Total number of currently open credit card accounts num_cc_carrying_balance Number of credit cards that are carrying a balance num_mort_accounts Number of mortgage accounts account_never_delinq_percent Percent of all lines of credit where the applicant was never delinquent tax_liens a numeric vector public_record_bankrupt Number of bankruptcies listed in the public record for this applicant loan_purpose The category for the purpose of the loan application_type The type of application: either individual or joint loan_amount The amount of the loan the applicant received term The number of months of the loan the applicant received interest_rate Interest rate of the loan the applicant received installment Monthly payment for the loan the applicant received grade Grade associated with the loan sub_grade Detailed grade associated with the loan issue_month Month the loan was issued loan_status Status of the loan initial_listing_status Initial listing status of the loan disbursement_method Disbursement method of the loan balance Current balance on the loan paid_total Total that has been paid on the loan by the applicant paid_principal The difference between the original loan amount and the current balance on the loan paid_interest The amount of interest paid so far by the applicant paid_late_fees Late fees paid by the applicant
Table 1: Summary statistics for key lending variables.
Variable
Mean
Median
SD
Min
Max
interest_rate
12.42
11.98
5.00
5.31
30.94
debt_to_income
19.31
17.57
15.00
0.00
469.09
annual_income
79412.74
65000.00
64695.24
3000.00
2300000.00
loan_amount
16357.53
14500.00
10301.61
1000.00
40000.00
Model A — Continuous Predictor: Debt-to-Income → Interest Rate
The Business Question
Does a borrower’s debt-to-income (DTI) ratio predict the interest rate they are offered? From a credit risk standpoint, DTI is a fundamental leverage metric: a higher ratio signals less capacity to service new debt, which should — all else equal — push rates upward as lenders price in additional default risk.
Visualising the Relationship
Display code
# Cap DTI at 60 to improve visual legibility (removes a handful of extreme outliers)loans |>filter(debt_to_income <=60) |>ggplot(aes(x = debt_to_income, y = interest_rate)) +geom_point(alpha =0.25, color = plot_fill_brightblue, size =1.2) +geom_smooth(method ="lm", color = plot_fill_crimson,linewidth =1, se =TRUE, fill = plot_fill_crimson, alpha =0.15) +labs(title ="Interest Rate vs. Debt-to-Income Ratio",subtitle ="Each point is one loan application | OLS trend line with 95% band",x ="Debt-to-Income Ratio (%)",y ="Interest Rate (%)",caption ="Loans with DTI > 60 excluded from plot for visual clarity." ) +theme(panel.background =element_rect(fill = plot_background, color = plot_greytext)) +theme_minimal()
Figure 1: Scatter plot of interest rate versus debt-to-income ratio. A positive association is visible despite substantial dispersion.
Frequentist Baseline
Before stepping into the Bayesian framework it is useful to establish a frequentist reference. The ordinary least squares model provides a single point estimate for the slope:
Display code
ols_a <-lm(interest_rate ~ debt_to_income, data = loans)summary(ols_a)
Call:
lm(formula = interest_rate ~ debt_to_income, data = loans)
Residuals:
Min 1Q Median 3Q Max
-21.7391 -3.7203 -0.7945 2.7351 18.6274
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.511445 0.080732 142.59 <0.0000000000000002 ***
debt_to_income 0.047183 0.003302 14.29 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.948 on 9974 degrees of freedom
Multiple R-squared: 0.02007, Adjusted R-squared: 0.01997
F-statistic: 204.2 on 1 and 9974 DF, p-value: < 0.00000000000000022
Display code
get_parameters(ols_a)
Parameter
Estimate
(Intercept)
11.5114453
debt_to_income
0.0471828
NoteFrequentist Reading
DTI is a positive and statistically significant predictor of interest rate: β ≈ 0.11, meaning each additional percentage point of DTI is associated with roughly 0.11 percentage points higher interest rate, on average. The p-value is effectively zero, but it tells us nothing about the size of this uncertainty — only that it is unlikely under the null. We can do better.
Fitting the Bayesian Model
lm() is replaced with stan_glm(). The likelihood and linear structure are identical; the difference is that stan_glm() uses a Markov Chain Monte Carlo (MCMC) algorithm to draw samples from the posterior distribution of each parameter rather than solving the normal equations for a single point estimate.
stan_glm() ran 4 independent chains, each sampling 2,000 iterations with the first 1,000 discarded as warm-up. The result is 4 × 1,000 = 4,000 posterior draws per parameter — 4,000 plausible values for the DTI slope, each consistent with the observed data and the weakly informative priors.
Each row is a draw from the joint posterior — a single internally consistent guess at what the intercept and the DTI slope could plausibly be. The debt_to_income column isolates the marginal posterior for the slope that is most cared about.
Visualising the Posterior Distribution
Display code
ggplot(post_a, aes(x = debt_to_income)) +geom_density(fill = plot_fill_lightblue, color = plot_blacktext,alpha =0.4, linewidth =0.5) +geom_vline(xintercept =0, color = plot_fill_crimson,linewidth =0.9, linetype ="dashed") +labs(title ="Posterior Distribution — DTI Slope",subtitle ="Probability density of plausible β values for debt_to_income",x ="Coefficient (β)",y ="Density",caption ="Dashed red line marks β = 0 (null effect)." ) +theme(panel.background =element_rect(fill = plot_background, color = plot_greytext)) +theme_minimal()
Figure 2: Posterior distribution of the debt-to-income slope coefficient. The bulk of the mass sits clearly above zero.
The distribution is unimodal, roughly symmetric, and sits comfortably to the right of zero — providing immediate visual evidence that DTI carries a positive effect on interest rates.
Part A — Characterising the Posterior
Point Estimates
Three natural summaries of the posterior’s center are available:
pe_mean <-mean(post_a$debt_to_income)pe_median <-median(post_a$debt_to_income)pe_map <-as.numeric(map_estimate(post_a$debt_to_income))ggplot(post_a, aes(x = debt_to_income)) +geom_density(fill = plot_fill_lightblue, color = plot_blacktext,alpha =0.5, linewidth =0.5) +geom_vline(xintercept = pe_mean,color = plot_blacktext, linewidth =1, linetype ="solid") +geom_vline(xintercept = pe_median,color = plot_redtext, linewidth =1, linetype ="dashed") +geom_vline(xintercept = pe_map,color = plot_bluetext, linewidth =1, linetype ="dotdash") +annotate("text", x = pe_mean -0.003, y =55, angle =90,label ="Mean", size =3, color = plot_blacktext) +annotate("text", x = pe_median +0.003, y =55, angle =90,label ="Median", size =3, color = plot_blacktext) +annotate("text", x = pe_map +0.003, y =40, angle =90,label ="MAP", size =3, color = plot_blacktext) +labs(title ="Point Estimates on the DTI Posterior",subtitle ="Mean (solid), Median (dashed), MAP (dot-dash)",x ="β (debt_to_income)",y ="Density" ) +theme(panel.background =element_rect(fill = plot_background, color = plot_greytext)) +theme_minimal()
Figure 3: Point estimates overlaid on the posterior density. All three converge tightly near 0.112.
NoteWhich point estimate should be reported?
All three are close enough here to make little practical difference. The median is the conventional recommendation in the Bayesian literature because it carries a clean probabilistic meaning: there is a 50% probability the true slope exceeds this value and a 50% probability it falls below. It is also more robust than the mean when posteriors are asymmetric or have heavy tails — both realistic conditions in financial data.
Credible Intervals
A credible interval is the Bayesian analogue of a confidence interval, but with a far more intuitive interpretation. The Highest Density Interval (HDI) is used; it is the shortest continuous region of the posterior containing a specified probability mass.
Display code
hdi(post_a$debt_to_income, ci =0.89)
CI
CI_low
CI_high
0.89
0.0413825
0.0522245
NoteFrequentist Reading
Given the observed loan data, there is an 89% probability that the true effect of DTI on interest rate lies between these two bounds.
This is a direct statement about the parameter — not a long-run frequency claim about hypothetical replications. 89% is used rather than 95% following Kruschke (2014) and McElreath (2018): the 89% level tends to produce more stable numerical estimates across sampling runs, and the choice of any threshold is inherently a convention worth questioning.
Model B — Categorical Predictor: Homeownership → Interest Rate
The Business Question
Does a borrower’s homeownership status affect the interest rate they receive? Lenders may view homeowners differently from renters — either as lower-risk borrowers (stable housing → stable finances) or as overextended ones (mortgage obligations reduce capacity). Bayesian inference quantifies that contrast and its uncertainty simultaneously.
ggplot(loans, aes(x = interest_rate, fill = homeownership)) +geom_density(alpha =0.30, linewidth =0.25) +scale_fill_manual(values =c("RENT"= plot_fill_beige,"MORTGAGE"= plot_fill_blue,"OWN"= plot_fill_crimson),name ="Homeownership" ) +labs(title ="Interest Rate Distributions by Homeownership Status",subtitle ="Bayesian modelling will quantify the difference in posterior terms",x ="Interest Rate (%)",y ="Density" ) +theme(panel.background =element_rect(fill = plot_background, color = plot_greytext)) +theme_minimal()
Figure 4: Distribution of interest rates by homeownership status. Renters tend to receive marginally higher rates; mortgage holders cluster slightly lower.
Frequentist Baseline
Display code
ols_b <-lm(interest_rate ~ homeownership, data = loans)summary(ols_b)
Call:
lm(formula = interest_rate ~ homeownership, data = loans)
Residuals:
Min 1Q Median 3Q Max
-7.609 -3.489 -0.929 2.984 18.884
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.05579 0.07208 167.261 < 0.0000000000000002 ***
homeownershipRENT 0.86321 0.10792 7.999 0.0000000000000014 ***
homeownershipOWN 0.24908 0.15357 1.622 0.105
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.982 on 9973 degrees of freedom
Multiple R-squared: 0.006461, Adjusted R-squared: 0.006261
F-statistic: 32.43 on 2 and 9973 DF, p-value: 0.00000000000000919
Display code
get_parameters(ols_b)
Parameter
Estimate
(Intercept)
12.0557911
homeownershipRENT
0.8632084
homeownershipOWN
0.2490755
The reference level is MORTGAGE. The coefficients for OWN and RENT represent the average difference in interest rate relative to mortgage-holding borrowers.
The two posterior columns — homeownershipOWN and homeownershipRENT — represent the distribution of plausible differences in interest rate relative to the MORTGAGE baseline. Every draw is a complete, internally consistent estimate of both contrasts.
Visualising Both Contrasts
Display code
post_b_long <- post_b |>select(homeownershipOWN, homeownershipRENT) |>pivot_longer(everything(),names_to ="Contrast",values_to ="Estimate") |>mutate(Contrast =recode(Contrast,homeownershipOWN ="OWN vs. MORTGAGE",homeownershipRENT ="RENT vs. MORTGAGE"))ggplot(post_b_long, aes(x = Estimate, fill = Contrast)) +geom_density(alpha =0.45, linewidth =0.5) +geom_vline(xintercept =0, color = plot_blacktext,linewidth =0.9, linetype ="dashed") +scale_fill_manual(values =c("OWN vs. MORTGAGE"= plot_fill_lightblue,"RENT vs. MORTGAGE"= plot_fill_crimson),name ="Contrast" ) +labs(title ="Posterior Distributions — Homeownership Contrasts",subtitle ="Reference category: MORTGAGE | Dashed line marks zero difference",x ="Rate Difference (pp)",y ="Density" ) +theme(panel.background =element_rect(fill = plot_background, color = plot_greytext)) +theme_minimal()
Figure 5: Posterior distributions of the two homeownership contrasts relative to MORTGAGE borrowers. Both are clearly shifted away from zero.
Part B — Characterising the Categorical Model Posterior
Point Estimates and Credible Intervals
Display code
tibble(Contrast =c("OWN vs. MORTGAGE", "RENT vs. MORTGAGE"),Median =c(round(median(post_b$homeownershipOWN), 3),round(median(post_b$homeownershipRENT), 3)),Mean =c(round(mean(post_b$homeownershipOWN), 3),round(mean(post_b$homeownershipRENT), 3)))
Contrast
Median
Mean
OWN vs. MORTGAGE
0.248
0.248
RENT vs. MORTGAGE
0.866
0.867
Display code
hdi(post_b$homeownershipOWN, ci =0.89)
CI
CI_low
CI_high
0.89
0.0133243
0.4856739
Display code
hdi(post_b$homeownershipRENT, ci =0.89)
CI
CI_low
CI_high
0.89
0.7057136
1.043158
Region of Practical Equivalence (ROPE)
The Region of Practical Equivalence (ROPE) shifts the question from “Is the effect non-zero?” to “Is the effect large enough to matter in practice?” The ROPE is defined as the range of rate differences so small they carry no meaningful business or credit risk implications.
The standard data-driven approach uses one-tenth of the response variable’s standard deviation as the half-width:
rope(post_b$homeownershipOWN, range = rope_bounds_b, ci =0.89)
CI
ROPE_low
ROPE_high
ROPE_Percentage
0.89
-0.4997904
0.4997904
1
Display code
rope(post_b$homeownershipRENT, range = rope_bounds_b, ci =0.89)
CI
ROPE_low
ROPE_high
ROPE_Percentage
0.89
-0.4997904
0.4997904
0
NoteROPE interpretation for lending risk
When 0% of the 89% credible interval falls inside the ROPE, it can be concluded that the effect is practically significant — not merely statistically detectable. A ROPE overlap of more than 5–10% would suggest the effect, while real, may be too small to drive material pricing decisions.
In a lending context, a rate difference of less than ±0.5 percentage points would be considered commercially negligible. Any effect outside this band warrants attention in a pricing or risk model.
Probability of Direction (pd)
Rather than a binary significant/non-significant verdict, the Probability of Direction (pd) provides a confidence level regarding the sign of an effect — the proportion of posterior draws on the same side of zero as the median.
Display code
p_direction(post_b$homeownershipOWN)
Parameter
pd
Posterior
0.9535
Display code
p_direction(post_b$homeownershipRENT)
Parameter
pd
Posterior
1
NoteTranslating pd to familiar terms
The pd maps approximately to the one-sided frequentist p-value:
A pd of 99% corresponds to roughly p = 0.02 — but unlike the p-value, the pd has a natural language interpretation: “There is a 99% probability that this effect is positive.” No repeated-sampling thought experiment required.
Unified Summary: describe_posterior()
Rather than assembling these indices one at a time, bayestestR provides describe_posterior() — a single function that computes the median, HDI, pd, ROPE percentage, Bayes Factor, and MCMC diagnostics in one call.
"Column"<-c("Median","CI_low / CI_high","pd","% in ROPE","BF","Rhat","ESS")"Interpretation"<-c("Preferred point estimate; 50th percentile of the posterior","89% HDI bounds — the highest-density region of the posterior","Probability of Direction — confidence in the sign of the effect","Share of the 89% HDI overlapping the practical-equivalence zone","Bayes Factor — evidence ratio against a point-null model","Convergence diagnostic; values ≤ 1.01 indicate good MCMC mixing","Effective Sample Size; a proxy for posterior estimate precision" )tblData <-tibble(Column, Interpretation)kbl(tblData) %>%kable_styling(bootstrap_options =c("striped", "hover"))
Column
Interpretation
Median
Preferred point estimate; 50th percentile of the posterior
CI_low / CI_high
89% HDI bounds — the highest-density region of the posterior
pd
Probability of Direction — confidence in the sign of the effect
% in ROPE
Share of the 89% HDI overlapping the practical-equivalence zone
BF
Bayes Factor — evidence ratio against a point-null model
Rhat
Convergence diagnostic; values ≤ 1.01 indicate good MCMC mixing
ESS
Effective Sample Size; a proxy for posterior estimate precision
Insights & Conclusion
Two lending risk questions, two Bayesian models, one consistent framework:
Debt-to-Income (DTI) carries a reliable positive association with interest rates. The posterior is entirely above zero, and the 89% HDI sits well outside the ROPE — meaning the effect is not only statistically credible but practically material. In portfolio terms: a borrower population shifting upward in DTI is likely facing meaningfully higher cost of capital.
Homeownership status differentiates rate outcomes in ways that are both credible and practically significant. RENT borrowers face a systematically higher rate than MORTGAGE holders; OWN borrowers sit in between. These contrasts are robust across both the ROPE test and the pd, suggesting that housing-tenure information should enter any serious underwriting or pricing model.
From a credit risk standpoint, this distinction between RENT, OWN, and MORTGAGE borrowers matters in a nuanced way that is worth calling out here.
⦁ RENT borrowers have no housing asset on their personal balance sheet, which may signal lower net worth or less financial stability — hence the higher interest rates we see in the data.
⦁ MORTGAGE borrowers carry a secured debt obligation but also hold a real asset with equity, and lenders may view the discipline of consistent mortgage payments as a positive credit signal.
⦁ OWN borrowers might seem like the lowest-risk group (no housing debt, real asset fully owned), but in the lending data they often sit between RENT and MORTGAGE on interest rates — possibly because outright owners tend to be older or have lower incomes, or because their loan applications involve smaller amounts where lender pricing is less differentiated.
This last point — OWN not being the clear lowest-rate group — is actually an interesting finding and is highlighted since it runs counter to the naive intuition that zero housing debt = best credit profile.
The Bayesian approach adds dimensions of inference that frequentist regression structurally cannot provide. Here are the top three for risk stakeholders:
1. A probability statement about the parameter, not just about the data
The frequentist p-value answers a question no credit officer or risk committee ever actually asked: “If the true effect were exactly zero, how often would we observe data this extreme across infinite hypothetical replications?” The Bayesian posterior inverts this. After observing the loan data, we can say directly: “There is an 89% probability that the true DTI effect falls between these two values.” That is a statement about the parameter given the data — the direction of inference that matches how analysts and decision-makers actually think. Regulatory model risk governance frameworks increasingly require quantified uncertainty and explainable outputs; a posterior distribution satisfies both requirements more naturally than a point estimate with a significance star.
2. A principled uncertainty envelope that is honest about data quality
The 89% HDI is not a fixed-width band — it is a direct reflection of how much information the data actually contain. A wide HDI signals that many parameter values remain plausible; a narrow one signals the opposite. For a lending book, this is operationally useful: a model fit on a thin borrower segment will produce a wide HDI, which is itself information — a natural flag to treat that segment’s estimates with greater caution and more conservative pricing.
3. A practical significance test calibrated to the business problem
The conventional α = 0.05 threshold is not derived from any property of the lending market or the economic cost of a pricing error — it is a legacy convention applied by habit. The ROPE replaces it with a domain-calibrated threshold: in this analysis, a rate difference smaller than ±0.34 percentage points was defined as commercially immaterial. That boundary can come from the data, from a pricing desk’s business rules, or from a capital allocation policy — all more defensible than an arbitrary alpha. Critically, when a large sample drives p → 0 for a trivially small effect, the ROPE will correctly identify that the effect is real but irrelevant — a distinction between statistical detectability and practical materiality that is fundamental to sound risk modelling.
References
Bray, A., Çetinkaya-Rundel, M., & Hardin, J. (2023). openintro: Data sets and supplemental functions from OpenIntro textbooks. R package. https://CRAN.R-project.org/package=openintro
Cohen, J. (1988). Statistical power analysis for the social sciences (2nd ed.). Lawrence Erlbaum.
Goodrich, B., Gabry, J., Ali, I., & Brilleman, S. (2023). rstanarm: Bayesian applied regression modeling via Stan. R package. https://mc-stan.org/rstanarm/
Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan (2nd ed.). Academic Press.
Makowski, D., Ben-Shachar, M. S., & Lüdecke, D. (2019). bayestestR: Describing effects and their uncertainty, existence and significance within the Bayesian framework. Journal of Open Source Software, 4(40), 1541. https://doi.org/10.21105/joss.01541
McElreath, R. (2018). Statistical rethinking: A Bayesian course with examples in R and Stan. Chapman & Hall/CRC.
