Cross-family dependence: statistical foundations

Status — Experimental (engine), paper-section reference

The cross-family route is Experimental and not release-ready: the evidence is a set of simulated fixtures covering five family pairs (test/test_mixed_family.jl) and there is no interval-coverage study. The maths below is what the code evaluates; the release status is what limits its use. See Evidence & limits. This page is the methods reference for the cross-family bivariate model implemented by fit_mixed_family in src/mixed_family.jl, with the link-scale standardisation in src/link_residual.jl. It documents the shared-latent construction, the Gauss–Hermite marginal likelihood, the link-residual variances , the latent-scale correlation and its identifiability, the three confidence intervals, and the exact reduction to the bivariate rho12 Gaussian model. It is written to paper-section standard; every formula here is the one the code evaluates.

1. Motivation and scope

The bivariate Gaussian model (bf(mu1=…, mu2=…, sigma1=…, sigma2=…, rho12=…)) couples two responses through a residual correlation on a common identity-link, real-valued scale. That construction is only available when both responses are Gaussian: is the correlation of two residuals that live on the same scale, so it is undefined the moment one response is a count, a proportion, or a strictly positive quantity.

The cross-family model removes that restriction. It couples two responses that may come from different families — Gaussian × Poisson, Gamma × Beta, negative-binomial × Binomial, and so on — through a single shared per-observation latent variable, and reports the induced dependence as a correlation on the latent (link) scale, . When both families happen to be Gaussian the model reduces exactly to the rho12 bivariate Gaussian (Section 7), so it is a strict generalisation rather than a separate object.

Throughout, the two axes are indexed  ;   indexes the paired observations. We follow the package convention of naming the scale (never ). The residual-correlation symbol is reserved for the same-scale Gaussian construction; the latent-scale cross-family correlation is written without a subscript.

2. The shared-latent construction

Each observation carries one scalar standard-normal latent variable , shared by both responses. Conditional on the two responses are independent, each drawn from its own family with a family-specific linear predictor:

Here is the fixed-effect linear predictor for axis (design matrices , coefficients ), and is a scalar factor loading that controls how strongly the shared latent feeds into axis . All dependence between and flows through : it is a one-factor model with a single common factor per observation. This is the canonical generalised-linear latent-variable (GLLVM) device — a shared latent on the linear predictor inducing cross-response dependence — specialised to two responses and a single factor (Skrondal & Rabe-Hesketh 2004; Niku et al. 2019).

Each axis uses its family's canonical link, so is on the link scale: identity for Gaussian, for Poisson / negative-binomial / Gamma, for Binomial / Beta. The conditional mean on the response scale is (identity, , or logistic).

Dispersion sub-model

Dispersion-carrying families also estimate a per-observation log-native dispersion through its own linear predictor,

with design matrix (default a single intercept column, which collapses to a constant). The natural dispersion plays the family-specific role given in the table below. Poisson and Binomial are dispersionless — they are pinned by alone and carry no .

Family	is	Conditional density at node
Gaussian	residual SD
Poisson	—
Binomial	—
NegBinomial2	size	,
Beta	SD , precision	,
Gamma	SD (CV), shape	,

The full parameter vector is

where is the log of the first loading (Section 5) and the blocks are present only for dispersion-carrying axes.

3. The Gauss–Hermite marginal likelihood

Because is one-dimensional, the marginal likelihood is a product of scalar integrals — one per observation — that we evaluate exactly (to quadrature order) rather than by Laplace. Marginalising the shared latent,

with    and the standard normal density. The conditional log-densities are the family terms tabulated in Section 2 (and coded in _mf_obs_ll).

Change of variables to Gauss–Hermite form

Gauss–Hermite quadrature approximates      with nodes/weights (Golub–Welsch). Writing   and substituting    (so    and  ) turns each observation's integral into

The prefactor is exactly the   correction subtracted per observation in the code. The negative marginal log-likelihood is therefore

evaluated with the numerically stable (subtract the per-observation maximum, exponentiate, sum, re-log). The default node count is  .

Why quadrature, not Laplace

For this model the integrand is genuinely one-dimensional, so a -node Gauss–Hermite rule is both cheap ( per likelihood evaluation) and far more accurate than a single mode-plus-curvature Laplace match — it captures skewed and heavy-tailed conditional posteriors that Laplace cannot see (the shape/dispersion bias documented in Laplace vs variational marginals). The summand is a smooth, deterministic function of (the nodes are constants), so the objective is differentiable by forward-mode automatic differentiation and is optimised with L-BFGS. A finite-penalty guard maps any non-finite objective value (e.g. an overflow from a large Poisson loading mid-line-search) to a large finite number of the same dual type, so the optimiser backtracks rather than aborting; the well-conditioned path — hence every converged estimate and interval — is untouched.

