Using the Julia engine • drmTMB

The Julia engine is available now as an optional backend for drmTMB. You keep the R formula and data interface and select engine = "julia" when fitting a supported model. The default, engine = "tmb", does not require Julia.

Julia is useful for repeated fitting and for models that can use DRM.jl’s sparse structured solvers. It also fits ordinary regression models: a large phylogeny is not a prerequisite. Which backend is faster depends on the model, requested inference, data size, and whether Julia has already started and compiled the code.

This guide starts with an ordinary regression, then shows group-level and phylogenetic scale regression. The examples are not run during the standard R package build, so that installing or reading drmTMB does not require Julia.

Setup and prerequisites

Install Julia 1.10 or later and the optional R package JuliaCall. DRM.jl currently needs a local source checkout; installing a registered package named DRM is not the setup route used by this bridge.

Run once in a terminal, choosing a directory where you keep source packages:

git clone https://github.com/itchyshin/DRM.jl.git
julia --project=DRM.jl -e 'using Pkg; Pkg.instantiate()'

Then, in R, set the checkout path before the first Julia call:

install.packages("JuliaCall")
Sys.setenv(DRM_JL_PATH = "/absolute/path/to/DRM.jl")
library(drmTMB)

The same setup works in ordinary noninteractive scripts run with Rscript; no test opt-in is required. Julia startup remains disabled in marked R package check processes unless explicitly enabled for repository testing. Do not enable live Julia tests on CRAN or win-builder.

Replace the path with the directory containing DRM.jl’s Project.toml and src/. The bridge initializes Julia when needed. The first call includes startup and compilation; time later calls separately when comparing repeated workflows. Keep the checkout revision with your analysis so results can be reproduced; a package version alone may describe more than one source build.

Start with an ordinary regression

These data are generated in R and use the same formula for both backends:

set.seed(42)
dat <- data.frame(x = seq(-1, 1, length.out = 120))
dat$y <- 1 + 0.8 * dat$x + rnorm(nrow(dat), sd = exp(-0.4 + 0.2 * dat$x))
form <- bf(y ~ x, sigma ~ x)

fit_r <- drmTMB(form, data = dat, family = gaussian(), engine = "tmb")
fit_julia <- drmTMB(form, data = dat, family = gaussian(), engine = "julia")

summary(fit_julia)
coef(fit_julia)
predict(fit_julia, dpar = "mu")
predict(fit_julia, dpar = "sigma", type = "response")
logLik(fit_r)
logLik(fit_julia)

The mu coefficients describe the conditional mean. The sigma coefficients are on the log standard-deviation scale: exponentiating a coefficient gives a multiplicative change in residual standard deviation, not residual variance. Compare the same estimator and requested outputs when checking the two fits.

Model variation between groups

In a location-scale-scale model, predictors can enter the mean, residual standard deviation, and group-level standard deviation. Here sex is constant within each individual, and each individual has repeated observations:

set.seed(43)
individuals <- data.frame(
  individual = factor(seq_len(40)),
  sex = factor(rep(c("female", "male"), each = 20))
)
individuals$b <- rnorm(40, sd = exp(-0.3 + 0.25 * (individuals$sex == "male")))
personality_data <- individuals[rep(seq_len(40), each = 4), ]
personality_data$exploration_score <-
  1 + 0.5 * (personality_data$sex == "male") + personality_data$b +
  rnorm(nrow(personality_data), sd = 0.5)

formula_lss <- bf(
  exploration_score ~ sex + (1 | individual),
  sigma ~ sex,
  sd(individual) ~ sex
)
fit_lss <- drmTMB(
  formula_lss, family = gaussian(), data = personality_data,
  REML = TRUE, engine = "julia"
)
summary(fit_lss)
names(coef(fit_lss))

Use the returned coefficient-block names when selecting a block with coef(); formula text such as sd(individual) need not be the stored block name. Gaussian REML is available on supported routes. It accounts for estimating specified fixed effects, but does not guarantee unbiased estimates in every model or sample. Keep ML/REML choices explicit when comparing fits.

Add a phylogeny

The Julia bridge accepts rooted polytomies with positive branch lengths. Tip labels remain literal: spaces, apostrophes and Unicode do not need to be renamed. The bridge quotes labels for transport and restores observation order. Zero-length branches and unary nodes remain unsupported; this bridge uses ultrametric trees on the correlation scale.

