Structured random slopes for non-Gaussian traits

Developer Note – under audit. This page is buildable and linked from Developer Notes so maintainers can audit the structured-slope surface. It documents validated syntax and point-recovery evidence for structured single-slope models, but it reports point estimates and recovery bands, not calibrated confidence intervals, and is not yet a public worked example. It is not a first-stop tutorial. Start with Get started, Formula keyword grid, and Convergence and start values before using this page.

A reaction norm is a slope: how a trait responds to an environmental gradient. The ordinary individual-level behavioural reaction-norm article is now an internal Gaussian draft for default latent(1 + x | individual, d = K) with its diagonal Psi_B,aug companion, while plain non-structured bare-bar (1 + x | individual) random slopes remain reserved. This article covers a different, internal technical surface: structured random slopes, where the slope varies across species along a phylogeny or across sites in space. Those are random regressions with an explicit dependence structure.

This article covers two directions that matter for ecological and evolutionary data:

Non-Gaussian responses. Counts, proportions, bounded scores, and ordinal categories are the rule, not the exception, in trait data. The structured random-slope grid now reaches them: recovery is established for the diagonal (indep) mode under poisson, nbinom1/nbinom2, Gamma, and Beta, and for the latent and dep modes across the full family list including binomial and ordinal_probit. Binomial and ordinal_probit under the diagonal indep mode remain structural-contract-only, with recovery not yet established.
Structured intercept + slope covariance. Beyond a single shared slope variance, the augmented forms phylo_dep(1 + x | species) and spatial_indep(1 + x | coords) estimate a per-unit random intercept and random slope together, and let you read the intercept-slope covariance back out.

The full grid – every core family (poisson, nbinom1, nbinom2, Gamma, Beta, binomial, ordinal_probit) crossed with the phylogenetic and spatial sources, across the indep, latent, and dep correlation modes – is now largely reachable: the latent and dep modes recover across the full family list, and the diagonal indep mode recovers for poisson, nbinom2, Gamma, and Beta. Binomial and ordinal_probit under indep remain structural-contract-only, with recovery evidence not yet established. This article is the worked tour of that grid, flagged where evidence stops short of recovery.

The scope is deliberately single-slope:

Claim	Status
Structured random regression with one covariate, `1 + x` (`s = 1`)	Covered for `latent`/`dep` modes across the core grid, and for `indep` mode under poisson/nbinom/Gamma/Beta; `indep`-mode binomial and ordinal_probit remain structural-only, their recovery evidence still pending
Gaussian `phylo_dep(1 + x1 + x2 \| species)` (`s = 2`)	Covered for the structured phylogenetic dependent path
Non-Gaussian `phylo_dep(..., s >= 2)`	Blocked (reserved); the fit is refused with a clear error until the evidence sweep clears
Delta, hurdle, and two-stage zero-inflated families	Blocked / out of scope for latent-scale slope covariance
Confidence intervals on slope variances	Not calibrated here; this article reports point estimates and recovery bands, not coverage guarantees

Reader path and scope

This draft is a structured-dependence map rather than a first-use tutorial:

Reader question	Model object	Code section	Readout
How does the augmented-slope grammar work in long and wide data?	`phylo_indep(1 + x \| species)`	`fit-small-long`, `fit-small-wide`	long/wide log-likelihood difference
Can a small structured Gaussian fit return usable slope scales?	diagonal phylogenetic slope fit	`recover-small`, `check-small`	fitted intercept/slope SDs and `check_gllvmTMB()` rows
Does the same grammar compose with counts?	Poisson `phylo_indep(1 + x \| species)`	`poisson-small`	convergence, Hessian flag, fitted slope scale
Where is the correlated structured slope path?	`phylo_dep(1 + x \| species)`	`dep-poisson`	exact shipped syntax and validation-cell recovery prose
Where is the spatial twin?	`spatial_indep(1 + x \| coords)`	`spatial-poisson`	exact shipped syntax and validation-cell recovery prose

Scope boundary: IN rows are the one-covariate structured random-slope grid for phylogenetic, spatial, and animal paths and fit-health diagnostics. PARTIAL rows include non-Gaussian s >= 2 dependent slopes and interval coverage for slope summaries, which is not yet calibrated. BLOCKED or out-of-scope rows are delta, hurdle, and zero-inflated slope covariance families. Treat the build-time examples as syntax and point-fit demonstrations; the heavy validation prose names recovery bands, not calibrated uncertainty.

The augmented-slope grammar

The random-slope term sits inside the RHS of the formula and follows the same dual-format grammar as every other covariance keyword in the package.

