Model specification

Status — Reference

Mirrors drmTMB's Model specification (13 items in drmTMB). A model is one bf formula bundle (one linear predictor per distributional parameter) plus a response family. All 13 drmTMB families are available.

Formula bundle

DRM.bf Function

julia

bf(response_formula, dpar_formulas...)
drm_formula(response_formula, dpar_formulas...)

Bundle one formula per distributional parameter, exactly as drmTMB. The first formula y ~ … sets the response and the μ predictor; each later formula param ~ … (e.g. sigma ~ …) sets that parameter's predictor. sigma defaults to ~ 1 when omitted.

The secondary parameter on each formula's left-hand side must be one of sigma, nu, zi, hu, zoi, coi. Mirroring drmTMB, bf rejects reserved / mis-typed syntax with a clear error: tau (the scale is sigma), mu as a separate formula (μ comes from the response), the two-response parameters mu1/mu2/sigma1/sigma2/rho12 in this positional form (use the keyword form), unknown parameter names, and a parameter given more than once.

julia

bf(; mu1, mu2, sigma1=…, sigma2=…, nu=…, rho12=…)

Bivariate formula bundle, mirroring drmTMB. mu1 = y1 ~ … and mu2 = y2 ~ … set the two responses and their mean predictors; sigma1, sigma2 (log σ) and rho12 (atanh ρ) default to ~ 1. For one-sided predictors give the parameter name as a placeholder LHS, e.g. sigma1 = @formula(sigma1 ~ x).

