Comparing drmTMB with other packages

library(drmTMB)

If you already fit mixed models in lme4, meta-analyses in metafor, proportional-odds models in ordinal, or location-scale models in glmmTMB, this article answers one question in each of those packages’ own vocabulary: on a dataset you may already know, does drmTMB land on the same fit?

Eight comparisons follow, each built from a real dataset, a real drmTMB fit, and a real fit from the comparator package alongside it. These eight models were chosen because a package this article already depends on fits them too; that set is small next to what drmTMB implements as a whole, and the package’s design documents record the rest of that surface.

Every comparison here is single-seed and single-dataset. Each one shows that drmTMB’s likelihood and optimizer reach the same optimum as another implementation of the same model, on the data shown. It is not evidence about interval calibration, coverage, bias, or small-sample behaviour, and no standard error or confidence interval is compared anywhere in this article. Agreement licenses the model both packages fit on the dataset shown; for any other model, drmTMB’s evidence is its own recovery and coverage studies, recorded per capability in the package’s design documents and capability ledger.

What this article compares

drmTMB’s location submodel is compared in every one of the eight fits below. Its distributional submodel — the sigma formula that lets residual spread depend on covariates — is compared in some of them, and its random-effect form is not compared here. The table states only what this article did, comparison by comparison.

drmTMB capability Comparisons Compared in this article? Independence
Location (mu) submodel, fixed effects 1–8 (all) yes, eight times 1–3 STRONG · 4–5 unclassified · 6–8 WEAK
Location (mu) submodel, ordinary random intercept / slope 1, 5, 6, 8 yes, four times 1 STRONG · 5 unclassified · 6, 8 WEAK
Distributional (sigma) submodel, fixed effects 2, 3, 6, 7, 8 yes, five times 2, 3 STRONG (a between-study heterogeneity SD, not a residual scale; 2’s sigma formula is intercept-only, 3’s carries a covariate) · 6, 7, 8 WEAK
Distributional (sigma) submodel, ordinary random effect none not compared here nothing to classify
Residual correlation rho12 with a predictor none not compared here nothing to classify
The rest of the implemented model surface none not compared here; the design documents record what it contains nothing to classify

Comparisons 6, 7, and 8 fit mu and sigma jointly in one model. The agreement reported for each is a property of that joint fit; it does not license either submodel taken on its own.

Three independence labels appear throughout: STRONG means the comparator shares no estimation code with drmTMB; WEAK means it is built on the same TMB automatic-differentiation stack and outer optimizer, so agreement is a consistency check between related implementations rather than a cross-implementation one; unclassified means the comparator is a separate engine whose independence drmTMB’s design documents have not yet classified.

Checked against a separate estimation engine

The three comparisons in this section are against lme4 and metafor, which share no estimation code with drmTMB (docs/design/242-external-comparator-evidence-class.md). Agreement here is a genuine cross-implementation check.

Comparison 1: contagious bovine pleuropneumonia across periods and herds

If you already fit lme4::glmer() on lme4::cbpp while learning mixed models, this is the same model in drmTMB: infection risk by observation period, with a random herd intercept.

data(cbpp, package = "lme4")
fit1 <- drmTMB(
  bf(mu = cbind(incidence, size - incidence) ~ period + (1 | herd)),
  data = cbpp, family = binomial()
)
coef(fit1, "mu")
#> (Intercept)     period2     period3     period4 
#>  -1.3985281  -0.9923386  -1.1286497  -1.5803225
cmp1 <- lme4::glmer(
  cbind(incidence, size - incidence) ~ period + (1 | herd),
  family = binomial, data = cbpp
)
lme4::fixef(cmp1)
#> (Intercept)     period2     period3     period4 
#>   -1.398343   -0.991925   -1.128216   -1.579745

Both sides fit the same logit-scale coefficients directly — no conversion — and the herd random-intercept SD compares to lme4’s the same way:

