DRModels.jl

Appearance

Response families ​

Status — Stable

Every family below is fitted by the verified drm front end, and every worked snippet on this page runs. This is a reference catalogue with one minimal fit per family; for the modelling rationale and side-by-side R parity see Choosing response families and the per-topic tutorials.

You pick a family from the shape of your response, then give each of its parameters a formula with bf. The mean μ always comes from the response formula (y ~ …); the family's second parameter lives in the sigma slot (sigma ~ …) with a family-specific link, and some families add further parameters (nu, zoi, coi, ordinal cutpoints).

Response looks like…FamilyMean linksigma slot
Real-valued, symmetricGaussian()identityresidual SD σ (log)
Real-valued, heavy tailsStudent()identityscale σ (log) + d.o.f. ν = 2 + exp(η) (not log ν)
Strictly positive, multiplicativeLogNormal()identity on log ySD of log y (log)
Strictly positive, right-skewedGamma()logshape via α = 1/σ²
Positive with exact zerosTweedie()log√dispersion σ + power ν ∈ (1,2)
Counts (variance ≈ mean)Poisson()log—
Counts, overdispersedNegBinomial2()logoverdispersion scale σ (log); size θ = 1/σ²
Positive counts (zero-truncated)TruncatedNegBinomial2()logoverdispersion scale σ (log); size θ = 1/σ²
Positive counts, no overdispersionTruncatedPoisson()log—
Proportions in (0,1)Beta()logitprecision via φ = 1/σ²
Successes out of trialsBinomial()logit— (cbind(s,f) or 0/1)
Successes, overdispersedBetaBinomial()logitprecision via φ = 1/σ²
Proportions on [0,1] incl. 0/1 (needs both exact 0s and 1s for coi)ZeroOneBeta()logitprecision φ + zoi / coi
Ordered categories 1..KCumulativeLogit()logit cutpoints— (K−1 cutpoints)

A handful of families share the name of a Distributions.jl distribution (Gaussian excepted: Poisson, Binomial, Beta, Gamma, LogNormal). DRModels exports its own family object of that name; when a snippet also needs the distribution to simulate, it qualifies it as Distributions.Poisson etc.

Continuous, real-valued ​

Gaussian — Gaussian() ​

Identity link on the mean, log link on the residual SD σ. The textbook location–scale model: both the mean and the spread can depend on covariates.

The code below also loads the Distributions package to simulate data. If you have not added it to your Julia environment, run using Pkg; Pkg.add("Distributions") once first.

julia
using DRModels, Random
import Distributions          # qualified, for the simulating distributions
Random.seed!(1)

n = 600
x = randn(n)
y = (1.0 .+ 0.5 .* x) .+ exp.(-0.3 .+ 0.4 .* x) .* randn(n)
fit = drm(bf(@formula(y ~ x), @formula(sigma ~ x)), Gaussian(); data = (; y, x))
(mu = coef(fit, :mu), logsigma = coef(fit, :sigma))   # ≈ (1.0, 0.5) and (-0.3, 0.4)
(mu = [1.0012462425015338, 0.5235831241758869], logsigma = [-0.3630675863006263, 0.3958446953329363])

Student-t — Student() ​

The robust sibling of Gaussian: identity μ, log σ, plus degrees of freedom ν that govern the tail weight, on the link ν = 2 + exp(η) (so ν > 2). Small ν downweights outliers; ν → ∞ returns to Gaussian. Recover ν with 2 + exp(coef(fit, :nu)), never exp(...) alone.

julia
Random.seed!(2)
yt = 0.5 .+ 0.8 .* x .+ 2.0 .* rand(Distributions.TDist(4), n)
fitt = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(nu ~ 1)), Student(); data = (; y = yt, x))
(slope = coef(fitt, :mu)[2], nu = 2 + exp(coef(fitt, :nu)[1]))   # d.o.f., ν = 2 + exp(η)
(slope = 0.6711518143080886, nu = 5.16048805149628)

LogNormal — LogNormal() ​

