library(glmbayesCore)

Core 20-S03: From calibrated priors to `pfamily` objects

Prior_Setup_GLMM() (Core 20-S02) returns calibrated prior settings. Samplers instead take pfamily objects, the same prior-family objects rglmb() and rlmb() use (see ?pfamily). Each one bundles a prior_list with the simulation function that knows how to use it. Two functions do the conversion:

pfamily_list(): one pfamily per random coefficient, for the Block 2 regressions of $b_{\cdot k}$ on $W_k$. This is the pfamily_list argument of lmerb() and glmerb().
dGamma_list(): one dGamma() per group, for the Block 1 measurement variance when each group has its own $\sigma^2_j$. This is the dispersion_ranef argument of lmerb().

1. User interface

1.1 `pfamily_list()`

pfamily_list(object, ptypes = NULL, ...)
print(x, components = NULL, ...)

pfamily_list() is an S3 generic with two methods:

Method	One `pfamily` per	Allowed `ptypes`	Documented in
`pfamily_list.Prior_Setup_GLMM()`	random coefficient $k$	`"dNormal"` (default), `"dIndependent_Normal_Gamma"`	this vignette
`pfamily_list.Prior_SetupGroup()`	row group	`"dNormal"` (default), `"dNormal_Gamma"`, `"dIndependent_Normal_Gamma"`	Core 20-S04

ptypes can be a single string applied to every element, or one string per element (named or in order). The result has class c("pfamily_list", "list") and an attribute "ptypes" recording the resolved choice. print() shows all elements by default. Use components = (names or indices) to show a subset.

For a Prior_Setup_GLMM object, component $k$ is built from pop.prior_list[[k]] and pop.ing_prior[[k]]:

`ptypes[k]`	Constructor call	Effect in the sampler
`"dNormal"`	`dNormal(mu, Sigma, dispersion = tau2_k)`	$\tau^2_k$ fixed at the calibrated centre
`"dIndependent_Normal_Gamma"`	`dIndependent_Normal_Gamma(mu, Sigma, shape, rate, disp_lower, disp_upper)`	$\tau^2_k$ also sampled, with a truncated Gamma prior on $1/\tau^2_k$

The independent Normal–Gamma choice requires $n_0 \le J$ (pop.dispersion.nprior at most the number of groups, or equivalently pop.dispersion.pwt at most 0.5).

1.2 `dGamma_list()`

dGamma_list(object, ...)
## S3 method for class 'Prior_Setup_GLMM'
dGamma_list(object, max_disp_perc_measurement = NULL, ...)
print(x, groups = NULL, ...)

Requires a Gaussian Prior_Setup_GLMM() object fitted with dispformula = ~<group>. Otherwise it stops with a message saying what to change. Group $j$ gets

dGamma(shape = shape_ING_j, rate = rate_j, beta = <placeholder>,
       Inv_Dispersion = TRUE, max_disp_perc = mdp_j,
       disp_lower = lower_j, disp_upper = upper_j)

with values from group.ing_prior[[j]]. max_disp_perc_measurement = NULL reuses each group’s stored window. A scalar or a vector with one value per group recomputes the windows without re-running Prior_Setup_GLMM(). The result has class c("dGamma_list", "list") and three attributes:

Attribute	Content
`window_diagnostics`	Data frame, one row per group: `n_j`, `sigma2_hat`, `shape_ING`, `rate`, `max_disp_perc`, `disp_lower`, `disp_upper`
`measurement_prior_group`	The same `shape`/`rate`/`disp_lower`/`disp_upper` values as named vectors, in the form the per-group Gaussian engines take
`group.dispersion.fit`, `calibration_source`	Carried over from the setup so `lmerb()` can reuse the `glmmTMB` fit

2. Example

Same data and formula as Core 20-S02, with one pooled-dispersion setup and one per-group setup.

