library(glmbayesCore)

Core 20-S04: Row-group priors and samplers

In the two-block Gibbs samplers, Block 1 draws a coefficient vector $b_j$ for every group $j$. Given Block 2, those draws are conditionally independent: each group is a separate Bayesian regression on its own rows, with its own prior. glmbayesCore provides the pieces for that step:

normalize_group(): turns a group specification into a standard row partition;
Prior_SetupGroup(): runs Prior_Setup() (Core 10) separately on each group’s rows;
pfamily_list(): converts those priors into one pfamily per group;
rNormal_reg_group() and rNormalGLM_reg_group(): draw every group’s coefficients in one call. A C++ loop runs the same samplers as rNormal_reg() and rNormalGLM_reg() (Parts 2 and 3) on each group.

These functions also work on their own, for any model that fits a separate regression per group with no pooling. Unlike Prior_Setup_GLMM() (Core 20-S02), Prior_SetupGroup() shares no $\gamma$ or $\tau^2$ across groups.

1. User interface

1.1 Group specifications and `normalize_group()`

normalize_group(group, l2)

Every function in this vignette accepts the grouping (group) in one of three forms, where $l_2$ is the number of rows:

Form	Example	Notes
factor or vector of length $l_2$	`iris$Species`	Factor level order is kept; other vectors are converted with `factor()`
group sizes (length $k < l_2$)	`c(2, 3)`	Rows must already be stored group by group; sizes must sum to $l_2$
list of row-index vectors	`list(a = c(1, 3), b = c(2, 4, 5))`	Must not overlap and must together cover `1:l2`

Prior_SetupGroup() also accepts a column name in data. The return value is a list with the number of groups k, labels ids, sizes l2_blocks, contiguous start positions starts and row indices rows.

1.2 `Prior_SetupGroup()`

Prior_SetupGroup(formula, group, family = gaussian(), data = NULL,
                 weights = NULL, subset = NULL, na.action = na.fail,
                 offset = NULL, contrasts = NULL,
                 pwt = NULL, pwt_default_low = 0.01, pwt_default_high = 0.05,
                 n_prior = NULL, sd = NULL, dispersion = NULL,
                 intercept_source, effects_source, mu = NULL, k = 1, ...)
print(x, blocks = NULL, ...)

Families: gaussian(), poisson(), binomial().
The six calibration arguments pwt, n_prior, sd, mu, dispersion and k take either the usual Prior_Setup() form (applied to every group) or a list named by group with one value per group. A NULL element keeps the default for that group.
Returns a named list of class "Prior_SetupGroup", one Prior_Setup() result per group. Its attributes are call, formula, group and group_info (the normalize_group() result).
print() shows a header. blocks = (names or indices) also prints the selected groups’ setups.

pfamily_list(<Prior_SetupGroup>, ptypes) builds one pfamily per group:

`ptypes`	Constructor call	Requires
`"dNormal"` (default)	`dNormal(mu, Sigma, dispersion)` (dispersion 1 if absent)	any family
`"dNormal_Gamma"`	`dNormal_Gamma(mu, Sigma_0, shape, rate)`	Gaussian
`"dIndependent_Normal_Gamma"`	`dIndependent_Normal_Gamma(mu, Sigma, shape_ING, rate)`	Gaussian

1.3 `rNormal_reg_group()` and `rNormalGLM_reg_group()`

rNormal_reg_group(n, y, x, group, prior_list = NULL, prior_lists = NULL,
                  offset = NULL, weights = 1, Gridtype = 2L)

rNormalGLM_reg_group(n, y, x, group, prior_list = NULL, prior_lists = NULL,
                     offset = NULL, weights = 1, family = gaussian(),
                     Gridtype = 2L, n_envopt = NULL, use_parallel = TRUE,
                     use_opencl = FALSE, verbose = FALSE, progbar = FALSE)

Argument	Role
`n`	Draws per group. Use `n = 1` (one Gibbs step); the returned matrices hold one draw per group
`y`, `x`, `group`	Response, design (same columns in every group) and grouping
`prior_list`	One prior for all groups (`mu`, `P` or `Sigma`, `dispersion`). Group-specific variants: `mu` as an $l_1 \times k$ matrix, `P`/`Sigma` as a list, `dispersion` with one value per group, or a `blocks` sublist
`prior_lists`	Alternatively, a list with one prior per group, in group order (`ids` from `normalize_group()`). Length 1 is recycled
`offset`, `weights`	Length 1 or $l_2$; split by group like `y`
`family`	`rNormalGLM_reg_group()` only: `poisson`, `binomial`, `quasipoisson`, `quasibinomial` or `Gamma` (not `gaussian`)
`Gridtype`, `n_envopt`, `use_parallel`, `use_opencl`, `verbose`	Envelope settings passed to each group’s sampler (Core 06–08)

