Chapter 05: Reading the output

library(lmebayes)

Chapter 05: Reading the output

An lmerb fit is a list with about thirty elements. This chapter is the map: which accessor gives which quantity, what the raw containers hold, and where the naming departs from lme4 in ways that will surprise you.

data(sleepstudy, package = "lme4")
dat <- sleepstudy
dat$Days_c <- dat$Days - mean(dat$Days)
form <- Reaction ~ Days_c + (1 + Days_c || Subject)

ps <- Prior_Setup_GLMM(form, data = dat, pop.pwt = 0.01)

set.seed(1)
fit <- lmerb(
  form, data = dat,
  pfamily_list     = pfamily_list(ps),
  dispersion_ranef = ps$group.dispersion,
  n                = 500L,
  progbar          = FALSE
)

Two blocks, two prefixes

Everything in the output is named after which block it belongs to.

Prefix Block Contents
popef.* Block 2 Population coefficients \(\gamma\) and their variance components
groupef.* Block 1 Per-group coefficients \(\beta_j\)
group.dispersion* — Observation dispersion \(\sigma^2\)

Within each prefix the suffix says what kind of quantity it is.

Suffix Meaning
none The raw draws
.mode Point estimate before sampling (exact posterior mean for Gaussian, mode otherwise)
.means / .mean Average over the stored draws
.iters Accept-reject candidates per draw (sampler efficiency) 
names(fit)
#>  [1] "call"                        "formula"                    
#>  [3] "lmer"                        "glmmTMB"                    
#>  [5] "dispformula"                 "dispersion_fit"             
#>  [7] "prior"                       "model_setup"                
#>  [9] "popef.mode"                  "popef.init"                 
#> [11] "groupef.mode"                "popef.means"                
#> [13] "popef"                       "groupef"                    
#> [15] "popef.dispersion"            "popef.dispersion.mean"      
#> [17] "popef.iters"                 "popef.iters.mean"           
#> [19] "groupef.iters"               "groupef.iters.mean"         
#> [21] "group.dispersion"            "group.dispersion.mean"      
#> [23] "group.dispersion.iters"      "group.dispersion.iters.mean"
#> [25] "draw_engine"                 "sim_method_used"            
#> [27] "m_convergence"               "pilot"                      
#> [29] "gap_tol"                     "mode_gap_max"               
#> [31] "convergence"                 "sweep_history"

print()

fit
#> Call:
#>   lmerb(formula = form, data = dat, pfamily_list = pfamily_list(ps), 
#>     dispersion_ranef = ps$group.dispersion, n = 500L, progbar = FALSE)
#> 
#> Bayesian linear mixed model  [500 draws, exact iid]
#> Formula: Reaction ~ Days_c + (1 + Days_c || Subject) 
#> 
#> Random effects (variance components fixed at lmer estimates):
#>  Groups    Name        Std.Dev.
#>  Subject   (Intercept) 37.533  
#>  Subject.1 Days_c       5.922  
#>  Residual              25.592  
#> Number of obs: 180,  groups: Subject, 18
#> 
#> --- Posterior means (ICM exact, under fixed variance components) ---
#> 
#>   RE component  parameter      fixef.mode   fixef.means    draws SD
#>   ------------  -----------  ------------  ------------  ----------
#>   (Intercept)   (Intercept)      298.5079      297.9952      9.0581
#>   Days_c        (Intercept)       10.3748       10.3326      1.4691

The header line carries three facts: how many draws, and which engine produced them. exact iid means the joint posterior was multivariate normal and was sampled directly. two-block Gibbs means sweeps were needed; Chapter 06 covers that case.

The table’s three columns are worth checking every time. fixef.mode is computed analytically before sampling; fixef.means is the Monte Carlo average. They should agree to within draws SD / sqrt(n). If they do not, something is wrong with the sampler, not with your model.

summary()

