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 the R ↔ Julia bridge and Rosetta; for the full capability map see What can I fit today?.

DRM.jl is distributional regression: instead of a single linear predictor for the mean, you give each parameter of the response distribution its own formula. 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

DRM.jl is pre-release, so develop it from a local checkout:

using Pkg
Pkg.develop(path = "/path/to/DRM.jl")   # or Pkg.add(url = "https://github.com/itchyshin/DRM.jl")
using DRM

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

• bf(...) — bundle one formula per distributional parameter (alias drm_formula), exactly like drmTMB / brms.

• 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 DRM, 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

Mean model (μ):
                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 σ):
                 Coef.  Std.Error    z  Pr(>|z|)
  (Intercept)  -0.4008     0.0354  NaN       NaN
  x             0.3181     0.0363  NaN       NaN

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, and 95% Wald intervals for every block at once. In this first-fit scale block, the z and p entries are currently unavailable (NaN); that is not evidence for a zero scale effect, and the estimate, standard error, and interval remain finite outputs. Supplying scale-block z/p results consistent with the R-facing post-fit contract remains an open parity obligation. 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  NaN       NaN     -0.470112  -0.331395
sigma: x             0.318075  0.0363008  NaN       NaN      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

• R ↔ Julia bridge — the experimental engine = "julia" route and the cells it admits today.

• 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.