Spelling	Meaning
`phylo_indep(1 + x`	species)
`phylo_latent(1 + x`	species)
`phylo_dep(1 + x`	species)
`spatial_indep(1 + x`	coords)
`spatial_dep(1 + x`	coords)

In long format, the wide shorthand (1 + x | species) expands to the explicit stacked-trait form 0 + trait + (0 + trait):x | species; the two are byte-identical (the recovery tests assert this). In wide format, traits(t1, t2, ...) on the LHS lets the RHS use the compact (1 + x | ...) shorthand directly.

library(gllvmTMB)
library(ape)

A small worked fit: Gaussian, both shapes

We start with the lightest possible illustrative fit so the dual-format machinery is concrete before the heavier grid. Per-species intercepts and slopes are drawn jointly from $\mathcal{N}(\mathbf 0,\, \boldsymbol\Sigma_b \otimes \mathbf{A})$ on a small tree, with a positive intercept-slope correlation ( $\rho = 0.5$ ). This is a syntax-and-recovery illustration on one small data set, not a simulation study.

set.seed(2026)
n_sp <- 30
tree <- ape::rcoal(n_sp)
tree$tip.label <- paste0("sp", seq_len(n_sp))
Cphy <- ape::vcv(tree, corr = TRUE)
Lphy <- t(chol(Cphy + diag(1e-8, n_sp)))

# True 2x2 (intercept, slope) covariance, shared across traits
s2_int <- 0.4; s2_slope <- 0.3; rho <- 0.5
cov_is <- rho * sqrt(s2_int * s2_slope)
Sigma_b <- matrix(c(s2_int, cov_is, cov_is, s2_slope), 2, 2)

ab <- (Lphy %*% matrix(rnorm(n_sp * 2), n_sp, 2)) %*% chol(Sigma_b)
colnames(ab) <- c("alpha", "beta")
rownames(ab) <- tree$tip.label

# Long-format frame: trait x species x rep, x shared within a (species, rep)
trait_lvls <- paste0("t", 1:3)
sr <- expand.grid(
  species = factor(tree$tip.label, levels = tree$tip.label),
  rep = seq_len(4)
)
sr$x <- rnorm(nrow(sr))
long <- merge(sr, data.frame(trait = factor(trait_lvls, levels = trait_lvls)),
              all = TRUE)
long <- long[order(long$species, long$rep, long$trait), ]
mu_t <- c(2, 1, 0.5)[as.integer(long$trait)]
long$value <- mu_t +
  ab[as.character(long$species), "alpha"] +
  ab[as.character(long$species), "beta"] * long$x +
  rnorm(nrow(long), sd = 0.3)
head(long)
#>     species rep         x trait    value
#> 1       sp1   1 0.2252053    t1 2.303619
#> 121     sp1   1 0.2252053    t2 1.710815
#> 241     sp1   1 0.2252053    t3 1.023804
#> 31      sp1   2 1.0019599    t1 3.417987
#> 151     sp1   2 1.0019599    t2 1.883815
#> 271     sp1   2 1.0019599    t3 1.341544

The long call uses the explicit stacked-trait grammar; the wide call uses the traits() LHS shorthand. The tree is passed once at the top level via phylo_tree =. We use the diagonal phylo_indep form here because it is the lightest structured cell that still estimates both an intercept variance and a slope variance.

fit_long <- gllvmTMB(
  value ~ 0 + trait + phylo_indep(1 + x | species),
  data = long, trait = "trait", phylo_tree = tree, unit = "species"
)

The wide form is the same model on one-row-per-(species, rep) data; we reshape with tidyr::pivot_wider and confirm the two surfaces agree.

library(tidyr)
wide <- as.data.frame(tidyr::pivot_wider(
  long, id_cols = c(species, rep, x),
  names_from = trait, values_from = value
))

fit_wide <- gllvmTMB(
  traits(t1, t2, t3) ~ 1 + phylo_indep(1 + x | species),
  data = wide, phylo_tree = tree, unit = "species"
)

cat("long vs wide logLik difference:",
    signif(abs(as.numeric(logLik(fit_long)) -
               as.numeric(logLik(fit_wide))), 3), "\n")
#> long vs wide logLik difference: 0

The two surfaces are the same likelihood problem; the difference is numerical zero. The rest of the article uses whichever shape reads more clearly for the point being made, and always shows both spellings for each structured cell.

Reading the intercept-slope covariance

