Getting started

DRModels.jl fits distributional regression models: the mean and the residual standard deviation (SD) of a response each get their own formula, so a covariate can change both the average and the spread. In formulas and in coef(fit, …) the mean μ is called mu and the residual SD σ is called sigma. Some response families have further parameters that can also take a formula; see Choosing response families.

What you will do

Fit a Gaussian model in which a predictor changes both an outcome's average and its residual spread. You will then read the estimated changes, check whether the optimiser converged, and see what an interval means. For the full set of currently supported models, see What can I fit today?.

This page takes you from a clean Julia session to a fitted model and shows how to turn its output into a scientific statement. The worked data are simulated so you can see the intended pattern clearly; a successful fit here is practice, not evidence about a real biological system.

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").

Some pages in this manual also load Distributions (to simulate data) or CairoMakie and AlgebraOfGraphics (to draw figures). Julia lets you load only packages you have added yourself, so add these once if a page uses them:

Pkg.add(["Distributions", "CairoMakie", "AlgebraOfGraphics"])

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

• bf(...) — bundle one formula per distributional parameter (alias drm_formula). In glmmTMB terms, bf(@formula(y ~ x), @formula(sigma ~ x)) is glmmTMB(y ~ x, dispformula = ~ x); drmTMB and brms use the same bf name.

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

Fit a first model where variability matters

This is the same example as the drmTMB package's first model, on the same data. 120 individuals were sampled in two habitats, forest and grassland. Habitat and temperature shift mean growth, and habitat also changes how much growth varies between individuals. A mean-only model would miss that second result.

The data are simulated, with these true values: mean growth 1 + 0.6·grassland + 0.4·temperature, and log residual SD -0.5 + 0.45·grassland. Julia's random numbers cannot reproduce R's set.seed(1), so the table is shipped with the package as test/data/habitat_growth.csv, generated from drmTMB's code (provenance in test/data/README.md). Reading it needs only base Julia:

using DRModels

path = pkgdir(DRModels, "test", "data", "habitat_growth.csv")
rows = split.(readlines(path)[2:end], ',')          # skip the header line
habitat     = String[r[1] for r in rows]
temperature = parse.(Float64, [r[2] for r in rows])
growth      = parse.(Float64, [r[3] for r in rows])
dat = (; growth, habitat, temperature)              # a named table of columns

fit = drm(bf(@formula(growth ~ habitat + temperature), @formula(sigma ~ habitat)),
          Gaussian(); data = dat)

Distributional regression fit (Gaussian)
  nobs = 120   logLik = -140.7119   convergence_status = converged   converged = true
  Residual SD (response scale) varies with the σ model: 0.6156 to 0.9925   dof_residual = 115

Mean model (μ):
  H0: coefficient = 0
                       Coef.  Std.Error        z   Pr(>|z|)
  (Intercept)         0.9703     0.0799  12.1420  6.322e-34
  habitat: grassland  0.5467     0.1508   3.6261  0.0002877
  temperature         0.4134     0.0781   5.2929  1.204e-07

Scale model (log σ):
  H0: coefficient = 0 on log σ
  intercept ⇔ σ = 1 (unit-dependent); slope ⇔ equal dispersion
                        Coef.  Std.Error        z   Pr(>|z|)
  (Intercept)         -0.4852     0.0913  -5.3147  1.068e-07
  habitat: grassland   0.4777     0.1291   3.6998  0.0002157

For your own data, a DataFrame (for example read from a CSV file with CSV.jl) works directly: data = df. Two pieces of Julia syntax appear above. (; growth, habitat, temperature) makes a named table from existing vectors, and @formula(…) writes a model formula, as ~ does in R.

bf(@formula(growth ~ habitat + temperature), @formula(sigma ~ habitat)) 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 residual SD — the ordinary regression that lm or lmer fits.)

Read the coefficients

coef(fit, :param) returns the coefficient block for one distributional parameter. The mean block is on the response scale (intercept, grassland contrast, temperature slope). The true values are (1.0, 0.6, 0.4); the estimates land within about one standard error of them:

coef(fit, :mu)

3-element Vector{Float64}:
 0.9703459148219472
 0.546744398029094
 0.41340703801448975

The scale block acts on log σ (intercept, grassland contrast). The true values are (-0.5, 0.45):

coef(fit, :sigma)

2-element Vector{Float64}:
 -0.4851668755620842
  0.4776545506953457