summary(fit1)$parameters[, c("term", "estimate")]
#>                        term  estimate
#> sd:mu:(1 | herd) (1 | herd) 0.6422603
attr(lme4::VarCorr(cmp1)$herd, "stddev")
#> (Intercept) 
#>   0.6420699
coef(fit1, "mu") - lme4::fixef(cmp1)
#>   (Intercept)       period2       period3       period4 
#> -0.0001852107 -0.0004136493 -0.0004334457 -0.0005771155
c(drmTMB = as.numeric(logLik(fit1)), glmer = as.numeric(logLik(cmp1)))
#>    drmTMB     glmer 
#> -92.02628 -92.02657

The two fits agree to about three decimal places on the coefficients and the log-likelihood, not to five or six. Both sides use a first-order (nAGQ = 1) Laplace approximation of the marginal likelihood; a gap of this size between two different Laplace solvers — TMB’s automatic-differentiation inner solve here, lme4’s PIRLS solve there — is the expected size for that comparison, not a sign that either optimizer stopped early.

binomial() here has no dispersion parameter, on either side. cbpp’s known extra-binomial variation is therefore unmodelled symmetrically, not a discrepancy between the packages; beta_binomial() adds a dispersion parameter to this same two-column response if that matters for your own data. (nbinom2() also carries a dispersion parameter, but it is a count family and takes a single-column count response, so it is not a substitute here.)

Independence: STRONG (lme4).

Comparison 2: the pooled BCG vaccine effect

If you run meta-analyses in metafor, this is a fixed pooled effect with a between-trial heterogeneity term, written as drmTMB gaussian regression with the sampling variance supplied through meta_V().

data(dat.bcg, package = "metadat")
dat_bcg <- metafor::escalc(
  measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg
)
fit2 <- drmTMB(
  bf(mu = yi ~ 1 + meta_V(V = vi), sigma = ~ 1),
  data = dat_bcg, family = gaussian()
)
coef(fit2, "mu")
#> (Intercept) 
#>   -0.711199

metafor::rma.uni()’s default estimator is REML, and it gives a visibly different fit:

cmp2_reml <- metafor::rma.uni(yi, vi, data = dat_bcg)
c(coef(cmp2_reml), tau2 = cmp2_reml$tau2)
#>    intrcpt       tau2 
#> -0.7145323  0.3132433

method = "ML" is the comparator that matches:

cmp2 <- metafor::rma.uni(yi, vi, data = dat_bcg, method = "ML")
c(coef(cmp2), tau2 = cmp2$tau2)
#>    intrcpt       tau2 
#> -0.7111991  0.2800282

drmTMB reports log(tau) as the sigma intercept; metafor reports tau^2 directly. Converting with tau^2 = exp(2 * coef(fit, "sigma")) and comparing the two log-likelihoods shows the two packages are maximising the same function, not merely landing close to each other:

c(drmTMB_tau2 = exp(2 * coef(fit2, "sigma")), metafor_tau2 = cmp2$tau2)
#> drmTMB_tau2.(Intercept)            metafor_tau2 
#>               0.2800281               0.2800282
c(drmTMB = as.numeric(logLik(fit2)), metafor = as.numeric(logLik(cmp2)))
#>    drmTMB   metafor 
#> -12.66508 -12.66508

BCG vaccination reduces tuberculosis risk on average, with real heterogeneity between the 13 trials, and drmTMB reproduces the estimate a metafor user already trusts before being asked to trust anything new.

Independence: STRONG (metafor, docs/design/242-external-comparator-evidence-class.md).

Comparison 3: does heterogeneity itself differ by allocation design?

The 13 BCG trials were allocated to treatment in three different ways (alternate, random, systematic). Rather than splitting the data into three separate meta-analyses, one drmTMB fit gives each allocation type its own heterogeneity tau through sigma ~ alloc.