rNormal_reg_group() needs dispersion in every group’s prior; it holds $\sigma^2$ fixed for the draw. For families without a dispersion parameter, rNormalGLM_reg_group() sets it to 1.

Both return a list (class "rNormal_reg_group" or "rNormalGLM_reg_group"):

Field	Content
`coefficients`	$k \times l_1$ matrix: one draw per group (rows named by group)
`coef.mode`	$k \times l_1$ matrix of posterior modes
`dispersion`	Dispersion used for each group
`group_info`	`normalize_group()` result
`group_results`	Each group’s full sampler output
`prior_lists`	Each group’s prior after normalization (`mu`, `Sigma`, `P`, …)
`n`, `k`, `l1`, `l2`, `y`, `x`, `offset`, `prior.weights`, `call` (and `family`)	Inputs and dimensions

2. Examples

2.1 Group partitions

str(normalize_group(factor(c("a", "a", "b", "b", "b")), l2 = 5L))
#> List of 5
#>  $ k        : int 2
#>  $ ids      : chr [1:2] "a" "b"
#>  $ l2_blocks: Named int [1:2] 2 3
#>   ..- attr(*, "names")= chr [1:2] "a" "b"
#>  $ starts   : Named int [1:2] 1 3
#>   ..- attr(*, "names")= chr [1:2] "" "a"
#>  $ rows     :List of 2
#>   ..$ a: int [1:2] 1 2
#>   ..$ b: int [1:3] 3 4 5
str(normalize_group(c(2L, 3L), l2 = 5L))
#> List of 5
#>  $ k        : int 2
#>  $ ids      : chr [1:2] "group1" "group2"
#>  $ l2_blocks: int [1:2] 2 3
#>  $ starts   : int [1:2] 1 3
#>  $ rows     :List of 2
#>   ..$ : int [1:2] 1 2
#>   ..$ : int [1:3] 3 4 5
str(normalize_group(list(g1 = c(1L, 3L), g2 = c(2L, 4L, 5L)), l2 = 5L))
#> List of 5
#>  $ k        : int 2
#>  $ ids      : chr [1:2] "g1" "g2"
#>  $ l2_blocks: Named int [1:2] 2 3
#>   ..- attr(*, "names")= chr [1:2] "g1" "g2"
#>  $ starts   : Named int [1:2] 1 3
#>   ..- attr(*, "names")= chr [1:2] "" "g1"
#>  $ rows     :List of 2
#>   ..$ g1: int [1:2] 1 3
#>   ..$ g2: int [1:3] 2 4 5

2.2 Independent priors per species

ps_block <- Prior_SetupGroup(
  Sepal.Length ~ Sepal.Width + Petal.Length,
  group  = "Species",
  data   = iris,
  family = gaussian()
)
#> Using default pwt = 0.01 (low-d default).
#> Using default pwt = 0.01 (low-d default).
#> Using default pwt = 0.01 (low-d default).
ps_block
#> 
#> Call:
#> Prior_SetupGroup(formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
#>     group = "Species", family = gaussian(), data = iris)
#> 
#> --- Row-block prior setup ---
#>   formula : Sepal.Length ~ Sepal.Width + Petal.Length
#>   blocks  : 3 (setosa, versicolor, virginica)
#>   Each element is a Prior_Setup() result; build pfamilies with
#>   pfamily_list(). Use print(x, blocks = ...) for details.
names(ps_block)
#> [1] "setosa"     "versicolor" "virginica"
ps_block$setosa$mu
#>                 mu
#> (Intercept)  5.006
#> Sepal.Width  0.000
#> Petal.Length 0.000
ps_block$setosa$dispersion
#> [1] 0.05640586

Use a different prior weight for each species, keeping the default for versicolor:

ps_block2 <- Prior_SetupGroup(
  Sepal.Length ~ Sepal.Width + Petal.Length,
  group  = "Species",
  data   = iris,
  family = gaussian(),
  pwt    = list(setosa = 0.05, versicolor = NULL, virginica = 0.2)
)
#> Using default pwt = 0.01 (low-d default).
sapply(ps_block2, function(p) p$PriorSettings$pwt[1])
#>     setosa versicolor  virginica 
#>       0.05       0.01       0.20

