Student-t distributed random effects

(x | gr(g, dist = "student")) gives a grouping term a Student-t latent instead of a gaussian one, which is brms's spelling (brms::gr(), argument dist). A group far from the others then costs the variance component much less than it does under a gaussian latent, because the t's tail can hold it.

Value

gr(g, dist = "student") is a formula term, not a free-standing function: bf() reads it at parse time, and the value it contributes is the heavy-tailed random-effect block of the model, reachable through ranef() and VarCorr(). This page itself documents the term grammar and returns nothing.

Parameterization

A level's coefficients are drawn as b_j = sqrt(nu * u_j) W z_j with u_j ~ inv-chi2(nu) and z_j ~ N(0, I), which is the multivariate t with nu degrees of freedom and scale matrix Sigma = W W'. The mixing variable u_j is per LEVEL and shared across that level's coefficients, so a correlated or a diag() block is one multivariate t and not several independent univariate ones. So is a || term, which is one diag() block. brms builds one term the same way.

All the Student-t terms of one grouping factor also share one u_j per level, as in brms. An example is (1 | gr(g, dist = "student")) + (0 + x | gr(g, dist = "student")). A term in mu and a term in sigma also share it. At each level, the coefficients of these terms are one multivariate t, with a block-diagonal scale matrix. So the example is the model (1 + x || gr(g, dist = "student")).

A gr(g, by = f, dist = "student") term has one block per level of f. Each block holds different levels of g, so alone the blocks share nothing. Beside a plain Student-t term on g, each level of g shares one u_j with the by-term's block that holds it.

A term on another grouping factor has its own u_j. A gaussian term on g has none. The Student-t terms of one factor take one dist_nu, because brms has one df per factor. simulate(re_formula = ) keeps or redraws the terms together at the levels they share. quadrature = TRUE refuses terms that share a level, because it integrates one random effect at a time.

Sigma is the SCALE, not the covariance. The variance is Sigma * nu / (nu - 2), so a standard deviation is scale * sqrt(nu / (nu - 2)). VarCorr() reports the scale in its sd rows, as brms's sd_ is; confint(), variables() and frm_simulate(newparams = ) speak of it as sd_<group>__<term>, which is the name brms gives the same quantity. The correlations are the same either way.

Why `nu` is fixed

brms estimates nu under a gamma(2, 0.1) prior truncated at 1, and that prior is carrying the parameter. Maximum likelihood has no such help: profiling a simulated fit with 20 groups leaves the whole grid from 2.1 to 500 inside the 95% profile interval, and joint ML sends nu to a boundary in 24% to 41% of replicates at 20 groups and still 5% to 14% at 100 (dev/tre-feasibility.md). So nu is a constant here, set by dist_nu and defaulting to 5. That is the frequentist analogue of brms's own prior(constant(3), class = "df"). To ask what the data say about it, fit two or three values and compare logLik().

dist_nu must exceed 2, so that the latent has a variance for VarCorr() to convert and for new-level prediction to use.

Accuracy

The t density is not log-concave, so the Laplace approximation is not exact over it as it is over a gaussian latent. Measured against adaptive quadrature (dev/tre-feasibility.md), the approximation pushes the estimated latent SCALE UPWARD, and how much depends on how much the data say about each level:

With the latent scale near the residual SD, the bias is under 2% of one standard error at 8 observations per group and 0.2% at 25, at every nu from 2.5 up. It reaches 12% of a standard error in the worst case tested, nu = 2.5 with 3 observations per group.
It becomes material only where the variance component is small AND the groups are tiny: at 2 observations per group and a true scale a quarter of the residual SD, the Laplace estimate of the scale came back three times the exact one.

Two checks, both already in the package. quadrature = TRUE marginalizes a scalar random intercept by Gauss-Kronrod quadrature instead, which over a t latent is EXACT, not merely better; it is the one to run when a t block's variance component matters and the groups are small. frmtmb.sample::check_laplace() measures the same thing without refitting, by NUTS on the objective: its z_shift for theta reproduced the displacement above to within a percentage point in the probe.