nu is the shared degrees-of-freedom formula for the bivariate Student-t family (drm(…, Student()), drmTMB's biv_student()), on the logm2 scale ν = 2 + exp(η). It is omitted entirely unless supplied, so the Gaussian and lognormal bundles are unchanged. Under the exact bivariate-t density a single scalar mixing variable governs both margins, so ν is structurally shared — there is no per-margin nu1/nu2.

Like the univariate form, bf rejects reserved / mis-typed syntax: a placeholder LHS that is not its own parameter name (e.g. sigma1 = @formula(tau ~ x) or a swapped sigma1/sigma2), and a two-column cbind(…) response on mu1/mu2.

DRM.DrmFormula Type

julia

DrmFormula

A bundle of one linear-predictor formula per distributional parameter, built by bf. response is the response column; forms is ordered :mu => rhs, :sigma => rhs, ….

DRM.BivariateDrmFormula Type

julia

BivariateDrmFormula

Two-response formula bundle (μ1, μ2, σ1, σ2, ρ12), built by the keyword form of bf.

Response families

DRM.Gaussian Type

julia

Gaussian()

Gaussian response family: identity link on the mean μ, log link on the scale σ (so σ coefficients act on log σ). Mirrors drmTMB::gaussian().

DRM.Student Type

julia

Student()

Student-t response family: identity link on the location μ, log link on the scale σ, and the degrees of freedom as ν = 2 + exp(η) (so ν > 2 and the variance is always finite; ν coefficients act on log(ν − 2)). Note that σ is the scale, not the standard deviation: for ν > 2, SD[y] = σ·sqrt(ν/(ν − 2)) (mirroring drmTMB). Robust sibling of Gaussian — heavy tails downweight outliers, and ν → ∞ tends to Gaussian. Mirrors drmTMB's student family.

julia

fit = drm(bf(y ~ x, sigma ~ 1, nu ~ 1), Student(); data = dat)
2 + exp(coef(fit, :nu)[1])  # estimated degrees of freedom (ν = 2 + exp(η))

DRM.Poisson Type

julia

Poisson()

Poisson response family for counts: log link on the mean μ (so μ coefficients act on log λ). No scale parameter. Mirrors drmTMB's poisson family.

Note

DRM.Poisson shadows Distributions.Poisson; if you need the distribution too (e.g. to simulate), qualify it as Distributions.Poisson.

julia

fit = drm(bf(y ~ x), Poisson(); data = dat)
fitted(fit)        # fitted counts λ = exp(Xβ̂), on the response scale

fit_phy = drm(bf(@formula(y ~ x + phylo(1 | species))), Poisson();
              data = dat, tree = tr, se = false)

# Experimental (#136 Arc 0): Poisson random-intercept variational (ELBO) marginal.
# Default remains Laplace (`marginal = :LA`). `loglik` on a VA fit is an ELBO.
fit_va = drm(bf(@formula(y ~ x + (1 | g))), Poisson(); data = dat, marginal = :VA)

# Opt-in Cox–Reid restricted estimation (#443 / #450). ML is the default.
fit_reml = drm(bf(@formula(y ~ x + (1 | g))), Poisson(); data = dat, method = :REML)
estimation_method(fit_reml)   # :REML

# Opt-in 1-D Liu–Pierce AGHQ for `(1 | g)` only (#448). Default `:LA` stays
# today's non-adaptive GHQ-32. k=1 ≡ 1-point Laplace plumbing, not a recovery
# headline. Capability row stays missing. `:REML` is not wired to `:AGHQ`.
fit_aghq = drm(bf(@formula(y ~ x + (1 | g))), Poisson();
               data = dat, marginal = :AGHQ, nAGQ = 5)
fit_aghq.marginal   # :AGHQ

fit_phy_reml = drm(bf(@formula(y ~ x + phylo(1 | species))), Poisson();
                   data = dat, tree = tr, se = false, method = :REML)

Restricted (Cox–Reid) estimation

method = :REML is opt-in — ML remains the default. It maximises the Cox–Reid adjusted profile likelihood ℓ_ML − ½·log|I_ββ| on two Poisson routes:

a scalar random intercept (1 | g) integrated by GHQ-32 (#443)
phylogenetic / relmat / animal / precomputed-spatial Laplace (_fit_poisson_general_laplace, #450)

On the Gaussian route this correction is exactly Patterson–Thompson REML, so the mechanism is anchored rather than ad hoc.

Probe Cell D is not a recovery result

A 16-tip / 12-seed Poisson phylo cell over-corrected under Cox–Reid (ML +8.18%, CR +17.41%). Do not read those percentages as a bias-sign headline or a reason to prefer :REML on trees. ADEMP on a larger tree is a follow-on. Evidence: docs/dev-log/evidence/2026-08-18-cox-reid-scoping-probe.md.

It over-corrects when clusters are plentiful

On a Poisson (1 | g) cell with true σ_b = 0.6 and 6 observations per cluster, ML was −12.4% at G = 10 where Cox–Reid was −1.8% — but by G = 40 ML was +1.4% and Cox–Reid +4.4%. Reach for it when clusters are few, not by habit. That asymmetry is why ML is the default.

Crossed intercepts, correlated slopes (1 + x | g), coordinate-spatial with estimated range ρ, fixed-effects-only, zi/hu, marginal = :VA, and marginal = :AGHQ still error on method = :REML rather than silently returning an ML fit. REML log-likelihoods are not comparable across different fixed-effect structures, so do not use them for model selection over μ. This slice does not flip a capability chip.

DRM.NegBinomial2 Type

julia

NegBinomial2()

Negative-binomial (NB2) family for overdispersed counts: log link on the mean μ, and log link on the scale σ (the sigma formula slot, so coef(fit, :sigma) is log σ). The NB2 size is θ = 1/σ² = exp(-2·coef(:sigma)), matching the Gamma/Beta convention and drmTMB (size = 1/σ²). Var = μ + μ²/θ; as θ → ∞ it tends to Poisson. Mirrors drmTMB's nbinom2 family. Crossed random intercepts on the mean, such as (1 | g) + (1 | h), use the sparse-Laplace engine when sigma ~ 1. A phylogenetic random intercept on the mean, phylo(1 | species), also uses the sparse-Laplace engine; here a covariate dispersion formula sigma ~ x is supported (a per-observation log σ, so a per-observation size θ = 1/σ²; #164), while the crossed-intercept route still requires sigma ~ 1.

julia

fit = drm(bf(y ~ x, sigma ~ 1), NegBinomial2(); data = dat)
fit = drm(bf(y ~ x + (1 | g) + (1 | h), sigma ~ 1), NegBinomial2(); data = dat)
fit_phy = drm(bf(@formula(y ~ x + phylo(1 | species)), @formula(sigma ~ 1)),
              NegBinomial2(); data = dat, tree = tr, se = false)
fit_disp = drm(bf(@formula(y ~ x + phylo(1 | species)), @formula(sigma ~ x)),
               NegBinomial2(); data = dat, tree = tr, se = false)  # log σ ~ x
# Coupled phylogenetic location–scale (#202): shared tag + phylo group on both axes.
fit_ls = drm(bf(@formula(y ~ x + (1 | p | phylo(species))),
                @formula(sigma ~ 1 + (1 | p | phylo(species)))),
             NegBinomial2(); data = dat, tree = tr, se = false)
exp(-2 * coef(fit, :sigma)[1])  # estimated size θ = 1/σ²

# Experimental (#136 Rung 1): NB2 random-intercept variational (ELBO) marginal.
# Requires `sigma ~ 1`. Default remains Laplace (`marginal = :LA`).
fit_va = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), NegBinomial2();
             data = dat, marginal = :VA)

DRM.TruncatedNegBinomial2 Type

julia

TruncatedNegBinomial2()

Zero-truncated negative-binomial (NB2) family for strictly-positive counts (≥ 1) — litter sizes, group sizes given presence. Same parameterisation as NegBinomial2 (log-link mean μ, log σ in the sigma slot, size θ = 1/σ²) but conditioned on y ≥ 1: P(k) = NB(k) / (1 − NB(0)). Mirrors drmTMB's truncated_nbinom2.

julia

fit = drm(bf(y ~ x, sigma ~ 1), TruncatedNegBinomial2(); data = dat)

Adding an hu part fits the HURDLE NB2 instead: P(0) = logistic(X_huᵀβ) and the positive counts follow this zero-truncated NB2. Zeros are then allowed in the response. This is drmTMB's spelling of the hurdle model — it has no hurdle_nbinom2() constructor — and it fits the same likelihood as NegBinomial2 with hu, which is what the call delegates to (so the returned fit reports NegBinomial2 as its family).

julia

fit = drm(bf(y ~ x, sigma ~ 1, hu ~ w), TruncatedNegBinomial2(); data = dat)

julia

Beta()

Beta response family for proportions in (0,1): logit link on the mean μ, and the sigma slot carries σ with the precision mapping φ = 1/σ² (so coef(fit, :sigma) is log σ; recover precision as exp(-2·log σ)). Likelihood Beta(μφ, (1-μ)φ). Mirrors drmTMB's beta_family. Crossed random intercepts on the mean, such as (1 | g) + (1 | h), use the sparse-Laplace engine when sigma ~ 1. A phylogenetic random intercept on the mean, phylo(1 | species), also uses the sparse-Laplace engine, as do general user-supplied PD-covariance intercepts — relmat(1 | id) with K = C, the animal(1 | id) (A = C) and spatial(1 | id) (K = C) aliases.

Note

DRM.Beta shadows Distributions.Beta; qualify the latter if you need it.

julia

fit = drm(bf(y ~ x, sigma ~ 1), Beta(); data = dat)
fit = drm(bf(y ~ x + (1 | g) + (1 | h), sigma ~ 1), Beta(); data = dat)
fit_phy = drm(bf(@formula(y ~ x + phylo(1 | species)), @formula(sigma ~ 1)),
              Beta(); data = dat, tree = tr, se = false)
exp(-2 * coef(fit, :sigma)[1])     # estimated precision φ

# Experimental (#136 Rung 1): Beta random-intercept variational (ELBO) marginal.
# Requires `sigma ~ 1`. Default remains Laplace (`marginal = :LA`).
fit_va = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), Beta();
             data = dat, marginal = :VA)

