# Chapter 04: Tailoring priors — leveraging the Prior_Setup function

#### Kjell Nygren

#### 2026-06-21

## 1. Introductory Discussion

In this chapter, we explore how the Prior_Setup function can be used to specify multivariate normal priors for both Bayesian linear models and the generalized linear models (GLMs) introduced in later chapters. In the previous chapter, we covered the minimal arguments required to invoke the function: a model formula and a family (with an optional link function). These define the structure of the model to which the prior will be applied—namely, the variables being modeled, the assumed distribution of the response, and the link function connecting the linear predictor to the data. Multivariate normal and g-prior structure connects to standard GLM treatments (McCullagh and Nelder 1989; Venables and Ripley 2002; Zellner 1986). Prior–data conflict diagnostics (when used elsewhere in the package) relate to (Evans and Moshonov 2006). Envelope sampling builds on (Nygren and Nygren 2006).

For full technical derivations of the priors and calibrated quantities returned by `Prior_Setup()` (including Gaussian Normal–Gamma calibration), see Chapter A12 (https://knygren.r-universe.dev/articles/glmbayes/Chapter-A12.html).

In many cases, it is also helpful to provide a data frame to Prior_Setup, allowing the function to correctly identify and extract the relevant variables from the model environment.

Now, let’s examine the remaining key arguments that enable the specification of multivariate normal priors. We begin by reviewing the function signature and a representative call—similar to the one used in the previous chapter, but now with visibility into the default settings:

    ps = Prior_Setup(
      weight ~ group, data = dat, family = gaussian(), pwt = 0.001,
      n_prior = NULL, sd = NULL,
      intercept_source = c("null_model", "full_model"),
      effects_source = c("null_effects", "full_model"),
      mu = NULL
    )

    #mu       = ps$mu
    #V        = ps$Sigma
    #disp_ML  = ps$dispersion

As shown above, we introduce six additional arguments beyond the core model specification:

- pwt, n_prior, and sd: These control the strength of the prior—i.e., how much information the prior contributes relative to the likelihood.
- intercept_source, effects_source, and mu: These govern the location or centering of the prior distribution.

In the next section, we discuss the default settings and their implications for model interpretation. Section 3 focuses on the parameters that control prior strength, while Section 4 addresses how the prior is centered.

# 2. The default prior specification

When you call `Prior_Setup()` with a single `pwt` value, the prior on the coefficient vector
\(\beta = (\beta_0,\beta_1,\dots,\beta_{k})^\top\)
is set up to be

\[ \beta \sim \mathcal{N}\bigl(\mu_{Prior},\;\Sigma_{Prior}\bigr), \]

where \(\mu_{Prior}\) is the derived prior mean and \(\Sigma_{Prior}\) the derived prior covariance.

In this section,we will frst focus on how the `pwt` argument is used to set \(\Sigma_{Prior}\) and how some of the other defaults help set \(\mu_{Prior}\). We also take a look at how these together help determine the Posterior mean in Bayesian linear models.

## 2.1 Identifiability and Full–Rank Likelihood Covariance

Before constructing the prior, `Prior_Setup()` fits the full GLM and extracts
the maximum likelihood estimate \(\widehat\beta_{\mathrm{MLE}}\) and the variance–covariance of the maximum‐likelihood estimator. \[ V_{MLE} \;=\;\mathrm{vcov}(\widehat\beta_{\mathrm{MLE}}). \]

The matrix \(V_{MLE}\) is checked to verify that it is a symmetric positive–definite (full rank). If \(V_{MLE}\) is rank‐deficient, the classical GLM lacks identifiability and `Prior_Setup()` aborts with an error. While it is possible to estimate Bayesian models even if the classical GLM lacks identifiability, it is not possible to specificy a g-prior for such a likelihood function and we urge the user to use caution with such models.

## 2.2 Canonical Precision Definitions and Scaling

For enhanced understanding and to simplify derivations, let us now define the Likelihood precision, the Prior Precision, and the Posterior Precision as follows:

\[ P_{MLE} = V_{MLE}^{-1}, \quad P_{Prior} = \Sigma^{-1}, \quad P_{Post} = P_0 + P_{Prior}. \]

For `length(pwt)=1`, the `Prior_Setup()` uses the \(pwt\) arguments to set the prior variance-covariance matrix as an inflated version of the likelihood variance–covariance matrix

\[ \Sigma_{Prior} = \frac{1 - pwt}{pwt}\;V_{MLE}, \quad \Longrightarrow P_{Prior} = \frac{pwt}{1 - pwt}\;P_{MLE}, \]

and therefore

\[ P_{\mathrm{Post}} = \Bigl(1 + \frac{pwt}{1 - pwt}\Bigr)\,P_0 = \frac{1}{1 - pwt}\;P_0. \]

Thus both \(P_{\mathrm{Prior}}\) and \(P_{\mathrm{Post}}\) are simple scalar multiples of \(P_0\). The scaling factor being applied to set \(\Sigma_{Prior}\) is actually a Zellner g-prior:

\[ g \;=\;\frac{1 - \text{pwt}}{\text{pwt}}. \]

It is worth noting that the default for `Pwt`is 0.05 which means that g=19 and the prior Variance is thus 19 times the likelihood variance (and the standard deviation is \(\sqrt{19}\) times larger than the likelihood standard deviations).

## 2.3 Posterior Mean via Precision Matrices

In the Gaussian‐data case, the posterior is

\[ \beta\mid y \;\sim\; \mathcal{N}\bigl(\mu_{\mathrm{Post}},\,\Sigma_{\mathrm{Post}}\bigr), \]

with

\[ \Sigma_{\mathrm{Post}} = P_{\mathrm{Post}}^{-1}, \quad \mu_{\mathrm{Post}} = P_{\mathrm{Post}}^{-1} \bigl(P_0\,\widehat\beta_{\mathrm{MLE}} \;+\;P_{\mathrm{Prior}}\,\mu\bigr). \]

Plugging in the scalar‐multiple forms:

\[ P_{\mathrm{Post}}^{-1} = (1 - pwt)\,P_{MLE}^{-1}, \quad P_{\mathrm{Prior}} = \frac{pwt}{1 - pwt}\,P_{MLE}, \]

we recover the familiar weighted‐average form:

\[ \mu_{\mathrm{Post}} = (1 - pwt)\,\widehat\beta_{\mathrm{MLE}} \;+\;pwt\,\mu_{Prior}. \]

Hence we can see that the pwt argument determine the relative weight on the prior in the formula for the posterior mean in these models. While technically, we don’t have the same closed form solution for other models, the pwt still plays a role in how close the posterior will be to the prior mean.

Since the default value for pwt is 0.01, we have a 1% weight on the prior for the default and the prior (when different from the Maximum likelihood estimate) hence pulls the posterior only minimally away from the maximum likelihood estimate.

## 2.4 The impact of the “null_model” and “null_effects” defaults for the interpretation of model diagnostics

In addition to the pwt arguments with its default, we also see from above that the default for the `intercept_source` is the “null_model” and the default for the `effects_source` is “null_effects”. The former means that the default prior for the intercept defaults to the estimate from the intercept only model while the latter means that the default prior for all effects variables is set to 0. This is in many ways consistent with the reference model utilized in most implementations for classical models. In particular, in classical models:

-The R^2 is defined relative to the null model (i.e., the R^2 is always zero for the null model) -The default F-test in a classical model is usally relative to the null model -The t-statistics usually are designed to measure how far estimates are from 0 and this is the default for how p-values are calculated.

It is worthwhile to now take a look at what this default specification means for the posterior mean when `pwt` is a scalar and the posterior mean is simplified as in the earlier section. In that case, we have:

- \(\mu_{\mathrm{Post},1} = pwt \cdot \bar{y} + (1 - pwt) \cdot \widehat{\beta}_{\mathrm{MLE},\mathrm{Intercept}}\)
- \(\mu_{\mathrm{Post},\mathrm{Effects}} = (1 - pwt) \cdot \widehat{\beta}_{\mathrm{MLE},\mathrm{Effects}}\)

Let us now re-examine the the Plant data from the previous chapter.

    lm.D9 <- lm(weight ~ group)
    sumlm<-summary(lm.D9)
    sumlm$coefficients
    #>             Estimate Std. Error  t value     Pr(>|t|)
    #> (Intercept)    5.032  0.2202177 22.85012 9.547128e-15
    #> groupTrt      -0.371  0.3114349 -1.19126 2.490232e-01

The groupTrt estimate from the classical model is -0.371 with a two sided p-value of 0.2490 (which implies that the one sided p-value is 0.1245).

    lmb.D9=lmb(weight ~ group,dNormal(mu=ps$mu,Sigma=ps$Sigma,dispersion=ps$dispersion),n=1000)
    sumlmb<-summary(lmb.D9)
    sumlmb$coefficients
    #>             Post.Mode  Post.Mean   Post.Sd    MC Error Pr(Null_tail)   SE(tail)
    #> (Intercept)  5.031814  5.0311363 0.2172248 0.006869250         0.202 0.01269630
    #> groupTrt    -0.370629 -0.3641917 0.3092222 0.009778466         0.120 0.01027619
    #>             Pr(Prior_tail)
    #> (Intercept)          0.202
    #> groupTrt             0.120

We can now see that the corresponding output from the Bayesian run in many ways mirror the output from the classical model. In particular, four out of the columns have outputs analogous to corresponding outputs from the classical model. The Post.Mode and Post.Mean, both mirrors the Estimate column while the Post.sd column mirrors the Std. Error column. The Pr(tail), meanwhile is close to half of the value for Pr(>|t|) column in the classical output.

In addition to this, we have a column with MC Error which represents the error connected to the estimated Post. mean columns. This reflects the fact that we are using Monte Carlo simulation (in this case based on n=1000 draws) to estimate the mean. In the gaussian() case, the actual Post.Mean is the same as the Post.Mode so if we increase the number of draws, we are increasingly likely to observe a Post.Mean that is close to the Post.Mode column and the MC Error column shrinks (see below).

Finally, we also have a column SE(tail). This reflects the Standard error connected to the Pr(tail) column. Just like with the MC Error, we can reduce the error associated with this estimate by generating a larger number of draws (again, see below).

    lmb.D9=lmb(weight ~ group,dNormal(mu=ps$mu,Sigma=ps$Sigma,dispersion=ps$dispersion),n=100000)
    sumlmb<-summary(lmb.D9)
    sumlmb$coefficients
    #>             Post.Mode  Post.Mean   Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept)  5.031814  5.0315824 0.2200432 0.0006958378       0.19915
    #> groupTrt    -0.370629 -0.3707955 0.3107300 0.0009826145       0.11643
    #>                SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.001262891        0.19915
    #> groupTrt    0.001014268        0.11643

For the purposes of this vignette, we simply note that under our default assumptions for the intercept and effects, the Pr(tail) for slope parameters in Bayesian linear models tend to be similar to the one-sided p-value in the classical models when we combineex them with our default value for our pwt (i.e., 0.01). While we don’t explore it fully here, there are certain aspects in which these become more fully equivalent in the limit as pwt is pushed further towards zero. A later vignette provides more details on exactly how these tail probabilities relate to their classical counterparts, while also demonstrating some properties as pwt, n_prior, and the sample size approach certain limits.

# 3. Adjusting the strength of priors using pwt, n_prior, and sd

In setting the strength of priors, pwt serves as the default mechanism, blending the prior mean with the maximum likelihood estimates. When pwt is a scalar, it uniformly adjusts this blend across parameters. Alternatively, a vector form of pwt allows for more nuanced adjustments, deviating from a Zellner g-prior. n_prior offers an alternative approach, influencing the strength of priors based on sample size equivalence. Finally, sd directly sets the standard deviation, offering another layer of control over the priors’ influence. Each of these parameters plays a crucial role in tailoring the prior to reflect the analyst’s beliefs and the context of the data.

To illustrate the role of these prior specifications, we will work with a model where the intercept dimension is orthogonal to the predictions as it makes interpretations clearer.

## 3.1 Data Setup and Visualization

    data(mtcars)

    # Create gallons per mile
    mtcars$gpm <- 1 / mtcars$mpg

The mtcars dataset contains fuel efficiency and engine characteristics for 32 vehicles. We begin by preparing the mtcars dataset for modeling. To analyze fuel consumption, we define gallons per mile (gpm) as the reciprocal of miles per gallon (gpm = 1/mpg). This transformation allows us to treat fuel use as a direct cost per unit distance, which aligns more naturally with physical energy expenditure.

    # Mean-center predictors
    #mtcars$c_hp <- scale(mtcars$hp, center = TRUE, scale = FALSE)
    #mtcars$c_drat <- scale(mtcars$drat, center = TRUE, scale = FALSE)
    #mtcars$c_qsec <- scale(mtcars$qsec, center = TRUE, scale = FALSE)

    mtcars$c_wt  <- as.numeric(scale(mtcars$wt, center = TRUE, scale = FALSE))
    mtcars$c_cyl <- as.numeric(scale(mtcars$cyl, center = TRUE, scale = FALSE))

Next, we mean-center the predictors of interest—vehicle weight (wt) and number of cylinders (cyl)—to improve interpretability and reduce collinearity in the regression model. Mean-centering ensures that the intercept corresponds to the expected fuel consumption for a vehicle with average weight and cylinder count.

    # Fit linear model
    fit <- lm(gpm ~ c_wt, data = mtcars)

    # Create scatter plot
    plot(mtcars$c_wt, mtcars$gpm,
         pch = 19, col = "steelblue",
         xlab = "Mean-Centered Weight (c_wt)",
         ylab = "Gallons per Mile (gpm)",
         main = "Fuel Consumption vs. Mean-Centered Weight")

    # Add regression line
    abline(fit, col = "darkred", lwd = 2)

    # Optional: Add confidence band manually
    x_vals <- seq(min(mtcars$c_wt), max(mtcars$c_wt), length.out = 100)
    pred <- predict(fit, newdata = data.frame(c_wt = x_vals), interval = "confidence")

    lines(x_vals, pred[, "lwr"], col = "darkred", lty = 2)
    lines(x_vals, pred[, "upr"], col = "darkred", lty = 2)

To assess the plausibility of a linear relationship, we first plot gpm against mean-centered weight (c_wt). The scatter plot reveals a clear positive trend: heavier vehicles tend to consume more fuel per mile. The fitted regression line aligns well with the data, and the symmetric confidence band suggests that linearity and homoscedasticity are reasonable assumptions. This relationship is also physically motivated, as rolling resistance and energy demand scale approximately linearly with mass.

    # Ensure gpm is defined
    mtcars$gpm <- 1 / mtcars$mpg

    # Create boxplot
    boxplot(gpm ~ factor(cyl), data = mtcars,
            col = "lightblue", border = "darkblue",
            xlab = "Number of Cylinders",
            ylab = "Gallons per Mile (gpm)",
            main = "Fuel Consumption by Number of Cylinders")

    # Add jittered points
    set.seed(123)  # for reproducibility
    points(jitter(as.numeric(factor(mtcars$cyl)), amount = 0.2), mtcars$gpm,
           pch = 19, col = rgb(70/255, 130/255, 180/255, alpha = 0.6))

We then examine fuel consumption across cylinder categories using a boxplot. Treating cyl as a factor, we observe that vehicles with more cylinders tend to have higher median fuel consumption, with 8-cylinder cars showing the greatest fuel use. The jittered points reveal the distribution within each group and highlight the presence of outliers and some overlap acorss groups. This plot supports the inclusion of cyl as a categorical or mean-centered numeric predictor, depending on modeling goals.

Together, these visualizations provide both empirical and theoretical justification for modeling gpm as a function of c_wt and c_cyl. The linear trend with weight is physically motivated, while the categorical differences in cylinder count suggest potential predictive value. These insights set the stage for exploring how prior strength influences inference in Bayesian regression

## 3.2 Using the pwt argument to adjust the Zellner g-prior

To illustrate the impact of prior strength on marginal predictors, we fit a sequence of Bayesian models using gpm as the response and two mean-centered predictors: c_wt and c_cyl.

To ensure reproducibility of posterior estimates and tail probabilities, we set a fixed random seed before simulation. This allows consistent interpretation across models and avoids misleading comparisons due to Monte Carlo variability.

    # Use Prior_Setup to illustrate prior strength
    pscars  <- Prior_Setup(formula = gpm ~  c_wt+c_cyl,family=gaussian(),data = mtcars,pwt=0.01) ## pwt=0.01-->n_prior= 32/99
    pscars2 <- Prior_Setup(formula = gpm ~  c_wt+c_cyl,family=gaussian(),data = mtcars,pwt=0.5)  ## pwt=0.5-->n_prior= 32
    pscars3 <- Prior_Setup(formula = gpm ~  c_wt+c_cyl,family=gaussian(),data = mtcars,pwt=0.55) ## pwt=0.55-->n_prior= 39.111

    pscars
    #>
    #> Call:  Prior_Setup(formula = gpm ~ c_wt + c_cyl, family = gaussian(),
    #>     data = mtcars, pwt = 0.01)
    #>
    #> Setting up a Zellner g-type prior:
    #>   pwt = 0.01
    #>   g   = (1 - pwt)/pwt = 99
    #>
    #> Note: n_prior was computed as (pwt / (1 - pwt)) * n_likelihood:
    #>   n_prior      = 0.3232
    #>   n_likelihood = 32
    #>
    #> Prior Setup Summary
    #> ====================
    #>
    #> Prior Estimates with 95% Confidence Intervals
    #>             Prior.Mean Prior.SD  CI.Lower CI.Upper  pwt
    #> (Intercept)   0.054227 0.012701  0.029335 0.079120 0.01
    #> c_wt          0.000000 0.021180 -0.041512 0.041512 0.01
    #> c_cyl         0.000000 0.011604 -0.022743 0.022743 0.01
    #>
    #> Prior Correlation Matrix
    #>             (Intercept)    c_wt   c_cyl
    #> (Intercept)           1  0.0000  0.0000
    #> c_wt                  0  1.0000 -0.7825
    #> c_cyl                 0 -0.7825  1.0000
    #>
    #> Conditional Dispersion (Gaussian family):  1e-04
    #>
    #> Gamma Prior on Residual Precision:
    #>   shape = 0.6616
    #>   rate = 6.05648e-05
    #>   Expected precision (inverse variance) = 1.09241e+04 which implies 1/Expected precision  = 9.15408e-05
    #>
    #>   Applicable to gaussian models with compound pfamilies (e.g., dNormal_Gamma, dIndependent_Normal_Gamma),
    #>   as well as for Gamma regression, quasipoisson, and quasibinomial models.

With pwt = 0.01, the prior exerts minimal influence, resulting in wide prior standard deviations and broad credible intervals. For example, the prior SD for c_wt is approximately 0.0207, yielding a 95% credible interval of ±0.0406 around the prior mean. This reflects substantial uncertainty and allows the posterior to be largely driven by the data. The intercept’s prior interval spans nearly 0.05 units, indicating weak regularization and high flexibility in model fitting

    pscars2
    #>
    #> Call:  Prior_Setup(formula = gpm ~ c_wt + c_cyl, family = gaussian(),
    #>     data = mtcars, pwt = 0.5)
    #>
    #> Setting up a Zellner g-type prior:
    #>   pwt = 0.5
    #>   g   = (1 - pwt)/pwt = 1
    #>
    #> Note: n_prior was computed as (pwt / (1 - pwt)) * n_likelihood:
    #>   n_prior      = 32
    #>   n_likelihood = 32
    #>
    #> Prior Setup Summary
    #> ====================
    #>
    #> Prior Estimates with 95% Confidence Intervals
    #>             Prior.Mean Prior.SD  CI.Lower CI.Upper pwt
    #> (Intercept)   0.054227 0.002298  0.049722 0.058732 0.5
    #> c_wt          0.000000 0.003833 -0.007513 0.007513 0.5
    #> c_cyl         0.000000 0.002100 -0.004116 0.004116 0.5
    #>
    #> Prior Correlation Matrix
    #>             (Intercept)    c_wt   c_cyl
    #> (Intercept)           1  0.0000  0.0000
    #> c_wt                  0  1.0000 -0.7825
    #> c_cyl                 0 -0.7825  1.0000
    #>
    #> Conditional Dispersion (Gaussian family):  2e-04
    #>
    #> Gamma Prior on Residual Precision:
    #>   shape = 16.5
    #>   rate = 0.0029
    #>   Expected precision (inverse variance) = 5.74117e+03 which implies 1/Expected precision  = 1.74181e-04
    #>
    #>   Applicable to gaussian models with compound pfamilies (e.g., dNormal_Gamma, dIndependent_Normal_Gamma),
    #>   as well as for Gamma regression, quasipoisson, and quasibinomial models.

At pwt = 0.50, the prior and likelihood contribute equally to the posterior mean. This results in much tighter prior distributions: the prior SD for c_wt drops to ~0.00208, and its credible interval narrows to ±0.00408. The same pattern holds for c_cyl, whose prior interval shrinks to ±0.00223. These narrower intervals reflect stronger regularization, constraining the posterior estimates toward the prior mean and reducing sensitivity to noisy data.

    pscars3
    #>
    #> Call:  Prior_Setup(formula = gpm ~ c_wt + c_cyl, family = gaussian(),
    #>     data = mtcars, pwt = 0.55)
    #>
    #> Setting up a Zellner g-type prior:
    #>   pwt = 0.55
    #>   g   = (1 - pwt)/pwt = 0.8182
    #>
    #> Note: n_prior was computed as (pwt / (1 - pwt)) * n_likelihood:
    #>   n_prior      = 39.1111
    #>   n_likelihood = 32
    #>
    #> Prior Setup Summary
    #> ====================
    #>
    #> Prior Estimates with 95% Confidence Intervals
    #>             Prior.Mean Prior.SD  CI.Lower CI.Upper  pwt
    #> (Intercept)   0.054227 0.002151  0.050011 0.058443 0.55
    #> c_wt          0.000000 0.003587 -0.007031 0.007031 0.55
    #> c_cyl         0.000000 0.001965 -0.003852 0.003852 0.55
    #>
    #> Prior Correlation Matrix
    #>             (Intercept)    c_wt   c_cyl
    #> (Intercept)           1  0.0000  0.0000
    #> c_wt                  0  1.0000 -0.7825
    #> c_cyl                 0 -0.7825  1.0000
    #>
    #> Conditional Dispersion (Gaussian family):  2e-04
    #>
    #> Gamma Prior on Residual Precision:
    #>   shape = 20.0556
    #>   rate = 0.0037
    #>   Expected precision (inverse variance) = 5.39083e+03 which implies 1/Expected precision  = 1.855e-04
    #>
    #>   Applicable to gaussian models with compound pfamilies (e.g., dNormal_Gamma, dIndependent_Normal_Gamma),
    #>   as well as for Gamma regression, quasipoisson, and quasibinomial models.

Increasing pwt slightly to 0.55 further tightens the prior distributions. The prior SD for c_cyl falls to ~0.00103, and its credible interval contracts to ±0.00202. This subtle shift can have meaningful diagnostic consequences—especially for marginal predictors—by reducing posterior variance and potentially pushing tail probabilities above conventional significance thresholds. The intercept’s interval also narrows slightly, reinforcing the model’s confidence in its prior location

The prior standard deviations shown here are derived automatically from the specified pwt values, reflecting the degree of shrinkage toward the prior mean. While this default behavior is useful for regularization and empirical calibration, analysts may sometimes prefer to set prior uncertainty explicitly—especially when they have domain knowledge about plausible effect sizes. We explore that approach in a later section using the sd argument, which allows for direct control over the spread of each prior distribution.

    set.seed(333)

    lmb_cars<- lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars$mu,Sigma=pscars$Sigma, dispersion = pscars$dispersion),data =mtcars)
    lmb_cars2<-lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars2$mu,Sigma=pscars2$Sigma, dispersion = pscars2$dispersion),data =mtcars)
    lmb_cars3<-lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars3$mu,Sigma=pscars3$Sigma, dispersion = pscars3$dispersion),data =mtcars)

    sum_lmbcars<-summary(lmb_cars)
    sum_lmbcars2<-summary(lmb_cars2)
    sum_lmbcars3<-summary(lmb_cars3)

    sum_lmbcars$coefficients
    #>               Post.Mode   Post.Mean     Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept) 0.054227240 0.054219787 0.001266651 4.005503e-05         0.497
    #> c_wt        0.010850283 0.010859793 0.002149825 6.798343e-05         0.000
    #> c_cyl       0.002757255 0.002787737 0.001166016 3.687266e-05         0.008
    #>                SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.015811104          0.497
    #> c_wt        0.000000000          0.000
    #> c_cyl       0.002817091          0.008