summary(fit)
#> Call:
#>   lmerb(formula = form, data = dat, pfamily_list = pfamily_list(ps), 
#>     dispersion_ranef = ps$group.dispersion, n = 500L, progbar = FALSE)
#> 
#> Bayesian linear mixed model fit  [500 draws, exact iid]
#> Formula: Reaction ~ Days_c + (1 + Days_c || Subject) 
#> 
#> Random effects (variance components fixed at lmer estimates):
#>  Groups    Name        Std.Dev.
#>  Subject   (Intercept) 37.533  
#>  Subject.1 Days_c       5.922  
#>  Residual              25.592  
#> Number of obs: 180,  groups: Subject, 18
#> 
#> === Block 2 dispersion (RE variance tau^2_k) ===
#> 
#> Prior and lmer reference:
#> 
#>               Prior 1/E[1/tau2]   E[tau2] sqrt(E[tau2]) disp_lower disp_upper
#> (Intercept) dNormal   1346.1456 1346.1456       36.6899         NA         NA
#> Days_c      dNormal     29.3722   29.3722        5.4196         NA         NA
#> Residual      fixed    709.8260  709.8260       26.6426         NA         NA
#>                  lmer lmer SD
#> (Intercept) 1408.7541 37.5334
#> Days_c        35.0716  5.9221
#> Residual     654.9399 25.5918
#> 
#> Overview:
#>             Post.Mode Post.Mean Post.Sd MC Error Mean SD Pr(Prior_tail)
#> (Intercept)   1346.15   1346.15    0.00     0.00   36.69             NA
#> Days_c          29.37     29.37    0.00     0.00    5.42             NA
#> Residual       709.83    709.83      NA       NA   26.64             NA
#>             Cand/draw
#> (Intercept)         1
#> Days_c              1
#> Residual           NA
#> 
#> Distribution percentiles (tau^2):
#> 
#>                1.0%    2.5%    5.0%  Median   95.0%   97.5%   99.0%
#> (Intercept) 1346.15 1346.15 1346.15 1346.15 1346.15 1346.15 1346.15
#> Days_c        29.37   29.37   29.37   29.37   29.37   29.37   29.37
#> Residual     709.83  709.83  709.83  709.83  709.83  709.83  709.83
#> 
#> SD credible interval (sqrt(tau^2) draws; lmer SD for reference):
#> 
#>              2.5% Median 97.5% lmer SD
#> (Intercept) 36.69  36.69 36.69  37.533
#> Days_c       5.42   5.42  5.42   5.922
#> Residual    26.64  26.64 26.64  25.592
#> 
#> === Block 2: Level-2 fixed effects (hyperparameters) ===
#> 
#> Prior and lmer reference:
#> 
#>                          Prior Mean Prior.sd   lmer lmer.se
#> (Intercept)::(Intercept)     298.51    90.05 298.51   9.050
#> Days_c::(Intercept)            0.00    15.38  10.47   1.546
#> 
#> Overview:
#>                          fixef.mode fixef.means Post.Sd MC.Error Pr(Prior_tail)
#> (Intercept)::(Intercept)    298.508     297.995   9.058    0.405          0.438
#> Days_c::(Intercept)          10.375      10.333   1.469    0.066         <2e-16
#>                             
#> (Intercept)::(Intercept)    
#> Days_c::(Intercept)      ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Distribution percentiles:
#> 
#>                             1.0%    2.5%    5.0%  Median   95.0%   97.5%  99.0%
#> (Intercept)::(Intercept) 278.379 281.610 284.103 297.280 313.238 316.400 319.07
#> Days_c::(Intercept)        7.147   7.554   8.022  10.323  12.810  13.233  13.71
#> 
#> === Block 1: Random effects (group-level) ===
#> 
#> Summary of posterior mode (ranef.mode) across groups:
#> 
#>                 Mean       SD      Min       Q1   Median       Q3      Max
#> (Intercept) 298.5079  36.4724 219.4041 290.8999 304.3215 314.6673 371.8535
#> Days_c       10.4463   5.0724   0.1222   7.1416  10.4194  12.8294  19.1842
#>             MCMC.mean
#> (Intercept)    298.37
#> Days_c          10.46
#> 
#> Per-group random effects: inspect fit$groupef.mode or fit$groupef,
#>   or call summary(fit, groups = <level ids>) for selected groups.