DRM.BetaBinomial Type

julia

BetaBinomial()

Beta-binomial response family — successes out of known trials with extra-binomial overdispersion. Logit link on the mean success probability μ; the sigma slot carries σ with precision φ = 1/σ² (so coef(fit, :sigma) is log σ). Likelihood BetaBinomial(n, μφ, (1-μ)φ). Requires a two-column response via cbind. Mirrors drmTMB's beta_binomial. Crossed random intercepts on the mean, such as (1 | g) + (1 | h), use the sparse-Laplace engine when sigma ~ 1. A phylogenetic random intercept on the mean, phylo(1 | species), also uses the sparse-Laplace engine (#166); both routes are constant-σ (overdispersion) only for now.

julia

fit = drm(bf(cbind(successes, failures) ~ x, sigma ~ 1), BetaBinomial(); data = dat)
fit = drm(bf(cbind(successes, failures) ~ x + (1 | g) + (1 | h), sigma ~ 1), BetaBinomial(); data = dat)
fit_phy = drm(bf(@formula(cbind(successes, failures) ~ x + phylo(1 | species)), @formula(sigma ~ 1)),
              BetaBinomial(); data = dat, tree = tr, se = false)

DRM.Binomial Type

julia

Binomial()

Binomial response family — successes out of known trials (logistic regression). Logit link on the mean success probability μ; no scale/dispersion parameter (mean-only, like Poisson). Accepts either a two-column response via cbind (trials = successes + failures) or a plain 0/1 Bernoulli vector. Likelihood Binomial(n, μ) with μ = logistic(η). Mirrors drmTMB's binomial family. A random intercept (1 | g) on the mean fits a logistic GLMM; crossed intercepts such as (1 | g) + (1 | h) use the sparse-Laplace engine. A phylogenetic random intercept on the mean, phylo(1 | species), also uses the sparse-Laplace engine.