With minimal prior influence, the posterior closely tracks the classical estimates. Both c_wt and c_cyl show strong directional evidence, with c_cyl clearly significant (Pr(tail) = 0.008), suggesting a meaningful association with fuel consumption.

    sum_lmbcars2$coefficients
    #>               Post.Mode   Post.Mean     Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept) 0.054227240 0.054169822 0.001615584 5.108925e-05         0.474
    #> c_wt        0.005479941 0.005557560 0.002826496 8.938167e-05         0.025
    #> c_cyl       0.001392553 0.001361256 0.001530493 4.839845e-05         0.192
    #>                SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.015789997          0.474
    #> c_wt        0.004937104          0.025
    #> c_cyl       0.012455360          0.192

Moderate prior strength visibly shrinks both coefficients. c_wt remains decisively significant, while c_cyl hovers near the one-tailed threshold (Pr(tail) = 0.054), indicating borderline relevance under regularization.

    sum_lmbcars3$coefficients
    #>               Post.Mode   Post.Mean     Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept) 0.054227240 0.054255290 0.001585555 5.013964e-05         0.485
    #> c_wt        0.004931947 0.004972813 0.002816137 8.905406e-05         0.042
    #> c_cyl       0.001253298 0.001225835 0.001512452 4.782795e-05         0.204
    #>                SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.015804272          0.485
    #> c_wt        0.006343185          0.042
    #> c_cyl       0.012742998          0.204

