demo/Ex_10_Prior_Setup_gaussian_calibration.R

In glmbayes: Bayesian Generalized Linear Models (IID Samples)

## Prior_Setup Gaussian calibration: dNormal, dNormal_Gamma, dIndependent_Normal_Gamma
##
## `Prior_Setup()` uses a single Gamma(shape, rate) calibration from `n_prior`
## (shape = (n_prior+k)/2, default k=1; see ?compute_gaussian_prior). `dNormal()` uses
## `ps$Sigma` and `dispersion = ps$dispersion`. The same returned `rate`, `mu`,
## `Sigma`, and `Sigma_0` apply to both `dNormal_Gamma()` and
## `dIndependent_Normal_Gamma()`; for the latter, pass `shape = ps$shape_ING`
## (see ?Prior_Setup).
##
## Run: demo(Ex_10_Prior_Setup_gaussian_calibration, package = "glmbayes")

library(glmbayes)

vcov_congruence <- function(V_post, V_base) {
  U <- chol(V_base)
  Ui <- backsolve(U, diag(nrow(U)))
  M <- crossprod(Ui, V_post %*% Ui)
  ev <- eigen(M, symmetric = TRUE, only.values = TRUE)$values
  list(M = M, eigenvalues = ev)
}

congruence_line <- function(V_post, label, V_base, base_label) {
  cq <- vcov_congruence(V_post, V_base)
  ev <- cq$eigenvalues
  sprintf(
    "%s vs %s: eigenvalues min / mean / max = %s / %s / %s",
    label,
    base_label,
    format(min(ev), digits = 6),
    format(mean(ev), digits = 6),
    format(max(ev), digits = 6)
  )
}

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)

n_mc <- 100000L
pwt <- 0.001

fit_lm <- lm(weight ~ group, x = TRUE, y = TRUE)
V_lm <- vcov(fit_lm)
V_lm_shrunk <- (1 - pwt) * V_lm

ps <- Prior_Setup(
  weight ~ group,
  pwt = pwt,
  intercept_source = "full_model",
  effects_source = "full_model"
)

fit_ng <- lmb(
  weight ~ group,
  pfamily = dNormal_Gamma(
    ps$mu,
    Sigma_0 = ps$Sigma_0,
    shape = ps$shape,
    rate = ps$rate
  ),
  n = n_mc
)

fit_ing <- lmb(
  weight ~ group,
  pfamily = dIndependent_Normal_Gamma(
    ps$mu,
    ps$Sigma,
    shape = ps$shape_ING,
    rate = ps$rate
  ),
  n = n_mc
)

fit_dn <- lmb(
  weight ~ group,
  pfamily = dNormal(
    mu = ps$mu,
    Sigma = ps$Sigma,
    dispersion = ps$dispersion
  ),
  n = n_mc
)

cat("\n======== One Prior_Setup(); NG uses ps$shape; ING uses ps$shape_ING =========\n")
cat("shape (NG) =", ps$shape, "  shape_ING =", ps$shape_ING, "  rate =", ps$rate, "\n")

cat("\n======== Classical vs. Posterior vcov(lmb) =========\n")
cat("Classical Scaled:\n")
print((1 - pwt) * vcov(fit_lm))
cat("\ndNormal (dispersion = ps$dispersion):\n")
print(vcov(fit_dn))
cat("\ndNormal_Gamma:\n")
print(vcov(fit_ng))
cat("\ndIndependent_Normal_Gamma:\n")
print(vcov(fit_ing))

cat("\n======== Congruence vs (1-pwt)*vcov(lm) =========\n")
cat("pwt =", pwt, "\n")
cat(congruence_line(vcov(fit_dn), "dNormal", V_lm_shrunk, "(1-pwt)*vcov(lm)"), "\n", sep = "")
cat(congruence_line(vcov(fit_ng), "NG", V_lm_shrunk, "(1-pwt)*vcov(lm)"), "\n", sep = "")
cat(congruence_line(vcov(fit_ing), "ING", V_lm_shrunk, "(1-pwt)*vcov(lm)"), "\n", sep = "")

cat("\n======== Dispersion (classical vs posterior means) =========\n")
disp_classical <- summary(fit_lm)$sigma^2
cat("lm sigma^2 (RSS/(n-p))   =", disp_classical, "\n")
cat("ps$dispersion (dNormal)  =", ps$dispersion, "\n")
cat("mean(lmb dNormal dispersion) =", mean(fit_dn$dispersion), "\n")
cat("mean(lmb NG dispersion)  =", mean(fit_ng$dispersion), "\n")
cat("mean(lmb ING dispersion) =", mean(fit_ing$dispersion), "\n")

invisible(NULL)