Cross-family bivariate dependence

Status — Experimental

A first slice of bivariate modelling for two responses from different families (e.g. Gaussian × Poisson), via a shared per-observation latent. The matrix-level entry point is DRM.fit_mixed_family (not exported); the bf(...) front end is wired as drm(bf(...), (Gaussian(), Poisson()); data = …) (src/mixed_family.jl:440). Gaussian, Poisson, Binomial, NegBinomial2, Beta and Gamma axes are supported (src/mixed_family.jl:32-66); the route takes family instances, e.g. (Gaussian(), Poisson()). The dependence is reported on the link/latent scale; the residual-correlation rho12 model is the Gaussian × Gaussian special case.

The residual-correlation model (Changing residual coupling with rho12) couples two Gaussian responses through a residual correlation ρ12. But often the two responses are not from the same family — a continuous trait and a count, a count and a proportion. There is no single residual covariance matrix to write down, because the two responses live on different scales with different mean–variance relationships. DRM.fit_mixed_family handles that case.

The model

Each observation i carries a single shared scalar latent uᵢ ~ N(0, 1). The two linear predictors load on it with their own loadings λ₁, λ₂:

and each response is then drawn from its own family, conditional on its own linear predictor:

The shared uᵢ is what induces the dependence: a large draw pushes both η₁ᵢ and η₂ᵢ in directions set by the signs and sizes of λ₁ and λ₂. With λ₁ λ₂ > 0 the two responses move together; with λ₁ λ₂ < 0, they oppose.

Because uᵢ is a single scalar per observation, the marginal likelihood is a one-dimensional integral, evaluated by Gauss–Hermite quadrature with K nodes:

The quadrature makes the marginal smooth and ForwardDiff-friendly, so the fit is a plain L-BFGS optimisation with forward-mode autodiff. For identifiability, λ₁ is constrained positive (λ₁ = exp(·)) to remove the u → -u sign flip.

Gaussian × Gaussian ≡ the rho12 model

When both axes are Gaussian the marginal is exactly bivariate normal, so the marginal log-likelihood and the reported correlation reduce to the residual-correlation rho12 model. In that case the two loadings are not separately identified (a flat λ₁ ridge), but the marginal covariance — hence the log-likelihood and ρ — is. For Gaussian × non-Gaussian (no free residual variance on the non-Gaussian axis), all parameters are identified.

The latent-scale correlation

The loadings λ₁, λ₂ are not directly comparable across families, because each family contributes its own observation-level noise on its link scale. To put the dependence on a single, interpretable scale we standardise it the way Nakagawa & Schielzeth (2010) standardise variance components — adding each family's link-scale residual variance v_k to the latent variance it already carries:

The numerator λ₁ λ₂ is the cross-axis covariance contributed by the shared latent; each denominator term λ_k² + v_k is the total latent-scale variance of axis k. v_k comes from DRM.link_residual:

• Gaussian (identity link): v = σ², the fitted residual variance.

• Poisson (log link): v = log(1 + 1/μ̄), evaluated at a representative fitted mean μ̄.

• Binomial (logit link): v = π²/3 (distribution-free).

• Beta (logit link): v = trigamma(μ̄ φ) + trigamma((1-μ̄) φ), φ the precision.

• Gamma (log link): v = trigamma(1/φ).

• NegBinomial2 (log link): v = trigamma(θ), θ the size.

Because ρ is a ratio of the same latent variance on top and bottom, it is rotation-invariant and lies in (-1, 1). It is a reporting quantity: it is computed from the fitted parameters and never enters the objective.

Why v = σ² for the Gaussian axis here

In this shared-latent parameterisation the Gaussian residual lives in v directly (there is no separate Ψ as in GLLVM/gllvmTMB, which report v = 0 for Gaussian because the residual sits elsewhere). Using v = σ² is what makes the Gaussian × Gaussian case collapse onto the rho12 model; using 0 would force ρ = 1.

Confidence intervals for ρ

fit_mixed_family can return three intervals for ρ, each on the correlation scale:

Field	Method	When to use
rho_ci_wald	Fisher-z (delta method on atanh ρ)	Always computed (confint = true); cheapest.
rho_ci_profile	Profile likelihood (profile = true)	Recommended — best calibrated near the boundary.
rho_ci_boot	Parametric bootstrap (B = … refits)	Most robust; most expensive.

