1.1 Formula rules

model_setup() accepts the formula language of lme4::lmer(), with three structural rules:

1. Exactly one grouping factor, for example school_id.
2. || for more than one random coefficient. The sampler treats the random-effect covariance as diagonal, so write (1 + x || school_id) rather than (1 + x | school_id). A random intercept alone must be written (1 | school_id); lme4 does not accept (1 || school_id).
3. Each fixed effect must be one of three types:
   • a group-level (level-2) covariate: constant within each group, such as private_school;
   • a population mean slope: an observation-level covariate that also has a random slope, such as distracted_ppvt;
   • a cross-level interaction between a level-2 covariate and a random slope, such as free_reduced_lunch:distracted_a1.

Any other fixed term (for example an observation-level covariate with no matching random slope) is rejected with an error naming the offending columns.

These rules let model_setup() split the lme4 fixed effects into one group-level regression per random coefficient:

1.2 Arguments

1.3 Returned fields

1.4 check_identifiability()

check_identifiability(y, D, group, W, family = gaussian(), group_name, verbose = FALSE) runs the identifiability checks on any D/group/W triple. model_setup() calls it internally, but it is exported so that hand-built designs can be checked too (Section 2.5). group must be a factor, colnames(D) must equal names(W), and each W[[k]] must have one row per level of group, in level order. With verbose = TRUE it reports progress through message().

2.1 The sampler design

dim(design$D)
#> [1] 551   3
colnames(design$D)
#> [1] "(Intercept)"     "distracted_ppvt" "distracted_a1"
design$group_name
#> [1] "school_id"
nlevels(design$group)
#> [1] 47

2.2 The hyper-design matrices

lapply(design$W, colnames)
#> $`(Intercept)`
#> [1] "(Intercept)"        "private_school"     "title1"            
#> [4] "free_reduced_lunch"
#> 
#> $distracted_ppvt
#> [1] "(Intercept)"
#> 
#> $distracted_a1
#> [1] "(Intercept)"        "free_reduced_lunch"
design$popef.moderation
#>                    interaction_col          moderator  random_slope
#> 1 free_reduced_lunch:distracted_a1 free_reduced_lunch distracted_a1

These follow the table in Section 1.1. The random intercept is regressed on all three school-level covariates. The distracted_a1 slope has its population mean plus the moderator free_reduced_lunch. The distracted_ppvt slope has only its population mean.

2.3 The reference fit

class(design$lmer)
#> [1] "lmerMod"
#> attr(,"package")
#> [1] "lme4"
design$Psi
#>     (Intercept) distracted_ppvt   distracted_a1 
#>       199.28938        22.53904        42.06903
design$dispersion
#> [1] 357.2463

Psi holds the classical between-school variance \(\hat\tau^2_k\) of each random coefficient, and dispersion holds the residual \(\hat\sigma^2\). Set fit_mer = FALSE to skip the reference fits:

design_only <- model_setup(form_lmer, data = dat, fit_mer = FALSE)
is.null(design_only$lmer)
#> [1] TRUE

2.4 Identifiability at two levels

Level 1 (within each group). Can the coefficient vector \(b_j\) be estimated from group \(j\)’s data alone?

• groupef.rank: TRUE when the rows of D for group \(j\) have full column rank.
• groupef.estimable: full rank and a family-specific check passes. For Gaussian models the group needs more observations than coefficients (\(n_j > p_{\mathrm{re}}\)) so its residual variance is estimable. For binomial, Poisson and Gamma models a classical glm() is fitted to the group and must give finite estimates. This catches, for example, complete separation or an all-zero count response.
• groupef.glm_check: a data frame with a short note for each group that failed the family-specific check.

Level 2 (across groups). For each random coefficient \(k\), popef.rank[k] is TRUE when W[[k]], restricted to the Level 1 estimable groups, has full column rank. popef.rank_ok is TRUE when this holds for every \(k\).

table(design$groupef.rank)
#> 
#> FALSE  TRUE 
#>    13    34
table(design$groupef.estimable)
#> 
#> FALSE  TRUE 
#>    13    34
design$popef.rank
#>     (Intercept) distracted_ppvt   distracted_a1 
#>            TRUE            TRUE            TRUE
design$popef.rank_ok
#> [1] TRUE

Groups that fail Level 1 stay in the reference fit. Prior_Setup_GLMM() needs popef.rank_ok = TRUE to derive default hyperpriors. The Part 4 examples keep only full-rank groups before calibrating priors:

full_rank <- names(design$groupef.rank)[design$groupef.rank]
dat_fr <- subset(dat, school_id %in% full_rank)
dat_fr$school_id <- droplevels(dat_fr$school_id)

design_fr <- model_setup(form_lmer, data = dat_fr, fit_mer = FALSE)
all(design_fr$groupef.rank)
#> [1] TRUE
nlevels(design_fr$group)
#> [1] 34

2.5 Calling check_identifiability() directly