A last consequence of the constant: logLik(), AIC() and BIC() carry roughly G * c(nu) where G is the number of levels and c(nu) = lgamma((nu+1)/2) - lgamma(nu/2) + log(2/(nu+1))/2 (-0.226 at nu = 3, -0.141 at nu = 5). Comparing two t fits with the same nu and the same grouping is fine, because the offset cancels; comparing a t fit against a gaussian one by AIC is not.

What is refused

gr(cov = ) / gr(prec = ): A relationship matrix correlates the LEVELS, and the t's mixing variable is per level, so the joint density over the field is not a multivariate t and has no closed form to hand the Laplace machinery. brms writes this combination, but as a hierarchical construction Stan samples rather than a density.
Other covariance structures: us (the default) and diag only. ar1(), cs(), toep() and the rest describe a covariance over the block's LEVELS, which the per-level mixing variable does not compose with.
mm(): A multi-membership row loads several levels at once, so the per-level mixing variable has no single value on it.
|ID| keys: Merged blocks are assembled as gaussian ones. Write the merged coefficients as one term instead - (x1 + x2 | gr(g, dist = "student")) is the same multivariate-t block.

Downstream

ranef() and its conditional variances, sdreport() standard errors, REML = TRUE and frm_bootstrap() all work: the Laplace machinery does not care which density the latent has. simulate() and frm_simulate() draw a multivariate t with one chi-square per level. predict(allow_new_levels = TRUE) inflates the unseen level's variance by nu / (nu - 2), so the interval has the right variance around a heavier-tailed truth. It is still built as a gaussian interval, so it is not the right quantile far into the tail. Use simulate() for that.

Examples

set.seed(1)
n <- 12
d <- data.frame(x = rnorm(20 * n), g = factor(rep(1:20, each = n)))
b <- rnorm(20)
b[20] <- b[20] + 6          # one outlying group
d$y <- 1 + 0.5 * d$x + b[d$g] + rnorm(20 * n)

fit_t <- frm(bf(y ~ x + (1 | gr(g, dist = "student"))),
             family = gaussian(), data = d)
fit_n <- frm(bf(y ~ x + (1 | g)), family = gaussian(), data = d)

# the gaussian latent has to widen to cover the outlying group
VarCorr(fit_t)
#> $g
#> $g$sd
#>            Estimate Est.Error      Q2.5    Q97.5
#> Intercept 0.8432083 0.2066928 0.4380979 1.248319
#> 
#> 
#> $residual__
#> $residual__$sd
#>  Estimate  Est.Error      Q2.5    Q97.5
#>  1.074199 0.05103885 0.9741644 1.174233
#> 
#> 
VarCorr(fit_n)
#> $g
#> $g$sd
#>           Estimate Est.Error     Q2.5    Q97.5
#> Intercept 1.500653 0.2475548 1.015454 1.985851
#> 
#> 
#> $residual__
#> $residual__$sd
#>  Estimate  Est.Error      Q2.5    Q97.5
#>  1.075971 0.05129669 0.9754312 1.176511
#> 
#> 

# heavier tails, at the cost of a fixed nu
frm(bf(y ~ x + (1 | gr(g, dist = "student", dist_nu = 3))),
    family = gaussian(), data = d)
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + (1 | gr(g, dist = "student", dist_nu = 3)) 
#>    Data: d (Number of observations: 240) 
#>  Method: ML   logLik: -384.804   AIC: 777.607   BIC: 791.53 
#> 
#> Multilevel Hyperparameters:
#> ~g (Number of levels: 20) 
#>               Estimate Est.Error l-95% CI u-95% CI
#> sd(Intercept)     0.67      0.18     0.39     1.14
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.28      0.20      0.9     1.66    6.56  5.5e-11
#> x             0.65      0.08      0.5     0.80    8.56  < 2e-16
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     1.07      0.05     0.98     1.18