• Fisher-z Wald applies the delta method to atanh(ρ(θ)) using the observed information (the Hessian of the marginal negative log-likelihood), then maps the symmetric atanh-scale interval back through tanh. The atanh transform keeps the interval inside (-1, 1). Returns (NaN, NaN) if the Hessian is not invertible.

• Profile likelihood (recommended) re-optimises the model under a soft constraint that fixes ρ(θ) at a trial value on the atanh scale, and bisects for the values where the deviance rises by the χ²(1) quantile. It is better calibrated than Wald when ρ is near ±1, and cheaper than the bootstrap.

• Parametric bootstrap resamples the shared latent and the per-family draws at the fitted parameters, refits each replicate, and takes the percentile interval of the resulting ρ̂. It makes the fewest distributional assumptions but costs B extra fits.

Worked example: Gaussian × Poisson

We simulate a continuous response (Gaussian) and a count (Poisson) that share a latent, recover the structure with fit_mixed_family, and read off the latent-scale correlation with its profile-likelihood interval.

using DRM, Random, Statistics
using DRM: fit_mixed_family
Random.seed!(2024)

n = 800
u = randn(n)                       # shared per-observation latent, u ~ N(0,1)
x = randn(n)                       # one shared covariate, on both axes
X1 = hcat(ones(n), x)
X2 = hcat(ones(n), x)

# True parameters.
β1, β2 = [1.0, 0.4], [0.2, 0.3]    # fixed effects (Gaussian, Poisson)
λ1, λ2 = 0.8, 0.5                  # latent loadings

η1 = X1 * β1 .+ λ1 .* u
η2 = X2 * β2 .+ λ2 .* u
y1 = η1 .+ randn(n)                                            # Gaussian, σ = 1
y2 = [Float64(rand(DRM.Distributions.Poisson(exp(η2[i])))) for i in 1:n]  # Poisson

fit = fit_mixed_family(; y1 = y1, X1 = X1, fam1 = Gaussian(),
                         y2 = y2, X2 = X2, fam2 = Poisson(),
                         K = 32, profile = true)

fit.converged        # did L-BFGS converge?

true

The fitted loadings and the link-scale residual variances feed the latent-scale correlation:

(λ1 = round(fit.λ1, digits = 3), λ2 = round(fit.λ2, digits = 3),
 v1 = round(fit.v1, digits = 3), v2 = round(fit.v2, digits = 3))

(λ1 = 0.791, λ2 = 0.487, v1 = 0.84, v2 = 0.564)

v1 is the fitted Gaussian residual variance (the true value is σ² = 1); v2 = log(1 + 1/μ̄) is the Poisson axis' link-scale variance at its representative mean. The latent-scale correlation and its profile-likelihood interval are:

round(fit.rho_latent, digits = 3)

0.355

round.(fit.rho_ci_profile, digits = 3)        # recommended interval for ρ

(0.284, 0.422)

The point estimate matches the value implied by the simulation truth, ρ = λ₁λ₂ / sqrt((λ₁²+v₁)(λ₂²+v₂)) ≈ 0.36, and the profile interval covers it. The cheaper Fisher-z Wald interval is also available:

round.(fit.rho_ci_wald, digits = 3)           # Fisher-z delta-method interval

(0.287, 0.42)

For the bootstrap interval, pass B (number of refits); it is omitted here to keep the page fast to build:

fit_boot = fit_mixed_family(; y1, X1, fam1 = Gaussian(),
                              y2, X2, fam2 = Poisson(),
                              K = 32, B = 500)
fit_boot.rho_ci_boot       # percentile interval from 500 parametric-bootstrap refits

Returned fields

fit_mixed_family returns a NamedTuple. The dependence-related fields are:

Field	Meaning
rho_latent	Latent-scale correlation ρ (the headline estimate).
rho_ci_wald	Fisher-z Wald interval (lo, hi) for ρ (or (NaN, NaN)).
rho_ci_profile	Profile-likelihood interval (recommended); (NaN, NaN) unless profile = true.
rho_ci_boot	Parametric-bootstrap interval; (NaN, NaN) unless B > 0.