2.3 Converting to `pfamily` objects

pf <- pfamily_list(ps_block)
attr(pf, "ptypes")
#>     setosa versicolor  virginica 
#>  "dNormal"  "dNormal"  "dNormal"
str(pf$setosa$prior_list)
#> List of 4
#>  $ mu        : num [1:3, 1] 5.01 0 0
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:3] "(Intercept)" "Sepal.Width" "Petal.Length"
#>   .. ..$ : chr "mu"
#>  $ Sigma     : num [1:3, 1:3] 14.892 -2.343 -4.616 -2.343 0.819 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:3] "(Intercept)" "Sepal.Width" "Petal.Length"
#>   .. ..$ : chr [1:3] "(Intercept)" "Sepal.Width" "Petal.Length"
#>  $ dispersion: num 0.0564
#>  $ ddef      : logi FALSE
#>  - attr(*, "Prior Type")= chr "dNormal"

pf_ng <- pfamily_list(ps_block, ptypes = "dNormal_Gamma")
attr(pf_ng, "ptypes")
#>          setosa      versicolor       virginica 
#> "dNormal_Gamma" "dNormal_Gamma" "dNormal_Gamma"

2.4 Gaussian group draws with one shared prior

set.seed(42)
n_groups <- 3L
n_per    <- 10L
grp      <- rep(seq_len(n_groups), each = n_per)
x        <- cbind(1, rnorm(n_groups * n_per))
colnames(x) <- c("(Intercept)", "X1")
b_true   <- matrix(c(5, 0.5, 3, -0.2, 7, 0.3), nrow = n_groups, byrow = TRUE)
sigma2   <- 1.5
y        <- rowSums(x * b_true[grp, ]) + rnorm(nrow(x), sd = sqrt(sigma2))

prior_list <- list(
  mu         = rep(0, ncol(x)),
  Sigma      = diag(100, ncol(x)),
  dispersion = sigma2
)

out <- rNormal_reg_group(n = 1L, y = y, x = x, group = grp,
                         prior_list = prior_list)
out$coefficients
#>   (Intercept)          X1
#> 1    4.712982 -0.06290662
#> 2    3.230819  0.01360927
#> 3    6.811915  0.75763105
out$coef.mode
#>   (Intercept)         X1
#> 1    4.904492 -0.1532808
#> 2    2.948469 -0.3367292
#> 3    7.011645  0.2979490

Each row is one group. Compare coef.mode with b_true.

2.5 Gaussian group draws with per-group priors from `Prior_SetupGroup()`

The prior_list inside each dNormal element of pfamily_list(ps_block) has the form prior_lists expects. Because prior_lists is matched by position, order it by the group levels:

X_iris <- model.matrix(~ Sepal.Width + Petal.Length, data = iris)
prior_lists <- lapply(pf[levels(iris$Species)], function(p) p$prior_list)

out_iris <- rNormal_reg_group(
  n           = 1L,
  y           = iris$Sepal.Length,
  x           = X_iris,
  group       = iris$Species,
  prior_lists = prior_lists
)
out_iris$coef.mode
#>            (Intercept) Sepal.Width Petal.Length
#> setosa       2.3307608   0.6607420    0.2805851
#> versicolor   2.1546271   0.2451657    0.7282309
#> virginica    0.6844146   0.2573544    0.9254707
out_iris$dispersion
#> [1] 0.05640586 0.11713771 0.10376182

2.6 Non-Gaussian group draws

set.seed(1)
grp2 <- factor(rep(c("A", "B"), each = 15L))
x2 <- cbind(1, rnorm(length(grp2)))
colnames(x2) <- c("(Intercept)", "x1")
eta <- x2 %*% c(-0.5, 0.3) + ifelse(grp2 == "A", 0, 0.4)
y2 <- rpois(length(grp2), exp(eta))

out_glm <- rNormalGLM_reg_group(
  n            = 1L,
  y            = y2,
  x            = x2,
  group        = grp2,
  prior_list   = list(mu = rep(0, 2), P = diag(10, 2), dispersion = 1),
  family       = poisson(),
  use_parallel = FALSE
)
out_glm$coefficients
#>   (Intercept)        x1
#> A  -0.5460683 0.1416992
#> B   0.2841374 0.3334142
out_glm$coef.mode
#>   (Intercept)           x1
#> A -0.35332535 -0.001877103
#> B  0.04112814  0.349076535