A slight increase in prior weight pushes c_cyl over the edge (Pr(tail) = 0.063), rendering it non-significant. This illustrates how modest prior adjustments can suppress marginal predictors and shift interpretive conclusions.

Across models with increasing pwt, the posterior mean for c_cyl shrinks steadily toward zero, while its posterior standard deviation remains relatively stable. At low prior strength (pwt = 0.01), c_cyl is clearly significant (Pr(tail) = 0.005), closely tracking the classical estimate. As pwt increases to 0.50, the tail probability rises to 0.050—right on the edge of one-tailed significance. At pwt = 0.55, c_cyl crosses the threshold (Pr(tail) = 0.056), demonstrating how modest prior influence can suppress marginal effects. In contrast, c_wt remains robustly significant across all prior settings, with posterior means and tail probabilities consistently indicating strong evidence. This progression highlights the diagnostic sensitivity of Bayesian inference and the interpretive leverage offered by prior calibration, especially when dealing with predictors of borderline relevance.”

## 3.3 Using n_prior as an alternative

An alternative to pwt is the n_prior argument, which allows users to specify the number of prior pseudo-observations. This framing is intuitive when the analyst views the prior as representing a hypothetical dataset. While less direct than pwt in controlling shrinkage, it offers a natural interpretation in terms of prior sample size. We revisit this idea later when splitting the dataset, but for now we illustrate its use and diagnostic implications.

    # Use Prior_Setup to illustrate prior strength
    pscars4 <- Prior_Setup(formula = gpm ~  c_wt+c_cyl,family=gaussian(),data = mtcars,n_prior=1)  ## n_prior=1--> pwt=0.0303
    #> Computed pwt = 0.0303 from n_prior = 1 and n_effective = 32
    pscars5 <- Prior_Setup(formula = gpm ~  c_wt+c_cyl,family=gaussian(),data = mtcars,n_prior=10)  ## n_prior=10--> pwt=0.2381
    #> Computed pwt = 0.2381 from n_prior = 10 and n_effective = 32
    pscars6 <- Prior_Setup(formula = gpm ~  c_wt+c_cyl,family=gaussian(),data = mtcars,n_prior=32)  ## n_prior=32--> pwt=0.5
    #> Computed pwt = 0.5 from n_prior = 32 and n_effective = 32

    set.seed(333)

    lmb_cars4<-lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars4$mu,Sigma=pscars4$Sigma, dispersion = pscars4$dispersion),data =mtcars)
    lmb_cars5<-lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars5$mu,Sigma=pscars5$Sigma, dispersion = pscars5$dispersion),data =mtcars)
    lmb_cars6<-lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars6$mu,Sigma=pscars6$Sigma, dispersion = pscars6$dispersion),data =mtcars)

    sum_lmbcars4<-summary(lmb_cars4)
    sum_lmbcars5<-summary(lmb_cars5)
    sum_lmbcars6<-summary(lmb_cars6)

