Bivariate non-Gaussian models: choose a joint or staged association
Source:vignettes/bivariate-nongaussian.Rmd
bivariate-nongaussian.RmdWith two non-Gaussian outcomes, first decide what kind of association answers the scientific question. A direct joint likelihood estimates both outcomes and their residual association together. A frozen-margin association first fits the two outcomes separately and then estimates a latent-normal association while holding the fitted margins fixed. These are different estimands.
This article helps you choose a route for two outcomes measured on
the same rows. The mixed-outcome route is a beta interface with
alpha-scale Godambe standard errors and Wald intervals for every
admitted pair class. The Bernoulli x ordinary-NB2 intercept route has
coverage evidence and is inference-ready with caveats; the other
admitted routes are interval-feasible and warn that coverage remains
uncalibrated. Random effects and a generic family-pair interface remain
outside this staged route. The exact direct lognormal route also has
callable uncertainty methods for a constant rho12, with a
retained fixed-effect calibration ladder for its tested domain; do not
treat a finite interval as a validated general non-Gaussian claim. For a
released direct Gaussian joint model, use biv_gaussian()
and Changing residual coupling with
rho12.
Start by choosing the association scale
| Your paired outcomes | Route | What the association means |
|---|---|---|
| Two Gaussian traits | Direct biv_gaussian() model |
Gaussian residual rho12
|
| Two positive traits plausibly lognormal | Post-0.6 direct biv_lognormal() model |
Log-response residual rho12
|
| Two heavy-tailed real-valued traits | Post-0.6 direct biv_student() model |
Shared-Student-t residual/scatter rho12
|
| A reviewed mixed or discrete pair |
biv_associate() staged route |
Latent-normal copula eta after frozen margins; an
intercept-bearing fixed-effect association formula is interval-feasible
for Bernoulli x NB2 |
Neither rho12 nor eta is automatically the
Pearson correlation of the two raw response columns. The scale is part
of the scientific interpretation, not just a computational detail.
Direct joint model: two positive traits
For two positive lognormal outcomes, the model is joint on the log-response scale:
Here rho12 is the correlation of the two log-response
residuals. A positive value means that units above their fitted
log-response expectation on one trait also tend to be above it on the
other. It is not automatically the correlation of the two traits on
their original scale.
Real example: penguin size after measured composition differences
The openly licensed palmerpenguins data set records body
size in three Palmer Archipelago penguin species. It is distributed
under CC0 and collates data from the Palmer Station LTER study (data policy;
Gorman et
al. 2014).
Our question is conditional: after allowing both margins to vary by species, sex, and centred year, do birds with longer-than-expected flippers also have greater-than-expected body mass? This is not a causal claim that flippers determine mass, nor a replacement for a raw scatterplot.
penguins <- subset(
palmerpenguins::penguins,
complete.cases(species, sex, year, flipper_length_mm, body_mass_g) &
flipper_length_mm > 0 & body_mass_g > 0
)
penguins$year_c <- penguins$year - mean(penguins$year)
nrow(penguins)
#> [1] 333
plot(
log(body_mass_g) ~ log(flipper_length_mm), data = penguins,
pch = 16, col = grDevices::adjustcolor("#087f8c", alpha.f = 0.45),
xlab = "log flipper length (mm)", ylab = "log body mass (g)"
)
abline(stats::lm(log(body_mass_g) ~ log(flipper_length_mm), data = penguins),
col = "#14344a", lwd = 2
)
Raw log-scale flipper length and body mass for analysed complete rows. This descriptive relationship is not the fitted residual association.
fit_log <- drmTMB::drmTMB(
drmTMB::bf(mu1 = flipper_length_mm ~ species + sex + year_c,
mu2 = body_mass_g ~ species + sex + year_c,
sigma1 = ~ 1, sigma2 = ~ 1, rho12 = ~ 1),
family = drmTMB::biv_lognormal(), data = penguins
)
drmTMB::rho12(fit_log)[1]
#> [1] 0.3583776
drmTMB::check_drm(fit_log)
#> check status
#> optimizer_convergence ok
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> hessian_conditioning ok
#> standard_errors_finite ok
#> standard_errors_inflated ok
#> dropped_rows ok
#> positive_scale ok
#> rho12_boundary ok
#> fixed_effect_design_size ok
#> value
#> 0
#> iterations=13; function=30; gradient=13
#> 3403.
#> <NA>
#> max=0.0000000001547; component=beta_mu1[1]
#> ok
#> TRUE
#> min_eig=298.3; cond=2683.
#> range=[0.001821,0.05480]
#> n_inflated=0; max_se=0.05480; median_se=0.007649
#> nobs=333; dropped=0
#> min=0.02695
#> 0.3584
#> total_mb=0.1401; max_cols=5; largest=mu1; largest_class=matrix; largest_density=0.6132
#> message
#> nlminb convergence code is 0.
#> Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#> Objective and log-likelihood are finite.
#> The log(sigma) clamp is not active at the optimum.
#> Maximum absolute fixed gradient is <= 0.001; largest component is beta_mu1[1].
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> Minimum eigenvalue and condition number of TMB's sdreport() fixed-effect covariance (sdr$cov.fixed), inverted. These are a genuinely different read of the fit's conditioning than TMB's internal pdHess flag -- comparable across fits, not claimed to be numerically identical to any raw TMB gradient or Hessian quantity. This fit's Hessian conditioning is within the requested threshold.
#> All fixed-effect standard errors are finite.
#> No fixed-effect standard error is inflated relative to the others.
#> No rows were dropped by model-frame or known-covariance filtering.
#> All fitted scale values are finite and positive.
#> All fitted residual correlations have absolute value <= 0.98.
#> Dense fixed-effect design matrices are modest for this fit.
# All three target the log-residual rho12, not a raw-scale correlation.
ci_wald <- confint(fit_log, parm = "rho12", method = "wald")
ci_profile <- confint(
fit_log, parm = "rho12", method = "profile", profile_engine = "endpoint"
)
ci_bootstrap <- confint(
fit_log, parm = "rho12", method = "bootstrap", R = 99, seed = 20260724
)
interval_row <- function(x, label) {
data.frame(
label = label,
estimate = unname(drmTMB::rho12(fit_log)[1]),
lower = x$lower,
upper = x$upper,
status = x$conf.status,
note = x$profile.message,
stringsAsFactors = FALSE
)
}
intervals <- rbind(
interval_row(ci_wald, "Wald"),
interval_row(ci_profile, "Profile"),
interval_row(
ci_bootstrap,
sprintf(
"Bootstrap (%d/%d retained)",
ci_bootstrap$bootstrap.n,
ci_bootstrap$bootstrap.n + ci_bootstrap$bootstrap.failed
)
)
)
bootstrap_total <- ci_bootstrap$bootstrap.n + ci_bootstrap$bootstrap.failed
cat(sprintf(
"Bootstrap diagnostic: %d/%d full refits retained. Failed refits remain a diagnostic, not hidden precision.\n",
ci_bootstrap$bootstrap.n, bootstrap_total
))
#> Bootstrap diagnostic: 99/99 full refits retained. Failed refits remain a diagnostic, not hidden precision.
Reported 95% confidence intervals for the fitted direct log-residual association in the penguin model. Each eye spans the reported interval and is centred on its estimate; the taper is a visual interval cue, not a likelihood, sampling density, or posterior distribution. The retained bootstrap count is a diagnostic, not hidden precision.
fitted(fit_log) returns arithmetic marginal means on the
original response scale. The mu1 and mu2
predictors remain log-scale locations. These direct models estimate the
two margins and rho12 together; they are not staged
frozen-margin fits. Wald is a fast link-scale approximation; profile
likelihood checks the exact likelihood geometry; the bootstrap refits
both margins and rho12. Inspect
check_drm(fit_log) and bootstrap diagnostics, especially
near
.
In this data set, do not hide failed bootstrap refits: the printed
bootstrap.n and bootstrap.failed fields are
part of the result. The external data example illustrates interpretation
and diagnostics; it does not itself establish coverage calibration or
prove lognormal adequacy. Inspect residual diagnostics and compare
scientifically defensible alternatives before reporting a substantive
association.
If a fit, covariance, or interval is withheld, follow Errors, warnings, and convergence rather than repairing, clipping, or hiding the diagnostic result.
Direct joint model: two heavy-tailed traits
For two real-valued traits with occasional extreme observations, the
bivariate Student-t model uses a shared degrees-of-freedom parameter,
nu:
fit_t <- drmTMB(
bf(mu1 = activity ~ habitat, mu2 = boldness ~ habitat,
sigma1 = ~ 1, sigma2 = ~ 1, nu = ~ 1, rho12 = ~ 1),
family = biv_student(), data = behaviour
)
rho12(fit_t)sigma1 and sigma2 are Student-t scales, not
necessarily marginal standard deviations. The model has one shared
nu, which controls tail heaviness. At a finite
nu, rho12 = 0 means zero residual correlation,
but not necessarily complete independence because the two outcomes share
the same row-level heavy-tail mechanism.
Staged association: different kinds of outcome
If one response is binary and the other is a count, there is no
shared direct Gaussian residual scale. The frozen-margin route instead
fits a binary probability model and a count model, then estimates
latent-normal eta from their fitted distributions.
paired_data <- na.omit(data.frame(
bred, offspring, habitat, habitat_score, season
))
assoc <- biv_associate(
bf(mu = bred ~ habitat),
bf(mu = offspring ~ habitat, sigma = ~ season),
family = list(binomial(), nbinom2()),
data = paired_data,
association = ~ 1
)
association(assoc)
vcov(assoc)
confint(assoc)
confint(assoc, type = "eta")This is one convenient R call, but not one jointly fitted model.
biv_associate() fits the two margins, freezes them, and
estimates only eta. The reported value is a latent-normal
copula correlation conditional on those fitted margins. It is not
rho12, a logit coefficient, an odds ratio, or an
observed-scale correlation. For every admitted association route,
vcov() and plain confint() report two-stage
Godambe uncertainty for the unbounded association-link coefficients
alpha whenever the fit-specific covariance diagnostics
pass. confint(assoc, type = "eta") transforms the constant-
association interval to bounded eta, while predict()
supplies delta-method eta standard errors and pointwise transformed
intervals.
For a binary response, the model observes whether a hidden tendency crossed a threshold. It does not observe an ordinary numeric residual. Read Frozen-margin association for mixed outcome pairs for the binary-threshold explanation, the reviewed pair classes, and the rule for withheld estimates.
For the literal-Bernoulli x ordinary-NB2 beta route, an
intercept-bearing fixed-effect formula such as
association = ~ x, ~ habitat, or
~ x * habitat can model a row-specific latent-normal
eta_i after both margins have been fitted and frozen.
Multiple predictors, factors, interactions, and explicit transformations
use the ordinary fixed-effect model matrix. See Association between mixed outcome pairs for
the exact syntax, output, and limits.
What is currently outside these routes
The direct same-family first slices keep sigma1,
sigma2, rho12, and nu constant
across rows. The staged route usually keeps eta constant;
its beta exception is an intercept-bearing fixed-effect association
formula for literal-Bernoulli x ordinary-NB2. Its reviewed fixed-effect
margins may include their explicitly supported covariates before they
are frozen. Neither route currently adds random effects or missing-data
support, and the staged formula rejects offsets, missing predictors,
aliased columns, dot expansion, and random effects. Other pair classes
remain intercept-only. Neither route is a general non-Gaussian bivariate
claim. The staged route’s alpha-scale standard errors and Wald intervals
are interval-feasible for every admitted pair class and for the
Bernoulli x ordinary-NB2 association-regression formula. Routes without
coverage calibration emit an experimental-interval warning. Eta
uncertainty inherits that warning and evidence tier; simultaneous eta
bands and profiles remain unavailable. The direct lognormal
constant-rho12 profile, Wald, and bootstrap methods have
direct calibration evidence only for the fixed-effect DGP in the Arc 6
coverage artifact; they do not validate the staged route or a general
non-Gaussian association claim. Choose the route whose stated
association scale matches the biological question, and treat every
boundary as part of the model definition.