The augmented intercept + slope covariance comes back through extract_Sigma(). For the diagonal phylo_indep form the off-diagonal is pinned to zero by construction, so it reports the two marginal variances; the correlated phylo_dep form (below) reports the full block. The fitted slope standard deviation for this small fit:

# sd_b carries (sd_intercept, sd_slope); cor_b carries rho (0 here for indep)
sd_b <- as.numeric(fit_long$report$sd_b)
data.frame(
  component = c("sd_intercept", "sd_slope"),
  true      = signif(sqrt(c(s2_int, s2_slope)), 3),
  estimated = signif(sd_b, 3)
)
#>      component  true estimated
#> 1 sd_intercept 0.632     0.610
#> 2     sd_slope 0.548     0.534
#> 3 sd_intercept 0.632     0.641
#> 4     sd_slope 0.548     0.451
#> 5 sd_intercept 0.632     0.687
#> 6     sd_slope 0.548     0.566

Always check the fit, do not just trust convergence

A non-positive-definite Hessian (pdHess = FALSE) is an uncertainty warning, not model death: the point estimates can still be sensible, but the standard errors around them are not trustworthy. gllvmTMB surfaces this in fit$fit_health, and check_gllvmTMB() rolls the convergence, gradient, and Hessian checks into one table you should read before interpreting any structured-slope fit.

isTRUE(fit_long$fit_health$pd_hessian)
#> [1] TRUE
check_gllvmTMB(fit_long)
#>                component status     value         threshold
#> 1  optimizer_convergence   PASS         0                 0
#> 2           max_gradient   PASS 0.0008213              0.01
#> 3               sdreport   PASS      TRUE              TRUE
#> 4             pd_hessian   PASS      TRUE              TRUE
#> 5           hessian_rank   PASS     13/13         full rank
#> 6           max_fixed_se   PASS    0.2578               100
#> 7        restart_history   PASS         1              >= 1
#> 8       selected_restart   PASS         1 finite restart id
#> 9         boundary_flags   PASS      none              none
#> 10   rotation_convention   PASS      none              none
#> 11    boundary_sigma_eps   PASS    0.3003             1e-04
#>                                                        message
#> 1                               optimizer reported convergence
#> 2  largest absolute gradient component at the selected optimum
#> 3                                           sdreport available
#> 4      positive-definite Hessian for curvature-based inference
#> 5  rank of the fixed-parameter covariance matrix from sdreport
#> 6                          largest fixed-effect standard error
#> 7               number of optimizer starts recorded on the fit
#> 8                        restart selected by minimum objective
#> 9                            no simple boundary flags detected
#> 10                    no fitted latent loading matrix detected
#> 11                         estimated continuous residual scale
#>                                                                                    action
#> 1            try multiple starts, stronger starts, rescaling, or an alternative optimizer
#> 2                    tighten optimization, rescale predictors, or inspect weak components
#> 3                   use point summaries cautiously and prefer profile/bootstrap intervals
#> 4        check gradients, boundary variances, rank, starts, and profile/bootstrap targets
#> 5                                    treat rank loss as a Hessian/identifiability warning
#> 6                         check collinearity, scaling, or weakly identified fixed effects
#> 7                                    refit with current gllvmTMB if provenance is missing
#> 8                                 inspect restart_history for competing likelihood basins
#> 9                     still inspect profile/bootstrap output for target-specific weakness
#> 10                                                  no loading rotation diagnostic needed
#> 11 if estimated near zero, check row-level unique terms or residual-scale identifiability

If pd_hessian is FALSE, the usual remedies are more units (more species, more sites), more replication per unit, a different seed for the starting values, or stepping down to a lighter correlation mode (dep $\to$ latent $\to$ indep). Report the point estimate with the caveat; do not silently drop the fit.

A non-Gaussian reaction norm: Poisson counts

Count traits are where random regression earns its keep in ecology – fecundity, abundance, and visit counts all vary in their environmental response across species. The slope term composes with the Poisson family exactly as it does with Gaussian: the augmented (intercept, slope) contribution enters the linear predictor before the family dispatch, so swapping family = is the only change.

Here is a small illustrative Poisson fit that runs at build time. We keep $n$ small and use the diagonal phylo_indep cell so it converges in a few seconds.

set.seed(5640)
n_sp2 <- 25
tree2 <- ape::rcoal(n_sp2)
tree2$tip.label <- paste0("sp", seq_len(n_sp2))
C2 <- ape::vcv(tree2, corr = TRUE)
L2 <- t(chol(C2 + diag(1e-8, n_sp2)))