fit3 <- drmTMB(
  bf(mu = yi ~ 1 + meta_V(V = vi), sigma = ~ alloc),
  data = dat_bcg, family = gaussian()
)
coef(fit3, "sigma")
#>     (Intercept)     allocrandom allocsystematic 
#>       -1.676592        1.175489        1.047993
dat_bcg$id <- seq_len(nrow(dat_bcg))
cmp3 <- metafor::rma.mv(
  yi, vi, random = ~ alloc | id, struct = "DIAG", data = dat_bcg, method = "ML"
)
cmp3$tau2
#> [1] 0.03497276 0.36706788 0.28444971

drmTMB’s sigma linear predictor is on the log(tau) scale, and each allocation level’s value is the intercept plus that level’s contrast, with alternate as the reference. rma.mv() under struct = "DIAG" has no reference level: it estimates one free tau^2 per level directly, so its four parameters are one beta plus three tau^2. What makes the two line up is therefore level order, not a shared reference — rma.mv() reports in the order of g.levels.f, and the comparison below relies on that. Exponentiating drmTMB’s side gives per-level tau, matching rma.mv()’s response-scale tau:

sigma3 <- coef(fit3, "sigma")
tau_by_level <- exp(c(
  alternate = sigma3[1],
  random = sigma3[1] + sigma3[2],
  systematic = sigma3[1] + sigma3[3]
))
tau_by_level
#>  alternate.(Intercept)     random.(Intercept) systematic.(Intercept) 
#>              0.1870101              0.6058616              0.5333383
sqrt(cmp3$tau2)
#> [1] 0.1870100 0.6058613 0.5333383
c(drmTMB = as.numeric(logLik(fit3)), metafor = as.numeric(logLik(cmp3)))
#>    drmTMB   metafor 
#> -11.99732 -11.99732

The two fits’ per-level tau values and log-likelihood agree closely, and that agreement is not a coincidence: with one row per study, rma.mv()‘s ~ alloc | id under struct = "DIAG" puts one random effect per study with variance tau^2 for that study’s allocation level, so the marginal variance is v_i + tau^2_{alloc(i)} on both sides. These are the same model, written in two packages’ syntax, not two models that happen to agree.

Read the three fitted tau values as a point-agreement result only, not as a finding about which allocation design is more heterogeneous: the alternate level carries two studies, random seven, and systematic four, and drmTMB’s own recorded evidence for this route says heterogeneity intervals are not usable at this few studies per level. Report the estimates this comparison licenses — that the two packages reach the same optimum — and no more.

Independence: STRONG. This and Comparison 2 are the two comparisons in this article that check a sigma linear predictor against a separate estimation engine — Comparison 2’s sigma formula is intercept-only, and this one gives it a covariate. In both, sigma is a between-study heterogeneity SD under meta_V(), not a residual scale.

Checked against ordinal

ordinal is a separate estimation engine, but docs/design/242-external-comparator-evidence-class.md has not yet classified its independence strength the way it classifies lme4, metafor, and glmmTMB, so the two comparisons in this section are reported without one.

Comparison 4: does serving temperature shift wine bitterness ratings?

If you fit proportional-odds models with ordinal::clm(), this is the same grammar in cumulative_logit(): five ordered bitterness ratings, from Randall (1989)’s wine-tasting data, as a function of serving temperature and skin contact.

data(wine, package = "ordinal")
fit4 <- drmTMB(
  bf(mu = rating ~ temp + contact), data = wine, family = cumulative_logit()
)
summary(fit4)$coefficients
#>               estimate std_error
#> mu:tempwarm   2.503102 0.5286802
#> mu:contactyes 1.527797 0.4766226
summary(fit4)$ordinal$cutpoints
#>       1|2       2|3       3|4       4|5 
#> -1.344385  1.250807  3.466888  5.006405
cmp4 <- ordinal::clm(rating ~ temp + contact, data = wine)
coef(cmp4)
#>        1|2        2|3        3|4        4|5   tempwarm contactyes 
#>  -1.344383   1.250809   3.466887   5.006404   2.503102   1.527798
c(drmTMB = as.numeric(logLik(fit4)), clm = as.numeric(logLik(cmp4)))
#>    drmTMB       clm 
#> -86.49192 -86.49192