The same scale-regression idea applies to a phylogenetic group effect. This small simulated example creates an ultrametric tree and matches data rows to its tip labels explicitly:

install.packages("ape") # once, if not already installed
set.seed(44)
bird_tree <- ape::rcoal(40)
bird_data <- data.frame(
  species = bird_tree$tip.label,
  habitat = rnorm(40),
  latitude = runif(40, -1, 1)
)
A <- ape::vcv(bird_tree, corr = TRUE)[bird_data$species, bird_data$species]
sd_phylo <- exp(-0.3 + 0.2 * bird_data$habitat)
u <- sd_phylo * drop(t(chol(A)) %*% rnorm(40))
bird_data$log_body_mass <- 1 + 0.4 * bird_data$habitat + u +
  rnorm(40, sd = exp(-0.7 + 0.1 * bird_data$latitude))

formula_phylo_lss <- bf(
  log_body_mass ~ habitat + phylo(1 | species, tree = bird_tree),
  sigma ~ latitude,
  sd(species, level = "phylogenetic") ~ habitat
)
fit_phylo_lss <- drmTMB(
  formula_phylo_lss, family = gaussian(), data = bird_data,
  REML = TRUE, engine = "julia"
)
summary(fit_phylo_lss)

Here habitat predicts both the mean and the phylogenetic standard deviation; latitude predicts residual standard deviation. The dense matrix above is only for generating this small example, not a recipe for simulating very large trees.

Missing responses

For supported Gaussian routes, response = "include" retains the full tree and group design while conditioning the likelihood on observed responses. This is not predictor imputation. Proper tree pruning can preserve covariances between observed tips; dropping a response does not inherently change those covariances.

bird_data_with_missing <- bird_data
bird_data_with_missing$log_body_mass[c(3, 9)] <- NA_real_
fit_missing <- drmTMB(
  formula_phylo_lss, family = gaussian(), data = bird_data_with_missing,
  missing = miss_control(response = "include"),
  REML = TRUE, engine = "julia"
)
summary(fit_missing)

Check each prediction method’s documented target and row behavior. Retaining unobserved tips in the model does not by itself establish support for every prediction target, such as predict(..., dpar = "sd(species)").

Inference and post-fit methods

The bridge has methods for summary(), coef(), vcov(), logLik(), predict() and confint(). Availability depends on the fitted route and target; a drmTMB_julia object does not automatically inherit every native-TMB diagnostic.

vcov(fit_julia)
confint(fit_julia, method = "wald")
confint(fit_julia, method = "profile", parm = "fixef:mu:x")

With matching development versions of drmTMB and DRM.jl, transformed terms retain their public formula names. For a model containing I(x^2), select parm = "fixef:mu:I(x^2)"; factor and interaction selectors likewise use the names returned by coef(fit_julia). Do not substitute generated Julia column numbers. Supported fixed-effect predictions on new data retain the training centering, scaling and polynomial basis. Older fitted objects without this label metadata keep their original coefficient names.

These interval examples use the ordinary ML fit above. Check the documented objective and target before comparing ML and REML inference. Agreement between engines does not establish nominal interval coverage.

Inspect conf.status and profile.message, not just the limits. A profile_failed result has an invalid endpoint; signed infinity can be a failure placeholder, and a transformed SD lower limit of zero does not clear that failure. A profile that did not cross the threshold within the searched range is a separate diagnostic. Generic Julia profiles now report failed nuisance solves, but successful optimizer termination does not by itself prove a global optimum or a practical large-tree runtime.

Which routes have bridge-side profile and bootstrap intervals

r_bridge_status in inst/extdata/julia-capabilities.tsv marks whether engine = "julia" profile and bootstrap intervals have been checked against engine = "tmb" on the same fitted target, not just whether the family fits. Two routes carry r_bridge_status "supported" today: base_gaussian_location_scale (the fixed-effect Gaussian location-scale model) and plain_binomial_nonphylo (a stats::binomial() fit without phylo()). For each, Wald and profile confidence intervals on a fixed-effect mean coefficient agree with engine = "tmb" within 1e-4, and both engines converge. Bootstrap intervals (R = 99) overlap between engines, but this is a weaker claim than profile or Wald agreement: engine = "tmb" and engine = "julia" draw bootstrap replicates from independent random-number streams, so the same seed does not reproduce the same replicates across engines, and there is no same-seed design to compare against. Treat a bootstrap-interval comparison as overlap evidence, not agreement within tolerance, even where both engines report 0/99 failed replicates.