# Independent intercept and slope fields (rho = 0), sd_int = 0.6, sd_slope = 0.5
ab2 <- L2 %*% cbind(rnorm(n_sp2, sd = 0.6), rnorm(n_sp2, sd = 0.5))
rownames(ab2) <- tree2$tip.label

dfp <- merge(
  data.frame(species = factor(tree2$tip.label, levels = tree2$tip.label),
             x = rnorm(n_sp2)),
  data.frame(trait = factor(paste0("t", 1:3), levels = paste0("t", 1:3))),
  all = TRUE
)
# replicate to give the counts something to estimate from
dfp <- dfp[rep(seq_len(nrow(dfp)), each = 4), ]
dfp$x <- dfp$x + rnorm(nrow(dfp), sd = 0.2)
eta <- c(1.4, 1.2, 1.0)[as.integer(dfp$trait)] +
  ab2[as.character(dfp$species), 1] +
  ab2[as.character(dfp$species), 2] * dfp$x
dfp$value <- rpois(nrow(dfp), exp(eta))

fit_pois <- gllvmTMB(
  value ~ 0 + trait + phylo_indep(1 + x | species),
  data = dfp, trait = "trait", phylo_tree = tree2, unit = "species",
  family = poisson(link = "log")
)

c(convergence = fit_pois$opt$convergence,
  pd_hessian  = isTRUE(fit_pois$fit_health$pd_hessian))
#> convergence  pd_hessian 
#>           0           1
sd_b_pois <- as.numeric(fit_pois$report$sd_b)
data.frame(
  component = c("sd_intercept", "sd_slope"),
  true      = c(0.6, 0.5),
  estimated = signif(sd_b_pois, 3)
)
#>      component true estimated
#> 1 sd_intercept  0.6     0.451
#> 2     sd_slope  0.5     0.475
#> 3 sd_intercept  0.6     0.408
#> 4     sd_slope  0.5     0.456
#> 5 sd_intercept  0.6     0.368
#> 6     sd_slope  0.5     0.333

The slope variance recovers under the count family. The structured-slope grid is validated across nbinom1, nbinom2, Gamma, and Beta under this same diagonal (indep) route; binomial and ordinal_probit are validated under the latent and dep modes but remain structural-contract-only under indep – recovery there is not yet established. The response-families article is the family catalogue these slope terms compose with.

The correlated form: `phylo_dep(1 + x | species)`

The diagonal phylo_indep cell above pins the intercept-slope correlation to zero. The full unstructured phylo_dep cell estimates it. This is the canonical evolutionary reaction-norm question: do species that start high (large intercept) also respond more steeply (large slope)? The estimand is the $2T \times 2T$ covariance $\boldsymbol\Sigma_b$ over the per-trait (intercept, slope) columns, and its intercept-slope blocks carry that correlation.

This cell is heavy at the validation $n$ ( $n_{\text{species}}$ in the low hundreds for the count families). The chunk below shows the exact shipped syntax for both shapes and the extractor call, with eval = FALSE; the recovery numbers quoted afterward come from the validation tests, not from a build-time fit.

# Long format: explicit stacked-trait grammar
fit_dep_long <- gllvmTMB(
  value ~ 0 + trait + phylo_dep(1 + x | species),
  data = long_counts, trait = "trait", phylo_tree = tree, unit = "species",
  family = poisson(link = "log")
)

# Wide format: traits() LHS shorthand, same model
fit_dep_wide <- gllvmTMB(
  traits(t1, t2) ~ 1 + phylo_dep(1 + x | species),
  data = wide_counts, phylo_tree = tree, unit = "species",
  family = poisson(link = "log")
)

# Check before interpreting
check_gllvmTMB(fit_dep_long)

# The full augmented (intercept, slope) covariance, interleaved per trait:
#   intercept.t1, slope.t1, intercept.t2, slope.t2
Sigma_aug <- extract_Sigma(fit_dep_long, level = "phy")$Sigma
round(Sigma_aug, 3)

# Slope variances sit on the interleaved diagonal positions 2, 4, ...;
# the intercept-slope correlation is the (1,2) / (3,4) block off-diagonal.
slope_var <- diag(Sigma_aug)[c(2, 4)]

In the validation cell (Poisson, phylo_dep, $n_{\text{species}} = 150$ , 2 traits, 10 reps), the engine converges with a positive-definite Hessian and recovers the per-trait slope variances inside the mean-dependent-family band (within a factor of about 4 of truth – the honest tolerance for a count family at this fixture). The extract_Sigma(level = "phy")$Sigma matrix carries interleaved dimnames intercept.<t>, slope.<t>, so the intercept-slope correlation reads directly off the $2 \times 2$ block for each trait.