No conversion is needed: both sides write the same cumulative-logit convention, logit P(Y <= j) = alpha_j - x'beta, and the two log-likelihoods agree to five decimal places with matching-sign coefficients. Warm serving and skin contact both push ratings toward “more bitter.”

The comparison that makes this pair valuable runs the other way. ordinal can add a scale formula to this same model, and it fits:

cmp4_scale <- ordinal::clm(rating ~ temp + contact, scale = ~ temp, data = wine)
as.numeric(logLik(cmp4_scale))
#> [1] -86.43946

cumulative_logit() does not accept the equivalent formula in drmTMB 0.7.0:

drmTMB(
  bf(mu = rating ~ temp + contact, sigma = ~ temp),
  data = wine, family = cumulative_logit()
)
#> Error in `drm_build_cumulative_logit_spec()`:
#> ! `cumulative_logit()` models currently support only a `mu` location
#>   formula.
#> ✖ Unsupported parameter: "sigma".
#> ℹ Ordinal scale/discrimination formulas are planned after the identifiability
#>   contract is finalized.

cumulative_logit() currently fits only a mu location formula; a scale or discrimination formula for this family is not yet implemented. If you need that, ordinal::clm(scale = ~ ...) is the tool that already has it. This asymmetry is worth sitting with: the two packages’ surfaces do not coincide, and they fail to coincide in the direction that costs drmTMB something, which is the direction worth believing.

Independence: ordinal is not yet classified in drmTMB’s design documents.

Comparison 5: do judges differ, and does that change the temperature effect?

Nine judges each rated eight bottles. Adding a judge random intercept asks whether that matters, and whether drmTMB still matches ordinal::clmm() once it is added.

fit5 <- drmTMB(
  bf(mu = rating ~ temp + contact + (1 | judge)),
  data = wine, family = cumulative_logit()
)
summary(fit5)$coefficients
#>               estimate std_error
#> mu:tempwarm   3.063001 0.5953910
#> mu:contactyes 1.834900 0.5125582
summary(fit5)$ordinal$cutpoints
#>       1|2       2|3       3|4       4|5 
#> -1.623669  1.513357  4.228520  6.088769
cmp5 <- ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine)
coef(cmp5)
#>        1|2        2|3        3|4        4|5   tempwarm contactyes 
#>  -1.623667   1.513365   4.228527   6.088773   3.062997   1.834885

For the random effect, compare SD to SD, not clmm’s printed variance:

summary(fit5)$parameters[, c("term", "estimate")]
#>                          term estimate
#> sd:mu:(1 | judge) (1 | judge) 1.131139
attr(ordinal::VarCorr(cmp5)$judge, "stddev")
#> (Intercept) 
#>    1.131133
c(drmTMB = as.numeric(logLik(fit5)), clmm = as.numeric(logLik(cmp5)))
#>    drmTMB      clmm 
#> -81.56541 -81.56541

Every quantity here — slopes, cutpoints, the random-intercept SD, and the log-likelihood — agrees closely between the two fits. Ignoring judge identity attenuates the temperature effect: Comparison 4’s marginal slope was smaller than the conditional slope shown here, which is the expected conditional-versus-marginal shift, and having a working comparator is what makes it checkable rather than merely plausible.

This licenses point-fit parity only, on this dataset. It does not license any REML or interval claim for cumulative_logit() random effects in general. The capability-ledger cell behind the model fitted here is mc-0225 (cumulative_logit, mu, ordinary random intercept, ML), recorded at interval_feasible. That tier is narrow: it establishes computational interval feasibility for that one cell’s direct intercept-SD target on a frozen low-rung fixture, it was not evaluated for coverage or calibration, and it carries no public interval guidance, no broader family or formula support, and nothing about correlated or labelled effects. A random slope on this family is a different cell (mc-0227, recorded one tier lower at point_fit_recovery, the rung directly beneath interval_feasible) and is not what this comparison fits.