4. Link-residual standardisation

The loadings and the marginal log-likelihood are well defined, but and are not directly comparable across families: a loading of means something different on a scale (Poisson) than on a scale (Binomial) than on an identity scale (Gaussian). To report a single interpretable dependence we standardise to the latent (link) scale, following Nakagawa & Schielzeth (2010).

The key quantity is the family's own observation-level variance on its link scale, the distribution-specific variance . On the link scale the total latent variance of axis decomposes as

so the latent-scale variance of axis is   and the cross-axis latent covariance is (the only shared term). is computed by link_residual and is reporting-only: it never enters the fit objective.

Per-family derivations

The link-residual map is exactly Nakagawa & Schielzeth's (2010, Table 2; see also Nakagawa, Johnson & Schielzeth 2017) distribution-specific variance for each family's canonical link, evaluated at a representative fitted mean  .

Family	Link		Origin
Gaussian	identity		residual variance enters the latent scale directly
Poisson			delta-method / lognormal variance of
Binomial			variance of a standard logistic; distribution-free
NegBinomial2			lognormal-approx. dispersion term, = trigamma
Gamma			variance of under
Beta			variance of ,

A few derivations worth spelling out, because they are the ones a reader will want to check:

• Gaussian,  . With an identity link the residual is on the link scale, so its variance enters   directly. This differs from the GLLVM / gllvmTMB convention, which reports   for a Gaussian trait because there the residual variance lives in a separate covariance block. In this shared-latent parameterisation there is no separate : setting   would force   for Gaussian × Gaussian.   is the value verified by the exact Gaussian × Gaussian ≡ rho12 parity (Section 7).

• Poisson,   . Under the link a Poisson mean with   has, to the lognormal / delta-method order used throughout Nakagawa & Schielzeth, an observation-level variance on the log scale of   , refined to   (the exact lognormal-matching expression that reproduces for large while staying finite and positive for small ). It is evaluated at the representative mean .

• Binomial,  . On the link the observation variance is taken as the variance of the standard logistic distribution, . This is the standard latent-threshold value (the variance of the logistic latent in a logit model) and is distribution-free — it does not depend on or any dispersion.

• Gamma,  . With a link and shape  , for    has variance (the trigamma function   is exactly the variance of of a Gamma deviate). Here is the squared coefficient of variation, the package's Gamma sigma-slot convention.

• NegBinomial2,  . The lognormal-approximation distribution-specific variance contributed by the size/dispersion on the link, matching gllvmTMB's nbinom2 entry.

• Beta,    . For    the variance of     is the sum of the trigamma terms of the two Beta shape parameters (precision  ), evaluated at .

5. The latent-scale correlation and identifiability

With the latent variances   and the shared covariance in hand, the latent-scale correlation is the natural ratio:

This is the cross-family generalisation of a residual correlation: it is the correlation of the two latent linear predictors after each axis is standardised by its own total link-scale variance.    by construction (Cauchy–Schwarz, since  ), with the sign carried by .

Identifiability

Three identifiability facts govern how is parameterised and what is estimable: 2. Sign of the shared factor (  ). The likelihood is invariant under flipping the latent,   , which sends     and leaves both distributions unchanged. This is the familiar factor-loading sign indeterminacy. We break it by constraining the first loading positive,

so is unconstrained and cannot change sign. With   pinned, the sign of is the sign of , which is now identified.

3. Gaussian × non-Gaussian: fully identified. When at least one axis is non-Gaussian, that axis has no free residual variance competing with the loading on the same scale (its observation variance is fixed by the mean / dispersion, not by ). All parameters — , both loadings, the dispersions, and hence — are then separately identified.

4. Gaussian × Gaussian: only and identified. When both axes are Gaussian the model has a flat ridge in the individual loadings: the marginal is exactly bivariate normal (Section 7), and a bivariate normal is fully described by its means, the two marginal variances  , and the covariance . The split of   between the loading and the residual is not separately identified (any on the ridge gives the same marginal variance). What is identified is precisely the marginal-covariance summary: the log-likelihood and the correlation . This is not a defect — is the only cross-axis quantity the data can speak to in the Gaussian × Gaussian case, and it is exactly the bivariate rho12 (Section 7). Reported individual loadings on this ridge should not be over-interpreted; the log-likelihood and are the trustworthy outputs.

