--- title: "Chapter 15: Estimating models with unknown dispersion parameters" output: rmarkdown::html_vignette bibliography: REFERENCES.bib reference-section-title: References vignette: > %\VignetteIndexEntry{Chapter 15: Estimating models with unknown dispersion parameters} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning= FALSE, message=FALSE ) ``` ```{r setup} library(glmbayes) ``` ## 1. Introduction Earlier chapters developed Bayesian linear and generalized linear models under the assumption that the dispersion parameter, denoted by \( \phi \) (or equivalently the residual variance \( \sigma^{2} \)), was known. For Gaussian models this corresponds to fixing the variance of the error term, and for quasi‑families it corresponds to fixing the scale parameter. This chapter extends the framework to unknown dispersion parameters---a standard topic in Bayesian linear and generalized linear models [@OHaganForster2004; @Gelman2013]---showing how glmbayes supports: (1) Priors on the dispersion only (2) Joint conjugate Normal–Gamma priors (3) Independent Normal–Gamma priors • iid sampling via truncated Gamma priors • two‑block Gibbs sampling Full derivations for Gaussian dispersion, shape, and rate calibration from `Prior_Setup()` / `compute_gaussian_prior()`—and how those objects connect to conjugate and independent Normal--Gamma priors—are given in **Chapter A12** (). All Gaussian examples use the **Dobson** plant‑weight data [@Dobson1990], first introduced in Chapter 03. ```{r dobson} ## Annette Dobson (1990) "An Introduction to Generalized Linear Models". ## Page 9: Plant Weight Data. ctl <- c(4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14) trt <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69) group <- gl(2, 10, 20, labels = c("Ctl", "Trt")) weight <- c(ctl, trt) ``` ## 2. Exponential-Family Background and Log-Concavity Many likelihoods used in generalized linear models belong to the exponential family. With observation weights \(w_1,\dots,w_n \ge 0\), the weighted exponential-family density is \[ f(y \mid \theta,\phi,w) = \exp\left\{ \sum_{i=1}^{n} w_i \left[ \frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + c(y_i,\phi) \right] \right\}. \] Here: - \( \theta_i \) is the canonical parameter - \( b(\theta_i) \) is the cumulant function - \( a(\phi) \) is the dispersion function - \( c(y_i,\phi) \) collects terms not involving \( \theta_i \) - \( w_i \) is the observation weight The **weighted log-likelihood** is therefore \[ \ell(\beta,\phi) = \sum_{i=1}^{n} w_i \left[ \frac{y_i\theta_i - b(\theta_i)}{a(\phi)} + c(y_i,\phi) \right]. \] This is the form used by both `glm()` and the Bayesian functions `glmb()`, `lmb()`, `rglmb()`, and `rlmb()`. --- ## 2.2 Special Case: Weighted Gaussian Linear Regression Consider the weighted Gaussian model \[ y_i \mid \beta,\phi \sim N(x_i^\top\beta,\;\phi/w_i), \qquad i=1,\dots,n. \] Let \(W = \mathrm{diag}(w_1,\dots,w_n)\) and define the precision \[ \tau = \frac{1}{\phi}. \] ### **Weighted log-likelihood** The log-likelihood (up to constants) is \[ \ell(\beta,\phi) = -\frac{1}{2\phi} \sum_{i=1}^{n} w_i (y_i - x_i^\top\beta)^2 \;-\; \frac{1}{2} \sum_{i=1}^{n} w_i \log(2\pi\phi). \] Equivalently, in precision form: \[ \ell(\beta,\tau) = \frac{\tau}{2} \sum_{i=1}^{n} w_i (y_i - x_i^\top\beta)^2 \;-\; \frac{1}{2} \sum_{i=1}^{n} w_i \log\left(\frac{2\pi}{\tau}\right). \] ### **Exponential-family identification** This matches the weighted exponential-family form with \[ \theta_i = \mu_i = x_i^\top\beta, \qquad b(\theta_i) = \tfrac{1}{2}\theta_i^2, \qquad a(\phi) = \phi, \] and \[ c(y_i,\phi) = -\frac{1}{2\phi}y_i^2 -\frac{1}{2}\log(2\pi\phi). \] ### **Weighted least squares representation** Base R implements weighted least squares via \[ y^\ast = W^{1/2} y, \qquad X^\ast = W^{1/2} X. \] The glmbayes functions use the same weighted likelihood, so all theoretical expressions above apply directly. ## Log-Concavity Properties in the Weighted Gaussian Case 1. **Concavity in \( \beta \)** For fixed \( \phi \), the negative weighted log-likelihood is \[ -\ell(\beta \mid \phi) = \frac{1}{2\phi} (y - X\beta)^{T} W (y - X\beta) + \text{const}, \] which is strictly convex whenever \( W \) has positive diagonal entries. Therefore the weighted log-likelihood is strictly concave in \( \beta \). 2. **Concavity in precision \( \tau = 1/\phi \)** Rewriting the weighted likelihood in terms of \( \tau \): \[ \ell(\beta,\tau) = -\frac{\tau}{2} (y - X\beta)^{T} W (y - X\beta) - \frac{1}{2} \left( \sum_{i=1}^{n} w_i \right) \log \tau + \text{const}, \] which is concave in \( \tau \) because the first term is linear in \( \tau \) and the second is concave. These properties guarantee: - unique posterior modes for both \( \beta \) and \( \tau \) - efficient iid sampling under Normal, Normal-Gamma, and Independent Normal-Gamma priors - stable two-block Gibbs sampling when the dispersion is updated separately The weighted Gaussian model thus provides the cleanest illustration of how weights enter the exponential-family structure and how they preserve log-concavity. ### 2.3 Special Case: Gamma Regression With a Log Link Consider the weighted Gamma regression model with log link: \[ y_i \mid \beta,\phi \;\sim\; \text{Gamma}\!\left( \text{shape}=\frac{1}{\phi},\; \text{scale}=\phi\,\mu_i \right), \qquad \mu_i = \exp(x_i^\top\beta), \] with observation weights \(w_i \ge 0\). Let the **precision** be \[ \tau = \frac{1}{\phi}. \] --- #### **Weighted log-likelihood** Up to an additive constant, the weighted log-likelihood is \[ \ell(\beta,\tau) = \sum_{i=1}^n w_i\left[ \tau\log\tau - \tau\log \mu_i + (\tau-1)\log y_i - \frac{\tau y_i}{\mu_i} \right]. \] Using \(\mu_i = \exp(x_i^\top\beta)\), this becomes \[ \ell(\beta,\tau) = \sum_{i=1}^n w_i\left[ \tau\log\tau - \tau x_i^\top\beta + (\tau-1)\log y_i - \tau y_i e^{-x_i^\top\beta} \right]. \] --- #### **Log‑concavity** For fixed \(\tau\): - the log-likelihood is **concave in \(\beta\)** because \(\mu_i = \exp(x_i^\top\beta)\) enters through a convex function and the Gamma log-likelihood is concave in the canonical parameter. For fixed \(\beta\): - the log-likelihood is **concave in \(\tau\)** because \(\tau\log\tau\), \((\tau-1)\log y_i\), and \(-\tau y_i/\mu_i\) combine into a concave function. These properties ensure: - unique posterior modes - stable envelope construction to support accept-reject sampling for \(\beta\) - valid accept–reject sampling for \(\tau\) ### Implications for Bayesian GLMs The Gamma log–link model thus shares the same structural advantages as the Gaussian and Poisson cases: * concavity in both \( \beta \) and \( v \), * unique posterior modes, * stable envelope construction, * efficient iid sampling for the dispersion parameter, * compatibility with two–block Gibbs sampling when desired. These properties make the Gamma log–link model an especially clean example of how exponential–family structure supports fast, reliable posterior simulation in glmbayes. ## 3. Gaussian Model With Unknown Dispersion For the Gaussian model \[ y = X\beta + \varepsilon,\qquad \varepsilon \sim N(0,\phi I_{n}), \] the likelihood is \[ p(y \mid \beta,\phi) \propto \phi^{-n/2}\exp\!\left(-\frac{S(\beta)}{2\phi}\right), \] where \[ S(\beta) = (y - X\beta)^{T}(y - X\beta) \] is the residual sum of squares. To support the different prior specification below, we use the Prior_Setup function to set and extract default settings for the needed parameters. ```{r Prior_Setup,results = "hide"} ps <- Prior_Setup(weight ~ group) x <- ps$x mu <- ps$mu V <- ps$Sigma y <- ps$y shape <- ps$shape rate <- ps$rate ``` ## 3. Priors for Gaussian Models With Unknown Dispersion We now specify priors for \((\beta,\tau)\) and derive the corresponding log-priors and log-posteriors. Throughout: - dispersion: \( \phi \) - precision: \( \tau = 1/\phi \) - weighted residual sum of squares: \[ \mathrm{RSS}(\beta) = \sum_{i=1}^{n} w_i (y_i - x_i^\top\beta)^2. \] --- ## 3.1 Prior on the Dispersion Only (Gamma Prior on Precision) We place a Gamma prior on the precision: \[ \tau \sim \mathrm{Gamma}(a_0,b_0), \] with log-prior \[ \log p(\tau) = (a_0 - 1)\log\tau - b_0\tau + \text{const}. \] Given fixed \(\beta\), the conditional posterior is \[ \tau \mid \beta,y \sim \mathrm{Gamma}\left( a_0 + \frac{n_{w}}{2}, \; b_0 + \frac{\mathrm{RSS}(\beta)}{2} \right), \] where \(n_{w} = \sum_i w_i\). --- To illustrate this, we fix \(\beta\) at the full-model coefficient vector from `Prior_Setup()` (`ps$coefficients`, the internal GLM maximum-likelihood estimate). We can then generate conditional Bayesian draws for the dispersion as below. This is essentially the equivalent of the classical function **gamma.dispersion()** function provided by the **MASS** package. ```{r glm_call_gamma_prior} out_rlmb <- rlmb( n = 1000, y = y, x = x, pfamily = dGamma(shape = shape, rate = rate, beta = ps$coefficients) ) ``` The posterior mean of the dispersion is then: ```{r glm_call_gamma_mean} mean(out_rlmb$dispersion) ``` ## 3.2 Joint Conjugate Normal–Gamma Prior The conjugate Normal--Gamma prior couples \(\beta\) and \(\tau\) [@LindleySmith1972; @OHaganForster2004]: \[ \beta \mid \tau \sim N(\mu_0,\;(\tau\Sigma_0)^{-1}), \qquad \tau \sim \mathrm{Gamma}(a_0,b_0). \] ### **Log-prior** \[ \log p(\beta,\tau) = (a_0 - 1)\log\tau - b_0\tau -\frac{\tau}{2}(\beta-\mu_0)^\top\Sigma_0^{-1}(\beta-\mu_0) +\text{const}. \] ### **Log-posterior** Combining with the Gaussian log-likelihood: \[ \log p(\beta,\tau \mid y) = \frac{\tau}{2}\mathrm{RSS}(\beta) -\frac{n_{w}}{2}\log\left(\frac{2\pi}{\tau}\right) +(a_0 - 1)\log\tau - b_0\tau -\frac{\tau}{2}(\beta-\mu_0)^\top\Sigma_0^{-1}(\beta-\mu_0) +\text{const}. \] This yields closed-form conditional posteriors: - \( \tau \mid y \) is Gamma - \( \beta \mid \tau,y \) is Normal allowing pure i.i.d. Monte Carlo sampling. --- The defining property of conjugacy is that the **posterior has the same posterior functional form**. In glmbayes, this model can be estimated as below. ```{r Conjugate_prior} ## Conjugate Normal_Gamma Prior lmb.D9_v2 <- lmb( weight ~ group, pfamily = dNormal_Gamma( ps$mu, Sigma_0 = ps$Sigma_0, shape = ps$shape, rate = ps$rate ) ) summary(lmb.D9_v2)$dispersion ``` ## 3.3 Independent Normal–Gamma Prior Here the prior factorizes: \[ \beta \sim N(\mu_0,\Sigma_0), \qquad \tau \sim \mathrm{Gamma}(a_0,b_0), \] so \(\beta\) and \(\tau\) are *a priori independent*. ### **Log-prior** \[ \log p(\beta,\tau) = -\frac{1}{2}(\beta-\mu_0)^\top\Sigma_0^{-1}(\beta-\mu_0) + (a_0 - 1)\log\tau - b_0\tau +\text{const}. \] ### **Log-posterior** \[ \log p(\beta,\tau \mid y) = \frac{\tau}{2}\mathrm{RSS}(\beta) -\frac{n_{w}}{2}\log\left(\frac{2\pi}{\tau}\right) -\frac{1}{2}(\beta-\mu_0)^\top\Sigma_0^{-1}(\beta-\mu_0) + (a_0 - 1)\log\tau - b_0\tau +\text{const}. \] A variation on this prior specification—using a **truncated Gamma prior** for the dispersion—is the model developed in **Chapter A07**, where: - the conditional posterior in \( \beta \) is *not conjugate*, - the conditional posterior in \( \tau = 1/\phi \) is *not conjugate*, and - the joint posterior remains **log‑concave**, enabling the **accept–reject envelope sampler**. This truncated‑Gamma, non‑conjugate Normal–Gamma model is therefore a direct extension of the independent Normal–Gamma prior described above, but requires the more advanced envelope‑based sampling strategy documented in Chapter A07. ### 3.3.1 Independent Normal Gamma Prior in glmbayes This prior is implemented in glmbayes using the below specification (but where a truncated gamma replaces the full gamma). The chunk below is **not evaluated** when the vignette is built. Precomputed results are read from **`inst/extdata/`** file **`Chapter11_Dobson_two_sampler.rds`**. To regenerate that file, run **`make_Chapter11_Dobson_two_sampler_output.R`** under **`data-raw/`**, from the package source root. ```{r chapter11-load-two-sampler} ch11_path <- system.file( "extdata", "Chapter11_Dobson_two_sampler.rds", package = "glmbayes" ) stopifnot(nzchar(ch11_path), file.exists(ch11_path)) ch11_saved <- readRDS(ch11_path) stopifnot( ncol(ch11_saved$gibbs_two_block$beta_out) == ncol(ps$x), nrow(ch11_saved$indep_norm_gamma$coefficients) == ch11_saved$indep_norm_gamma$n_draws ) ``` ```{r Independent_normal_gamma_prior, eval = FALSE} lmb.D9_v3 <- lmb(n = 10000, weight ~ group, dIndependent_Normal_Gamma( ps$mu, ps$Sigma, shape = ps$shape_ING, rate = ps$rate, max_disp_perc = 0.99, disp_lower = NULL, disp_upper = NULL ) ) summary(lmb.D9_v3)$coefficients summary(lmb.D9_v3)$dispersion sd(lmb.D9_v3$dispersion) ``` ```{r Independent_normal_gamma_loaded} inc <- ch11_saved$indep_norm_gamma coef_summ_iid <- data.frame( posterior_mean = colMeans(inc$coefficients), posterior_SD = apply(inc$coefficients, 2, sd), row.names = inc$coef_colnames ) coef_summ_iid cat("Dispersion: posterior mean =", mean(inc$dispersion), " SD =", sd(inc$dispersion), "\n") ``` ### 3.3.2 Two-Block Gibbs Sampling using glmbayes If iid sampling is not desired, or if the user wants to inspect mixing behavior, a simple two-block Gibbs sampler [@RobertCasella2004] may be used: 1. Sample \( \beta \mid \tau, y \) from a Normal distribution 2. Sample \( \tau \mid \beta, y \) from a Gamma distribution The following chunk is **not evaluated** during the vignette build (`eval = FALSE`). Stored draws are loaded from the same RDS as above. ```{r two_block_Gibbs_sampler, eval = FALSE} set.seed(180) dispersion2 <- ps$dispersion ## Burn-in for (i in 1:1000) { out1 <- rlmb( n = 1, y = y, x = x, pfamily = dNormal(mu = mu, Sigma = V, dispersion = dispersion2) ) out2 <- rlmb( n = 1, y = y, x = x, pfamily = dGamma(shape = shape, rate = rate, beta = out1$coefficients[1, ]) ) dispersion2 <- out2$dispersion } ## Sampling beta_out <- matrix(0, nrow = 10000, ncol = 2) disp_out <- rep(0, 10000) for (i in 1:10000) { out1 <- rlmb( n = 1, y = y, x = x, pfamily = dNormal(mu = mu, Sigma = V, dispersion = dispersion2) ) out2 <- rlmb( n = 1, y = y, x = x, pfamily = dGamma(shape = shape, rate = rate, beta = out1$coefficients[1, ]) ) beta_out[i, ] <- out1$coefficients[1, 1:2] disp_out[i] <- out2$dispersion } colMeans(beta_out) sd(beta_out[, 1]) sd(beta_out[, 2]) mean(disp_out) sd(disp_out) ``` ```{r two_block_Gibbs_loaded} gb <- ch11_saved$gibbs_two_block beta_out <- gb$beta_out disp_out <- gb$disp_out colMeans(beta_out) sd(beta_out[, 1]) sd(beta_out[, 2]) mean(disp_out) sd(disp_out) ``` ### 3.3.3 Comparison of the Two Samplers Under default settings, both the truncated‑Gamma iid sampler and the two‑block Gibbs sampler produce nearly identical posterior summaries for the regression coefficients and the dispersion parameter. To make this explicit, the table below reports posterior means and standard deviations from: * the iid sampler using the independent Normal–Gamma prior (same specification as `lmb(..., dIndependent_Normal_Gamma(...), n = 10000)` above), and * the two‑block Gibbs sampler constructed from alternating Normal and Gamma updates. These values illustrate that both samplers target essentially the same posterior distribution (although the former has a truncated prior), with only minor Monte Carlo variation. ```{r chapter11-sampler-comparison-table, echo = FALSE} inc <- ch11_saved$indep_norm_gamma gb <- ch11_saved$gibbs_two_block cn <- inc$coef_colnames iid_m <- colMeans(inc$coefficients) iid_s <- apply(inc$coefficients, 2, sd) gib_m <- colMeans(gb$beta_out) gib_s <- apply(gb$beta_out, 2, sd) cmp <- data.frame( Parameter = c(cn, "Dispersion"), iid_Mean = c(iid_m, mean(inc$dispersion)), iid_SD = c(iid_s, sd(inc$dispersion)), Gibbs_Mean = c(gib_m, mean(gb$disp_out)), Gibbs_SD = c(gib_s, sd(gb$disp_out)) ) knitr::kable( cmp, digits = 4, caption = paste( "Dobson plant weight: independent Normal-Gamma (iid) vs", "two-block Gibbs" ) ) ``` The difference in the iid and Gibbs SD for the dispersion is largely due to the difference between a truncated and an unrestricted gamma prior. However, the iid sampler has two clear advantages: * **Independent draws**: no autocorrelation and no loss of effective sample size * **No tuning or burn‑in**: each draw is a direct sample from the posterior The Gibbs sampler mixes well for this model class, but it inevitably shows mild autocorrelation and convergence diagnostics are beyond the scope of this vignette. For most users, the iid sampler is easier to work with and safer as no convergence diagnostics are needed. For readers who want the full mathematical details of the iid sampler, see **Chapter A07**, which documents the complete envelope construction and accept‑reject algorithm for the independent Normal–Gamma prior. ## 4. Gamma Regression Models With Unknown Dispersion We now specify priors for \((\beta,\tau)\) and derive the corresponding log‑priors and log‑posteriors, using the **carinsca** data [@BaileySimon1960; @statsciCarinsca; @McCullagh1989]. Throughout: - dispersion: \(\phi\) - precision: \(\tau = 1/\phi\) - mean: \(\mu_i = \exp(x_i^\top\beta)\) - weighted log-likelihood: \(\ell(\beta,\tau)\) from Section 2.3 ```{r Carinsca} data(carinsca) carinsca$Merit <- ordered(carinsca$Merit) carinsca$Class <- factor(carinsca$Class) oldopt <- options(contrasts = c("contr.treatment", "contr.treatment")) Claims=carinsca$Claims Insured=carinsca$Insured Merit=carinsca$Merit Class=carinsca$Class Cost=carinsca$Cost Claims_Adj<-Claims/1000 glm.carinsca <- glm(Cost / Claims ~ Merit + Class, family = Gamma(link = "log"), weights = Claims_Adj, x = TRUE) ``` ```{r Prior_Setup_gamma,results = "hide"} ps <- Prior_Setup( Cost / Claims ~ Merit + Class, family = Gamma(link = "log"), weights = Claims_Adj ) mu=ps$mu V=ps$Sigma shape <- ps$shape rate <- ps$rate x <- ps$x y <- ps$y m<-nrow(x) p<-ncol(x) ## Starting dispersion for beta | tau in the two-block Gibbs sampler: same quasi-likelihood ## estimate as used elsewhere for this Carinsca Gamma GLM (not the Dobson Gaussian ps). dispersion2 <- gamma.dispersion(glm.carinsca) options(oldopt) ``` --- ### 4.1 Prior on the Dispersion Only (Gamma Prior on Precision) We place a Gamma prior on the precision: \[ \tau \sim \text{Gamma}(a_0,b_0), \] with log‑prior \[ \log p(\tau) = (a_0 - 1)\log\tau - b_0\tau + \text{const}. \] Given fixed \(\beta\), the conditional posterior is \[ \log p(\tau \mid \beta,y) = \ell(\beta,\tau) + (a_0 - 1)\log\tau - b_0\tau + \text{const}. \] This posterior is **log‑concave in \(\tau\)**, enabling envelope‑based accept–reject sampling in the Gamma regression case. ```{r glm_call_gamma_prior2} ## Carinsca Gamma GLM (already fitted in the Carinsca chunk for gamma.dispersion etc.) gamma.dispersion(glm.carinsca) out2 <- rglmb( n = 1000, y = y, x = x, family =Gamma(link="log"), pfamily = dGamma(shape = shape, rate = rate, beta = ps$coefficients), weights = Claims_Adj ) mean(out2$dispersion) ``` ### 4.2 Independent Normal–Gamma Prior Here the priors factorize: \[ \beta \sim N(\mu_0,\Sigma_0), \qquad \tau \sim \text{Gamma}(a_0,b_0), \] so \(\beta\) and \(\tau\) are **a priori independent**. --- #### **Log‑prior** \[ \log p(\beta,\tau) = -\tfrac12(\beta-\mu_0)^\top\Sigma_0^{-1}(\beta-\mu_0) + (a_0 - 1)\log\tau - b_0\tau + \text{const}. \] --- #### **Log‑posterior** \[ \log p(\beta,\tau \mid y) = \ell(\beta,\tau) -\tfrac12(\beta-\mu_0)^\top\Sigma_0^{-1}(\beta-\mu_0) + (a_0 - 1)\log\tau - b_0\tau + \text{const}. \] In Gamma regression: - the conditional posterior in \(\beta\) is **not conjugate** - the conditional posterior in \(\tau\) is **not conjugate** - the joint posterior remains **log‑concave** This in theory enables two sampling strategies for the Gamma regression model: 1. **Two‑block Gibbs sampler** - \( \beta \mid \tau, y \): sampled using a Normal approximation (via `rglmb`) - \( \tau \mid \beta, y \): sampled using a Gamma‑based quasi‑likelihood update This is the **only fully implemented method** for Gamma regression in the current version of *glmbayes*. 2. **Accept–reject envelope sampler (Chapter A07)** Chapter A07 develops an accept–reject sampler using a **truncated Gamma prior** for \( \tau \) and a joint envelope over \((\beta,\tau)\). **However, this method is implemented only for the Gaussian family.** Extending the envelope‑based sampler to the Gamma regression model would require additional mathematical development and corresponding modifications to the code. The truncated‑Gamma independent Normal–Gamma model in **Chapter A07** therefore applies exclusively to the Gaussian case, while Gamma regression currently relies on the two‑block Gibbs sampler described above. Implementation details: alternate `rglmb` draws for \(\beta \mid \tau, y\) using `dNormal(..., dispersion = dispersion2)` with draws for \(\tau \mid \beta, y\) using `dGamma(...)` on the Gamma regression likelihood, updating `dispersion2` after each \(\tau\) step. A **1000‑iteration burn‑in** is followed by **1000 stored iterations** (coefficient matrix `beta_out`, dispersion vector `disp_out`, and accept–reject counts `iters_out`). The same logic is written to **`inst/extdata/Chapter11_Carinsca_gamma_gibbs.rds`** by **`data-raw/make_Chapter11_Carinsca_gamma_gibbs_output.R`** (**`set.seed(190)`**). The chunk **`Block_Gibbs_gamma_Regression`** below is **not evaluated** when the vignette is built (`eval = FALSE`). Stored draws are loaded from **`inst/extdata/Chapter11_Carinsca_gamma_gibbs.rds`** in **`chapter11-load-carinsca-gamma`**, and summaries are produced in **`Block_Gibbs_gamma_loaded`**. ```{r chapter11-load-carinsca-gamma} ch11_cg_path <- system.file( "extdata", "Chapter11_Carinsca_gamma_gibbs.rds", package = "glmbayes" ) stopifnot(nzchar(ch11_cg_path), file.exists(ch11_cg_path)) ch11_cg <- readRDS(ch11_cg_path) stopifnot( ncol(ch11_cg$gibbs_gamma$beta_out) == ncol(ps$x), length(ch11_cg$gibbs_gamma$disp_out) == nrow(ch11_cg$gibbs_gamma$beta_out) ) ``` ```{r Block_Gibbs_gamma_Regression, eval = FALSE} set.seed(190) dispersion2 <- gamma.dispersion(glm.carinsca) suppressWarnings( suppressMessages( for(i in 1:1000){ ## --- Block 1: Regression coefficients --- out1 <- rglmb( n = 1, y = y, x = x, family = Gamma(link="log"), pfamily = dNormal(mu = mu, Sigma = V, dispersion = dispersion2), weights = Claims_Adj ) ## --- Block 2: Dispersion (quasi-likelihood sampler) --- out2 <- rglmb( n = 1, y = y, x = x, family = Gamma(link="log"), pfamily = dGamma(shape = shape, rate = rate, beta = out1$coefficients[1,]), weights = Claims_Adj ) ## --- SCALE dispersion for the next beta update --- ## Convert quasi (from out2) to MLE for consistency dispersion2 <- out2$dispersion ##* ((m - p)/m) } ) ) ## Run 1000 additional iterations and store output beta_out <- matrix(0, nrow = 1000, ncol = ncol(x)) disp_out <- rep(0, 1000) iters_out <- rep(0, 1000) suppressWarnings( suppressMessages( for (i in 1:1000) { out1 <- rglmb( n = 1, y = y, x = x, family = Gamma(link="log"), pfamily = dNormal(mu = mu, Sigma = V, dispersion = dispersion2), weights = Claims_Adj ) out2 <- rglmb( n = 1, y = y, x = x, family = Gamma(link="log"), pfamily = dGamma(shape = shape, rate = rate, beta = out1$coefficients[1,]), weights = Claims_Adj ) dispersion2 <- out2$dispersion ##* ((m - p) / m) beta_out[i, ] <- out1$coefficients[1, seq_len(ncol(x))] disp_out[i] <- out2$dispersion ##* ((m - p) / m) iters_out[i]<-out2$iters } ) ) ## Review output beta_mean <- colMeans(beta_out) beta_sd <- apply(beta_out, 2, sd) beta_tlike <- beta_mean / beta_sd # analogous to GLM t-values bayes_coef_table <- data.frame( Estimate = beta_mean, Std.Error = beta_sd, t.like = beta_tlike ) bayes_coef_table mean(disp_out) mean(iters_out) ``` ```{r Block_Gibbs_gamma_loaded} cg <- ch11_cg$gibbs_gamma beta_out <- cg$beta_out disp_out <- cg$disp_out iters_out <- cg$iters_out beta_mean <- colMeans(beta_out) beta_sd <- apply(beta_out, 2, sd) beta_tlike <- beta_mean / beta_sd bayes_coef_table <- data.frame( Estimate = beta_mean, Std.Error = beta_sd, t.like = beta_tlike ) bayes_coef_table mean(disp_out) mean(iters_out) ``` ### 4.3 Conditional Sampling for the Dispersion Parameter Chapter A06 provides a detailed derivation of the conditional posterior for the precision parameter \( v = 1/\phi \) in Gamma regression. The posterior can be written in the form \[ \ell(v) = ( \text{shape2} - 1 ) \log v - \text{rate1} \, v + \text{testfunc}(v) + \text{const}, \] where \( \text{testfunc}(v) \) is a concave remainder term capturing the curvature contributed by the Gamma likelihood. This decomposition enables a **tangent–envelope accept–reject sampler**, in which: 1. A Gamma proposal distribution is constructed using the linear part \[ ( \text{shape2} - 1 ) \log v - \text{rate2} \, v, \] where \( \text{rate2} = \text{rate1} - g^{*} \) and \( g^{*} \) is the derivative of \( \text{testfunc}(v) \) at the anchor point \( v^{*} \). 2. A proposed value \( v_{\text{prop}} \) is accepted with probability \[ \exp\!\left[ \text{testfunc}(v_{\text{prop}}) - \text{testfunc}(v^{*}) - g^{*}(v_{\text{prop}} - v^{*}) \right]. \] Because \( \text{testfunc}(v) \) is concave, the tangent line at \( v^{*} \) always lies above it, ensuring a valid rejection sampler. ## 8. Summary This chapter introduced Bayesian estimation of Gaussian models with unknown dispersion [@Gelman2013; @OHaganForster2004], including: • Gamma priors on the precision \( \tau \) • Joint Normal–Gamma priors • Independent Normal–Gamma priors • iid sampling via envelope methods • two‑block Gibbs sampling The log‑concavity of the Gaussian likelihood in both \( \beta \) and \( \tau \) ensures that all samplers are stable, efficient, and theoretically well‑behaved. ## A1: Posterior Distribution Details for Conjugate Normal-Gamma prior ### A1.1 Marginal Posterior for \( \tau \) Integrating out \( \beta \) yields a Gamma posterior **independent of \( \beta \)**: \[ \tau \mid y \sim \mathrm{Gamma}\!\left( a_{0} + \frac{n}{2},\; b_{0} + \frac{1}{2} S_{\text{marg}} \right), \] where \[ S_{\text{marg}} = (y - X\mu_{0})^{T} \left(I + X\Sigma_{0}X^{T}\right)^{-1} (y - X\mu_{0}). \] Here \( \Sigma_{0} \) is the same matrix as in Section 3.2 and in A1.2 (so that \( \Sigma_{0}^{-1} \) is the prior precision factor in \( -\frac{\tau}{2}(\beta-\mu_{0})^{\mathsf T}\Sigma_{0}^{-1}(\beta-\mu_{0}) \)). Let \( \hat\beta = (X^{\mathsf T}X)^{-1} X^{\mathsf T} y \) be the ordinary least squares estimator and \( \mathrm{RSS} = (y - X\hat\beta)^{\mathsf T}(y - X\hat\beta) \). Write \( G = X^{\mathsf T}X \) (invertible in the usual full-rank setup). The **same** scalar \( S_{\text{marg}} \) can be written as \[ S_{\text{marg}} = \mathrm{RSS} + (\hat\beta - \mu_{0})^{\mathsf T} \left(\Sigma_{0} + G^{-1}\right)^{-1} (\hat\beta - \mu_{0}). \] Equivalently, expanding \( (I + X\Sigma_{0}X^{\mathsf T})^{-1} \) with the Sherman--Morrison--Woodbury formula gives the same middle matrix in the precision form \[ G - G\left(\Sigma_{0}^{-1} + G\right)^{-1}G = \left(\Sigma_{0} + G^{-1}\right)^{-1}. \] In particular, if \( \mu_{0} = \hat\beta \) then the second term vanishes and \( S_{\text{marg}} = \mathrm{RSS} \). **Step-by-step (Woodbury + least-squares split).** Set \(e = y - X\hat\beta\), \(\delta = \hat\beta - \mu_{0}\), so \(y - X\mu_{0} = e + X\delta\). Normal equations give \(X^{\mathsf T}e = 0\), hence \(\|e + X\delta\|^{2} = \mathrm{RSS} + \delta^{\mathsf T}G\delta\) and \(X^{\mathsf T}(e + X\delta) = G\delta\). Woodbury says \((I + X\Sigma_{0}X^{\mathsf T})^{-1} = I - X(\Sigma_{0}^{-1}+G)^{-1}X^{\mathsf T}\). Substitute \(y-X\mu_{0}=e+X\delta\) into the quadratic form for \(S_{\text{marg}}\); cross terms vanish because \(e^{\mathsf T}X=0\), and the algebra reduces to \(\mathrm{RSS} + \delta^{\mathsf T}\bigl(G - G(\Sigma_{0}^{-1}+G)^{-1}G\bigr)\delta\), which equals the display above. For observation weights \( w_{i} \) (Section 2.2), define \(G = X^{\mathsf T}WX\), \(\hat\beta = G^{-1} X^{\mathsf T}Wy\), and \(\mathrm{RSS} = (y - X\hat\beta)^{\mathsf T}W(y - X\hat\beta)\). The **RSS + \((\hat\beta-\mu_{0})^{\mathsf T}(\Sigma_{0}+G^{-1})^{-1}(\hat\beta-\mu_{0})\)** identity is unchanged (same proof in whitened coordinates \(\tilde y = W^{1/2}y\), \(\tilde X = W^{1/2}X\)). The **first** display becomes \((y-X\mu_{0})^{\mathsf T}W^{1/2}(I + W^{1/2}X\Sigma_{0}X^{\mathsf T}W^{1/2})^{-1}W^{1/2}(y-X\mu_{0})\) for any factor \(W^{1/2}\) with \(W^{1/2}(W^{1/2})^{\mathsf T}=W\) (e.g.\ diagonal \(\sqrt{w_i}\)). This closed‑form marginal is what makes the Normal–Gamma prior so computationally attractive: sampling \( \tau \) requires **only a single Gamma draw**. ### A1.2 Conditional Posterior for \( \beta \) Given \( \tau \), the posterior for \( \beta \) is Gaussian: \[ \beta \mid \tau,y \sim N(m_{\text{post}}, V_{\text{post}}), \] with \[ V_{\text{post}} = \left(\Sigma_{0}^{-1} + \tau X^{T}X\right)^{-1}, \qquad m_{\text{post}} = V_{\text{post}} \left(\Sigma_{0}^{-1}\mu_{0} + \tau X^{T}y\right). \] ### A1.3 Sampling the Posterior Because the Normal--Gamma prior is *conjugate* to the Gaussian likelihood [@OHaganForster2004], the **joint posterior factorizes into two standard distributions**: \[ p(\tau,\beta \mid y) = p(\tau \mid y)\, p(\beta \mid \tau, y). \] This structure means **no MCMC is required**. Posterior draws are obtained by **independent Monte Carlo sampling**, not by a Gibbs chain. --- #### Step 1: Draw the marginal posterior for \( \tau \) The dispersion (precision) parameter has a closed‑form marginal posterior: \[ \tau \mid y \sim \mathrm{Gamma}\!\left( a_{0} + \frac{n}{2},\; b_{0} + \frac{1}{2} S(\hat\beta_{\text{GLS}}) \right), \] where \[ S(\hat\beta_{\text{GLS}}) = (y - X\hat\beta_{\text{GLS}})^{T}(y - X\hat\beta_{\text{GLS}}) \] and \( \hat\beta_{\text{GLS}} \) is the generalized least‑squares estimator under the prior precision. Each draw of \( \tau \) is **i.i.d.** from this Gamma distribution. --- #### Step 2: Draw \( \beta \) conditionally on each sampled \( \tau \) Given a sampled value \( \tau^{(s)} \), the regression coefficients have a Normal posterior: \[ \beta \mid \tau^{(s)}, y \sim N\!\left(m_{\text{post}}^{(s)},\, V_{\text{post}}^{(s)}\right), \] with \[ V_{\text{post}}^{(s)} = \left(\Sigma_{0}^{-1} + \tau^{(s)} X^{T}X\right)^{-1}, \] \[ m_{\text{post}}^{(s)} = V_{\text{post}}^{(s)} \left(\Sigma_{0}^{-1}\mu_{0} + \tau^{(s)} X^{T}y\right). \] Each draw of \( \beta \) is **conditionally independent** given its corresponding \( \tau^{(s)} \). --- ### Result: Pure i.i.d. Monte Carlo Sampling The conjugate posterior allows: - **direct sampling** of \( \tau \) from its marginal Gamma distribution - **direct sampling** of \( \beta \) from its conditional Normal distribution No Markov chain is constructed, and **no Gibbs or Metropolis steps** are involved. All posterior samples are **independent**, making the conjugate Normal–Gamma model one of the simplest Bayesian regression frameworks for simulation. ## A2. Detailed sampling procedured for the Independent Normal-Gamma prior ### A2.1 iid Sampling Under Truncated Gamma Priors When iid sampling is requested, glmbayes uses a truncated Gamma prior for the precision parameter \( \tau \). The truncation is introduced for numerical stability and does not change the basic Normal-Gamma structure. The iid sampler for this prior family is based on the envelope accept-reject method described in: **Chapter A07 — Accept-Reject Sampling for Gaussian Regression Models with Independent Normal-Gamma Priors** That vignette explains the full simulation procedure, including: * how the dispersion parameter is truncated, * how the envelope is constructed across dispersion values, * how the mixture of truncated normal proposals for \( \beta \) is formed, * how the Gamma proposal for \( \tau \) is calibrated. The key point for users is that the posterior remains log-concave, so the accept-reject sampler produces independent draws for \( (\beta, \tau) \). There is no Markov chain, no burn-in, and no autocorrelation. ## A3. Posterior mean of \( \beta \) and marginal covariance (Zellner-type prior) This subsection records closed-form **posterior moments** for \( \beta \) in the **conjugate Normal-Gamma** Gaussian regression model, in the same weighted likelihood notation as Section 2.2 and in the **prior weight / effective prior sample size** notation used by `Prior_Setup()`. ### A3.1 Weighted likelihood and prior precision Assume \[ y \mid \beta,\phi \sim N\!\left(X\beta,\;\phi W^{-1}\right), \qquad \tau = 1/\phi, \] with \(W = \mathrm{diag}(w_{1},\dots,w_{n})\) and \(w_i \ge 0\). Define the **effective sample size** \[ n_{w} = \sum_{i=1}^{n} w_i \] (the quantity `sum(wt)` in `rNormalGammaReg()`). Write the **weighted Gram matrix** as \[ G = X^{\mathsf T} W X. \] Let \(\hat\beta\) denote the **weighted least squares** estimator \[ \hat\beta = G^{-1} X^{\mathsf T} W y \] and the **weighted residual sum of squares** \[ \mathrm{RSS} = (y - X\hat\beta)^{\mathsf T} W (y - X\hat\beta). \] Take a **Zellner-type** conditional Normal prior on \(\beta\) given \(\tau\) that is proportional to the data precision: \[ \beta \mid \tau \sim N\!\left(\mu_{0},\; \tau^{-1}\left(\frac{pwt}{1-pwt}G\right)^{-1}\right), \qquad 0 < pwt < 1, \] equivalently prior **precision** \(\tau\, (pwt/(1-pwt))\, G\). This is the structure implied by `Prior_Setup()` together with `dNormal_Gamma(..., Sigma_0, ...)` when \(\Sigma_0 = ((1-pwt)/pwt)\, G^{-1}\) (Zellner \(g\)-prior relative to \(G^{-1}\); see Section 3.2). Define the **effective prior sample size** \(n_{\mathrm{prior}} > 0\) by \[ pwt = \frac{n_{\mathrm{prior}}}{n_{\mathrm{prior}} + n_{w}}, \qquad 1 - pwt = \frac{n_{w}}{n_{\mathrm{prior}} + n_{w}}. \] This is the relationship implemented when `Prior_Setup()` is called with scalar `n_prior` (`n_prior` and `n_effective` in the function). ### A3.2 Posterior mean of \( \beta \) Combining likelihood precision \(\tau G\) with prior precision \(\tau (pwt/(1-pwt)) G\) gives posterior precision **given \(\tau\)**: \[ \tau G + \tau \frac{pwt}{1-pwt} G = \frac{\tau}{1-pwt}\, G. \] Hence \[ \mathrm{Cov}(\beta \mid \tau, y) = \frac{1-pwt}{\tau}\, G^{-1} = \frac{n_{w}}{n_{\mathrm{prior}} + n_{w}}\,\frac{1}{\tau}\, G^{-1}. \] The **posterior mean** \(m_{\mathrm{post}} = E(\beta \mid \tau,y)\) solves the same normal equations as the **precision-weighted average** of \(\hat\beta\) and \(\mu_{0}\), and simplifies to \[ \boxed{ E(\beta \mid y) = (1-pwt)\,\hat\beta + pwt\,\mu_{0} = \frac{n_{w}}{n_{\mathrm{prior}} + n_{w}}\,\hat\beta + \frac{n_{\mathrm{prior}}}{n_{\mathrm{prior}} + n_{w}}\,\mu_{0}. } \] The right-hand side **does not depend on \(\tau\)** (for this prior), so \(E(\beta \mid y) = E(\beta \mid \tau,y)\). ### A3.3 Marginal covariance of \( \beta \) and the Gamma parameters Marginalize \(\tau\) with \[ \tau \mid y \sim \mathrm{Gamma}(a_{n}, b_{n}) \] in the **shape-rate** parameterization \(p(\tau) \propto \tau^{a_{n}-1} e^{-b_{n}\tau}\) used in `rNormalGammaReg()` (`R::rgamma(shape, 1/rate)`), so \(E(\tau^{-1} \mid y) = b_{n}/(a_{n}-1)\) for \(a_{n} > 1\). Let \[ V_{n} = (1-pwt)\, G^{-1} = \frac{n_{w}}{n_{\mathrm{prior}} + n_{w}}\, G^{-1}, \] so \(\mathrm{Cov}(\beta \mid \tau,y) = \tau^{-1} V_{n}\). Because \(E(\beta \mid \tau,y)\) is constant in \(\tau\), the **law of total covariance** gives \[ \boxed{ \mathrm{Cov}(\beta \mid y) = E(\tau^{-1} \mid y)\, V_{n} = \frac{b_{n}}{a_{n}-1}\, \frac{n_{w}}{n_{\mathrm{prior}} + n_{w}}\, G^{-1}, \qquad a_{n} > 1. } \] In the sampler, \[ a_{n} = a_{0} + \frac{n_{w}}{2}, \qquad b_{n} = b_{0} + \frac{1}{2} S_{\mathrm{marg}}, \] with \(S_{\mathrm{marg}}\) the marginal sum of squares from Section A1.1. In the notation of this section (\(G = X^{\mathsf T}WX\), weighted least squares \(\hat\beta = G^{-1} X^{\mathsf T}Wy\), weighted \(\mathrm{RSS} = (y - X\hat\beta)^{\mathsf T}W(y - X\hat\beta)\), and \(\Sigma_{0}\) as in Section 3.2), \[ S_{\mathrm{marg}} = \mathrm{RSS} + (\hat\beta - \mu_{0})^{\mathsf T} \left(\Sigma_{0} + G^{-1}\right)^{-1} (\hat\beta - \mu_{0}). \] This is the same scalar as the augmented residual sum of squares `S` produced by `rNormal_reg.wfit()` (Section A1). Under the Zellner \(g\)-prior of A3.1, the prior **covariance** is \(\Sigma_{0} = ((1-pwt)/pwt)\,G^{-1}\), so \(\Sigma_{0} + G^{-1} = G^{-1}/pwt\) and \((\Sigma_{0} + G^{-1})^{-1} = pwt\, G\). Equivalently, \[ S_{\mathrm{marg}} = \mathrm{RSS} + pwt\, (\hat\beta - \mu_{0})^{\mathsf T} G (\hat\beta - \mu_{0}). \] The **structural** point is that the prior--data mismatch enters only through the scalar factor **\(pwt\)** multiplying a nonnegative quadratic in \(\hat\beta-\mu_{0}\). Hence, holding the data and \(\mu_{0}\) fixed so that \((\hat\beta-\mu_{0})^{\mathsf T}G(\hat\beta-\mu_{0})\) stays finite, \[ pwt \to 0^{+} \quad\Longrightarrow\quad pwt\,(\hat\beta-\mu_{0})^{\mathsf T}G(\hat\beta-\mu_{0}) \to 0, \] and \(S_{\mathrm{marg}} \to \mathrm{RSS}\): weak prior weight **washes out** the gap even when \(\mu_{0} \neq \hat\beta\). Equivalently, \(pwt\to 0\) whenever \(n_{\mathrm{prior}}/(n_{\mathrm{prior}}+n_{w})\to 0\)---for example \(n_{\mathrm{prior}}\to 0\) with \(n_{w}\) held fixed, or \(n_{w}\to\infty\) with \(n_{\mathrm{prior}}\) held fixed. The **exact** cancellation \(S_{\mathrm{marg}} = \mathrm{RSS}\) also holds for any \(pwt\) when \(\mu_{0} = \hat\beta\) (for example `Prior_Setup(..., intercept_source = "full_model", effects_source = "full_model")`). For fixed \(pwt>0\) and \(\mu_{0} \neq \hat\beta\), the second term is strictly positive unless \(\hat\beta-\mu_{0}\) lies in the null space of \(G\). Under standard regularity, with fixed \(n_{\mathrm{prior}}\) and \(n_{w}\to\infty\), one has \(pwt = O(1/n_{w})\) while \((\hat\beta-\mu_{0})^{\mathsf T}G(\hat\beta-\mu_{0})\) is typically \(O_{p}(n_{w})\), so the mismatch term is often \(O_{p}(1)\): the **additive** difference \(S_{\mathrm{marg}}-\mathrm{RSS}\) need not vanish at finite \(n_{w}\), but \(S_{\mathrm{marg}}/n_{w}\) and \(\mathrm{RSS}/n_{w}\) share the same leading order. **Limit \(pwt\to 0\) and coefficient covariance.** Under the Zellner decomposition above, **\(pwt\to 0^{+}\)** forces **\(S_{\mathrm{marg}}\to \mathrm{RSS}\)** (for fixed data and \(\mu_{0}\), with a finite mismatch quadratic), not only when \(\mu_{0}=\hat\beta\). Since \(b_{n} = b_{0} + \tfrac{1}{2} S_{\mathrm{marg}}\), the marginal Gamma rate then limits to **\(b_{0} + \tfrac{1}{2}\mathrm{RSS}\)** whenever \(b_{0}\) is held fixed along the path; \(E(\tau^{-1}\mid y)=b_{n}/(a_{n}-1)\) therefore ties to **RSS** through \(b_{n}\) in the same limit. At the same time, **\((1-pwt)\,G^{-1}\to G^{-1}\)** from A3.2--A3.3. Thus \(\mathrm{Cov}(\beta \mid \tau,y) = \tau^{-1}(1-pwt)\,G^{-1}\) tends to \(\tau^{-1}G^{-1}\) (likelihood-only WLS covariance given \(\tau\)), and \[ \mathrm{Cov}(\beta \mid y) = \frac{b_{n}}{a_{n}-1}\,(1-pwt)\,G^{-1}, \qquad a_{n}>1, \] has both **\(b_{n}\)** driven by **\(S_{\mathrm{marg}}\to \mathrm{RSS}\)** and the matrix prefactor \((1-pwt)\,G^{-1}\to G^{-1}\): weak \(pwt\) removes prior shrinkage of precision **and** drives the residual sum of squares in the \(\tau\) update back to the **data** \(\mathrm{RSS}\). The usual \(\sigma^{2}(X^{\mathsf T}WX)^{-1}\) limit is recovered when, in addition, \(a_{n}\) and \(b_{0}\) follow the weak-prior path discussed below (e.g.\ formal \(n_{\mathrm{prior}}\to 0\) with \(n_{w}\) fixed, holding \(a_{n}>1\)). If \(pwt\) is driven to zero by \(n_{w}\to\infty\) at fixed \(n_{\mathrm{prior}}\), then \((1-pwt)\to 1\) **and** \(G^{-1}\) typically scales as \(O(1/n_{w})\), so \(\mathrm{Cov}(\beta\mid y)\) shrinks with more weighted data even though the prior contributes negligibly to precision. With `Prior_Setup()` defaults for the Gamma on \(\tau\), \(a_{0} = n_{\mathrm{prior}}/2\) in the pre-calibration Step 9 rate scaling, so \(a_{n} = (n_{\mathrm{prior}} + n_{w})/2\); the prior rate is \(b_{0} = (n_{\mathrm{prior}}/2)\, d\) with `dispersion` \(d\) returned by `Prior_Setup()` (see Section 3.2 and `?Prior_Setup`). The **marginal** posterior of \(\beta\) is **multivariate \(t\)**; the matrix above is its **covariance** (not the scale matrix in every textbook parameterization of the \(t\) density). #### `Prior_Setup()` dispersion and the Gamma rates \(b_{0}\), \(b_{n}\) For Gaussian models with scalar `n_prior`, `Prior_Setup()` sets the prior Gamma **rate** to \(\mathrm{rate} = (n_{\mathrm{prior}}/2)\, d\), where \(d\) is the scalar `dispersion` returned by the function. In the notation above, \(b_{0} = \mathrm{rate}\). The **posterior** Gamma rate is \(b_{n} = b_{0} + \frac{1}{2} S_{\mathrm{marg}}\). Thus a larger returned \(d\) increases \(b_{0}\) and, holding \(S_{\mathrm{marg}}\) fixed, would increase \(b_{n}\). In practice, for **`dNormal_Gamma`** the matrix passed to the sampler is **`Sigma_0`** (dispersion-free prior covariance from `Prior_Setup()`), while the returned **`dispersion`** calibrates the Gamma rate and the coefficient-scale **`Sigma`**; changing \(d\) rescales the Normal--Gamma coupling of \(\beta\) given \(\tau\); the augmented least-squares step then produces a **different** \(S_{\mathrm{marg}}\) and posterior mean \(m_{\mathrm{post}}\) as well. #### Special case: \(b_{0} = \frac{1}{2}(n_{\mathrm{prior}}/n_{w})\,S_{\mathrm{marg}}\) \[ \frac{S_{\mathrm{marg}}}{n_{w}} \] This choice is **not** the `Prior_Setup()` default; it is a **hypothetical** reparameterization that sets the prior Gamma rate proportional to the **same** marginal sum of squares that enters the likelihood fragment, with weight \(n_{\mathrm{prior}}/n_{w}\) on the prior piece relative to the data piece (see the discussion in the main text). Because \(S_{\mathrm{marg}}\) depends on \(y\) (and on \(\mu_{0}\) versus \(\hat\beta\)), such a \(b_{0}\) is **data-dependent** unless used only in stylized derivations. Assume \(a_{0} = n_{\mathrm{prior}}/2\) (nominal rate scaling). Then \(a_{n} = (n_{\mathrm{prior}} + n_{w})/2\) and, provided \(a_{n} > 1\) (equivalently \(n_{\mathrm{prior}} + n_{w} > 2\)), \[ a_{n} - 1 = \frac{n_{\mathrm{prior}} + n_{w} - 2}{2}. \] Take \[ b_{0} = \frac{1}{2}\,\frac{n_{\mathrm{prior}}}{n_{w}}\, S_{\mathrm{marg}}, \qquad b_{n} = b_{0} + \frac{1}{2} S_{\mathrm{marg}} = \frac{1}{2}\, S_{\mathrm{marg}}\left(1 + \frac{n_{\mathrm{prior}}}{n_{w}}\right) = \frac{S_{\mathrm{marg}}}{2}\,\frac{n_{w} + n_{\mathrm{prior}}}{n_{w}}. \] Then the posterior expectation of \(\tau^{-1} = \phi = \sigma^{2}\) is \[ \frac{b_{n}}{a_{n}-1} = \frac{S_{\mathrm{marg}}(n_{w} + n_{\mathrm{prior}})/n_{w}} {(n_{\mathrm{prior}} + n_{w} - 2)}. \] Under the Zellner prior of A3.1, \(V_{n} = (1-pwt)\, G^{-1} = \dfrac{n_{w}}{n_{\mathrm{prior}} + n_{w}}\, G^{-1}\). One may write \(\mathrm{Cov}(\beta \mid y)\) as the product of \(b_{n}/(a_{n}-1)\) and \(V_{n}\) **without** combining the \((n_{\mathrm{prior}} + n_{w})\) terms: \[ \boxed{ \mathrm{Cov}(\beta \mid y) = \frac{b_{n}}{a_{n}-1}\, V_{n} = \dfrac{S_{\mathrm{marg}}\,(n_{\mathrm{prior}} + n_{w})/n_{w}} {n_{\mathrm{prior}} + n_{w} - 2} \cdot \dfrac{n_{w}}{n_{\mathrm{prior}} + n_{w}}\, G^{-1}. } \] The two fractions cancel to \(S_{\mathrm{marg}} / (n_{\mathrm{prior}} + n_{w} - 2)\) times \(G^{-1}\). So with this \(b_{0}\), the marginal coefficient covariance matches that heuristic pattern---analogous in structure to \(\mathrm{RSS}/(n-p)\) times \((X^{\mathsf T}WX)^{-1}\) when \(S_{\mathrm{marg}} = \mathrm{RSS}\) (e.g.\ \(\mu_{0} = \hat\beta\)) and one reads \(n_{\mathrm{prior}} + n_{w} - 2\) as an effective residual degrees of freedom for the **joint** Normal--Gamma fragment. The factor \(a_{n}-1\) in the denominator of \(b_{n}/(a_{n}-1)\) is exactly half that count, and the factor \((1-pwt) = n_{w}/(n_{\mathrm{prior}}+n_{w})\) from \(V_{n}\) supplies the remaining alignment. **Limit \(n_{\mathrm{prior}} \to 0\) (formal comparison to classical).** Holding \(n_{w} > 2\) and the data fixed, \(pwt = n_{\mathrm{prior}}/(n_{\mathrm{prior}}+n_{w}) \to 0\), \(b_{0} = \frac{1}{2}(n_{\mathrm{prior}}/n_{w})S_{\mathrm{marg}} \to 0\), \(a_{n} = (n_{\mathrm{prior}}+n_{w})/2 \to n_{w}/2\), and \(a_{n}-1 \to (n_{w}-2)/2\). Along the same path, \(S_{\mathrm{marg}} = \mathrm{RSS} + pwt\,(\hat\beta - \mu_{0})^{\mathsf T}G(\hat\beta-\mu_{0}) \to \mathrm{RSS}\) because \(pwt\to 0\) (Section A3.3): the prior--data mismatch vanishes **even when** \(\mu_{0} \neq \hat\beta\). Also \(b_{n} = b_{0} + \frac{1}{2}S_{\mathrm{marg}} \to \frac{1}{2}\mathrm{RSS}\), so \[ \frac{b_{n}}{a_{n}-1} \to \frac{\mathrm{RSS}}{n_{w}-2}, \qquad V_{n} = (1-pwt)\,G^{-1} \to G^{-1}. \] Hence \[ \mathrm{Cov}(\beta \mid y) \to \frac{\mathrm{RSS}}{n_{w}-2}\,(X^{\mathsf T}WX)^{-1}, \] the usual weighted least squares covariance under the Gaussian `glm.fit` dispersion convention of `Prior_Setup()` (see `?Prior_Setup`, Details). **Before** taking \(n_{\mathrm{prior}}\to 0\), a fixed \(n_{\mathrm{prior}}>0\) with \(\mu_{0} \neq \hat\beta\) typically has \(S_{\mathrm{marg}} > \mathrm{RSS}\) and hence a **larger** \(b_{n}\) and posterior spread than this classical limit suggests; only in the limit does \(S_{\mathrm{marg}}\) coincide with \(\mathrm{RSS}\) for every \(\mu_{0}\). At \(n_{\mathrm{prior}} = 0\) the Zellner \(g\)-prior weight and the Gamma shape \(a_{0}\) are degenerate, so this limit is best read as **asymptotic intuition** for weak prior weight, not as a literal prior specification. **Growth of \(n_{w}\) at fixed \(n_{\mathrm{prior}}\).** The posterior mean of \(\tau^{-1}=\sigma^{2}\) in the special case can be written with the uncancelled factors as \[ \frac{b_{n}}{a_{n}-1} = \frac{S_{\mathrm{marg}}\,(n_{\mathrm{prior}} + n_{w})/n_{w}} {n_{\mathrm{prior}} + n_{w} - 2} = \frac{S_{\mathrm{marg}}}{n_{w}} \cdot \frac{n_{\mathrm{prior}} + n_{w}}{n_{\mathrm{prior}} + n_{w} - 2}. \] As \(n_{w} \to \infty\) with \(n_{\mathrm{prior}}\) fixed, \((n_{\mathrm{prior}} + n_{w})/(n_{\mathrm{prior}} + n_{w} - 2) \to 1\) and \((n_{\mathrm{prior}} + n_{w})/n_{w} \to 1\), so asymptotically \((b_{n}/(a_{n}-1)) \sim S_{\mathrm{marg}}/n_{w}\): the leading behavior is **mean marginal residual energy per effective observation**. Moreover \(pwt = n_{\mathrm{prior}}/(n_{\mathrm{prior}}+n_{w}) \to 0\), so in the Zellner decomposition \(S_{\mathrm{marg}} = \mathrm{RSS} + pwt\,(\hat\beta-\mu_{0})^{\mathsf T}G(\hat\beta-\mu_{0})\) the mismatch term is multiplied by a coefficient tending to **zero**; under typical regularity \(pwt\,(\hat\beta-\mu_{0})^{\mathsf T}G(\hat\beta-\mu_{0}) = O_{p}(1)\) while \(\mathrm{RSS}\) grows with \(n_{w}\), so \(S_{\mathrm{marg}}/n_{w}\) and \(\mathrm{RSS}/n_{w}\) have the same leading order. The factor \((n_{\mathrm{prior}}+n_{w})/n_{w}\) is the only explicit **\(O(1/n_{w})\)** correction linking total sample size to effective data count before cancellation with \(V_{n}\). #### `dGamma` posterior with fixed \(\beta\): general form and special choices For a precision-only prior \[ \tau \sim \mathrm{Gamma}(a_0,b_0), \] and Gaussian weighted likelihood with \(\beta\) treated as fixed, define \[ \mathrm{RSS}_w(\beta)=\sum_i w_i(y_i-x_i^\top\beta)^2. \] Then \[ \tau\mid y,\beta \sim \mathrm{Gamma}\!\left(a_0+\frac{n_w}{2},\; b_0+\frac12\,\mathrm{RSS}_w(\beta)\right). \] Using the Chapter 15 calibration \[ a_0=\frac{n_{\mathrm{prior}}}{2}, \qquad b_0=\frac12\frac{n_{\mathrm{prior}}}{n_w}\,\mathrm{RSS}_w(\beta), \] we get \[ \tau\mid y,\beta \sim \mathrm{Gamma}\!\left( \frac{n_{\mathrm{prior}}+n_w}{2}, \frac12\left(1+\frac{n_{\mathrm{prior}}}{n_w}\right)\mathrm{RSS}_w(\beta) \right). \] Hence, for \(n_{\mathrm{prior}}+n_w>2\), \[ E(\tau^{-1}\mid y,\beta) = \frac{b_n}{a_n-1} = \frac{\left(\frac{n_{\mathrm{prior}}+n_w}{n_w}\right)\mathrm{RSS}_w(\beta)} {n_{\mathrm{prior}}+n_w-2}. \] This expression holds for any fixed \(\beta\). #### Choice \(\beta=\beta_\star\) (posterior mode from `dNormal_Gamma`) and weak-prior limit If \(\beta=\beta_\star\), where \(\beta_\star\) is the posterior mode from the `dNormal_Gamma` fit, the same formula applies with \(\mathrm{RSS}_w(\beta_\star)\). Using the weighted least-squares decomposition: \[ \mathrm{RSS}_w(\beta) = \mathrm{RSS}_w(\hat\beta) + (\beta-\hat\beta)^\top G(\beta-\hat\beta),\qquad G=X^\top W X, \] so \[ \mathrm{RSS}_w(\beta_\star) = \mathrm{RSS}_w(\hat\beta) + (\beta_\star-\hat\beta)^\top G(\beta_\star-\hat\beta). \] Under scalar Zellner weighting, \(\beta_\star=(1-pwt)\hat\beta+pwt\,\mu_0\), giving \[ \mathrm{RSS}_w(\beta_\star) = \mathrm{RSS}_w(\hat\beta) + pwt^2(\hat\beta-\mu_0)^\top G(\hat\beta-\mu_0). \] As \(n_{\mathrm{prior}}\to 0\) with \(n_w\) fixed, \(pwt\to 0\), so \[ \mathrm{RSS}_w(\beta_\star)\to \mathrm{RSS}_w(\hat\beta), \] and therefore \[ E(\tau^{-1}\mid y,\beta_\star) \to \frac{\mathrm{RSS}_w(\hat\beta)}{n_w-2}, \] the classical weighted residual-variance limit.