Getting started

Status — Stable

A standalone first-fit walkthrough. Everything on this page runs against the verified Gaussian front end (drm / bf) and the post-fit accessors (coef, loglik, confint, summary). For moving between R and Julia, see Coming from R and the Rosetta; for the full capability map see What can I fit today?.

DRModels.jl is distributional regression: instead of a single linear predictor for the mean, a chosen response family can let you give its supported parameters their own formulas. The simplest case puts a formula on the mean μ and a formula on the residual scale σ, so the spread of the data can change with covariates just like the mean does.

This page takes you from a clean Julia session to a fitted model and shows how to read what came back.

Install

This is the direct-Julia route. If Julia is new to you, first install Julia 1.10 or later and open its Julia prompt (the REPL). If you would rather keep working in R, you do not need this setup: use the optional drmTMB Julia-engine route instead.

Install the current package directly from GitHub:

using Pkg
Pkg.add(url = "https://github.com/itchyshin/DRModels.jl")
using DRModels

If you are developing DRModels.jl from a local clone, replace the installation line with Pkg.develop(path = "/absolute/path/to/your/DRModels.jl").

The two verbs you will use the most are exported at the top level:

• bf(...) — bundle one formula per distributional parameter (alias drm_formula), using the familiar drmTMB / brms pattern.

• drm(formula, family; data = ...) — fit the model by maximum likelihood.

Fit your first distributional regression

We simulate data whose mean and spread both depend on a covariate x, then recover that structure. The mean rises with x; the residual standard deviation also rises with x (heteroscedasticity), so a mean-only model would be mis-specified.

using DRModels, Random
Random.seed!(20260610)

n    = 400
x    = randn(n)
μ    = 1.0 .+ 0.5 .* x            # mean increases with x
logσ = -0.4 .+ 0.3 .* x           # log residual SD ALSO increases with x
y    = μ .+ exp.(logσ) .* randn(n)
dat  = (; y, x)

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

Distributional regression fit (Gaussian)
  nobs = 400   logLik = -412.5663   converged = true
  Residual SD (response scale) varies with the σ model: 0.2689 to 1.7955   dof_residual = 396

Mean model (μ):
  H0: coefficient = 0
                Coef.  Std.Error        z    Pr(>|z|)
  (Intercept)  1.0445     0.0358  29.1985  2.027e-187
  x            0.4995     0.0308  16.2226   3.492e-59

Scale model (log σ):
  H0: coefficient = 0 on log σ
  intercept ⇔ σ = 1 (unit-dependent); slope ⇔ equal dispersion
                 Coef.  Std.Error         z   Pr(>|z|)
  (Intercept)  -0.4008     0.0354  -11.3247  9.902e-30
  x             0.3181     0.0363    8.7622  1.915e-18

bf(@formula(y ~ x), @formula(sigma ~ x)) is the whole idea in one line: the first formula sets the response and the μ sub-model; the second sets the σ sub-model. (Drop the sigma formula and it defaults to sigma ~ 1, a constant scale — ordinary homoscedastic regression.)

Read the coefficients

coef(fit, :param) returns the coefficient block for one distributional parameter. The mean block is on the response scale and should recover (1.0, 0.5):

coef(fit, :mu)

2-element Vector{Float64}:
 1.0444504833162647
 0.49949236057252433

The scale block acts on log σ, so it should recover (-0.4, 0.3):

coef(fit, :sigma)

2-element Vector{Float64}:
 -0.400753568936077
  0.3180749930522409

A positive σ-slope means the residual spread grows with x. Because the σ link is logarithmic, exponentiating a slope turns it into a multiplicative effect on the residual SD per unit of x:

exp(coef(fit, :sigma)[2])      # residual-SD ratio per one-unit increase in x

1.3744793337631318

Calling coef(fit) with no parameter concatenates every block into one vector (μ first, then σ) — the raw parameter vector the optimiser returns.

Read the model fit

loglik(fit) is the maximised log-likelihood; aic / bic build on it for model comparison, and nobs / dof report the sample size and number of estimated parameters:

