--- title: "Bivariate non-Gaussian models: choose a joint or staged association" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Bivariate non-Gaussian models: choose a joint or staged association} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") has_penguins <- requireNamespace("palmerpenguins", quietly = TRUE) ``` With 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](bivariate-coscale.html). ## 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: \[ (\log Y_{1i}, \log Y_{2i}) \sim N_2\left( (\mu_{1i}, \mu_{2i}), \begin{bmatrix} \sigma_1^2 & \rho_{12}\sigma_1\sigma_2 \\ \rho_{12}\sigma_1\sigma_2 & \sigma_2^2 \end{bmatrix} \right). \] 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](https://allisonhorst.github.io/palmerpenguins/); [Gorman et al. 2014](https://doi.org/10.1371/journal.pone.0090081)). 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. ```{r penguin-data, eval = has_penguins} 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) ``` ```{r penguin-raw-scale, eval = has_penguins, fig.width = 6, fig.height = 4, fig.cap = "Raw log-scale flipper length and body mass for analysed complete rows. This descriptive relationship is not the fitted residual association."} 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 ) ``` ```{r penguin-direct-lognormal, eval = has_penguins, warning = FALSE, message = FALSE} 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] drmTMB::check_drm(fit_log) # 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 )) ``` ```{r penguin-intervals, eval = has_penguins, echo = FALSE, fig.width = 7, fig.height = 3.8, fig.cap = "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.", fig.alt = "Three 95 percent confidence intervals for the direct log-residual correlation rho12. Wald, profile, and bootstrap estimates are all about 0.35. The bootstrap interval is annotated with its retained full-refit count."} at <- rev(seq_len(nrow(intervals))) old_par <- par(no.readonly = TRUE) par(mar = c(4.2, 7.4, 0.8, 2.6), bg = "white") x_limits <- range(c(0, intervals$lower, intervals$upper)) + c(-0.08, 0.18) plot( NA, xlim = x_limits, ylim = c(0.5, 3.5), yaxt = "n", bty = "n", xlab = expression(rho[12]), ylab = "" ) abline(v = 0, lty = 2, lwd = 1, col = "grey65") draw_interval_eye <- function(lower, estimate, upper, y) { x <- seq(lower, upper, length.out = 101) left_width <- max(estimate - lower, .Machine$double.eps) right_width <- max(upper - estimate, .Machine$double.eps) taper <- ifelse( x <= estimate, (x - lower) / left_width, (upper - x) / right_width ) half_height <- 0.18 * sqrt(pmax(taper, 0)) polygon( c(x, rev(x)), c(y + half_height, rev(y - half_height)), col = grDevices::adjustcolor("#087f8c", alpha.f = 0.28), border = NA ) points( estimate, y, pch = 21, cex = 1.50, lwd = 1.80, col = "#087f8c", bg = "white" ) } for (i in seq_len(nrow(intervals))) { draw_interval_eye(intervals$lower[i], intervals$estimate[i], intervals$upper[i], at[i]) } axis(2, at = at, labels = c("Wald", "Profile", "Bootstrap"), las = 1, tick = FALSE) text( intervals$upper[3] + 0.025, at[3], labels = sprintf("%d/%d retained", ci_bootstrap$bootstrap.n, bootstrap_total), adj = c(0, 0.5), cex = 0.82, col = "#4b5563" ) par(old_par) ``` `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 \(|\rho_{12}|=1\). 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](convergence.html) 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`: ```{r student-example, eval = FALSE} 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. ```{r staged-example, eval = FALSE} 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](cross-family.html) 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](cross-family.html) 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.