Note

DRM.Binomial shadows Distributions.Binomial; if you need the distribution too (e.g. to simulate), qualify it as Distributions.Binomial.

julia

fit = drm(bf(cbind(successes, failures) ~ x), Binomial(); data = dat)   # logistic regression
fit = drm(bf(y ~ x + (1 | g)), Binomial(); data = dat)                  # 0/1 logistic GLMM
fit = drm(bf(cbind(successes, failures) ~ x + (1 | g) + (1 | h)), Binomial(); data = dat)
fit_phy = drm(bf(@formula(cbind(successes, failures) ~ x + phylo(1 | species))),
              Binomial(); data = dat, tree = tr, se = false)
fitted(fit)        # fitted success probabilities μ̂ = logistic(Xβ̂)

# Experimental (#136 Rung 1): Binomial random-intercept variational (ELBO) marginal.
fit_va = drm(bf(@formula(y ~ x + (1 | g))), Binomial(); data = dat, marginal = :VA)

julia

Gamma()

Gamma response family for positive continuous data: log link on the mean μ, and the sigma slot carries σ = the coefficient of variation, mapped to the shape α = 1/σ² (coef(fit, :sigma) is log σ; recover the shape as exp(-2·log σ)). Likelihood Gamma(α, μ/α); variance μ²σ². Mirrors drmTMB's Gamma family. Crossed random intercepts on the mean, such as (1 | g) + (1 | h), use the sparse-Laplace engine when sigma ~ 1. A phylogenetic random intercept on the mean, phylo(1 | species), also uses the sparse-Laplace engine.

Note

DRM.Gamma shadows Distributions.Gamma; qualify the latter if you need it.

julia

fit = drm(bf(y ~ x, sigma ~ 1), Gamma(); data = dat)
fit = drm(bf(y ~ x + (1 | g) + (1 | h), sigma ~ 1), Gamma(); data = dat)
fit_phy = drm(bf(@formula(y ~ x + phylo(1 | species)), @formula(sigma ~ 1)),
              Gamma(); data = dat, tree = tr, se = false)
exp(-2 * coef(fit, :sigma)[1])     # estimated shape α

# Experimental (#136 Rung 1): Gamma random-intercept variational (ELBO) marginal.
# Requires `sigma ~ 1`. Default remains Laplace (`marginal = :LA`).
fit_va = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), Gamma();
             data = dat, marginal = :VA)

DRM.LogNormal Type

julia

LogNormal()

LogNormal response family for positive continuous data: the mean formula μ is the mean of log y (identity link on the log scale), and σ (log link) is the SD of log y. The response-scale median is exp(μ). Mirrors drmTMB's lognormal family.

Note

DRM.LogNormal shadows Distributions.LogNormal; qualify the latter if needed.

A random intercept (1 | g) or a correlated random intercept+slope (1 + x | g) may be placed on the log-mean μ; the group effect is integrated out by Gauss–Hermite quadrature (re_sd(fit)[:g] / vc(fit)[:g]).

Structured markers (phylo/relmat)