Independence: ordinal is not yet classified in drmTMB’s design documents.

Checked against a package built on the same machinery

The three comparisons in this section are against glmmTMB, which is built on the same TMB automatic-differentiation stack and outer optimizer as drmTMB (docs/design/242-external-comparator-evidence-class.md). Agreement here is a consistency check between related implementations, not a cross-implementation check. It is also, uncomfortably, where drmTMB’s stated differentiator lives: no separate estimation engine checked a residual-scale submodel anywhere in this article.

glmmTMB is called as glmmTMB::glmmTMB() below and is never attached with library(). That is deliberate. glmmTMB shares several names with drmTMB – the family constructors lognormal() and nbinom2(), and the extractor ranef() – so attaching it silently changes which package’s function a bare call reaches. The drmTMB:: prefixes below are for the same reason, not style.

Comparison 6: does reaction-time variability grow with sleep deprivation?

If you already fit lme4::lmer(Reaction ~ Days + (Days | Subject)) on lme4::sleepstudy, this asks a further question: does the spread of reaction times widen with each day of deprivation, on top of the mean?

The scale conversion is checked first on a fixed-effect variant of the model, where the likelihood is exact and no Laplace step intervenes on either side:

data(sleepstudy, package = "lme4")
fit6_fe <- drmTMB(
  bf(mu = Reaction ~ Days, sigma = ~ Days), data = sleepstudy, family = gaussian()
)
mu_hat <- predict_parameters(fit6_fe, dpar = "mu", type = "response")$estimate
eta_sigma <- predict_parameters(fit6_fe, dpar = "sigma", type = "link")$estimate
c(
  reported = as.numeric(logLik(fit6_fe)),
  `hand, sd = exp(eta)` = sum(dnorm(sleepstudy$Reaction, mu_hat, exp(eta_sigma), log = TRUE)),
  `hand, sd = sqrt(exp(eta))` = sum(dnorm(sleepstudy$Reaction, mu_hat, sqrt(exp(eta_sigma)), log = TRUE))
)
#>                  reported       hand, sd = exp(eta) hand, sd = sqrt(exp(eta)) 
#>                 -938.7164                 -938.7164                -4647.8169

drmTMB’s sigma linear predictor is log(SD); no transform, no squaring. The wrong alternative misses drmTMB’s own reported log-likelihood by thousands of units, so the link is structural, not something that happens to hold for one fit.

On the displayed random-effect model, that same no-transform rule is the one to use, and the two packages’ sigma coefficients compare directly to glmmTMB’s dispformula coefficients:

fit6 <- drmTMB(
  bf(mu = Reaction ~ Days + (1 + Days | Subject), sigma = ~ Days),
  data = sleepstudy, family = gaussian()
)
coef(fit6, "mu")
#> (Intercept)        Days 
#>   252.85599    10.08607
coef(fit6, "sigma")
#> (Intercept)        Days 
#>  2.81184909  0.08463394
cmp6 <- glmmTMB::glmmTMB(
  Reaction ~ Days + (Days | Subject), dispformula = ~ Days, data = sleepstudy
)
glmmTMB::fixef(cmp6)$cond
#> (Intercept)        Days 
#>   252.85595    10.08606
glmmTMB::fixef(cmp6)$disp
#> (Intercept)        Days 
#>   2.8118492   0.0846339

With a random effect on mu, logLik() is a Laplace marginal on both sides, so no row-by-row density identity holds the way it did above. What the two packages share instead is the marginal log-likelihood itself:

c(drmTMB = as.numeric(logLik(fit6)), glmmTMB = as.numeric(logLik(cmp6)))
#>    drmTMB   glmmTMB 
#> -870.0003 -870.0003
abs(coef(fit6, "sigma") - glmmTMB::fixef(cmp6)$disp)
#>  (Intercept)         Days 
#> 1.569626e-07 3.767808e-08