data(big_word_club, package = "bayesrules")
dat <- big_word_club
dat$school_id <- factor(dat$school_id)
dat <- subset(
  dat,
  !is.na(score_ppvt) &
    !is.na(invalid_ppvt) & invalid_ppvt == 0L &
    complete.cases(dat[, c(
      "score_ppvt", "distracted_a1", "distracted_ppvt",
      "private_school", "title1", "free_reduced_lunch", "school_id"
    )])
)

form_lmer <- score_ppvt ~
  private_school + title1 + free_reduced_lunch +
  distracted_ppvt + distracted_a1 +
  free_reduced_lunch:distracted_a1 +
  (1 + distracted_ppvt + distracted_a1 || school_id)

design0 <- model_setup(form_lmer, data = dat, fit_mer = FALSE)
dat <- subset(dat, school_id %in% names(design0$groupef.rank)[design0$groupef.rank])
dat$school_id <- droplevels(dat$school_id)

ps_pooled <- Prior_Setup_GLMM(form_lmer, data = dat, pop.pwt = 0.01)

ps_group <- Prior_Setup_GLMM(
  form_lmer,
  data                 = dat,
  pop.pwt              = 0.01,
  dispformula          = ~school_id,
  group.max_disp_perc  = 0.8,
  group.dispersion.pwt = 0.1,
  group.alpha_target   = NULL
)

2.1 Normal priors with fixed $\tau^2_k$ (default)

pf <- pfamily_list(ps_pooled)
class(pf)
#> [1] "pfamily_list" "list"
attr(pf, "ptypes")
#>     (Intercept) distracted_ppvt   distracted_a1 
#>       "dNormal"       "dNormal"       "dNormal"
print(pf, components = "distracted_a1")
#> 
#> [[distracted_a1]]
#> 
#> Call:  dNormal(mu = mu_k, Sigma = Sig_k, dispersion = d_k)
#> 
#> Prior Family: dNormal 
#> 
#> Prior List:
#> 
#> $mu
#>                    [,1]
#> (Intercept)           0
#> free_reduced_lunch    0
#> 
#> $Sigma
#>                    (Intercept) free_reduced_lunch
#> (Intercept)         1715.54041        -20.7836363
#> free_reduced_lunch   -20.78364          0.3411801
#> 
#> $dispersion
#> [1] 11.47851
#> 
#> $ddef
#> [1] FALSE
#> 
#> attr(,"Prior Type")
#> [1] "dNormal"

Each element is an ordinary pfamily. Its prior_list holds the values the sampler reads:

str(pf[["distracted_a1"]]$prior_list)
#> List of 4
#>  $ mu        : num [1:2, 1] 0 0
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:2] "(Intercept)" "free_reduced_lunch"
#>   .. ..$ : NULL
#>  $ Sigma     : num [1:2, 1:2] 1715.54 -20.784 -20.784 0.341
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:2] "(Intercept)" "free_reduced_lunch"
#>   .. ..$ : chr [1:2] "(Intercept)" "free_reduced_lunch"
#>  $ dispersion: num 11.5
#>  $ ddef      : logi FALSE
#>  - attr(*, "Prior Type")= chr "dNormal"

2.2 Independent Normal–Gamma priors (sample $\tau^2_k$)

pf_ing <- pfamily_list(ps_pooled, ptypes = "dIndependent_Normal_Gamma")
attr(pf_ing, "ptypes")
#>                 (Intercept)             distracted_ppvt 
#> "dIndependent_Normal_Gamma" "dIndependent_Normal_Gamma" 
#>               distracted_a1 
#> "dIndependent_Normal_Gamma"
str(pf_ing[["(Intercept)"]]$prior_list[c("shape", "rate", "disp_lower", "disp_upper")])
#> List of 4
#>  $ shape     : num 2.67
#>  $ rate      : num 266
#>  $ disp_lower: num 91.6
#>  $ disp_upper: num 284

These are exactly ps_pooled$pop.ing_prior[["(Intercept)"]].