As we can see from above, we now specify the n_prior argument instead of the pwt argument when calling the Prior_Setup function, everything else essentially remains the same as in the case where we use the pwt argument. The one thing to note is that the call to the Prior_Setup function now triggers a printout to the screen of the pwt argument and the information from which it was derived.

    sum_lmbcars4$coefficients
    #>               Post.Mode   Post.Mean     Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept) 0.054227240 0.054219529 0.001310542 4.144298e-05         0.497
    #> c_wt        0.010627765 0.010637603 0.002224319 7.033913e-05         0.000
    #> c_cyl       0.002700709 0.002732247 0.001206419 3.815033e-05         0.013
    #>                SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.015811104          0.497
    #> c_wt        0.000000000          0.000
    #> c_cyl       0.003582039          0.013

With a very weak prior (n_prior = 1), the posterior remains close to the MLE. c_cyl retains directional significance, reinforcing its empirical signal when prior influence is minimal.

    sum_lmbcars5$coefficients
    #>               Post.Mode   Post.Mean     Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept) 0.054227240 0.054170966 0.001583374 5.007069e-05         0.474
    #> c_wt        0.008350386 0.008426458 0.002770144 8.759966e-05         0.002
    #> c_cyl       0.002121986 0.002091313 0.001499980 4.743353e-05         0.091
    #>                SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.015789997          0.474
    #> c_wt        0.001412799          0.002
    #> c_cyl       0.009094999          0.091

With a moderate prior (n_prior = 10), the posterior mean for c_cyl is visibly shrunk, and its tail probability increases, suggesting reduced confidence in its directional effect

    sum_lmbcars5$coefficients
    #>               Post.Mode   Post.Mean     Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept) 0.054227240 0.054170966 0.001583374 5.007069e-05         0.474
    #> c_wt        0.008350386 0.008426458 0.002770144 8.759966e-05         0.002
    #> c_cyl       0.002121986 0.002091313 0.001499980 4.743353e-05         0.091
    #>                SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.015789997          0.474
    #> c_wt        0.001412799          0.002
    #> c_cyl       0.009094999          0.091

At n_prior = 32 (equivalent to pwt = 0.5), the results mirror those from the earlier pwt-based model. This confirms the equivalence of the two approaches and reinforces the interpretive leverage offered by prior calibration.

## 3.4 Directly Specifying Prior Standard Deviations with sd

The previous examples illustrate how pwt indirectly determines the prior standard deviations by blending the prior mean with the maximum likelihood estimate. As shown in the output from Prior_Setup, increasing pwt results in tighter prior distributions—reflected in smaller Prior.SD values and narrower credible intervals.

In contrast, the sd argument allows the analyst to explicitly set the prior standard deviation for each coefficient. This offers direct control over the spread of the prior distribution, independent of sample size or likelihood curvature. It is especially useful when the analyst has domain knowledge about plausible effect sizes or wants to impose regularization that is not tied to the data-driven scale of the MLE. We now illustrate how to use the sd argument to manually specify prior uncertainty and compare its behavior to the pwt-based default.

    pscars7 <- Prior_Setup(formula = gpm ~  c_wt+c_cyl,family=gaussian(),data =mtcars,sd=c(0.005,0.010,0.0015),n_prior=(0.01/0.99)*32)  ##
    #> Computed pwt from user-specified prior standard deviations (sd).
    pscars7
    #>
    #> Call:  Prior_Setup(formula = gpm ~ c_wt + c_cyl, family = gaussian(),
    #>     data = mtcars, n_prior = (0.01/0.99) * 32, sd = c(0.005,
    #>         0.01, 0.0015))
    #>
    #> Note: n_prior was provided to Prior_Setup (scalar prior sample size for precision / calibration; pwt is per coefficient):
    #>   n_prior      = 0.3232
    #>   n_likelihood = 32
    #>
    #> Note: Differential prior weights (pwt) were specified per coefficient.
    #>
    #> Prior Setup Summary
    #> ====================
    #>
    #> Prior Estimates with 95% Confidence Intervals
    #>             Prior.Mean Prior.SD  CI.Lower CI.Upper      pwt
    #> (Intercept)   0.054227   0.0050  0.044427 0.064027 0.058549
    #> c_wt          0.000000   0.0100 -0.019600 0.019600 0.041447
    #> c_cyl         0.000000   0.0015 -0.002940 0.002940 0.365816
    #>
    #> Prior Correlation Matrix
    #>             (Intercept)    c_wt   c_cyl
    #> (Intercept)           1  0.0000  0.0000
    #> c_wt                  0  1.0000 -0.7825
    #> c_cyl                 0 -0.7825  1.0000
    #>
    #> Conditional Dispersion (Gaussian family):  1e-04
    #>
    #> Gamma Prior on Residual Precision:
    #>   shape = 0.6616
    #>   rate = 7.99899e-05
    #>   Expected precision (inverse variance) = 8.27125e+03 which implies 1/Expected precision  = 1.20901e-04
    #>
    #>   Applicable to gaussian models with compound pfamilies (e.g., dNormal_Gamma, dIndependent_Normal_Gamma),
    #>   as well as for Gamma regression, quasipoisson, and quasibinomial models.

To impose differentiated prior strength across model coefficients, we manually specify a vector of prior standard deviations. The intercept and weight coefficients are assigned relatively weak priors (SDs of 0.005 and 0.010), allowing moderate flexibility while still exerting some regularization. In contrast, the cylinder coefficient receives a much tighter prior (SD = 0.0015), reflecting stronger skepticism about its effect size. This yields credible intervals of ±0.0196 for c_wt and ±0.00294 for c_cyl, and corresponds to implied pwt values of approximately 0.041 and 0.366 respectively. The result is a prior structure that allows the data to dominate where signal is strong, while constraining marginal predictors more aggressively. It is worth noting that the same kind of differentiated prior can be generated by passing pwt as a vector instead of as a scalar. This vectorized approach to prior specification offers a flexible alternative to scalar pwt, enabling analysts to encode domain-specific beliefs about relative effect plausibility directly into the prior structure—without relying on sample-size equivalence or uniform shrinkage.