phylo(1 | group) (needs tree = ...) and relmat(1 | group) (needs K = ...) may be placed on the mean, exactly as drmTMB's lognormal family supports. log(y) is exactly Gaussian, so a structured call delegates WHOLESALE to drm(f, Gaussian(); data = data-with-logged-response, tree = ..., K = ...) — the same identity biv_lognormal() already uses (src/bivariate_lognormal.jl) — and shifts the reported log-likelihood by the parameter-free Jacobian -sum(log y). theta/vcov/ranef and the fitted objective's gradient carry over untouched from the Gaussian-on-log(y) fit; only loglik (and everything derived from it: aic/bic/deviance) shifts. animal/spatial markers and meta_V are not implemented for LogNormal() (no R parity cell) and stay refused, as does any marker on a formula other than the mean.

julia

fit = drm(bf(y ~ x, sigma ~ 1), LogNormal(); data = dat)
coef(fit, :mu)                 # on the log scale
exp(coef(fit, :sigma)[1])      # SD of log y

fit_phy = drm(bf(@formula(y ~ x + phylo(1 | species)), @formula(sigma ~ 1)),
              LogNormal(); data = dat, tree = tr)

DRM.ZeroOneBeta Type

julia

ZeroOneBeta()

Zero-one-inflated beta family for proportions on the closed interval [0,1] (values may be exactly 0 or 1). Parameters: mu (logit; beta mean on the interior), sigma (log; precision φ = 1/σ²), zoi (logit; probability the value is a boundary 0/1), coi (logit; probability of 1 given a boundary). Mirrors drmTMB's zero_one_beta. fitted returns the unconditional mean (1-zoi)·μ + zoi·coi.

julia

fit = drm(bf(y ~ x, sigma ~ 1, zoi ~ 1, coi ~ 1), ZeroOneBeta(); data = dat)

DRM.Tweedie Type

julia

Tweedie()

Tweedie response family (compound Poisson–Gamma, power 1 < p < 2) for semicontinuous data — an exact-zero mass plus a positive continuous part. Log link on the mean μ; sigma is the √dispersion (so φ = σ², coef(fit, :sigma) is log σ); nu is the power p on a logit-(1,2) link (p = 1 + logistic(coef(fit, :nu))). Var(y) = φ·μ^p. Density via the Dunn–Smyth series. Mirrors drmTMB's tweedie. An ordinary (1 | g) random intercept and an independent (0 + x | g) random slope on mu are both supported (32-node Gauss–Hermite marginal, #563); sigma ~ 1 and nu ~ 1 remain fixed-effect sub-models on either route. The CORRELATED random slope (1 + x | g), random effects on sigma/nu, and structured (phylo/relmat/animal/spatial) markers on mu are not implemented — matches drmTMB 0.7.0, which rejects the same forms ("Only independent tweedie() mu random intercepts and slopes are implemented in this slice").

julia

fit = drm(bf(y ~ x, sigma ~ 1, nu ~ 1), Tweedie(); data = dat)
exp(2 * coef(fit, :sigma)[1])               # dispersion φ
1 + 1 / (1 + exp(-coef(fit, :nu)[1]))       # power p ∈ (1,2)

fit_re = drm(bf(y ~ x + (1 | g), nu ~ 1), Tweedie(); data = dat)
exp(coef(fit_re, :resd)[1])                 # random-intercept SD

fit_slope = drm(bf(y ~ x + (0 + x | g), nu ~ 1), Tweedie(); data = dat)
exp(coef(fit_slope, :resd)[1])              # random-slope SD

DRM.CumulativeLogit Type

julia

CumulativeLogit()