2.3 Mixing prior types

Sample $\tau^2$ for the random intercept only:

re <- names(ps_pooled$pop.prior_list)
ptypes_mixed <- stats::setNames(
  c("dIndependent_Normal_Gamma", rep("dNormal", length(re) - 1L)),
  re
)
attr(pfamily_list(ps_pooled, ptypes = ptypes_mixed), "ptypes")
#>                 (Intercept)             distracted_ppvt 
#> "dIndependent_Normal_Gamma"                   "dNormal" 
#>               distracted_a1 
#>                   "dNormal"

2.4 Per-group measurement-variance priors

A pooled setup has no per-group calibration:

try(dGamma_list(ps_pooled))
#> Error : object's group.ing_prior is not per-group; call Prior_Setup_GLMM(..., dispformula = ~school_id) on a Gaussian model to calibrate per-group measurement-dispersion priors (dispformula = ~1, the default, skips this calibration).

With the per-group setup:

disp_pf <- dGamma_list(ps_group)
length(disp_pf)
#> [1] 34
print(disp_pf, groups = 1:2)
#> 
#> [[2]]
#> 
#> Call:  dGamma(shape = g$shape_ING, rate = g$rate, beta = matrix(0, 1, 
#>     1, dimnames = list("(Intercept)", NULL)), Inv_Dispersion = TRUE, 
#>     max_disp_perc = mdp_lev, disp_lower = disp_lower, disp_upper = disp_upper)
#> 
#> Prior Family: dGamma 
#> 
#> Prior List:
#> 
#> $shape
#> [1] 2.666667
#> 
#> $rate
#> [1] 1703.799
#> 
#> $beta
#>             [,1]
#> (Intercept)    0
#> 
#> $Inv_Dispersion
#> [1] TRUE
#> 
#> $max_disp_perc
#> [1] 0.8
#> 
#> $disp_lower
#> [1] 712.6046
#> 
#> $disp_upper
#> [1] 1275.768
#> 
#> attr(,"Prior Type")
#> [1] "dGamma"
#> 
#> [[3]]
#> 
#> Call:  dGamma(shape = g$shape_ING, rate = g$rate, beta = matrix(0, 1, 
#>     1, dimnames = list("(Intercept)", NULL)), Inv_Dispersion = TRUE, 
#>     max_disp_perc = mdp_lev, disp_lower = disp_lower, disp_upper = disp_upper)
#> 
#> Prior Family: dGamma 
#> 
#> Prior List:
#> 
#> $shape
#> [1] 2.722222
#> 
#> $rate
#> [1] 1503.164
#> 
#> $beta
#>             [,1]
#> (Intercept)    0
#> 
#> $Inv_Dispersion
#> [1] TRUE
#> 
#> $max_disp_perc
#> [1] 0.8
#> 
#> $disp_lower
#> [1] 616.871
#> 
#> $disp_upper
#> [1] 1084.198
#> 
#> attr(,"Prior Type")
#> [1] "dGamma"
head(attr(disp_pf, "window_diagnostics"), 3L)
#>   group n_j sigma2_hat shape_ING      rate max_disp_perc disp_lower disp_upper
#> 1     2  12  1022.2794  2.666667 1703.7989           0.8   712.6046  1275.7680
#> 2     3  13   872.8049  2.722222 1503.1639           0.8   616.8710  1084.1978
#> 3     5  11   447.6161  2.611111  721.1592           0.8   307.2302   561.3532

2.5 Changing the truncation windows

disp_pf_wide <- dGamma_list(ps_group, max_disp_perc_measurement = 0.95)
cbind(
  narrow = head(attr(disp_pf,      "window_diagnostics")$disp_upper, 3L),
  wide   = head(attr(disp_pf_wide, "window_diagnostics")$disp_upper, 3L)
)
#>         narrow      wide
#> [1,] 1275.7680 1759.1728
#> [2,] 1084.1978 1477.7743
#> [3,]  561.3532  784.1306