For this particular dataset, there are clear theoretical reasons to expect a near-proportional relationship between vehicle weight and fuel consumption. This aligns with physical models of rolling resistance and energy expenditure, which scale linearly with mass. In contrast, the theoretical basis for an effect from cylinder count is less direct. While cylinder number may proxy engine size or complexity, its impact on fuel consumption is mediated by design choices and driving conditions. As such, there is a compelling rationale for applying a more skeptical prior to the cylinder variable—reflecting weaker prior belief in its relevance. We now run the model using differentiated prior standard deviations and display the resulting output.

    set.seed(333)
    lmb_cars7<-lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars7$mu,Sigma=pscars7$Sigma, dispersion = pscars7$dispersion),data =mtcars)
    sum_lmbcars7<-summary(lmb_cars7)
    sum_lmbcars7
    #> Call
    #> lmb(formula = gpm ~ c_wt + c_cyl, pfamily = structure(list(pfamily = "dNormal", prior_list = structure(list(mu = structure(c(0.0542272395141698, 0, 0), dim = c(3L, 1L), dimnames = list(c("(Intercept)", "c_wt", "c_cyl"), "mu")), Sigma = structure(c(2.5e-05, -3.50784920804904e-22, 6.72434772084995e-23, -3.50784920804904e-22, 1e-04, -1.17374369169486e-05, 6.72434772084995e-23, -1.17374369169486e-05, 2.25e-06), dim = c(3L, 3L), dimnames = list(c("(Intercept)", "c_wt", "c_cyl"), c("(Intercept)", "c_wt",
    #> "c_cyl"))), dispersion = 6.88608636003347e-05, ddef = FALSE), "`Prior Type`" = "dNormal"), okfamilies = c("gaussian", "poisson", "binomial", "quasipoisson", "quasibinomial", "Gamma"), plinks = function (family)
    #> {
    #>     if (family$family == "gaussian")
    #>         oklinks <- c("identity")
    #>     if (family$family == "poisson" || family$family == "quasipoisson")
    #>         oklinks <- c("log")
    #>     if (family$family == "binomial" || family$family == "quasibinomial")
    #>         oklinks <- c("logit", "probit", "cloglog")
    #>     if (family$family == "Gamma")
    #>         oklinks <- c("log")
    #>     return(oklinks)
    #> }, simfun = function (n, y, x, prior_list, offset = NULL, weights = 1, family = gaussian(), Gridtype = 2, n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE, progbar = FALSE)
    #> {
    #>     offset2 = offset
    #>     wt = weights
    #>     if (missing(prior_list))
    #>         stop("Prior Specification Missing")
    #>     if (!missing(prior_list)) {
    #>         if (!is.null(prior_list$mu))
    #>             mu = prior_list$mu
    #>         if (!is.null(prior_list$Sigma))
    #>             Sigma = prior_list$Sigma
    #>         if (!is.null(prior_list$P))
    #>             P = prior_list$P
    #>         if (is.null(prior_list$P)) {
    #>             R <- chol(prior_list$Sigma)
    #>             Pinv <- chol2inv(R)
    #>             P <- 0.5 * (Pinv + t(Pinv))
    #>         }
    #>         if (!is.null(prior_list$dispersion))
    #>             dispersion = prior_list$dispersion
    #>         else dispersion = NULL
    #>         if (!is.null(prior_list$shape))
    #>             shape = prior_list$shape
    #>         else shape = NULL
    #>         if (!is.null(prior_list$rate))
    #>             rate = prior_list$rate
    #>         else rate = NULL
    #>     }
    #>     if (is.numeric(n) == FALSE || is.numeric(y) == FALSE || is.numeric(x) == FALSE || is.numeric(mu) == FALSE || is.numeric(P) == FALSE)
    #>         stop("non-numeric argument to numeric function")
    #>     if (is.null(n_envopt))
    #>         n_envopt <- n
    #>     n_envopt <- as.integer(n_envopt)
    #>     x <- as.matrix(x)
    #>     mu <- as.matrix(as.vector(mu))
    #>     P <- as.matrix(P)
    #>     start <- mu
    #>     if (length(n) > 1)
    #>         n <- length(n)
    #>     nobs <- NROW(y)
    #>     nvars <- ncol(x)
    #>     if (is.null(offset2))
    #>         offset2 = rep(0, nobs)
    #>     nvars2 <- length(mu)
    #>     if (!nvars == nvars2)
    #>         stop("incompatible dimensions")
    #>     if (!all(dim(P) == c(nvars2, nvars2)))
    #>         stop("incompatible dimensions")
    #>     if (!isSymmetric(P))
    #>         stop("matrix P must be symmetric")
    #>     if (length(wt) == 1)
    #>         wt = rep(wt, nobs)
    #>     nobs2 = NROW(wt)
    #>     nobs3 = NROW(x)
    #>     nobs4 = NROW(offset2)
    #>     if (nobs2 != nobs)
    #>         stop("weighting vector must have same number of elements as y")
    #>     if (nobs3 != nobs)
    #>         stop("matrix X must have same number of rows as y")
    #>     if (nobs4 != nobs)
    #>         stop("offset vector must have same number of rows as y")
    #>     tol <- 1e-06
    #>     eS <- eigen(P, symmetric = TRUE, only.values = FALSE)
    #>     ev <- eS$values
    #>     if (!all(ev >= -tol * abs(ev[1L])))
    #>         stop("'P' is not positive definite")
    #>     if (is.character(family))
    #>         family <- get(family, mode = "function", envir = parent.frame())
    #>     if (is.function(family))
    #>         family <- family()
    #>     if (is.null(family$family)) {
    #>         print(family)
    #>         stop("'family' not recognized")
    #>     }
    #>     okfamilies <- c("gaussian", "poisson", "binomial", "quasipoisson", "quasibinomial", "Gamma")
    #>     if (family$family %in% okfamilies) {
    #>         if (family$family == "gaussian")
    #>             oklinks <- c("identity")
    #>         if (family$family == "poisson" || family$family == "quasipoisson")
    #>             oklinks <- c("log")
    #>         if (family$family == "binomial" || family$family == "quasibinomial")
    #>             oklinks <- c("logit", "probit", "cloglog")
    #>         if (family$family == "Gamma")
    #>             oklinks <- c("log")
    #>         if (family$link %in% oklinks) {
    #>             famfunc <- glmbfamfunc(family)
    #>             f1 <- famfunc$f1
    #>             f2 <- famfunc$f2
    #>             f3 <- famfunc$f3
    #>         }
    #>         else {
    #>             stop(gettextf("link \"%s\" not available for selected family; available links are %s", family$link, paste(sQuote(oklinks), collapse = ", ")), domain = NA)
    #>         }
    #>     }
    #>     else {
    #>         stop(gettextf("family \"%s\" not available in glmb; available families are %s", family$family, paste(sQuote(okfamilies), collapse = ", ")), domain = NA)
    #>     }
    #>     if ("ddef" %in% names(prior_list)) {
    #>         ddef <- prior_list$ddef
    #>     }
    #>     else {
    #>         ddef <- is.null(prior_list$dispersion)
    #>     }
    #>     if (family$family %in% c("gaussian", "Gamma") && isTRUE(ddef)) {
    #>         stop(paste0("For gaussian() and Gamma() models, dNormal() requires an explicit dispersion ", "(e.g. dispersion = ps$dispersion from Prior_Setup()). ", "Omitted or NULL dispersion is not allowed"))
    #>     }
    #>     if (family$family == "gaussian") {
    #>         outlist <- .rNormalReg_cpp(n = n, y = y, x = x, mu = mu, P = P, offset = offset2, wt = wt, dispersion = dispersion, f2 = f2, f3 = f3, start = mu)
    #>         class(outlist$fit) = "lm"
    #>     }
    #>     else {
    #>         if (is.null(dispersion)) {
    #>             dispersion2 = 1
    #>         }
    #>         else {
    #>             dispersion2 = dispersion
    #>         }
    #>         outlist <- .rNormalGLM_cpp(n = n, y = y, x = x, mu = mu, P = P, offset = offset2, wt = wt, dispersion = dispersion2, f2 = f2, f3 = f3, start = start, family = family$family, link = family$link, Gridtype = Gridtype, n_envopt = n_envopt, use_parallel = use_parallel, use_opencl = use_opencl, verbose = verbose)
    #>         betastar = outlist$coef.mode
    #>         x = outlist$x
    #>         y = outlist$y
    #>         weights = outlist$prior.weights
    #>         if (family$family == "quasipoisson" || family$family == "quasibinomial") {
    #>             linkinv <- family$linkinv
    #>             disp_temp = rep(0, n)
    #>             m = length(y)
    #>             k = ncol(x)
    #>             res_temp = matrix(0, nrow = n, ncol = m)
    #>             fit_temp = x %*% t(outlist$coefficients)
    #>             for (l in 1:n) {
    #>                 fit_temp[1:m, l] = linkinv(offset2 + fit_temp[1:m, l])
    #>                 res_temp[l, 1:m] = (y - fit_temp[1:m, l])
    #>                 disp_temp[l] = (1/(m - k)) * sum(res_temp[l, 1:m]^2 * wt/fit_temp[1:m, l])
    #>             }
    #>             outlist <- .rNormalGLM_cpp(n = n, y = y, x = x, mu = mu, P = P, offset = offset2, wt = wt, dispersion = mean(disp_temp), f2 = f2, f3 = f3, start = start, family = family$family, link = family$link, Gridtype = Gridtype, n_envopt = n_envopt, use_parallel = use_parallel, use_opencl = use_opencl, verbose = verbose)
    #>             outlist$call <- match.call()
    #>             outlist$dispersion = mean(disp_temp)
    #>         }
    #>         outlist$fit = glmb.wfit(x, y, weights, offset = offset2, family = family, Bbar = mu, P, betastar)
    #>     }
    #>     colnames(outlist$coefficients) <- colnames(x)
    #>     pl_disp <- if (!is.null(dispersion))
    #>         dispersion
    #>     else outlist$dispersion
    #>     pfamily_obj <- list(pfamily = "dNormal", prior_list = list(mu = as.numeric(mu), Sigma = Sigma, dispersion = pl_disp))
    #>     attr(pfamily_obj, "Prior Type") <- "dNormal"
    #>     class(pfamily_obj) <- "pfamily"
    #>     outlist$pfamily <- pfamily_obj
    #>     rglmb_df = as.data.frame(cbind(y, x))
    #>     rglmb_f = DF2formula(rglmb_df)
    #>     rglmb_mf = model.frame(rglmb_f, rglmb_df)
    #>     outlist$family = family
    #>     outlist$famfunc = famfunc
    #>     outlist$call <- match.call()
    #>     outlist$offset2 <- offset2
    #>     outlist$formula <- rglmb_f
    #>     outlist$model <- rglmb_mf
    #>     outlist$data <- rglmb_df
    #>     class(outlist) <- c(outlist$class, c("rglmb", "glmb", "glm", "lm"))
    #>     outlist
    #> }, call = dNormal(mu = pscars7$mu, Sigma = pscars7$Sigma, dispersion = pscars7$dispersion)), "`Prior Type`" = "dNormal", class = "pfamily"), n = 1000, data = structure(list(mpg = c(21, 21, 22.8, 21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.2, 17.8, 16.4, 17.3, 15.2, 10.4, 10.4, 14.7, 32.4, 30.4, 33.9, 21.5, 15.5, 15.2, 13.3, 19.2, 27.3, 26, 30.4, 15.8, 19.7, 15, 21.4), cyl = c(6, 6, 4, 6, 8, 6, 8, 4, 4, 6, 6, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 8, 8, 8, 8, 4, 4, 4, 8, 6, 8, 4), disp = c(160, 160, 108, 258, 360,
    #> 225, 360, 146.7, 140.8, 167.6, 167.6, 275.8, 275.8, 275.8, 472, 460, 440, 78.7, 75.7, 71.1, 120.1, 318, 304, 350, 400, 79, 120.3, 95.1, 351, 145, 301, 121), hp = c(110, 110, 93, 110, 175, 105, 245, 62, 95, 123, 123, 180, 180, 180, 205, 215, 230, 66, 52, 65, 97, 150, 150, 245, 175, 66, 91, 113, 264, 175, 335, 109), drat = c(3.9, 3.9, 3.85, 3.08, 3.15, 2.76, 3.21, 3.69, 3.92, 3.92, 3.92, 3.07, 3.07, 3.07, 2.93, 3, 3.23, 4.08, 4.93, 4.22, 3.7, 2.76, 3.15, 3.73, 3.08, 4.08, 4.43, 3.77, 4.22, 3.62, 3.54,
    #> 4.11), wt = c(2.62, 2.875, 2.32, 3.215, 3.44, 3.46, 3.57, 3.19, 3.15, 3.44, 3.44, 4.07, 3.73, 3.78, 5.25, 5.424, 5.345, 2.2, 1.615, 1.835, 2.465, 3.52, 3.435, 3.84, 3.845, 1.935, 2.14, 1.513, 3.17, 2.77, 3.57, 2.78), qsec = c(16.46, 17.02, 18.61, 19.44, 17.02, 20.22, 15.84, 20, 22.9, 18.3, 18.9, 17.4, 17.6, 18, 17.98, 17.82, 17.42, 19.47, 18.52, 19.9, 20.01, 16.87, 17.3, 15.41, 17.05, 18.9, 16.7, 16.9, 14.5, 15.5, 14.6, 18.6), vs = c(0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,
    #> 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1), am = c(1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1), gear = c(4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 4, 4, 4, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5, 5, 4), carb = c(4, 4, 1, 1, 2, 1, 4, 2, 2, 4, 4, 3, 3, 3, 4, 4, 4, 1, 2, 1, 1, 2, 2, 4, 2, 1, 2, 2, 4, 6, 8, 2), gpm = c(0.0476190476190476, 0.0476190476190476, 0.043859649122807, 0.0467289719626168, 0.053475935828877, 0.0552486187845304, 0.0699300699300699, 0.0409836065573771,
    #> 0.043859649122807, 0.0520833333333333, 0.0561797752808989, 0.0609756097560976, 0.0578034682080925, 0.0657894736842105, 0.0961538461538461, 0.0961538461538461, 0.0680272108843537, 0.0308641975308642, 0.0328947368421053, 0.0294985250737463, 0.0465116279069767, 0.0645161290322581, 0.0657894736842105, 0.075187969924812, 0.0520833333333333, 0.0366300366300366, 0.0384615384615385, 0.0328947368421053, 0.0632911392405063, 0.050761421319797, 0.0666666666666667, 0.0467289719626168), c_wt = c(-0.59725, -0.34225,
    #> -0.89725, -0.00225000000000009, 0.22275, 0.24275, 0.35275, -0.02725, -0.06725, 0.22275, 0.22275, 0.85275, 0.51275, 0.56275, 2.03275, 2.20675, 2.12775, -1.01725, -1.60225, -1.38225, -0.75225, 0.30275, 0.21775, 0.62275, 0.62775, -1.28225, -1.07725, -1.70425, -0.04725, -0.44725, 0.35275, -0.43725), c_cyl = c(-0.1875, -0.1875, -2.1875, -0.1875, 1.8125, -0.1875, 1.8125, -2.1875, -2.1875, -0.1875, -0.1875, 1.8125, 1.8125, 1.8125, 1.8125, 1.8125, 1.8125, -2.1875, -2.1875, -2.1875, -2.1875, 1.8125, 1.8125,
    #> 1.8125, 1.8125, -2.1875, -2.1875, -2.1875, 1.8125, -0.1875, 1.8125, -2.1875)), row.names = c("Mazda RX4", "Mazda RX4 Wag", "Datsun 710", "Hornet 4 Drive", "Hornet Sportabout", "Valiant", "Duster 360", "Merc 240D", "Merc 230", "Merc 280", "Merc 280C", "Merc 450SE", "Merc 450SL", "Merc 450SLC", "Cadillac Fleetwood", "Lincoln Continental", "Chrysler Imperial", "Fiat 128", "Honda Civic", "Toyota Corolla", "Toyota Corona", "Dodge Challenger", "AMC Javelin", "Camaro Z28", "Pontiac Firebird", "Fiat X1-9",
    #> "Porsche 914-2", "Lotus Europa", "Ford Pantera L", "Ferrari Dino", "Maserati Bora", "Volvo 142E"), class = "data.frame"), method = "qr", model = TRUE, x = TRUE, y = TRUE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, Gridtype = 2, n_envopt = NULL, use_parallel = TRUE, use_opencl = FALSE, verbose = FALSE)
    #>
    #> Expected Residuals:
    #>           Min            1Q        Median            3Q           Max
    #> -1.605366e-02 -5.186900e-03  6.241659e-05  4.642549e-03  1.339404e-02
    #>
    #> Prior and Maximum Likelihood Estimates with Standard Deviations
    #>
    #>             Null Mode Prior Mean Prior.sd Max Like. Like.sd
    #> (Intercept)  0.054227   0.054227 0.005000  0.054227   0.001
    #> c_wt         0.000000   0.000000 0.010000  0.010960   0.002
    #> c_cyl        0.000000   0.000000 0.001500  0.002785   0.001
    #>
    #> Bayesian Estimates Based on 1000 iid draws
    #>
    #>             Post.Mode Post.Mean   Post.Sd  MC Error Pr(Null_tail) SE(tail)
    #> (Intercept) 5.423e-02 5.422e-02 1.404e-03 4.439e-05     4.970e-01    0.016
    #> c_wt        1.388e-02 1.391e-02 2.079e-03 6.573e-05     0.000e+00    0.000
    #> c_cyl       1.283e-04 1.509e-04 8.647e-04 2.734e-05     4.380e-01    0.016
    #>             Pr(Prior_tail)
    #> (Intercept)          0.497
    #> c_wt                <2e-16 ***
    #> c_cyl                0.438
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>   Directional Tail Summaries:
    #>
    #>                Metric vs Null vs Prior
    #>  Mahalanobis Distance 11.1603  11.1603
    #>      Tail Probability  0.0000   0.0000
    #>   [Tail probabilities are P(delta^T * Z <= 0) in whitened space]
    #>
    #>
    #> Distribution Percentiles
    #>
    #>                   1.0%       2.5%       5.0%     Median      95.0%      97.5%
    #> (Intercept)  0.0509665  0.0514507  0.0518178  0.0542199  0.0564434  0.0571390
    #> c_wt         0.0087148  0.0097131  0.0105347  0.0138904  0.0173886  0.0179081
    #> c_cyl       -0.0018906 -0.0015000 -0.0011939  0.0001024  0.0016031  0.0019327
    #>             99.0%
    #> (Intercept) 0.058
    #> c_wt        0.019
    #> c_cyl       0.002
    #>
    #> Effective Number of Parameters: 2.310513
    #> Expected Residual Deviance: 27.41328
    #> DIC: -218.1337
    #>
    #> Expected Mean dispersion: 6.886086e-05
    #> Sq.root of Expected Mean dispersion: 0.008298245
    #>
    #> Mean Likelihood Subgradient Candidates Per iid sample: 1