The sigma slope on Days is positive: residual spread grows with each day of deprivation. A mean-only model cannot show that, because it carries a single residual SD for every observation; this is an addition to the mean model you already trust, not a replacement for it.

Independence: WEAK. The two fits are a consistency check between related implementations, not an independent confirmation.

Comparison 7: are male penguins more variable in body mass than females?

palmerpenguins::penguins lets one model ask two things at once: does mean body mass differ by species, and is one sex more variable than the other?

pen <- palmerpenguins::penguins[
  stats::complete.cases(
    palmerpenguins::penguins[, c("species", "sex", "body_mass_g")]
  ),
]
fit7 <- drmTMB(
  bf(mu = body_mass_g ~ species, sigma = ~ sex),
  data = pen, family = drmTMB::lognormal()
)
coef(fit7, "mu")
#>      (Intercept) speciesChinstrap    speciesGentoo 
#>       8.17272037       0.02447963       0.32491983
coef(fit7, "sigma")
#> (Intercept)     sexmale 
#>  -2.4243950   0.4374102
cmp7 <- glmmTMB::glmmTMB(
  log(body_mass_g) ~ species, dispformula = ~ sex, data = pen, family = gaussian()
)
glmmTMB::fixef(cmp7)$cond
#>      (Intercept) speciesChinstrap    speciesGentoo 
#>       8.17272036       0.02447967       0.32491980
glmmTMB::fixef(cmp7)$disp
#> (Intercept)     sexmale 
#>  -2.4243951   0.4374103

Compare on the log scale, coefficient for coefficient, with sigma unsquared. Two things matter here that a coefficient-by-coefficient eyeball can miss. First, drmTMB’s mu is E[log y], not log E[y]:

mu_hat <- predict_parameters(fit7, dpar = "mu", type = "link")$estimate
eta_sigma <- predict_parameters(fit7, dpar = "sigma", type = "link")$estimate
y <- pen$body_mass_g
c(
  reported = as.numeric(logLik(fit7)),
  `hand, meanlog = mu` = sum(dlnorm(y, mu_hat, exp(eta_sigma), log = TRUE)),
  `hand, meanlog = mu - sigma^2/2` = sum(dlnorm(y, mu_hat - exp(eta_sigma)^2 / 2, exp(eta_sigma), log = TRUE))
)
#>                       reported             hand, meanlog = mu 
#>                      -2511.448                      -2511.448 
#> hand, meanlog = mu - sigma^2/2 
#>                      -2517.472

Second, the two packages’ raw log-likelihoods are not comparable directly: drmTMB reports it on the original body_mass_g scale, including the log-Jacobian for the log transform, while glmmTMB reports it for log(body_mass_g). The gap between them is exactly that Jacobian term:

c(drmTMB = as.numeric(logLik(fit7)), glmmTMB = as.numeric(logLik(cmp7)))
#>     drmTMB    glmmTMB 
#> -2511.4482   261.3321
sum(-log(y))
#> [1] -2772.78

Coefficient agreement, once the two traps above are accounted for, is close on every term. Gentoo penguins are heavier than Adelie on the log scale, and males are noticeably more variable in log body mass than females — the sigma ~ sex submodel is what makes that second finding representable at all; a model with one residual SD has no parameter for it.

Independence: WEAK.

Comparison 8: are satiated owl broods more variable, not just quieter?

glmmTMB::Owls records sibling negotiation calls by nest, food treatment, and parent sex. Beyond the mean call rate, does food treatment change the overdispersion too?

The scale conversion is again checked first on a fixed-effect variant, where the likelihood is exact:

