What can I fit today?
Use this page when you know what you measured but are deciding what kind of model to fit. DRModels.jl is for questions where predictors may change not only the average response, but also its variability, its zeros, or the way two responses vary together.
Start with the question closest to yours:
| Your response or question | A good first page |
|---|---|
| A continuous measurement, such as a trait, concentration, or score | Fit your first model |
| Counts, including overdispersed counts or extra zeros | Count abundance and extra zeros |
| Proportions, success rates, or binomial data | Proportions and success rates |
| Whether predictors change both the average and the spread | When variance carries signal |
| Two responses that may be associated | Changing residual coupling |
| Observations linked by a phylogeny, space, pedigree, or relatedness matrix | Biological examples |
The first examples are motivated by ecology, evolution, and environmental science, but the models are also useful wherever a response can vary in more than its average. Read on for the model choices, or use the which page next table to move directly to a guide.
The modelling idea: a formula per parameter
DRModels.jl is distributional regression — it lets you model more than the mean of a response. Each parameter supported by the chosen response family can have its own formula, bundled together with bf:
| Parameter | What it controls | Formula |
|---|---|---|
| μ (mean / location) | where the response sits | the response formula, y ~ … |
sigma (scale / dispersion) | how spread out it is | sigma ~ … |
family extras — nu, zi, hu, zoi, coi | shape, zero-inflation, hurdle, boundary mass | nu ~ …, zi ~ …, … |
rho12 (bivariate only) | residual correlation between two responses | rho12 ~ … |
Always write sigma (never tau) for the scale and rho12 for residual correlation. The same machinery drives every parameter: a formula → a linear predictor → a family-specific link. Which parameters are available depends on the family (Gaussian has no nu; only counts take zi / hu).
See Which scale are you modelling? for the difference between the residual sigma, a group-level SD, and a known sampling variance — they are distinct quantities DRModels.jl keeps separate.
The bf(...) front end
bf (alias drm_formula) collects one formula per parameter, following the same formula-per-parameter pattern as drmTMB and brms. It has two forms:
Univariate — positional. The first formula is the mean; the rest name their parameter on the left-hand side:
bf(@formula(y ~ x), @formula(sigma ~ x)) # mean + scale
bf(@formula(y ~ x), @formula(sigma ~ 1), @formula(nu ~ 1)) # + shape (Student)Bivariate — keyword. Name each of the two responses' parameters explicitly, including the residual correlation rho12:
bf(mu1 = @formula(y1 ~ x), mu2 = @formula(y2 ~ x),
sigma1 = @formula(sigma1 ~ 1), sigma2 = @formula(sigma2 ~ 1),
rho12 = @formula(rho12 ~ x))You pass the bundle and a family to drm: drm(bf(...), Gaussian(); data = dat). ML is the default (REML likelihoods are not comparable across fixed-effect structures, so ML is what model selection needs). For the full grammar — including cbind(successes, failures) ~ … for binomial-type responses — see the model specification reference.
Supported families
Pick the family from the shape of the response; the Choosing response families guide has the full decision table and worked examples. DRModels.jl implements drmTMB's complete family set:
| Family | Response | Mean link | sigma slot / extras |
|---|---|---|---|
Gaussian | real-valued, symmetric | identity | residual SD σ (log) |
Student | real-valued, heavy tails | identity | scale σ + d.o.f. nu |
LogNormal | positive, multiplicative | identity on log y | SD of log y |
Gamma | positive, continuous | log | CV → shape α = 1/σ² |
Tweedie | positive with exact zeros | log | √dispersion σ + power nu ∈ (1,2) |
Poisson | counts (var ≈ mean) | log | — (+ zi / + hu) |
NegBinomial2 | overdispersed counts (NB2) | log | dispersion θ (+ zi / + hu) |
TruncatedNegBinomial2 | positive counts (≥ 1) | log | dispersion θ |
Beta | proportions in (0,1) | logit | precision φ = 1/σ² |
BetaBinomial | successes / trials, overdispersed | logit | φ = 1/σ² (cbind) |
Binomial | successes / trials | logit | — (cbind or 0/1) |
ZeroOneBeta | proportions on [0,1] incl. 0 and 1 | logit | φ + zoi / coi |
CumulativeLogit | ordered categories 1..K | logit (cutpoints) | K−1 ordered cutpoints |
The count modifiers — zi (zero-inflation, ZIP / ZINB), hu (hurdle), and the beta boundary modifiers zoi / coi — are themselves formulas you add to the bundle, e.g. bf(@formula(y ~ x), @formula(zi ~ 1)).
Structured and random effects
On top of fixed effects, DRModels.jl carries ordinary random effects and several structured effects whose covariance comes from a known matrix or geometry. Write them as terms in the mean formula:
| Effect | Term | What it encodes |
|---|---|---|
| Random intercept / slope | (1 | g), (0 + x | g), (1 + x | g) | exchangeable group variation; correlated slopes |
| Crossed / nested RE | (1 | g) + (1 | h) | multiple grouping factors |
phylo | phylo(1 | species) | covariance from a phylogenetic tree |
spatial | spatial(1 | site) | exponential kernel over coordinates, estimated range |
animal | animal(1 | id) | additive-genetic covariance from a pedigree A |
relmat | relmat(1 | id) | a user-supplied relatedness matrix K |
meta_V | meta_V(v) | known sampling (co)variances for meta-analysis |
Which families route which effects. Gaussian gets the full set — random intercept/slope, correlated and crossed RE, phylo / spatial / animal / relmat structure, meta_V, and a random effect on sigma — fit in closed form or via a sparse augmented-state Laplace approximation. The non-Gaussian families (Poisson, NB2, Binomial, Gamma, Beta, Student-t, LogNormal, Beta-binomial, Tweedie and CumulativeLogit) carry random intercepts via Gauss–Hermite marginals, and random slopes for every one of them except Binomial(), which supports (1 | g) on the mean only. The capability matrix records which slope form — independent or correlated — each family admits. Phylogenetic (phylo) effects go via a sparse Laplace path for Poisson, NB2, Binomial, Gamma, Beta, Beta-binomial and CumulativeLogit. On that path NB2, Gamma and Beta accept a covariate dispersion formula sigma ~ x (a per-observation log σ), while BetaBinomial() requires a constant sigma and Binomial() carries no dispersion parameter at all; Student() rejects meta_V and every structured marker. LogNormal() is the exception: because log y is exactly Gaussian, its phylo/relmat structured markers on the mean delegate WHOLESALE to Gaussian() on log y (exact, not a Laplace approximation) rather than the shared non-Gaussian sparse path; animal/spatial are not implemented for LogNormal(). So put predictors on sigma for Gaussian freely. On the non-Gaussian phylogenetic path, follow the family-specific boundary: NB2, Gamma, and Beta accept sigma ~ x; BetaBinomial() requires constant dispersion; and Binomial() has no dispersion formula.
CumulativeLogit (ordinal) carries an ordinary random intercept (1 | g) or an independent random slope (0 + x | g) on mu via the same Gauss–Hermite scheme, and an intercept-only phylo(1 | species) through the same sparse Laplace engine; the correlated form (1 + x | g) and the other structured effects (relmat / animal / spatial) are not implemented yet.
For the verified engine behind the phylogenetic models — the q=4 phylogenetic bivariate location–scale model, which matches drmTMB's fit and still returns usable Wald and bootstrap intervals where drmTMB's Hessian is singular — see Large data for the scoped performance evidence.
After the fit
The common starting points are coef, summary / coeftable, fitted, response residuals, and aic / bic for appropriate ML comparisons. The rest depends on the fitted route: standard errors and Wald intervals require an available covariance estimate; profile intervals require a stored or precomputed profile target; and bootstrap, prediction, simulation, and variance-component summaries apply only where that model implements them. Use profile_targets(fit) to see which profile intervals are ready for a particular fit, and consult the capability matrix before planning an analysis around a post-fit method. re_sd and vc describe mixed or structured models, not every fit.
See Checking and using fitted models for a worked Gaussian workflow, and Which scale are you modelling? for honest inference at a variance boundary.
Which page next
| If you want to… | Go to |
|---|---|
| Fit your first model, end to end | Get started |
| Choose the right response family | Choosing response families |
Tell residual σ, group SD, and known V apart | Which scale are you modelling? |
| Extract coefficients, CIs, predictions | Checking and using fitted models |
| Diagnose convergence | Convergence |
| Scale to large data | Large data |
| Choose the marginal method (Laplace vs VA) | Marginal: LA vs VA |
| Model variability as signal (location–scale) | When variance carries signal |
| Predictors on a random effect's SD (location–scale–scale, incl. climate-dependent phylogenetic SD) | Part 2: location–scale–scale |
Change residual coupling with rho12 | Changing residual coupling with rho12 |
| Robust continuous responses | Robust continuous responses |
| Counts and extra zeros | Count abundance and extra zeros |
| Proportions and success rates | Proportions and success rates |
| Phylogenetic / spatial / animal models | Phylogenetic · Spatial · Animal |
| Meta-analysis with known variances | Meta-analysis |
| The full API reference | Model specification · Fitting & post-fit |