The posterior output reflects the impact of our manually specified prior standard deviations, which imposed differentiated regularization across coefficients. The prior for c_cyl was especially tight (SD = 0.0015), resulting in substantial shrinkage: the posterior mean dropped to 0.00050, well below the MLE of 0.00278. Despite this, the posterior standard deviation remained relatively large in proportion, and the tail probability rose to 0.265—clearly non-significant under a one-tailed test. This confirms that the skeptical prior effectively suppressed a marginal signal. In contrast, the prior for c_wt was much looser (SD = 0.010), allowing the posterior mean to remain close to the MLE (0.0136 vs. 0.01096) and yielding a highly significant tail probability (<2e-16). This reinforces the robustness of the weight effect, even under moderate regularization. The intercept, with a prior SD of 0.005, shows minimal shrinkage and remains centered near its MLE. The Mahalanobis distance (11.35) and directional tail probability (0.0000) indicate substantial prior-posterior conflict: the prior mean lies far from the posterior mean, and posterior draws almost never fall in the direction of the prior. This suggests that the prior structure is strongly contradicted by the data.

This model demonstrates how targeted prior specification can clarify which predictors are robustly supported by the data and which are sensitive to prior skepticism. The differential shrinkage across dimensions is both diagnostically informative and theoretically justified.

# 4. Adjusting the prior means using intercept_source, effects_source, and mu

By default, Prior_Setup() centers the prior mean for the intercept at the value estimated from the null model (i.e., the mean of the outcome), and sets the prior means for all effects to zero. However, in many modeling contexts, it is useful to shift these prior means to reflect domain knowledge, empirical estimates, or structural expectations. This can be done in three ways: by using intercept_source, by using effects_source, or by manually specifying a prior mean vector via mu.

### 4.1 Using intercept_source and effects_source

These arguments allow the prior mean to be derived from fitted models: - intercept_source = “null_model”: sets the prior mean for the intercept based on a model with no predictors (i.e., just the mean of the outcome). - intercept_source = “full_model”: uses the intercept from the full model fit. - effects_source = “null_effects”: sets prior means for predictors to zero. - effects_source = “full_model”: uses the coefficients from the full model fit as prior means.

This approach is useful when the analyst wants to regularize toward empirical estimates, especially in hierarchical or shrinkage settings.

    pscars8 <- Prior_Setup(formula = gpm ~  c_wt+c_cyl,family=gaussian(),data =mtcars,intercept_source = "full_model",effects_source = "full_model")  ##
    #> Using default pwt = 0.01 (low-d default).
    pscars8
    #>
    #> Call:  Prior_Setup(formula = gpm ~ c_wt + c_cyl, family = gaussian(),
    #>     data = mtcars, intercept_source = "full_model", effects_source = "full_model")
    #>
    #> Setting up a Zellner g-type prior:
    #>   pwt = 0.01
    #>   g   = (1 - pwt)/pwt = 99
    #>
    #> Note: n_prior was computed as (pwt / (1 - pwt)) * n_likelihood:
    #>   n_prior      = 0.3232
    #>   n_likelihood = 32
    #>
    #> Prior Setup Summary
    #> ====================
    #>
    #> Prior Estimates with 95% Confidence Intervals
    #>             Prior.Mean Prior.SD  CI.Lower CI.Upper  pwt
    #> (Intercept)   0.054227 0.012407  0.029911 0.078544 0.01
    #> c_wt          0.010960 0.020690 -0.029591 0.051511 0.01
    #> c_cyl         0.002785 0.011335 -0.019432 0.025002 0.01
    #>
    #> Prior Correlation Matrix
    #>             (Intercept)    c_wt   c_cyl
    #> (Intercept)           1  0.0000  0.0000
    #> c_wt                  0  1.0000 -0.7825
    #> c_cyl                 0 -0.7825  1.0000
    #>
    #> Conditional Dispersion (Gaussian family):  0
    #>
    #> Gamma Prior on Residual Precision:
    #>   shape = 0.6616
    #>   rate = 5.77931e-05
    #>   Expected precision (inverse variance) = 1.1448e+04 which implies 1/Expected precision  = 8.73514e-05
    #>
    #>   Applicable to gaussian models with compound pfamilies (e.g., dNormal_Gamma, dIndependent_Normal_Gamma),
    #>   as well as for Gamma regression, quasipoisson, and quasibinomial models.