Here the prior is given by its precision P. Sigma works too.

3. Lineage

3.1 `Prior_SetupGroup()` and its `pfamily_list()` method

Prior_SetupGroup(formula, group, family, data, ..., pwt, n_prior, sd, mu,
                 dispersion, k)
 ├─ family check (gaussian / poisson / binomial)
 ├─ .prior_setup_group_formula_block_meta()        model frame + partition
 │   ├─ stats::model.frame()
 │   ├─ .prior_setup_group_resolve_block()          column name / vector / list → group
 │   ├─ normalize_group()                           exported
 │   └─ stats::model.matrix()                       predictor names
 ├─ .prior_setup_group_resolve_calibration_arg()    x6: shared value or named per-group list
 └─ for each group b:
     ├─ .prior_setup_group_rows_to_data_subset()     model-frame rows → data rows
     └─ Prior_Setup(formula, subset = rows_b, ...)   exported (Core 10)

print.Prior_SetupGroup(x, blocks)
 ├─ .prior_setup_group_select_named_list_keys()
 └─ print.PriorSetup() for each selected group

pfamily_list.Prior_SetupGroup(object, ptypes)
 └─ for each group:
     ├─ "dNormal"                    → dNormal(mu, (Sigma + t(Sigma))/2, dispersion)
     ├─ "dNormal_Gamma"              → dNormal_Gamma(mu, Sigma_0, shape, rate)
     └─ "dIndependent_Normal_Gamma"  → dIndependent_Normal_Gamma(mu, Sigma, shape_ING, rate)

3.2 The groupwise samplers (R to C++)

rNormal_reg_group()                    rNormalGLM_reg_group()
  R/simfunction_group.R                  R/simfunction_group.R
  ├─ recycle offset / weights            ├─ family + link checks, recycle offset / weights
  ├─ glmbfamfunc(gaussian())             ├─ glmbfamfunc(family)       f2 / f3 callbacks
  └─ .group_rNormalReg_cpp()             └─ .group_rNormalGLM_cpp()   R/rcpp_wrappers.R
        │  .Call                                │  .Call
        ▼                                       ▼
  group_rNormalReg_cpp_export()          group_rNormalGLM_cpp_export()     src/group_utils.cpp
  ├─ normalize_group_cpp()               ├─ normalize_group_cpp()
  │   └─ split_factor_rows()             │   └─ split_factor_rows()
  ├─ normalize_prior_for_groups_cpp()    ├─ normalize_prior_for_groups_cpp()
  │   ├─ base_prior_block()              │   ├─ base_prior_block()
  │   │   └─ prior_list_to_P_Sigma()     │   │   └─ prior_list_to_P_Sigma()
  │   └─ rep_list_blocks()               │   └─ rep_list_blocks()
  ├─ require dispersion in each group    │
  ├─ prior_payload_from_groups()         ├─ prior_payload_from_groups()
  └─ rNormalRegGroups()                  └─ rNormalGLMGroups()
       src/rNormalRegGroups.cpp               src/rNormalGLMGroups.cpp
       for each group b:                      for each group b:
       ├─ slice y, x, offset, wt              ├─ slice y, x, offset, wt
       ├─ mu / P / dispersion for b           ├─ mu / P / dispersion for b
       └─ glmbayes::sim::rNormalReg()         └─ glmbayes::sim::rNormalGLM()
            (Gaussian conjugate draw)               (envelope accept-reject; Parts 2-3)

3.3 Helper reference: R

Helper	File	Role
`.prior_setup_group_formula_block_meta()`	`Prior_SetupGroup.R`	Builds the model frame once, resolves `group` against it and partitions the rows
`.prior_setup_group_resolve_block()`	`Prior_SetupGroup.R`	Converts a column name, vector or list into a group vector aligned with model-frame rows
`.prior_setup_group_resolve_calibration_arg()`	`Prior_SetupGroup.R`	Turns one calibration argument into a per-group list; checks names exactly match the group IDs
`.prior_setup_group_rows_to_data_subset()`	`Prior_SetupGroup.R`	Maps model-frame row positions back to `data` rows for `Prior_Setup(subset = )`
`.prior_setup_group_select_named_list_keys()`	`Prior_SetupGroup.R`	Resolves `blocks =` for `print.Prior_SetupGroup()`
`.group_rNormalReg_cpp()`, `.group_rNormalGLM_cpp()`	`rcpp_wrappers.R`	`.Call` wrappers for the C++ entry points
`.rNormalRegGroups_cpp()`, `.rNormalGLMGroups_cpp()`	`rcpp_wrappers.R`	Lower-level wrappers that call the per-group loops directly with an already normalized partition and prior; not called by the exported functions in this vignette
`glmbfamfunc()`	`simulationpipeline.R`	Exported; returns the family’s log-posterior (`f2`) and gradient (`f3`) callbacks