For strictly positive, multiplicative data. The μ formula is the mean of log y; σ (log link) is the SD of log y. The response-scale median is exp(μ). A random intercept (1 | g) or correlated slope (1 + x | g), or a phylo(1 | group) (needs tree = ...) / relmat(1 | group) (needs K = ...) structured marker, may be placed on the mean; the structured routes delegate to Gaussian() on log y (exact — no approximation) and shift the reported log-likelihood by the parameter-free Jacobian. animal/ spatial markers are not implemented for LogNormal().

julia
Random.seed!(3)
yln = exp.(0.5 .+ 0.3 .* x .+ 0.4 .* randn(n))
fitln = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), LogNormal(); data = (; y = yln, x))
(intercept = coef(fitln, :mu)[1], sd_logy = exp(coef(fitln, :sigma)[1]))   # ≈ (0.5, 0.4)
(intercept = 0.5051299105689142, sd_logy = 0.39721037135420106)

Gamma — Gamma() ​

Strictly positive, right-skewed (durations, sizes, concentrations). Log link on the mean; the sigma slot maps to the shape α = 1/σ², so the coefficient of variation is 1/√α. Recover the shape as exp(-2·log σ).

julia
Random.seed!(4)
αsh = 8.0
μg  = exp.(0.5 .+ 0.4 .* x)
yg  = Float64.([rand(Distributions.Gamma(αsh, μi / αsh)) for μi in μg])
fitg = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), Gamma(); data = (; y = yg, x))
(mean_ratio = exp(coef(fitg, :mu)[2]), shape = exp(-2 * coef(fitg, :sigma)[1]))   # shape ≈ 8
(mean_ratio = 1.4766839765751336, shape = 7.712957344355953)

Tweedie — Tweedie() ​

Positive data with exact zeros (biomass, rainfall, total loss): a point mass at 0 plus a positive continuous part. Log link on the mean; sigma is the √dispersion (φ = σ²); nu is the power p ∈ (1,2) on a logit link (p = 1 + logistic(coef(:nu))).

julia
Random.seed!(5)
function rtw(μi, φ, p)                              # compound Poisson–Gamma draw
    λ = μi^(2 - p) / (φ * (2 - p)); γ = φ * (p - 1) * μi^(p - 1)
    N = rand(Distributions.Poisson(λ))
    N == 0 ? 0.0 : rand(Distributions.Gamma(N * (2 - p) / (p - 1), γ))
end
ytw = [rtw(exp(0.5 + 0.3 * xi), 2.0, 1.5) for xi in x]
fittw = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(nu ~ 1)), Tweedie(); data = (; y = ytw, x))
(zeros = count(==(0.0), ytw), power = 1 + 1 / (1 + exp(-coef(fittw, :nu)[1])))   # power ≈ 1.5
(zeros = 186, power = 1.498574838587078)

An ordinary (1 | g) random intercept, and an independent (0 + x | g) random slope, on mu are both supported. Each group's effect is integrated out by per-group adaptive Gauss–Hermite quadrature with 21 nodes; drmTMB fits the same tweedie() mu ~ x + (1 | g) / (0 + x | g) models with a one-point Laplace approximation, so log-likelihoods can differ slightly between the two packages (see Laplace vs variational marginals). The group SD is exp(coef(fit, :resd)[1]) either way. The CORRELATED random slope (1 + x | g), random effects on sigma/nu, and structured (phylo/relmat) markers on mu remain unimplemented (drmTMB rejects the same forms).

Counts ​

Poisson — Poisson() ​

Counts with variance ≈ mean. Log link on the mean; no scale parameter. Slopes are log rate ratios (exp(β) is the multiplicative effect on the expected count).

julia
Random.seed!(6)
λp = exp.(0.3 .+ 0.5 .* x)
yp = Float64.([rand(Distributions.Poisson(λi)) for λi in λp])
fitp = drm(bf(@formula(y ~ x)), Poisson(); data = (; y = yp, x))
(logmu = coef(fitp, :mu), rate_ratio = exp(coef(fitp, :mu)[2]))   # rate ratio ≈ 1.65
(logmu = [0.28367487517333523, 0.518705073083966], rate_ratio = 1.6798509564552755)

Negative-binomial — NegBinomial2() ​

Overdispersed counts. Log link on the mean; the sigma slot carries log σ (not log θ), with size θ = 1/σ²; variance = μ + μ²/θ, so θ → ∞ recovers Poisson.