A wider window covers more of each $\sigma^2_j$ posterior but usually lowers the sampler’s acceptance rate. The recomputed window uses the posterior shape $\text{shape}_{ING} + n_j/2$ with a rate matched to sigma2_hat. The prior handed to the sampler (shape_ING, rate) does not change.

2.6 Pooled measurement dispersion

With dispformula = ~1 there is a single $\sigma^2$. Its Gamma calibration (shape, rate, disp_lower, disp_upper) is in ps_pooled$group.ing_prior, so no list is needed.

2.7 Putting it together

fit <- lmebayes::lmerb(
  form_lmer, data = dat,
  pfamily_list     = pfamily_list(ps_group),
  dispersion_ranef = dGamma_list(ps_group),
  dispformula      = ~school_id
)

dispformula must be the same in Prior_Setup_GLMM() and lmerb(). The two arguments are independent, so the caller has to keep them consistent.

3. Lineage

3.1 Call trees

pfamily_list(object, ptypes)                    generic (R/pfamily_list.R)
 └─ pfamily_list.Prior_Setup_GLMM()             R/pfamily_list_Prior_Setup_GLMM.R
     ├─ validate / recycle / reorder ptypes     against names(pop.prior_list)
     ├─ n_prior_for() (local)                   pop.dispersion.nprior[k], or J w/(1-w)
     └─ for each RE coefficient k:
         ├─ "dNormal"
         │    └─ dNormal(mu, Sigma, dispersion)                      exported
         └─ "dIndependent_Normal_Gamma"
              ├─ stop if n_prior_k > J
              ├─ pop.ing_prior[[k]]   (or, if absent, rebuild with
              │                        .lmebayes_ing_limiting_posterior_window())
              └─ dIndependent_Normal_Gamma(mu, Sigma, shape, rate,
                                           disp_lower, disp_upper)   exported

print.pfamily_list(x, components)
 ├─ .pfamily_list_select_keys()                 names / indices → keys
 └─ print.pfamily() for each selected element

dGamma_list(object, ...)                        generic (R/dGamma_list.R)
 └─ dGamma_list.Prior_Setup_GLMM()              R/dGamma_list_Prior_Setup_GLMM.R
     ├─ stop unless gaussian and
     │   .lmebayes_ing_prior_is_grouped(group.ing_prior)
     ├─ .lmebayes_expand_scalar_or_vector()     max_disp_perc_measurement → per group
     ├─ for each group j:
     │   ├─ [window override] .lmebayes_ing_prior_quantile_window(
     │   │        shape_ING + n_j/2, sigma2_hat * (shape_ING + n_j/2 - 1), mdp_j)
     │   └─ dGamma(shape_ING, rate, Inv_Dispersion = TRUE, ...)      exported
     └─ attach window_diagnostics, measurement_prior_group,
        group.dispersion.fit, calibration_source

print.dGamma_list(x, groups)
 ├─ .lmebayes_select_named_list_keys()
 └─ print.pfamily() for each selected group

3.2 Helper reference

Helper	File	Role
`pfamily_list.Prior_Setup_GLMM()`	`pfamily_list_Prior_Setup_GLMM.R`	Block 2 method; one `pfamily` per random coefficient
`pfamily_list.Prior_SetupGroup()`	`Prior_SetupGroup.R`	Row-group method (Core 20-S04)
`.pfamily_list_select_keys()`	`pfamily_list.R`	Resolves `components =` for `print.pfamily_list()`
`dGamma_list.Prior_Setup_GLMM()`	`dGamma_list_Prior_Setup_GLMM.R`	Block 1 method; one `dGamma()` per group
`.lmebayes_select_named_list_keys()`	`dGamma_list.R`	Resolves `groups =` for `print.dGamma_list()`
`.lmebayes_ing_prior_is_grouped()`	`Prior_Setup_GLMM.R`	Tells the per-group and pooled `group.ing_prior` shapes apart
`.lmebayes_expand_scalar_or_vector()`	`mer_prior_helpers.R`	Recycles or reorders `max_disp_perc_measurement` to one value per group
`.lmebayes_ing_prior_quantile_window()`	`mer_prior_helpers.R`	Gamma-quantile window when a group’s tail probability changes
`.lmebayes_ing_limiting_posterior_window()`	`mer_prior_helpers.R`	Fallback $\tau^2_k$ window for setups without `pop.ing_prior`
`dNormal()`, `dIndependent_Normal_Gamma()`, `dGamma()`	`pfamily.R`	Exported constructors; see `?pfamily` and Core 02