(loglik = loglik(fit), aic = aic(fit), bic = bic(fit), nobs = nobs(fit), dof = dof(fit))

(loglik = -412.56627778015184, aic = 833.1325555603037, bic = 849.0984137487357, nobs = 400, dof = 4)

Always check that the optimiser actually converged before trusting any of the above:

is_converged(fit)

true

summary(fit) prints estimates, standard errors, Wald intervals, and its coefficient-level z and p summaries for every block at once. Treat those summaries as evidence conditional on this fitted model, not as proof that a biological effect is real; check diagnostics and the size of the estimated scale change as well. Row names are prefixed with the parameter (mu: …, sigma: …) so they stay unique:

summary(fit)

────────────────────────────────────────────────────────────────────────────────
                     Estimate  Std.Error       z  Pr(>|z|)  Lower 95%  Upper 95%
────────────────────────────────────────────────────────────────────────────────
mu: (Intercept)      1.04445   0.0357707   29.20    <1e-99   0.974341   1.11456
mu: x                0.499492  0.0307899   16.22    <1e-58   0.439145   0.55984
sigma: (Intercept)  -0.400754  0.0353876  -11.32    <1e-29  -0.470112  -0.331395
sigma: x             0.318075  0.0363008    8.76    <1e-17   0.246927   0.389223
────────────────────────────────────────────────────────────────────────────────

Confidence intervals

confint(fit) returns one row (param, coef, estimate, lower, upper) per coefficient, on each block's working scale (μ on the response scale, σ on log σ). The default is a 95% Wald interval:

confint(fit)

4-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
 (param = :mu, coef = "(Intercept)", estimate = 1.0444504833162647, lower = 0.9743412045976196, upper = 1.1145597620349097)
 (param = :mu, coef = "x", estimate = 0.49949236057252433, lower = 0.43914519895183113, upper = 0.5598395221932175)
 (param = :sigma, coef = "(Intercept)", estimate = -0.400753568936077, lower = -0.4701119295150904, upper = -0.3313952083570636)
 (param = :sigma, coef = "x", estimate = 0.3180749930522409, lower = 0.24692670465969427, upper = 0.38922328144478757)

For a coefficient near a boundary, or whenever you want intervals that do not assume a quadratic log-likelihood, ask for the profile-likelihood interval instead. It re-optimises the other parameters at each fixed value and uses a likelihood-ratio calibration. That calibration is asymptotic, so a profile interval does not have universal exact coverage:

confint(fit; method = :profile)

4-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
 (param = :mu, coef = "(Intercept)", estimate = 1.0444504833162647, lower = 0.97413231982766, upper = 1.1147107123562001)
 (param = :mu, coef = "x", estimate = 0.49949236057252433, lower = 0.43905965187659424, upper = 0.5601510661608071)
 (param = :sigma, coef = "(Intercept)", estimate = -0.400753568936077, lower = -0.4685422102581028, upper = -0.32975066456783414)
 (param = :sigma, coef = "x", estimate = 0.3180749930522409, lower = 0.24711546066516699, upper = 0.38945046057662325)

Beyond a first fit

The same front end also provides these next steps; the capability map records their current boundaries:

• Per-parameter prediction — predict_parameters (fitted μ/σ/… on new data), marginal_parameters (population-averaged), and prediction_grid for building a swept newdata grid from a reference table.

• Auditable profile-likelihood CIs — profile_result returns the full profile object behind confint(fit; method = :profile).

• Post-fit accessors — summary, family, is_converged, deviance, dof_residual, and rho12 (bivariate residual correlation).

• Non-Gaussian phylogenetic random effects — phylo(1 | species, tree) on the mean for Poisson, NegBinomial2, Gamma, Beta, and Binomial families (constant σ), via a sparse Laplace approximation, plus crossed intercepts (1 | g) + (1 | h) for the same families.

Where to go next

• Coming from R — choosing native R, direct Julia, or the optional engine = "julia" route.

• Rosetta (R ↔ Julia) — vocabulary and workflow translation.

• Choosing response families — the full list of response families and how to fit each one.

• What can I fit today? — the live capability map.

• When variance carries signal — a deeper location–scale walkthrough.