A positive grassland coefficient in the σ block means residual growth is more variable in grassland than in forest. Because the σ link is logarithmic, exponentiating a coefficient turns it into a multiplicative effect on the residual SD:

exp(coef(fit, :sigma)[2])      # grassland / forest residual-SD ratio (true value 1.57)

1.6122884232564298

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 = -140.711871958826, aic = 291.423743917652, bic = 305.36120263156226, nobs = 120, dof = 5)

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

is_converged(fit)

true

true is necessary but not sufficient. It says the optimiser stopped by its own rules; it does not prove the estimates sit at a well-defined maximum. A fit can stop on a flat stretch of the likelihood, or with a parameter at its limit (a variance at zero, a Student-t nu running to infinity), and still report true. Its stricter checks for a collapsed residual SD apply to Gaussian fits only. check_drm adds the size of the gradient (max_abs_grad) and whether the covariance is usable; read its grad_source before trusting its ok summary, and see Improving convergence.

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)             0.970346  0.0799162  12.14    <1e-33   0.813713   1.12698
mu: habitat: grassland      0.546744  0.150779    3.63    0.0003   0.251224   0.842265
mu: temperature             0.413407  0.0781055   5.29    <1e-06   0.260323   0.566491
sigma: (Intercept)         -0.485167  0.0912877  -5.31    <1e-06  -0.664087  -0.306246
sigma: habitat: grassland   0.477655  0.129101    3.70    0.0002   0.224621   0.730688
──────────────────────────────────────────────────────────────────────────────────────

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)

5-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
 (param = :mu, coef = "(Intercept)", estimate = 0.9703459148219472, lower = 0.8137131074081216, upper = 1.1269787222357728)
 (param = :mu, coef = "habitat: grassland", estimate = 0.546744398029094, lower = 0.2512238323010047, upper = 0.8422649637571833)
 (param = :mu, coef = "temperature", estimate = 0.41340703801448975, lower = 0.26032305900709707, upper = 0.5664910170218824)
 (param = :sigma, coef = "(Intercept)", estimate = -0.4851668755620842, lower = -0.6640874333558606, upper = -0.3062463177683078)
 (param = :sigma, coef = "habitat: grassland", estimate = 0.4776545506953457, lower = 0.22462105424341355, upper = 0.7306880471472779)

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)

5-element Vector{@NamedTuple{param::Symbol, coef::String, estimate::Float64, lower::Float64, upper::Float64}}:
 (param = :mu, coef = "(Intercept)", estimate = 0.9703459148219472, lower = 0.81120019816929, upper = 1.1295128210579561)
 (param = :mu, coef = "habitat: grassland", estimate = 0.546744398029094, lower = 0.24835203643241344, upper = 0.8451363828395108)
 (param = :mu, coef = "temperature", estimate = 0.41340703801448975, lower = 0.25886714240741604, upper = 0.5678713915900959)
 (param = :sigma, coef = "(Intercept)", estimate = -0.4851668755620842, lower = -0.6540233422693458, upper = -0.2949080754885855)
 (param = :sigma, coef = "habitat: grassland", estimate = 0.4776545506953457, lower = 0.22326862025589667, upper = 0.7320398583325971)

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 (the stored in-sample values), and prediction_grid for building a swept newdata grid from a reference table. predict(...; type = :response) returns invlink(Xβ) with random effects set to 0: an average group, not an average over groups.

• 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 sigma), fitted with a Laplace approximation (the method glmmTMB also uses to integrate out random effects), plus crossed intercepts (1 | g) + (1 | h) for the same families.

Random effects and the default approximation

A model with random effects, such as y ~ x + (1 | g), has a likelihood that integrates over them, and that integral has to be approximated. For an ordinary random intercept on the mean of most non-Gaussian families, DRModels.jl's default (marginal = :LA) uses adaptive Gauss–Hermite quadrature with 5 nodes per group, whereas drmTMB uses the Laplace approximation. The two packages therefore report different log-likelihoods for the same model, by about 0.001 to 0.24 units on the count, Beta, Gamma and Student datasets compared so far. To reproduce drmTMB's numbers for Poisson, Binomial or Student, or for NegBinomial2, Gamma or Beta with sigma ~ 1, pass marginal = :Laplace. That option provides a matched approximation for package comparisons. The listed log-likelihood differences do not establish which approximation is more accurate. Compare loglik, lrtest, aic and bic only between fits that use the same approximation. The drmTMB parity page lists the defaults for the main random-effect routes and the models that accept :Laplace.

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.