set.seed(1)
J <- 8L; n_j <- 6L
group <- factor(rep(paste0("g", seq_len(J)), each = n_j))
x1 <- rnorm(J * n_j)
D <- cbind("(Intercept)" = 1, x1 = x1)
y <- 1 + x1 + rnorm(J * n_j)
W <- list(
  "(Intercept)" = matrix(1, J, 1, dimnames = list(NULL, "(Intercept)")),
  x1            = matrix(1, J, 1, dimnames = list(NULL, "(Intercept)"))
)

ident <- check_identifiability(
  y = y, D = D, group = group, W = W,
  family = gaussian(), group_name = "group", verbose = TRUE
)
#> check_identifiability: Level 1 (rank) -- 8 of 8 group full column rank.
#> check_identifiability: Level 1 (estimability) -- 8 of 8 full-rank group are estimable.
#> check_identifiability: Level 2 -- all hyper-design matrices are full rank (8 estimable group).
ident$popef.rank_ok
#> [1] TRUE

Now make one group’s covariate constant. Its slope can no longer be estimated from that group alone, so it fails Level 1:

D_bad <- D
D_bad[group == "g3", "x1"] <- 0.5
ident_bad <- check_identifiability(
  y = y, D = D_bad, group = group, W = W, group_name = "group"
)
ident_bad$groupef.rank
#>    g1    g2    g3    g4    g5    g6    g7    g8 
#>  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE

3.1 Call tree

model_setup(formula, data, family, ..., fit_mer, dispformula)
 │
 ├─ .lmebayes_normalize_family()              family string/function → family object
 │
 ├─ extract_re_hyper_matrices()               ── STAGE A: sampler design
 │   ├─ is_single_factor_model()              exactly one grouping factor?
 │   ├─ .lmebayes_validate_uncorrelated_re_formula()
 │   │   └─ is_fixed_effects_only()           reject (1 || g) and correlated (x | g)
 │   ├─ lme4::lFormula()                      model frame, X, reTrms
 │   ├─ classify_lme4_fixed_columns()         split X into level-1 / level-2 columns
 │   ├─ classify_crosslevel_re_moderation()   level-1 columns → mean slopes,
 │   │                                        cross-level moderators, or disallowed
 │   ├─ get_lme4_components()                 X, sparse Z and its row/column maps
 │   │   └─ .lme4_Z_random_column_map(), .lme4_Z_random_row_map(),
 │   │      .lme4_label_Z_random_sparse()
 │   ├─ extract_lme4_fixed_group_matrix()     one row per group: level-2 X  → W
 │   │   └─ classify_lme4_fixed_columns()
 │   ├─ extract_re_Z_obs()                    sparse Z → dense n x p_re matrix D
 │   └─ .lmebayes_weights_offset_from_frame() weights/offset aligned with y
 │       └─ .lmebayes_normalize_weights(), .lmebayes_normalize_offset()
 │
 ├─ .lmebayes_dispformula_kind()              "pooled" (~1) or "group" (~<group>)
 │
 ├─ [fit_mer = TRUE]                          ── STAGE B: reference fits
 │   ├─ .lmebayes_mer_optional_args()         forward only supplied args
 │   ├─ lme4::lmer() | lme4::glmer()          always (backs design$lmer / $glmer)
 │   ├─ lme4::isSingular()                    message() if singular
 │   ├─ extract_mer_variance_components()     VarCorr → Psi, dispersion
 │   └─ [gaussian & dispformula = ~<group>]
 │       └─ .lmebayes_fit_glmmtmb_reference() glmmTMB(..., dispformula)
 │
 └─ check_identifiability()                   ── STAGE C: diagnostics
     ├─ Level 1 rank:   Matrix::rankMatrix() on each D_j
     ├─ Level 1 estimability:
     │   └─ .lmebayes_block_glm_estimable()
     │       └─ .lmebayes_glm_estimable_precheck()   quick family-specific rejects
     │          then stats::glm() per group           (binomial/poisson/Gamma)
     └─ Level 2 rank:   Matrix::rankMatrix() on W[[k]][estimable groups, ]

3.2 Helper reference

Stage A — design extraction (R/lme4_design_utilities.R unless noted)

Stage B — reference fits

Stage C — identifiability (R/check_identifiability.R)

3.3 Exploring the internals

The non-dot design utilities are documented internally and can be called with :::. For example, to see how the fixed columns were classified:

fr <- lme4::lFormula(form_lmer, data = dat_fr)
cls <- glmbayesCore:::classify_lme4_fixed_columns(fr$X, fr$reTrms$flist[[1]])
str(cls, max.level = 1)
#> List of 3
#>  $ level1_cols : chr [1:3] "distracted_ppvt" "distracted_a1" "free_reduced_lunch:distracted_a1"
#>  $ level2_cols : chr [1:4] "(Intercept)" "private_school" "title1" "free_reduced_lunch"
#>  $ level2_flags: Named logi [1:7] TRUE TRUE TRUE TRUE FALSE FALSE ...
#>   ..- attr(*, "names")= chr [1:7] "(Intercept)" "private_school" "title1" "free_reduced_lunch" ...