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 intercept for site
fit <- tratio(
count | effort ~ depth + (1 | site),
data = df,
family = ratiod_poisson_gamma(),
control = ctl
)
summary(fit)This estimates a site-specific intercept \(b_{\text{site}[i]}\) for each observation \(i\):
\[ \log(\mu_i) = \beta_0 + \beta_1 \cdot \text{depth}_i + b_{\text{site}[i]} \]
where \(b_j \sim \text{Normal}(0, \sigma_{\text{site}}^2)\).
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: - \(b^{\text{site}}_j \sim \text{Normal}(0, \sigma_{\text{site}}^2)\) - \(b^{\text{year}}_k \sim \text{Normal}(0, \sigma_{\text{year}}^2)\)
The linear predictor becomes:
\[ \log(\mu_i) = \beta_0 + \beta_1 \cdot \text{depth}_i + b^{\text{site}}_{\text{site}[i]} + b^{\text{year}}_{\text{year}[i]} \]
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.
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 allow the effect of a predictor to vary by group.
# 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
)tulpaRatio supports lme4-style formula expansion for random slopes.
tulpaRatio cleans parameter names for readability:
| Formula Term | Parameter Name |
|---|---|
x:z |
x_z |
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))
# 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)tulpaRatio uses penalized complexity (PC) priors that favor simpler models:
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.
# 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))Latent factors are observation-level scores on which each process loads with its own loading:
\[ \log(\mu_{\text{count}}) = \beta_0 + \beta_1 \cdot \text{depth} + b_{\text{site}} + \lambda_1 f_i \] \[ \log(\mu_{\text{effort}}) = \alpha_0 + b_{\text{site}} + \lambda_2 f_i \]
where \(f_i \sim \text{Normal}(0, 1)\) with a sum-to-zero constraint and \(\lambda_1 > 0\). The factor adds the covariance \(\lambda_1 \lambda_2\) between an observation’s two linear predictors and moves the ratio by \((\lambda_1 - \lambda_2) f_i\).
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).
| 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
vignette("getting-started") - Basic usagevignette("philosophy") - Why shared structure
mattersvignette("spatial-temporal") - Spatial and temporal
extensions