julia
Random.seed!(7)
θnb = 2.5
μnb = exp.(0.4 .+ 0.5 .* x)
ynb = Float64.([rand(Distributions.NegativeBinomial(θnb, θnb / (θnb + μi))) for μi in μnb])
fitnb = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), NegBinomial2(); data = (; y = ynb, x))
exp(-2 * coef(fitnb, :sigma)[1])      # estimated dispersion θ = 1/σ² ≈ 2.5
2.4831294961225976

The sigma coefficient is on the log-σ scale, and the NB2 dispersion is θ = 1/σ² = exp(-2·coef(:sigma)) (not exp(coef(:sigma))).

No underdispersion

NB2 cannot represent dispersion below the Poisson limit: its variance is never below the mean. On underdispersed counts the size θ goes to infinity (σ → 0) 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 does not offer a Conway–Maxwell–Poisson or generalized-Poisson family (#976).

TruncatedNegBinomial2() is the same family conditioned on y ≥ 1, for counts that are positive by construction (litter sizes, group sizes given presence). TruncatedPoisson() is the Poisson conditioned on y ≥ 1; adding an hu ~ ... part to either truncated family (or to Poisson() / NegBinomial2()) fits the hurdle model.

Differs from drmTMB: hu on Poisson()

DRModels.jl accepts hu on a plain Poisson() and fits the hurdle Poisson. drmTMB currently refuses it: poisson() supports only mu and an optional zi, and hu raises "Poisson models only support mu and optional zi". There is therefore no drmTMB spelling for this model, and no R-parity comparison is available for it (twin track, #726). For a hurdle count model that exists on both sides, use NegBinomial2() + hu, which agrees numerically across the two packages. Whether the two should align by drmTMB adding hu to poisson() or DRModels.jl refusing it is an open design decision; this note changes no behaviour.

Proportions and binary data ​

Binomial — Binomial() ​

Successes out of known trials (logistic regression). Logit link on the success probability; no dispersion parameter. The response is either a cbind(successes, failures) two-column term or a plain 0/1 Bernoulli vector.

julia
Random.seed!(8)
ntrials = fill(20, n)
pb = 1 ./ (1 .+ exp.(-(0.2 .+ 0.7 .* x)))
succ = Float64.([rand(Distributions.Binomial(ntrials[i], pb[i])) for i in 1:n])
fail = ntrials .- succ
fitb = drm(bf(@formula(cbind(succ, fail) ~ x)), Binomial(); data = (; succ, fail, x))
coef(fitb, :mu)      # logit-scale intercept & slope ≈ (0.2, 0.7)
2-element Vector{Float64}:
 0.18193662500113625
 0.7023511831871687

Beta — Beta() ​

Proportions strictly inside (0,1). Logit link on the mean; the sigma slot carries σ with precision φ = 1/σ² (recover precision as exp(-2·log σ)).

julia
Random.seed!(9)
μbe = 1 ./ (1 .+ exp.(-(0.3 .+ 0.6 .* x))); φbe = 12.0
ybe = Float64.([rand(Distributions.Beta(μi * φbe, (1 - μi) * φbe)) for μi in μbe])
fitbe = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), Beta(); data = (; y = ybe, x))
(slope = coef(fitbe, :mu)[2], precision = exp(-2 * coef(fitbe, :sigma)[1]))   # φ ≈ 12
(slope = 0.5645152894705681, precision = 11.356120815443854)

Beta-binomial — BetaBinomial() ​

Successes out of trials with extra-binomial overdispersion. Logit mean; precision φ = 1/σ² in the sigma slot. Requires a cbind(successes, failures) response.

julia
Random.seed!(10)
ntr = fill(30, n)
μbb = 1 ./ (1 .+ exp.(-(0.1 .+ 0.5 .* x))); φbb = 5.0
sbb = Float64.([rand(Distributions.BetaBinomial(ntr[i], μi * φbb, (1 - μi) * φbb)) for (i, μi) in enumerate(μbb)])
fbb = ntr .- sbb
fitbb = drm(bf(@formula(cbind(s, f) ~ x), @formula(sigma ~ 1)), BetaBinomial(); data = (; s = sbb, f = fbb, x))
(slope = coef(fitbb, :mu)[2], precision = exp(-2 * coef(fitbb, :sigma)[1]))   # φ ≈ 5
(slope = 0.5632915727351466, precision = 4.944685720682151)