Cumulative-logit (proportional-odds) ordinal family. The response is an ordered category coded 1, 2, …, K. Pr(y ≤ k) = logistic(θ_k − η) with ordered cutpoints θ_1 < … < θ_{K-1}; the single linear predictor η comes from the mu formula with its intercept dropped (cutpoints absorb it). No scale parameter. coef(fit, :mu) are the slopes; coef(fit, :cutpoints) are the raw increment parameters (θ_1 = δ_1, θ_k = θ_{k-1} + exp(δ_k)). fitted returns the expected ordered-category score Σ_k k·Pr(y=k). Mirrors drmTMB's cumulative_logit. An ordinary random intercept (1 | g) or an independent random slope (0 + x | g) on mu integrates the group-level term out by 32-node Gauss–Hermite quadrature; coef(fit, :resd) is the log random-effect SD. An unlabelled, intercept-only phylogenetic random intercept phylo(1 | species) on mu is fit via the sparse-Laplace GLMM route instead (needs tree = …); it cannot be combined with an ordinary random effect. For the phylo route, re_sd(fit)[:group] is on the RAW branch-length scale (tip variance equals the tree height h, not 1) — this differs from drmTMB's sdpars$mu[["phylo(1 | species)"]], which is on the CORRELATION scale (ape::vcv(tree, corr = TRUE), tip variance 1 regardless of h). Convert with re_sd(fit)[:group] * sqrt(phylo_tree_height(augmented_phy(tree))) to compare against drmTMB's number (the same convention used throughout the Gaussian phylo-mean route, e.g. test_parity_gaussian_phylo_mean.jl). Correlated slopes (1 + x | g), crossed/multiple random effects, and relmat/animal/spatial structured markers are not implemented.

julia

fit = drm(bf(y ~ x), CumulativeLogit(); data = dat)          # y coded 1..K
fit_re = drm(bf(y ~ x + (1 | g)), CumulativeLogit(); data = dat)
fit_slope = drm(bf(y ~ x + (0 + x | g)), CumulativeLogit(); data = dat)
fit_phylo = drm(bf(y ~ x + phylo(1 | species)), CumulativeLogit(); data = dat, tree = tr)

Advanced family type

SkewNormal is a Julia family type with its own docstring. Its presence here does not widen the supported-model matrix; consult the capability matrix before selecting a family and structure.

The source docstring below uses a shorthand formula sketch. For a call with prepared data dat, use explicit Julia formula objects:

julia

fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(nu ~ 1)),
          SkewNormal(); data = dat)

DRM.SkewNormal Type

julia

SkewNormal()

Skew-normal response family for continuous, asymmetric data: identity link on the location μ, log link on the scale σ (the sigma slot), and identity link on the slant ν (the nu slot). The public parameterisation is the moment form mirroring drmTMB's skew_normal: μ is the mean of y, σ is the SD of y, and ν is Azzalini's slant α. These are mapped internally to the location–scale–slant (ξ, ω, α) of the density f(y) = 2·φ((y−ξ)/ω)·Φ(α(y−ξ)/ω)/ω.

ν = 0 recovers the symmetric Gaussian; ν > 0 skews right, ν < 0 skews left.

julia

fit = drm(bf(y ~ x, sigma ~ 1, nu ~ 1), SkewNormal(); data = dat)
coef(fit, :mu)[1]           # mean of y (identity)
exp(coef(fit, :sigma)[1])   # SD of y
coef(fit, :nu)[1]           # estimated slant α (identity)

Note

Near ν = 0 the model is (locally) symmetric Gaussian and the slant is only weakly identified — its standard error inflates and recovery is loose. The fitter starts ν at a small nonzero value seeded by the sample-skewness sign to break this symmetry.

Two-column response

DRM.cbind Function

julia

cbind(successes, failures)

Formula marker for a two-column count response, e.g. bf(cbind(successes, failures) ~ x, sigma ~ z) with BetaBinomial. Trials are successes + failures. Mirrors drmTMB's cbind(...) response. Only inspected structurally on a formula left-hand side.

Modelled missing predictors

Experimental

Exported for evaluation; fenced for v1.0 (D-181). API and numerics may change; not covered by the R-parity scoreboard.

For a Gaussian response, mi(x) marks an additive predictor whose missing values are integrated out under a joint model. Observed predictor values inform the predictor distribution; missing predictors are not filled in before fitting. This route supports one Gaussian or Bernoulli predictor, or two independent Gaussian predictor models, with complete remaining exogenous fixed-effect covariates. Neither the response nor either modelled predictor may appear in those fixed designs. Other response families, three or more predictors, mixed predictor families, random or structured effects, REML, and profile/bootstrap intervals are still outside this admission. The same narrow route is available as development through drmTMB(..., engine = "julia"); it does not establish full native fitted-result parity.

julia