3.3 Where the values come from

`pfamily` field	Source in `Prior_Setup_GLMM()` (stage numbers from Core 20-S02)
Block 2 `mu`, `Sigma`	`pop.prior_list[[k]]` (Stage 7)
Block 2 `dispersion` (dNormal)	`pop.prior_list[[k]]$dispersion` (Stage 8)
Block 2 `shape`, `rate`, `disp_lower`, `disp_upper`	`pop.ing_prior[[k]]` (Stage 9)
Block 1 `shape`, `rate`, window	`group.ing_prior[[j]]` (Stages 10–11)

4. Next

Core 20-S04 covers the row-group samplers that draw the Block 1 group coefficients given Block 2.

- Core 20-S03: From calibrated priors to pfamily objects

Core 20-S03: From calibrated priors to pfamily objects — pfamily_list() and dGamma_list()

Core 20-S03: From calibrated priors to `pfamily` objects

1. User interface

1.1 `pfamily_list()`

1.2 `dGamma_list()`

2. Example

2.1 Normal priors with fixed \(\tau^2_k\) (default)

2.2 Independent Normal–Gamma priors (sample \(\tau^2_k\))

2.3 Mixing prior types

2.4 Per-group measurement-variance priors

2.5 Changing the truncation windows

2.6 Pooled measurement dispersion

2.7 Putting it together

3. Lineage

3.1 Call trees

3.2 Helper reference

3.3 Where the values come from

4. Next

Method	One `pfamily` per	Allowed `ptypes`	Documented in
`pfamily_list.Prior_Setup_GLMM()`	random coefficient \(k\)	`"dNormal"` (default), `"dIndependent_Normal_Gamma"`	this vignette
`pfamily_list.Prior_SetupGroup()`	row group	`"dNormal"` (default), `"dNormal_Gamma"`, `"dIndependent_Normal_Gamma"`	Core 20-S04

`ptypes[k]`	Constructor call	Effect in the sampler
`"dNormal"`	`dNormal(mu, Sigma, dispersion = tau2_k)`	\(\tau^2_k\) fixed at the calibrated centre
`"dIndependent_Normal_Gamma"`	`dIndependent_Normal_Gamma(mu, Sigma, shape, rate, disp_lower, disp_upper)`	\(\tau^2_k\) also sampled, with a truncated Gamma prior on \(1/\tau^2_k\)

Core 20-S03: From calibrated priors to pfamily objects — pfamily_list() and dGamma_list()

Core 20-S03: From calibrated priors to pfamily objects

1. User interface

1.1 pfamily_list()

1.2 dGamma_list()

2. Example

2.1 Normal priors with fixed \(\tau^2_k\) (default)

2.2 Independent Normal–Gamma priors (sample \(\tau^2_k\))

2.3 Mixing prior types

2.4 Per-group measurement-variance priors

2.5 Changing the truncation windows

2.6 Pooled measurement dispersion

2.7 Putting it together

3. Lineage

3.1 Call trees

3.2 Helper reference

3.3 Where the values come from

4. Next

Core 20-S03: From calibrated priors to `pfamily` objects

1.1 `pfamily_list()`

1.2 `dGamma_list()`