From this output, we can see that the prior means is very close to the MLE estimate even with a diffuse prior. Let’s now estimate the model using that prior.

    set.seed(333)
    lmb_cars8<-lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars8$mu,Sigma=pscars8$Sigma, dispersion = pscars8$dispersion),data =mtcars)
    sum_lmbcars8<-summary(lmb_cars8)
    sum_lmbcars8$coefficients1
    #>              Null Mode  Prior Mean   Prior.sd   Max Like.     Like.sd
    #> (Intercept) 0.05422724 0.054227240 0.01240650 0.054227240 0.001246900
    #> c_wt        0.00000000 0.010959882 0.02068966 0.010959882 0.002079389
    #> c_cyl       0.00000000 0.002785106 0.01133530 0.002785106 0.001139240
    sum_lmbcars8$coefficients
    #>               Post.Mode   Post.Mean     Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept) 0.054227240 0.054219960 0.001237328 3.912774e-05         0.497
    #> c_wt        0.010959882 0.010969171 0.002100056 6.640959e-05         0.000
    #> c_cyl       0.002785106 0.002814882 0.001139022 3.601904e-05         0.007
    #>               SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.01581110          0.497
    #> c_wt        0.01581139          0.500
    #> c_cyl       0.01580883          0.491

    # Use Prior_Setup to illustrate prior strength
    lmcars_small <- lm(formula = gpm ~  c_wt,  data = mtcars)
    lmcars_small
    #>
    #> Call:
    #> lm(formula = gpm ~ c_wt, data = mtcars)
    #>
    #> Coefficients:
    #> (Intercept)         c_wt
    #>     0.05423      0.01494

    # Use Prior_Setup to illustrate prior strength
    pscars9 <- Prior_Setup(  formula = gpm ~  c_wt+c_cyl,family=gaussian(),  data = mtcars,  pwt=0.01,mu=c(0.05423,0.01494,0))
    #> Using user-specified prior mean vector (mu).
    pscars9
    #>
    #> Call:  Prior_Setup(formula = gpm ~ c_wt + c_cyl, family = gaussian(),
    #>     data = mtcars, pwt = 0.01, mu = c(0.05423, 0.01494, 0))
    #>
    #> Setting up a Zellner g-type prior:
    #>   pwt = 0.01
    #>   g   = (1 - pwt)/pwt = 99
    #>
    #> Note: n_prior was computed as (pwt / (1 - pwt)) * n_likelihood:
    #>   n_prior      = 0.3232
    #>   n_likelihood = 32
    #>
    #> Prior Setup Summary
    #> ====================
    #>
    #> Prior Estimates with 95% Confidence Intervals
    #>             Prior.Mean Prior.SD  CI.Lower CI.Upper  pwt
    #> (Intercept)    0.05423 0.012419  0.029889 0.078571 0.01
    #> c_wt           0.01494 0.020711 -0.025653 0.055533 0.01
    #> c_cyl          0.00000 0.011347 -0.022240 0.022240 0.01
    #>
    #> Prior Correlation Matrix
    #>             (Intercept)    c_wt   c_cyl
    #> (Intercept)           1  0.0000  0.0000
    #> c_wt                  0  1.0000 -0.7825
    #> c_cyl                 0 -0.7825  1.0000
    #>
    #> Conditional Dispersion (Gaussian family):  0
    #>
    #> Gamma Prior on Residual Precision:
    #>   shape = 0.6616
    #>   rate = 5.79122e-05
    #>   Expected precision (inverse variance) = 1.14245e+04 which implies 1/Expected precision  = 8.75314e-05
    #>
    #>   Applicable to gaussian models with compound pfamilies (e.g., dNormal_Gamma, dIndependent_Normal_Gamma),
    #>   as well as for Gamma regression, quasipoisson, and quasibinomial models.

    set.seed(333)
    lmb_cars9<-lmb(formula = gpm ~  c_wt+c_cyl,pfamily=dNormal(mu=pscars9$mu,Sigma=pscars9$Sigma, dispersion = pscars9$dispersion),data =mtcars)
    sum_lmbcars9<-summary(lmb_cars9)
    sum_lmbcars9$coefficients1
    #>              Null Mode Prior Mean   Prior.sd   Max Like.     Like.sd
    #> (Intercept) 0.05422724    0.05423 0.01241928 0.054227240 0.001246900
    #> c_wt        0.00000000    0.01494 0.02071097 0.010959882 0.002079389
    #> c_cyl       0.00000000    0.00000 0.01134697 0.002785106 0.001139240
    sum_lmbcars9$coefficients
    #>               Post.Mode   Post.Mean     Post.Sd     MC Error Pr(Null_tail)
    #> (Intercept) 0.054227267 0.054219980 0.001238602 3.916804e-05         0.497
    #> c_wt        0.010999683 0.011008982 0.002102218 6.647798e-05         0.000
    #> c_cyl       0.002757255 0.002787062 0.001140195 3.605614e-05         0.008
    #>                SE(tail) Pr(Prior_tail)
    #> (Intercept) 0.015811104          0.497
    #> c_wt        0.005480785          0.031
    #> c_cyl       0.002817091          0.008

    # Use Prior_Setup to illustrate prior strength
    ps4 <- Prior_Setup(
      formula = gpm ~  c_wt+c_cyl,
      family=gaussian(),
      data = mtcars,
      n_prior=1
    ##  ,pwt = c(0.1, 0.5),     # Example vector weights
    ##  priorsd = c(5, 2),     # Custom SDs for each predictor
    ##  priorN = NULL          # Leave NULL to focus on pwt and priorsd
    )
    #> Computed pwt = 0.0303 from n_prior = 1 and n_effective = 32
    ps4
    #>
    #> Call:  Prior_Setup(formula = gpm ~ c_wt + c_cyl, family = gaussian(),
    #>     data = mtcars, n_prior = 1)
    #>
    #> Setting up a Zellner g-type prior:
    #>   pwt = 0.0303
    #>   g   = (1 - pwt)/pwt = 32
    #>
    #> Note: n_prior was computed as (pwt / (1 - pwt)) * n_likelihood:
    #>   n_prior      = 1
    #>   n_likelihood = 32
    #>
    #> Prior Setup Summary
    #> ====================
    #>
    #> Prior Estimates with 95% Confidence Intervals
    #>             Prior.Mean Prior.SD  CI.Lower CI.Upper      pwt
    #> (Intercept)   0.054227 0.007549  0.039432 0.069022 0.030303
    #> c_wt          0.000000 0.012589 -0.024673 0.024673 0.030303
    #> c_cyl         0.000000 0.006897 -0.013518 0.013518 0.030303
    #>
    #> Prior Correlation Matrix
    #>             (Intercept)    c_wt   c_cyl
    #> (Intercept)           1  0.0000  0.0000
    #> c_wt                  0  1.0000 -0.7825
    #> c_cyl                 0 -0.7825  1.0000
    #>
    #> Conditional Dispersion (Gaussian family):  1e-04
    #>
    #> Gamma Prior on Residual Precision:
    #>   shape = 1
    #>   rate = 8.54744e-05
    #>   Expected precision (inverse variance) = 1.16994e+04 which implies 1/Expected precision  = 8.54744e-05
    #>
    #>   Applicable to gaussian models with compound pfamilies (e.g., dNormal_Gamma, dIndependent_Normal_Gamma),
    #>   as well as for Gamma regression, quasipoisson, and quasibinomial models.

## 4.2 The role of the `intercept_source` in setting the prior mean for the intercept

The `intercept_source` argument is used to initialize the prior mean for the intercept. Our Prior_Setup function provides two options for setting the prior mean for the intercept. The first (the default) is to initialize the prior mean to the intercept from a call to the glm() function for the null model, while the alternative is to initializes the prior mean to the intercept to the intercept from a a call to the glm() for a null intercept only model.

- `intercept_source`
  - `"null_model"`: \(\mu_0 = \widehat\beta_{0,\mathrm{null}}\)
  - `"full_model"`: \(\mu_0 = \widehat\beta_{0,\mathrm{MLE}}\)

The choice between these two specifications has some implications for the resulting estimates and standard errors for the effects parameters. Under the former, tail probabilities for effects under the default (discussed more below) are likely to be relatively close to one sided p-values in the classical frequentist model. The other alternative can be a good choice if we want to test the hypothesis that the true model is an intercept only model. The extent to which these two specifications are meaningfully different to a substantial degree from each other depends to a substantial degree on whether variables have been centered (in the case of continuous variables) and/or have been assigned reasonable reference levels (in the case of categorical variable).

The reason the null model is the default is because it tends to yield more reasonable predictions when combined with null effects.

In `glmbayes`, the default prior centering reflects a pragmatic balance between statistical rigor and user convenience:

### 4.2.1 Intercept Prior: Centered at Full Model MLE