The fit also returns the fixed effects β1/β2, the loadings λ1/λ2, the Gaussian residual SDs σ1/σ2 (NaN on non-Gaussian axes), the link-scale variances v1/v2, the loglik, converged, and iterations.

See also

• Changing residual coupling with rho12 — the Gaussian × Gaussian residual-correlation model that this generalises.

• DRM.link_residual — the per-family link-scale variance v_k.

API

DRM.link_residual Function

link_residual(fam, μ̂; dispersion) -> Float64

Distribution-specific observation variance v of fam on its link (latent) scale.

• Gaussian → dispersion (residual variance σ²); identity link.

• Poisson → log(1 + 1/μ̂); log link (μ̂ a representative fitted mean).

• Binomial → π²/3; logit link (distribution-free).

• Beta → trigamma(μ̂ φ) + trigamma((1-μ̂) φ); logit link, dispersion = φ.

• Gamma → trigamma(1/φ); log link, dispersion = φ (the variance σ²).

• NegBinomial2 → trigamma(θ); log link, dispersion = θ (the size/dispersion).

Feeds the cross-family latent correlation ρ = λ1 λ2 / sqrt((λ1²+v1)(λ2²+v2)).

source

Post-fit accessors

DRM.mf_coef Function

mf_coef(fit) -> NamedTuple

Tidy table of the cross-family fit's point estimates as three equal-length vectors (; axis, term, estimate). Rows, in order:

• β1/β2: fixed-effect coefficients per axis, labelled b1[1], b1[2], ….

• bσ1/bσ2: dispersion sub-model coefficients (log-native scale), labelled bsig1[…] (omitted for dispersionless axes — Poisson/Binomial).

• lambda: the two latent loadings λ1/λ2 (axis = :shared).

• rho: the latent-scale correlation rho_latent (axis = :shared).

axis is :y1 / :y2 for the per-axis rows and :shared for the loadings and ρ. Estimates are on the model's native parameter scale (link scale for β, log-native for βσ).

source

DRM.mf_summary Function

mf_summary(fit; nobs = nothing, io = stdout)

Print a human-readable summary of a fit_mixed_family fit: the coefficient table from mf_coef, the latent correlation ρ with its available CIs (Wald / profile / bootstrap, whichever are finite), the log-likelihood, the AIC (and BIC when nobs is supplied), and the convergence flag. Returns fit invisibly. Pass nobs (the number of observation pairs) to include BIC.

source

DRM.mf_fitted Function

mf_fitted(fit, X1, X2) -> NamedTuple

Per-axis fitted means on the RESPONSE scale, evaluated at the latent u = 0 (the population/marginal-mode convention): μ_k = g_k⁻¹(X_k β_k) with the shared random effect set to zero, i.e. the fitted value for a "typical" individual whose latent draw is at the mean. Returns (; mu1, mu2), each a length-size(X_k, 1) vector. The inverse link g_k⁻¹ is the family's own (identity for Gaussian, exp for Poisson/NB2/Gamma, logistic for Binomial/Beta — see _mf_mean).

NOTE: this is the conditional mean at u = 0, NOT the marginal mean E[y_k] (which for a non-identity link differs from g_k⁻¹(Xβ) by a Jensen term in the shared latent). The u = 0 convention matches the link-scale linear predictor used everywhere else in the model.

source

mf_fitted(fit)

Fitted means for a cross-family fit produced by the formula route (drm(bf(...), (Fam1(), Fam2()); data = …)), which carries its own design matrices. A fit built from raw matrices via fit_mixed_family does not, so there they must still be supplied: mf_fitted(fit, X1, X2).

Same u = 0 convention as the three-argument method above.

source

DRM.mf_aic Function

mf_aic(fit) -> Float64

Akaike information criterion, -2·loglik + 2·k, where k = _mf_nparams(fit) is the number of free parameters. fit_mixed_family fits by ML, so this is directly comparable across mean/dispersion structures. Lower is better.

source

DRM.mf_bic Function

mf_bic(fit; nobs) -> Float64

Bayesian (Schwarz) information criterion, -2·loglik + k·log(nobs), with k = _mf_nparams(fit). The fit NamedTuple does not carry the sample size, so nobs — the number of observation PAIRS n — is a required keyword. Lower is better; comparable across structures (ML fit).

source

References

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