Three sections:

  1. Block 2 dispersion — the variance components \(\tau^2_k\), and the residual variance, each against its lmer reference. When a component is dNormal its percentiles are constant, because it is fixed rather than sampled.
  2. Block 2 level-2 fixed effects — the population parameters, with prior means, lmer reference values, posterior SD, Monte Carlo error, and Pr(Prior_tail).
  3. Block 1 random effects — a distributional summary across groups rather than a per-group listing.

summary() takes a groups argument when you want specific levels:

summary(fit, groups = c("308", "309"))

Pr(Prior_tail)

This is a directional tail probability, not a p-value. It measures how far the posterior has moved into the tail of the prior, in the direction of the shift. Near 0.5 means the data barely moved that coefficient away from its prior mean. Near zero means the posterior sits far out in the prior’s tail — the data disagrees with the prior, decisively. The theory is in vignette("Chapter-A04", package = "glmbayes").

Accessors

fixef(), ranef(), and coef() are lme4 generics. lmebayes defines methods for them but does not re-export the generics, so call them from lme4.

ranef()

re <- lme4::ranef(fit)
head(re$Subject)
#>     (Intercept)    Days_c
#> 308    339.7664 19.173801
#> 309    219.3786  4.153587
#> 310    234.4354  7.161023
#> 330    302.3252  4.789262
#> 331    308.6786  6.538323
#> 332    307.4449  9.699681

These are not lme4 random effects. They are the full per-subject coefficients \(\beta_j\) — subject 308’s intercept is their expected reaction time, not an offset from the population mean. See Chapter 02. The lme4 quantity is \(\beta_j\) minus the population mean:

head(re$Subject) - rep(unlist(fit$popef.means), each = 6)
#>     (Intercept)     Days_c
#> 308   41.771257  8.8412503
#> 309  -78.616555 -6.1789628
#> 310  -63.559785 -3.1715271
#> 330    4.330046 -5.5432878
#> 331   10.683417 -3.7942276
#> 332    9.449768 -0.6328696

Two arguments matter. type chooses between the chain mean (default) and the ICM point estimate, and condVar attaches posterior variances:

re_mode <- lme4::ranef(fit, type = "mode")
head(re_mode$Subject, 3)
#>     (Intercept)    Days_c
#> 308    339.9486 19.184176
#> 309    219.4041  4.099875
#> 310    234.3826  7.080020

re_cv <- lme4::ranef(fit, condVar = TRUE)
head(as.data.frame(re_cv), 4)
#>    grpvar        term grp  condval   condsd
#> 1 Subject (Intercept) 308 339.7664 8.237316
#> 2 Subject (Intercept) 309 219.3786 7.741607
#> 3 Subject (Intercept) 310 234.4354 8.783966
#> 4 Subject (Intercept) 330 302.3252 8.297803

The long format has condsd filled in when condVar = TRUE, which is what you want for a caterpillar plot.

coef()

head(coef(fit)$Subject)
#>     (Intercept)    Days_c
#> 308    339.7664 19.173801
#> 309    219.3786  4.153587
#> 310    234.4354  7.161023
#> 330    302.3252  4.789262
#> 331    308.6786  6.538323
#> 332    307.4449  9.699681

For merMod objects coef() adds fixef to ranef. Here ranef() already returns full coefficients, so coef() returns the same values — kept for interface compatibility rather than because it computes anything new.

fixef()

lme4::fixef(fit)                        # type = "mer": the reference lmer fit
#> (Intercept)      Days_c 
#>   298.50789    10.46729
lme4::fixef(fit, type = "hyper")        # Block 2 posterior mode
#> (Intercept).(Intercept)      Days_c.(Intercept) 
#>               298.50789                10.37476
lme4::fixef(fit, type = "hyper.mean")   # Block 2 posterior mean
#> (Intercept).(Intercept)      Days_c.(Intercept) 
#>               297.99517                10.33255

The default is the lmer reference so that fixef() on this object behaves the way it does on a merMod. The Bayesian quantities are behind type.

Names differ between the two: the lmer version uses model-matrix names (Days_c), while the hyper versions use component.parameter (Days_c.(Intercept)), reflecting that each population parameter belongs to a specific random coefficient’s level-2 row.

The raw draws

Accessors summarize. When you want the draws themselves:

