Skip to contents

Random Effects

Gilles Colling

2026-10-13

Source: vignettes/random-effects.Rmd
random-effects.Rmd

Installation 

install.packages("pak")                   # development version
pak::pak("gcol33/tulpaRatio")

Overview

This vignette covers random effects in tulpaRatio, including random intercepts, random slopes, crossed effects, and nested effects. tulpaRatio uses lme4-style formula syntax with extensions for correlated vs. uncorrelated slopes.

The examples use one simulated data set: counts and sampling effort from three plots at each of 16 sites in four regions, surveyed in six years. Region, site, plot and observation-level effects enter both the count and the effort; year and observer effects and the site-specific slopes of the covariates enter only the count.

library(tulpaRatio)

set.seed(42)
df <- expand.grid(plot = 1:3, site = 1:4, region = 1:4, year = 1:6)
df$site <- factor(paste(df$region, df$site, sep = "_"))
df$plot <- factor(df$plot)
df$region <- factor(df$region)
df$year <- factor(df$year)
n <- nrow(df)
df$depth <- rnorm(n)
df$temp <- rnorm(n)
df$x <- rnorm(n)
df$z <- rnorm(n)
df$observer <- factor(sample(1:6, n, replace = TRUE))

b_region <- rnorm(4, 0, 0.4)[df$region]
b_site <- rnorm(16, 0, 0.5)[df$site]
b_plot <- rnorm(48, 0, 0.3)[interaction(df$site, df$plot)]
b_year <- rnorm(6, 0, 0.4)[df$year]
b_observer <- rnorm(6, 0, 0.3)[df$observer]
slopes <- matrix(rnorm(16 * 5, 0, 0.2), nrow = 16)[df$site, ]
u <- rnorm(n, 0, 0.4)
df$effort <- rgamma(n, shape = 5, rate = 5 / exp(0.5 + b_region + b_site + b_plot + u))
df$count <- rpois(n, exp(
  1 + b_region + b_site + b_plot + b_year + b_observer + u +
    (0.4 + slopes[, 1]) * df$depth + (-0.2 + slopes[, 2]) * df$temp +
    (0.3 + slopes[, 3]) * df$x + (-0.2 + slopes[, 4]) * df$z +
    (0.1 + slopes[, 5]) * df$x * df$z
))

ctl <- list(seed = 1, adapt_delta = 0.95)

Random Intercepts

Single Grouping Factor 