data(Owls, package = "glmmTMB")
fit8_fe <- drmTMB(
  bf(mu = SiblingNegotiation ~ FoodTreatment, sigma = ~ FoodTreatment),
  data = Owls, family = drmTMB::nbinom2()
)
mu_hat <- predict_parameters(fit8_fe, dpar = "mu", type = "response")$estimate
eta_sigma <- predict_parameters(fit8_fe, dpar = "sigma", type = "link")$estimate
y <- Owls$SiblingNegotiation
c(
  reported = as.numeric(logLik(fit8_fe)),
  `hand, size = 1/sigma^2` = sum(dnbinom(y, mu = mu_hat, size = 1 / exp(eta_sigma)^2, log = TRUE)),
  `hand, size = 1/sigma` = sum(dnbinom(y, mu = mu_hat, size = 1 / exp(eta_sigma), log = TRUE)),
  `hand, size = sigma` = sum(dnbinom(y, mu = mu_hat, size = exp(eta_sigma), log = TRUE))
)
#>               reported hand, size = 1/sigma^2   hand, size = 1/sigma 
#>              -1725.172              -1725.172              -1735.690 
#>     hand, size = sigma 
#>              -1822.905

drmTMB’s sigma maps to size = 1 / sigma^2; glmmTMB’s dispformula predictor is log(size) directly. The conversion is therefore log(theta) = -2 * log(sigma), and only the first hand-built likelihood above reproduces drmTMB’s reported value — the other two miss it by a wide margin, in an order of magnitude that could otherwise look plausible if you did not check it.

The drmTMB:: prefix on nbinom2() above is mandatory here, not stylistic. glmmTMB also exports a family constructor named nbinom2(). If you attach glmmTMB with library(), a bare family = nbinom2() silently hands drmTMB glmmTMB’s family object instead of its own, and the call is rejected. This article never attaches glmmTMB, so the clash is shown explicitly:

drmTMB(
  bf(mu = SiblingNegotiation ~ FoodTreatment, sigma = ~ FoodTreatment),
  data = Owls, family = glmmTMB::nbinom2()
)
#> Error in `drm_family_type()`:
#> ! Currently supported families are `gaussian()`, `student()`,
#>   `skew_normal()`, `lognormal()`, `biv_lognormal()`, `biv_student()`,
#>   `Gamma(link = "log")`, `tweedie()`, `beta()`, `zero_one_beta()`,
#>   `beta_binomial()`, `binomial(link = "logit"/"probit"/"cloglog")`,
#>   `cumulative_logit()`, `poisson(link = "log")`, `nbinom2()`,
#>   `truncated_nbinom2()`, `biv_gaussian()`, `c(gaussian(), gaussian())`, and
#>   `list(gaussian(), gaussian())`. Zero-inflated Poisson and NB2 models use the
#>   same family route plus a `zi ~ ...` formula; hurdle NB2 models use
#>   `truncated_nbinom2()` plus a `hu ~ ...` formula.

On the displayed random-effect model, the joint fits agree on the marginal log-likelihood, and the same -2 * conversion carries over to the sigma coefficients:

fit8 <- drmTMB(
  bf(
    mu = SiblingNegotiation ~ FoodTreatment * SexParent + (1 | Nest),
    sigma = ~ FoodTreatment
  ),
  data = Owls, family = drmTMB::nbinom2()
)
coef(fit8, "mu")
#>                         (Intercept)               FoodTreatmentSatiated 
#>                          2.03752206                         -0.63582963 
#>                       SexParentMale FoodTreatmentSatiated:SexParentMale 
#>                          0.03176616                          0.10467627
coef(fit8, "sigma")
#>           (Intercept) FoodTreatmentSatiated 
#>            -0.1887961             0.6040870
cmp8 <- glmmTMB::glmmTMB(
  SiblingNegotiation ~ FoodTreatment * SexParent + (1 | Nest),
  dispformula = ~ FoodTreatment, family = glmmTMB::nbinom2, data = Owls
)
glmmTMB::fixef(cmp8)$cond
#>                         (Intercept)               FoodTreatmentSatiated 
#>                          2.03752367                         -0.63582849 
#>                       SexParentMale FoodTreatmentSatiated:SexParentMale 
#>                          0.03176311                          0.10467526
glmmTMB::fixef(cmp8)$disp
#>           (Intercept) FoodTreatmentSatiated 
#>             0.3775933            -1.2081753
c(drmTMB = as.numeric(logLik(fit8)), glmmTMB = as.numeric(logLik(cmp8)))
#>    drmTMB   glmmTMB 
#> -1716.218 -1716.218
sigma8 <- coef(fit8, "sigma")
abs(-2 * sigma8 - glmmTMB::fixef(cmp8)$disp)
#>           (Intercept) FoodTreatmentSatiated 
#>          1.153512e-06          1.315022e-06