using DRM, LinearAlgebra
BLAS.set_num_threads(1)
n = 32
z = collect(range(-1.2, 1.2; length=n))
x_full = 0.15 .+ 0.65 .* z .+ 0.15 .* sin.(1:n)
y_full = 0.3 .+ 0.4 .* z .+ 0.7 .* x_full .+ 0.18 .* cos.(1:n)
data = (y = Union{Missing,Float64}[i in (4, 24) ? missing : y_full[i] for i in 1:n],
        x = Union{Missing,Float64}[i in (8, 18, 24) ? missing : x_full[i] for i in 1:n], z=z)
fit = drm(bf(@formula(y ~ z + mi(x)), @formula(sigma ~ 1)), Gaussian();
    data=data, impute=(x=@formula(x ~ z),),
    missing=miss_control(response="include", predictor="model"))
@assert is_converged(fit)
@assert isposdef(Symmetric(vcov(fit)))
(coef(fit, :mu), coef(fit, :mi_x), coef(fit, :sigma_mi_x), imputed(fit))

([0.2956124851028067, 0.37501846625780827, 0.7485514011793478], [0.15689535402141302, 0.6447678352387548], [0.10189508204961784], (variable = ["x", "x", "x"], original_row = [8, 18, 24], model_row = [8, 18, 24], observed = Bool[0, 0, 0], estimate = [-0.24247894915511675, 0.2428178772879237, 0.5312766777084319], std_error = Union{Missing, Float64}[0.09074205154086845, 0.08920306612733002, 0.10470157533534986], source = ["conditional_mode", "conditional_mode", "conditional_mode"], uncertainty_status = ["ok", "ok", "ok"]))

For a binary predictor use impute=(x=impute_model(@formula(x ~ z); family=Binomial()),). miss_control(predictor="model") retains the default response="fail" when the response is complete; request response="include" for missing responses. Unmarked incomplete predictors and unsupported options produce an error.

imputed(fit; rows=:all) retains original row numbers and includes observed values, whose SEs are missing. Gaussian imputation SEs combine conditional variance with a first-order parameter-uncertainty correction. Bernoulli summaries return conditional probabilities and their Bernoulli SDs. These are neither multiple-imputation draws nor interval-coverage guarantees. Check uncertainty_status; se=false hides SEs without erasing a failure status. coef(fit, :sigma_mi_x) returns natural predictor SD; unqualified coef(fit) and vcov(fit) use the raw coordinates, including its log SD. Complete fixed-effect interactions/transforms are allowed; interactions involving mi(x) are not yet admitted.

For the R bridge, write the corresponding R formula and use engine = "julia", impute = list(x = x ~ z) (or impute = list(x = impute_model(x ~ z, family = binomial())) for a binary predictor), and missing = miss_control(response = "drop", predictor = "model") or miss_control(response = "include", predictor = "model"). The R response-drop preprocessing differs from native-TMB behaviour and is not native parity. Its profile/bootstrap methods are explicitly unsupported. The Gaussian predictor-SD Wald interval is a natural-scale delta interval, can cross zero, and is neither native-interval parity nor coverage evidence.

With two Gaussian predictors, mark each separately and provide both predictor models. The impute entry order does not change which model belongs to which variable. Their prior distributions are conditionally independent given their covariates; conditioning on an observed response can correlate their missing values. The prepared fit retains that full conditional covariance.

julia

using DRM, Random, LinearAlgebra
BLAS.set_num_threads(1)
rng = MersenneTwister(9302)
n = 64
z = collect(range(-1, 1; length=n))
x1 = 0.2 .+ 0.4z .+ 0.5randn(rng, n)
x2 = -0.1 .+ 0.3z .+ 0.6randn(rng, n)
y = 0.3 .+ 0.6x1 .- 0.4x2 .+ 0.2z .+ 0.5randn(rng, n)
data = (y=Union{Missing,Float64}[i % 13 == 0 ? missing : y[i] for i in 1:n],
        x1=Union{Missing,Float64}[i % 5 == 0 ? missing : x1[i] for i in 1:n],
        x2=Union{Missing,Float64}[i % 7 == 0 ? missing : x2[i] for i in 1:n], z=z)
fit = drm(bf(@formula(y ~ z + mi(x1) + mi(x2))), Gaussian(); data=data,
    impute=(x2=@formula(x2 ~ z), x1=@formula(x1 ~ z)),
    missing=miss_control(response="include", predictor="model"))
@assert is_converged(fit)
(coef(fit, :sigma_mi_x1), coef(fit, :sigma_mi_x2),
 imputed(fit; variable=:x2))

