Model specification
Status — Reference
Mirrors drmTMB's Model specification. A model is one bf formula bundle (one linear predictor per distributional parameter) plus a response family. DRModels.jl has fifteen response-family constructors: the fourteen under Response families and SkewNormal under Advanced family type.
Formula bundle
DRModels.bf Function
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.
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.
DRModels.drm_formula Function
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.
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.
DRModels.DrmFormula Type
DrmFormulaA bundle of one linear-predictor formula per distributional parameter, built by bf. response is the response column; forms is ordered :mu => rhs, :sigma => rhs, ….
DRModels.BivariateDrmFormula Type
BivariateDrmFormulaTwo-response formula bundle (μ1, μ2, σ1, σ2, ρ12), built by the keyword form of bf. The schema_cache field holds the training-data factor schema written at fit time, so predict / predict_parameters keep the fitted contrasts when newdata is a level subset (#1000).
DRModels.offset Function
offset(x)Formula marker for a known offset in the mean formula, for example y ~ x + offset(log_exposure). The term enters the mean's linear predictor with its coefficient fixed at 1 (η = Xμβ + offset); it is not estimated, as in drmTMB. Every univariate family applies it on the fixed-effects mean (for Gaussian() that is exactly the fit of y - offset). Some families also apply it with random effects: Poisson() with one ordinary (1 | g) (marginal = :LA, :AGHQ or :Laplace); NegBinomial2() with zi/hu and with one ordinary (1 | g) (:LA or :Laplace); Gamma() and Beta() with one (1 | g) under marginal = :Laplace. Every other route raises an ArgumentError that names it. An offset in any formula other than the mean (sigma, nu, zi, ...) is refused, as drmTMB refuses it. At most one offset(...) term is allowed. predict(fit, newdata) needs the offset column in newdata. Outside a formula, offset(x) returns x unchanged.
Other names used in formulas
These exported names have no separate entry below:
@formula(...)is re-exported from StatsModels.jl; it writes one formula, as~does in R. See the StatsModels.jl documentation.drm_formulais another name forbf(drmTMB's spelling); the two are the same function.offset(x)inside the mean formula addsxto the linear predictor with its coefficient fixed at 1, as in R, for example@formula(y ~ x + offset(log_exposure)); at most one offset term is allowed. Every univariate family applies it on the fixed-effects mean. With random effects it is applied forPoisson()with one ordinary(1 | g), forNegBinomial2()with one ordinary(1 | g)(:LAor:Laplace), and forGamma()andBeta()with one(1 | g)undermarginal = :Laplace. Fixed-effectsNegBinomial2()routes withziorhualso support offsets; they cannot be combined with a random effect. Other routes refuse offsets by name (seeoffset). An offset insigmaor any other non-mean formula is refused, as in drmTMB.
Response families
DRModels.Gaussian Type
Gaussian()Gaussian response family: identity link on the mean μ, log link on the scale σ (so σ coefficients act on log σ). Mirrors drmTMB::gaussian().
# Growth curves: a correlated random intercept + slope on the mean AND a random
# intercept on sigma, same grouping factor (#960; per-group 3-D Laplace, as drmTMB):
fit = drm(bf(@formula(weight ~ 1 + time + diet + (1 + time | chick)),
@formula(sigma ~ 1 + time + diet + (1 | chick))), Gaussian(); data = dat)
vc(fit)[:chick] # 2×2 intercept/slope covariance
exp(coef(fit, :resd)[1]) # SD of the sigma random interceptDRModels.Student Type
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.
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(nu ~ 1)), Student(); data = dat)
2 + exp(coef(fit, :nu)[1]) # estimated degrees of freedom (ν = 2 + exp(η))
predict_parameters(fit, dat)[:nu] # same quantity (response-scale ν, not exp(η))
# Opt-in TMB-convention Laplace for one ordinary `(1 | g)` on the mean (what drmTMB fits):
fit_lap = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), Student(); data = dat,
marginal = :Laplace)
# A group-level covariate can model an ordinary mean intercept's SD:
# `dat` must include numeric `z` values that are constant within each `g` group.
fit_sd = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1),
@formula(nu ~ 1), @formula(sd(g) ~ z)), Student(); data = dat)
# It can also model the coefficient SD of one independent numeric mean slope:
fit_slope_sd = drm(bf(@formula(y ~ x + (0 + x | g)), @formula(sigma ~ 1),
@formula(nu ~ 1), @formula(sd(g, coef(x)) ~ z)),
Student(); data = dat)
# Independent random intercept and slope on one factor (#960, as drmTMB fits it; Laplace):
fit_ind = drm(bf(@formula(y ~ x + (1 | g) + (0 + x | g)), @formula(sigma ~ 1)), Student(); data = dat)
exp.(coef(fit_ind, :resd)) # SDs of the intercept and slope effectsStudent-t SD regression currently supports one ordinary mean intercept or one independent numeric mean slope, after the family's existing observed-response filtering, with predictors that are constant within each retained group. For multiple possible coefficients, use the matching coef(intercept) or coef(x) selector. A slope formula models the SD of its random coefficient; the observation-level contribution also depends on the slope predictor. This route requires fixed sigma ~ 1 and nu ~ 1 submodels. It uses adaptive Gauss–Hermite quadrature (marginal = :LA); correlated blocks, structured effects, and marginal = :Laplace are not admitted on this route yet.
DRModels.Poisson Type
Poisson()Poisson response family for counts: log link on the mean μ (so μ coefficients act on log λ). No residual scale parameter. Mirrors drmTMB's poisson family.
Note
DRModels.Poisson shadows Distributions.Poisson; if you need the distribution too (e.g. to simulate), qualify it as Distributions.Poisson.
fit = drm(bf(@formula(y ~ x)), Poisson(); data = dat)
fitted(fit) # fitted counts λ = exp(Xβ̂), on the response scale
# Regress the SD of an ordinary Poisson random intercept on a group-level
# predictor. Each value of z must be constant within its group g.
fit_sd = drm(bf(@formula(y ~ x + (1 | g)), @formula(sd(g) ~ z)),
Poisson(); data = dat)
length(predict_parameters(fit_sd, dat)[:sd]) # one value per observation
# A known-exposure offset (#727, drmTMB twin gap): `offset(log_exposure)`'s
# coefficient is fixed at 1, so λ = exposure · exp(Xβ̂). Fixed-effects-only mean, or
# (#941) with one ordinary random intercept — the standard rate GLMM:
fit_off = drm(bf(@formula(y ~ x + offset(log_exposure))), Poisson(); data = dat)
fit_rate = drm(bf(@formula(y ~ x + offset(log_exposure) + (1 | g))), Poisson(); data = dat)
fit_rate_tmb = drm(bf(@formula(y ~ x + offset(log_exposure) + (1 | g))), Poisson();
data = dat, marginal = :Laplace) # drmTMB's Laplace objective
fit_phy = drm(bf(@formula(y ~ x + phylo(1 | species))), Poisson();
data = dat, tree = tr, se = false)
# A predictor-dependent SD for one animal slope-only mean coefficient.
# `sd(animal, animal) ~ habitat` models the coefficient SD by animal;
# the row-level random contribution is `x * b_animal`.
fit_animal_sd = drm(
bf(@formula(y ~ x + animal(0 + x | animal)),
@formula(sd(animal, animal) ~ habitat)),
Poisson(); data = dat, A = A)
# Experimental: Poisson random-intercept variational (ELBO) marginal.
# Default remains `marginal = :LA` (adaptive GHQ, 5 nodes per group, on an ordinary `(1 | g)`). `loglik` on a VA fit is an ELBO.
fit_va = drm(bf(@formula(y ~ x + (1 | g))), Poisson(); data = dat, marginal = :VA)
# TMB-convention Laplace on an ordinary `(1 | g)`, as native drmTMB fits it (ML only).
fit_lap = drm(bf(@formula(y ~ x + (1 | g))), Poisson(); data = dat, marginal = :Laplace)
# Opt-in Cox–Reid restricted estimation. 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, with `nAGQ` nodes. The default
# `:LA` is already per-group adaptive GHQ with 5 nodes. k=1 ≡ 1-point Laplace; this
# identity does not establish parameter recovery. `:REML` is not available with `: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)For an ordinary random-intercept route, the group SD follows log(s_g) = Z_g α and uses adaptive Gauss–Hermite quadrature. One independent numeric random slope is also supported. Sparse-Laplace routes support a predictor-dependent SD for one phylogenetic or related-matrix mean intercept, and an animal mean intercept or slope-only coefficient. For an animal slope, sd(animal, animal) ~ z models the SD of the coefficient multiplying x; each row contributes xᵢ b_g to the log mean. The SD predictors must be constant within animal. The animal SD multiplies the supplied covariance A: the marginal animal effect covariance is D A D, so non-unit diagonal values in A affect the marginal SD. Ordinary IID intercept and slope SD routes, and structured phylogenetic, related-matrix, and animal SD routes, accept a fixed offset(...) in the mean formula. For an animal intercept, callers may instead supply a positive-definite precision matrix with Ainv; Julia uses it directly as Q in the latent covariance D Q⁻¹ D. Its rows and columns must follow the first-seen animal order. The non-SD Poisson animal-slope route uses the same covariance scaling. Correlated blocks, zero-inflation/hurdle terms, and REML for animal slopes are not supported with these routes. Profile and bootstrap methods are available as fit-level plumbing; their interval coverage has not been calibrated.
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 per-group adaptive GHQ (5 nodes)phylogenetic /
relmat/animal/ precomputed-spatial Laplace (_fit_poisson_general_laplace)
On the Gaussian route this correction is exactly Patterson–Thompson REML, so the mechanism is anchored rather than ad hoc.
Limited evidence for phylogenetic models
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. Recovery on larger trees has not been established. This is a small exploratory result, not a general performance claim.
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 μ.
DRModels.NegBinomial2 Type
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. NB2 cannot represent dispersion below the Poisson limit (variance below the mean): on underdispersed counts θ → ∞ and the fit sits at the Poisson log-likelihood. For underdispersed counts use glmmTMB's compois or genpois (as drmTMB advises, drmTMB#1466); DRModels.jl has no such family yet (#976). 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/σ²), while the crossed-intercept route still requires sigma ~ 1.
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), NegBinomial2(); data = dat)
fit = drm(bf(@formula(y ~ x + (1 | g) + (1 | h)), @formula(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: 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: NB2 random-intercept variational (ELBO) marginal.
# Requires `sigma ~ 1`. Default remains `marginal = :LA` (adaptive GHQ, 5 nodes per group, on an ordinary `(1 | g)`).
fit_va = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), NegBinomial2();
data = dat, marginal = :VA)
# TMB-convention Laplace on an ordinary `(1 | g)`, as native drmTMB fits it.
# Requires `sigma ~ 1`; ML only.
fit_lap = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), NegBinomial2();
data = dat, marginal = :Laplace)
# Predictor-dependent SDs are available for ordinary intercepts and
# independent numeric slopes, and for an animal-model mean intercept with a
# supplied covariance matrix `A` or precision matrix `Ainv`. The animal route
# uses sparse Laplace and requires `sigma ~ 1`. `A` and `Ainv` are used as
# supplied (raw kernel, as in native drmTMB and the Poisson and Gamma animal
# routes), as is a `relmat` `K`.
fit_animal_sd = drm(bf(@formula(y ~ x + animal(1 | id)),
@formula(sigma ~ 1),
@formula(sd(id, animal) ~ habitat)),
NegBinomial2(); data = dat, A = A, se = false)
# `coef(fit_animal_sd, :sd)` models the log SD multiplier by habitat. The
# reported `:sd_group` values multiply the raw supplied `A` kernel. Interval coverage for predictor-dependent SDs is not yet calibrated.
# A predictor-dependent SD can also describe one ordinary random coefficient
# on the log-`sigma` predictor. The mean is fixed-effect-only, and the `sd`
# predictors are fixed and constant within group:
fit_sigma_sd = drm(bf(@formula(y ~ x),
@formula(sigma ~ 1 + (1 | g)),
@formula(sd(g, sigma) ~ habitat)),
NegBinomial2(); data = dat)
fit_sigma_slope_sd = drm(bf(@formula(y ~ x),
@formula(sigma ~ 1 + (0 + z | g)),
@formula(sd(g, sigma, coef(z)) ~ habitat)),
NegBinomial2(); data = dat)
# The sigma-side route uses 11-node adaptive GHQ by default; `nAGQ` overrides it.
# `sigma(fit_sigma_sd)` reports conditional scales at group modes.
# Interval coverage is not calibrated.
# Known exposure: `offset(z)` enters log μ with its coefficient fixed at 1 (#938).
# Applied on the fixed-effects mean, with `zi`/`hu`, and with one ordinary `(1 | g)`
# (`:LA` or `:Laplace`); other routes, and an offset in `sigma`, are refused.
fit_off = drm(bf(@formula(y ~ x + offset(log_exposure)), @formula(sigma ~ 1)),
NegBinomial2(); data = dat)DRModels.TruncatedNegBinomial2 Type
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.
fit = drm(bf(@formula(y ~ x), @formula(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).
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(hu ~ w)), TruncatedNegBinomial2(); data = dat)DRModels.TruncatedPoisson Type
TruncatedPoisson()Zero-truncated Poisson family for strictly-positive counts (≥ 1) — litter sizes, group sizes given presence. Log link on the mean λ, no dispersion parameter, conditioned on y ≥ 1: P(k) = Poisson(k; λ) / (1 − e^{-λ}). Mirrors drmTMB's truncated_poisson().
fit = drm(bf(@formula(y ~ x)), TruncatedPoisson(); data = dat)For one ordinary random intercept or independent numeric mean slope, sd(group) ~ z models its log SD using predictors constant within group. Select the random-slope coefficient explicitly with coef(x):
fit = drm(bf(@formula(y ~ x + (1 | site)), @formula(sd(site) ~ habitat)),
TruncatedPoisson(); data = dat)
fit_slope = drm(bf(@formula(y ~ x + (0 + x | site)),
@formula(sd(site, coef(x)) ~ habitat)),
TruncatedPoisson(); data = dat)This route uses adaptive one-dimensional quadrature. Correlated, structured and other random-effect blocks are not yet supported for TruncatedPoisson SD regression. For a slope, the returned SD describes the random coefficient multiplied by the row-level value of x in the linear predictor. Missing or NaN responses are dropped with a warning, following the package's observed-row policy; missing predictors still error.
Adding an hu part fits the HURDLE Poisson instead: P(0) = logistic(X_huᵀβ) and the positive counts follow this zero-truncated Poisson. Zeros are then allowed in the response. This is the shared drmTMB / DRModels.jl spelling of the hurdle Poisson, and it fits the same likelihood as Poisson with hu, which is what the call delegates to (so the returned fit reports Poisson as its family). Poisson() + hu keeps working.
fit = drm(bf(@formula(y ~ x), @formula(hu ~ w)), TruncatedPoisson(); data = dat)DRModels.Beta Type
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. An ordinary mean random intercept or one independent numeric mean slope also supports a predictor-dependent SD through sd(group) ~ predictors. A slope targets its coefficient SD with sd(group, coef(x)) ~ predictor; the slope loading is applied at each row. A supplied-A animal random intercept also supports sd(group, animal) ~ predictors with sigma ~ 1; the animal covariance is used as supplied (raw), as native drmTMB does. These SD routes require complete responses and fixed sigma ~ 1. Correlated SD blocks, phylogenetic SDs, and calibrated interval coverage remain unsupported. 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 (used as supplied, not rescaled; the precomputed spatial kernel keeps its rescaling) — relmat(1 | id) with K = C, the animal(1 | id) (A = C) and spatial(1 | id) (K = C) aliases.
Note
DRModels.Beta shadows Distributions.Beta; qualify the latter if you need it.
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), Beta(); data = dat)
fit = drm(bf(@formula(y ~ x + (1 | g) + (1 | h)), @formula(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)
fit_animal_sd = drm(bf(@formula(y ~ x + animal(1 | id)),
@formula(sigma ~ 1),
@formula(sd(id, animal) ~ z)),
Beta(); data = dat, A = A)
exp(-2 * coef(fit, :sigma)[1]) # estimated precision φ
# Experimental: Beta random-intercept variational (ELBO) marginal.
# Requires `sigma ~ 1`. Default remains `marginal = :LA` (adaptive GHQ, 5 nodes per group, on an ordinary `(1 | g)`).
fit_va = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), Beta();
data = dat, marginal = :VA)
# TMB-convention Laplace on an ordinary `(1 | g)`, as native drmTMB fits it.
# Requires `sigma ~ 1`; ML only.
fit_lap = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), Beta();
data = dat, marginal = :Laplace)DRModels.BetaBinomial Type
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; both routes require constant σ (overdispersion). One ordinary mean random intercept or one independent numeric mean slope also supports a predictor-dependent SD with sd(g) ~ predictors, provided sigma ~ 1; the SD predictors must be constant within groups. For a random slope, the SD is the coefficient SD and the observation-level slope predictor multiplies that coefficient in the mean linear predictor. These routes use adaptive Gauss–Hermite quadrature and do not add a free Bernoulli-style residual scale.
One ordinary random intercept on sigma also supports a predictor-dependent SD with sd(group, sigma) ~ predictors. Here the second selector names the distributional parameter. The mean formula must be fixed-effect only, and the SD predictors must be constant within groups. This models between-group variation in log overdispersion separately from variation in the logit mean.
fit = drm(bf(@formula(cbind(successes, failures) ~ x), @formula(sigma ~ 1)), BetaBinomial(); data = dat)
fit = drm(bf(@formula(cbind(successes, failures) ~ x + (1 | g) + (1 | h)), @formula(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)
fit_sd = drm(bf(@formula(cbind(successes, failures) ~ x + (1 | g)),
@formula(sigma ~ 1), @formula(sd(g) ~ treatment)),
BetaBinomial(); data = dat)
fit_sigma_sd = drm(bf(@formula(cbind(successes, failures) ~ x),
@formula(sigma ~ 1 + (1 | g)),
@formula(sd(g, sigma) ~ treatment)),
BetaBinomial(); data = dat)DRModels.Binomial Type
Binomial(; link = :logit)Binomial response family — successes out of known trials. The default :logit link can be changed to :probit or :cloglog; 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. With inverse link h, the likelihood is Binomial(n, μ) with μ = h(η). Non-logit links are available for fixed effects and ordinary adaptive-quadrature random-effect and SD-regression routes. Structured, crossed-Laplace, variational, and TMB-convention Laplace routes currently require :logit, except for :probit with one phylogenetic SD-regression target. A random intercept (1 | g) on the mean fits a GLMM by per-group ADAPTIVE Gauss–Hermite quadrature (nq nodes on the mode- centred scale, #712/#713); 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. A correlated random intercept + slope (1 + x | g) is fit by the same per-group adaptive quadrature helper, nq nodes per axis (the same scheme as BetaBinomial); it needs within-group variation in x — a slope predictor that is constant within any level of g leaves the slope SD and group-level correlation unidentified, and is refused with an informative error (mirrors drmTMB's drm_validate_q2_slope_variation). An independent random slope (0 + x | g) and marginal = :VA on (1 + x | g) remain out of scope. Independent scalar ordinary effects can each have a group-level SD regression, integrated by per-group adaptive Gauss–Hermite quadrature. For an intercept, write sd(g) ~ z; with an independent slope on the same group, select each coefficient in its own formula, for example sd(g, coef(intercept)) ~ z0 and sd(g, coef(x)) ~ z1. The predictors must be constant within group. The ordinary independent-effect route requires marginal = :LA; provider-specific structured SD support is documented in the family guide. Missing responses are refused. Ordinary independent SD fits use local quadrature state and can use threaded profiling; structured Laplace SD fits use shared warm-start state and require threads = false. predict_parameters(fit, newdata)[:sd] returns the fitted random-effect SD on each row, repeating the value for rows in the same group; it accepts predictor-only new data and checks it against the fitted SD design schema. With se = true it returns the prediction and delta-method SE. In-sample this agrees with marginal_parameters(fit)[:sd].
For a sparse mesh field, a single spatial random intercept can use mesh = (; projection, precision, alignment_ids), where projection maps observations to mesh vertices and precision is the sparse vertex precision matrix. Optional alignment_ids must match data.row_id and are used to reorder projection rows after data reordering. The mesh route fits spatial(1 | g) with sd(g, spatial) ~ z by sparse Laplace approximation and supports seeded marginal simulation. It currently requires complete responses, one spatial intercept and one spatial SD formula; additional random effects and other structured inputs are refused. This mesh implementation is Julia-only; native TMB parity and interval-coverage calibration have not been established.
Separation: detect and warn (fixed-effects-only fits)
A fixed-effects-only Bernoulli/Binomial fit is screened at fit time for (quasi-)complete separation: the maximum-likelihood estimate does not exist when some direction d makes the signed design X̃ d ≥ 0 on every row and > 0 on at least one (Albert & Anderson 1984; Konis 2007), decided here by a non-negative least-squares (Farkas) solve. The fit is still returned, a warning names the affected coefficients, their standard errors are Inf (their vcov rows/columns are NaN with an Inf diagonal) and their Wald confint / coeftable intervals are (-Inf, Inf). A coefficient is flagged when ANY separating direction moves it (which can list more coefficients than R's detectseparation, which reports one direction); an aliased column that is a multiple of a flagged one is named too. The check is invariant to column units and robust to a single far outlier. When it proves there is no separation the fit is left alone, even if a Wald |z| is tiny. Only when it cannot decide (iteration budget exhausted) is a scale-free near-separation rule consulted: a fitted probability within 1e-8 of 0/1, a coefficient whose own contribution |βⱼ|·range(xⱼ) spans that whole probability range (> 36.84 on the logit scale), and Wald |z| < 0.05; a failed (non-finite) SE is never read as near separation. If that rule does not fire, it warns that the check was inconclusive rather than reporting no separation. Bootstrap replicate refits do not repeat the warning. No refusal and no penalised estimator by default. drmTMB runs the same check with the same constants, so the two packages flag the same coefficients (the coefficient values themselves are arbitrary stopping points and are not expected to agree). Random-intercept routes are not screened.
Crossed intercepts: Laplace is biased low for Bernoulli data
The default crossed fit is the Laplace approximation drmTMB and lme4 use. For Bernoulli / small-n Binomial data with a large random-intercept SD and few observations per group it underestimates the log-likelihood (≈9.4 nat on a 300 × 4 Bernoulli design with σ_g = 2.5) and shrinks σ. When one grouping has at most 8 levels, marginal = :AGHQ integrates the same model by nested adaptive Gauss–Hermite quadrature and corrects this (#761). Its loglik is not comparable with a Laplace fit's, so lrtest refuses mixed marginals.
Note
DRModels.Binomial shadows Distributions.Binomial; if you need the distribution too (e.g. to simulate), qualify it as Distributions.Binomial.
fit = drm(bf(@formula(cbind(successes, failures) ~ x)), Binomial(); data = dat) # logistic regression
fit = drm(bf(@formula(y ~ x + (1 | g))), Binomial(); data = dat) # 0/1 logistic GLMM
fit = drm(bf(@formula(cbind(successes, failures) ~ x + (1 | g) + (1 | h))), Binomial(); data = dat)
fit_q = drm(bf(@formula(cbind(successes, failures) ~ x + (1 | g) + (1 | h))), Binomial();
data = dat, marginal = :AGHQ) # accurate crossed integral, h ≤ 8 levels
fit_phy = drm(bf(@formula(cbind(successes, failures) ~ x + phylo(1 | species))),
Binomial(); data = dat, tree = tr, se = false)
fit_phy_probit_sd = drm(
bf(@formula(cbind(successes, failures) ~ x + phylo(1 | species)),
@formula(sd(species, phylogenetic) ~ z)),
Binomial(link = :probit); data = dat, tree = tr, se = false)
fit_sd = drm(bf(@formula(y ~ x + (1 | g)), @formula(sd(g) ~ z)),
Binomial(); data = dat)
predict_parameters(fit_sd, dat)[:sd] # one fitted random-effect SD per row
fit_relmat_sd = drm(bf(@formula(y ~ x + relmat(1 | id)),
@formula(sd(id, relmat) ~ z)),
Binomial(); data = dat, K = K)
fitted(fit) # fitted success probabilities on the response scale
# Experimental: Binomial random-intercept variational (ELBO) marginal.
fit_va = drm(bf(@formula(y ~ x + (1 | g))), Binomial(); data = dat, marginal = :VA)
# TMB-convention Laplace on an ordinary `(1 | g)`, as native drmTMB fits it (ML only).
fit_lap = drm(bf(@formula(y ~ x + (1 | g))), Binomial(); data = dat, marginal = :Laplace)DRModels.Gamma Type
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. One ordinary random intercept or independent numeric slope on sigma can have a predictor-dependent SD, while mu stays fixed-effect only. For a slope, use sigma ~ 1 + (0 + x | g) and sd(g, sigma, coef(x)) ~ z; the SD describes the random coefficient on log(sigma), and each row's x multiplies that effect. Here sigma is the Gamma coefficient of variation. These routes use per-group adaptive Gauss-Hermite quadrature and require complete responses and group-constant SD predictors. sigma(fit) includes the fitted conditional group modes; predict_parameters reports the population-level sigma predictor alongside the separate SD prediction.
Note
DRModels.Gamma shadows Distributions.Gamma; qualify the latter if you need it.
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), Gamma(); data = dat)
fit = drm(bf(@formula(y ~ x + (1 | g) + (1 | h)), @formula(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 α
fit_sigma_sd = drm(bf(@formula(y ~ x), @formula(sigma ~ 1 + (1 | g)), @formula(sd(g, sigma) ~ z)),
Gamma(); data = dat)
fit_sigma_slope_sd = drm(
bf(@formula(y ~ x), @formula(sigma ~ 1 + (0 + x | g)), @formula(sd(g, sigma, coef(x)) ~ z)),
Gamma(); data = dat)
# Experimental: Gamma random-intercept variational (ELBO) marginal.
# Requires `sigma ~ 1`. Default remains `marginal = :LA` (adaptive GHQ, 5 nodes per group, on an ordinary `(1 | g)`).
fit_va = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), Gamma();
data = dat, marginal = :VA)
# TMB-convention Laplace on an ordinary `(1 | g)`, as native drmTMB fits it.
# Requires `sigma ~ 1`; ML only.
fit_lap = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), Gamma();
data = dat, marginal = :Laplace)
# Independent random intercept and slope on one factor, `(1 | g) + (0 + x | g)` (#960,
# what drmTMB fits; Laplace). A correlated `(1 + x | g)` is DRModels.jl's own GHQ route
# (drmTMB refuses it for Gamma). `sigma` may carry covariates here.
fit_ind = drm(bf(@formula(y ~ x + (1 | g) + (0 + x | g)), @formula(sigma ~ 1)), Gamma(); data = dat)
exp.(coef(fit_ind, :resd)) # SDs of the intercept and slope effectsDRModels.LogNormal Type
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
DRModels.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]). Crossed/multiple random intercepts, (1 | g) + (1 | h) (#736), delegate the same way as the structured markers below: log(y) is exactly Gaussian, so the fit runs WHOLESALE through drm(f, Gaussian(); data = data-with-logged-response) (_fit_multi_ranef_gaussian, the closed-form exact crossed marginal, not a Laplace approximation), with the reported log-likelihood shifted by -sum(log y).
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.
fit = drm(bf(@formula(y ~ x), @formula(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)DRModels.ZeroOneBeta Type
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. predict / predict_parameters(..., type = :response)[:mu] return the interior beta mean logistic(Xβ̂) (the same quantity drmTMB returns for the mu dpar), not the mixture mean.
A random intercept on the mean, (1 | g), is admitted the same way Beta() admits it: b_g ~ N(0, σ_b²) on the logit mean is integrated out per group by 32-node Gauss–Hermite quadrature (mirrors drmTMB's ordinary RI on zero_one_beta()'s mu, #723); sigma, zoi, coi stay fixed-effects-only.
A phylogenetic/relatedness random intercept on the mean, phylo(1 | species) / relmat(1 | id) (with K = C) / animal(1 | id) (A = C) / spatial(1 | id) (K = C), uses the sparse-Laplace mean-phylo engine (#739): sigma, zoi ~ 1, coi ~ 1 (intercept-only) stay fixed effects; zoi/coi are estimated by their closed-form Bernoulli MLEs (the atom/coi mixture is completely separable from the mean/dispersion/phylo likelihood — see _zeroonebeta_laplace_setup), and the Beta interior likelihood on the mean uses the same verified :beta_fixed-derived phylo/relmat spine as Beta().
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(zoi ~ 1), @formula(coi ~ 1)), ZeroOneBeta(); data = dat)
fit = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1), @formula(zoi ~ 1), @formula(coi ~ 1)), ZeroOneBeta(); data = dat)
fit_phy = drm(bf(@formula(y ~ x + phylo(1 | species)), @formula(sigma ~ 1),
@formula(zoi ~ 1), @formula(coi ~ 1)),
ZeroOneBeta(); data = dat, tree = tr, se = false)DRModels.Tweedie Type
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 (per-group adaptive Gauss–Hermite marginal); sigma ~ 1 and nu ~ 1 remain fixed-effect sub-models on either route. Crossed/multiple random intercepts, (1 | g) + (1 | h), are also supported, via the sparse augmented-state Laplace GLMM engine shared with Gamma/NB2/Beta's crossed routes; sigma ~ 1 and nu ~ 1 are REQUIRED (not just default) on the crossed route, and p is found by an outer 1-D profile search over the one dispersion-like nuisance slot the shared engine does not also give to p. The CORRELATED random slope (1 + x | g), random effects on sigma/nu, and phylogenetic, animal, or spatial markers on mu are not implemented. The related-matrix intercept is supported only on the predictor-dependent SD route described below; this does not establish general structured-effect support.
Predictor-dependent SD regression is available for one ordinary mean random intercept or one independent numeric mean slope with a group-constant SD design, for example sd(group) ~ z or sd(group, coef(x)) ~ z. A single related-matrix mean intercept is also available with relmat(1 | group), sd(group, relmat) ~ z, and K = K. For ordinary effects, the mean, dispersion, and power retain their existing submodels. The related-matrix route currently requires sigma ~ 1 and nu ~ 1. Ordinary effects use adaptive Gauss–Hermite quadrature; the related-matrix route uses a Laplace approximation. Correlated ordinary blocks and other structured providers are not implemented, and interval coverage is not calibrated.
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(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(@formula(y ~ x + (1 | g)), @formula(nu ~ 1)), Tweedie(); data = dat)
exp(coef(fit_re, :resd)[1]) # random-intercept SD
fit_slope = drm(bf(@formula(y ~ x + (0 + x | g)), @formula(nu ~ 1)), Tweedie(); data = dat)
exp(coef(fit_slope, :resd)[1]) # random-slope SD
fit_crossed = drm(bf(@formula(y ~ x + (1 | g) + (1 | h))), Tweedie(); data = dat)
re_sd(fit_crossed) # Dict(:g => ..., :h => ...)DRModels.CumulativeLogit Type
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 per-group adaptive Gauss–Hermite quadrature (41 nodes centred on each group's conditional mode); 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). Crossed/multiple random intercepts, (1 | g) + (1 | h), are also supported, via the same sparse-Laplace engine Poisson/Binomial's crossed routes reuse; the cutpoints are found by an outer profile search since that engine has no cutpoint slot. Correlated slopes (1 + x | g) and relmat/animal/spatial structured markers are not implemented.
fit = drm(bf(@formula(y ~ x)), CumulativeLogit(); data = dat) # y coded 1..K
fit_re = drm(bf(@formula(y ~ x + (1 | g))), CumulativeLogit(); data = dat)
fit_slope = drm(bf(@formula(y ~ x + (0 + x | g))), CumulativeLogit(); data = dat)
fit_phylo = drm(bf(@formula(y ~ x + phylo(1 | species))), CumulativeLogit(); data = dat, tree = tr)
fit_crossed = drm(bf(@formula(y ~ x + (1 | g) + (1 | h))), CumulativeLogit(); data = dat)
re_sd(fit_crossed) # Dict(:g => ..., :h => ...)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:
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(nu ~ 1)),
SkewNormal(); data = dat)DRModels.SkewNormal Type
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.
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(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)
predict_parameters(fit, dat)[:nu] # same slant (identity; not exp(η))
# Predictor-dependent SD for an ordinary random intercept
fit_sd = drm(bf(@formula(y ~ x + (1 | site)), @formula(sigma ~ 1),
@formula(nu ~ 1), @formula(sd(site) ~ habitat)),
SkewNormal(); data = dat)
predict_parameters(fit_sd, dat)[:sd] # fitted random-intercept SD by site
# Predictor-dependent SD for one independent mean slope
fit_slope_sd = drm(bf(@formula(y ~ x + (0 + x | site)), @formula(sigma ~ 1),
@formula(nu ~ 1), @formula(sd(site, coef(x)) ~ habitat)),
SkewNormal(); data = dat, nAGQ = 9)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. The information matrix is singular at ν = 0, so the Wald SE for nu is unreliable there. check_drm reports skew_normal_nu_wald and warns when any fitted |ν̂| < 0.5. This cutoff is an uncalibrated heuristic of this package; no Wald test is run (#1060). Predictor-dependent SD supports one ordinary mean intercept or independent numeric slope and uses one-dimensional adaptive Gauss–Hermite quadrature. It requires fixed sigma ~ 1, nu ~ 1, complete responses, and a full-rank group-level SD design. Correlated and structured random effects are not supported on this route; interval coverage calibration is deferred.
Two-column response
DRModels.cbind Function
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, not yet stable. API and numerics may change. The bounded R bridge workflow is shown below.
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 not supported here. The same narrow model is available experimentally through drmTMB(..., engine = "julia"); it does not yet reproduce every fitted result and post-fit method from drmTMB's default R/TMB engine.
using DRModels, 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.29561248510280674, 0.37501846625780844, 0.7485514011793476], [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.24247894915511678, 0.2428178772879237, 0.5312766777084319], std_error = Union{Missing, Float64}[0.09074205154086845, 0.08920306612733003, 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 supported.
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"). Dropping responses before fitting behaves differently from drmTMB's default R/TMB engine. Profile and bootstrap intervals are not available here. The Gaussian predictor-SD Wald interval is a natural-scale delta interval, can cross zero, and has not been shown to have reliable coverage across datasets or sample sizes.
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.
using DRModels, 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.4135447546956656, -0.46248422948769086, 0.41618947536448614, 0.06592534270473177, -0.03289463233738368, -0.18891817752408835, 0.01648601894117721], 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").
DRModels.mi Function
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.
DRModels.miss_control Function
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.
DRModels.impute_model Function
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.
DRModels.imputed Function
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.
sourceConditional imputation for one explicitly selected two-predictor marker variable.
sourceConditional finite-state imputation for the marker variable in a formula-based fit.
sourceDRModels.JointTwoDrmFit Type
JointTwoDrmFitFormula-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>).
DRModels.JointMissingControl Type
Validated controls for the joint missing-predictor formula admission.
sourceDRModels.JointImputeModel Type
Formula and family for the one predictor distribution in a joint fit.
source