biv_gaussian_residual (the residual-only bivariate Gaussian route, fitting rho12 without a structured or phylogenetic covariance term) is expected to join those two once drmTMB#1187 lands; as of this vignette it still carries r_bridge_status "partial" because the bridge has no profile- or bootstrap-ready target for that route on any parameter. #1187 adds one. Until it merges, use engine = "tmb" for residual bivariate Gaussian intervals.

REML standard errors are not directly comparable across engines

Where REML is admitted on the bridge (see docs/design/261-reml-by-route.md), a fixed-effect mean-block Wald standard error is not expected to match between engines by construction, not error: engine = "tmb"’s sdreport() moves the mean coefficients into TMB’s random-effect block under REML and propagates variance-parameter uncertainty into their reported standard errors, while DRM.jl profiles the mean coefficients out exactly at each variance-parameter value and reports the canonical GLS covariance (Xmu' Vhat^-1 Xmu)^-1 with no such propagation (drmTMB#1201). Two more REML routes are landing after this vignette: the mean-only phylogenetic Gaussian cell (drmTMB#1199) and the residual-only bivariate Gaussian cell (drmTMB#1197).

One modelled missing predictor (development route)

The Julia engine also has a narrow, development route for one modelled missing predictor. It is limited to a Gaussian identity-link response, one bare additive mi(x) term, complete fixed-effect exogenous designs, and a Gaussian or Bernoulli fixed-effect predictor model. Missing x values are integrated in the joint likelihood; they are not filled before fitting.

set.seed(45)
n <- 80
joint_data <- data.frame(z = seq(-1, 1, length.out = n))
joint_data$x <- 0.2 + 0.7 * joint_data$z + rnorm(n, sd = 0.25)
joint_data$y <- 0.4 + 0.5 * joint_data$z + 0.8 * joint_data$x + rnorm(n, sd = 0.3)
joint_data$x[c(8, 21, 53)] <- NA_real_
joint_data$y[c(15, 53)] <- NA_real_

fit_joint <- drmTMB(
  bf(y ~ z + mi(x), sigma ~ 1), family = gaussian(), data = joint_data,
  impute = list(x = x ~ z),
  missing = miss_control(response = "include", predictor = "model"),
  engine = "julia"
)
coef(fit_joint)
imputed(fit_joint, rows = "all")

For a binary predictor, use impute = list(x = impute_model(x ~ z, family = binomial())). The R development route also accepts response = "drop"; it drops missing-response rows during preparation. That preprocessing differs from native-TMB behaviour, so it is not a response-policy parity claim.

This route excludes other response families, additional or interacted mi() terms, random or structured effects, offsets, weights, non-default controls, and REML. summary() and Wald confint() require a usable returned covariance; profile and bootstrap intervals are unsupported. For a Gaussian predictor, the predictor-SD Wald interval is transformed to the natural SD scale, can cross zero, and is not a native-interval-parity or coverage result. Two public bridge adapter cases pass. Native training prediction and binary newdata handling have been repaired and checked independently. Full numerical parity remains open: small differences in optimizer stopping affect coefficients and predictions beyond the declared tolerance. This development route makes no speed or interval-coverage claim.

Species counts, sparse solvers and parallel work

There is no universal 5,000-species limit. The limit belongs to particular dense Gaussian location-scale-scale (LSS) routes, and is a limit on observations, not a general Julia limit.

For a single phylogenetic LSS component, DRM.jl automatically chooses its sparse tree engine once there are more than 500 species. In direct Julia use, you can also request that route explicitly with sparse = true or algorithm = :sparse_lbfgs. The sparse route avoids constructing the dense species-by-species covariance matrix. A deliberately forced dense fallback, and the current multi-component LSS route, still stop at 5,000 observations. Repeated observations per species can therefore reach the dense limit before a dataset has 5,000 species.

Other sparse phylogenetic routes have different practical limits. Species count, observation count, number of latent effects, and requested uncertainty all affect cost.

A dense double-precision 5,000 by 5,000 matrix alone occupies about 200 MB, before factorizations and temporary storage. Its storage grows quadratically and a general dense factorization grows cubically. Tree-based augmented-state solvers can avoid the dense covariance, but the precision is tree-sparse, not generally tridiagonal. Crossed effects and other structures can introduce additional fill.

Sparse likelihood calculations, model construction, standard errors, profiles, and bootstrap refits are different workloads. A fast likelihood evaluation is not evidence that a complete fit or interval calculation takes the same time. There is no single speedup factor for all Julia-backed models.

To make multiple Julia threads available, set JULIA_NUM_THREADS before Julia starts (for example in the shell launching R). A multi-threaded runtime does not make every fit parallel automatically. Use the documented options for the specific workflow, and avoid combining many Julia workers with many BLAS threads.

Current boundaries

The Julia backend is implemented, but it is not yet a replacement for every native-TMB workflow. Check the capability guide for the model you intend to fit. In particular:

The one-modelled-predictor development route above is available through engine = "julia"; other missing-predictor models remain native-TMB workflows.
Missing-response inclusion and REML depend on the model route. Most non-Gaussian bridge routes fit ML only, but large-p Poisson phylo() is a genuine exception: the bridge fits REML = TRUE there by Cox-Reid Laplace, a route native engine = "tmb" does not have. See the REML column in the table below and docs/design/261-reml-by-route.md before assuming ML-only.
Accepted families do not imply that every combination of structured effects, weights, controls and post-fit methods is admitted. phylo(), relmat(), animal() and spatial() terms on the fixed-effect-only families below are already refused by existing gates before Julia starts. An ordinary random-effect term ((1 | g), a random slope) or an sd()/sd_phylo() scale submodel on one of these families lands with drmTMB#1196: for student(), lognormal(), truncated_nbinom2(), zero_one_beta(), tweedie(), beta_binomial(), cumulative_logit(), and skew_normal() once #1176 lands, such a term is refused before Julia starts, naming the family and pointing to engine = "tmb". Until #1196 merges, that combination is not refused by drmTMB before Julia starts; the call round-trips to DRM.jl and DRM.jl refuses it in its own words instead (for example, TruncatedNegBinomial2() currently supports fixed effects only, pinned in tests/testthat/test-julia-family-truncated_nbinom2.R). Treat any such combination as unsupported regardless of which side reports the refusal.
predict() on new data currently has narrower Julia support than in-sample prediction; do not assume that every distributional parameter is available.

If the bridge rejects a combination, use the native engine for that analysis and retain the warning or error when reporting the limitation.

Family and route reference

Every family below has a same-target parity receipt against engine = "tmb" on a committed fixture (coefficients and log-likelihood within 1e-4, Wald SEs within 1e-3 relative), except the one row marked “not yet admitted”. The “Receipt” column names the capability-ledger row (inst/extdata/julia-capabilities.tsv, r_bridge_status column) where one exists; three fixed-effect-only families were admitted without a ledger row yet, so their receipt lives only in the family’s own test file.

Family	`engine = "julia"` route	Receipt
`gaussian()`	fixed effect; mean-only `phylo()`; coupled mean+`sigma` phylo; ordinary random intercept/slope on `mu` and `sigma`; `relmat()`/`spatial()`/`animal()`	`base_gaussian_location_scale`, `gaussian_phylo_mean`, `gaussian_random_intercept_mu`, `gaussian_random_slope_mu`, `gaussian_sigma_random_intercept`, `general_covariance_structured`, `location_scale_scale` (partial/experimental; see the TSV for which)
`biv_gaussian()`	fixed effect (residual `rho12`); coupled q = 4 phylo	`biv_gaussian_residual`, `biv_q4_phylo_reml`
`student()`	fixed effect only	`fe_student`
`lognormal()`	fixed effect only	`fe_lognormal`
`poisson()`	fixed effect; large-p `phylo()`; `zi` dpar on the fixed-effect route	`fe_poisson`, `phylo_count_large_p`, `zi_poisson`
`nbinom2()`	fixed effect; large-p `phylo()`; coupled mean+`sigma` phylo; `zi` dpar; `hu` dpar (bridge spelling only – see note below)	`fe_nbinom2`, `phylo_count_large_p`, `zi_nbinom2`, `hurdle_nbinom2`
`gamma()`	fixed effect; large-p `phylo()`; coupled mean+`sigma` phylo	`fe_gamma`, `phylo_gamma_beta_binomial`
`beta()`	fixed effect; large-p `phylo()`; coupled mean+`sigma` phylo	`fe_beta`, `phylo_gamma_beta_binomial`
`stats::binomial(link = "logit")`	fixed effect; large-p mean-only `phylo()`	`plain_binomial_nonphylo`, `phylo_gamma_beta_binomial`
`truncated_nbinom2()`	fixed effect only	no ledger row yet; `tests/testthat/test-julia-family-truncated_nbinom2.R`
`zero_one_beta()`	fixed effect only	no ledger row yet; `tests/testthat/test-julia-family-zero_one_beta.R`
`tweedie()`	fixed effect only	no ledger row yet; `tests/testthat/test-julia-family-tweedie.R`
`beta_binomial()`	fixed effect only (`phylo()` refuses: DRM.jl’s phylo `BetaBinomial` route is constant-`sigma` only, no receipt)	`fe_beta_binomial`
`cumulative_logit()`	fixed effect only	no ledger row yet; `tests/testthat/test-julia-family-cumulative_logit.R`
`skew_normal()`	not yet admitted (open PR)	none
any other family (`zi_poisson()`/`zi_nbinom2()`/`hurdle_nbinom2()` as their own family constructors, mixture families, etc.)	refused	none

The zi and hu dpars are not separate families: they are formula terms (bf(y ~ x, zi ~ x) or bf(y ~ x, hu ~ x)) added to poisson() or nbinom2(), and the bridge admits them through those families’ fixed-effect route rather than through a family-specific ledger row. hu has a spelling mismatch between engines today: native engine = "tmb" fits the hurdle model as truncated_nbinom2() with hu ~ ..., while engine = "julia" fits it as nbinom2() with hu ~ ... (DRM.jl reads hu on nbinom2() as the hurdle model). The same call does not fit on both engines yet; see hurdle_nbinom2 in the capability table for the measured cross-spelling receipt and the tracked defect. ### Julia optimizer controls

For the base Julia bridge, use the Julia-native controls through the existing drm_control() interface:

control = drm_control(optimizer = list(g_tol = 1e-6, algorithm = "lbfgs"))

These settings control DRM.jl’s gradient tolerance and solver choice. They do not make its optimizer equivalent to TMB’s nlminb() presets or iteration budgets.

The forwarded set is a closed whitelist, and it is short:

`drm_control()` setting	reaches DRM.jl as	where
`optimizer$g_tol`	`drm(; g_tol =)`	every base / phylo / bivariate route
`optimizer$algorithm`	`drm(; algorithm =)`, one of `"auto"`, `"gls"`, `"lbfgs"`, `"em"`, `"sparse"`, `"sparse_lbfgs"`	as above, except q = 4
`optimizer$q4_vcov`	`drm(; q4_vcov =)`	bivariate q = 4 phylogenetic only

On the bivariate q = 4 phylogenetic route optimizer$g_tol is forwarded as DRM.jl’s q4_g_tol (that route’s outer-gradient tolerance) and optimizer$algorithm is refused, because its optimiser has no solver-selection setting.

Everything else drm_control() carries is rejected before Julia starts, with an error naming the setting – it is never silently ignored. That includes the storage flags (se, keep_data, keep_model_frame, keep_tmb_object), the SE-shaping flags (se_report_covariance, se_skip_delta_method, se_group_sd), sparse_fixed, aggregate_gaussian, logsigma_clamp, logsigma_clamp_margin, optimizer_preset, newton_polish, multi_start, fallback_optimizer, start, and any nlminb() name inside optimizer (iter.max, eval.max, rel.tol, …). Each of these describes an nlminb()/TMB program that DRM.jl does not run, so there is nothing to forward it to: use engine = "tmb" when you need one of them.

The structured, bivariate q2 structured, and cross-family Julia routes are narrower again: they accept only a default control, so even optimizer$g_tol is refused there.