str(fit$popef)
#> List of 2
#>  $ (Intercept): num [1:500, 1] 293 308 295 297 296 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : chr "(Intercept)"
#>  $ Days_c     : num [1:500, 1] 10.64 11.48 10.38 9.81 9.03 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : chr "(Intercept)"

fit$popef is a list with one matrix per random coefficient: rows are draws, columns are that coefficient’s level-2 parameters. Because the draws are independent, ordinary Monte Carlo standard errors apply directly:

draws_int <- fit$popef[["(Intercept)"]][, 1]
c(mean = mean(draws_int),
  sd   = sd(draws_int),
  mcse = sd(draws_int) / sqrt(length(draws_int)))
#>        mean          sd        mcse 
#> 297.9951696   9.0581192   0.4050914

No effective-sample-size correction, no thinning. That is the payoff of the replicate-chain design described in vignette("Chapter-B01", package = "lmebayesCore").

head(fit$groupef)
#>   draw Subject (Intercept)    Days_c
#> 1    1     308    332.8099 23.359618
#> 2    1     309    221.8328  2.043515
#> 3    1     310    238.1080  9.044810
#> 4    1     330    307.4362  3.949451
#> 5    1     331    321.0254  7.489291
#> 6    1     332    301.4832  4.096830
nrow(fit$groupef)
#> [1] 9000

fit$groupef is long format: n draws times J groups. Filter it directly, or use the printer:

lmebayesCore::print_groupef(fit, draws = 1:3, groups = 1:3)
#> 
#> Simulated group coefficients:
#>  draw Subject (Intercept) Days_c
#>     1     308       332.8 23.360
#>     1     309       221.8  2.044
#>     1     310       238.1  9.045
#>     2     308       339.1 18.781
#>     2     309       225.6  5.786
#>     2     310       229.2  5.505
#>     3     308       340.4 17.664
#>     3     309       214.5  3.752
#>     3     310       243.9  3.150

print_groupef() lives in lmebayesCore and is not re-exported here, so it needs the namespace prefix. It takes the fit, not the data frame.

Credible intervals

Nothing special is required — these are ordinary independent draws.

t(apply(fit$popef[["(Intercept)"]], 2, quantile,
        probs = c(0.025, 0.5, 0.975)))
#>                 2.5%    50%    97.5%
#> (Intercept) 281.6104 297.28 316.4002
t(apply(fit$popef[["Days_c"]], 2, quantile,
        probs = c(0.025, 0.5, 0.975)))
#>                 2.5%      50%    97.5%
#> (Intercept) 7.554183 10.32334 13.23322

Per-group intervals come from the long-format container:

by_subject <- split(fit$groupef$Days_c, fit$groupef$Subject)
head(t(sapply(by_subject, quantile, probs = c(0.025, 0.5, 0.975))), 5)
#>           2.5%       50%     97.5%
#> 308 13.9003503 19.147798 24.252926
#> 309 -0.6780557  4.285271  8.954318
#> 310  2.2318144  7.001911 12.211377
#> 330 -0.3958314  4.878103  9.915853
#> 331  1.8899615  6.538472 10.979420

Dispersion and convergence metadata

fit$group.dispersion       # sigma^2 (fixed here)
#> [1] 709.826
fit$m_convergence          # inner Gibbs sweeps per stored draw
#> [1] 1
fit$sim_method_used        # which route ran
#> [1] "DEFAULT"

summary_sigma2() gives a full report, but only when the observation dispersion is per-group — with a single fixed sigma^2 there is nothing to report and it stops with a message. Chapter 08 covers that case.

Quick reference

You want Use
Population posterior mean fit$popef.means or fixef(fit, "hyper.mean")
Population draws fit$popef[[k]]
Per-group coefficients lme4::ranef(fit) or fit$groupef.mode
Per-group draws fit$groupef (long format)
Observation variance fit$group.dispersion
Variance components summary(fit), Block 2 dispersion section
Sweeps per draw fit$m_convergence
Reference lmer fit fit$lmer
The parsed design fit$model_setup

See also

  • Chapter 06 — m_convergence and the sweep history.
  • Chapter 08 — summary_sigma2() and per-group dispersion.
  • ?lmerb, ?ranef.lmerb, ?fixef.lmerb