# Random intercept for site
fit <- tratio(
  count | effort ~ depth + (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
summary(fit)
#> ratio model summary
#> ===================
#> 
#> Inference: Exact (Tier 1) via hmc
#> Family: poisson_gamma 
#> 
#> Fixed effects (numerator):
#>        parameter  mean    sd  2.5% 97.5%  rhat ess_bulk
#>  num:(Intercept) 0.958 0.221 0.503 1.376 1.017      200
#>        num:depth 0.303 0.035 0.230 0.370 1.002     1102
#> 
#> Fixed effects (denominator):
#>          parameter   mean    sd   2.5% 97.5%  rhat ess_bulk
#>  denom:(Intercept)  0.346 0.221 -0.109 0.775 1.016      200
#>        denom:depth -0.046 0.039 -0.121 0.030 1.006     1049
#> 
#> Variance components:
#>  parameter  mean    sd  2.5% 97.5%  rhat ess_bulk
#>      shape 2.421 0.195 2.049 2.809 1.007      870
#>   sigma_re 0.837 0.173 0.576 1.255 1.006      294
#> 
#> Diagnostics:
#>   Divergences: 0 
#>   Avg. acceptance: 0.954 
#>   Chains: 4 
#>   Iterations: 2000 
#>   Warmup: 1000 
#>   R-hat > 1.01 for 2 parameter(s), largest 1.017 (num:(Intercept)).
#>   Bulk ESS < 400 for 3 parameter(s), smallest 200 (denom:(Intercept)).
#>   Run longer chains or reparameterise before using these estimates.

This estimates a site-specific intercept bsite[i]b_{\text{site}[i]} for each observation ii:

log⁡(μi)=β0+β1⋅depthi+bsite[i] \log(\mu_i) = \beta_0 + \beta_1 \cdot \text{depth}_i + b_{\text{site}[i]}

where bj∼Normal(0,σsite2)b_j \sim \text{Normal}(0, \sigma_{\text{site}}^2).

Multiple (Crossed) Grouping Factors

When observations are classified by multiple independent grouping factors:

# Crossed random effects: site and year
fit <- tratio(
  count | effort ~ depth + (1 | site) + (1 | year),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

This estimates: - bjsite∼Normal(0,σsite2)b^{\text{site}}_j \sim \text{Normal}(0, \sigma_{\text{site}}^2) - bkyear∼Normal(0,σyear2)b^{\text{year}}_k \sim \text{Normal}(0, \sigma_{\text{year}}^2)

The linear predictor becomes:

log⁡(μi)=β0+β1⋅depthi+bsite[i]site+byear[i]year \log(\mu_i) = \beta_0 + \beta_1 \cdot \text{depth}_i + b^{\text{site}}_{\text{site}[i]} + b^{\text{year}}_{\text{year}[i]}

Nested Random Effects

When groups are hierarchically nested (plots within sites):

# Nested: plots within sites
fit_nested <- tratio(
  count | effort ~ depth + (1 | site/plot),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

This expands to:

# Equivalent explicit form
fit <- tratio(
  count | effort ~ depth + (1 | site) + (1 | site:plot),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

The site:plot term groups by each plot-within-site combination.

Deep Nesting

For deeper hierarchies:

# region/site/plot expands to:
# (1 | region) + (1 | region:site) + (1 | region:site:plot)
# Every region, site and plot holds several observations, so each effect
# is well informed by the data and the centred parameterisation suits it
fit <- tratio(
  count | effort ~ x + (1 | region/site/plot),
  data = df,
  family = ratiod_poisson_gamma(),
  control = modifyList(ctl, list(re_param = "centered"))
)

Random Slopes

Random slopes allow the effect of a predictor to vary by group.

Correlated Intercept and Slope

Use single | for correlated random effects:

# Intercept and depth effect vary by site, with correlation
fit_slope <- tratio(
  count | effort ~ depth + (1 + depth | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

This estimates:

(b0jb1j)∼MVN((00),(σ02ρσ0σ1ρσ0σ1σ12)) \begin{pmatrix} b_{0j} \\ b_{1j} \end{pmatrix} \sim \text{MVN}\left(\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_0^2 & \rho\sigma_0\sigma_1 \\ \rho\sigma_0\sigma_1 & \sigma_1^2 \end{pmatrix}\right)

The correlation ρ\rho captures whether sites with higher intercepts also tend to have steeper slopes. The summary reports it as cor_re[1,intercept:depth] beside the two standard deviations:

coefs <- coef(summary(fit_slope))
coefs[grepl("^(sigma_re|cor_re)", coefs$parameter), c("parameter", "mean", "sd")]
#>                   parameter        mean         sd
#> 5 cor_re[1,intercept:depth] -0.07041415 0.27823696
#> 7         sigma_re[1,depth]  0.18419592 0.05326567
#> 8     sigma_re[1,intercept]  0.85023296 0.18423371

Uncorrelated Intercept and Slope

Use double || for independent random effects:

# Intercept and depth effect vary by site, independently
fit <- tratio(
  count | effort ~ depth + (1 + depth || site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

This is computationally simpler and appropriate when you don’t expect correlation between intercepts and slopes.

Multiple Slopes 

# Multiple correlated slopes
fit <- tratio(
  count | effort ~ depth + temp + (1 + depth + temp | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

# Multiple uncorrelated slopes
fit <- tratio(
  count | effort ~ depth + temp + (1 + depth + temp || site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

Formula Expansion

tulpaRatio supports lme4-style formula expansion for random slopes.

Interactions 

# x*z expands to x + z + x:z
fit <- tratio(
  count | effort ~ x*z + (1 + x*z | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)
# Random slopes for x, z, and their interaction

Polynomials 

# Quadratic random slope
fit <- tratio(
  count | effort ~ x + I(x^2) + (1 + x + I(x^2) | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

Parameter Naming

tulpaRatio cleans parameter names for readability:

Formula Term Parameter Name
x:z x_z

Shared Random Effects

By default, tulpaRatio shares random effects between numerator and denominator when the same grouping factor appears in both. This is the core of the tulpaRatio philosophy (see vignette("philosophy")).

Default: Shared Structure 

# Site effect shared between count and effort processes
fit <- tratio(
  count | effort ~ (1 | site),
  formula_num = ~ depth,
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

The site effect enters both linear predictors identically:

log⁡(μcount)=β0+β1⋅depth+𝐛𝐬𝐢𝐭𝐞 \log(\mu_{\text{count}}) = \beta_0 + \beta_1 \cdot \text{depth} + \mathbf{b_{\text{site}}} log⁡(μeffort)=α0+𝐛𝐬𝐢𝐭𝐞 \log(\mu_{\text{effort}}) = \alpha_0 + \mathbf{b_{\text{site}}}

Shared effects cancel in the ratio, removing spurious variation from common unmeasured drivers.

Explicit Shared Structure 

# Only site is shared, observer is process-specific
fit <- tratio(
  count | effort ~ (1 | site) + (1 | observer),
  shared = ~ (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = ctl
)

Independent Structure (Use with Caution) 

set.seed(7)
trials_df <- data.frame(
  block = factor(rep(1:8, each = 6)),
  treatment = factor(rep(c("control", "treated"), 24)),
  trials = 20
)
b_block <- rnorm(8, 0, 0.5)[trials_df$block]
trials_df$successes <- rbinom(
  48, size = 20,
  prob = plogis(-0.5 + 0.8 * (trials_df$treatment == "treated") + b_block)
)

# Request independence (triggers warning)
fit <- tratio(
  successes | trials ~ treatment + (1 | block),
  shared = ~ 0,
  data = trials_df,
  family = ratiod_binomial(),
  control = ctl
)
#> Warning: You have specified `shared = ~ 0`, assuming independence between numerator and denominator processes.
#> 
#> This assumption is often violated in practice and can lead to spurious ratio effects when both processes share unmeasured drivers.
#> 
#> Consider whether shared random effects would be more appropriate.

The warning is issued every time shared = ~ 0 is requested, because independence between the two processes is an assumption the data rarely support.

Extracting Random Effects

Random-effect draws are named re[term,group], with terms numbered in formula order (for site/plot, term 1 is site and term 2 is site:plot) and a third index naming the component of a random slope. as_draws() returns them in the posterior format.

# Posterior summaries for random effects
re <- posterior::subset_draws(as_draws(fit_nested), variable = "re")
head(posterior::summarise_draws(re))
#> # A tibble: 6 × 10
#>   variable   mean median    sd   mad      q5   q95  rhat ess_bulk ess_tail
#>   <chr>     <dbl>  <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 re[1,1]   0.786  0.788 0.310 0.311  0.280  1.28   1.01     437.    1077.
#> 2 re[1,2]   0.343  0.345 0.313 0.307 -0.169  0.867  1.01     559.    1084.
#> 3 re[1,3]   0.418  0.419 0.312 0.317 -0.0928 0.928  1.01     474.    1070.
#> 4 re[1,4]   0.364  0.357 0.315 0.301 -0.151  0.892  1.01     578.     946.
#> 5 re[1,5]   0.789  0.780 0.301 0.299  0.314  1.30   1.01     561.    1157.
#> 6 re[1,6]  -0.395 -0.393 0.317 0.316 -0.920  0.118  1.00     621.    1377.

# For random slopes: intercept and depth components of the first site
re_slope <- posterior::subset_draws(
  as_draws(fit_slope),
  variable = c("re[1,1,intercept]", "re[1,1,depth]")
)
posterior::summarise_draws(re_slope)
#> # A tibble: 2 × 10
#>   variable       mean median     sd    mad     q5    q95  rhat ess_bulk ess_tail
#>   <chr>         <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
#> 1 re[1,1,inte…  0.867  0.876 0.242  0.231   0.451  1.25   1.03     184.     473.
#> 2 re[1,1,dept… -0.278 -0.276 0.0979 0.1000 -0.441 -0.120  1.00    1129.    1796.

Prior Specification

Default PC Priors

tulpaRatio uses penalized complexity (PC) priors that favor simpler models:

# Default: P(sigma > 1) = 0.01
fit <- tratio(
  count | effort ~ depth + (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  priors = ratiod_priors(),
  control = ctl
)

Custom Priors 

# Expect smaller random effects
fit <- tratio(
  count | effort ~ depth + (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  priors = ratiod_priors(
    sigma = prior_half_cauchy(0.5)  # Half-Cauchy with scale 0.5
  ),
  control = ctl
)

Latent Factors for Unmeasured Confounders

When you suspect unmeasured confounders are driving both numerator and denominator but don’t have a natural grouping structure, latent factors provide an alternative to random effects.

Single Latent Factor 

# Capture shared unmeasured variation
# Observation-level factors mix slowly, so the chains run longer
fit_latent <- tratio(
  count | effort ~ (1 | site),
  formula_num = ~ depth,
  data = df,
  family = ratiod_poisson_gamma(),
  latent = latent_factor(n_factors = 1),
  control = list(seed = 1)
)
head(latent_factors(fit_latent))
#>   observation factor        mean        sd        q2.5         q50      q97.5
#> 1           1      1 -1.44020897 0.6014386 -2.67137054 -1.41283988 -0.3321069
#> 2           2      1  1.09839661 0.4033631  0.27574066  1.10558310  1.8680457
#> 3           3      1 -0.05163295 0.4786275 -1.04969767 -0.03596029  0.8303736
#> 4           4      1 -1.11649373 0.6688172 -2.49970142 -1.11253748  0.1268144
#> 5           5      1  0.82943143 0.4523220 -0.09131511  0.84582517  1.6877243
#> 6           6      1 -0.35728172 0.6157591 -1.63180607 -0.33230960  0.7811181

Latent factors are observation-level scores on which each process loads with its own loading:

log⁡(μcount)=β0+β1⋅depth+bsite+λ1fi \log(\mu_{\text{count}}) = \beta_0 + \beta_1 \cdot \text{depth} + b_{\text{site}} + \lambda_1 f_i log⁡(μeffort)=α0+bsite+λ2fi \log(\mu_{\text{effort}}) = \alpha_0 + b_{\text{site}} + \lambda_2 f_i

where fi∼Normal(0,1)f_i \sim \text{Normal}(0, 1) with a sum-to-zero constraint and λ1>0\lambda_1 > 0. The factor adds the covariance λ1λ2\lambda_1 \lambda_2 between an observation’s two linear predictors and moves the ratio by (λ1−λ2)fi(\lambda_1 - \lambda_2) f_i.

Two Latent Factors

Two factors with a lower-triangular loading matrix give the two linear predictors an unrestricted residual covariance. A third factor would have nothing left to identify, so latent_factor() stops at two (one when shared = FALSE).

fit <- tratio(
  count | effort ~ depth,
  data = df,
  family = ratiod_poisson_gamma(),
  latent = latent_factor(n_factors = 2),
  control = list(seed = 1)
)

Latent Factors vs Random Effects

Feature Random Effects Latent Factors
Grouping Required (site, year, etc.) None needed
Level Group-level Observation-level
Use case Known grouping structure Residual confounding
Dimensionality 1 per group K factors × N observations

Use latent factors when: - You suspect unmeasured confounders but don’t know the grouping structure - Residual variation after random effects is still correlated - You want to capture observation-level shared structure

Prior Specification 

# Stronger regularization on latent factors
fit <- tratio(
  count | effort ~ depth,
  data = df,
  family = ratiod_poisson_gamma(),
  latent = latent_factor(
    n_factors = 1,
    prior = prior_pc(U = 0.5, alpha = 0.01)  # P(loading > 0.5) = 0.01
  ),
  control = list(seed = 1)
)

See Also