6. Confidence intervals for {#6.-Confidence-intervals-for-\rho}

Three intervals for are provided, in increasing order of cost and robustness. All operate on Fisher's   where a transformation is used, because the sampling distribution of a correlation is far closer to normal on the scale.

6.1 Fisher- Wald (delta-method)

The default interval. Let be the maximum-likelihood estimate and   the inverse observed-information matrix, where    is the Hessian of the negative log-likelihood at (computed by automatic differentiation). Define   and its gradient    (also by automatic differentiation through and ). The delta-method variance of is

and the   interval is formed on the scale and mapped back by :

with the standard normal quantile. The back-transform guarantees the endpoints stay in  . If is not invertible or is non-positive, the interval is returned as (NaN, NaN) rather than a spurious number — the same honesty discipline the engine applies whenever the observed information is singular.

6.2 Profile likelihood

The recommended interval when the Wald approximation is doubted (small , or near the boundary), better calibrated than Wald and cheaper than the bootstrap. Fix the correlation at a trial value and re-optimise all other parameters subject to that constraint, enforced with a quadratic penalty on the scale:

The profile deviance at is

and the profile-likelihood interval is the set   , i.e. the two where crosses the quantile at level  . The two endpoints are found by bisection on each side of (within ); a side that never reaches the threshold inside   is reported at that boundary.

6.3 Parametric bootstrap

The most robust and most expensive interval, useful as an external check on the two analytic intervals. At the fitted simulate replicate data sets from the model itself — draw a fresh shared latent   per observation, form   , and draw from the corresponding family at (the per-family samplers in _mf_rand) — refit the model to each replicate, and collect the bootstrap correlations . The interval is the percentile interval

the empirical and   quantiles of the bootstrap distribution. Non-converged replicates are discarded; the interval is returned only if at least   replicates succeed, otherwise (NaN, NaN).

7. Exact reduction to the rho12 bivariate model

When both axes are Gaussian the cross-family model is not an approximation of the bivariate Gaussian — it is algebraically identical to it, and equals the residual correlation .

Take    with residual SDs . Conditional on ,

and marginalising the Gaussian keeps the pair Gaussian (a linear combination of independent normals). The marginal moments are

so

This is exactly the bivariate Gaussian fitted by bf(mu1=…, mu2=…, sigma1=…, sigma2=…, rho12=…), whose residual correlation is

Substituting the Gaussian link-residual   (Section 4) into the latent-scale correlation gives

The marginal model is the same as the rho12 bivariate Gaussian model. A finite Gauss–Hermite rule is still a numerical quadrature approximation here: its polynomial-degree exactness does not make a Gaussian-density integrand exact. Thus fitted likelihood values need numerical agreement checks at a stated tolerance. The mathematical equivalence identifies the same and reported correlation: the loading split is the unidentified ridge, while   is the trustworthy cross-axis summary. This equivalence is used as a correctness anchor for the engine, with numerical parity checked separately.

8. Summary of the construction

Object	Expression	Where
Linear predictor		§2
Shared latent		§2
Dispersion sub-model		§2
Marginal likelihood	via -node Gauss–Hermite	§3
Link-residual variance	(Nakagawa & Schielzeth 2010, per family)	§4
Latent correlation		§5
Sign fix		§5
CIs	Fisher- Wald · profile likelihood · parametric bootstrap	§6
Gaussian × Gaussian	rho12 bivariate Gaussian,	§7

References

• Nakagawa, S., & Schielzeth, H. (2010). Repeatability for Gaussian and non-Gaussian data: a practical guide for biologists. Biological Reviews, 85(4), 935–956. — distribution-specific link-scale variances .

• Nakagawa, S., Johnson, P. C. D., & Schielzeth, H. (2017). The coefficient of determination and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. Journal of the Royal Society Interface, 14(134), 20170213. — the observation/link-scale variance decomposition used here.

• Skrondal, A., & Rabe-Hesketh, S. (2004). Generalized Latent Variable Modeling: Multilevel, Longitudinal, and Structural Equation Models. Chapman & Hall/CRC. — shared-latent (common-factor) construction for cross-response dependence.

• Niku, J., Hui, F. K. C., Taskinen, S., & Warton, D. I. (2019). gllvm: Fast analysis of multivariate abundance data with generalized linear latent variable models in R. Methods in Ecology and Evolution, 10(12), 2173–2182. — GLLVM latent-variable dependence and link-scale standardisation (the gllvmTMB link_residual_per_trait lineage this implementation mirrors).

• Golub, G. H., & Welsch, J. H. (1969). Calculation of Gauss quadrature rules. Mathematics of Computation, 23(106), 221–230. — the Gauss–Hermite nodes/weights.

See also

• Choosing response families — the families that can appear on either axis.

• Which scale are you modelling? — link vs response scale.

• Laplace vs variational marginals — why a 1-D quadrature beats a Laplace match for the shape/dispersion parameters.

• fit_mixed_family (in src/mixed_family.jl) — the fitting function this page documents.