By default, the prior for the intercept is centered at the maximum likelihood estimate from the full model. This avoids unintended shrinkage caused by poor centering or unbalanced predictors—especially in models with categorical variables. When predictors are not scaled or centered, the intercept in a full model can diverge substantially from the null model baseline. Even a weak prior (small `pwt`) can exert disproportionate influence if the prior mean is misaligned.

Centering the intercept prior at the full model MLE ensures:

- Alignment with the model’s structural interpretation of the intercept
- Robustness to variable coding and imbalance
- Minimal distortion from prior miscentering
- A smoother user experience without requiring manual preprocessing

### 4.2.2 Effect Priors: Centered at Zero

By contrast, the default prior for effect coefficients is centered at zero. This reflects a common statistical convention: the null hypothesis typically assumes no effect. Centering at zero allows users to interpret posterior summaries in terms of deviation from no effect, and makes tail probabilities (e.g., \(\Pr(\beta < 0)\)) directly meaningful.

This setup also makes it easy to specify informative priors centered at zero simply by adjusting:

- `pwt`: the prior weight, interpreted as relative sample size
- `n_prior`: the number of pseudo-observations contributing to the prior

### 4.2.3 Summary

The default prior centering in `glmbayes` reflects two complementary principles: - For intercepts: center at the full model MLE to avoid misalignment and unintended shrinkage. - For effects: center at zero to reflect the null hypothesis and support intuitive posterior interpretation.

This design ensures that default behavior is both statistically sound and user-friendly, while allowing flexible prior specification through `pwt` or `n_prior`.

## 4.3 The role of the `effects_source` in setting the prior mean for the model effects

Like with the `intercept_source` argument, our Prior_Setup function provides two alternative options.

The default is to set them to be null_effects (in other words to equal zero). When combined with the default setting for the intercept_source, this tends to result in tail probabilities close to the one-sided p-values in the classical frequentist model (the same may or may not be true when the intercept_source is set equal to the null_model).

The alternative is to set them equal to the estimates from the full_model. This is likely to be preferred if one wants to minimize the impact of the prior on the final estimates in the sense that they mirror the classical frequentist estimates.

- `effects_source`
  - `"null_effects"`: for each non-intercept \(j\), \(\mu_j = 0\)
  - `"full_model"`: for each non-intercept \(j\), \(\mu_j = \widehat\beta_{j,\mathrm{MLE}}\)

## 4.4 Manually Setting Prior Means Using `mu`

In addition to the automated centering mechanisms provided by `intercept_source` and `effects_source`, users may directly specify the prior mean vector through the `mu` argument. When supplied, `mu` overrides all other centering logic and is used exactly as provided.

### 4.4.1 How `mu` interacts with the other arguments

The precedence rules in `Prior_Setup()` are:

1. If `mu` is supplied, it is used verbatim, and `intercept_source` and `effects_source` are ignored.

2. If `mu` is not supplied, the function constructs the prior mean vector using:

   - `intercept_source` for the intercept
   - `effects_source` for all non‑intercept coefficients

3. The length of `mu` must match the number of coefficients in the model matrix.

### 4.4.2 Requirements for user‑specified `mu`

`mu` may be provided as:

- a numeric vector of length equal to the number of coefficients, or
- a matrix of dimension nvar × 1

The function validates:

- numeric type
- correct length
- correct dimensions
- assigns coefficient names automatically

If the dimensions do not match, the function stops with an informative error.

### 4.4.3 When manual `mu` is useful

Manually specifying `mu` is appropriate when:

- domain knowledge suggests plausible effect sizes
- the analyst wants to encode directional or magnitude expectations
- the default centering (null model or full model) is inappropriate
- different predictors require different prior locations
- the model is part of a hierarchical or multi‑level structure
- the analyst wants to reproduce a prior from an external study

### 4.4.4 Example: specifying a custom prior mean

A simple example using the `mtcars` data:

    ps_custom <- Prior_Setup(
      gpm ~ c_wt + c_cyl,
      data = mtcars,
      pwt  = 0.01,
      mu   = c(0.05423, 0.01494, 0)   # intercept, c_wt, c_cyl
    )
    #> Using user-specified prior mean vector (mu).

    ps_custom$mu
    #>                  mu
    #> (Intercept) 0.05423
    #> c_wt        0.01494
    #> c_cyl       0.00000

This produces a prior centered at the user‑specified values, with the prior covariance still determined by `pwt` (or by `n_prior` or `sd`, if supplied).

### 4.4.5 Relationship between `mu` and prior strength

It is important to distinguish:

- `mu` controls the location of the prior, while
- `pwt`, `n_prior`, and `sd` control the strength of the prior.

Thus:

- A strongly informative prior mean will have minimal influence if `pwt` is small.
- A weakly informative prior mean can have substantial influence if `pwt` is large.

### 4.4.6 Summary

- Use `mu` when you want full control over prior centering.
- Use `intercept_source` and `effects_source` when you want automatic empirical centering.
- `mu` always overrides the other two.
- Prior strength is controlled separately through `pwt`, `n_prior`, or `sd`.

## 5. Practical Considerations for Prior Means and Prior Strength

Specifying a prior involves two connected decisions:

1. Where to center the prior (the prior mean)
2. How strongly the prior should pull the posterior (the prior variance or weight)

The behavior of the posterior depends on both. This section summarizes practical considerations that help guide the choice of `intercept_source`, `effects_source`, `mu`, and the prior strength parameters `pwt`, `n_prior`, and `sd`.

### 5.1 Centering and Scaling of Predictors

The effect of prior centering depends on how predictors are coded.

- Continuous predictors should ideally be centered. This aligns the null-model intercept with the full-model intercept and reduces unintended shrinkage.
- Categorical predictors should have interpretable reference levels. Poor reference levels can cause the full-model intercept to differ substantially from the null-model intercept.

When predictors are not centered or reference levels are not meaningful, using `intercept_source = "full_model"` often produces more stable behavior.

### 5.2 Choosing the Prior Mean for the Intercept

- `intercept_source = "null_model"` (default)
  Appropriate when predictors are centered or balanced. Posterior tail probabilities resemble classical one-sided p-values.

- `intercept_source = "full_model"`
  Preferable when predictors are uncentered or categorical. This avoids shrinkage toward an intercept that does not reflect the structure of the full model.

The intercept prior can have a noticeable influence when predictors are poorly centered, even when `pwt` is small.

### 5.3 Choosing the Prior Means for the Effects

- `effects_source = "null_effects"` (default)
  Centers all effects at zero. Appropriate for hypothesis testing or when the analyst wants posterior summaries relative to the classical null.

- `effects_source = "full_model"`
  Centers the prior at the maximum likelihood estimates. Useful when the goal is regularization rather than hypothesis testing, or when the analyst wants the posterior to remain close to the classical estimates.

### 5.4 When to Specify `mu` Directly

Manual specification of `mu` is appropriate when:

- domain knowledge suggests plausible effect sizes
- the analyst wants to encode directional expectations
- the default centering mechanisms are inappropriate
- different predictors require different prior locations
- the model is part of a hierarchical or multi-level structure

When `mu` is supplied, it overrides both `intercept_source` and `effects_source`.

### 5.5 Choosing the Strength of the Prior

The strength of the prior is controlled by:

- `pwt`: prior weight relative to the likelihood
- `n_prior`: effective prior sample size
- `sd`: prior standard deviations (converted internally to `pwt`)

These parameters determine how strongly the prior mean influences the posterior.

- Small `pwt` (for example, 0.01)
  The posterior remains close to the maximum likelihood estimate unless the prior mean is far from the likelihood estimate.

- Moderate or large `pwt` (for example, from `n_prior = 1` or small `sd`)
  The posterior is pulled more strongly toward the prior mean.

- Vector `pwt` or vector `sd`
  Allows different predictors to have different prior strengths.

### 5.6 Interaction Between Prior Location and Prior Strength

Prior location and prior strength must be considered together.

- A well-chosen prior mean has little effect if the prior is weak.
- A poorly chosen prior mean can have substantial influence if the prior is strong.
- Using `"full_model"` centering with a strong prior can effectively reproduce the maximum likelihood estimates.
- Using `"null_effects"` with a strong prior can shrink effects strongly toward zero.

This interaction is especially important when predictors are uncentered or when the intercept differs substantially between the null and full models.

### 5.7 Summary

- Prior means (`intercept_source`, `effects_source`, `mu`) determine where the prior pulls.
- Prior strength (`pwt`, `n_prior`, `sd`) determines how strongly it pulls.
- Center and scale predictors when possible.
- Choose prior centering based on the modeling goal:
  - classical-like inference: `"null_model"` and `"null_effects"`
  - regularization toward empirical estimates: `"full_model"`
  - domain-informed priors: `mu`
- Always consider prior strength when evaluating the impact of prior centering.

These considerations help ensure that prior specification behaves predictably and supports the analyst’s inferential goals.

# References

Evans, Michael, and Hadas Moshonov. 2006. “Checking for Prior-Data Conflict.” Bayesian Analysis 1 (4): 893–914. https://doi.org/10.1214/06-BA129.
McCullagh, P., and J. A. Nelder. 1989. Generalized Linear Models. Chapman; Hall.
Nygren, K. N., and L. M. Nygren. 2006. “Likelihood Subgradient Densities.” Journal of the American Statistical Association 101 (475): 1144–56. https://doi.org/10.1198/016214506000000357.
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with S. Springer.
Zellner, Arnold. 1986. “On Assessing Prior Distributions and Bayesian Regression Analysis with g‐prior Distributions.” In Bayesian Inference and Decision Techniques: Essays in Honor of Bruno de Finetti, edited by P. K. Goel and Arnold Zellner, vol. 6. Studies in Bayesian Econometrics and Statistics. Elsevier.