3.4 Helper reference: C++ (`glmbayes::sim::group`)

Function	File	Role
`group_rNormalReg_cpp_export()`, `group_rNormalGLM_cpp_export()`	`group_utils.cpp`	Entry points: normalize the partition and priors, then call the per-group loop
`normalize_group_cpp()`, `split_factor_rows()`	`group_utils.cpp`	C++ equivalent of `normalize_group()`
`normalize_prior_for_groups_cpp()`	`group_utils.cpp`	Expands `prior_list` / `prior_lists` / `blocks` / matrix `mu` into one prior per group
`prior_list_to_P_Sigma()`, `base_prior_block()`	`group_utils.cpp`	Check symmetry and positive definiteness; fill in whichever of `P` and `Sigma` is missing; check dimensions
`rep_list_blocks()`	`group_utils.cpp`	Copies one prior to all groups
`prior_payload_from_groups()`	`group_utils.cpp`	Packs per-group priors into the `mu` matrix / `P_blocks` / `dispersion` layout; detects when all groups share one prior
`rNormalRegGroups()`, `rNormalGLMGroups()`	`rNormalRegGroups.cpp`, `rNormalGLMGroups.cpp`	Per-group loops: slice data, select that group’s prior, call the single-model sampler, and collect draws and modes
`slice_`, `_for_block`, `validate_row_blocks`, `coef_row_from_block*`	same files	Small slicing and validation utilities used by the loops

4. Role in the mixed-model samplers

In lmerb() and glmerb(), each Gibbs iteration draws Block 1 with these row-group samplers. Each group’s prior mean comes from the current Block 2 draw ($W_{k,j}^\top \gamma_k$ for each coefficient $k$), and its prior covariance from the current $\tau^2_k$. The Gibbs engines that build these per-sweep priors are in lmebayesCore. The per-group draw itself is the glmbayesCore code shown here.

- Core 20-S04: Row-group priors and samplers

Form	Example	Notes
factor or vector of length \(l_2\)	`iris$Species`	Factor level order is kept; other vectors are converted with `factor()`
group sizes (length \(k < l_2\))	`c(2, 3)`	Rows must already be stored group by group; sizes must sum to \(l_2\)
list of row-index vectors	`list(a = c(1, 3), b = c(2, 4, 5))`	Must not overlap and must together cover `1:l2`

Field	Content
`coefficients`	\(k \times l_1\) matrix: one draw per group (rows named by group)
`coef.mode`	\(k \times l_1\) matrix of posterior modes
`dispersion`	Dispersion used for each group
`group_info`	`normalize_group()` result
`group_results`	Each group’s full sampler output
`prior_lists`	Each group’s prior after normalization (`mu`, `Sigma`, `P`, …)
`n`, `k`, `l1`, `l2`, `y`, `x`, `offset`, `prior.weights`, `call` (and `family`)	Inputs and dimensions

Core 20-S04: Row-group priors and samplers

Core 20-S04: Row-group priors and samplers

1. User interface

1.1 Group specifications and normalize_group()

1.2 Prior_SetupGroup()

1.3 rNormal_reg_group() and rNormalGLM_reg_group()

2. Examples

2.1 Group partitions

2.2 Independent priors per species

2.3 Converting to pfamily objects

2.4 Gaussian group draws with one shared prior

2.5 Gaussian group draws with per-group priors from Prior_SetupGroup()

2.6 Non-Gaussian group draws

3. Lineage

3.1 Prior_SetupGroup() and its pfamily_list() method

3.2 The groupwise samplers (R to C++)

3.3 Helper reference: R

3.4 Helper reference: C++ (glmbayes::sim::group)

4. Role in the mixed-model samplers

1.1 Group specifications and `normalize_group()`

1.2 `Prior_SetupGroup()`

1.3 `rNormal_reg_group()` and `rNormalGLM_reg_group()`

2.3 Converting to `pfamily` objects

2.5 Gaussian group draws with per-group priors from `Prior_SetupGroup()`

3.1 `Prior_SetupGroup()` and its `pfamily_list()` method

3.4 Helper reference: C++ (`glmbayes::sim::group`)