Skip the conversion and the mismatch is large enough to notice — but it is signed the opposite way for nbinom2() than it was for the identity rule in Comparison 6, which is exactly the kind of error that survives a casual eyeball because it still lands in a plausible-looking range. Applied correctly, satiated broods have higher sigma, meaning lower size, meaning more overdispersion than food-deprived broods, and the two packages agree once, and only once, the conversion is applied.

Independence: WEAK.

What these eight comparisons show

Three of the comparisons above — the lme4 and metafor ones — are against separate estimation engines. Three are against glmmTMB, which is built on the same TMB automatic-differentiation stack and outer optimizer as drmTMB, so they check consistency between related implementations rather than across independent ones. Two are against ordinal, whose independence drmTMB’s design documents have not yet classified. It is worth naming where the distributional-parameter comparisons fall among those three groups: every comparison of a residual-scale submodel in this article — the capability drmTMB exists for — is in the glmmTMB group. The separate-engine checks on a sigma linear predictor are the two meta-analysis comparisons, where sigma is a between-study heterogeneity SD, not a residual scale.

These eight models were compared against packages listed in drmTMB’s DESCRIPTION under Suggests, because those are the packages this article can call. drmTMB implements a great deal more than eight models; what it implements, and what evidence stands behind each capability, is recorded in the package’s design documents and capability ledger rather than here.

Reproducing these numbers

Every figure above is a fitted result, so it is tied to the versions that produced it. Two comparisons are quoted to a precision — agreement to 1e-3, and a log-likelihood gap of 1.63e-11 — that a different optimizer build can move. The block below records what this rendering used, so a reader who gets a different number can tell whether the package changed or their environment did.

pkgs <- c(
  "drmTMB", "TMB", "Matrix", "lme4", "metafor", "metadat", "ordinal", "glmmTMB",
  "palmerpenguins"
)
data.frame(
  package = pkgs,
  version = vapply(
    pkgs,
    function(p) {
      if (requireNamespace(p, quietly = TRUE)) {
        as.character(utils::packageVersion(p))
      } else {
        NA_character_
      }
    },
    character(1)
  ),
  row.names = NULL
)
#>          package    version
#> 1         drmTMB      0.7.0
#> 2            TMB     1.9.21
#> 3         Matrix      1.7.5
#> 4           lme4      2.0.1
#> 5        metafor      5.0.1
#> 6        metadat      1.6.0
#> 7        ordinal 2025.12.29
#> 8        glmmTMB     1.1.14
#> 9 palmerpenguins      0.1.1
c(
  R = R.version.string,
  platform = R.version$platform,
  BLAS = basename(extSoftVersion()[["BLAS"]])
)
#>                              R                       platform 
#> "R version 4.6.0 (2026-04-24)"       "aarch64-apple-darwin23" 
#>                           BLAS 
#>             "libRblas.0.dylib"

No comparison on this page draws random numbers, so none of them carries a seed: each fit is a deterministic optimisation of a fixed dataset shipped by the package named beside it. The two quantities that do vary between machines are the optimiser’s stopping point and the BLAS in use, which is why both are recorded above rather than the seed that does not exist.