Differs from drmTMB: structured effects on BetaBinomial()

DRModels.jl admits structured random effects (e.g. tree= phylogenetic terms) on BetaBinomial(). drmTMB refuses structured markers on beta_binomial (planned there), so a phylogenetic beta-binomial fit has no drmTMB comparator and no R-parity claim (twin track, #740). Plain fixed-effect beta-binomial fits are available on both sides.

Zero-one-inflated beta — ZeroOneBeta() ​

Proportions on the closed interval [0,1] (values may be exactly 0 or 1). Beta mean (logit) and precision φ on the interior, plus zoi (logit probability the value is a boundary) and coi (logit probability of a 1 given a boundary).

Data with only zeros, or only ones

coi is estimated from the boundary values alone. If the data contain exact zeros but no exact ones (or ones but no zeros), the best estimate of coi is exactly 0 (or 1): the parameter is at its boundary and not identified. The fit still reports converged = true, and the standard error and interval for coi are meaningless. The mu, sigma and zoi estimates are not affected, because coi ~ 1 enters the likelihood separately; report those and not coi. coi cannot currently be dropped from the model, and an ordered-beta family, the usual alternative for this case, is not available (DRModels.jl#975).

julia
Random.seed!(11)
function rzob(xi)
    μi = 1 / (1 + exp(-(0.2 + 0.5 * xi))); φ = 10.0; zoi = 0.15; coi = 0.4
    rand() < zoi ? (rand() < coi ? 1.0 : 0.0) : rand(Distributions.Beta(μi * φ, (1 - μi) * φ))
end
yzob = [rzob(xi) for xi in x]
fitzob = drm(bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(zoi ~ 1), @formula(coi ~ 1)), ZeroOneBeta(); data = (; y = yzob, x))
(boundaries = count(yi -> yi == 0 || yi == 1, yzob), zoi = 1 / (1 + exp(-coef(fitzob, :zoi)[1])))
(boundaries = 87, zoi = 0.145)

Ordered categories ​

Cumulative-logit — CumulativeLogit() ​

Ordered categories coded 1, 2, …, K (Likert ratings, severity classes): Pr(y ≤ k) = logistic(θ_k − η) with ordered cutpoints. The single linear predictor comes from the mu formula with its intercept dropped (the cutpoints absorb the level), so coef(fit, :mu) is the slope.

julia
Random.seed!(12)
θcut = [-1.0, 0.0, 1.2]; K = 4
yo = map(x) do xi
    u = rand(); η = 0.8 * xi; cat = K
    for j in 1:(K-1)
        (u < 1 / (1 + exp(-(θcut[j] - η)))) && (cat = j; break)
    end
    cat
end
fito = drm(bf(@formula(y ~ x)), CumulativeLogit(); data = (; y = Float64.(yo), x))
coef(fito, :mu)[1]      # slope ≈ 0.8
0.8846602950829421

An ordinary random intercept (1 | g) or an independent random slope (0 + x | g) on mu is integrated out by per-group adaptive Gauss–Hermite quadrature with 41 nodes (the same scheme as the Poisson/Gamma/Tweedie GLMM routes, with more nodes because ordinal group posteriors can be far from normal); the correlated form (1 + x | g) is not implemented.

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

An unlabelled, intercept-only phylogenetic random intercept phylo(1 | species) on mu is fit via the sparse augmented-state Laplace GLMM route (the same engine as the Poisson/Gamma/Binomial/Beta phylo(1 | g) cells), matching drmTMB 0.7.0's validate_ordinal_phylo_mu_structured_term(). It cannot be combined with an ordinary (1 | g)/(0 + x | g) random effect on mu.

julia
fit_phylo = drm(bf(@formula(y ~ x + phylo(1 | species))), CumulativeLogit(); data = dat, tree = tr)
re_sd(fit_phylo)[:species]   # phylogenetic SD, raw branch-length scale

See also ​