Why the band, not a point claim? Mean-dependent families (Poisson, negative binomial, Gamma, Beta) carry an implicit, mean-varying latent-scale residual, so the recoverable precision on a slope variance is looser than for a Gaussian or fixed-scale family. The validation suite uses a 3x-4x variance band for these families and a tighter (2.5x-3x) band for the fixed-scale families (binomial, ordinal probit). These are recovery tolerances, not calibrated confidence intervals – see the CI note below.

The spatial twin: `spatial_indep(1 + x | coords)`

When the structured units are locations rather than species, the same random-regression idea runs over a spatial field. The slope on $x$ varies smoothly across space, modelled as a second SPDE field alongside the intercept field. The spatial_indep form treats the two fields as independent; spatial_dep estimates their cross-covariance.

This needs a mesh (make_mesh()), which makes it heavier than the phylo cells, so it is shown eval = FALSE. The marginal field standard deviations come back through extract_Sigma(level = "spatial").

library(fmesher)

# A mesh over the site coordinates
mesh <- make_mesh(sites, c("lon", "lat"), cutoff = 0.1)

# Long format
fit_sp_long <- gllvmTMB(
  value ~ 0 + trait + spatial_indep(1 + x | coords),
  data = sites_long, trait = "trait", mesh = mesh,
  family = poisson(link = "log")
)

# Wide format
fit_sp_wide <- gllvmTMB(
  traits(t1, t2, t3) ~ 1 + spatial_indep(1 + x | coords),
  data = sites_wide, mesh = mesh,
  family = poisson(link = "log")
)

check_gllvmTMB(fit_sp_long)

# Marginal field SDs for the intercept and slope surfaces (2x2 cross-field)
S_sp <- extract_Sigma(fit_sp_long, level = "spatial")$Sigma
sqrt(diag(S_sp))   # intercept-field SD, slope-field SD

In the structural cell (Poisson, spatial_indep, 150 sites), the fit routes to the correct per-trait engine (2T augmented columns, 3T free block-diagonal parameters, correctly sized sd_spde_b) – this is wiring evidence, not recovery evidence. Convergence, Hessian positive-definiteness, and BLUP correlation with simulated truth are not yet measured for this cell (the cited test’s own header records full non-Gaussian recovery as ‘a follow-up’). The Multivariate spatial models article develops the spatial machinery (mesh construction, range and field-SD interpretation) in full.

Where to simulate from

For a Gaussian DGP without phylogeny or space, the shipped simulators give you a faithful stacked-trait long frame to attach a slope to. simulate_unit_trait() builds the generic (unit, observation, trait) cube and simulate_site_trait() the domain-specific (site, species, trait) cube with phylogenetic and spatial structure baked in. For the slope term specifically you add the per-unit slope draw and the x covariate by hand, exactly as the worked chunks above do – the simulators give you the intercept-level scaffold, and the reaction-norm slope is the piece you layer on top.

The full grid at a glance

Every cell below now fits; recovery within its honest band is established for the named validated cells, while binomial and ordinal traits under the diagonal mode remain structural-only. “Indep” pins the intercept-slope (or cross-field) correlation to zero; “latent” is the reduced-rank block; “dep” is the full unstructured covariance.

Source	`indep`	`latent`	`dep`
`phylo_*(1 + x`	species)	validated	validated
`spatial_*(1 + x`	coords)	validated	validated

Families covered for each cell: gaussian, poisson, nbinom1, nbinom2, Gamma, Beta, binomial (incl. multi-trial), and ordinal_probit. The deprecated single-variance slope keywords are not used in new public examples. indep: validated except binomial and ordinal_probit, which remain structural-contract-only under indep with recovery not yet established.

This table covers one random-slope covariate. Gaussian phylo_dep(1 + x1 + x2 | species) has its own recovery evidence, but the non-Gaussian s >= 2 path remains unavailable because its recovery evidence is not yet sufficient.

A note on confidence intervals

Everything here reports point estimates and recovery bands, not calibrated confidence intervals. The coverage study attained nominal coverage only for narrow Gaussian cases; for the non-Gaussian structured slopes it did not. Treat the recovery bands above as evidence that the point estimate lands near truth, not as a coverage guarantee. For inference on a slope variance, prefer external validation (a Bayesian refit in MCMCglmm or brms) until coverage is established. See profile-likelihood-ci for what the interval tools currently do and do not guarantee.