([0.5109814606591342], [0.5462545159810994], (variable = ["x2", "x2", "x2", "x2", "x2", "x2", "x2", "x2", "x2"], original_row = [7, 14, 21, 28, 35, 42, 49, 56, 63], model_row = [7, 14, 21, 28, 35, 42, 49, 56, 63], observed = Bool[0, 0, 0, 0, 0, 0, 0, 0, 0], estimate = [-0.29345510184587387, 0.375875419892618, 0.41354475469566565, -0.46248422948769086, 0.41618947536448614, 0.06592534270473177, -0.03289463233738368, -0.18891817752408835, 0.016486018941177183], std_error = Union{Missing, Float64}[0.44479454549246544, 0.44707555287542833, 0.44011705349308594, 0.44332052884703266, 0.4676349885435535, 0.43420219264865206, 0.43703430213853767, 0.444213997579822, 0.44476566844270216], source = ["conditional_mode", "conditional_mode", "conditional_mode", "conditional_mode", "conditional_mode", "conditional_mode", "conditional_mode", "conditional_mode", "conditional_mode"], uncertainty_status = ["ok", "ok", "ok", "ok", "ok", "ok", "ok", "ok", "ok"]))

With two predictors, imputed requires variable. Direct Julia returns raw coefficients and covariance in the prepared order: exogenous mean coefficients, the two marker slopes in formula order, residual scale coefficients, then each predictor's coefficients and log SD. Natural-SD accessors do not change that raw covariance convention. The R bridge instead restores native model-matrix column order and transforms both predictor SDs and their full covariance to public natural scales. In R use impute=list(x1=x1~z, x2=x2~z) and imputed(fit, variable="x2").

DRM.mi Function

julia

mi(x)

Marker for exactly one additive predictor whose missing values are modelled by the joint missing-predictor frontend. It is valid only inside @formula; direct evaluation is an error.

DRM.miss_control Function

julia

miss_control(; response = "fail", predictor = "fail")

Declare missing-data handling for the joint frontend. Its currently verified admission requires predictor = "model" and, when responses are missing, response = "include"; the defaults fail rather than silently changing a conventional drm() fit.

DRM.impute_model Function

julia

impute_model(formula; family = Gaussian(), levels = nothing)

Wrap a predictor model when its distribution is not the default Gaussian. The bounded frontend admits Gaussian, Bernoulli (Binomial()), ordinal (CumulativeLogit()), and categorical (CategoricalLogit()) predictors. levels is required for textual ordinal data, recommended for categorical data, and rejected for Gaussian/Bernoulli predictors.

DRM.imputed Function

julia

imputed(fit::PreparedJointFit; variable = :x, rows = :missing, se = true)

Return a Tables-compatible column table of modelled predictor summaries, with variable, original_row, model_row, observed, estimate, std_error, source, and uncertainty_status. rows=:all includes observed predictors, whose standard errors are missing; the default returns missing-predictor rows.

Gaussian means are conditional modes. Their standard errors add first-order parameter uncertainty to conditional variance, matching TMB's prediction-error approximation. Bernoulli estimates are conditional probabilities and use the conditional Bernoulli standard deviation. Neither route returns multiple imputations, posterior intervals, or exact integration over parameter uncertainty.

se=false omits all standard errors but retains fit-level uncertainty failures. The prepared-array interface names its predictor x; formula fits retain their user-supplied predictor name.

Conditional imputation for the marker variable in a formula-based joint fit.

Conditional imputation for one explicitly selected two-predictor marker variable.

Conditional finite-state imputation for the marker variable in a formula-based fit.

DRM.JointDrmFit Type

User-facing wrapper around the prepared exact joint fit.

DRM.JointTwoDrmFit Type

julia

JointTwoDrmFit

Formula-facing wrapper for the prepared two-independent-Gaussian-predictor kernel. Raw coefficients and vcov preserve the kernel order beta, b1, b2, delta, alpha1, logtau1, alpha2, logtau2; natural predictor SDs are exposed only through coef(fit, :sigma_mi_<variable>).

DRM.JointMissingControl Type

Validated controls for the joint missing-predictor formula admission.

DRM.JointImputeModel Type

Formula and family for the one predictor distribution in a joint fit.