Changelog

Source: NEWS.md

drmTMB 0.7.1

Independent-evaluation fixes (wave 3)

• weights() composed with mi() now leaves the maximum-likelihood estimate unchanged under a constant reweighting for EVERY imputation family, not only the Bernoulli one repaired in wave 1. The eleven mi_family quadrature blocks in src/drmTMB.cpp (ordinal, categorical, beta, zero-one-beta, beta-binomial, Poisson, NB2, truncated NB2, lognormal, gamma, Tweedie) multiplied weights(i) into each quadrature node before logspace_add() combined them and then subtracted the combined log-density with no outer weight at all; the observed-row imputation-prior term was unweighted too. weights(i) now multiplies the combined mixture log-density and the observed-row prior, so weights = c is literal row duplication and the weighted objective is exactly c times the unweighted one (asserted through the fitted TMB objective in the new per-family test, which fails 43 times on the previous kernel). Three neighbours are recorded, not fixed. (1) The Gaussian latent mi() route (mi_family == 0, two blocks, src/drmTMB.cpp near lines 1215 and 4429) and the second-covariate has_mi2 prior (near line 1257) leave the covariate-model density unweighted: for OBSERVED rows that is the same one-line change as the template; for MISSING rows the covariate is a Laplace-integrated latent, and weighting both terms inside the integral is not the same as weighting the marginal, so that half needs its own design (the exact Gaussian marginal, then weighted, is the clean route). (2) The Tweedie imputation model sizes its fixed 35-node quadrature support from a start-value dispersion built with stats::var() (n - 1), so duplicating rows moves the support and the same integral evaluates about 1.3 nats apart; more importantly the support is frozen at the start value, so the integral is silently truncated whenever the fitted imputation scale outgrows it. Its duplication test arm is therefore not a negative control for Tweedie; the objective identity is. Credit: the independent evaluation by Russell Dinnage (finding M1, #1307) and the review that named the eleven sites.

• sigma(), predict(dpar = "sigma") (both type = "link" and "response"), residuals(), fitted() and simulate() now report the soft-clamped scale the likelihood actually evaluated for every family in drm_clamped_scale_families(). Previously the accessors returned exp(eta) on the raw linear predictor while the TMB objective used the clamped value, so on a fit where the clamp binds sigma() could exceed the fitted scale many-fold, a hand-recomputed log-likelihood under sigma(fit) missed logLik(fit) by hundreds of nats, and the Pearson residual spread read about two rather than one. One internal helper now applies the clamp at the two live prediction paths (predict.drmTMB() and the marginal predictor simulate() uses); the fix is a no-op whenever the clamp is inactive. The independent review of that change found two more raw-scale consumers, repaired in follow-up commits: predict_parameters() built its Wald interval from the raw predictor while its estimate was clamped, so 60 of 200 rows on a clamp-active fit showed an estimate outside its own interval; clamping the endpoints instead collapsed the interval to zero width at saturated rows (measured coverage 0), so clamp-bent rows now carry NA endpoints flagged conf.status = "clamp_limited", as profile() already does for direct SD targets, and in-band rows keep the ordinary Wald interval. That refusal selects on the estimate, not the truth: near the band edge the intervals that are returned are conditioned on the estimate having stayed in band, and the review measured their coverage falling from 0.985 to 0.882 as the true predictor approached the edge (0.98 unconditional), so treat any scale-parameter interval on a clamp-active fit as provisional; and predict(dpar = "sd(group)") for a modelled random-effect scale returned exp(raw eta) while the kernel and the fit’s own sdpars use the clamped value. Not yet aligned, and recorded as follow-ups: the response-scale profile() output for sigma, the Julia-bridge scale target, summary_parameter_delta_derivative(), and check_drm()’s clamp detector, which does not read a modelled sd(group) scale. Reporting the clamped value makes such a fit honest, not correct: when check_drm() says the clamp is active, rescale the response and refit. Credit: the independent evaluation by Russell Dinnage (finding M2, #1308).

• summary()$derived and the heritability()/icc()/repeatability() accessors now use the marginal residual variance when sigma carries random intercepts. drm_constant_residual_sigma() returned exp(b0), the median scale, so with a random intercept on log(sigma) of SD omega the variance-share denominator was too small by a factor exp(2 omega^2) (Russell measured +60.6% on the ratio at omega = 0.8). The residual variance is now exp(2 b0 + 2 sum(omega_k^2)), exactly E[sigma^2] for Gaussian random intercepts on the log scale, including a phylogenetic random intercept on sigma under the default unit-diagonal correlation; the delta-method standard error of the accessors carries the extra dependence on each omega_k (hand-derived and package SEs agree to 1e-11). That change corrects WHICH quantity the interval brackets; it does not make the bracket nominal: the accessors’ interval is still a Wald interval on a variance ratio, and the independent review measured 95% coverage of 0.910 with a random intercept on sigma and 0.928 for a sigma ~ 1 control over 500 replicates each, a shortfall that predates this change (method = "profile" where available is the safer choice). Random slopes on sigma, a structured effect whose correlation diagonal (measured on the rows the design uses, e.g. the tips of a tree, never the latent internal nodes) is not one, or a fit that carries a sigma random effect while the log(sigma) soft clamp bent the assembled predictor for at least one observation (the kernel clamps the assembled predictor, so the closed form is no longer the moment the likelihood uses; a clamp-active sigma ~ 1 fit still returns its constant clamped scale) give NA with a message naming why instead of the silently wrong median. summary()$derived$residual_sd is the marginal sqrt(E[sigma^2]) while the sigma parameter row remains the median exp(b0); they coincide only when sigma carries no random effect. Both repeatability loci are Gaussian-only, so the three-scale question of de Villemereuil et al.

  2016. does not arise for this number; the audit note docs/design/275-repeatability-scale-and-residual-variance.md records that non-Gaussian fits are refused by an error rather than mislabelled, and that latent-scale support for them is a feature decision. Credit: Russell Dinnage (finding S2, #1301) and the independent review that found the accessors rebuild the residual themselves.

• ?drm_phylo_penalty now describes the estimator the package actually reports. The prior on each phylogenetic SD is the documented exponential (PC) prior with P(sd > sd_u) = sd_alpha holding exactly; the penalty is evaluated on log(sd) with the change-of-variables Jacobian, so the reported penalised sd_phylo is the mode in log(sd), is never zero, and sits at 1/rate (0.334 at the defaults) under a flat likelihood. A penalised fit must not be used to test a null of no phylogenetic signal. No code or default changed (decision D-266). Credit: Russell Dinnage (finding S3, #1312) and the independent review that overturned the first proposed wording.

Independent-evaluation fixes (wave 2)

• simulate() now returns NA at masked missing-response rows for every family, matching residuals(). Twelve of thirteen families previously returned the internal missing-response sentinel (0, or 1 for the positive families) at those rows, so posterior-predictive tools such as DHARMa saw fabricated observations. The zero-inflated, hurdle and truncated count families were masked in a second commit after review; two tests that had pinned finite draws at masked rows now assert the masking invariant. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding M4.
• fitted_distribution()’s $p() and $d() recycle a scalar threshold across every row for the zero-inflated, hurdle and truncated count families; they previously returned row 1’s value only. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding S1.
• vcov(fit, type = "robust") (and robust = TRUE) now aborts with class drmTMB_vcov_robust_unsupported and says what to try instead; it used to return the model-based matrix silently. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding S4.
• AIC() and BIC() called with a drmTMB fit and a foreign model (for example an lm fit) return the standard one-row-per-model data frame; they used to drop the foreign model and return a bare scalar. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding S5.
• Finding M3 (the fixed_gradient row of check_drm() firing on correct fits at large n) no longer reproduces: the Newton polish added in #1130 already drives every correct fit’s gradient far below the tolerance. A regression test now locks that behaviour at n = 2000. Finding S3 (the phylogenetic SD penalty’s documented prior) stays open as a design decision after review; no change ships for it. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval).

Independent-evaluation fixes (wave 1)

• drm_logsigma_clamp_active() (and check_drm()’s logsigma_clamp_active row) now detects the LOWER log(sigma) clamp arm, not just the upper one. Previously a scale-model fit whose raw log(sigma) predictor sat on the lower bound (e.g. a response on a small numeric scale) was reported as “the clamp is not active”, even though the clamp had changed what the likelihood evaluated. The fit-time cli_warn() still fires only for the upper (runaway-scale) arm, since the lower arm is often a legitimate variance-zero boundary (meta-analysis tau = 0). Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding C1. check_drm() now reports the lower arm as a note whose text names the lower bound and the legitimate tau = 0 case, and the upper arm as a warning (Fisher review of the fix).
• drm_clamped_scale_families() (and hence check_drm()’s logsigma_clamp_active row) now names biv_lognormal and biv_student, not just biv_gaussian. src/drmTMB.cpp clamps log_sigma1/log_sigma2 identically for all three bivariate families in one shared branch, but the R-side list previously named only biv_gaussian, so check_drm() printed the false sentence “The log(sigma) clamp does not apply to this family” for the other two. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding Md-A.
• weights() composed with mi() (missing-predictor imputation) no longer moves the maximum-likelihood estimate under a constant reweighting. Ten duplicate call sites in src/drmTMB.cpp (one per response family sharing the Bernoulli-imputed, mi_family == 1 two-point mixture) multiplied weights(i) into each mixture leaf BEFORE logspace_add() combined them, computing log(p1*f1^w + p0*f0^w) instead of the correct w*log(p1*f1 + p0*f0); a constant weight therefore shifted the mi(x) coefficient. weights(i) now multiplies the combined mixture log-density and the (previously unweighted) imputation-prior term for observed rows, matching the general contract that weights = c is identical to literal row duplication. The drm_response_log_density() contract comment in src/drm_response_kernels.h is updated to spell out that “outside the leaf” means outside the whole mixture, not just outside each leaf call. The mi_family quadrature blocks for other imputed-covariate families (ordinal, categorical, beta, Poisson, lognormal, gamma, NB2, Tweedie, zero-one-beta, truncated-NB2, beta-binomial) share a structurally similar pattern and were fixed in wave 3 (above). Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding M1.
• mi() is now rejected with an error on every formula parameter except mu. The public mi() marker is an identity stub (function(x) x), and no non-mu formula path extracted or rejected it, so e.g. sigma ~ mi(z) was silently parsed as an ordinary covariate and gave a bit-identical logLik to sigma ~ z – no error, and no missing-data handling for z. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding Md-D.
• An unused factor level (for example, left over after subset() without droplevels()) no longer zeroes out every fixed-effect standard error. The empty level’s design-matrix column of all zeros made the fit’s Hessian singular, so estimates came back exactly right but every SE (not just the empty level’s) came back NA with a sdreport_non_pd_hessian status, and refitting with drm_control(se = TRUE) – the message check_drm() printed – was a no-op, since se = TRUE is already the default. drmTMB() now drops unused levels from the factor columns that enter a fixed-effect design matrix, once, up front; responses, mi() predictors and random-effect or structured-marker grouping variables keep their declared level sets, so an ordinal response with an empty category still reaches its own refusal. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding Md-E.
• check_drm()’s dropped_rows row now reflects rows the experimental MSPL estimator discarded because of a zero frequency weight, not just rows dropped by complete-case or known-covariance filtering. MSPL filters data before any family builder computes its keep vector, so a builder’s keep was always relative to the already-MSPL-filtered data – MSPL-discarded rows were invisible to it, and the row printed “no rows were dropped” even when MSPL had discarded some. mspl_frequency_rows$kept is now threaded through to re-express model$keep relative to the original input data. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding Md-N.
• skew_normal() now uses the package’s own tail-safe drm_log_pnorm() (already used by the binomial probit link) for its skew-CDF factor, instead of flooring pnorm(...) with + 1e-300. The floor saturated the far-tail log-density to a constant (log(1e-300) ~= -690.8) regardless of how far in the tail a point actually was, giving a gradient wrong by orders of magnitude and a plateau a maximiser could sit on. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding Md-M.
• The default: branch of drm_response_log_density() (src/drm_response_kernels.h) now calls error() instead of silently returning Type(0.0) (a likelihood contribution of 1) for an unhandled model_type. Every current call site is reachable-safe, so this is a guard against the next family wired into an mi() two-point sum before its case is added here, not a fix to an observable behaviour today. Credit: the independent evaluation by Russell Dinnage (rdinnager/drmTMB_eval), finding Mi-9.

Version bump only – tagging, release and CRAN submission remain the maintainer’s ceremonies. This heading summarizes, at a glance, the engine = "julia" bridge work already recorded in detail under the 0.7.0 section below:

• engine = "julia" masked-response fits: the is_converged() and bootstrap defects are fixed upstream in DRM.jl (#646).
• engine = "julia" bridge-side profile/bootstrap inference is qualified (G3) on two routes: base_gaussian_location_scale and plain_binomial_nonphylo.
• engine = "julia" admits beta_binomial() and cumulative_logit() on fixed-effect routes, and predict() on cumulative_logit() Julia fits now matches engine = "tmb".
• REML support is tabled by route and measured across both engines (#1142), and drm_julia_reml_supported() now also covers Location-Scale-Scale sd(...) models, promoting Capability Row 12 (location_scale_scale) from partial to covered (DRM.jl #558).

See the 0.7.0 section immediately below for the full, measured detail on every item above.

drmTMB 0.7.0

engine = "julia" control surface: no silent drops, boundary made permanent (leaf-engine-control-surface)

• drmTMB(..., engine = "julia") no longer silently ignores a drm_control() setting. Seven fields – se_report_covariance, se_skip_delta_method, se_group_sd, logsigma_clamp, logsigma_clamp_margin, newton_polish and fallback_optimizer – reached the bridge, were dropped, and the fit came back as if nothing had been asked for. Measured live at DRM.jl pin 430ef64cc: drmTMB(bf(y ~ x, sigma ~ x), data, engine = "julia", control = drm_control(newton_polish = FALSE)) returned a fit byte-identical to the default one (logLik -199.0299845089, same four coefficients) with no error. They now abort with an error naming the setting, like the nine settings that already did. The refused set is derived from drm_control() itself rather than hand-listed, so a control added to drm_control() is refused on the Julia path without a matching bridge edit; a totality test asserts the classification covers names(drm_control()) exactly.

• Behaviour change. A script that passed one of those seven settings with engine = "julia" used to run and now errors. It was never honoured, so no result changes – what changes is that the mismatch is now visible. Use engine = "tmb" for those controls.

• The forwarded set is documented as a closed whitelist in ?drm_control and the Julia-engine vignette: optimizer$g_tol and optimizer$algorithm on the base bridge, plus optimizer$q4_vcov on the bivariate q = 4 phylogenetic route, where optimizer$g_tol is forwarded as DRM.jl’s q4_g_tol and optimizer$algorithm is refused. Both admitted knobs were confirmed to reach the engine in the same run: g_tol = 1 moved the fit to logLik -199.0364652191, and algorithm = "em" reached DRM.jl and threw its own ArgumentError from gaussian_core.jl.

• The engine_control_surface capability row is no longer ambiguous. Its next_action used to read “Design engine_control explicitly before relaxing the gate”, advertising work that cannot be finished: most of drm_control() describes an nlminb()/TMB program that DRM.jl does not run, so there is nothing to forward it to and no native comparator against which a parity claim could be measured. The row now carries a permanent claim boundary naming every setting that does not cross and what to do instead, and claim_status stays experimental permanently – it records a deliberately narrow Julia-native surface, not unfinished work. ## Julia routes that refuse a whole control now name the offending settings (#1108)

• The four engine = "julia" routes that accept only a default drm_control() – structured, bivariate q2 structured, cross-family, and the joint missing-predictor adapter – refused without saying WHICH setting they refused. A caller who passed drm_control(newton_polish = FALSE, se_group_sd = TRUE, optimizer = list(iter.max = 500)) to a cross-family model was told only that “cross-family models currently accept only default control” and had to bisect their own call to find out why. Those routes were already fail-CLOSED – nothing was silently dropped on them – so this is a message defect, not a behaviour one, and the fitting behaviour is unchanged.

• The refusal now adds a line naming every non-default setting, derived from drm_control() itself (drm_julia_nondefault_control_fields()) rather than hand-listed, and reports optimizer entries the way a user wrote them (optimizer$iter.max). Measured on this branch, the same call now reports Non-default control settings: "se_group_sd", "newton_polish", and "optimizer$iter.max".

• optimizer_preset is a macro that also rewrites the nlminb iteration budgets, so it reports optimizer_preset, optimizer$iter.max and optimizer$eval.max together. That is measured behaviour and is pinned by a test rather than hidden.

• tests/testthat/test-julia-control-refusal-names.R enumerates names(drm_control()), so a new control field that the naming layer does not cover fails the suite instead of going unnamed.

engine = "julia" refuses a factor design it cannot reproduce, before Julia starts (DRM.jl #467, #609)

• The A6 guard already compared CODING SCHEMES (ordered factors, an explicit contrasts attribute, a non-default options("contrasts")). Two shapes agree on the scheme and disagree about the LEVELS it is applied to, so they went straight past it. Both were measured this run, driving 58 formula-construct cases end to end through engine = "julia" against engine = "tmb" and probing the marshalled columns in the live Julia session: 46 faithful, 10 refused, and 2 SILENTLY MISLABELLED.

• A character column whose R level order is not code-point order now refuses. R builds levels with sort() under the session collation, so c("a", "B", "c") gives levels a, B, c and codes against "a". A character column crosses the bridge as a plain Julia Vector{String} – no pool, no order – and DRM.jl sorts it by code point, giving levels B, a, c and coding against "B". Same column COUNT, so no count check fires. Measured at the standing pin DRM.jl 430ef64cc: the fit SUCCEEDED, both engines returned IDENTICAL coefficient names, and max|coef diff| was 0.1785 (n = 120) and 1.1036 (n = 150) – engine = "julia"’s intercept was the mean of the wrong baseline group. factor(<character column>) is the same hazard, because the bridge materialises the ORIGINAL values and sorts those. A FACTOR column is immune: its level order crosses intact in the CategoricalVector pool (probed for levels = c("c", "b", "a")), which is why the refusal tells you to store the column as a factor with an explicit level order.

• A factor level that no row uses now refuses. model.matrix() gives every DECLARED level a column, including an all-zero one; DRM.jl codes only the levels it observes, so it builds one column fewer. This was already fail-closed, but only deep inside Julia and with a message that named neither the column nor the fix (the R side must send exactly one name per column, plus a Julia stack trace). drmTMB now names the unused level and points at droplevels().

• Newer DRM.jl builds carry their own echo check (_bridge_check_coef_labels_fidelity) that refuses the first shape from the Julia side. That check does NOT exist at the pin above, which is why this guard is upstream of it: a refusal that lives only in the engine still lets an older build report the wrong parameter under the right name. The engine check remains the second line of defence, and DRM.jl’s count message gained the unused-level explanation in a companion change.

• Scope, measured, NOT derived: the other 46 cases in the battery (factor(), a bare factor column, reversed declared levels, mixed-case factor LEVELS, factor() over an integer column with 2/9/10 ordering, character columns whose order already agrees, logical columns, factor:factor and numeric:factor interactions, I(x^2), I(x*z), I(x^2 + z), poly(x, 2), poly(x, 3), poly(x, 2) * g, scale(x), (x + z)^2, (x + z + g)^3, x - 1, 0 + x, x:z, x + z + x:z - z, (x + z)^2 - x:z, log(x + 2), sqrt(x + 2), and the sigma-side repeats) were FAITHFUL: identical coefficient names and max|coef diff| <= 2.675e-10 against engine = "tmb". An independent design oracle went further and compared VALUES: for 37 constructs DRM.jl’s own design matrix is element-wise identical to R’s model.matrix() (31 at exactly 0; poly(x, 2) 8.33e-16, poly(x, 3) 8.26e-16, scale(x) 4.44e-16 – so DRM.jl reproduces R’s ORTHOGONAL raw = FALSE polynomial basis, #467’s flagged high-risk case). I(log(x + 2)), x * z - x:z (a - removal over an unexpanded *) and a single-level factor keep their existing refusals. This is one Gaussian location-scale fixture per construct, not a coverage study; the receipts are in docs/dev-log/evidence/julia-r-parity/formula-construct-fidelity/. ## Hurdle NB2 through engine = "julia": ONE call now fits on both engines

• drmTMB(bf(y ~ x, sigma ~ z, hu ~ w), family = truncated_nbinom2()) – drmTMB’s spelling of the hurdle negative binomial, since there is no hurdle_nbinom2() constructor – previously fitted on engine = "tmb" and ABORTED on engine = "julia", while the bridge instead accepted family = nbinom2() with hu, which the native engine refuses. A user could not switch engine = on one call. DRM.jl PR #662 makes TruncatedNegBinomial2() accept an hu part (delegating to its existing NegBinomial2 hurdle kernel), and the bridge fit now reports the native model_type "hurdle_nbinom2" instead of "truncated_nbinom2", so predict(fit, dpar = "hu") resolves its logit link the way it does natively. Measured on the native test suite’s own hurdle fixture (n = 1800, 486 zeros), both engines on the same call: coefficients agree to 9.4e-12 over 8 coefficients, logLik -2941.45558666655 vs -2941.45558666657, Wald SEs to 1.1e-06 relative, estimator "ML" on both. Fixed effects only – no phylogenetic, random-effect or structured hurdle route, and no interval claim. Known gap, declared not fixed: fitted()/residuals() on this route return the untruncated count mean through the bridge rather than the hurdle mean; every dpar is correct, so hurdle_nbinom2_mean() reproduces the native fitted() from them exactly. ## biv_lognormal() on the engine = "julia" fixed-effect route

• drmTMB(bf(mu1 = y1 ~ x, mu2 = y2 ~ x, sigma1 = ~ 1, sigma2 = ~ 1, rho12 = ~ 1), family = biv_lognormal(), engine = "julia") now fits instead of refusing with “currently supports Workflow G fixed-effect families”. One registry row (spec("biv_lognormal", fe = TRUE)) admits it, and the bridge’s default label branch now matches the biv_ prefix rather than the single literal "biv_gaussian" – without that, the registry row alone left BOTH the full and the short bf(mu1 = y1 ~ x, mu2 = y2 ~ x) form aborting inside DRM.jl (coef_labels supplies names for unknown dpar "sigma"), because a scalar sigma label was being sent to a model whose blocks are sigma1/sigma2/rho12. The scale contract is unchanged from the native engine and was measured, not assumed: both engines take y1/y2 on the raw positive scale, mu1/mu2 are means of log y, sigma1/sigma2 are SDs of log y, and both carry the change-of-variables Jacobian, so both match the independent raw-scale oracle in tests/testthat/test-biv-lognormal.R (logLik -164.673669588 on both, agreeing with the oracle to 4e-13). Measured same-target receipts on an n = 600 draw at DRM.jl pin 430ef64cc: coefficients 9.288083e-07 (7/7 name-matched), logLik -1124.208196846248 vs -1124.208196846242, Wald SE 1.707768e-06 relative. Scope: fixed-effect mu1/mu2 with intercept-only sigma1/sigma2/rho12 only – the same cell the native engine admits; no random effects, no phylogenetic or structured route, no interval-coverage claim.

• The admission needed a scope fence to hold that boundary, and this is the part worth reading twice. Bivariate families are exempt from the A4.G17 fixed-effect-only random-effect fence because biv_gaussian legitimately fits predictor-driven sigma1/sigma2/rho12 – so the new registry row inherited an exemption biv_lognormal has not earned. Measured at the pin before the fence existed: sigma1 = ~ x FIT through engine = "julia" (logLik -71.4056477) and so did rho12 = ~ x (logLik -70.64289338), while engine = "tmb" refused both. engine = "julia" now refuses a predictor on sigma1/sigma2/rho12 and an ordinary random-effect bar for this family on the R side, before Julia starts, in the same words the native engine uses. biv_gaussian() is untouched. ## Bootstrap replicates keep a masked fit’s response mask (#1188)

• confint(method = "bootstrap") on a fit made with missing = miss_control(response = "include") drew every replicate over the full design and refitted on all rows, regardless of how many rows the seed fit actually observed. A bootstrap whose replicates are richer than the original understates uncertainty, and the shortfall grew with the missing fraction. bootstrap_response_data() now re-applies the seed fit’s response NA mask to each simulated draw, and the replicate is refitted under the seed fit’s own response policy, so every replicate uses the same rows the seed fit did.

• Measured this run on y = 0.3 + 0.5x + N(0,1) exp(0.1x), n = 60, bf(y ~ x, sigma ~ x), Gaussian, target fixef:mu:x with truth 0.5; S = 200 datasets per cell, B = 99 replicates each, nominal 95% percentile interval, both arms on the same datasets and bootstrap seeds. Coverage at 10% / 30% / 50% of responses masked was 0.895 / 0.820 / 0.720 before and 0.910 / 0.895 / 0.910 after, against a Wald reference of 0.920 / 0.925 / 0.915 on the same datasets. Monte Carlo standard errors are 0.020-0.032; the paired improvement is significant at 30% masking (p = 2.75e-04) and 50% masking (p = 7.28e-12) but not at 10% (p = 0.375), so this release does not claim the 10% cell. Receipt: docs/dev-log/evidence/julia-r-parity/p2-g3/1188-bootstrap-mask-receipt.md.

• The default missing = miss_control(response = "drop") policy is untouched: a drop fit stores complete-cased data, so it never had a mask to lose and the restoration is a measured no-op there. Post-fix bootstrap coverage is still slightly below nominal on this n = 60 location-scale fixture, but so is Wald on the same datasets; that residual is small-sample behaviour, not a missing-data effect, and is not claimed fixed.

• The matching DRM.jl change (_bootstrap_data) lands separately, so engine = "julia" and engine = "tmb" keep the same replicate semantics. ## engine = "julia" admits REML on the bivariate q = 2 structured routes

• drmTMB(bf(mu1 = y1 ~ x + phylo(1 | p | id, tree = tree), mu2 = y2 ~ x + phylo(1 | p | id, tree = tree), sigma1 = ~ 1, sigma2 = ~ 1, rho12 = ~ 1), biv_gaussian(), engine = "julia", REML = TRUE) now fits instead of refusing, and likewise for matching relmat(1 | p | id, K = K) and spatial(1 | p | id, coords = coords) markers. The bridge previously refused this cell for EVERY provider – through two different branches, one per marker family – while DRM.jl fits it (fit_coevolution_q2_reml) and native engine = "tmb" has always admitted it. Measured at DRM.jl pin 430ef64cc against engine = "tmb" REML = TRUE on the same fixture: max |d coef| 5.15e-05 (phylo), 5.29e-05 (relmat), 4.20e-05 (spatial); |d logLik| 1.80e-04, 3.70e-07, 4.49e-08; estimator and DRM.jl’s own estim_method read REML on all three, with the ML and REML log-likelihoods about 6 units apart. phylo’s 1.80e-04 is the one number above 1e-4, and it is a property of that route rather than of REML: the already-shipped ML fit on the same fixture disagrees by 5.42e-04, three times wider.

• Point estimates only on these routes. No standard error, confidence interval, or coverage claim is made here. When the receipt was measured DRM.jl reported an all-NaN covariance for the bivariate q = 2 structured cell; DRM.jl has since started returning a finite one, confirmed live at pin 0edb916a5 alongside an exact reproduction of every point number above. The receipt stays point-only regardless: that covariance is the ML observed-information curvature evaluated at the REML point, which is not the quantity engine = "tmb" reports for a REML fit, and no comparison of the two has been measured. Use engine = "tmb" when you need uncertainty on this model.

• animal() q = 2 stays refused. DRM.jl fits it, but native engine = "tmb" still refuses bivariate animal() q = 2 REML, so there is no same-target comparator to measure against; the route is admitted only once that receipt exists. Every other bivariate structured shape keeps its existing refusal unchanged. Evidence: docs/dev-log/evidence/julia-r-parity/reml/biv-q2-bridge-receipt.md.

The formula-construct battery extended off its one Gaussian fixture, and a silently mislabelled sd(<group>) block found there (DRM.jl #467/#609/#730)

• PR #1227 ran 58 formula constructs through engine = "julia" on one Gaussian location-scale fixture and recorded two of them as silently mislabelled at the standing DRM.jl pin 430ef64cc: a character column whose R locale-collated level order is not Julia’s code-point order, and a factor level no row uses. Re-measured at drmTMB 2fcbb0fbf against a live DRM.jl (aee371cc9), neither is silent any more on a mu or sigma block – DRM.jl’s own _bridge_check_coef_labels_fidelity refuses both, on every family tried (poisson, nbinom2, binomial, gamma, cumulative_logit) and on both dpar sides, while a properly declared factor stays faithful to 1.05e-11 (poisson), 1.07e-11 (nbinom2) and 5.58e-12 (cumulative_logit). The silence was a property of the dead pin.
• It is not gone everywhere. That check iterates _bridge_rendered_regression_blocks, which skips every location-scale-scale sd_<group> / sdphy_<group> block by construction, so the group-level SD formula was echoed with R’s names and never compared against the design DRM.jl built. Measured at 2fcbb0fbf + aee371cc9: bf(y ~ x + (1 | study), sigma ~ z, sd(study) ~ s_chr) converged on both engines to an identical logLik (-69.917488, diff 2.98e-13) under identical coefficient names, with mu and sigma faithful to 2.12e-11 and the sd block off by 1.3853 – s_chrBeta reported as 0.692648 by engine = "tmb" and -0.692648 by engine = "julia", the baseline having moved from alpha to Beta. Declaring the same column as a factor in R makes the identical model faithful to 1.46e-10, which identifies the level order rather than the fit as the mechanism.
• tests/testthat/test-julia-formula-constructs-nongaussian.R pins the contract that matters here: a construct must be either FAITHFUL (identical base-R names AND coefficients within tolerance) or REFUSED, never the third thing – converging, returning identical names, and reporting a coefficient that means something other than its label says. It sets a Latin collation explicitly, because testthat’s default LC_COLLATE = "C" is exactly Julia’s code-point order and every case would otherwise skip while reading green. 31 pass / 0 fail / 0 skip live.
• Two independent fixes close the sd() cell and both were verified here: drmTMB’s own drm_julia_check_factor_level_fidelity() (PR #1227) refuses it before Julia starts – a route #1227 never claimed – and DRM.jl #730’s _bridge_check_lss_coef_labels_fidelity refuses it in the engine. NOT covered: the multi-IID sd route, sdphy_<group> live, and random-effect routes generally, where DRM.jl supplies no bridge_formula_labels_v1 at all and so refuses every factor, declared or not – honest, but a separate gap.

engine = "julia" default coefficient labels widened for the A4 family admissions

• The Julia bridge’s default coefficient labeller (drm_julia_bridge_default_dpar_labels()) previously defaulted a nu label only for student(); every other family with a native dpar beyond mu/sigma – tweedie()’s and skew_normal()’s nu, zero_one_beta()’s zoi/coi – was left unlabelled, so the SAME short formula forms engine = "tmb" accepts (a bare bf(y ~ x), or bf(y ~ x, sigma ~ z)) aborted at DRM.jl’s echo with coef_labels is missing an entry for dpar "nu". The defaulter now reads the extra dpars off each family’s own constructor. Measured live at DRM.jl pin 430ef64cc: tweedie()’s bare bf(y ~ x) and bf(y ~ x, sigma ~ z) now reach the SAME logLik as engine = "tmb" exactly (-479.8586143 and -463.2274318); a bare zero_one_beta() formula that omits zoi/coi now fits instead of aborting. A nu formula that native TMB refuses but DRM.jl fits (tweedie’s nu ~ z) is neither silently accepted nor pre-refused: engine = "julia" now emits one informational message naming the shape and the native refusal it has no cross-engine receipt against. This completes the fixed-effect (Workflow G) admission of the six A4 families landed across separate PRs: tweedie(), zero_one_beta(), beta_binomial(), truncated_nbinom2(), cumulative_logit(), and skew_normal().

engine = "julia" target discovery closes over what confint() accepts (#1156)

• profile_targets() on a Julia fit and confint() on the same fit disagreed in two ways that #1156’s own fix did not reach. Measured at DRM.jl pin 430ef64cc on one Gaussian 32-tip fit of bf(y ~ x + phylo(1 | species, tree), sigma ~ 1, sd(species, level = "phylogenetic") ~ z) fitted on both engines. (1) The two engines name the same estimand differently: native TMB reports fixef:sd_phylo(species):z, the bridge reports fixef:sd_phylo:z, the Wald intervals agree to six significant figures ([0.1548584, 0.5961240] vs [0.1548585, 0.5961239]), and each engine refused the other’s spelling – so a script could not be moved between engines by changing engine = alone. The native spelling is now documented as CANONICAL (docs/design/258-coefficient-naming-contract.md section 9) and the bridge ACCEPTS it, resolving it to the row it reports; the bridge’s own short form is still what coef(), vcov() and profile_targets() print. (2) A listed target was called unknown: sigma is listed by profile_targets() (profile_ready = FALSE, note missing_tmb_parameter) but confint(fit, "sigma", method = "profile") answered Unknown confidence-interval target: "sigma". It now says the row is listed, gives its inventory note, names the profile-ready alias to use instead, and lists the profile-ready targets. The union #1156 introduced was measured CLOSED at this pin (nothing confint() accepted was unlisted) and is now pinned by tests that fail if it regresses.

engine = "julia" admits REML = TRUE for the residual-only bivariate Gaussian cell

• drmTMB(bf(mu1 = y1 ~ x, mu2 = y2 ~ x, sigma1 = ~1, sigma2 = ~1, rho12 = ~1), family = biv_gaussian(), REML = TRUE, engine = "julia") now fits instead of refusing. Native engine = "tmb" has always fitted this cell by REML, and DRM.jl (PR itchyshin/DRM.jl#652) now fits the same restricted likelihood in closed form: both engines integrate out exactly beta_mu1 and beta_mu2 and both report the normalised Patterson-Thompson log-likelihood, so there is no data-independent constant to remove before comparing. Same-target on the committed fixture (n = 60, seed 1): logLik -97.021205818372 on both engines (difference 0.0), 7/7 coefficient names identical with max scaled difference 4.33e-07, Wald SE max relative difference 6.54e-07; the Julia fit’s estimator and DRM.jl’s estim_method both read "REML". The ML route is unchanged (-90.202703298791 on both engines). The rho12 max coefficient gap (4.33e-07) is entirely a link-guard convention difference, not model disagreement: native TMB bounds rho12 via 0.999999 * tanh(eta) (src/drmTMB.cpp) while DRM.jl uses 0.99999999 * tanh(eta) (src/sparse_aug_plsm.jl), so rho12:(Intercept) is not the same parameter on the two engines; on the natural (bounded) rho scale the two engines agree to 3.95e-12. The gap grows with |rho| and reaches 4.95e-04 at rho = 0.999, which would break the 1e-4 coefficient bar with both engines exactly correct (drmTMB#1190). This shape only: the gate requires all five bivariate dpars, no structured marker, no meta_V(), no random bar, and intercept-only sigma1, sigma2 and rho12. A covariate-carrying sigma/rho12 design keeps refusing even though DRM.jl’s closed form covers it, because nothing has measured it against a native comparator. Receipt: docs/dev-log/evidence/julia-r-parity/reml/reml-biv-residual-receipt.md; docs/design/261-reml-by-route.md row biv_gaussian_residual now reads FITS / FITS / FITS (drmTMB #1142, DRM.jl #624). ## Bridge-side profile and bootstrap inference qualified on the masked-response Julia route (#544)

• gaussian_response_mask is promoted partial -> supported on the r_bridge_status axis. On a Gaussian location-scale fit with missing = miss_control(response = "include"), confint() through engine = "julia" now agrees with engine = "tmb" on the same fit and the same target: Wald to 7.9e-08, profile to 5.1e-06/7.2e-06 (against a 1e-4 bar stated before measuring), and a 99-replicate parametric bootstrap completes with 0 failures on both engines. This closes the “bridge-side inference remains unqualified (G3)” fence for this route.

• This requires DRM.jl at or after #646. Against an older DRM.jl the same fit reports is_converged() FALSE and its bootstrap loses every replicate; both were DRM.jl defects, not drmTMB ones.

• Disclosed, and it affects engine = "tmb" equally: a parametric bootstrap on a masked-response fit draws each replicate response over the FULL design and refits on every row, so the interval is calibrated to the complete-data sample size and is narrower than the observed-data Wald interval by an amount that grows with the missing fraction (measured bootstrap/Wald width ratio at 10%/30%/50% masked: 0.82/0.67/0.65 on engine = "julia", 0.83/0.68/0.68 on engine = "tmb"). If you bootstrap a heavily masked fit on either engine, treat the interval as anti-conservative. A cross-engine fix is tracked as #1188. ## Ordered cutpoints through engine = "julia": discoverable, and refused by name (#1144)

• #1144 polished the constrained stats::nlminb() solve behind the native ordinal cutpoint profile, so engine = "tmb" reports honest "ordinal:cutpoint:<label>" intervals. cumulative_logit() is an admitted bridge family, so the same question needed an answer on the other engine. Measured on the committed fixture at DRM.jl pin 430ef64cc: the Julia fit’s fit$ordinal$cutpoints already agreed with the native slot to 8.98e-13, but profile_targets() listed no ordinal row at all, and confint(fit, parm = "ordinal:cutpoint:low|medium") answered Unknown confidence-interval target for method = "wald", "profile" and "bootstrap" alike – a typo diagnosis for a target plainly on the fit. DRM.jl’s bridge inference accepts only fixed-effect and random-effect-SD targets at that pin, so the cutpoint interval genuinely cannot be routed. profile_targets() now lists the ordinal rows for a Julia-engine cumulative_logit() fit – the same names, in the same order, as the native fit – with the public cutpoint rows profile_ready = FALSE and note julia_ordinal_cutpoint_native_only, and confint() refuses a cutpoint target for every method, naming the cutpoint, engine = "tmb", and fit$ordinal$cutpoints. No numbers changed on either engine: this is a discovery-and-diagnosis fix plus an explicit engine boundary. ## Same-target REML receipts for three capabilities the parity scoreboard read UNCITED (#1142)

• docs/design/parity-scoreboard.md reported 23 of 45 drmTMB-native capabilities UNCITED on its bridge axis – reachable through engine = "julia", but with no receipt a reader could check. Three were REML cells, and no REML row existed in any of DRM.jl’s four evidence tables at all. All three now carry a same-target comparison against engine = "tmb", measured at DRM.jl aee371cc9: the fixed-effect Gaussian location-scale cell (coefficients 9.47e-12, REML logLik 7.11e-14), the Gaussian mean ordinary random intercept (1 | g) (coefficients 9.36e-11, logLik 1.07e-12, SEs 3.33e-07 relative), and the bivariate dense q4 phylogenetic location-scale cell (coefficients 7.47e-04, logLik 2.02e-02, inside that row’s own recorded atol_loglik of 0.03). Estimator honesty is read from both sides on every fit – drmTMB’s fit$estimator and DRM.jl’s own estim_method – and the ML/REML logLik gap is asserted non-zero, so an ML fit wearing a REML label would fail the tests rather than pass them. engine = "julia" gains no new route here; only evidence.

• Two boundaries are recorded rather than papered over. (a) On the heteroscedastic fixed-effect cell the two mean-block Wald SEs differ by 3.09e-03 relative, past the 1e-3 bar, while the scale block agrees to 3.1e-07; a control refitting the same data with sigma ~ 1 makes the mean block agree exactly, so the gap is carried by the sigma covariate, not by REML. (b) A block-diagonal q4 call – two distinct phylo labels, one on the means and one on the scales – is silently fitted as the dense q4 model by engine = "julia": on one fixture the native block-diagonal REML fit gives logLik -934.738 at df = 11 while the bridge returns -930.165 at df = 15, the dense answer. Use engine = "tmb" for a block-diagonal q4 REML fit until the bridge refuses that layout. Both are in docs/design/261-reml-by-route.md, and the ledger note claiming the engine = "julia" path for the q4 cell was “halted by design” is corrected – it fits.

biv_student() admitted through engine = "julia" (leaf fam-biv-student)

• drmTMB(bf(mu1 = y1 ~ x, mu2 = y2 ~ x, sigma1 = ~1, sigma2 = ~1, nu = ~1, rho12 = ~1), family = biv_student(), engine = "julia") now fits instead of refusing. DRM.jl needed no change: at pin 430ef64cc its _bridge_family("biv_student") already routed the tag to the bivariate Student-t model in src/bivariate_student.jl, so the admission is one row in the Julia family registry plus the retirement of a family-specific abort in drmTMB() that fired before the registry was ever consulted. Measured this run against engine = "tmb" on the tests/testthat/test-biv-student.R draw (n = 400, seed 6401): coefficients agree to 3.771e-07 (8/8 matched by name), logLik -928.707976 on both engines (diff 2.569e-11), and per-coefficient Wald standard errors agree to 9.013e-07 relative, with the comparator’s own negative control still failing as it must. sigma1/sigma2 are Student-t scales (marginal SD = sigma * sqrt(nu / (nu - 2))), nu is one shared degrees-of-freedom parameter, and rho12 is the scatter correlation – zero rho12 is not independence at finite nu – identically on both engines.
• The Julia route is fenced to exactly the shape native engine = "tmb" fits. Random-effect bars and non-intercept sigma1, sigma2, nu or rho12 formulas are refused with the native wording before Julia starts, and confint() is deferred for this family on both engines as it always was natively. Each of these fitted or returned intervals through engine = "julia" before this release while engine = "tmb" refused them; a shape the native engine refuses has no same-target comparator and so can carry no parity receipt. No interval-coverage claim is made.

predict(type = "quantile") now works through engine = "julia" (#1198)

• predict() on an engine = "julia" fit accepted type = c("response", "link") only, so predict(fit, type = "quantile") was refused by match.arg() for EVERY family while engine = "tmb" has accepted it since the distributional-output layer landed. The bridge method now takes type = "quantile" and a prob argument, and hands the fit to the SAME native quantile code the TMB engine uses: the per-family density/CDF/quantile registry, the prob validation, the percentage column labels, and the calibrated/prob/label attributes all keep one source of truth, with only the per-row distributional parameters coming from the bridge’s own reconstruction. All 14 families the bridge admits on its fixed-effect route qualify: gaussian, biv_gaussian, student, lognormal, poisson, nbinom2, gamma, beta, binomial, truncated_nbinom2, zero_one_beta, tweedie, beta_binomial and cumulative_logit. Measured against engine = "tmb" on the same fit at DRM.jl pin 430ef64cc (tests/testthat/test-julia-predict-quantile.R, prob = 0.1 / 0.5 / 0.9), identical dimensions, column names and attributes, and max |d quantile| at most 4.93e-11 on stored rows and 2.34e-11 on fresh newdata – exactly 0 for the five discrete families. NOT covered: a meta_V() fit refuses type = "quantile" on this engine, because a Julia-bridge fit does not retain the per-row known sampling variance and its quantiles would silently use sigma alone (measured 3.485e-01 too narrow against engine = "tmb"); use engine = "tmb" there. This is a distributional plug-in interval on both engines – attr(., "calibrated") is FALSE – and it makes no interval-coverage claim. ## Zero-inflated Poisson through engine = "julia": focused tests, and a corrected registry note (leaf-fam-zi-poisson)

• The zi_poisson capability row has carried a same-target parity receipt since 2026-09-05, but no test drove a zi ~ formula part through engine = "julia" – measured that day, of the tests/testthat files matching engine = "julia", none also matched a zi ~ formula. Of design 168’s four limbs of “covered”, the focused-test limb was missing. tests/testthat/test-julia-zi-poisson.R now supplies it ([ FAIL 0 | WARN 0 | SKIP 0 | PASS 39 ], 1 live test): the admission mechanism, the payload and coefficient labels for both dpar blocks, the live same-target round trip, and a plain-Poisson contrast on the same data proving the zi part reaches the engine rather than being dropped in marshalling. Re-measured in that run on its own fixture (n = 600, seed 20260905, bf(count ~ x, zi ~ z), DRM.jl pin 430ef64cc, comparison code from DRM.jl tools/parity_numeric.R and tools/parity_se.R): coefficients 5.776934e-12 (4/4 matched by name), logLik -786.1016601045 on both engines (diff 1.136868e-12), Wald SE 1.317571e-08 absolute / 2.373401e-07 relative over 4 SEs, negative control NEGATIVE_CONTROL_OK at rel 9.090911e-02, estimator ML on both. This closes a limb; it does not move r_bridge_status, which stays partial pending a bridge-side inference (G3) receipt.

• R/julia-family-registry.R listed zi_poisson, zi_nbinom2 and hurdle_nbinom2 among families “the Julia bridge has NO case yet” for. That was wrong on both halves, and the note is corrected. They are model_type values, not family_type values: drm_family_type() never returns any of them, a zero-inflated Poisson is spelled family = poisson() plus a zi ~ formula part, and the bridge already routes all three through the poisson / nbinom2 registry rows plus the zi / hu entries in julia_bridge_supported_dpars(). No registry row is added: one keyed on those names would admit a family tag drmTMB never emits. ## Bivariate animal() q2 REML admitted (leaf-biv-animal-reml)

• drmTMB(bf(mu1 = y1 ~ x1 + animal(1 | p | id, A = A), mu2 = y2 ~ x2 + animal(1 | p | id, A = A), sigma1 = ~1, sigma2 = ~1, rho12 = ~1), family = biv_gaussian(), REML = TRUE) now fits instead of refusing. The bivariate q2 exact-covariance REML route was already provider-agnostic in TMB’s C++ (no branch on marker identity) and in DRM.jl’s own REML implementation (animal and relmat markers call the identical make_coevo_problem_from_covariance() path, bit-identical loglik measured on a matched fixture); the prior refusal was a recorded scope decision, not a mathematical one. Measured this run: the new animal() REML fit agrees with the relmat() control on the same matrix to machine precision (identical mathematics), with the dense restricted-likelihood oracle, and with a direct DRM.jl method = :REML call (|d logLik| < 1e-4, coefficients to 1e-7); see docs/dev-log/evidence/julia-r-parity/reml/biv-animal-q2-receipt.md. That DRM.jl comparison is point-estimate and logLik only: DRM.jl’s bivariate q2 structured route returns vcov = NaN for every provider and both estimators, so there is no Julia SE to compare. engine = "julia" itself is untouched by this change and still refuses REML = TRUE on every bivariate q2 structured route (phylo, spatial, relmat, animal alike); the fit above is native TMB (engine = "tmb") only. Precision (Ainv) and pedigree-built animal representations, slopes, q4+, and scale-side bivariate relmat/animal REML routes remain deferred and refused. ## engine = "julia" zero-inflated NB2: fitted() and residuals() now agree with the native engine (DRM.jl bridge fix)

• A zero-inflated negative-binomial model (family = nbinom2() with a zi ~ formula part, native model_type "zi_nbinom2") already routed through engine = "julia", but nothing tested or documented it – and it carried a silent disagreement. fitted() and residuals() differed between the two engines on the SAME converged fit, because DRM.jl’s fitted() is the count-component mean mu while drmTMB’s is the unconditional mean (1 - zi) * mu. Measured on the package’s own tests/testthat/test-zi-nbinom2.R fixture (n = 1800): coefficients agreed to 4.56745752330789e-13 and logLik to 1.72803993336856e-11, while fitted() disagreed by 1.3665229755584 – invisible to any coefficient or likelihood check. Fixed on the DRM.jl bridge boundary (_bridge_fitted_marginal), which leaves DRM.jl’s own fitted/simulate/marginal_parameters and the bridge’s mu dpar untouched; after the fix the difference is 4.83169060316868e-13. The same repair covers zero-inflated Poisson; hurdle fits are deliberately not repaired (their mean also divides by 1 - P(0) and has no bridge receipt). A new tests/testthat/test-julia-family-zi_nbinom2.R pins the route: coef, logLik and Wald-SE parity, the coefficient-naming contract for all three dpar blocks (design 258 section 8.10), the engine’s own refusals, and two KNOWN GAPS that are pinned rather than claimed – predict(dpar = "zi") is refused on a Julia-engine fit (the bridge tags it "nbinom2", whose link table has no zi row), and sigma() returns a list rather than the numeric vector the native engine returns. No registry row and no family admission: the route was already reachable through the existing nbinom2 row.

engine = "julia" scope fence for the fixed-effect-only family cohort (A4.G17)

• Admitting a family on the engine = "julia" fixed-effect route (fe = TRUE in the family registry, R/julia-family-registry.R) had a gap: nothing pre-Julia checked an ordinary random-effect bar ((1 | g), including random slopes) or an sd()/sd_phylo() scale submodel against that family for a family with none of phylo_only/locscale_phylo/slope_phylo/structured set. Measured on origin/main with DRM_JL_PATH unset: student() y ~ x + (1 | g), beta_binomial() cbind(s, f) ~ x + (1 | g), lognormal() sigma ~ z + (1 | g), zero_one_beta() y ~ x + (1 | g), and truncated_nbinom2() y ~ x + (1 + x | g) all reached drm_julia_setup() – five of six fe-only families with no receipt. A new drm_julia_refuse_fe_only_random_effects() gate, registry-driven (no hard-coded family list, so a later fe-only admission is covered without an edit here), now refuses these before Julia starts, naming the family, the fe registry column that admitted it, and the offending term. phylo() and relmat()/animal()/spatial() markers are untouched by this gate – both are already refused pre-Julia by existing gates. No family admission or removal, no DRM.jl change, no comparison-row change: this closes a scope gap only. ## Profile and bootstrap intervals on the residual bivariate Gaussian Julia route (#544)

• confint(fit, parm = "fixef:mu1:x", method = "profile" | "bootstrap") now works for a residual-only bivariate Gaussian fit made with engine = "julia" – bf(mu1 = y1 ~ x, mu2 = y2 ~ x, sigma1 = ~ 1, sigma2 = ~ 1, rho12 = ~ 1), family = biv_gaussian(). Previously profile_targets() reported profile_ready = FALSE for EVERY row of EVERY bivariate fit, so no parameter of this route – including rho12 – had a profile or bootstrap target, and the call was refused R-side before any Julia call. The bivariate rule is now a SPLIT rather than a blanket refusal: a bivariate fit carrying a covariance provider (a phylo tree, or a K / A / coords matrix) is the structured route, whose inference target is the four among-axis SDs, and its fixed-effect rows stay not-ready; a bivariate fit carrying none is this residual route and is ready on the same precondition every univariate route uses. Needs DRM.jl PR #647 for the bootstrap half; profile needed no engine change. Measured on the committed gaussian-bivariate-rho12 fixture (n = 180): 7 profile-ready fixed-effect targets where there were none, all seven profiling to a finite interval; same-target agreement against engine = "tmb" of Wald 3.4e-14 / 3.0e-07 and profile 2.1e-06 / 6.5e-06 on fixef:mu1:x and fixef:rho12:(Intercept); bootstrap R = 99 with 0/99 failed on both engines. The capability row biv_gaussian_residual moves partial -> supported on the r_bridge_status axis. NOT an interval-coverage claim (one fixture, one seed), and the bootstrap comparison is distributional overlap only – the two engines draw from independent RNG streams, so the same seed does not reproduce the same replicates. The fixef:rho12:(Intercept) agreement is NOT solver agreement of the kind fixef:mu1:x shows: TMB guards rho12 with 0.999999 and DRM.jl with 0.99999999, a deterministic guard-constant reparameterisation of about 9.9e-07 * rho / (1 - rho^2) on the linear-predictor coefficient, predicted at 3.77138e-07 and measured at 3.77151e-07 for this fixture’s |rho12| ~ 0.34. The 1e-6 Wald bar on rho12 is therefore CONDITIONAL on this fixture’s weak correlation, not a general parity claim: the offset crosses 1e-6 near |rho12| ~ 0.62 and reaches ~4.9e-05 at rho12 = 0.99. Tracked cross-engine at itchyshin/drmTMB#1190; aligning the three guard constants (TMB, DRM.jl, and the bridge’s own atanh(rho) back-transform in drm_julia_residual_rho12_corpair()) changes numerics on every bivariate receipt and is a deliberate cross-engine decision, not part of this leaf. Receipts under docs/dev-log/evidence/julia-r-parity/p2-g3/. ## engine = "julia" masked-response fits: convergence flag and bootstrap fixed upstream (DRM.jl #646)

• A Gaussian fit with missing = miss_control(response = "include") through engine = "julia" reported is_converged() FALSE and failed all bootstrap replicates, even though its coefficients and profile interval agreed with engine = "tmb" to ~7e-06 on the same fixture. Both defects were DRM.jl’s and are fixed there (DRM.jl #646): the degeneracy check took std() of a response vector still carrying NaN in the masked rows (so its bar was NaN and every comparison against it false), and the replicate simulator drew fit.nobs values against full-design mean and scale vectors, throwing DimensionMismatch on every replicate. MEASURED on the fixture: raw Optim.converged TRUE, |grad|inf 6.41e-12 against a g_tol of 1e-8, and a parameter vector bit-identical to the complete-case fit – a wrong return-code mapping, not a hard optimisation surface.

• fit$opt for a Julia-engine fit is no longer a bare list(convergence = <0/1>). It now also carries iterations (the optimiser’s own count, NA_integer_ on a route that does not record one – never 0) and a message describing what DRM.jl reported. DRM.jl sends no optimiser message string, so message is composed from the facts that did cross the bridge and says so; it is not presented as a verbatim optimiser message.

• tests/testthat/test-julia-missing.R gains a live assertion block on the masked fixture: is_converged() TRUE, opt$convergence 0, opt carrying iterations and message, and confint(method = "bootstrap", R = 19) with 0 failed replicates. Receipts, including the per-hunk red controls and the re-qualification against the fixed DRM.jl, are under docs/dev-log/evidence/julia-r-parity/p2-g3/.

• NOT promoted here. The gaussian_response_mask capability row stays partial: the two defects that blocked its G3 qualification are fixed and re-measured (wald delta 7.86e-08, profile delta 7.18e-06, bootstrap 0/99 failed on both engines), but the status move belongs on top of #1184, which is concurrently rewriting that row and the guard that pins it.

• MEASURED: the parametric bootstrap on a masked-response fit is anti-conservative, and the narrowing GROWS with the missing fraction – masking 10% / 30% / 50% of the rows gave relative narrowing of roughly -7% / -27% / -40% against the seed fit’s Wald SE, and an implied nominal-95 interval covering roughly 93% / 85% / 76% of the time. The narrowing grows with the missing fraction because replicates are drawn over the FULL design regardless of how many rows were observed; this is shared by both engines. The fix (a mask-preserving bootstrap) is a cross-engine statistical calibration change, out of scope here, and tracked as #1188. ## engine = "julia" bridge-side profile/bootstrap inference qualified on two routes (G3)

• base_gaussian_location_scale and plain_binomial_nonphylo promoted r_bridge_status partial -> supported: profile and bootstrap confidence intervals through engine = "julia" on a real fixed-effect target agree with engine = "tmb" within 1e-4 (profile deltas measured 2.8e-06/5.9e-07 and 9.4e-08/2.3e-06), both engines converge, and bootstrap (R = 99) intervals overlap with 0/99 failed replicates on either side. plain_binomial_nonphylo’s bootstrap is enabled by drmTMB#1123 (already fixed). Two related rows do NOT promote: biv_gaussian_residual has no profile/bootstrap-ready target on the Julia bridge for any parameter (a structural gap), and gaussian_response_mask’s Julia bootstrap fails all 99 replicates while its underlying fit’s own optimizer convergence flag reads FALSE on the missing-response fixture – both newly measured, both left unqualified rather than rounded up. A purpose-built quasi-complete-separation binomial cell exercises DRM.jl’s #631 profile-endpoint-failure backstop through the public confint() entry point for the first time: engine = "julia" refuses rather than returning an infinite bound (docs/dev-log/evidence/julia-r-parity/p2-g3/). ## engine = "julia" fits the mean-only phylogenetic Gaussian cell by REML (#1142)

• drmTMB(bf(y ~ x + phylo(1 | species, tree = tree), sigma ~ 1), family = gaussian(), REML = TRUE, engine = "julia") now fits instead of refusing. DRM.jl #624 item (c) is closed on the engine side, so the bridge’s REML gate admits this one cell. Both engines integrate out the same set – the phylogenetic field and beta_mu – and both carry the same additive 0.5 * p_mu * log(2 * pi) from the Patterson-Thompson restriction, so the two REML log-likelihoods are directly comparable with no constant removed: drmTMB’s REML construction folds beta_mu into TMB’s random= vector, and DRM.jl profiles it out exactly by GLS and adds 0.5 * logdet(Xmu' V^-1 Xmu) - 0.5 * p_mu * log(2 * pi). Measured on the committed fixture (tests/testthat/test-reml-phylo-location.R draw, n = 90, 30 tips) against native engine = "tmb" with REML = TRUE: logLik -76.000977125761 vs -76.000977125105 (6.56e-10), max scaled coefficient difference 3.71e-08, identical coefficient names, identical nobs (90) and df (4); the fit’s estimator is "REML" and equals what DRM.jl reports as estim_method. One documented convention difference: drmTMB’s REML standard errors for beta_mu come from TMB’s sdreport over a random set that contains it, so they propagate variance-parameter uncertainty, while DRM.jl reports the canonical (Xmu' Vhat^-1 Xmu)^-1 – a hand-built GLS oracle confirms DRM.jl to 4.5e-08, and the two engines’ SEs differ by 4.84e-05, 1.50e-03 and 1.42e-07 relative on the three coefficients. This is a construction difference, not an error on either side: the test asserts DRM.jl-vs-oracle agreement at rtol ~5e-7 and bounds the cross-engine gap at rtol 2.5e-3, not a bare rtol-1e-3 cross-engine match. Tracked in #1201. Scope: phylo() with an intercept-only sigma only. A sigma predictor, an ordinary (1 | g) bar alongside the phylo term, relmat() / animal() / spatial(), and a non-default missing response engine all still refuse before Julia starts, with drmTMB’s own message. Receipt: docs/dev-log/evidence/julia-r-parity/reml/reml-phylo-mean-receipt.md; docs/design/261-reml-by-route.md row gaussian_phylo_mean now reads FITS / FITS / FITS. fit$bridge$gradient is deliberately absent on a Julia REML fit of this cell (it is present, and near zero, on the ML fit): the route’s analytic score belongs to the ML marginal and is (1.01, 0.99) on the variance parameters at the REML optimum, so reporting it would have read as “not converged” for a converged fit.

engine = "julia" admits beta_binomial() (fixed effects)

• drmTMB(bf(cbind(successes, failures) ~ x, sigma ~ z), family = beta_binomial(), engine = "julia") now routes to DRM.jl’s BetaBinomial instead of refusing with the Workflow G message. Both engines fit the same target with the same sigma parameterisation (precision phi = 1 / sigma^2); the per-row denominator crosses as trials context, not as a distributional parameter. Measured at DRM.jl pin 430ef64cc on the committed fixture (tests/testthat/test-beta-binomial.R draw, n = 1200): max |d coef| 7.77e-15, |d logLik| 2.57e-10, per-coefficient Wald SE max relative difference 1.73e-07; the fit’s estimator ("ML") equals what DRM.jl reports as estim_method. Fixed effects only: a phylo() term with this family is still refused before Julia starts, and no (1 | g) route, structured marker, or interval-coverage claim is made; use engine = "tmb" for those.

engine = "julia" admits cumulative_logit() on the fixed-effect route

• cumulative_logit() (dpar mu; ordered-factor response) now routes through engine = "julia" for fixed-effect models – one row in the Julia family registry plus a small family file that integer-codes the ordered response, sends mu without its intercept, labels DRM.jl’s cutpoints block with drmTMB’s own "low|medium" spelling, and moves that block into fit$ordinal (the native engine’s slot; cutpoints never appear in coef() or vcov()). Same target as engine = "tmb" on the committed test-cumulative-logit.R fixture at DRM.jl pin 430ef64cc: max |d coef| 2.20e-14, |d logLik| 1.13e-11, cutpoints within 8.98e-13, Wald SE within 5.76e-09 relative; estimator reads "ML" and equals DRM.jl’s estim_method. A sigma ~ formula and an unordered factor are refused before Julia starts with the native engine’s own reasons. Known gap, not fixed here: predict() on the Julia object aborts for this family (its design rebuild keeps the mu intercept); fitted() is unaffected. No phylogenetic, structured, or random-effect ordinal routes are admitted; (1 | g) fails closed at DRM.jl’s label echo.

engine = "julia" admits skew_normal() on the fixed-effect route

• skew_normal() (dpars mu, sigma, nu) now routes through engine = "julia" for fixed-effect models – one row in the Julia family registry on the R side, plus a five-line _bridge_family() case in DRM.jl (SkewNormal() existed there but the R bridge had no tag for it; that case is DRM.jl PR #641, so a DRM.jl checkout without it still aborts at the Julia boundary with drm_bridge: unsupported family). Same target as engine = "tmb" on the committed test-skew-normal-location-scale.R fixture at DRM.jl pin 430ef64cc + that case: max |d coef| 1.89e-11, |d logLik| 2.16e-12, per-coefficient Wald SE within 1.04e-06 relative; estimator reads "ML" and equals DRM.jl’s estim_method. The public moment parameterisation (mu = mean, sigma = SD, nu = Azzalini slant) is the same on both engines, so the bridged coefficients are the native ones; a predictor-dependent nu ~ z agrees to the same tolerance. Write

predict() on cumulative_logit() Julia fits now matches engine = "tmb"

• Closes the gap the cumulative_logit() admission above recorded: predict() on an engine = "julia" cumulative_logit() fit no longer aborts. The reconstructed mu design dropped the fitted coefficient block’s intercept but kept restoring “(Intercept)” via stats::model.matrix(); the design rebuild now drops it too, scoped to cumulative_logit()’s mu dpar so no other family’s Julia prediction path is touched. mu’s link is identity, so type = "response" and type = "link" are identical on both engines. Measured on the committed tests/testthat/test-julia-cumlogit-predict.R fixture (n = 900) at DRM.jl pin 430ef64cc: max |d prediction| for stored data and fresh newdata, both types, all below 1.5e-13. Thresholds are not read by predict() on either engine for this family – they live in fit$ordinal, not in the linear predictor. type = "quantile" on engine = "julia" fits was unavailable when this landed; it is added by the entry at the top of this file (drmTMB#1198).

REML support tabled by route, measured across both engines (#1142)

• docs/design/261-reml-by-route.md, generated by tools/write-reml-route-table.R, tables REML support for every engine = "julia"-admitted route (the 21-row capability TSV plus A5’s three ordinary-random-effect shapes) across native engine = "tmb", native DRM.jl, and the bridge – 30 rows, each FITS/REFUSES/N/A with a measured or cited receipt, no cell says “expected”. No REML behaviour changed: this is measurement only. Surfaces two new capability gaps where native TMB is ahead of DRM.jl (plain fixed-effect bivariate Gaussian REML, and mean-only Gaussian relmat() REML), one new bridge under-admission (DRM.jl’s own #450 Poisson+relmat() Cox-Reid REML is unreachable through engine = "julia" because the structured-term gate fires before the Poisson-specific gate), and one four-cell honesty-of-interface gap (fe_poisson, zi_poisson, Poisson+relmat(), and Gaussian random-slope REML requests reach a raw Julia stack trace instead of the bridge’s intended refusal message). None of these are fixed here; see the table’s Gaps section.

Coevolution accessors ported from DRM.jl (#1118)

• New coevolution_cor(), coevolution_vc(), and coevolution_summary() read the among-axis structure of a q = 4 structured bivariate location-scale fit – the 4 x 4 correlation matrix over mu1, mu2, sigma1, sigma2, the per-axis variance components, and a tidy long form over the six axis pairs – from the covariance the fit already stores, for both the native engine and engine = "julia" (whose raw-branch-length Sigma_a is rescaled to the unit-height convention, as profile_targets() already does). Ported term-for-term from DRM.jl src/coevo_accessors.jl and checked live against DRM.jl’s own accessors at pin 430ef64cc. ## Boundary-corrected likelihood-ratio test for variance components (DRM.jl #1116 port)

• New chibar_pvalue(statistic, q) and lrt_boundary(full, reduced, q): the chi-bar-square boundary-corrected likelihood-ratio test for q (1 or

  2. variance components tested at zero, ported term-for-term from DRM.jl src/chibar.jl (pin 430ef64cc). Dropping a random effect tests a variance on the boundary of its space, where the naive chisq(q) p-value is conservative; lrt_boundary() reports both the corrected and the naive p-value, refuses ML-vs-REML pairs, REML pairs with different mean structures, and penalized (MAP) fits. Measured against DRM.jl’s native lrt_boundary on the same data (four fixtures, q = 1 and q = 2): the statistic agrees to 3.8e-10 (strong random intercept, stat 670.27), 2.2e-11 (moderate, stat 2.2012, p 0.06895 vs 0.06895 to 9.7e-13) and 4.5e-11 (two crossed intercepts, q = 2, stat 94.418); at a variance estimated at zero both engines return the boundary point mass 0.5 (to 2.3e-7). chibar_pvalue() agrees with DRM.chibar_pvalue to 1e-12 relative on identical statistics (tests/testthat/test-lrt-boundary.R). ## engine = "julia" admits tweedie() on the fixed-effect route

• tweedie() (dpars mu, sigma, nu) now routes through engine = "julia" for fixed-effect models – one row in the Julia family registry, no bridge code. Same target as engine = "tmb" on the committed test-tweedie-location-scale.R fixture at DRM.jl pin 430ef64cc: max |d coef| 2.77e-11, |d logLik| 0, per-coefficient Wald SE within 3.28e-06 relative; estimator reads "ML" and equals DRM.jl’s estim_method. Write nu ~ 1 explicitly for now: a formula that omits nu aborts at DRM.jl’s label echo (coef_labels is missing an entry for dpar "nu") because the bridge’s label defaulter fills nu only for student(); that one-line fix sits outside this change. No phylogenetic, structured, or random-effect skew-normal routes are admitted; (1 | g) fails closed at DRM.jl. tweedie routes are admitted; (1 | g) fails closed at the same echo. ## engine = "julia" admits zero_one_beta() (fixed effects)

• drmTMB(bf(prop ~ x, sigma ~ z, zoi ~ w, coi ~ v), family = zero_one_beta(), engine = "julia") now routes to DRM.jl’s ZeroOneBeta instead of refusing with the Workflow G message. Both engines fit the same three-part mixture — mu (logit, the interior beta mean), sigma (log, phi = 1 / sigma^2), zoi (logit), coi (logit) — and all four dpars cross the design-258 label contract. Measured 2026-09-05 at DRM.jl pin 430ef64cc on the committed fixture (the tests/testthat/test-zero-one-beta.R draw, n = 1600): max |d coef| 3.98e-11, |d logLik| 1.93e-12, max |d fitted| 9.08e-12, per-coefficient Wald SE max relative difference 9.76e-07 over all 8 coefficients; the fit’s estimator ("ML") equals what DRM.jl reports as estim_method. Fixed effects only, and all four formula parts must be written: a (1 | g) term is refused by DRM.jl (ZeroOneBeta() currently supports fixed effects only), and a formula that omits zoi/coi (which native engine = "tmb" fits with intercept-only parts) aborts at the label echo — use engine = "tmb" for either. Not a phylogenetic, structured, or interval-coverage claim. ## Model comparison: aicc() and the nested likelihood-ratio test, ported from DRM.jl (#1117)

• New aicc() returns the small-sample corrected AIC, AIC + 2k(k + 1) / (n - k - 1), for a drmTMB fit on either engine (Inf when n - k - 1 <= 0, as in DRM.jl), and the internal drm_lrtest(reduced, full) computes DRM.jl’s lrtest() (statistic, df, p-value) with its REML, MAP and boundary-variance guards. Both are ports of DRM.jl src/comparison.jl at pin 430ef64cc, measured against DRM.jl’s own values on two committed fixtures: every quantity agrees to 1e-8 on the engine = "julia" object and to 2e-12 across engines. anova.drmTMB() still refuses; wiring it to drm_lrtest() is a follow-up. ## check_drm() now reads DRM.jl’s route-aware gradient/convergence diagnostics (#1108 / DRM.jl #569)

• DRM.jl #632 attaches the objective gradient at the optimum (gradient, index-aligned gradient_names) to the bridge payload, but only for routes whose internal fitter carries one; other routes correctly omit the field rather than filling zeros/NaN. check_drm() on a engine = "julia" fit previously errored with “no applicable method” even when that gradient was sitting unread in fit$bridge$gradient. check_drm.drmTMB_julia() now dispatches through the SAME generic a native TMB fit uses, reporting optimizer_convergence and a route-aware fixed_gradient row: numeric max|gradient| and its largest component when the route carries one, a NOTE naming the route (never a fabricated number) when it does not. Verified against DRM.jl 430ef64cc: the bivariate structured q2/q4 route and the sparse location-scale-scale ML route attach a gradient; the base Gaussian/GLMM route and the non-Gaussian phylo Laplace route do not — both directions are exercised by a live test, not just a mocked one.

Pinned DRM.jl clone moved from e0a65f96b to 430ef64cc

• The pin advances 18 commits, carrying DRM.jl #630 (sparse LSS gradient O(n+G) rather than O(G·n) — every warm-timing receipt taken at the old pin is now pessimistic, not wrong), #631 (profile endpoints never return an infinite bound; the sparse-LSS data race behind it is fixed), #632 (the bridge forwards optimizer options and exposes the stored gradient) and #633. Measured on the full live-Julia suite, all 37 test-julia*.R files in one session: 1440 passed, 0 failed, 0 errors, 0 skipped. Dependency versions were held fixed (Project.toml identical across the refs, so the previous pin’s Manifest.toml was carried over) and the clone was precompiled before measuring — the two controls whose absence invalidated an earlier probe.
• Compared file by file against the previous pin’s sweep, nothing changed — the one differing file is a test converted by the previous re-pin, not an engine change. In particular #631’s change to profile endpoint numbers reaches no test in this suite, which is recorded as a gap for the bridge-side inference qualification to close rather than as a pass.
• Totoro’s full-suite runner, which existed only on that host, is now vendored in DRM.jl as tools/totoro_run_suite.sh (DRM.jl #638).

engine = "julia" admits truncated_nbinom2() (fixed effects)

• drmTMB(bf(y ~ x, sigma ~ 1), family = truncated_nbinom2(), engine = "julia") now routes to DRM.jl’s TruncatedNegBinomial2 instead of refusing with the Workflow G message. Both engines fit the same zero-truncated NB2 target with the same size = 1 / sigma^2 parameterisation. Measured at DRM.jl pin 430ef64cc on the committed fixture (tests/testthat/test-family-dpq-batchC.R draw, n = 300): max |d coef| 8.81e-11, |d logLik| 2.84e-12, per-coefficient Wald SE max relative difference 2.71e-07; the fit’s estimator ("ML") equals what DRM.jl reports as estim_method. Fixed effects only: a (1 | g) term or a hu hurdle formula is refused (by DRM.jl, and by the design-258 label echo respectively), never silently dropped; use engine = "tmb" for those. Not a phylogenetic, structured, or interval coverage claim. ## engine = "julia": a structured marker with a non-intercept left side of the bar is refused before Julia boots

• phylo(1 + x | g), phylo(0 + x | g), and the same shapes on relmat(), animal(), and spatial() are now refused by an R-side capability gate the moment drmTMB(..., engine = "julia") is called, rather than reaching the JuliaCall subprocess and being refused there. At the pinned DRM.jl (430ef64cc, carrying DRM.jl #621’s _check_phylo_re_lhs), every such construct is already refused live — verified for a Gamma phylo(1 + x | species) slope, which DRM.jl throws as phylo(1 + x | species) is not implemented on the univariate routes -- only phylo(1 | species) (intercept) is after booting Julia. This gate moves the same refusal in front of that boot (defense-in-depth, drmTMB#1146), so it holds even if a future DRM.jl refactor narrows one of the more specific admission checks that also catch some of these shapes today. The gate is capability-gated, not a permanent ban: drm_julia_marker_slope_pin_supports() is the single switch to flip when a future pin fits one of these constructs. engine = "tmb" is unaffected and keeps fitting, for example, the two-SD Gaussian phylogenetic random slope.

engine = "julia": the ENGINE now decides which estimator ran, and two Poisson REML cells are admitted

• drmTMB labelled a Julia fit "REML" purely from its own support gate. Nothing checked what DRM.jl actually did, so a gate that was too permissive would have returned a maximum-likelihood fit labelled REML, with nothing to reveal it. DRM.jl #625 now reports estim_method on every fit, so the engine is the authority: if it says it fitted ML while REML was asked for and believed supported, the fit aborts and names it as a drmTMB gate defect rather than mislabelling the result. This closes the one failure direction a user cannot recover from.
• drm_julia_reml_supported() admits Poisson (1 | g) and Poisson phylo(1 | species), which it previously refused. Since the refusal landed these raised an avoidable error even though DRM.jl can restrict them. Widened on MEASUREMENT, not on report: each was verified with the estim_method oracle at pin e0a65f96b — poisson_random_intercept reports :REML with ml -123.1282 / reml -128.6623; poisson_phylo_intercept reports :REML with ml -74.6002 / reml -79.3865. Both restrictions are several log-likelihood units, not convergence noise.
• Not widened: the bivariate structured q=2 cell that DRM.jl’s refusal text also names. It was not verified, because drmTMB refuses that shape earlier and for an unrelated reason — the Julia route supports one phylo() term, and a q=2 spec carries one per mean axis. The REML question never arises there, so admitting it would have meant trusting the engine’s description of itself, which is exactly what the measurement above exists to avoid.

Pinned DRM.jl clone moved from 77513aa0 to e0a65f96b

• The engine = "julia" evidence and live tests are measured against a pinned DRM.jl checkout. That pin moved. Measured on the full live-Julia suite (all 37 test-julia*.R files, one session): 1433 passed, 0 failed, 1 error, where the single error is the one intended behaviour change below. The clone’s Project.toml is identical across the two refs, so the previous pin’s Manifest.toml was carried over and dependency versions are held fixed — only engine code differs between the measurements.
• phylo(1 + x | g) on a univariate route is now REFUSED by DRM.jl rather than silently fitted. At the old pin, DRM.jl’s parser kept only the grouping symbol, so it fitted an intercept-only model and quietly discarded x — the label was correct for the model it actually fitted, not for the two-SD random slope the user wrote. DRM.jl now fails closed. tests/testthat/test-julia-slope-nongaussian.R asserts that specific refusal, with negative controls so that a coef_labels regression, a precompile failure, or a missing-package error cannot pass as the refusal. A genuine two-SD phylogenetic random slope on non-Gaussian families remains unimplemented on both engines.
• The new pin also carries DRM.jl’s REML honesty work: the bridge now reports estim_method on every fit, ml_loglik always, and reml_loglik only when the fit really was restricted. A silent ML downgrade is therefore now directly observable rather than inferred from a likelihood gap.
• Reading that work surfaced a gap in the other direction: DRM.jl’s rewritten refusal message enumerates the cells it can fit by REML, and that list is wider than drm_julia_reml_supported() in three places — bivariate structured q=2, Poisson (1 | g), and Poisson phylo(1 | species). Since the Julia route now refuses rather than downgrading, those three raise an avoidable error. Tracked as #1152; deliberately not fixed here, because DRM.jl’s message is its own claim about itself and widening the gate on it, unverified, risks the more dangerous opposite fault.
• Provenance comments throughout the source and design notes that read “measured against DRM.jl 77513aa0” are left standing: they are dated records of what was verified then, and rewriting them would assert measurements never made at the new pin. docs/design/258 §7.8 records which of those statements a reader must carry forward differently.

Breaking: engine = "julia" now REFUSES unsupported REML = TRUE instead of quietly fitting ML

• Asking for REML = TRUE on a cell the DRM.jl bridge cannot fit by restricted maximum likelihood used to emit a warning and then fit by maximum likelihood anyway. The same script on two engines therefore returned two different estimators, separated only by a warning, and it did so on exactly the quantities REML exists to correct: variance components, and every ratio built from them (heritability(), repeatability(), icc()). Nothing in the returned fit announced the substitution beyond effective_REML, which nobody reads. It is now an error.
• The error names the cell, says the Julia engine cannot fit it by REML, and gives the two ways forward: REML = FALSE to ask for maximum likelihood explicitly, or a documented Gaussian REML cell. The engine = "tmb" pointer is deliberately bounded — TMB fits Gaussian cells by REML and has a separate diagnostic-only binomial route, but it does not offer a general REML fit for every cell the bridge refuses, and the error says so rather than sending users to a route that will also decline.
• The refusal fires before anything is marshalled to Julia, so an unsupported request costs no fit.
• Why this is not gated on finishing the support census: the two failure directions are not symmetric. A gate that is too narrow produces a visible error next to a route that works — recoverable in one line. A gate that is too wide hands back ML labelled REML — not recoverable, because nothing tells you to look. The refusal only ever affects the first direction; widening the gate where the census shows DRM.jl is more capable than drmTMB assumes is a separate, later change. Census batch 1 (4 cells) agrees with the current gate on every cell and found no instance of the dangerous wide direction.

First CRAN-targeted release candidate; not yet submitted to or accepted by CRAN. drmTMB fits distributional regression models – location, scale, shape, zero inflation, and residual correlation – for one or two responses, using Template Model Builder.

Regression cover for the non-interactive engine = "julia" abort

• engine = "julia" used to abort in an ordinary Rscript session with “Live Julia setup is disabled on the CRAN check lane”, because the CRAN-lane predicate inferred “package check” from !interactive(). A scripted analysis is a supported way to use the package, so this blocked real non-interactive use (reported by an outside user working on a 10,970-species tree). The behaviour was fixed on 2026-08-30 by requiring R’s own _R_CHECK_PACKAGE_NAME_ marker; this entry adds the regression test the fix shipped without. Reverting the fix left the existing CRAN-lane test file fully green (32 passed, 0 failed), so the bug could have been reintroduced without CI noticing; the new tests/testthat/test-julia-noninteractive-lane.R fails against the pre-fix predicate.
• The marker choice is now backed by measurement rather than inference (docs/dev-log/evidence/julia-r-parity/check-lane-markers/): probe packages under R CMD check show _R_CHECK_PACKAGE_NAME_ set in the examples lane, in tests/testthat.R, and inside test_check(), but absent from the vignette rebuild subprocess, which carries no check marker at all. TESTTHAT_IS_CHECKING would have been the wrong marker: it is missing from the examples lane, the lane the guard was originally written for. Because the vignette lane is invisible to any marker, the new test file also scans every vignette and fails if an evaluated chunk ever calls engine = "julia", with a positive control proving the scanner detects a planted chunk.

Renamed: summary()’s derived “repeatability”/“phylogenetic_signal” rows

• summary()’s derived component printed rows labelled "repeatability" and "phylogenetic_signal" for structured mu random-effect components. These rows divide by the TOTAL variance (every mu random-effect variance in the fit, summed, plus the residual variance). The icc() / repeatability() accessors, added separately, use a different focal-vs-residual denominator (that one component’s variance over itself plus the residual only) under the same word “repeatability” – the two agreed on a fit with one structured mu component and silently disagreed on a fit with two or more. summary()’s rows are renamed to "total_variance_share" and "phylo_total_variance_share", which state their denominator; their arithmetic is unchanged. icc(), repeatability(), and heritability() are unaffected – they keep their existing names and values. If you matched on the row name "repeatability" or "phylogenetic_signal" in summary(fit)$derived, or on the profile-target parm strings "derived:repeatability(<group>)" / "derived:phylogenetic_signal(<group>)", update to the new names. See docs/design/259-heritability-icc-repeatability.md section 3 item 5 for which function owns which denominator.

Fixed: profile_targets() on an engine = "julia" fit listed one target where confint() accepts five (#1156)

• profile_targets() returned early for a Julia fit with only the phylogenetic SD inventory, so bf(y ~ x + phylo(1 | species), sigma ~ 1) listed 1 target (sd:mu:phylo(1 | species)) while the same TMB fit listed 6 – even though confint(fit, parm = "fixef:mu:x", method = "profile") on the Julia fit already worked when the name was typed in full. Discovery now returns the union the engine accepts, in the native row order: fixef:mu:(Intercept), fixef:mu:x, fixef:sigma:(Intercept), sigma, sd:mu:phylo(1 | species) (measured live against the pinned DRM.jl clone 430ef64cc; each of the four profile_ready rows returns a finite profile interval). The sigma response-scale alias is listed as it is on the Wald path – profile_ready = FALSE, note missing_tmb_parameter – and the native fit’s derived phylo_total_variance_share row has no bridge counterpart, so it remains TMB-only.

Fixed: nlminb could report convergence short of the true optimum on flat surfaces (#1130)

• nlminb()’s own rel.tol/x.tol stopping rules trigger on a relative change in the objective or step, not on the gradient norm, so on a flat or weakly-curved surface (for example a variance component identified by few groups) it could report convergence == 0 while the exact TMB gradient at that point was still far from zero – invisible at the training rows but amplified at extrapolated newdata. drmTMB() now applies a Newton polish by default after every nlminb() fit: a few Newton steps on the exact TMB gradient and a finite-differenced Hessian of it, accepting a step only if it does not increase the objective, until the gradient’s max absolute component is small. This does not change nlminb’s own tolerances (tightening them was tried and found unreliable: see docs/design/260-nlminb-newton-polish-optimizer-tolerance.md), costs a handful of extra gradient evaluations at the reported optimum (skipped entirely once the gradient is already small; measured negligible on a small Gaussian, a location-scale, and a phylogenetic fit), and honestly recomputes the reported convergence code and message from the polished gradient rather than trusting nlminb’s own diagnostic. Opt out with drm_control(newton_polish = FALSE).
• confint(..., method = "profile", profile_engine = "endpoint") compares the free fit’s objective against a constrained endpoint solve; the polish above originally ran on the free fit only, so a constrained endpoint could stop short of its own minimum relative to the now-lower free objective and get rejected by the existing gradient guard, narrowing or invalidating the interval. The same polish now runs on the constrained endpoint solve too, keeping the comparison symmetric (and skipped on both sides together under drm_control(newton_polish = FALSE)).
• The constrained ordinal cutpoint profile (confint(fit, parm = "ordinal:cutpoint:<label>", method = "profile")) had the same asymmetry (#1144): its endpoint solve stopped wherever nlminb alone stopped. On the committed random-intercept cumulative_logit() fixture that was a gradient of 3.6e-4 to 1.9e-3 at the six endpoints (two above the 1e-3 acceptance guard) against 8.5e-10 at the free optimum. Both endpoints of every cutpoint are now polished the same way (measured 4e-14 to 7e-9 afterwards); the gradient guard PROFILE_ENDPOINT_GRADIENT_TOL is unchanged at 1e-3, and on that fixture the reported intervals move by less than 1e-9 – the fix is to the solve’s consistency, not to a visibly wrong interval.

Fixed: bootstrap confint() failed for cbind(successes, failures) binomial fits (#1123)

• confint(fit, method = "bootstrap") errored with “Bootstrap confidence intervals require a stored response column in the fitted data” for binomial fits whose response used trial-denominator syntax – bf(cbind(successes, failures) ~ x) – even though method = "wald" and method = "profile" worked on the same fit. The internal resampler only knew how to write a single simulated response column back into the refit data, so a two-column cbind() response tripped its guard. Bootstrap replicates now rebuild the two-column response from the simulated successes and each row’s original trial size; Bernoulli 0/1 and proportion-with-weights binomial fits are unaffected.

Fitted objects now store the final gradient (DRM.jl #569, R-side)

• drmTMB() fits now carry $gradient (the final outer gradient at opt$par) and $gradient_max_component (the label of its largest |component|), computed once at fit time. Previously this quantity was only available lazily, inside check_drm()’s fixed_gradient row, and only when the fit retained its TMB object (drm_control(keep_tmb_object = TRUE), the default). The stored value agrees exactly with what check_drm() computes live, and survives keep_tmb_object = FALSE.

Public start contract: drm_control(start = list(...))

• drm_control() gains a start argument implementing the public start contract from design 35 (“Public Start Contract”): a named list keyed by "fixef:<dpar>:<column>", "sd:<dpar>:<term>", or "cor:<dpar>:<term>" labels. Labels are validated against the parsed formula and family before optimization (unknown labels error before the fit runs, never during it). Values are transformed to the internal unconstrained scale before TMB::MakeADFun() sees them: sd: starts are natural-scale and log()-transformed, cor: starts are natural (-1, 1)-scale and atanh()-transformed, and fixef: starts are already on the model’s internal link scale. A partial start updates only the named targets; every other component keeps the ordinary family-builder default. Latent random-effect (u) values are not addressable through this contract. This is a validated translation layer onto the existing private drm_apply_start_override() hook – nothing changes at the TMB boundary, and a start changes only where the optimizer begins, never what is reported. start_from = <a fitted model> remains reserved and unimplemented.
• Two defects found by an adversarial pass are fixed. (1) A fixef:sigma: (or sigma1:/sigma2:) start whose implied log(sigma) linear predictor falls outside the configured logsigma_clamp band now warns (drmTMB_start_clamp_saturated_warning) before the fit that produces the bad result: the softclamp derivative is ~0 in that region, so nlminb can see an already-flat gradient at the start and report a spuriously clean convergence without ever moving. The start itself is never silently moved.
  2. A fixef:mu: (or scale-side fixef:sigma: when a sigma variance component makes REML fold it in too) start under REML = TRUE now errors instead of being silently accepted with no effect: REML integrates those fixed effects into the Laplace random block, so there is no free coordinate for the start to seed. This matches objective_at()’s existing refusal of the same labels under REML.

objective_at(): evaluate the fitted objective at a supplied point

• New exported S3 generic objective_at(object, at = list(...)), with a drmTMB method, implementing design 35’s “Objective At A Point”: it shares the start= label vocabulary ("fixef:<dpar>:<column>", "sd:<dpar>:<term>", "cor:<dpar>:<term>") to evaluate a fitted model’s TMB objective (negative log-likelihood) at a supplied point, without refitting. It reuses the existing R/profile.R evaluation pattern (substituting into a copy of fit$opt$par and calling fit$obj$fn()) rather than a second implementation. This is a diagnostic: it selects nothing, reports no uncertainty, and does not mutate the fitted object (the TMB object is re-pinned to its optimum after evaluation). Requires drm_control(keep_tmb_object = TRUE) (the default). Unknown labels error before evaluation, on the same rule as start=. For penalized (MAP) fits the return value is on the same unpenalized convention as logLik() (any penalty is subtracted back out, so objective_at(fit, <own optimum>) == -logLik(fit) holds for every fit type), and it errors for experimental MSPL fits exactly as logLik() does.
• start=/objective_at() labels now reach biv_gaussian’s rho12 fixed effect (via the existing generic fixef:rho12:<column> label – beta_rho12 is now named at start-build time) and the phylogenetic location/scale covariance block via three new families, phylo_sd:<axis>, phylo_cor:<axis1>:<axis2>, and phylo_theta:<axis1>:<axis2> (design 35, “Phylo Covariance Block”). phylo_cor: always denotes a correlation strictly in (-1, 1) (q = 2 or block-diagonal q > 2); on a DENSE q > 2 block, where the working parameter is a raw UNSTRUCTURED_CORR_t Cholesky-space entry with no correlation meaning in isolation, phylo_cor: refuses and names phylo_theta: – the honestly-named, unconstrained-scale counterpart – as the alternative, rather than returning a value on a scale its own name does not promise. Closes a gap the A5 cross-engine receipt had to work around by addressing beta_rho12/log_sd_phylo/theta_phylo by internal TMB parameter name.

Build provenance (drm_provenance())

• New exported drm_provenance() (DRM.jl#473) answers a question packageVersion("drmTMB") cannot: which BUILD produced this attachment, not just which release. Two builds can report the same version string while differing by commits under R/, src/, or NAMESPACE. drm_provenance() returns package_version, git_sha, git_dirty, build_time, source, reason, and queried_at. Every fitted drmTMB object now carries this record at fit$provenance.
• The git SHA and working-tree dirty flag are captured once, at build time, by ./configure/configure.win (plain git rev-parse/git status --porcelain, non-fatal) and baked into inst/build-provenance.dcf. An installed package has no .git to consult, so drm_provenance() never attempts a live git lookup at call time; when the baked file is absent (e.g. devtools::load_all(), or a build whose configure step did not run) it returns source = "unavailable" with a stated reason, and source = "baked-without-git" when configure ran but git itself was unavailable (no git binary, or no .git, as when installing from a released tarball).
• tools/drmtmb_provenance.R is a separate, unshipped developer/CI tool (excluded from the built tarball via the pre-existing .Rbuildignore tools/ rule, same as every other script in that directory) for stamping a build anchor into Julia-side fixture receipts. It duplicates rather than shares the few git-capture lines with configure, deliberately: a shipped configure must not depend on an unshipped file.
• This is drmTMB-side only. DRM.jl has no equivalent provenance surface yet.

Student-t response + one binary mi() predictor

• A student() response can now carry one binary mi() predictor (y ~ z + mi(treatment) with impute_model(treatment ~ z, family = binomial())). This is drmTMB #962 / S6 A7: C++ has_mi wiring via drm_student_log_density (identity location; nu = 2 + exp(eta_nu)). The shared 7-arg leaf ABI is not extended — see LOOP/notes/A7-student-nu-abi.md. Ledger cell mp-student-bernoulli records MCAR + MAR point recovery. This is not FIML, not impute_joint, not k ≥ 2, and not a continuous missing predictor under student. zi-* responses remain gated. Capability stays partial.

nbinom2 response + one Gaussian mi() predictor

• An nbinom2() response can now carry one Gaussian mi() predictor (y ~ z + mi(x) with impute_model(x ~ z, family = gaussian())). This is drmTMB #962 / S6 A7: C++ has_mi for mi_family == 0 inside model_type == 7, not a whitelist-only edit (nbinom2 was already on drm_missing_predictor_families()). Ledger cell mp-nbinom2-gaussian records a joint-density identity plus MCAR + MAR point recovery. This is not FIML, not impute_joint, not k ≥ 2, and not a grouped or structured predictor model. Poisson / binomial / beta responses still take a binary predictor only. zi-* responses remain gated.

Zero-inflated Poisson + one binary mi() predictor

• A zero-inflated Poisson response (family = poisson() plus zi ~ 1) can now carry one binary mi() predictor in mu only (y ~ z + mi(treatment), zi ~ 1 with impute_model(treatment ~ z, family = binomial())). This is drmTMB #962 / S6 A7 / D-23: C++ has_mi wiring that inlines the ZIP mixture in model_type == 8. It does not reuse the plain Poisson leaf for structural zeros, and eta_zi comes from observed-only predictors. Ledger cell mp-zi-poisson-bernoulli records MCAR + MAR point recovery. This is not FIML, not impute_joint, not k ≥ 2, not mi() on zi, and not zi_nbinom2. Student remains gated. Capability stays partial.

Beta-binomial response + one binary mi() predictor

• A beta_binomial() response can now carry one binary mi() predictor (cbind(success, failure) ~ z + mi(treatment) with impute_model(treatment ~ z, family = binomial())). This is drmTMB #962 / S6 A7: C++ has_mi wiring and a drm_response_log_density beta-binomial leaf (logit success probability), not a whitelist-only edit. Ledger cell mp-beta-binomial-bernoulli records MCAR + MAR point recovery. This is not FIML, not impute_joint, not k ≥ 2, and not a continuous missing predictor under beta_binomial. Student and zi-* responses remain gated.

Lognormal response + one binary mi() predictor

• A lognormal() response can now carry one binary mi() predictor (y ~ z + mi(treatment) with impute_model(treatment ~ z, family = binomial())). This is drmTMB #962 / S6 A7: C++ has_mi wiring and a drm_response_log_density lognormal leaf (identity log-location), not a whitelist-only edit. Ledger cell mp-lognormal-bernoulli records MCAR + MAR point recovery. This is not FIML, not impute_joint, not k ≥ 2, and not a continuous missing predictor under lognormal. Student and zi-* responses remain gated.

Gamma response + one binary mi() predictor

• A Gamma(link = "log") response can now carry one binary mi() predictor (y ~ z + mi(treatment) with impute_model(treatment ~ z, family = binomial())). This is drmTMB #962 / S6 A7: C++ has_mi wiring and a drm_response_log_density Gamma leaf, not a whitelist-only edit. Ledger cell mp-gamma-bernoulli records MCAR + MAR point recovery. This is not FIML, not impute_joint, not k ≥ 2, and not a continuous missing predictor under Gamma. Lognormal, student, beta_binomial, and zi-* responses remain gated.

Two independent Gaussian mi() terms

• A Gaussian response can now carry two independent mi() terms with a named impute list (y ~ mi(m1) + mi(m2) + x). This is #963 option (b): product of margins, no residual correlation. It is not impute_joint and not FIML. Non-Gaussian responses still require exactly one mi() term. Ledger cell mp-gaussian-gaussian-k2-indep records MCAR + MAR point recovery.

Lognormal ordinary correlated slope

• Complete-data ordinary lognormal() y ~ x + (1 + x | id) now has a point_fit_recovery ledger cell (mc-0720) for sd0, sd1, and group-level rho_re under the same design-17 map as Poisson (ρ = 0.999999 tanh(η)). This is not the independent-slope cell mc-0380 and not Wave 2.5 NB2 mc-0719. A slope predictor that is constant within every group now aborts before the fit (unidentified sd1 / rho_re). Gamma neighbour, REML, missing-response, labelled blocks, mixed (1 | g) + (1 + x | g), intervals, and coverage remain rejected.

• CRAN / win-builder Julia hang (post-#1061). The ^julia invert filter correctly excluded test-julia-*.R, but test-binomial-response.R still called drmTMB(..., engine = "julia") inside expect_error(). Workflow G admits fixed-effect binomial into the Julia bridge, so that call reached JuliaCall::julia_setup() and hung Ligges R-release for ~10448s. drm_julia_setup() now hard-blocks the non-interactive CRAN lane unless DRMTMB_JULIA_TESTS=true, drm_skip_live_julia() matches the same predicate, and the obsolete binomial expect_error(engine = "julia") is replaced with a pure-R tag check.

Binomial ordinary correlated slope

• Complete-data binomial logit cbind(success, failure) ~ x + (1 + x | id) now has a point_fit_recovery ledger cell (mc-0717) for sd0, sd1, and group-level rho_re under the log-sech Cholesky. This is not the independent-slope cell mc-0061. A slope predictor that is constant within every group now aborts before the fit (unidentified sd1 / rho_re). REML, missing-response, labelled blocks, mixed (1 | g) + (1 + x | g), intervals, and coverage remain rejected.

Poisson ordinary correlated slope

• Complete-data Poisson log count ~ x + (1 + x | id) now has a point_fit_recovery ledger cell (mc-0718) for sd0, sd1, and group-level rho_re under the design-17 map ρ = 0.999999 tanh(η). This is not the independent-slope cell mc-0431 and not Wave 1 binomial mc-0717. REML, missing-response, labelled blocks, mixed (1 | g) + (1 + x | g), intervals, and coverage remain rejected.

NB2 ordinary correlated slope

• Complete-data ordinary nbinom2() count ~ x + (1 + x | id) now has a point_fit_recovery ledger cell (mc-0719) for sd0, sd1, and group-level rho_re under the same design-17 map as Poisson (ρ = 0.999999 tanh(η)). This is not the independent-slope cell mc-0402 and not Wave 2 Poisson mc-0718. A slope predictor that is constant within every group now aborts before the fit (unidentified sd1 / rho_re). Zero-inflated or truncated NB2, REML, missing-response, labelled blocks, mixed (1 | g) + (1 + x | g), intervals, and coverage remain rejected.

Binomial responses accept a phylogenetic random effect

• binomial() was the only common response family that could not take a structured random effect – gaussian, poisson, nbinom2, Gamma and beta all accepted the identical phylo(1 | id, tree = tree) term while a binomial model aborted at the structured-effect gate (#1048). Phylogenetic logistic regression (a binary trait on a tree) is the canonical comparative-methods use of a binary response, so the hole was conspicuous.

  The first slice is deliberately narrow, matching how beta and zero_one_beta grew provider by provider: one unlabelled q1 phylo() intercept on mu, in either the Bernoulli or the two-column cbind(successes, failures) form. Phylogenetic slopes, labelled covariance blocks, spatial/animal/relmat providers, combination with ordinary random effects, and combination with missing-predictor mi() all refuse with explicit messages rather than fitting silently narrower models. Supporting recovery evidence (slope and sd_phylo essentially unbiased by 160 tips, Laplace small-sample attenuation ~9% at 40-80 tips) is recorded on #1048.

Fixed: penalized fits reported phylo_penalty and logLik off the optimum

• A penalized (MAP) phylogenetic fit read its penalty from a bare obj$report(), which TMB evaluates at obj$env$last.par. After TMB::sdreport() that is a finite-difference step away from the optimum, so fit$phylo_penalty – and therefore fit$logLik, which is -opt$objective + phylo_penalty – were both slightly wrong. The error scaled with the penalty (order 1e-3 at sd_u = 0.5) and, being a diagnostic-time artefact, made both values depend on whether standard errors had been requested at all: the same fit at the same optimum reported different numbers under se = TRUE and se = FALSE.

  The practical consequence was in drm_phylo_penalty_sweep(), whose logLik column is the basis for deciding whether a coupling is data-informed or prior-shaped. The error differed per row, so it did not cancel in exactly the comparison the sweep exists to make.

  Estimates were never affected – opt$par and sdreport() were always at the optimum. check_drm()’s log(sigma) clamp row had the same defect and is fixed alongside. The existing regression test could not catch this because it derived its “expected” penalty from the same bare report() call, so both sides of the assertion moved together; it now derives the expectation from parList(opt$par), and a new test asserts that se cannot move a reported estimate.

Experimental MSPL accepts probit and complementary log-log

• estimator = "mspl" previously required binomial(link = "logit") exactly. It now also accepts binomial(link = "probit") and binomial(link = "cloglog"). Everything else about the route is unchanged: still one complete Bernoulli or grouped-binomial model, one ordinary q = 1 or correlated q = 2 grouping block, engine = "tmb", no REML, no intervals.
• The basis is measured, not assumed. Kosmidis and Firth (2021, Theorem 1 with Section 3.1) prove the Jeffreys penalty yields finite estimates for any link whose working weight vanishes in both tails, which logit, probit and cloglog all do. drmTMB’s own TMB-Laplace evidence: no non-finite estimate in 43,972 completed fits, and Wald standard errors calibrated in the identified regime (probit mean(SE)/sd(beta) in [0.946, 1.008], cloglog in [0.957, 1.027]). Artifacts under docs/dev-log/simulation-artifacts/2026-08-11-mspl-nonlogit-links/. Scope of that measurement: Bernoulli responses with two fixed-effect columns. Grouped-binomial and wider designs inherit the n_eff extrapolation noted below without direct measurement – which was already true of the logit route and is not made worse here, only left unresolved for two more links.
• The soft-penalty scale is unchanged and is a logit constant. c_n = 2 * sqrt(p / n_eff) comes from a delta-method argument at beta = 0 for the logit link; the same argument gives about 1.25 * sqrt(p / n_eff) for probit and 1.31 for cloglog. Using the logit constant for all three was measured to move the estimate by roughly 1% of one standard error (at q = 1 with two fixed-effect columns; not measured at q = 2), so it is kept rather than made link-specific – a per-link constant would define a different estimator.
• Know this before using it on rare events. Under deep separation with a random slope, the standard error is frequently unavailable: vcov() and summary() return NA with a drmTMB_mspl_wald_unavailable warning rather than a fabricated number, and in the most extreme cells measured this affected the large majority of converged fits. The point estimate remains finite – that is what the penalty is for – but inference is not available there for any link.
• MSPL start values now place the intercept at the link of the observed event rate instead of at zero. At beta = 0 cloglog implies an event rate of 1 - exp(-1) = 0.632, so a rare-event design previously began several log units from its own intercept; this removes a class of avoidable optimizer failures. Slopes still start at zero and the start remains deterministic and finite.

Boundary intervals are now flagged on every route

• confint(method = "bootstrap") warns at a variance-component or correlation boundary (class drmTMB_bootstrap_boundary_warning, conf.status = "bootstrap_at_boundary"). The Wald and profile routes already flagged their own boundary cases; bootstrap did not, which made it the one route that could return a clean-looking interval for a target the other two warn about. Resampling does not repair a boundary — a percentile interval whose draws pile up at zero is reporting the constraint, not the sampling distribution.
• The flag fires when at least 5% of retained resamples land on the bound, and requires at least 20 retained draws, because a share computed from a handful of resamples is noise rather than evidence. Calibration at R = 200: a true SD of 0 put 43% of draws on the bound; a true SD of 0.25 put 5% there, with a lower endpoint of exactly zero; a true SD of 0.9 put none there.
• check_drm() now states what it does not check. It reads the fit, so a random-effect SD comfortably clear of zero passes every fit-level check while confint() still warns about that same target’s interval — measured on a 10-group fit with an SD estimate of 0.1936, every check ok, and 43% of bootstrap draws on the bound. A new interval_reliability_scope note records that interval reliability is not assessed there and points at conf.status. Previously that all-clear could reasonably be read as permission to report.

Offsets in mu for every univariate family

• offset() in the mu formula now works for all univariate families, not only Poisson, nbinom2(), and Bernoulli/binomial. Newly accepted: Gaussian, student(), skew_normal(), lognormal(), Gamma(), tweedie(), beta(), beta_binomial(), zero_one_beta(), and cumulative_logit().

• Read the offset on the family’s link scale. For log-link families it is the usual offset(log(exposure)) rate model. For identity-link families it is a known additive shift of the mean. For logit-link families it is a known log-odds shift — a calibration term, not an exposure; offset(log(effort)) on a logit-link family does not give a rate model.

• A zero_one_beta() offset shifts only the interior beta component and leaves zoi/coi untouched. A cumulative_logit() offset shifts the latent location against fixed cutpoints rather than an intercept.

• Still rejected, deliberately: truncated_nbinom2() and its hurdle path (both renormalise the observed mean over a restricted support, so an exposure term would not scale the reported mean — extending them needs a decision on whether the offset targets the latent untruncated rate or the observed mean), every bivariate family (no per-response offset contract yet), Gaussian sufficient-statistic aggregation, and every distributional parameter other than mu.

• Correctness is pinned by a link-agnostic identity: a constant offset c must lower the fitted intercept by exactly c and leave every slope unchanged. tests/testthat/test-offset-families.R checks that for each family.

• Fixes a documentation defect: the drmTMB() help previously claimed offset(log(exposure)) support for zero-truncated negative-binomial mu formulas, which the code rejected. ## Live Workflow G engine = "julia" FE gate (#499)

• Skip-safe live tests now round-trip DRM.jl’s eleven admitted Workflow G fixtures through drmTMB(..., engine = "julia") against committed expected.toml numbers (drmTMB 0.6.0 pin). Measured locally: 11/11 pass when JuliaCall + DRM.jl are available.

• The R bridge now admits those Workflow G fixed-effect families (Gaussian, bivariate Gaussian, Student-t, lognormal, Poisson, NB2, Gamma, Beta, Binomial) without requiring a phylo() term; large-p phylo routes for the count / non-Gaussian set remain. Beta-binomial stays gated.

• Marshalling fixes that the live gate surfaced: expand cbind(successes, failures) into data columns, and rewrite meta_V(V = v) / drmTMB::meta_V(...) to positional meta_V(v) for DRM.jl’s StatsModels parser.

• Claim fence: experimental live parity for the eleven cells only — not a CRAN-default Julia engine, not closure of all of #499 (FIML / remaining honest scope stay open).

Binomial models accept probit and cloglog links (new in 0.7.0)

• binomial(link = "probit") and binomial(link = "cloglog") now fit, alongside the existing binomial(link = "logit"). predict(), summary(), and predict_parameters() back-transform through the fitted link rather than assuming logit. Other binomial links (for example cauchit) are still rejected.
• The link is evaluated on the log scale throughout (a tail-safe log-normal CDF for probit; a log1mexp form for cloglog), not by computing a probability and clamping it. This keeps accuracy in the extreme tails; see docs/design/252-binomial-link-generalisation.md §3 and inst/COPYRIGHTS.
• This is a fitting capability, not an interval or coverage claim. The new links inherit binomial’s existing capability-ledger evidence; no new recovery or coverage campaign was run for them, and the census is unchanged.
• The DRM.jl bridge continues to reject probit and cloglog: DRM.jl implements the logit mean only, so engine = "julia" errors rather than silently fitting a different model.
• Superseded later in 0.7.0: the maximum-softly-penalized-likelihood (MSPL) entry point was initially logit-only. It now also admits probit and complementary log-log links under the explicitly bounded, experimentally supported contract described in “Experimental MSPL accepts probit and complementary log-log” above.

Default uncertainty story

• Fixed effects / routine Wald targets: start with confint(fit) (method = "wald"). It is the fastest fitted-object route when TMB::sdreport() succeeded.
• Random-effect SDs / variance components: prefer confint(fit, parm = ..., method = "profile") after profile_targets(fit) lists the row as profile-ready.
• Boundary rows: read profile.boundary and conf.status. A usable profile interval at a variance-component boundary warns with class drmTMB_profile_boundary_warning. That warning does not repair coverage; treat the interval as indicative of scale, not as a calibrated level interval (see ?confint.drmTMB Boundary intervals; D-117 evidence).
• No nominal-coverage-everywhere claim. A computable interval is not coverage certification. Report a route as inference-ready only when the capability guide names that exact cell.

See also the vignette First-week intervals: fit, profile, and boundary.

REML capability wording now matches the callable surface

• The capability guide now records the bounded O2 binomial REML route: one ordinary unlabelled location (mu) random intercept (mc-0060) or independent slope (mc-0062). This route is diagnostic-only despite deterministic agreement with glmmTMB(REML = TRUE); it has no recovery or coverage claim. Fixed-only, multiple-term, correlated, labelled, structured, and missing-response binomial REML remain unavailable.
• README and ?drmTMB now state the same public non-Gaussian REML boundary without requiring the capability guide: binomial O2 only, both diagnostic-only; O3 remains package-private; mc-0227 stays public ML point_fit_recovery. Other non-Gaussian families reject REML = TRUE.
• The retained cumulative-logit AGHQ plus Cox-Reid campaign is now identified as evidence for a package-private estimator. It does not grant a public drmTMB() fit or reader-facing reporting permission.
• The guide no longer describes zero-one-beta structured effects as wholly absent. It names only the exact q1 structured-intercept gates and their earned point-recovery or interval-feasible tiers.

---

drmTMB 0.6.0

Five Prong B routes now have interval-feasible profile evidence

• A Totoro 135-trace campaign (source SHA 6618e4b30, GNU parallel, ≤100 cores) retained five-seed stats::profile() / TMB::tmbprofile receipts for fourteen Prong B Tier-1 candidates. Five cells cleared the preregistered ten-clause contract on every seed and are promoted point_fit_recovery → interval_feasible: mc-0568 and mc-0576 (ordinary zero_one_beta() sigma intercept and slope), mc-0595 and mc-0596 (zero_one_beta() sigma under relmat() and spatial()), and mc-0653 (zi_nbinom2() sigma under phylo_interaction(), 8×8 campaign DGP). Census move: model_surface interval_feasible 182→187; frozen point_fit_recovery 59→54.

• Nine siblings stay at point_fit_recovery after honest WITHHOLD (mc-0593, mc-0594, mc-0597, the five labelled count-mu q2 cells, and mc-0425 at 4/5). Per the preregistration, 4/5 truth-bracketing is a block, not an 80% pass. This is still not coverage or calibration.

• Structured-sigma promotions (mc-0595, mc-0596, mc-0653) name in claim_boundary the documented ML sigma-axis low bias and that REML is unavailable for these families. Evidence: docs/dev-log/simulation-artifacts/2026-08-05-135-trace-campaign/.

Profile intervals now warn at a variance-component boundary

• confint(fit, method = "profile") now warns when it returns a usable interval whose profile.boundary flag is TRUE, with condition class drmTMB_profile_boundary_warning. Until now only the Wald path warned at a boundary (drmTMB_wald_boundary_warning), and it steers the user to method = "profile" – into a regime the package had measured as worse, with no signal. Conditional on the boundary flag, the D-117 10-group random-effect SD gate measured profile coverage at 0.1021, 0.2387, and 0.8683 against a nominal 0.95 – and 0 of 89 in a fourth cell where boundary hits are rare (0.09% of fits) – driven by a random-effect SD point estimate biased 8.3%-15.8% low, which anchors the interval low and makes it miss from above. Unconditionally, over all fits, the same gate measured 0.9248 against a nominal 0.95: a flagged interval is the bad case, not the typical one, and the flag tells you which you have. These figures are from a 400,000-attempt campaign and supersede the 1,000-replicate figures quoted in earlier releases.

• This is not a drmTMB defect, and the warning does not report one. lme4::lmer on the same data-generating process and the same seeds agreed on boundary incidence for 4000/4000 replicates and matched the conditional coverage to four decimal places; in the single divergence, lme4 returned an interval that excluded its own maximum likelihood estimate. Sub-nominal coverage here is a property of profile intervals near a variance boundary, not of this implementation. The comparator ran at 1,000 replicates per cell and was not re-run at 400,000. Evidence: docs/dev-log/simulation-artifacts/2026-08-04-d117-10group-profile-gate/ and .../2026-08-09-d117-100k-regate/.

• Scope. All of the above was measured on one design: Gaussian, a single random intercept on the mean, 10 groups, 4 or 10 observations per group, maximum likelihood. It is not a general claim about other families, providers, or group counts, and 0.9248 against a nominal 0.95 is real undercoverage – the package’s small-sample floor is deliberately tapered with the group count for that reason.

• REML = TRUE improves this design without repairing it – now measured. Over 400,000 paired replicates on the same design and seeds, refitting with REML = TRUE moved profile coverage from 0.9248 to 0.9463 against a nominal 0.95, roughly halved the SD point estimate’s downward bias (pooled -10.9% under maximum likelihood to -4.6% under REML), and brought the upper-to-lower miss asymmetry from 5.7:1 to 2.0:1 – a better-centred interval, not merely a wider one. Coverage conditional on the boundary flag improved but stayed well below nominal (0.74 to 0.83), so the boundary warning applies under either estimator. The maximum likelihood control arm of this campaign reproduced the 400,000-attempt gate above to five decimal places. The default estimator is unchanged, and this measurement covers the same single design as the figures above. Evidence: docs/dev-log/simulation-artifacts/2026-08-15-d117-reml-arm/.

• The interval is still returned – a boundary is a warning, not an auto-discard, matching the Wald path. Only usable intervals are flagged: profile_failed and clamp_limited rows also carry profile.boundary = TRUE but return missing endpoints, and they already report themselves through conf.status, so warning about the coverage of an interval that was never returned would be noise.

Direct profile targets for fourteen count and zero-one-beta routes

• confint(fit, method = "profile") and profile(fit) now reach the random-effect SD targets of fourteen routes that previously returned profile_ready = FALSE with a point_fit_only_* note: the zero_one_beta() ordinary sigma intercept and slope routes; the same family’s sigma routes under phylo(), animal(), relmat(), spatial(), and phylo_interaction(); the labelled intercept-slope covariance blocks of poisson() and nbinom2() mu, including their correlation target; and the nbinom2() and zi_nbinom2() sigma routes under phylo_interaction(). The retired notes are point_fit_only_count_q2, point_fit_only_count_sigma_interaction, point_fit_only_zi_nbinom2_sigma_interaction, and point_fit_only_zero_one_beta_sigma_q1.

• This change makes a profile computable on these routes. It is not a claim that the resulting interval attains nominal coverage. Five of the fourteen routes later cleared a seeded Totoro campaign and are now interval_feasible (see the Prong B section above); the other nine remain point_fit_recovery. Read an interval from an unpromoted route as “the profile traversed and returned two finite ordered endpoints”, not as calibrated inference.

• Seven of the fourteen profile a sigma-axis random-effect SD under a structure: zero_one_beta() with phylo(), animal(), relmat(), spatial(), or phylo_interaction(), and nbinom2() or zi_nbinom2() with phylo_interaction(). For this class of cell the underlying maximum likelihood point estimate is documented as biased low — eleven of twelve retained ML estimates fell below truth in the sibling nbinom2 provider cells, with a fit-level one-sided sign test at p = 0.0032 and a cell-level figure of p = 0.0625 once the shared data-generating process and shared random-number stream are respected. Neither existing correction reaches these routes: bias_correct shifts only the method = "wald" centre, never a profile endpoint, and native scale-side REML is unavailable because drm_validate_reml_spec() admits only Gaussian and binomial models. A profile interval built around a low-biased point estimate inherits that bias uncorrected; interpret these seven routes accordingly.

• Reaching a profile target is not the same as obtaining an interval from it. The change makes these targets eligible for stats::profile(); whether a given dataset yields two finite endpoints still depends on that fit. The zi_nbinom2() phylo_interaction() sigma route is a known case in point: on the fixture used for its point-fit evidence the variance component collapses to the lower boundary (estimate 5e-05 against a generating value of 0.60), and its profile returns conf.status = "profile_failed" at some ystep settings and a lower endpoint of zero with a near_sd_boundary message at others. That degeneracy predates this change and is a property of the fixture, not of the profile machinery, but it is the reason this entry claims reachability rather than intervals.

• zi_nbinom2() ordinary sigma q1, zero_one_beta() structured mu, and every zero_one_beta() zoi and coi route remain fenced and unchanged.

Fixed-kappa Gaussian mesh intercept at point-fit recovery

• spatial(1 | site, mesh = mesh) now has authenticated current-source point-recovery evidence for its exact tested fixed-domain n = 128 and n = 256 designs. A fresh 100-fit Totoro campaign retained 50/50 usable fits per rung; relative bias was -0.0301 and -0.0096, and log-scale RMSE was 0.1444 and 0.1009, with both Monte Carlo uncertainty bounds inside the frozen gates. The earlier n = 64 failure remains the lower tested boundary, so this is not a universal n >= 128 guarantee. kappa remains fixed; field-scale intervals, coverage, projected marginal-SD inference, range, slopes, non-Gaussian or bivariate meshes, anisotropy, barriers, replicated fields, and spatiotemporal fields remain unavailable.

Zero-one-beta ordinary coi random effects at point-fit recovery

• zero_one_beta() now admits the exact complete-response ML-Laplace routes bf(y ~ x, sigma ~ 1, zoi ~ 1, coi ~ 1 + (1 | id)) and bf(y ~ x, sigma ~ 1, zoi ~ 1, coi ~ x + (0 + x | id)). The slope route requires the same untransformed raw symbol in its fixed and random terms. Independent objective and gradient checks, current-source compatibility tests, and retained M = 64, 50-observations-per-group recovery rungs support population-level point recovery only. Sparse observed zeroes or ones, or weak boundary-row predictor spread, can weakly identify individual conditional modes; inspect both before interpreting them. Direct profiles, intervals, coverage, transformed or mismatched slopes, simultaneous or structured atom effects, missing responses, REML, and AGHQ remain unavailable.

Association alpha and eta intervals

• vcov() and confint() now expose two-stage Godambe-Wald uncertainty for the association-link coefficients alpha from every admitted fixed-effect, complete-pair association route. These routes are interval-feasible whenever their fit-specific covariance diagnostics pass; routes without coverage calibration warn that their intervals remain experimental. The literal-Bernoulli x ordinary-NB2 intercept route is inference-ready with caveats: a retained 16-cell high-information campaign passed its bias, availability, SE-calibration, and coverage gates. Lower-information fits warn. confint(assoc, type = "eta") transforms a constant-association interval to the bounded eta scale. predict(assoc, newdata = ..., type = "eta", se.fit = TRUE, interval = "confidence") returns row-specific eta estimates, delta-method standard errors, and transformed pointwise limits. These derived results inherit the alpha route’s evidence tier; simultaneous eta bands, profiles, random effects, missingness, weights, offsets, and REML remain outside this method.

simulate() redraws random effects (re.form) — corrected before first release

• Earlier 0.6.0 development builds held random effects frozen at their fitted values in simulate(). Every replicate reused the same û; only residual noise varied. Simulated data therefore under-represented between-group variability, and — because confint(method = "bootstrap") is driven by simulate() — parametric bootstrap intervals for models with random effects were anticonservative (too narrow). This is corrected here, within 0.6.0, before the first release: no released version ever shipped the frozen behaviour.

• simulate() gains re.form. re.form = NULL (the default) draws a fresh random-effect realisation for every replicate; re.form = NA reproduces the earlier conditional behaviour. confint() exposes the same choice for the parametric bootstrap as bootstrap_re_form, also defaulting to marginal.

• If you installed a 0.6.0 development build from GitHub, your results change. Any code depending on the earlier behaviour needs re.form = NA — and should first ask whether it wanted frozen random effects. Reproducing a fitted dataset is a legitimate use; a parametric bootstrap or a posterior predictive check is not, which is why the default is marginal rather than the argument merely being added. Bootstrap intervals computed on an earlier development build are anticonservative and should be recomputed.

• Marginal draws support ordinary grouped intercepts and slopes, intra-block correlation ((1 + x | g)), labelled cross-parameter correlation ((1 | p | id)), and phylo()/spatial()/relmat()/animal()/ phylo_interaction() structured mu effects at q = 1. Cross-trait q > 1 structured effects (including a multi-endpoint phylo_interaction()), correlated covariance blocks, corpair() regression, and modelled random-effect scale abort with an informative message rather than silently falling back to conditional simulation — a silent fallback would reproduce the defect this change fixes. Those structures can still be simulated with re.form = NA, and any bootstrap interval obtained that way is anticonservative and must not be used as coverage evidence.

• predict() is deliberately unchanged: conditional prediction at the original data is correct, and was not the defect.

• Found by the Arc B numerical audit’s score-consistency check, not by reading the code — the first Bartlett identity E[score] = 0 failed on the random-effect variance component at z = 5.36 (60 replicates) growing to 10.59 (200), and now measures z = -0.205. No certified capability-ledger cell relied on bootstrap intervals, so no evidence was retracted. See docs/design/243-marginal-simulation-and-re-form.md.

Legacy Julia cross-family extractor repair

• Legacy drmTMB_julia_xfam objects no longer let vcov(), fitted(), residuals(), stored predict(), or summary() fall through to absent fields and silently return NULL or an empty coefficient table. The bridge now reconstructs per-axis u = 0 response means and response residuals, reports point-only coefficients plus ML information criteria, and clearly errors for unavailable covariance/Wald inference. This is a compatibility repair only: Julia cross-family fitting remains deferred, rho_latent is not rho12, and no cross-family inference claim is added.

Unreleased beta association-gradient extension (superseded for intervals)

• Historical point-estimate note, superseded by Association alpha intervals above: the frozen-margin literal-Bernoulli × ordinary-NB2 beta route now also accepts an intercept-bearing fixed-effect association formula with multiple predictors, factors, interactions, and explicit transformations. It fits a row-specific latent-normal association after the two margins have been fitted and frozen. The coefficients now have interval-feasible alpha-scale Godambe-Wald uncertainty; fitted eta values and new-data predictions now have derived pointwise uncertainty as described above. This is not a direct rho12 model and supplies no simultaneous eta bands, profiles, calibrated association-regression coverage, random effects, offsets, missing association predictors, aliases, or dot expansion.

Arc 6.5 Bernoulli × Bernoulli development slice (superseded for intervals)

• Historical point-estimate note, superseded by Association alpha intervals above: associate_pairs() admits two frozen literal Bernoulli-logit margins and estimates a latent-normal eta using deterministic bivariate-normal rectangle probabilities. Its alpha-scale Godambe-Wald interval is now interval-feasible when diagnostics pass, while its retained recovery campaign remains HOLD and coverage remains uncalibrated. It is neither rho12 nor an observed-scale correlation or odds ratio.

Arc 6.4 exact bivariate Student-t development slice

• biv_student() now fits one bounded exact bivariate Student-t likelihood: fixed-effect mu1/mu2, constant Student-t scales sigma1/sigma2, one shared constant nu > 2, and constant scatter/residual correlation rho12. At finite nu, rho12 = 0 means uncorrelated residuals, not independent margins, because each pair shares the same scale-mixture draw. The route is source-tested only; random/structured effects, scale/shape/correlation predictors, partial pairs, offsets, weights, meta_V, mi(), REML, Julia, intervals, coverage, smoke, recovery, and capability promotion remain deferred.

Arc 6.3 exact bivariate lognormal development slice

• biv_lognormal() now fits one bounded exact bivariate-lognormal likelihood: fixed-effect mu1/mu2, constant log-response SDs, and a constant rho12. rho12 is the within-row residual correlation on the log-response scale, not the frozen-margin eta or a raw-scale correlation. The route accepts only complete, finite, positive pairs with unit weights; all random/structured effects, sigma/rho predictors, offsets, meta_V, mi(), REML, Julia, intervals, coverage, and capability claims remain deferred.

Cross-family association first implementation (superseded for intervals)

• Historical point-estimate note, superseded by Association alpha intervals above: post-0.6 development implements the bounded post-fit associate_pairs() interface with kernel = latent_normal() for fixed-effect Gaussian × literal-Bernoulli, Gaussian × ordinary-NB2, literal-Bernoulli × ordinary-NB2, and ordinary-NB2 × ordinary-NB2 complete pairs. It freezes stage-1 margins (including NB2 mu and sigma), then estimates latent-normal association (eta). Public alpha-scale Godambe-Wald intervals are now interval-feasible for every admitted route, with the exact Bernoulli × ordinary-NB2 intercept domain inference-ready with caveats. The staged interface still does not add mixed-family rho12; eta-scale uncertainty is a derived transformation of the alpha covariance.

First-impression formula surface (issue #776)

• (1 + x || g), the lme4/brms spelling for uncorrelated random effects, is now accepted and desugars to (1 | g) + (0 + x | g). Previously R parsed || as its own operator, the random-term parser recognised only |, and the term survived into the fixed-effect design matrix, where it aborted with 'length = N' in coercion to 'logical(1)' — a message whose N is the row count, so it read as a size problem and was not. The two spellings are now the same model: on a 20-group fixture they agree in log-likelihood and in every fixed-effect coefficient, while the correlated (1 + x | g) block stays distinct.
• This is a formula rewrite onto a route drmTMB already fits and has certified. No likelihood, no TMB change, and no new capability claim: no ledger cell changes tier, and the desugaring is inert on formulas that contain no ||.
• A categorical slope under || is rejected rather than silently accepted. In lme4 || splits by formula term rather than by design-matrix column, so a factor slope keeps its within-factor correlations and is not in fact uncorrelated; the error names the explicit two-term form instead of copying that behaviour.
• Smooth terms (s(), te(), ti(), t2()) are rejected by name before the model frame is evaluated. A reader arriving from mgcv or gamlss previously got R’s own could not find function "s"; the message now points at poly() and splines::ns() and states that penalised smooths are not implemented.
• A user-facing error no longer leaks internal roadmap vocabulary: the univariate Gaussian parameter check named “Phase 1” instead of naming the supported distributional parameters.
• The simple-grouping rule ((1 | g1/g2) and (1 | g1:g2) are not implemented) is now locked by tests on both the mu and sigma parse paths, where the guard is duplicated.

Reader-facing plotting surface complete (issue #58)

• The figure gallery now demonstrates all six public plotting functions. Alongside the existing plot_parameter_surface(), it adds worked examples of worm_plot() and qq_plot() (on a correctly specified fit paired with a deliberately mis-specified one, so a reader can see what misfit looks like), centile_chart(), plot.profile.drmTMB() via profile(), and a fitted plot_corpairs() correlation row with a computed profile interval.
• The gallery’s hand-typed illustrative figures now declare themselves as fixtures in their own titles, not only in captions, so a screenshot cannot be mistaken for a fitted result. A long-standing inverted axis label on the simulation bias panel is corrected. Documentation only; no code change.

Arc 4c ordinary mu random-slope profile coverage

• The exact independent mu random-slope cells for skew_normal() (mc-0464), tweedie() (mc-0539), and zero_one_beta() (mc-0575) are now inference_ready_with_caveats for the standard ML-Laplace profile interval. A retained 1,200-attempt-per-M Fir campaign supports a deployment floor of M=16 for each family at true slope SD 0.50. The ledger estimator remains ML; no cell earns supported status.
• This narrow result does not cover other SDs, observation counts, group grids, correlated or labelled slopes, scale/shape random effects, structured effects, REML, or AGHQ. Skew-normal retains slant-identification risk, Tweedie retains small-M zero-boundary profiles, and zero-one beta’s coverage is generator-qualified — a disclosed leak of machine-exact ones from the nominally-interior sampler means the interval claim holds only under the generator as executed; a strictly-interior rerun was specified, reviewed, and decided against (no faithful strictly-interior rendering of the intended generator exists), so this is the adjudicated terminal characterization, not a pending fix. Campaign point-bias and Wald coverage are unavailable because of a disclosed reporting defect; profile coverage is unaffected and the prospective extractor is repaired.

Beta q1 phylogenetic location intercept

• Native univariate ML beta() now admits one unlabelled intercept-only phylo(1 | species, tree = tree) term in mu. Family sigma remains a fixed-effect formula and controls phi = sigma^(-2); it is distinct from the constant latent phylogenetic location-effect SD.
• An independent dbeta() plus augmented-GMRF joint-likelihood oracle and central-difference gradient check cover the exact implementation. A fresh, predeclared 800-fit Totoro campaign retained 400 attempts at each of g = 512, m = 4 and g = 1024, m = 4; all fits converged with pdHess = TRUE. The log-latent-SD recovery gate held at g = 512 and passed only at the exact tested g = 1024, m = 4 cell.
• This is a point_fit_recovery admission only for that exact tested regime, not g >= 1024 or a universal minimum species count. Moderate-information g = 256 and g = 512 results remain explicit HOLDs. REML, q2/q4, labels, phylogenetic slopes, phylogeny in family sigma, direct sd() regression, zero_one_beta(), missing/external data, intervals, coverage, and broader Beta or all-family claims remain outside this PR.

Exact supplied-relatedness q2 REML intercept (Arc 1b-S2R)

• drmTMB(..., REML = TRUE) now admits one exact bivariate-Gaussian supplied-relatedness location cell: matching labelled relmat(1 | p | id, K = K) intercepts in mu1 and mu2. Both formulas must use the same label, group ordering, and supplied covariance matrix K; sigma1, sigma2, and rho12 must be intercept-only; response pairs must be complete; weights must equal one; and no known meta_V(), additional random effect, direct-SD formula, or corpair() regression may be present.
• An independent dense restricted-likelihood oracle matches the native TMB objective at the optimum and two displaced parameter vectors. A deliberately wrong precision orientation fails materially. The predeclared 2,400-fit Totoro campaign retained every attempt; every fit converged with pdHess = TRUE, and all structured-SD, structured-correlation, and RMSE gates passed.
• This is a point_fit_recovery admission only. Supplied precision Q, animal(), slopes, scale-side blocks, q4 or larger blocks, non-Gaussian families, intervals, coverage, AI-REML, and supported claims remain outside this arc.

Exact bivariate-spatial q2 REML intercept (Arc 1b-S1)

• drmTMB(..., REML = TRUE) now admits one exact bivariate-Gaussian coordinate-spatial location cell: matching labelled spatial(1 | p | site, coords = coords) terms in mu1 and mu2, with intercept-only sigma1, sigma2, and rho12, complete response pairs, unit weights, no known meta_V() covariance, and no additional ordinary random effect, direct-SD formula, or corpair() regression. The spatial covariance is fixed by the coordinates; range estimation is not part of this slice.
• An independent dense restricted-likelihood oracle matches the native TMB objective at the optimum and displaced parameter vectors. A predeclared 1,200-fit Totoro campaign completed without fit or convergence failures; 1,198 fits had positive-definite Hessians, and the high-information cells passed all structured-SD and latent-correlation recovery gates.
• This is a point_fit_recovery admission only. Spatial slopes, range estimation, animal-model bivariate REML, supplied-Q relatedness REML, scale-side q2, q4 or larger blocks, intervals, coverage, AI-REML, and supported claims remain outside this spatial arc. Arc 1b-S2R subsequently admits only the exact matching supplied-K relmat() q2 location-intercept cell at the same evidence tier.
• A subsequent prospective Confidence Eye campaign promoted the two direct spatial-SD cells and their joint latent spatial correlation target to inference_ready_with_caveats for exact tested M = 36 sites x 3 and H = 36 x 8 baseline-ring configurations. All 1,500 datasets and 4,500 target outcomes remained in the denominator. M and H passed jointly; L = 12 x 3 failed. Mesh intervals, estimated range, slopes, q4+, non-Gaussian spatial models, spatial scale models, derived observed correlations, other geometries or information configurations, and supported remain outside the claim.

Positive-continuous q1 structured location intercepts (Arc 3a)

• Native univariate ML now fits one unlabelled q1 structured intercept in mu for Gamma-phylo(), lognormal-phylo(), and lognormal-relmat() using K or Q. The effect enters the Gamma log-mean or lognormal log-response location predictor. Gamma-relmat() retains its existing intercept and independent one-slope route.
• Focused likelihood, extractor, prediction-decomposition, K/Q parity, and rejection-neighbour tests are green. A 6,000-fit primary Totoro campaign certified lognormal-relmat() and retained the predeclared phylogenetic intercept-RMSE HOLD. A separate, freshly seeded 2,400-fit phylogenetic addendum then passed its design-conditioned GLS-oracle and structured-field projection gates without relaxing the original threshold. All three new cells are therefore point_fit_recovery; no interval tier is implied.
• New-route slopes, labels/q2+, sigma structure, joint mu/sigma, simultaneous structured providers, spatial/animal, bivariate responses, REML, intervals, coverage, and supported claims remain outside this arc.

Exact-Gaussian REML for mean-side structured providers (Arc 1a)

• drmTMB(..., REML = TRUE) now admits pure-mu, univariate Gaussian spatial(), animal(), and relmat() terms as an unlabelled intercept or an independent intercept plus one numeric slope. These routes require sigma ~ 1 with no sigma random effect.
• Independent restricted-likelihood oracles, representation-parity fixtures, and Totoro recovery and profile campaigns support the three cells at inference_ready_with_caveats. The campaigns used spatial coordinates, animal A, and relmat K: spatial and relmat cover exactly M={8,16,32}; the animal A campaign covers one fixed M=8; all use n_each=20, where M is the number of structured levels (and hence the matrix dimension) and n_each is the number of observations per structured level. Pedigree and Ainv animal inputs and relmat Q have deterministic representation-parity evidence only, not multi-seed campaign coverage. Coverage clears the pre-specified small-sample floors but is not nominal-exact, so supported is withheld.
• The fitted structured SD scale s_j gives latent-field covariance s_j^2 K_h; a node’s marginal SD is s_j sqrt(K_h[ii]), so it equals s_j only when that diagonal entry is one. Slope-only, labelled or multiple slopes, sigma random effects, matched mu+sigma, bivariate and non-Gaussian routes remain outside this Arc 1a claim.

Residual-scale random intercepts for lognormal and Gamma (Arc 2c)

• predict(..., dpar = "sigma") now includes the fitted residual-scale random intercept for lognormal() and Gamma(link = "log") models; sigma(), printed fit summaries, and the emmeans preflight use the same corrected capability detection.
• A residual-scale (sigma) random intercept sigma ~ ... + (1 | id) is now accepted for lognormal() and Gamma(link = "log"), joining gaussian() (full) and nbinom2() (intercept-only) as the families that allow a random effect on a dispersion parameter.
• As with the mean random effects, the sigma-SD is fit by maximum likelihood with the Laplace approximation and can be biased downward when the number of groups or the per-group replication is small. Point recovery at 40 groups has -3% to -4% relative bias. The separate Arc 4a iid campaign promotes only the lognormal route to inference_ready_with_caveats for true SD 0.4, n_each=12, and exactly M={16,32,64}; coverage is mildly anti-conservative, not nominal. Gamma retains point-recovery evidence only. See docs/dev-log/simulation-artifacts/2026-07-12-arc2c-sigma-recovery/ and docs/dev-log/simulation-artifacts/2026-07-12-dg3-re-sd-coverage/README-profile-iid-v2.md. Sentinels remain in tests/testthat/test-arc2c-sigma-random-intercept.R.
• Scope (first gate): one independent sigma random intercept only. A sigma random slope, labelled covariance blocks, and combining a sigma random effect with a mu random effect in the same model remain rejected until joint recovery tests exist. The other non-Gaussian families still reject sigma random effects.

Random slopes for the intercept-only families (Arc 2b)

• One independent mu random slope (0 + x | id) is now accepted for the five families that gained a random intercept in Arc 2a: binomial(), cumulative_logit(), skew_normal(), tweedie(), and zero_one_beta(). Combined with Arc 2a, every fitted univariate family now supports both a mean random intercept and an independent mean random slope.
• Random-effect standard deviations for these families are fit by maximum likelihood with the Laplace approximation and can be biased downward when the number of groups or the per-group replication is small. All five have point recovery from a 60-seed sweep (per-family relative slope-SD bias of -2% to -9% at 40 groups). A separate Arc 4a iid campaign promoted the binomial route to inference_ready_with_caveats at true SD 0.6, 12 observations and 12 trials per observation, and exactly M={32,64}; it is coverage-backed but mildly anti-conservative rather than certified nominal. Point-recovery evidence is in docs/dev-log/simulation-artifacts/2026-07-12-arc2b-slope-recovery/, with single-seed DG2 sentinels in tests/testthat/test-arc2b-mu-random-slope.R; corrected binomial coverage evidence is in docs/dev-log/simulation-artifacts/2026-07-12-dg3-re-sd-coverage/README-profile-iid-v2.md. Later campaigns promoted cumulative-logit (mc-0227) and the three Arc 4c cells above under their own exact design-specific caveats.
• Scope: one independent mu slope only. Correlated intercept-slope blocks (1 + x | id), labelled covariance blocks (0 + x | p | id), and sigma/shape/inflation-dpar random effects remain rejected for these families.

Random intercepts for every family (Arc 2a)

• An ordinary mu random intercept (1 | group) is now accepted for the five families that previously rejected all random effects: binomial(), cumulative_logit(), skew_normal(), tweedie(), and zero_one_beta(). Every fitted family now supports at least a mean random intercept, joining the families (gaussian(), poisson(), nbinom2(), gamma(), lognormal(), beta(), beta_binomial(), student(), truncated_nbinom2()) that already did.
• Fits are by maximum likelihood (Laplace). Random-effect standard deviations are recovered at a known data-generating point (per-family DG2 sentinels in tests/testthat/test-arc2a-mu-random-intercept.R). With few or small clusters the random-effect standard deviation can be biased low under the Laplace approximation; adaptive Gauss-Hermite quadrature is the standard remedy for the non-Gaussian families and remains planned.
• Scope: mu-side intercepts only. Random slopes, sigma/shape/inflation-dpar random effects, labelled covariance blocks, and (for cumulative_logit) combining a phylogenetic effect with an ordinary intercept remain rejected in this slice.

Distributional output & adequacy layer (#747, #748)

• Every fitted family now exposes a shared fitted-distribution foundation: drm_family_dpq() (internal registry) and fitted_distribution() return per-row density (d), CDF (p), and quantile (q) closures at the fitted, fixed-effect distributional parameters. All 18 fitted model_type values are promoted to status = "reference", including bivariate biv_gaussian (marginal-only: response = 1 or 2 selects which response’s N(mu_k, sigma_k) marginal is returned; the joint distribution and rho12 are out of scope).
• residuals(fit, type = "quantile") returns Dunn-Smyth (1996) randomized quantile residuals, qnorm(F(y; theta_hat)), for every family; worm_plot() and qq_plot() draw the corresponding detrended and ordinary QQ diagnostics, with an optional nsim multi-realization seed envelope so a single randomized draw is not over-read.
• predict(fit, type = "quantile", prob = ) returns conditional response quantiles; exceedance(fit, threshold, newdata) returns Pr(Y > threshold) (or its complement) as a thin wrapper over the shared CDF; centile_chart() draws model-conditional centile curves against one covariate. All three, and the plug-in prediction intervals, carry attr(., "calibrated") <- FALSE; none of these outputs propagate theta_hat uncertainty.
• What the diagnostic detects, and what it does not (see docs/dev-log/simulation-artifacts/2026-07-12-dg3-power-arm-gated/, 400-seed gated campaign across all 18 families; tweedie: 99 of 400 seeds locally, 66/99 dispersion-arm non-convergence, full run deferred to Totoro). Under a correctly specified fixed-effect model, type-I error stays near or below the nominal rate (Type-I 0.0025-0.025 across families at alpha = 0.05; the KS+PIT statistic is conservative, so power is understated, not overstated). Under a genuine distributional shape/atom mis-specification that a family cannot reabsorb through its own free parameters – heavy tails fit as Gaussian (power 0.925-0.995), overdispersion or zero-inflation ignored by a family with no free dispersion parameter (Poisson: 0.9625-0.9825), truncation ignored and a plain count model fit instead (truncated_nbinom2 vs plain nbinom2: 1.0), a zero/one atom ignored and a plain beta fit instead (zero_one_beta/tweedie: 0.99-1.0) – power is high, typically >= 0.8 at n = 300-400 per arm. There is a genuine structural blind spot, not a bug: a mis-specification that a fitted family’s own free nuisance or dispersion parameter can absorb leaves the fitted-model residual marginally N(0,1) and is not detectable by this diagnostic – e.g. heteroscedasticity absorbed by Student-t nu (power 0.035 at n = 300, versus 1.0 for the same heteroscedasticity under Gaussian, which has no absorbing parameter), missing zero-inflation absorbed by nbinom2 sigma (power 0.035, versus 0.9625 for the same missing zero-inflation under Poisson), and fitting a plain nbinom2/beta-binomial/Tweedie to data whose TRUE dispersion actually varies with a covariate, which its own constant dispersion parameter partially soaks up (power 0.01-0.14, versus 0.81-1.0 for the same mis-specification in Gamma/beta/lognormal, which lack a matching absorbing structure). Detecting an absorbed mis-specification is a mean-structure diagnostic’s job, not this one’s. Zero-inflation/hurdle/ zero-one-inflation mechanism mis-specification (a constant inflation probability fit when it truly varies with a covariate) splits into two patterns under the n-ladder (tested to n = 3000): for hurdle_nbinom2/ zero_one_beta power stays flat at or below about 0.01 at every n – a genuine structural blind spot; for zi_nbinom2/zi_poisson power rises with n (to about 0.11/0.06 at n = 3000), so the marginal is not identical under the mechanism mis-spec, but power stays far below the >= 0.8 detectable benchmark even at n = 3000, so it remains impractical to detect at realistic sample sizes. Neither pattern should be relied on as an adequacy check for the mechanism. gamma-vs-lognormal wrong-family detection is sample-size limited rather than structurally blind: power rises from about 0.19 at n = 300 to 0.79 at n = 1000 and 1.0 at n = 3000, so this specific mis-specification needs n well above 1000 to be reliably caught.
• Honest-scope invariants apply throughout: intervals are labelled calibrated = FALSE; a pass is worded “no detectable departure”, never “adequate” or “the model is correct”; and a distributional-output/adequacy (DG) tick on a family never changes or implies anything about that family’s own inference-tier status (e.g. skew-normal’s diagnostic_hold fit-quality status is unaffected by its DG2/DG3 promotion) – see tests/testthat/test-dg-firewall.R.
• This is fixed-effect adequacy only: for random-effect or structured fits, quantile residuals are conditional on the fixed-effect prediction, not marginal, so a departure (or its absence) is evidence about fixed-effect adequacy only. Calibrated coverage (DG4/DG5), uncertainty beyond theta_hat, random-effect/structured residual adequacy, and bivariate joint (non-marginal) outputs remain separately authorized future work.

Missing responses: MR-T7 certification

• The generated capability ledger and live runtime oracle now reconcile all 18 fitted response routes at G3 recovery-verified, with zero G0 routes. The capability page, missing-data article, design inventory, NEWS, roadmap, and machine-readable evidence are regenerated from the same route state. This closeout adds no family, formula grammar, estimator, interval, or coverage claim: each tick remains bounded to its documented fixed, random, or structured route, and G4/G5 remain separate future evidence tiers.

Missing responses: MR-T6 count mixtures

• response = "include" now masks fixed-effect zero-inflated Poisson, zero-inflated NB2, and hurdle NB2 responses. Each route guards its complete zero-or-positive mixture contribution before response classification, uses observed-only starts, and has separate missing-zero and missing-positive parity tests. Zero-versus-positive sentinel retapes, row/extractor contracts, and exact fixed-seed 25% MCAR recovery promote all three routes to G3. All 18 fitted response routes are now G3 recovery-verified for their documented masking slice; random/structured mixture routes, response plus mi(), REML, intervals, and coverage remain outside this arc.

Missing responses: MR-T5 truncated counts

• response = "include" now masks positive-count responses for the non-hurdle truncated_nbinom2() route. The complete NB2 density and zero-truncation normalization are skipped together for masked rows; positive sentinels, observed-only starts, row/extractor contracts, and exact fixed-seed 25% MCAR recovery promote the ordinary mu random-intercept route to G3. Hurdle, sigma-random, structured, response-plus-mi(), interval, and coverage claims remain outside this tranche.

Missing responses: MR-T4 encoded responses

• response = "include" now masks beta-binomial and cumulative-logit responses. A missing success or failure count masks the entire beta-binomial row, including its derived trials; ordered-factor levels remain declared and any observed subset with an empty category rejects before cutpoints are built. Coordinated encoded-sentinel retapes, row/extractor contracts, and exact fixed-seed 25% MCAR recovery promote both routes to G3. Integer ordinal masking, broader random/structured routes, response plus mi(), intervals, and coverage remain outside this tranche.

Missing responses: MR-T3 atom and boundary families

• response = "include" now masks missing Tweedie and zero-one beta responses. Tweedie tests retape masked rows as a zero atom and a positive continuous value; zero-one beta tests zero and one atoms against an interior value. Observed-only starts, full mixture guards, row/extractor contracts, and exact fixed-seed 25% MCAR recovery promote both fixed-effect routes to G3. Random effects, structured effects, response plus mi(), intervals, and coverage remain outside this tranche.

Missing responses: MR-T2 continuous families

• response = "include" now masks missing Student-t, skew-normal, lognormal, and Gamma responses. Plain data-time likelihood guards prevent masked values from reaching density or positive-support transformations; observed-only starts, direct sentinel invariance, row/extractor contracts, and fixed-seed 25% MCAR recovery promote all four routes to G3. Student-t, lognormal, and Gamma are verified through their ordinary random-intercept routes; skew-normal remains fixed-effect only. This does not promote structured effects, intervals, or coverage.

Missing responses: MR-T1 verification

• The six previously admitted response = "include" routes—univariate and bivariate Gaussian, binomial, Poisson, NB2, and beta—now share direct retaped-sentinel tests, original-row and extractor contracts, and fixed-seed 25% MCAR recovery tests. Univariate residuals are NA on masked response rows while fitted values retain the original row length. These tests promote the six routes to the capability ledger’s G3 recovery-verified tier; they do not claim interval calibration or coverage.

drmTMB 0.5.0

drmTMB 0.5.0 was numbered as the intended first CRAN release (not 1.0) — historical note, superseded twice: 0.5.0 was never submitted to or accepted by CRAN and was retired in favour of a more capable 0.6.0-class cut; and 0.6 in turn is the development cycle, never submitted. The first CRAN submission is 0.7.0 (decided 2026-07-25). The honest version number reflects that much of the family and inference surface is still scaffolded or recovery-grade. Throughout the dev-log and the “Q-Series v1.0” ledger, “v1.0” is reserved for the later maturity milestone that 0.5.0 deliberately does not yet claim. This entry accumulates the 0.4.x development cycle (docs/ledger alignment, non-Gaussian coverage validation, the missing-data non-Gaussian arc) into the release; earlier tagged development lines appear below.

Missing data: non-Gaussian responses and predictors

The likelihood-based missing-data layer now extends beyond Gaussian responses. Both modes are validated per family against single sources of truth (drm_missing_response_families(), drm_missing_predictor_families()), and an anti-drift test asserts that every family outside those allow-lists still rejects loudly, so an unsupported request never silently degrades to a wrong likelihood. See vignette("missing-data") for the full capability matrix.

• Missing-response masking (FIML) for non-Gaussian responses. missing = miss_control(response = "include") now marginalises missing responses out of the joint likelihood for binomial(), poisson(), nbinom2(), and beta() fits, in addition to the existing univariate and bivariate Gaussian routes. Masked rows keep their complete predictors and row identity but contribute no response density (a plain data guard in the TMB kernel, so the placeholder is never taped). Valid under ignorable (MCAR/MAR) missingness.

• Missing-predictor mi() for non-Gaussian responses. missing = miss_control(predictor = "model") with an impute model now supports one binary (Bernoulli/logit) missing predictor on binomial(), nbinom2(), and beta() responses, joining the existing Poisson route. The missing predictor is marginalised by an exact 2-point sum inside the same joint likelihood, with the response density carrying its family dispersion (nbinom2() size = exp(-2*log_sigma); beta() phi = exp(-2*log_sigma)). Point-fit recovery of the mean, dispersion, and predictor-model coefficients is tested at scale for each family.

• Pluggable response-density leaf. The mi() quadrature now routes each family’s response density through one shared kernel (drm_response_log_density), so a non-Gaussian response reuses the same integration loop. The Gaussian extraction was a byte-identical refactor (verified by golden capture on the log-likelihood, gradient, and objective), and the per-family leaves replicate their inline densities exactly, including the beta() boundary nudge and shape floor.

Bug fixes

• Structured sigma random effects for family = nbinom2() (sigma ~ phylo()/spatial()/animal()/relmat()) now correctly modify the scale predictor. They were previously applied to the mean predictor (the TMB kernel’s model_type == 7 branch lacked the scale-side dispatch the beta family already had), so a sigma ~ phylo(...) fit silently matched a mean-phylo fit while reporting a *_sigma SD. Point-fit recovery is now verified (tests/testthat/test-nbinom2-sigma-structured-recovery.R); intervals and coverage remain out of scope (recovery-grade).

Coverage validation

• A simulation coverage campaign (400 seeds, n-ladder 50–800; see docs/dev-log/simulation-artifacts/2026-07-09-nongaussian-unstructured-coverage-pilot/) confirms that unstructured (fixed-effect) non-Gaussian confidence intervals are calibrated. The mean coefficients of binomial(), poisson(), beta(), and nbinom2() — including rare-event and low-count stress — and the location-scale sigma coefficients of nbinom2() all show finite-rate ≈ 1.0 and near-nominal Wald coverage. beta() location-scale intervals are calibrated for interior proportions; exact 0/1 observations require zero_one_beta().

Documentation and release-ledger alignment

• README.md, ROADMAP.md, and docs/dev-log/known-limitations.md now state the exact REML structured-effect boundary shipped across 0.2.0/0.3.0: univariate Gaussian REML accepts phylogenetic mean-side, scale-side, and matched q2 mean-and-scale blocks, plus univariate spatial/animal/relmat scale-side blocks. Arc 1a additionally admits the exact pure-mu spatial/animal/relmat intercept and independent one-slope cells over the documented discrete domains; other non-phylogenetic mean-side or mixed mean+scale structured effects, sparse-fixed designs, Gaussian row aggregation, and ordinary direct-SD formulae remain rejected. Bivariate Gaussian REML accepts phylogenetic structured effects in every covariance layout, including the dense q4 block, and rejects spatial/animal/relmat entirely. REML remains rejected outright for every non-Gaussian family.
• ROADMAP.md corrects the Q-Series inference_ready anchor count from five rows to eight, adding the three q1 mu:(Intercept) anchors (phylo, spatial, relmat) that the release ledger already carried but the roadmap text had not listed. Two of the eight rows – the phylo and relmat q2 mu1:x/mu2:x slope-SD rows – are inference_ready only through the bias-corrected confint() channel; their raw uncorrected Wald intervals fail coverage. No structured row is supported, and non-Gaussian structured rows remain point-recovery evidence only, with no intervals, coverage, or supported claim.
• docs/dev-log/known-limitations.md records that nbinom2() structured sigma terms (phylo/spatial/animal/relmat) now correctly target the scale predictor log_sigma (the routing fix announced under Bug fixes above); earlier versions mis-targeted the mean predictor. These four rows are recovery-grade only – point-fit recovery is verified, but intervals and coverage remain out of scope.

Inference guidance

• Historical note, superseded by the cell-specific 0.6.0 guidance above. This release originally described profile likelihood as the headline method for structured covariance targets. The retained evidence does not support that blanket recommendation: q1 mu and the exact phylo/relmat slope-only q2 mu1:x/mu2:x SD rows use the default location-axis bias-corrected, small-sample-t Wald channel; q1 sigma uses raw uncorrected log-SD Wald-z evidence and its profile channel is diagnostic-only at g = 8; Arc 1a REML uses direct structured-SD profiles only over its tested discrete domains. A target appearing in profile_targets() means that it can be computed, not that its profile interval is validated for reporting.

drmTMB 0.3.0

Large direct-SD models: uncertainty no longer scales with the square of the group count

Reported by Ayumi Mizuno on a 10,440-tip bivariate phylogenetic fit, where TMB::sdreport() exhausted 48 GB of memory.

• Breaking (default change). A direct-SD surface (sd(group, level = "phylogenetic") ~ ., formerly sd_phylo()) previously ADREPORTed one standard deviation per group, so the joint ADREPORT covariance was n_group x n_group. Under REML = TRUE the fixed effects are integrated into the Laplace random block and vcov() reads exactly that covariance, so a bivariate fit at ten thousand tips needed roughly 14 GB for it alone. Those per-group standard errors are now opt-in via drm_control(se_group_sd = TRUE). The fitted per-group standard deviations themselves are unchanged and always available. Parameter standard errors, vcov(), summary(), and pdHess now work under REML at that scale.

• New drm_control(se_report_covariance = ) and drm_control(se_skip_delta_method = ) pass through to the getReportCovariance and skip.delta.method arguments of TMB::sdreport(), for further control over the cost of uncertainty on large models.

• REML = TRUE no longer rejects an explicitly-passed missing = control when the data contain no missing values. The gate tested the setting rather than whether the missing-data engine actually engages, and miss_control(response = "include") is an exact no-op on complete-case data. REML combined with a missing-data engine that genuinely engages is still rejected, as that combination remains unvalidated.

• sigma ~ z + (1 + x | p | id) again reports the specific “labelled residual-scale random-slope covariance blocks are not implemented yet” error rather than a generic shape error. The behaviour (rejection) is unchanged.

Unified sd(..., level = ) scale grammar

• sd(group, level = "phylogenetic") is the new generic spelling for the phylogenetic direct-SD targets (sd1(...) / sd2(...) for the bivariate endpoints). The legacy sd_phylo() / sd_phylo1() / sd_phylo2() spellings are soft-deprecated: they keep working and emit a one-time deprecation warning. Reserved level values ("spatial", "animal", "relmat") are parsed but not yet implemented.

More REML coverage (Gaussian location-scale)

Restricted maximum likelihood now covers substantially more of the location-scale family, debiasing scale-side variance components with adequate within-group replication. Every combination admitted under REML is also admitted under ML (docs/dev-log/ml-reml-coverage-2026-07-07.md).

• Matched mean-and-scale phylogenetic block (q2) under REML. A univariate mu + sigma model with a correlated phylo(1 | p | id) block is now admitted; a sample-size ladder shows REML is less biased than ML (N >= 250 to identify, N >= 1000 for the location-scale correlation). This supersedes the earlier small-N “REML degrades the mean” verdict.

• Block-diagonal bivariate location-scale phylogenetic layout under REML. A phylo mean block and a phylo scale block with distinct labels (1 | p | id on the means, 1 | ps | id on the scales) are admitted; the scale-side random phylo is identifiable with per-group replication (it collapses at one observation per species, where a fixed sd_phylo() scale should be used instead). At this intermediate point the dense block stayed rejected; the later dense-q4 entry below supersedes that state with recovery evidence at adequate information.

• Ordinary sigma random effects under REML. A residual-scale random intercept (1 | id), an independent random slope (0 + x | id), the correlated mean-scale block (1 | p | id), and a bivariate labelled scale-side block (1 | s | id) are now admitted; REML debiases the scale-side variance component with adequate within-group replication (at very low replication it can underperform ML).

• Dense (unstructured) q4 phylogenetic location-scale block under REML. The previous “sign-flip” verdict is superseded: the DGP-to-endpoint mapping is correct (a single nonzero simulated correlation lands on the right pair with the right sign), and the apparent flip was an under-powered fit whose variance component collapsed. With adequate information (roughly n_tip >= 200 and per-species replication n_each >= 10) the dense q4 converges and recovers, and REML is strictly better than ML there – higher convergence/pdHess rate and variance components debiased toward truth. At one observation per species it still collapses; use the block-diagonal layout or a fixed sd(level = "phylogenetic") scale.

• Bivariate mean-scale random-effect correlations and q > 2 labelled location covariance blocks under REML. Both are now admitted; REML is consistently less biased than ML on the block standard deviations. ML/REML parity is now complete for every implemented cell (docs/dev-log/ml-reml-coverage-2026-07-07.md).

• Scale-side spatial / animal / relatedness structured effects under REML. sigma ~ spatial(...), sigma ~ animal(...), and sigma ~ relmat(...) now fit under REML = TRUE. A recovery + coverage campaign shows REML debiases the scale-side intercept standard deviation in every cell (bias approaching zero as the group count grows) and profile-CI coverage clears the small-sample floor. Mean-side non-phylogenetic structured effects under REML were unvalidated and rejected at the time of this entry. Superseded in 0.6.0: the Arc 1a spatial, animal, and relmat() unlabelled mean-intercept and independent intercept-plus-one-numeric- slope REML cells are now admitted only over their recorded discrete recovery domains; slope-only, labelled, multiple-slope, q > 1, simultaneous-provider, and adjacent structured cells remain rejected.

• Degrees of freedom under REML now count the marginalised scale fixed effects. A scale-side REML fit marginalises beta_sigma as well as beta_mu; logLik()’s df (and therefore AIC() / BIC()) now counts both, matching the ML parameter count. Fits without a sigma variance component are unchanged.

New: check_drm() diagnostics for weak identification and direct-SD surfaces

• phylo_mu_diagnostics no longer reports a false error for a fitted sd(group, level = "phylogenetic") ~ . surface (which has no scalar phylogenetic standard deviation). It now summarises the fitted per-group SD surface and errors only on genuinely non-finite or non-positive fitted standard deviations.

• New standard_errors_inflated check flags a finite-but-inflated Wald standard error on a converged, positive-definite-Hessian fit – the signature of a weakly identified, near-flat direction such as a boundary correlation. The bivariate phylogenetic mean-mean boundary warning now names the same symptom in words. A clean pdHess is necessary, not sufficient.

• New “Choosing between maximum likelihood and REML” guidance in the Improving convergence article: ML is the default; REML’s p / n correction to variance components matters mainly at small group counts, and it leaves the mean coefficients essentially unchanged.

New: correlated residual-scale random slopes

• sigma ~ x + (1 + x | id) – a correlated residual-scale intercept-slope block – is now implemented (and the multi-slope (1 + x1 + x2 | id) generalisation). Previously only independent residual-scale slopes ((0 + x | id)) were supported. The univariate TMB likelihood now applies the same-dpar correlation conditioning to the sigma random effects, mirroring the mu side. Recovery of the intercept SD, slope SD, and their correlation is validated against known truth.

• Consequently the ordinary two-level DHGLM with correlated random slopes on BOTH the location and the scale – y ~ x + (1 + x | id) with sigma ~ x + (1 + x | id) – now fits, under ML and REML. The remaining piece of the full q12 is the labelled cross-formula mu-sigma slope block (the mean-scale slope cross-correlation), which is still planned.

drmTMB 0.2.0

REML for Gaussian and bivariate-Gaussian location-scale models

Restricted maximum likelihood (REML = TRUE) now covers more of the phylogenetic location-scale model family, debiasing the variance components and giving better-conditioned, honest scale-side standard errors. Validated by exact restricted-likelihood references and known-truth recovery ladders (docs/design/221-native-reml-finish.md); the native REML test suite is green.

• Bivariate Gaussian REML with phylogenetic / random location effects. drmTMB(..., REML = TRUE) now fits biv_gaussian() models whose means carry correlated phylo() (or ordinary) random effects (the “correlate the means” model), matching an exact bivariate restricted-likelihood reference. A sample-size recovery ladder shows REML is less downward-biased than ML on the variance components at every sample size, and its standard errors track ML’s. At the time of this 0.2.0 entry, scale-side random effects, matched mean-and-scale phylogenetic effects, and q > 2 labelled covariance blocks remained rejected. Later 0.3.0/current entries supersede that boundary with row-specific point-fit or recovery evidence.

• Phylogenetic direct-SD scale (sd_phylo(...) ~ predictors) under REML. The heteroscedastic phylogenetic-variance model – a predictor (e.g. climate) on the phylogenetic SD – is now admitted under REML for univariate and bivariate Gaussian models, matching an exact restricted-likelihood reference. This is the scale side of the corrected ecogeographic location-scale model.

• Correct REML standard errors for direct-SD coefficients. vcov() and summary() previously returned NA standard errors for the sd_phylo coefficients under REML (they are absent from the sdreport ADREPORT joint covariance); they now fall back to the fixed-parameter covariance and report finite Wald standard errors.

Native phylogenetic location-scale fits with debiased variance components reduce the need for an external Bayesian comparator (e.g. MCMCglmm) for this workflow.

drmTMB 0.1.4

Current development claims in this NEWS section follow the finish-plan claim registry in docs/design/168-r-julia-finish-capability-matrix.md; fitted, planned, unsupported, and release-gate language should not be read more broadly than that matrix.

• drmTMB() now fits three non-count family structured mu one-slope cells as native point-fit/extractor recovery-only routes: Gamma() with relmat(1 + x | id, K = K), student() with spatial(1 + x | id, coords = coords), and beta() with animal(1 + x | id, pedigree = ped). Each extends the existing structured intercept gate to an unlabelled intercept-plus-one-slope term with no compiled-code change. On a crossed n_lvl in {10,20,30} x 30-seed ladder, a null-slope separability control, and a non-identity AR(1) relatedness check, both variance components recover with RMSE falling as levels increase (Gamma and beta 90/90 converged with positive-definite Hessian; Student-t 83/90, so use n_levels >= 20). This is recovery-only: labelled or multiple structured slopes, scale/shape/zero-inflation structured slopes, other families, intervals, coverage, inference_ready, supported, REML, AI-REML, and bridge support remain planned.

• Simultaneous two-provider structured count mu is now admitted at point-fit/recovery for NB2: nbinom2() with spatial(1 | site, coords = coords) + relmat(1 | id, Q = Q) on a crossed site x id design now builds and surfaces both structured fields (ranef() shows spatial_mu and relmat_mu; both SDs are direct log_sd_phylo/log_sd_phylo2 profile targets). On the crossed ladder both fixed-covariance variance components recover with a positive-definite Hessian, and a non-crossed control shows the separability requirement (site and id must vary independently). This is recovery evidence only — it does not authorize interval reliability, coverage, inference_ready, STAN cross-check, REML, AI-REML, bridge parity, or supported. Joint identifiability rests on the crossed design.

• Structured q12 two-slope all-four covariance is now admitted at point-fit/recovery for phylo(), spatial(), animal(), and relmat(): (1 + x + z | p | id) on mu1/mu2/sigma1/sigma2 builds a twelve-endpoint (66-correlation) among-trait covariance that recovers a known covariance at adequate sample size. pdHess=FALSE is expected here (the 66-correlation block is weakly identified) and is not failure: the twelve SDs are direct profile targets and the 66 correlations route through profile/bootstrap (ELR excluded). This is recovery evidence only — it does not authorize interval reliability, coverage, STAN cross-check, REML, AI-REML, bridge parity, or supported. With this admission every Gaussian structured-random-effect row now holds a v1.0 basic-working-or-better row-accounting role; the remaining rows outside that practical surface are non-Gaussian.

• Structured q6 two-slope location covariance is now admitted at point-fit/recovery for phylo(), spatial(), animal(), and relmat(): bf(mu1 = y1 ~ x + z + phylo(1 + x + z | p | id, tree = tree), mu2 = ..., sigma1 = ~1, sigma2 = ~1, rho12 = ~1) builds a six-endpoint (15-correlation) among-trait covariance that recovers a known covariance with a positive-definite Hessian at adequate sample size. This is recovery evidence only: the six SDs are direct profile targets and the fifteen correlations are derived (no Wald interval), so it does not authorize interval reliability, coverage, STAN cross-check, REML, AI-REML, bridge parity, the structured q8 rows, or supported wording.

• The Q-Series v1.0 release status is now generated from the 104-row support-cell board and recorded in docs/dev-log/release-audits/q-series-v1-release-status.md. It separates implemented/basic-working Gaussian structured-effect rows, 27 non-Gaussian recovery rows, and 10 non-Gaussian diagnostic-only rows from post-v1.0 inference_ready and supported validation. This is release-planning evidence only; it does not authorize coverage, q4/q8 promotion, broad bridge support, REML, AI-REML, or public-support wording.

• The Q-Series v1.0 practical surface now includes ten row-specific diagnostic-only gates outside the ordinary mu lanes: Student-t intercept-only mu ~ spatial(1 | id, coords = coords), Student-t nu ~ phylo(1 | id, tree = tree), cumulative-logit ordinal mu ~ phylo(1 | id, tree = tree), truncated-NB2 hurdle hu ~ relmat(1 | id, Q = Q), zero-inflated Poisson zi ~ spatial(1 | id, coords = coords), zero-inflated Poisson fixed-zi mu ~ spatial(1 | id, coords = coords), zero-inflated NB2 fixed-zi mu ~ spatial(1 | id, coords = coords), Poisson slope-only mu ~ spatial(0 + x | site, coords = coords), Poisson labelled-scalar mu ~ spatial(1 | p | site, coords = coords), and Poisson mu ~ spatial(1 | site, coords = coords) + (1 | id). These rows establish fit/extractor feasibility but not point-estimate recovery and are not interval, coverage, inference_ready, supported, bridge, REML, AI-REML, broad shape/inflation/ordinal/structured non-Gaussian support evidence, or neighbouring-row evidence.

• truncated_nbinom2() hurdle models now fit the row-specific Q-Series v1.0 hu ~ relmat(1 | id, Q = Q) local gate. The fitted relatedness-field SD for the hurdle probability is exposed through sdpars$hu and ranef("relmat_hu"). This is local fit-only/extractor evidence; hurdle slopes, labelled covariance, broader hurdle structured effects, intervals, coverage, inference_ready, supported, REML, AI-REML, and bridge support remain closed.

• cumulative_logit() now fits the row-specific Q-Series v1.0 ordinal phylogenetic mu intercept gate, for formulas such as bf(score ~ x + phylo(1 | species, tree = tree)). The fitted phylogenetic ordinal location SD is exposed through sdpars$mu, ranef("phylo_mu"), and a direct profile_targets() row. This is local fit-only/extractor evidence; ordinal slopes, scale/discrimination formulas, bivariate ordinal models, intervals, coverage, inference_ready, and supported status remain planned.

• nbinom2() now fits the row-specific Q-Series v1.0 zero-inflated NB2 fixed-zi spatial mu intercept gate, for formulas such as bf(count ~ x + spatial(1 | site, coords = coords), sigma ~ 1, zi ~ 1). The fitted spatial mu SD is exposed through sdpars$mu, ranef("spatial_mu"), and a direct profile_targets() row. This is local fit-only/extractor evidence; zero-inflated NB2 structured slopes, labels, simultaneous providers, structured zi, structured sigma, intervals, coverage, inference_ready, supported, REML, AI-REML, and bridge support remain closed.

• confint() now applies the small-sample t(g - 1) width plus the simulation-calibrated log(g/(g - 1)) centre shift by default for location-axis structured random-effect SD targets. This moves only the phylo and relmat bivariate q2 mu1:x/mu2:x slope SD cells to inference_ready for interval and coverage status; supported remains withheld because the engine grids still measure right-tail miss asymmetry and g-dependence. Spatial q2, animal q2, q4/q8, count, and non-Gaussian structured rows remain separate future arcs.

• The exact Gaussian q1 sigma one-slope phylo(), animal(), and relmat() rows are now inference_ready under the raw uncorrected log-SD Wald-z interval channel. The Nibi top-up, banked SR475 slope grid, and local animal SR1000 reconciliation show 100% fit/pdHess pass rates, Wald finite rates at or above 0.953, and Wald MCSE at or below 0.01; the caveat is that one-sided misses are asymmetric and sigma slope SDs over-cover, so this is not supported. Profile intervals remain diagnostic-only at deployment g=8 for low-finite sigma targets, and the location-axis bias+t correction does not apply to sigma.

• drmTMB() now fits ordinary Poisson and NB2 structured mu one-slope count cells with unlabelled phylo(1 + x | ...), fixed-covariance spatial(1 + x | ...), animal(1 + x | ...), and relmat(1 + x | ...) terms. These are native TMB ML/Laplace point-fit and extractor cells for non-zero-inflated count means only. Exact q1 NB2 structured sigma intercept-plus-one-slope routes for those four providers are separately fitted at recovery grade; pure or multiple structured count slopes, labelled count covariance, zero-inflated structured effects beyond the exact Poisson spatial zi, Poisson fixed-zi spatial mu, and NB2 fixed-zi spatial mu local-fit gates, richer or labelled NB2 structured sigma, q2/q4 count covariance, bridge support, REML, AI-REML, intervals, coverage, and public-support promotion remain planned or unsupported.

• biv_gaussian() now fits the first phylogenetic all-four one-slope structured block when the same labelled phylo(1 + x | p | species, tree = tree) term appears in mu1, mu2, sigma1, and sigma2. The fitted block exposes eight endpoint SDs and 28 derived latent phylogenetic correlations through sdpars$mu, corpars$phylo, corpairs(level = "phylogenetic"), summary()$covariance, profile_targets(), and structured_effects(). This is native ML point-fit/extractor evidence plus deterministic same-target fixture parity for the exact shared-label phylo cell only; block-diagonal layouts, broad bridge support beyond the fixture, intervals, coverage, REML, AI-REML, and public-support promotion remain planned.

• biv_gaussian() now also fits the first fixed-covariance spatial all-four one-slope structured block when the same labelled spatial(1 + x | p | site, coords = coords) term appears in mu1, mu2, sigma1, and sigma2. The fitted block exposes the same eight endpoint SDs and 28 derived latent spatial correlations through the standard covariance extractors. This is native ML point-fit/extractor evidence plus deterministic same-target fixture parity for the exact fixed-covariance spatial cell only; range-estimating spatial support, block-diagonal layouts, broad bridge support beyond the fixture, intervals, coverage, REML, AI-REML, and public-support promotion remain planned.

• biv_gaussian() now fits the corresponding exact A-matrix animal and K/Q lower-level relatedness all-four one-slope blocks when the same labelled animal(1 + x | p | id, A = A) or relmat(1 + x | p | id, K/Q = ...) term appears in mu1, mu2, sigma1, and sigma2. These cells expose eight endpoint SDs and 28 derived latent animal or relatedness correlations through the standard covariance extractors. This is native ML point-fit/extractor evidence plus deterministic same-target fixture parity for the exact A-matrix animal and K-matrix relmat cells only; pedigree/Ainv bridge marshalling, relmat Q bridge marshalling, block-diagonal layouts, broad bridge support beyond the fixture, intervals, coverage, REML, AI-REML, and public-support promotion remain planned.

• The residual rho12 correlation now uses the same 0.999999 (six-nines) tanh guard as every other latent correlation in the package, instead of an inconsistent 0.99999999 (eight-nines) bound (flagged in review). The guard is far from any realistic correlation, so fitted rho12 values are unchanged to about seven decimal places; this is a numerical-consistency fix, not a behaviour change for interior fits.

• drm_phylo_penalty_sweep() runs a penalized (MAP) phylogenetic fit across a range of cor_sd correlation-penalty values and returns a tidy sensitivity summary (convergence, pdHess, logLik per cor_sd) plus the fitted objects for extracting the couplings. This turns the mandatory prior-sensitivity sweep – the check of whether a weakly identified coupling is data-informed (stable across cor_sd) or prior-shaped (tracks cor_sd) – into a single call; there is no universal cor_sd. (drm_phylo_penalty() and the new sweep are now both in the pkgdown reference.)

• check_drm() now reports a logsigma_clamp_active row that flags when the log(sigma) clamp is active at the optimum – the diagnostic-surface complement to the fit-time clamp-active warning – so a clamp-bound fit is visible in the standard diagnostic table (a note when the TMB object was dropped, an ok otherwise).

• check_drm() now warns when q2 random-effect covariance correlations are close to +/-1, including univariate mu/sigma covariance and bivariate q2 mu/mu, sigma/sigma, and same-response mu/sigma covariance rows. The diagnostic prints the fitted rho_abs and rho_boundary so a converged, positive-Hessian fit at the correlation guard is visible rather than silently labelled ok.

• check_drm() now reports fitted-boundary diagnostics for bivariate coordinate-spatial, animal(), and relmat() q2 location covariance rows (biv_spatial_q2_covariance, biv_animal_q2_covariance, and biv_relmat_q2_covariance). This changes diagnostic visibility only; it does not change fitting behaviour or promote structured q2 recovery, interval, or power claims.

• biv_gaussian() now fits the first structured slope-only q=2 mu1/mu2 covariance cells for phylo(0 + x | p | species, tree = tree), fixed-covariance spatial(0 + x | p | site, coords = coords), animal(0 + x | p | id, A/Ainv = ...), and relmat(0 + x | p | id, K/Q = ...). The fitted SDs and cor(mu1:x,mu2:x | p | group) row are coefficient-aware in sdpars$mu, corpars, corpairs(), summary()$covariance, profile_targets(), and structured_effects(). This is native point-fit/extractor evidence plus deterministic same-target fixture parity for the exact slope-only q2 cells; it does not by itself promote separate intercept-plus-slope structured q4/q8 covariance cells, broad bridge support, interval reliability, coverage, REML, or AI-REML.

• drmTMB() now rejects combining REML = TRUE with penalty = (a penalized / MAP fit): a restricted-likelihood estimator and a maximum-a-posteriori estimator are different estimators of the variance components, so the combination is undefined.

• The Student-t nu documentation now states explicitly that the nu > 2 (finite-variance) bound, required by the sigma = SD contract, means the family cannot represent the very heavy tails of nu <= 2 (e.g. Cauchy); check_drm() warns as nu approaches the boundary.

• drmTMB(..., REML = TRUE) now fits bivariate Gaussian fixed-effect location models (mu1/mu2), marginalising both mean blocks (beta_mu1, beta_mu2) for an unbiased residual covariance. Validated against an exact restricted-likelihood reference (the OLS-residual covariance with the n - p correction): sigma1, sigma2, rho12, and both mean-coefficient blocks match, and the degrees of freedom count both marginalised blocks. At the time of this entry, bivariate random-effect and structured (phylo) means under REML were a later slice and rejected; later 0.2.0/0.3.0 entries supersede that boundary with row-specific q2/q4 evidence.

• REML now supports a fixed-effect heteroscedastic residual (sigma ~ predictors), not just an intercept-only sigma. REML restricts the likelihood for the mean fixed effects regardless of the scale model, so a Gaussian with residual variance V = diag(sigma_i^2) + random-effect covariance has an exact restricted likelihood; drmTMB’s REML estimates match a hand-computed restricted-likelihood reference (random-effect SD, the sigma coefficients, and the mean coefficients), and the degrees of freedom count the marginalised mean fixed effects. At the time of this entry, scale-side random effects under REML remained rejected; later 0.3.0/current rows supersede that boundary with exact point-fit or recovery admissions.

• drmTMB(..., REML = TRUE) now fits mean-side phylogenetic location models – a phylo() term on mu with an intercept-only sigma – extending REML beyond the ordinary-random-effect slice. REML restricts the likelihood for the mean fixed effects (TMB marginalises beta_mu through its exact Gaussian Laplace step), giving a less downward-biased phylogenetic variance component. Estimates match a hand-computed restricted Gaussian likelihood (phylogenetic SD, residual sigma, and the mean coefficients). At the time of this historical entry, scale-side structured effects and non-phylogenetic structured effects (spatial, animal, relatedness) under REML remained rejected. Superseded in 0.6.0: the later scale-side work and Arc 1a admit only the exact tested cells and discrete domains named in the current capability guide; neighbouring slopes, q > 1 blocks, simultaneous providers, and untested cells remain rejected.

• drm_control(fallback_optimizer = ) adds an opt-in fallback optimizer (an [stats::optim()] method such as "BFGS") that is tried as a final attempt when no nlminb() preset converges. A different algorithm sometimes succeeds on a numerically awkward but identified problem. It is off by default (NULL), so the default fit uses only the nlminb() preset ladder; when enabled, the fallback attempt is recorded in fit$optimizer_attempts like any preset.

• drm_control(multi_start = K) enables multi-start fitting: each optimizer preset is run from K starting points – the principled start plus K - 1 reproducibly perturbed starts – and the lowest-objective result is kept. This helps weakly identified models escape poor local optima. It is opt-in: multi_start = 1 (the default) is the single-start fit and is unchanged, and the perturbations use a fixed internal seed with the caller’s random stream saved and restored, so fits stay reproducible.

• The optimizer now escalates its preset ladder (default -> careful -> robust) when a preset returns a non-converged result, not only when it throws an error. Previously drmTMB() accepted the first preset that did not error – even a false convergence (convergence != 0) or a non-finite objective – so the careful and robust presets were effectively unreachable for the very cases they exist for. Now a non-converged attempt is recorded and the next preset is tried; the first cleanly-converged attempt is returned, or the best (lowest-objective) attempt if none converge. A clean first attempt is unchanged (no escalation, no warning). The full ladder is recorded in fit$optimizer_attempts, and the fit-time convergence warning now points there instead of suggesting a manual robust refit (the ladder is tried automatically).

• The Gaussian sigma-slope starting values are more robust for strong scale-heterogeneity models. The start heuristic previously discarded all scale slopes whenever its log-absolute-residual regression looked too large, handing a legitimately steep sigma ~ x model a flat intercept-only starting point; it now shrinks an over-large slope start toward zero (direction preserved, magnitude bounded) instead. This only changes the optimizer’s starting point, not the objective, so converged fits are unchanged; it gives the optimizer a better start on hard scale models. Moderate and intercept-only sigma starts are unaffected.

• drmTMB() now warns at fit time if the optimized objective is not finite (NaN/Inf), instead of storing a non-finite log-likelihood and returning a broken fit silently. TMB normally returns a finite objective, so this is a defensive guard (class drmTMB_nonfinite_objective_warning); when it fires, the fit never reached a usable optimum and its estimates and standard errors are meaningless.

• The log(sigma) overflow soft-clamp now guards every scale-bearing family, not just Gaussian. Student, skew-normal, lognormal, gamma, Tweedie, beta, zero-one-beta, beta-binomial, and the negative-binomial family (NB2, truncated, hurdle, zero-inflated), plus the Gaussian row-aggregation path, all route log(sigma) through the same use_logsigma_clamp-gated soft-clamp before exponentiation. The clamp is exactly the identity inside the band, so every in-band fit is bit-identical to the unclamped fit (verified per family); the only change is that a runaway scale becomes a finite, clamp-flagged fit instead of an overflow to NaN. Previously this guard existed only for the Gaussian likelihood, leaving the other families’ dispersion parameters unprotected. The same configurable band (drm_control(logsigma_clamp = )) and the clamp-active warning apply across families.

• drmTMB() now warns at fit time when the log(sigma) soft-clamp is active at the optimum – the fitted log(sigma) reached or passed the identity band, so the clamp bent the scale. Such a fit can converge artificially: the saturated tail is flat, so the gradient vanishes and the optimizer may report convergence even though the scale ran to the bound, leaving estimates and standard errors near the clamp unreliable. The warning (class drmTMB_clamp_active_warning) names the value reached and the band and recommends rescaling the response, widening the band with drm_control(logsigma_clamp = ), adding within-group replication, or a penalized/MAP fit. It covers every clamp-guarded scale family (the detector reads the main log_sigma/log_sigma1/log_sigma2 scales, ignoring the unclamped missing-predictor imputation scales) and is filtered by the simulation harness, which tracks scale state separately. This closes the case where a clamp-bound fit returned looking fine under a false convergence = 0.

• confint(method = "wald") now flags a Wald interval that sits at a variance-component or correlation boundary instead of presenting it as an ordinary interval. When a random-effect or structured standard deviation is within sd_boundary of zero, or a correlation is within rho_boundary of +/-1, the row’s conf.status becomes "wald_at_boundary" and a warning (class drmTMB_wald_boundary_warning) points to method = "profile", because the symmetric Wald interval undercovers under boundary (chi-square-mixture) inference. The interval is still returned – a boundary is a warning, not an auto-discard – and a residual or distributional scale near zero is regular and is not flagged. The thresholds are exposed as confint(..., sd_boundary = 1e-4, rho_boundary = 0.98), matching the check_drm() defaults.

• AIC() and BIC() now have drmTMB methods that warn when the criterion is not a valid comparison. A REML fit warns that its restricted likelihood is comparable only across models with identical fixed effects (never ML versus REML, or different mean structures); a penalized (MAP) fit warns that the criterion is not standard, because logLik() returns the unpenalized data log-likelihood and a penalized parameter does not contribute a full degree of freedom. Previously these calls dispatched to stats::AIC.default, which reads the log-likelihood value and ignores the estimator, silently returning a meaningless number. For a plain maximum-likelihood fit the value is unchanged and no warning is emitted; the warnings carry condition classes drmTMB_ic_reml_warning and drmTMB_ic_map_warning.

• drmTMB() now warns at fit time when the optimizer reports non-convergence, instead of returning a non-converged fit that looks fine. The warning surfaces the nlminb code and message (for example “false convergence (8)”), points to check_drm(), and suggests the robust optimizer preset; print() annotates a non-zero convergence code with “(not converged; see check_drm())” rather than showing a bare integer next to a clean-looking coefficient table. The warning carries condition class drmTMB_convergence_warning so callers that do their own convergence bookkeeping (the simulation harness, or user code) can muffle just this signal while still seeing every other warning. This does not change any estimate; it makes a weak fit visible at the point of fitting.

• drm_control() now exposes the Gaussian log(sigma) overflow guard as a configurable knob: logsigma_clamp = c(lo, hi) sets the identity-in-band soft-clamp band (default c(-12, 12), unchanged) and logsigma_clamp_margin its saturation margin (default 3, saturating to [-15, 15]). logsigma_clamp = NULL disables the guard. Widen the band for legitimately huge-variance unstandardized responses, or disable it to inspect the raw overflow on a near-degenerate per-group scale model. The default is bit-identical to the previous fixed guard; the band is a numerical guard only and does not change identifiability (#570).

• drmTMB() now accepts an optional penalty = drm_phylo_penalty(sd_u, sd_alpha, cor_sd) argument that switches a phylogenetic fit to a penalized / maximum-a-posteriori (MAP) estimator: a penalised-complexity prior on each phylogenetic SD and an optional mean-zero normal on the phylogenetic correlation. This regularizes a weakly-identified phylogenetic variance or correlation (for example a scale-side phylogenetic field at about one observation per tip, or a coupled location-scale correlation pinned near +/-1) so the fit returns a finite, positive-definite estimate instead of stalling at a boundary. Plain maximum likelihood stays the default and bit-identical when penalty = NULL; a penalized fit is labelled estimator = "MAP", logLik() returns the unpenalized data log-likelihood (the penalty is stored in fit$phylo_penalty), and check_drm() adds a penalized_map note. The penalty does not manufacture identifiability: cor_sd has no universal value and must be chosen by a prior-sensitivity sweep, likelihood-ratio tests or AIC across penalized fits are not standard, and within-group replication remains the clean route to a fully identified coupled model. A known-truth coupled-q4 recovery simulation backs this – the penalty rescues the model to a positive-definite fit (which also restores Wald, profile, and bootstrap intervals) and recovers strong correlations near the right prior while over-shrinking weak ones. The Improving convergence article gains a penalized/MAP section, and docs/design/174-controls-and-convergence.md records the control catalog, the generalization-via-controls principle, and the interval-method guidance (#570).

• drmTMB() now fits the first primary Bernoulli/binomial response slice with native TMB family = stats::binomial(link = "logit"). Supported responses are explicit 0/1 event indicators and cbind(successes, failures) count responses; the fixed-effect likelihood includes the binomial normalizing constant so logLik(), AIC, and BIC match stats::glm() on overlapping logit models. Use beta_binomial() when successes out of known trials need extra-binomial variation through sigma. Non-logit links, factor-response ordering, proportions plus weights, weights = trials, sigma, random effects, structured effects, bivariate or mixed responses, and non-phylogenetic engine = "julia" binomial fits remain unsupported for this first slice (#569).

• confint(method = "bootstrap") now attaches a lightweight "bootstrap.diagnostics" attribute to returned bootstrap interval tables, with one row per bootstrap refit and target. The diagnostics record refit convergence, target availability, finite draw use, seed/backend/worker provenance, and refit-control flags while leaving the visible interval table at the usual target grain. The Ayumi q4 developer harness writes the same ledger to bootstrap-diagnostics.csv; this is troubleshooting evidence for failed or partial bootstrap runs, not a coverage or 10k-workflow claim (#555).

• confint() now warns when a default Wald interval is requested for the skew-normal slant nu, recommending method = "profile" (or method = "bootstrap") instead. An ADEMP pilot found the Wald slant interval over-rejects near nu = 0 (24-40% false positives versus the nominal 5%) because the Azzalini information is near-singular at alpha = 0. A later fixed-effect guard grid kept tail-floor exposure and fixed-gradient diagnostics visible rather than treating finite likelihoods as interval evidence. The warning is scoped to skew-normal nu only; Wald intervals for other families, including Student-t nu (tail shape) and Tweedie nu (power), are unchanged.

• confint(method = "profile") now accepts profile_endpoint_max_eval for direct scalar endpoint profiles, giving long variance-component or correlation diagnostics an explicit endpoint-evaluation budget. When the budget is reached, confint() returns an endpoint profile_failed row with missing endpoints instead of silently falling through to a full-profile fallback. The Ayumi q4 developer harness can pass the same budget through DRMTMB_AYUMI_Q4_PROFILE_ENDPOINT_MAX_EVAL and now separates returned-fit status from convergence/Hessian inference status (#555).

• confint(method = "profile") now returns an explicit row-level conf.status = "profile_failed" with missing endpoints when a direct numeric profile target fails during endpoint or tmbprofile evaluation, instead of aborting the whole interval request or labelling a non-finite interval as a successful profile. Focused regression tests keep bivariate q=4 phylogenetic location-scale sigma SDs visible as direct native-TMB ML profile targets, check weak-Hessian profile status for those sigma targets, and leave q4 phylogenetic correlations derived and not profile-ready (#551).

• Phase 18 now has a standalone fixed-effect skew-normal artifact lane (skew_normal_fixed_effect). The new DGP, summariser, smoke runner, grid writer, manual Actions task, and focused tests save aggregate, replicate-level, manifest, failure-ledger, fixed-effect Wald interval, optional profile, optional parametric-bootstrap, interval-evidence, interval-diagnostic, and interval-failure artifacts for bf(y ~ x, sigma ~ z, nu ~ w), family = skew_normal(). The default grid uses moderate shape-recovery sample sizes (n = 720 and 1440) because stochastic skewness recovery is sample-size dependent; this is repeatable smoke/grid infrastructure, not a formal 500- or 1000-replicate operating-characteristic result.

• skew_normal() now fits the first univariate fixed-effect skew-normal location-scale-shape route with public mu = E[y], public sigma = SD[y], and residual slant nu on the identity scale. The TMB likelihood transforms internally to native skew-normal xi, omega, and alpha = nu; focused tests cover density normalization, native-density comparison, Gaussian normal limit, positive and negative skew recovery, predictor-dependent nu, Gaussian false-positive behaviour, simulation, fixed-effect interval visibility, and malformed-neighbour rejection. Random effects, sd(group), known sampling covariance, structured effects, bivariate skew-normal models, residual rho12, latent skew(id), and skew aliases remain planned.

• drmTMB() now forwards REML = TRUE through the experimental engine = "julia" bridge for one route-specific bivariate q = 4 Gaussian phylogenetic location-scale DRM.jl cell when the installed DRM.jl build supports that Patterson-Thompson REML diagnostic. The bridge article now shows the glmmTMB-style top-level REML = TRUE/FALSE switch, the labelled four-axis phylo() syntax, and the current missing-response boundary. This bridge evidence does not establish native-TMB q4 REML; the native route has separate recovery evidence. It also does not establish HSquared AI-REML, non-Gaussian REML, broad R-to-Julia bridge support, public optimizer controls, q4 interval reliability, or q4 interval coverage; weights, missing-predictor imputation, non-default control, most non-Gaussian families, unsupported phylogenetic neighbours, corpair() entries, simulation, and persistent Julia handles remain native-TMB or future bridge work (#544).

• drmTMB() now records estimator, REML, requested_REML, and effective_REML on experimental Julia-engine fits, so downstream diagnostics can distinguish requested REML from the estimator actually fitted. Unsupported Julia REML requests now warn with the exact unsupported cell, fall back to ML, and state that native engine = "tmb" is only an REML fallback for its documented univariate Gaussian REML slice rather than for every rejected Julia cell (#555).

• biv_gaussian() now fits the first ordinary q8 all-endpoint location-scale slope covariance slice: matching labelled (1 + x | p | id) terms in mu1, mu2, sigma1, and sigma2. The fitted block has eight endpoint SDs and 28 latent group-level correlations across response-specific location intercepts, location slopes, scale intercepts, and scale slopes. The SDs appear in sdpars$mu and sdpars$sigma; the correlations appear in corpars$re_cov, corpairs(level = "group", block = "p"), summary()$covariance, profile_targets(), and check_drm(). The Phase 18 registry now exposes opt-in biv_gaussian_q8_endpoint and biv_gaussian_q8_endpoint_recovery Actions tasks; the recovery lane reports bias, RMSE, MCSE, and explicit interval unavailability. A 2026-06-07 local two-cell audit ran 20 replicates per cell and kept q8 at hold_diagnostic: 38/40 manifests completed, model-convergence rates were 0.263 and 0.158, positive-Hessian rates were 0 in both cells, two replicates failed with non-positive leading minors, and no Wald intervals were usable. Q8 still has no coverage result, power claim, predictor-dependent corpair() regression, random rho12, structured q8 sibling, or non-Gaussian q8 route.

• biv_gaussian() now fits the first same-response location-scale slope covariance slice: matching labelled (0 + x | p | id) terms in mu1/sigma1 or mu2/sigma2. The location-slope SD appears in sdpars$mu, the scale-slope SD appears in sdpars$sigma, and the group-level cor(mu1:x,sigma1:x | p | id) or cor(mu2:x,sigma2:x | p | id) row appears in corpars$mu_sigma, corpairs(class = "mean-scale-slope"), summary()$parameters, profile_targets(), and check_drm(). Cross-response pairs, mismatched coefficients, and univariate labelled sigma slopes remain closed; the all-endpoint q8 route is a separate source-tested slice.

• The Phase 18 same-response bivariate Gaussian q=2 mu/sigma slope covariance lane now has smoke and multi-replicate recovery artifacts (biv_gaussian_mu_sigma_slope and biv_gaussian_mu_sigma_slope_recovery). The lane reuses the fitted matching (0 + x | p | id) terms in mu1 and sigma1, reports 12 estimands, emits bias, RMSE, empirical SE, MCSE, and fixed-effect Wald coverage tables, and keeps the two slope SDs plus the derived mu_sigma correlation out of Wald interval claims. A local 2026-06-06 formal audit ran 500 replicates in each of the two default recovery cells and produced 1,000 ok manifest rows, but convergence/positive-Hessian rates were 0.856 and 0.884 and all-replicate fixed-effect Wald coverage was 0.796-0.850. A follow-up hardening audit regenerated and robust-refit the 130 weak replicates; none were rescued, all retained false-convergence and pdHess = FALSE, and estimates were unchanged. Among interval-available converged fits, fixed-effect Wald coverage was 0.930-0.972, and endpoint profiles succeeded on two clean representative fits for rho12, both slope SDs, and cor(mu1:x,sigma1:x | p | id). This is diagnostic evidence and profile feasibility, not power-grid support.

• The ordinary NB2 mu random-effect surface now has a standalone recovery artifact lane (nbinom2_mu_re_recovery), parallel to the Poisson one: it runs the already-recovery-capable smoke summary at recovery-scale n_rep and emits isolated bias/RMSE/MCSE, Wald-coverage, and profile-coverage CSVs through an opt-in Actions task, as a ready_grid random_slopes registry row. (The truncated-NB2 mu random-intercept surface already had an equivalent standalone coverage-emitting lane through its existing truncated_nbinom2_mu_random_intercept task, so it needed no new writer.)

• The ordinary Poisson mu random-effect surface now has a standalone, dispatchable recovery artifact lane (poisson_mu_re_recovery) — the first non-Gaussian recovery artifact lane. The recovery contract (bias, RMSE, MCSE, Wald coverage for the fixed mean coefficients, and profile coverage for the random-effect SD) was already computed by the smoke summary; the new opt-in lane runs it at recovery-scale n_rep and emits isolated CSV artifacts instead of only riding the combined first-wave summary. It is a ready_grid random_slopes registry row.

• The Phase 18 bivariate Gaussian slope-only mu1/mu2 covariance lane now has a multi-replicate recovery companion (biv_gaussian_mu_slope_recovery). It reuses the smoke DGP, fit, and runner for the matching (0 + x | p | id) block and reports bias, RMSE, empirical SE, Monte Carlo standard error, and Wald interval coverage across its 10 estimands. Wald coverage is reported only for the fixed mu1/mu2 endpoints; the two slope random-effect SDs and the derived slope-slope correlation stay derived_interval_unavailable. The lane is a ready_grid random_slopes registry row with its own opt-in Actions task and grid writer.

• The Phase 18 bivariate Gaussian q=6 mu1/mu2 location covariance lane now has a multi-replicate recovery companion (biv_gaussian_q6_location_recovery). It reuses the smoke DGP, fit, and runner for the matching (1 + x + z | p | id) block in both location formulas and reports bias, RMSE, empirical SE, Monte Carlo standard error, and Wald interval coverage across its 30 estimands. Wald coverage is reported only for the fixed mu1/mu2 endpoints; the six location random-effect SDs and the fifteen derived location-location correlations stay derived_interval_unavailable. The lane is a ready_grid random_slopes registry row with its own opt-in Actions task and grid writer.

• The Phase 18 bivariate Gaussian q=2 residual-scale intercept covariance lane now has a multi-replicate recovery companion (biv_gaussian_q2_scale_recovery). It reuses the smoke DGP, fit, and runner, runs at recovery-scale n_rep, and reports bias, RMSE, empirical SE, Monte Carlo standard error, and Wald interval coverage. Wald coverage is reported only for the fixed mu1/mu2 endpoints that carry a standard error; the random-effect scale SDs and the derived scale-scale correlation stay derived_interval_unavailable. The lane is registered as a ready_grid correlation_blocks row with its own opt-in Actions task and grid writer, and its design sheet is docs/design/156-phase-18-bivariate-scale-q2-recovery-ademp.md.

• The Phase 18 bivariate Gaussian q=2 residual-scale slope covariance lane now fits matching sigma1 = ~ x + (0 + x | p | id) and sigma2 = ~ x + (0 + x | p | id) blocks under biv_gaussian(). The two scale-slope SDs appear in sdpars$sigma, the group-level cor(sigma1:x,sigma2:x | p | id) row appears in corpars$sigma, corpairs(class = "scale-scale"), summary()$parameters, profile_targets(), and check_drm(), and residual rho12 remains a separate row-level correlation. The biv_gaussian_q2_scale_slope and biv_gaussian_q2_scale_slope_recovery Phase 18 tasks report smoke, bias, RMSE, MCSE, and fixed-effect Wald coverage evidence, while q8 now has separate diagnostic smoke/recovery artifact tasks that do not promote coverage or power (#483).

• The Phase 18 bivariate Gaussian q=4 mu1/mu2 location covariance lane now has a multi-replicate recovery companion (biv_gaussian_q4_location_recovery). It reuses the smoke DGP, fit, and runner for the matching (1 + x | p | id) block in both location formulas, runs at recovery-scale n_rep, and reports bias, RMSE, empirical SE, Monte Carlo standard error, and Wald interval coverage. Wald coverage is reported only for the fixed mu1/mu2 endpoints; the four location random-effect SDs and the six derived location-location correlations stay derived_interval_unavailable. The lane is a ready_grid random_slopes registry row with its own opt-in Actions task and grid writer.

• The Phase 18 simulation programme now has a bivariate Gaussian q=2 residual-scale intercept covariance smoke lane (biv_gaussian_q2_scale). It fits the already-supported matching sigma1 = ~ 1 + (1 | p | id) and sigma2 = ~ 1 + (1 | p | id) block under biv_gaussian(), summarising the two direct scale SDs from sdpars$sigma and the derived scale-scale correlation from corpars$sigma while keeping residual rho12 a separate layer. The lane is registered as a correlation_blocks row, has its own Actions task and grid writer, and provides one fittable scale-covariance prerequisite for the q8 endpoint gate.

• biv_gaussian() now fits matching ordinary mu1/mu2 location covariance blocks beyond the slope-only route, including q=4 (1 + x | p | id) and q=6 (1 + x + z | p | id) blocks in both location formulas with Phase 18 smoke artifact routing. The location SDs appear in sdpars$mu as direct log_sd_re_cov profile targets, and the group-level correlations appear in corpars$re_cov, corpairs(), and summary(fit)$covariance as derived-unavailable interval rows; same-response q2 location-scale slope covariance and q8 all-endpoint covariance now have separate smoke/recovery lanes, while predictor-dependent slope corpair() regressions, broad q > 2 recovery, coverage, power, and non-Gaussian correlated slopes remain planned (#440, #446).

• drmTMB() now accepts REML = TRUE for the first univariate Gaussian mixed-model slices: dense ordinary mu fixed effects, ordinary mu random intercepts or slopes, diagonal or dense known sampling covariance through meta_V(V = V), intercept-only sigma, complete responses, and no row aggregation, structured effects, or direct random-effect scale formulae. The ordinary mixed-model path is checked against lme4::lmer(..., REML = TRUE); the known-V path matches manual full restricted Gaussian likelihoods and metafor REML estimates, with the expected fixed-design determinant shift in reported metafor log likelihoods.

• drmTMB() now retries optimizer-call errors from the default deterministic nlminb() budget with the existing "careful" and "robust" optimizer presets when no explicit optimizer controls were supplied. Successful retries warn and record the selected preset in fit$optimizer_used and all attempted presets in fit$optimizer_attempts, while nonzero convergence-code fits still return for diagnostic inspection rather than being silently rerun (#506).

• drmTMB() now treats namespace-qualified formula markers such as drmTMB::phylo(...) and drmTMB::meta_V(...) as their unqualified equivalents during formula parsing, fixing the cryptic length-3 condition error triggered by drmTMB::phylo() in bivariate formulas (#504).

• is_converged() now returns a compact no-rerun convergence flag for drmTMB fits, with include_hessian = TRUE available when downstream workflows need successful TMB::sdreport() output and pdHess = TRUE before using Wald-style uncertainty (#317).

• structured_effects() now returns a stable post-fit metadata table for fitted phylo(), spatial(), animal(), relmat(), and phylo_interaction() structured-effect markers, so downstream tools can inspect grouping variables, matrix attachments, matrix slot/source/role IDs, compact precision fingerprints, provider and observed levels, level-alignment policy, input scale, bridge-marshalling boundary, fitted blocks, endpoint sets, coefficient sets, and random-effect block names without grepping formula text (#335).

• drmTMB() now fits the first sigma-only Gaussian structured one-slope cells for phylo(1 + x | species, tree = tree), fixed-covariance spatial(1 + x | site, coords = coords), A-matrix animal(1 + x | id, A = A), and relmat(1 + x | id, K/Q = ...) on the residual-scale formula. The same tranche opens matched mu+sigma one-slope native point-fit/extractor cells for those four providers by tracking mu:(Intercept), mu:x, sigma:(Intercept), and sigma:x as separate endpoint members, with deterministic same-target fixture parity banked for the sigma-only and matched cells. The fitted scale-side structured SDs appear in sdpars$sigma, ranef(), profile_targets(), and log-sigma predictions; broad bridge support beyond these fixtures, interval reliability, coverage, REML, and AI-REML remain planned.

• tweedie() now fits the first univariate fixed-effect route for non-negative semicontinuous responses with exact zeros and positive continuous values. The supported syntax is bf(y ~ x, sigma ~ z, nu ~ 1) with log(mu), log(sigma), public sigma = sqrt(phi), nu = 1 + plogis(eta_nu), E[y] = mu, and Var(y) = sigma^2 * mu^nu. Random effects, predictor-dependent nu, structured effects, bivariate or mixed-response Tweedie models, zero-inflation aliases, and hurdle aliases remain planned.

• student(), lognormal(), Gamma(link = "log"), beta(), beta_binomial(), and truncated_nbinom2() now support ordinary unlabelled independent numeric mu random slopes such as (0 + x | id) beside their first ordinary mu random-intercept slices. Historical boundary note, superseded for the exact zero-one-beta ordinary zoi and coi q1 intercept/same-raw-symbol slope gates: other zero-one-beta random effects, correlated slopes, labelled covariance, broad structured effects, and unsupported shape/inflation neighbours remain planned.

• The non-Gaussian tutorial route now connects count, beta-binomial, beta, and zero-one beta examples through the getting-started article, model map, implemented source map, worked-example inventory, and pkgdown article navigation. Historical boundary note, superseded for the exact zero-one-beta ordinary zoi q1 intercept and same-symbol slope gates: broader atom random effects, structured bounded responses, known covariance, ordered beta, beta-binomial zero inflation, and mixed bounded-response models remain planned.

• The model-guide route now includes a model-selection article for AIC/BIC comparisons. It shows Gaussian versus Student-t, NB2 versus ZINB2, and constant versus predictor-dependent sigma candidates, and it reads a 200-replicate seeded Phase 18 article-support summary that keeps MCSEs, convergence, Hessian, and warning rates beside AIC/BIC target-selection rates. The table is documentation evidence, not a formal power or operating-characteristic grid.

• zero_one_beta() now fits the first fixed-effect route for continuous proportions on [0, 1] with structural exact 0 or 1 outcomes. The interior beta component uses mu and public scale sigma; zoi models the probability of an exact boundary outcome; coi models the probability that a boundary outcome is exactly 1; and fitted() returns the unconditional mean (1 - zoi) * mu + zoi * coi. Historical boundary note, superseded for the exact ordinary zoi and coi q1 intercept/same-raw-symbol slope gates: structured effects, other atom random effects, denominator syntax, known covariance, and bivariate bounded-response models remain planned or blocked.

• gr() is now deprecated as a public formula marker. Existing direct calls warn and remain no-op placeholders for compatibility, while new known-relatedness formulas should use relmat() or the biological structured-effect markers animal(), phylo(), and spatial().

• meta_known_V() is now deprecated as a formula marker. Existing direct calls and formulas warn, while fitted Gaussian known-covariance models continue to use the same additive known-V likelihood path; new code should use meta_V(V = V).

• Historical 0.5.0 first-slice record (superseded for current response masking by MR-T2–MR-T6 and vignette("missing-data")): miss_control() exposed the initial missing-data control surface. The default response = "drop" kept existing complete-case behaviour, while response = "include" retained missing-response rows for univariate Gaussian models and independent-observation bivariate Gaussian models with complete predictors. Univariate masked responses contributed zero Gaussian likelihood; bivariate partial-response rows contributed the appropriate marginal Gaussian likelihood, both-missing rows contributed zero response likelihood, and fit$missing_data stored original-row accounting and response-pattern counts. The first missing-predictor slices also supported one numeric univariate Gaussian location term such as mi(x) with impute = list(x = x ~ z), impute = list(x = x ~ z + (1 | group)), or an explicit intercept-only structured covariate model such as impute = list(x = x ~ z + relmat(1 | line, Q = Q)) and missing = miss_control(predictor = "model"), integrating missing x values by TMB’s Laplace approximation under a fixed-effect, one random-intercept, or one structured-intercept Gaussian predictor model. The first family-specific non-Gaussian predictor slices supported one binary mi(treatment) term with impute = list(treatment = impute_model(treatment ~ z, family = binomial())), one ordered categorical mi(score) term with impute = list(score = impute_model(score ~ z, family = cumulative_logit())), one unordered categorical mi(habitat) term with impute = list(habitat = impute_model(habitat ~ z, family = categorical())), one strict proportion mi(cover) term with impute = list(cover = impute_model(cover ~ z, family = beta())), one boundary-proportion mi(cover) term with impute = list(cover = impute_model(cover ~ z, family = zero_one_beta())), one denominator-aware proportion mi(cover) term with impute = list(cover = impute_model(success ~ z, family = beta_binomial(), trials = trials)), one Poisson count mi(abundance) term with impute = list(abundance = impute_model(abundance ~ z, family = poisson())), one negative-binomial count mi(abundance) term with impute = list(abundance = impute_model(abundance ~ z, family = nbinom2())), one zero-truncated negative-binomial count mi(abundance) term with impute = list(abundance = impute_model(abundance ~ z, family = truncated_nbinom2())), one lognormal positive continuous mi(biomass) term with impute = list(biomass = impute_model(biomass ~ z, family = lognormal())), one Gamma positive continuous mi(biomass) term with impute = list(biomass = impute_model(biomass ~ z, family = Gamma(link = "log"))), and one Tweedie semi-continuous mi(biomass) term with impute = list(biomass = impute_model(biomass ~ z, family = tweedie())). Finite-state predictors sum exactly over possible missing states; strict beta/proportion predictors use deterministic quadrature over possible missing proportion values; zero-one beta boundary proportions use exact zero and one mass plus deterministic interior beta quadrature; beta-binomial denominator-aware proportions use deterministic finite summation over possible success counts; Poisson, negative-binomial, and zero-truncated negative-binomial count predictors use deterministic finite summation over count states; lognormal positive predictors use deterministic quadrature over log-scale states; Gamma positive predictors use deterministic quadrature under the Gamma mean-CV predictor model; Tweedie semi-continuous predictors use exact zero mass plus deterministic positive-support quadrature with fixed predictor-model power 1.5. imputed() reports fitted conditional modes for Gaussian missing predictors, fitted conditional probabilities for binary missing predictors, fitted conditional expected scores plus level probabilities for ordered missing predictors, fitted conditional modal categories plus level probabilities for unordered missing predictors, fitted conditional means for strict beta/proportion, boundary-proportion, beta-binomial, lognormal, Gamma, and Tweedie predictors, and fitted conditional expected counts for count predictors. Dense known-V partial-response slicing, structured covariate slopes, automatic response-structure inheritance, joint response-covariate structured correlations, multiple missing predictors, grouped or structured non-Gaussian predictor models, EM/profile engines, simulation-based imputed summaries, and measurement-error models remain planned.

• miss_control(predictor = "model") now has a first non-Gaussian response route: ordinary family = poisson() models can include one fixed-effect binary mi(treatment) predictor with impute = list(treatment = impute_model(treatment ~ z, family = binomial())). The route sums over the two missing treatment states with the Poisson response likelihood; missing Poisson responses, zero-inflated Poisson response models, Poisson response random or structured effects with mi(), non-binary missing predictors in Poisson response models, and multiple missing predictors remain planned.

• phylo_interaction() now fits the first q=1 pair-level phylogenetic interaction slice for univariate Gaussian mu and ordinary Poisson/NB2 mu models, using a sparse Kronecker precision from the two partner phylogenies. Use a precomputed pair column with ordinary (1 | pair_id) for independent pair effects, and keep additive partner main phylogenies, binary/Bernoulli incidence models, structured pair slopes, labelled count covariance, and simultaneous structured layers as planned follow-up work (#447).

• The status and learning docs now consistently lead with meta_V(V = V) for known sampling covariance, keep meta_known_V(V = V) as a deprecated compatibility alias, and refresh known-limitations wording so constant spatial, animal-model, and relmat() q=4 routes plus ordinary Poisson/NB2 q=1 structured mu routes are described as bounded first slices rather than broad planned or broad fitted support.

• Phase 18 now exposes a manual-only correlation_block_status Actions task that writes read-only CSV status artifacts for residual rho12, ordinary and structured q=2 corpairs() rows, and q=4 diagnostic rows. It removes the correlation-block wrapper-target gap without running new models, promoting q=4 intervals, or changing the fitted support boundary (#446).

• Phase 18 now exposes manual-only phylo_mu_slope, spatial_mu_slope, animal_mu_slope, and relmat_mu_slope Actions tasks for the Gaussian structured mu one-slope grid writers. Historical note, superseded by current 0.6.0 evidence: q1 structured sigma one-slope routes now fit for all four providers, and exact non-Gaussian provider gates are recorded in the live ledger. The manual Gaussian tasks remain excluded from task = "all"; mesh/SPDE, sparse large-pedigree speed claims, additional multiple or labelled structured-slope layouts outside the exact fitted ledger cells, slope correlations, and non-Gaussian structured slopes outside the exact later gates stay out of their scope.

• Phase 18 now has a local phylogenetic Gaussian mu one-slope artifact writer. Historical note, superseded by current 0.6.0 evidence: the exact q1 phylogenetic sigma one-slope route and exact non-Gaussian phylogenetic gates are now fitted at their recorded tiers. The manual Gaussian task remains excluded from task = "all"; additional multiple or labelled phylogenetic-slope layouts outside the exact fitted ledger cells, slope correlations, and non-Gaussian phylogenetic effects outside the exact later gates stay out of its scope.

• Phase 18 now has a local dense-pedigree animal() Gaussian mu one-slope artifact writer. The DGP, smoke runner, summary helper, grid writer, manual animal_mu_slope task, and focused tests save aggregate, replicate-level, manifest, and failure-ledger artifacts for animal(1 + x | id, pedigree = pedigree). Historical note, superseded by current 0.6.0 evidence: the exact A-matrix q1 sigma one-slope route is now fitted and inference-ready with caveats; pedigree/Ainv bridge marshalling, task = "all" inclusion, sparse large-pedigree speed claims, additional multiple or labelled animal-slope layouts outside the exact fitted ledger cells, and slope correlations remain out of scope.

• Phase 18 now has a local known-matrix relmat() Gaussian mu one-slope artifact writer. The DGP, smoke runner, summary helper, grid writer, manual relmat_mu_slope task, and focused tests save aggregate, replicate-level, manifest, and failure-ledger artifacts for relmat(1 + x | id, Q = Q). Historical note, superseded by current 0.6.0 evidence: the exact K/Q q1 sigma one-slope route is now fitted and inference-ready with caveats; broader bridge claims, task = "all" inclusion, additional multiple or labelled relmat()-slope layouts outside the exact fitted ledger cells, and slope correlations remain out of scope.

• drmTMB() now exposes the fitted TMB::sdreport() object as both $sdr and $sdreport, making Hessian checks easier to discover. Wald standard errors, vcov(), and Wald confidence intervals are unavailable when TMB::sdreport() returns but pdHess = FALSE, so non-positive-definite Hessian fits keep point estimates but no longer advertise Hessian-based intervals.

• The bivariate coscale tutorial and correlation-pair design note now sharpen the reader boundary between residual rho12, singular corpair() formula markers, and plural corpairs() extraction rows, including fitted structured extraction rows while keeping random effects in rho12 and unsupported corpair() regressions planned (#443).

• The ordinary Gaussian random-slope closeout now links the q=3 recovery, q=4 output-contract, extractor, corpairs(), summary(), profile_targets(), and independent log-sigma slope evidence before larger Phase 18 power simulations use those rows (#439).

• The random-slope support matrix now consistently separates fitted ordinary, bivariate slope-only, structured one-slope, selected non-Gaussian mu, and structured count q=1 routes from their planned neighbours (#438). Historical note, superseded by current 0.6.0 evidence: q1 structured sigma one-slope routes now fit for phylo/spatial/animal/relmat; phylo, A-matrix animal, and K/Q relmat are inference-ready with caveats, while spatial intervals remain blocked. Additional multiple or labelled structured-slope layouts outside the exact fitted ledger cells, slope correlations, and broader non-Gaussian neighbours remain planned.

• confint(method = "profile") now has profile_engine = c("auto", "endpoint", "tmbprofile"). The default auto route uses a faster endpoint-only scalar solver for direct constant scale, SD, and correlation targets, with curvature-seeded endpoint brackets and lower/upper endpoint splitting when a single endpoint target is profiled with Unix parallel = "multicore". Fixed-effect profiles, newdata profiles, linear combinations, and derived targets remain on the existing TMB::tmbprofile() or status-only paths. Profile rows now record the engine in profile.engine, parallel = "multicore" uses about half the detected CPU cores when workers = NULL, and bench/profile-scalar-endpoint.R records endpoint-versus-tmbprofile timing evidence for the phylogenetic SD target.

• confint() now defaults to fast Wald intervals for fixed effects plus direct fitted scale, random-effect SD, random-effect correlation, and constant rho12 targets when TMB::sdreport() is available. SD intervals use the fitted log-SD scale, correlation intervals use a guarded Fisher-z/atanh scale, profile_precision = "fast" supplies quicker profile controls for selected long-running targets, and method = "bootstrap" adds bounded simulate/refit percentile intervals with refit success/failure counts; positive scale and SD bootstrap intervals use fitted log-scale percentiles before exponentiating endpoints (reported from itchyshin/bergmann-drmTMB#2).

• confint(method = "bootstrap") now gives a direct-target-only error when a user requests a derived target such as a q4 unstructured correlation, modelled sd(group) surface, repeatability, or phylogenetic signal by exact name or by a broad alias such as parm = "correlations" or parm = "variance_components", instead of reporting the target as unknown or silently dropping unsupported rows.

• confint(method = "profile") now forwards parallel and workers to the target loop and accepts profile_maxit as an explicit per-target TMB::tmbprofile() budget guard. Serial remains the default, Unix multicore can split independent profile targets, and user-supplied maxit in ... is rejected when profile_maxit is also supplied.

• profile() now returns full profile-likelihood curve data for selected direct profile_targets() rows, and plot() draws the likelihood-ratio curve with the fitted estimate, likelihood-ratio cutoff, and profile confidence endpoints. The model-workflow article shows a 95% residual-sigma profile and the focused test suite now checks that the sampled curve extends beyond the cutoff on both sides of the interval.

• profile_targets() now reports fitted q=4 Julia bivariate phylogenetic SD estimates from the stored phylocov covariance instead of placeholder 0.5 values. The R-side Julia bridge target inventory now matches the fitted among-axis Sigma_a block before profile or bootstrap intervals are requested (#555).

• Univariate Gaussian models now fit residual-scale structured random intercepts with sigma ~ phylo(...), sigma ~ spatial(...), sigma ~ animal(...), and sigma ~ relmat(...). Matching intercept-only structured terms in mu and sigma estimate one latent structured mu-sigma correlation and report separate SD rows under sdpars$mu and sdpars$sigma, correlation rows under corpars, corpairs(), and profile_targets(). Historical note, superseded by the current 0.6.0 guidance: exact q1 sigma one-slope routes are now fitted for all four providers, with phylo, A-matrix animal, and K/Q relmat inference-ready with caveats; spatial sigma-slope intervals remain blocked. Additional multiple or labelled structured-sigma layouts outside the exact fitted ledger cells, direct-SD formulas combined with structured sigma, mesh/SPDE, and non-Gaussian residual-scale structured effects beyond the exact NB2 q1 recovery-grade routes remain planned.

• Bivariate sigma() output is now a roundable list for biv_gaussian() fits, so round(sigma(fit), digits) preserves the documented $sigma1 and $sigma2 components.

• The large phylogenetic benchmark runner can now include a real cell-level random intercept through --cell-random-effect true, giving speed checks a mixed-model path that is closer to applied repeated-cell data.

• plot_corpairs() now uses Confidence Eye regions by default for supported finite correlation intervals: the eye is a pale Fisher-z/atanh confidence region and the hollow circle is the point estimate. Conventional CI lines remain available with interval_style = "line" for diagnostic or reader-preference displays.

• The figure gallery now treats the default Confidence Eye as a stricter visual contract for selected row-wise interval summaries: pale finite confidence region plus hollow point-estimate circle, a dotted zero reference where zero is meaningful, a bottom scale axis, and no filled points, outlines, in-plot titles, or CI bars in default examples. Other figure classes keep purpose-specific grammar: raw-data displays, model surfaces, point summaries, simulation summaries, and support-boundary strips are judged case by case. Variance-component SD rows use log-SD Wald eyes, correlation rows use Fisher-z/atanh eyes, compact point summaries use point-interval displays, and the gallery surface example now labels the shared sigma ~ temperature curve explicitly rather than implying a habitat-specific scale effect.

• poisson() now fits the labelled-scalar spatial count route bf(count ~ x + spatial(1 | p | site, coords = coords)). The label is treated as a scalar covariance-block tag for the existing q1 spatial mu field and is exposed through sdpars$mu, ranef("spatial_mu"), and a direct profile_targets() row. This is local fit-only evidence only: q2/q4 count covariance, labelled slopes, simultaneous structured providers, intervals, coverage, inference_ready, supported, REML, AI-REML, bridge support, and public support remain out of scope.

• nbinom2() now fits the ordinary, non-zero-inflated q=1 phylogenetic mu intercept with syntax bf(count ~ x + phylo(1 | species, tree = tree), sigma ~ z). The fitted effect is on the log-mean scale while sigma remains fixed-effect overdispersion; sdpars$mu, ranef("phylo_mu"), profile_targets() as a direct log_sd_phylo target, and check_drm() phylogenetic diagnostics expose the route. Historical note, superseded by current recovery evidence: exact q1 structured sigma intercept-plus-one-slope routes now fit for phylo/spatial/animal/relmat. Labelled q2/q4 count blocks, zero-inflated NB2 phylogeny, richer structured sigma blocks, structured-sigma intervals/coverage, simultaneous structured count types, and count cross-parameter covariance remain planned.

• poisson() and nbinom2() now fit q=1 spatial(), animal(), and relmat() mu intercepts for ordinary non-zero-inflated count models, extending the existing q=1 phylo() count route. The fitted log-mean structured SD appears in sdpars$mu, marker-specific ranef() blocks such as ranef("spatial_mu"), direct profile_targets() rows through log_sd_phylo, and check_drm() structured diagnostics. Historical note, superseded by current recovery evidence: exact q1 structured NB2 sigma intercept-plus-one-slope routes now fit for phylo/spatial/animal/relmat. Pure or multiple structured count slopes, labelled q2/q4 count covariance beyond the exact Poisson scalar-label gate, zero-inflated structured effects beyond the exact local-fit gates, simultaneous structured count types, richer structured sigma blocks, and structured-sigma intervals/coverage remain planned.

• Phase 18 now has an opt-in count structured q1 artifact lane for ordinary Poisson/NB2 spatial(), animal(), and relmat() mu intercepts. The new DGP, summariser, smoke runner, summary helper, grid writer, manual count_structured_q1 Actions task, and focused tests save aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, Wald coverage, direct log_sd_phylo profile-target, optional profile-interval, interval-evidence, interval-diagnostic, and interval-failure artifacts without adding zero-inflated structure, structured slopes, labelled count covariance, structured NB2 sigma, task = "all" inclusion, or formal recovery claims.

• nbinom2() now fits the first ordinary log-sigma random-intercept gate for non-zero-inflated models, with syntax such as bf(count ~ x, sigma ~ z + (1 | id)). The fitted effect models grouped overdispersion on the log-sigma scale and is exposed through sdpars$sigma, random_effects$sigma, sigma(), predict(dpar = "sigma"), direct log_sd_sigma profile targets, and check_drm() replication diagnostics. Historical note, superseded by current recovery evidence: exact q1 structured NB2 sigma intercept-plus-one-slope routes now fit for phylo/spatial/animal/relmat. Ordinary NB2 sigma slopes, labelled covariance blocks, joint mu/sigma random effects, zero-inflated/truncated/hurdle NB2 scale random effects, richer structured sigma blocks, structured-sigma intervals/coverage, and Poisson scale random effects remain planned or inapplicable.

• beta() and beta_binomial() now support ordinary unlabelled mu random intercepts and independent numeric slopes such as bf(prop ~ x + (1 | id) + (0 + x | id), sigma ~ z) for strict (0, 1) responses and bf(cbind(success, failure) ~ x + (1 | id) + (0 + x | id), sigma ~ z) for counted successes out of known trials. The fitted logit-mean or logit-success-probability SD appears in sdpars$mu, random_effects$mu, direct profile_targets() rows, and check_drm() replication diagnostics; correlated bounded-response random slopes, labelled covariance blocks, sigma random effects, exact 0/1 boundary mass, zoi/coi, structured effects, known covariance, and bivariate or mixed bounded-response models remain planned.

• student() now supports ordinary unlabelled mu random intercepts and independent numeric slopes such as bf(y ~ x + (1 | id) + (0 + x | id), sigma ~ z, nu ~ 1). The fitted location SD appears in sdpars$mu, random_effects$mu, direct profile_targets() rows, and check_drm() replication diagnostics. Historical note: correlated Student-t random slopes, labelled covariance blocks, sigma random effects, nu random effects beyond the exact phylo local-fit gate, broad structured effects, and known covariance remain planned; the later Arc 6.4 exact biv_student() source slice is deliberately narrower and carries no recovery claim.

• lognormal() and Gamma(link = "log") now support ordinary unlabelled mu random intercepts and independent numeric slopes such as bf(y ~ x + (1 | id) + (0 + x | id), sigma ~ z). The fitted SDs appear in sdpars$mu, random_effects$mu, direct profile_targets() rows, and check_drm() replication diagnostics. This historical ordinary-effect entry is superseded in part by the exact Arc 3a q1 Gamma-phylo and lognormal-phylo/relmat intercept gates; correlated positive-continuous random slopes, labelled covariance blocks, sigma slopes, labelled or combined sigma random effects, other structured positive-continuous effects, known covariance, and bivariate or mixed positive-continuous models remain planned.

• Phase 18 now has a Student-t mu random-intercept artifact lane for student(). The new DGP, summariser, smoke runner, grid writer, first-wave runner inclusion, manual student_mu_random_intercept Actions task, and focused tests save aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, Wald coverage, direct-SD profile interval, and profile coverage artifacts for ordinary (1 | id) in mu with fixed-effect sigma and nu, while keeping correlated Student-t random slopes, labelled covariance blocks, sigma random effects, nu random effects beyond the exact phylo local-fit gate, broad structured effects, known covariance, and bivariate Student-t models out of scope.

• Phase 18 now has a zero-truncated NB2 mu random-intercept artifact lane for truncated_nbinom2(). The new DGP, summariser, smoke runner, grid writer, first-wave runner inclusion, manual truncated_nbinom2_mu_random_intercept Actions task, and focused tests save aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, Wald coverage, direct-SD profile interval, and profile coverage artifacts for ordinary (1 | id) in mu, while keeping correlated zero-truncated NB2 random slopes, labelled covariance blocks, sigma random effects, hurdle random effects, zero-inflated zero-truncated models, structured effects, and bivariate count models out of scope.

• Phase 18 now has a bounded-response mu random-intercept artifact lane for beta() and beta_binomial(). The new DGP, summariser, smoke runner, grid writer, first-wave runner inclusion, manual bounded_response_mu_random_intercept Actions task, and focused tests save aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, Wald coverage, direct-SD profile interval, and profile coverage artifacts for ordinary (1 | id) in mu. Historical boundary note, superseded for the exact zero-one-beta ordinary zoi q1 intercept and same-raw-symbol slope gates: other zero-one beta random effects, correlated bounded-response slopes, labelled covariance blocks, structured effects, known covariance, and mixed bounded-response models remain out of scope.

• Phase 18 now has a fixed-effect proportion artifact lane for beta() and beta_binomial(). The new DGP, summariser, smoke runner, grid writer, first-wave runner inclusion, manual proportion_fixed_effect Actions task, and focused tests save aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, and Wald coverage artifacts while keeping exact 0/1 boundary mass, zoi/coi, correlated bounded-response random slopes, labelled covariance blocks, sigma random effects, structured bounded responses, and mixed-response bounded models out of scope.

• Phase 18 now has a fixed-effect positive-continuous artifact lane for lognormal() and Gamma(link = "log"). The new DGP, summariser, smoke runner, grid writer, first-wave runner inclusion, manual positive_continuous_fixed_effect Actions task, and focused tests save aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, and Wald coverage artifacts while keeping Tweedie, generalized Gamma, positive-response random effects beyond the ordinary mu intercept slice, known-covariance positive responses, structured positive responses, and mixed-response positive-continuous models out of scope.

• Phase 18 now has a positive-continuous mu random-intercept artifact lane for lognormal() and Gamma(link = "log"). The new DGP, summariser, smoke runner, grid writer, first-wave runner inclusion, manual positive_continuous_mu_random_intercept Actions task, and focused tests save aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, Wald coverage, direct-SD profile interval, and profile coverage artifacts for ordinary (1 | id) in mu, while keeping correlated positive-continuous random slopes, labelled covariance blocks, sigma random effects, Tweedie, generalized Gamma, structured effects, known covariance, and mixed positive-continuous models out of scope.

• Phase 18 added a fixed-effect ordinal artifact lane for cumulative_logit(). Historical note, superseded in part by later ordinary mu intercept/slope recovery and an exact phylogenetic intercept gate: this artifact lane itself covers only aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, Wald coverage, cutpoint, and cutpoint-ordering outputs; correlated/labelled or other structured ordinal effects, ordinal scale/discrimination formulas, bivariate ordinal models, and mixed-response ordinal models remain out of scope.

• Phase 18 now has a fixed-effect zero-one beta artifact lane for zero_one_beta(). The new DGP, summariser, smoke runner, grid writer, first-wave runner inclusion, manual zero_one_beta_fixed_effect Actions task, and focused tests save aggregate, replicate, manifest, failure-ledger, fixed-effect Wald interval, and Wald coverage artifacts. Historical boundary note, superseded for the exact ordinary zoi q1 intercept and same-raw-symbol slope gates: other zero-one random effects, covariance blocks, denominator syntax, known covariance, structured bounded responses, and bivariate or mixed bounded-response models remain out of scope.

• Phase 18 Slice C now closes the count first-wave review lane as a documented evidence inventory rather than a new syntax lane. The new count-closure note ties together paired Poisson/NB2 mu random effects, NB2 log-sigma random intercepts, Poisson q1 phylo, NB2 q1 phylo, the NB2 hold_smoke_only formal gate, and the next Slice D choices while keeping COM-Poisson, generalized Poisson, Tweedie, zero-one beta, skew-normal, and new random-effect syntax out of scope.

• Phase 18 now has an overdispersion-aware NB2 phylogenetic q1 formal-admission lane for bf(count ~ x + phylo(1 | species, tree = tree), sigma ~ z). The new ADEMP sheet, DGP, target-plus-grouped-comparator fitter, summariser, smoke runner, summary helper, grid writer, formal-grid spec/read-back QA helpers, promotion-decision helper, focused tests, and manual nbinom2_phylo_q1_formal Actions task save aggregate, replicate, manifest, failure-ledger, Wald interval, Wald coverage, direct log_sd_phylo profile-target, optional profile-interval, interval-evidence, interval-diagnostic, and interval-failure artifacts. This does not create formal recovery or coverage claims until the 500-replicate grid is run and audited.

• Phase 18 Slices 541-555 now record a local NB2 phylogenetic q1 formal-audit pass. The all-cell sentinel ran 288 formal cells once and the representative audit ran 24 formal-shaped cells with five replicates each, both with direct log_sd_phylo profiles and grouped-comparator rows. The artifacts passed read-back QA, all sentinel fits and all replicate-audit rows converged, and the promotion helper correctly keeps the route at hold_smoke_only because the 500-replicate formal recovery gate was not run. Profile failures at true sd_phylo = 0 and fixed-sigma instability in low-count, low-overdispersion cells remain visible audit boundaries.

• Phase 18 formal phylogenetic q1 Actions tasks now accept one-based condition_shard and condition_shards inputs. The NB2 q1 full-grid singleton dispatch was cancelled after existing manifest timings implied a 27-31 hour run under optimistic ten-worker assumptions, so shard artifacts now record shard metadata and require a merged audit before any coverage claim.

• Phase 18 formal phylogenetic q1 Actions tasks now use shard-aware concurrency groups, so manually dispatched 16-shard NB2 q1 formal runs do not replace earlier pending shards. The supported non-Gaussian evidence goal is now recorded as an evidence closeout for fixed-effect non-Gaussian families plus first count mixed-model lanes, not as broad non-Gaussian random-effect or structured-effect parity.

• Phase 18 NB2 q1 formal shard artifacts have now been downloaded and audited together. The completed 16-shard x 500-replicate set has all 288 formal condition cells and 144,000 ok manifest rows, but the route remains hold_smoke_only because direct log_sd_phylo profile intervals are boundary-sensitive and low-count fixed-sigma recovery remains unstable.

• Phase 18 Slice D3 records the zero-one bounded-response design gate. The note separates strict beta(), denominator-aware beta_binomial(), and zero-one beta responses with exact 0/1 mass while keeping zero-one random effects, correlated or broader bounded-response random slopes, structured bounded responses, Tweedie, skew-normal, COM-Poisson, and generalized Poisson out of the fitted surface.

• poisson(link = "log") now fits the first structured non-Gaussian dependence slice: an ordinary, non-zero-inflated q=1 phylogenetic mu intercept with syntax bf(count ~ x + phylo(1 | species, tree = tree)). The fitted effect is on the log-mean scale and is exposed through sdpars$mu, ranef("phylo_mu"), profile_targets() as a direct log_sd_phylo target, and check_drm() phylogenetic diagnostics. This is intentionally not broad count parity: pure or multiple structured count slopes, labelled q=2/q=4 count blocks, zero-inflated structured count effects, simultaneous structured count types, and count cross-parameter covariance remain planned.

• biv_gaussian() now fits constant coordinate-spatial q=4 location-scale blocks when the same labelled spatial(1 | p | site, coords = coords) term appears in mu1, mu2, sigma1, and sigma2. The fitted route reports four spatial endpoint SDs and six derived latent spatial correlations through corpairs(level = "spatial"), summary()$covariance, profile_targets(), and check_drm(), with q=4 correlation intervals marked derived-unavailable. Historical note, superseded by the current 0.6.0 guidance: the q1 spatial sigma one-slope route now has point-fit/extractor evidence, although its interval gate remains blocked; mesh/SPDE, block-diagonal or broader intercept-plus-slope spatial covariance beyond the exact fixed-covariance all-four one-slope cell, additional multiple or labelled spatial-sigma layouts outside the exact fitted ledger cells, predictor-dependent spatial corpair() regression, direct spatial SD surfaces, and non-Gaussian spatial effects outside the exact ordinary Poisson/NB2 q1 spatial mu intercept-plus-one-slope, recovery-grade NB2 q1 spatial sigma, Student-t spatial mu, Poisson spatial zi, fixed-zi Poisson spatial mu, and fixed-zi NB2 spatial mu gates remain planned (#5).

• The implementation-map roadmap now records Slices 356-380 as the fitted spatial q4 closeout and Slices 381-388 as the first non-Gaussian structured front gate. The new Poisson q1 ADEMP sheet keeps the first structured non-Gaussian simulation gate to one non-zero-inflated Poisson mu phylogenetic intercept, while zi, hu, scale, shape, ordinal, bounded-response, mixed-response, spatial, animal, relmat(), slope, q2, and q4 structural layers remain planned rather than fitted.

• The implementation-map roadmap now records Slices 389-405 as the remaining non-Gaussian structured-dependence planning gates. These close the scale, shape, ordinal, known-covariance, extractor, diagnostic, simulation, interval, fallback, error-message, grammar, documentation, and issue-template contracts without adding new likelihood, TMB, or formula-grammar code.

• The implementation-map roadmap now records Slices 406-420 as route-specific non-Gaussian structured issue drafts. These add a Poisson q1 implementation issue body, Poisson q1 smoke-runner body, malformed-neighbour test body, documentation-sync body, NB2 q1 ADEMP skeleton, probability-component and scale/shape boundaries, and extractor/diagnostic name registries without opening new code.

• The implementation-map roadmap now records Slices 421-435 as the Poisson phylogenetic q1 runner contract. This names the direct log_sd_phylo target, sdpars$mu and ranef("phylo_mu") extractor checks, manifest and warning/error schemas, smoke and formal-grid gates, stale-doc corrections, malformed-neighbour error table, and focused test plan without adding new likelihood code.

• The implementation-map roadmap now records Slices 436-450 as Poisson phylogenetic q1 evidence-ledger synchronization. The source map, validation-debt register, Phase 18 programme, readiness matrix, and family registry now point to the runner contract before broad simulation claims.

• Phase 18 now has an opt-in Poisson phylogenetic q1 smoke surface for the ordinary non-zero-inflated phylo(1 | species, tree = tree) mu route. The new DGP, fitter, summariser, runner, summary helper, and focused tests return aggregate, replicate, manifest, failure-ledger, Wald fixed-effect interval, Wald coverage, and direct log_sd_phylo profile-target status tables; formal recovery grids remain future work.

• Phase 18 now has a repeatable grid-output writer for the Poisson phylogenetic q1 smoke surface, saving aggregate, replicate, manifest, failure-ledger, Wald interval, Wald coverage, and direct profile-target CSVs beside resumable per-replicate RDS files.

• Phase 18 Poisson phylogenetic q1 artifacts now include optional direct log_sd_phylo profile intervals, interval-evidence diagnostics, a formal-grid spec/read-back QA wrapper, a promotion-decision helper, and a manual poisson_phylo_q1_formal GitHub Actions task. The formal task is excluded from task = "all"; formal recovery or coverage claims still require the 500-replicate gate and artifact review.

• Phase 18 now has a separate ordinary NB2 log-sigma random-intercept smoke lane for bf(count ~ x, sigma ~ z + (1 | id)). The new DGP, summariser, runner, summary helper, and grid writer save aggregate, replicate, manifest, failure-ledger, Wald interval, Wald coverage, direct log_sd_sigma profile-target, optional profile-interval, interval-evidence, interval-diagnostics, and interval-failure artifacts beside resumable RDS results. Historical note, superseded by current 0.6.0 evidence: NB2 q1 structured sigma intercept-plus-one-slope routes for phylo/spatial/animal/relmat now have recovery-grade point-fit evidence. This is still not broad NB2 scale parity: ordinary NB2 sigma slopes, joint mu/sigma random effects, structured sigma intervals/coverage, richer structured blocks, and zero-inflated/truncated/hurdle scale random effects remain planned.

• The implementation-map roadmap now records Slices 526-540 as the NB2 phylogenetic q1 overdispersion-aware formal-admission lane. This adds the ADEMP sheet, DGP, target-plus-grouped-comparator fit, summariser, smoke runner, grid writer, formal-grid QA helpers, manual Actions task, tests, and docs sync while keeping formal recovery claims gated on a later 500-replicate run and audit.

• The implementation-map roadmap now records Slices 496-510 as the NB2 phylogenetic q1 implementation slice. This adds the ordinary non-zero-inflated NB2 mu phylogenetic intercept, TMB prior contribution, extractor/profile/diagnostic tests, and neighbouring-route guards. Historical note, superseded by current 0.6.0 evidence: NB2 q1 structured sigma intercept-plus-one-slope routes for phylo/spatial/animal/relmat now have recovery-grade point-fit evidence; their intervals/coverage, NB2 zi, richer structured slopes, and broader count covariance remain planned.

• The implementation-map roadmap now records Slices 341-355 as implementation-ready issue templates and acceptance gates for generic direct-SD syntax, p8/q8 slope covariance, spatial q4 parity, Poisson/NB2 q1 structured-count candidates, non-Gaussian structured ADEMP sheets, user documentation, review roles, and validation handoffs; these remain planning slices and do not add new fitted model surfaces.

• The implementation-map roadmap now records Slices 326-340 as pre-code specifications for generic direct-SD syntax, p8/q8 endpoint registries, spatial q4 parity, q4 diagnostics, Poisson/NB2 q1 structured-count candidates, user-route examples, and stale-claim checks; these remain planning slices and do not add new fitted model surfaces.

• The implementation-map roadmap now also records Slices 311-325 as design gates for generic structured direct-SD syntax, p8/q8 endpoint taxonomy, structured q=4 ordering, q=4 interval status, non-Gaussian structured-dependence candidate scoring, and user-route examples; these are planning slices, not new fitted likelihood claims.

• The implementation-map roadmap now records Slices 303-310 as planning and documentation gates rather than new likelihood work: generic sd*() design, p8/q8 location-scale planning, structured q=4 parity, q=4 interval policy, a historical no-fit decision for probability-component random effects, a non-Gaussian structured-dependence candidate map, maintenance scans, and user-route guidance. That no-fit decision is superseded only for the exact ordinary zero-one-beta zoi and coi q1 intercept/same-raw-symbol slope gates; broader zi, hu, zoi, and coi random effects remain outside the slice.

• The pkgdown site now includes an implementation map that separates fitted, first-slice, fixed-effect-only, planned, and blocked surfaces across families, random-effect layers, q, random slopes, corpairs(), zi, and hu, and uses it as the roadmap-facing ledger for future parity work.

• biv_gaussian() now fits the first ordinary bivariate random-slope covariance route: matching slope-only mu1/mu2 blocks such as (0 + x | p | id) in both location formulas. The fitted slope-slope row appears in sdpars$mu, corpars$mu, ranef(), corpairs(), summary()$covariance, profile_targets(), and check_drm(); matching q=4 and q=6 mu1/mu2 location blocks, same-response q2 location-scale slope covariance, and q8 all-endpoint covariance now have separate diagnostic artifact routes, while predictor-dependent slope corpair() regressions, q8 coverage/power evidence, and non-Gaussian structured slope covariance remain planned.

• The structural-dependence tutorial path now includes a focused coordinate-spatial page for fitted spatial(coords = coords) Gaussian mu intercept, residual-scale sigma intercept, one numeric mu slope, q=2 bivariate location-covariance, and constant q=4 location-scale routes. Historical note, superseded by the current 0.6.0 guidance: the q1 spatial sigma one-slope route now has point-fit/extractor evidence, although its interval gate remains blocked; mesh/SPDE inputs, additional multiple or labelled slope layouts outside the exact fitted ledger cells, direct spatial SD surfaces, spatial corpair() regressions, simultaneous phylo-plus-spatial layers, and non-Gaussian spatial effects outside the exact ordinary Poisson/NB2 q1 spatial mu intercept-plus-one-slope, recovery-grade NB2 q1 spatial sigma, Student-t spatial mu, Poisson spatial zi, fixed-zi Poisson spatial mu, and fixed-zi NB2 spatial mu gates remain planned.

• The structural-dependence tutorial path added a focused phylogenetic page for fitted Gaussian routes. Historical note, superseded in part by exact non-Gaussian gates: ordinary Poisson/NB2 q1 phylogenetic mu intercept-plus-one-slope, recovery-grade NB2 q1 phylogenetic sigma, Student-t q1 phylogenetic nu, and cumulative-logit q1 phylogenetic mu now fit at their recorded tiers; additional multiple or labelled slope layouts outside the exact fitted ledger cells, slope correlations, matrix-input phylogeny, combined phylo-plus-spatial layers, q4 corpair() regressions, and non-Gaussian phylogenetic neighbours outside those gates remain planned.

• The structural-dependence tutorial path now includes a focused relmat() page for fitted known-matrix Gaussian mu and sigma intercept slices, one-slope mu paths, matching q=2 bivariate location covariance, and constant q=4 location-scale covariance. Historical note, superseded by the current 0.6.0 guidance: the exact K/Q q1 sigma one-slope route is now fitted and inference-ready with caveats; broader bridge claims, additional multiple or labelled slope layouts outside the exact fitted ledger cells, slope correlations, predictor-dependent corpair() regression, and meta-analysis sampling covariance remain separate.

• The structural-dependence tutorial path now includes a focused animal-model page for the fitted animal(pedigree/A/Ainv) Gaussian mu and sigma intercept slices, one-slope mu paths, matching q=2 bivariate location covariance, and constant q=4 location-scale covariance. Historical note, superseded by the current 0.6.0 guidance: the exact A-matrix q1 sigma one-slope route is now fitted and inference-ready with caveats; pedigree/Ainv bridge marshalling, sparse pedigrees, additional multiple or labelled slope layouts outside the exact fitted ledger cells, slope correlations, and animal corpair() regression remain planned.

• The structural-dependence tutorial path now has a small overview article that helps readers choose between animal(), phylo(), coordinate spatial(), planned phylo-plus-spatial models, and relmat() before entering the longer technical tutorial.

• The figure gallery now shows spatial, animal, and relmat() q=2 fitted correlation rows beside residual rho12, ordinary group, and phylogenetic rows, and marks the constant spatial q=4 block as partly fitted beside still-planned structured correlation-regression and standalone scale extensions.

• animal() now fits a dense first pedigree route for Gaussian mu and sigma animal intercept effects: animal(1 | id, pedigree = pedigree) builds an additive relationship matrix from id, dam, and sire columns, one numeric animal(1 + x | id, pedigree = pedigree) mu slope fits as independent intercept and slope fields, and matching labelled animal(1 | p | id, pedigree = pedigree) terms work in the first bivariate q=2 location-covariance and constant q=4 location-scale paths. Historical note, superseded by the current 0.6.0 guidance: the exact A-matrix q1 sigma one-slope route is now fitted and inference-ready with caveats; pedigree/Ainv bridge marshalling, large-pedigree sparse precision construction, additional multiple or labelled structured-slope layouts outside the exact fitted ledger cells, slope correlations, predictor-dependent corpair() regressions, and generic direct-SD grammar remain planned (#147).

• biv_gaussian() now fits the first animal-model and lower-level relatedness q=2 known-matrix location covariance: matching animal(1 | p | id, A = A) / animal(1 | p | id, Ainv = Ainv) or relmat(1 | p | id, K = K) / relmat(1 | p | id, Q = Q) terms in mu1 and mu2. The fitted rows appear in sdpars$mu, corpars$animal or corpars$relmat, ranef("animal_mu") or ranef("relmat_mu"), corpairs(), summary()$covariance, profile_targets(), and check_drm(). Historical note, superseded by the current 0.6.0 guidance: univariate Gaussian sigma intercepts, one numeric mu slope, and the exact A-matrix animal and K/Q relmat q1 sigma one-slope routes are fitted; the sigma slopes are inference-ready with caveats. Pedigree/Ainv bridge marshalling, additional multiple or labelled structured-slope layouts outside the exact fitted ledger cells, slope correlations, predictor-dependent corpair() regressions, and generic direct-SD grammar remain planned (#147).

• biv_gaussian() now also fits constant all-four q=4 animal-model and relmat() location-scale blocks when the same labelled known-matrix term appears in mu1, mu2, sigma1, and sigma2. These rows reuse the structured covariance backend and report four endpoint SDs and six derived latent correlations through corpairs(), summary()$covariance, profile_targets(), and check_drm(). Historical note, superseded by the current 0.6.0 guidance: exact q1 A-matrix animal and K/Q relmat sigma one-slope routes are now fitted and inference-ready with caveats; additional multiple or labelled structured-slope layouts outside the exact fitted ledger cells, slope correlations, predictor-dependent corpair() regressions, and direct-SD grammar remain planned (#147).

• The pkgdown workflow now builds the advertised single-site URL rather than publishing only the dev/ subtree for development versions, so https://itchyshin.github.io/drmTMB/ remains the public entry point.

drmTMB 0.1.3 (2026-05-20)

• animal() and relmat() now fit the first known-relatedness Gaussian mu random-intercept slice: animal(1 | id, A = A), animal(1 | id, Ainv = Ainv), relmat(1 | id, K = K), and relmat(1 | id, Q = Q). The fitted latent scale appears in sdpars$mu, conditional effects appear in ranef("animal_mu") or ranef("relmat_mu"), direct scale targets appear in profile_targets(), and check_drm() reports replication and scale-ratio diagnostics. Historical note, superseded by current 0.6.0 evidence: one-slope mu, sigma intercept and exact q1 one-slope routes, and selected bivariate relatedness covariance are fitted where documented. Sparse large-pedigree construction, additional multiple or labelled structured-slope layouts outside the exact fitted ledger cells, broader covariance/bridge claims, and predictor-dependent corpair() regressions remain planned (#147).
• biv_gaussian() now fits the first coordinate-spatial q=2 mu1/mu2 location covariance through matching spatial(1 | p | site, coords = coords) terms. The fitted spatial SDs appear in sdpars$mu, conditional fields in ranef("spatial_mu"), the spatial mean-mean row in corpairs(level = "spatial") and summary()$covariance, and direct SD/correlation targets in profile_targets(); mesh/SPDE, spatial sigma, spatial q=4, direct spatial SD surfaces, and predictor-dependent spatial corpair() regression remain planned (#5).

drmTMB 0.1.2 (2026-05-16)

• drm_control() now has optimizer_preset = "careful" and "robust" for explicit nlminb() optimizer-budget presets. These expand to recorded iter.max and eval.max controls, keep the default fit fast, and can still be overridden with optimizer = list(...).
• drm_control() now reserves fallback-optimizer control names such as fallback_optimizer, fallback_optimizers, and optimizer_fallback while documenting the future selected-optimizer provenance contract. Fallback BFGS or L-BFGS-B refits remain planned, not automatic.
• drm_control() now reserves warm-start control names such as start_from, warm_start, and warm_start_from so simpler-fit starts cannot be silently passed to nlminb() before the source-fit contract is implemented.
• Phase 18 private parametric-bootstrap helpers now accept serial or Unix multicore execution, cap actual workers at 10, and record requested versus actual core counts in bootstrap draw and interval tables; PSOCK remains excluded until fitted TMB object rebuilds are explicit.
• Phase 18 replicate execution now has a private bounded runner helper, with the Gaussian location-scale, meta_V(V = V), Poisson and NB2 mu random-effect, Gaussian mu and sigma random-slope, coordinate spatial mu slope, Student-t shape, and bivariate residual rho12 smoke surfaces wired through serial or Unix multicore execution capped at 10 workers; closure-heavy runners use a per-replicate summary factory to preserve profile and bootstrap seeds. Higher-level grid and count-gallery wrappers now forward runner settings, and Student-t shape plus bivariate residual rho12 wrappers carry separate bootstrap backend settings with a guard against multicore replicate and multicore bootstrap layers running at the same time.
• The paired Poisson/NB2 mu random-effect Phase 18 lane now has a repeatable grid-output writer that saves aggregate, replicate, manifest, failure-ledger, Wald interval, Wald coverage, profile interval, and profile coverage CSV artifacts beside resumable per-replicate RDS files.
• The Phase 18 meta_V(V = V) lane now has a repeatable grid-output writer that saves aggregate, replicate, manifest, failure-ledger, Wald interval, and Wald coverage CSV artifacts beside resumable per-replicate RDS files.
• Phase 18 now also has repeatable simple grid-output writers for ordinary Gaussian mu random slopes, independent Gaussian sigma random slopes, and coordinate-spatial Gaussian mu slopes.
• Phase 18 grid writers now return an artifact manifest with path existence and CSV row counts, including zero-row handling for optional interval artifacts.
• Phase 18 artifact manifests can now be bound across grid-writer outputs and summarized by surface, giving report-staging code a compact check of present, missing, empty, and total CSV rows.
• Phase 18 first-wave report staging now has a private artifact-status writer that saves bound artifact-manifest and surface-status CSVs from multiple grid-writer outputs before a report reads the simulation tables.
• Phase 18 first-wave report staging now includes an artifact-status report template that reads the bound manifest and status CSVs first, renders a preflight page for complete outputs, and stops clearly when required artifacts are missing.
• Phase 18 first-wave report staging now has a private table-bundle writer that combines selected CSV artifacts across grid-writer outputs, preserving source surface and artifact names as leading columns while filling missing columns.
• Phase 18 first-wave report staging now includes a summary-report skeleton that reads artifact status, aggregate operating-characteristic rows, interval coverage, interval diagnostics, interval failures, manifests, and warning/error ledgers in one page, with priority columns, row caps, a compact aggregate-bias overview, compact interval-coverage summaries, run-manifest summaries, and compact warning/error summaries for table-first review.
• Phase 18 first-wave report staging now has a reusable private smoke runner that executes the Gaussian location-scale, meta_V(V = V), paired Poisson/NB2 mu random-effect, ordinary Gaussian mu random-slope, ordinary Gaussian sigma random-slope, and coordinate-spatial Gaussian mu slope grid writers, stages the combined first-wave summary report, and records requested versus actual worker counts.
• Phase 18 interval-heavy report staging now has a separate private smoke runner for Student-t shape and bivariate residual rho12 grid writers, keeping their Wald/profile/bootstrap interval artifacts separate from the baseline first-wave runner.
• Phase 18 first-wave report staging now has a render helper that writes artifact status, table bundles, and an optional HTML summary report from grid-writer outputs in one orchestration step.
• check_drm() now reports the largest fixed-gradient component label in the fixed_gradient diagnostic row, making non-converged or flat-surface fits easier to triage before Hessian or Wald inference is trusted.
• biv_gaussian() now fits multiple independent same-response labelled mu/sigma random-intercept covariance blocks in one two-response model, for example mu1/sigma1 with label p and mu2/sigma2 with label q plus residual rho12. corpairs(), profile_targets(), summary(), and check_drm() keep the two mean-scale rows separate from residual rho12 and from same-parameter mu1/mu2 or sigma1/sigma2 blocks.
• Interval documentation now separates fitted-model Wald and profile intervals from Phase 18 Fisher-z simulation helpers, and tests confirm Student-t nu fixed-effect interval targets plus Fisher-z-scale correlation-helper output.
• meta_V(V = V) now has explicit full-matrix alias coverage with Wald fixed-effect interval checks, documentation points to it as the preferred known-covariance spelling, and meta_V(V = V, scale = "exact") now errors with guidance that the exact additive route is already selected by meta_V(V = V).
• Bergmann-report follow-up now hardens invalid Wald standard-error rows, documents long-iteration triage, implements the labelled phylogenetic q4-to-two-q2 block-diagonal fallback, and records the Ayumi Mass + Beak fallback as a boundary/false-convergence diagnostic case. The earlier unsupported univariate sigma ~ phylo(...) boundary is now superseded by the fitted intercept-only structured sigma route described above.
• Structural-dependence docs now give a clearer user surface for planned animal() and relmat() models, including what fitted sensitivity model to use now and when a known matrix belongs to meta_V(V = V) instead of a future latent relatedness path.
• Structural-dependence design notes now separate dense covariance inputs (A, K) from sparse precision inputs (Ainv, Q) for future animal() and relmat() models, keeping large-pedigree or large-matrix speed claims blocked until sparse-precision recovery and scaling evidence exists.
• The family registry now has a Slice 283 family-and-parameter evidence map, listing each public family route, distributional-parameter link, shape or coscale slot, fitted random-effect allowance, and test evidence state before later count, proportion, shape, ordinal, or mixed-response hardening work expands those rows.
• Count-family tests now assert Wald fixed-effect interval rows for Poisson, NB2, zero-truncated NB2, zero-inflated Poisson, zero-inflated NB2, and hurdle NB2 dpars, and the count tutorial now names ordinary non-zero-inflated Poisson/NB2 mu random intercepts and independent slopes as the current fitted mixed-count route.
• Proportion-family tests now assert Wald fixed-effect interval rows for beta and beta-binomial mu and sigma coefficients, and the proportion tutorial now states that fixed-effect beta() and beta_binomial() plus ordinary unlabelled mu random intercepts and independent numeric slopes are fitted while zero-one inflation, correlated bounded-response random slopes, labelled covariance blocks, sigma random effects, and meta_V(V = V) bounded-response routes remain planned or blocked.
• Continuous-shape design notes now separate fitted fixed-effect Student-t nu, fitted fixed-effect skew-normal nu, planned skew-t nu/tau, and future latent-effect skew(id) ~ ..., keeping shape and skewness random effects out of Phase 18 simulation until likelihood, recovery, diagnostic, and interval evidence exists.
• Ordinal readiness docs recorded the original fixed-effect cumulative_logit() evidence ledger. Historical note, superseded in part by later ordinary mu intercept/slope recovery and an exact phylogenetic intercept gate: broader ordinal covariance and scale/discrimination formulas remain planned.
• Mixed-response bivariate family docs now keep Gaussian-count, Gaussian-proportion, count-proportion, ordinal mixed, and other two-response combinations planned until a joint likelihood or copula/latent-variable contract, prediction, simulation, extractors, intervals, examples, and comparator checks exist; tests cover mixed-family errors for c() and list() spellings.
• corpairs() now returns conf.status and interval_source columns by default, matching the prediction-table provenance contract. plot_corpairs() now draws finite bounds only when those columns mark a real interval source, so compatible pair tables cannot imply confidence intervals from bare numeric limits alone.
• User-facing docs now share a status vocabulary for stable, first slice, opt-in, planned or reserved, and unsupported or blocked surfaces across README, the model-map article, the package reference topic, getting-started article, source-map guidance, and pkgdown reference groups.
• The pre-simulation readiness matrix now includes a Slice 291 evidence-ledger gate: every public stable-core row is mapped to implementation evidence, tests or diagnostics, user-facing boundaries, and Phase 18 admission status before comprehensive simulation design can treat it as admitted.
• The Phase 18 simulation programme now has a Slice 292 comprehensive design map covering continuous, proportion, count, ordinal, meta-analysis, bivariate, random-slope, shape, phylogenetic, spatial, animal(), and relmat() lanes, while keeping planned or blocked lanes in the failure ledger instead of fitted grids.
• Phase 18 now has its first one-page ADEMP sheet for the admitted Gaussian location-scale lane, tying the existing phase18_dgp_gaussian_ls() helper to aims, DGP conditions, estimands, methods, performance measures, and Williams-style reporting checks before larger grids run.
• Phase 18 now has a one-page ADEMP sheet for the admitted Gaussian meta_V(V = V) lane, keeping known sampling covariance as input data and public residual sigma as the fitted heterogeneity estimand before vector or dense known-V grids expand.
• Phase 18 now has a one-page ADEMP sheet for the paired Poisson/NB2 mu random-effect lane, keeping the first count grid to ordinary non-zero-inflated mu random intercepts and independent numeric slopes while zero-inflated, hurdle, zero-truncated, structured, correlated-slope, and labelled covariance count models remain failure-ledger rows.
• Phase 18 now has a one-page ADEMP sheet for the fixed-effect proportion lane, separating strict continuous beta() responses from denominator-aware beta_binomial() success counts while keeping exact 0/1 boundary mass outside that earlier beta/beta-binomial artifact lane. Random effects beyond the beta and beta-binomial ordinary mu intercept/slope slices, structured effects, known sampling covariance, and mixed-response bounded models remain in the failure ledger.
• Phase 18 added a one-page ADEMP sheet for the fixed-effect ordinal lane. Historical note, superseded in part by later ordinary mu intercept/slope recovery and an exact phylogenetic intercept gate: correlated/labelled or other structured ordinal effects, scale/discrimination formulas, cutpoint-specific predictors, bivariate ordinal models, and mixed-response ordinal models remain in the failure ledger.
• Phase 18 now has a one-page ADEMP sheet for the bivariate Gaussian residual rho12 lane, defining the response-specific mean and scale DGP, residual covariance matrix, response-scale rho12 grids, and boundary diagnostics while keeping group-level corpairs(), structured correlations, known sampling covariance, random effects in rho12, mixed-response families, and bivariate random-slope covariance in separate design or failure-ledger lanes.
• Broader bivariate random-slope combination boundaries now have focused error coverage for unsupported residual-scale slope variants, cross-response or coefficient-mismatched same-response location-scale slope combinations, and all-four q=8-style slope requests while the matching q2 sigma1/sigma2 scale-slope route is fitted separately.
• Structured random-slope boundaries now have parser and fit-time audit coverage: one-slope animal() and relmat() markers are readable as planned grammar, and multiple structured slopes are rejected. Only coordinate spatial Gaussian mu one-slope models are fitted for Phase 18 admission.
• Shape and inflation random-effect boundaries now have random-slope-specific test coverage. Historical boundary note, superseded for the exact ordinary zero-one-beta zoi and coi q1 intercept/same-raw-symbol slope gates: Student-t nu, zero-inflation zi, hurdle hu, and all broader, transformed, mismatched, joint, or structured atom random effects still error before fitting.
• Gaussian sigma random-effect documentation now states the independent residual-scale slope boundary more explicitly: separate terms such as sigma ~ z + (0 + w_id | id) + (0 + w_site | site) are fitted with direct log_sd_sigma profile targets. Historical note, superseded by the current 0.6.0 guidance: unlabelled ordinary correlated blocks such as (1 + x | id) are now fitted; labelled residual-scale slope covariance remains planned.
• Ordinary Gaussian mu random-slope documentation now states the q > 2 boundary more explicitly: multi-slope blocks such as (1 + x1 + x2 + x3 | id) are fitted, their SDs are direct profile targets, and their block correlations are derived-unavailable for direct profile intervals until a dedicated interval method exists.
• Phase 18 now has a pre-simulation capability audit table that records implemented, tested, planned, and unsupported status for Gaussian, non-Gaussian, shape, inflation, bivariate, random-slope, meta-analysis, phylogenetic, spatial, animal, and relmat() model classes before broad simulation grids admit them.
• Phase 18 now has an optional resumable replicate runner under inst/sim/ that captures warnings, errors, elapsed time, session metadata, and optional RDS output for pilot simulation cells.
• Phase 18 now has a Simulation & Comparison plot-grammar article for bias, RMSE, coverage, power, convergence, runtime, and warning/error ledgers across continuous, proportion, count, and meta-analysis examples.
• plot_corpairs() now has a label argument so publication figures can use short row labels while keeping full correlation metadata in the source table.
• Phase 18 replicate results can now be reduced to a warning/error failure ledger so failed fits and warning-bearing replicates remain visible beside aggregate summaries.
• Phase 18 saved replicate directories can now be reloaded into result lists, allowing manifests and failure ledgers to be rebuilt from RDS output after a resumable simulation run.
• Phase 18 simulation results can now be reduced to a compact manifest table with cell id, replicate, seed, status, skipped/resumed flag, warning count, error message, and elapsed time.
• Phase 18 summary-smoke helpers now return run manifests and warning/error ledgers beside aggregate bias, RMSE, and MCSE tables.
• Developer design notes now include ASReml efficiency lessons for future animal() and relmat() work, emphasizing sparse precision matrices, explicit row-name matching, matrix-orientation metadata, and honest speed claims.
• nbinom2() now supports ordinary non-zero-inflated mu random intercepts and independent numeric random slopes such as bf(count ~ x + (1 | id) + (0 + x | id), sigma ~ z). The fitted SDs appear in sdpars$mu, random_effects$mu, and direct profile_targets() rows. Historical note, superseded by current 0.6.0 evidence: NB2 q1 structured sigma intercept-plus-one-slope routes for phylo/spatial/animal/relmat now have recovery-grade point-fit evidence, and one exact fixed-zi NB2 mu ~ spatial() intercept is diagnostic-only. Correlated or labelled NB2 mu slopes, joint mu/sigma random effects, ordinary NB2 sigma slopes, structured sigma intervals/coverage or richer blocks, and zero-inflated NB2 random effects outside that exact diagnostic gate remain planned.
• Phase 18 now includes a smoke simulation report template under inst/sim/reports/, giving pilot surfaces a reader-facing structure for purpose, aggregate summaries, reader checks, and interpretation boundaries.
• Phase 18 now has a Gaussian mu q=3 random-slope smoke surface under inst/sim/, covering seeded data generation, drmTMB() fitting for (1 + x1 + x2 | id), parameter summaries, aggregate output, manifests, and failure ledgers.
• Phase 18 now has a Gaussian sigma independent random-slope smoke surface under inst/sim/, covering seeded data generation, drmTMB() fitting for sigma ~ z + (0 + w | id), parameter summaries, aggregate output, manifests, and failure ledgers.
• Phase 18 now has a structured-slope parity gate: coordinate spatial one-slope Gaussian mu models are fitted enough for focused smoke grids, while phylogenetic, animal, and relmat() one-slope models remain planned until their implementation, diagnostics, profile targets, recovery tests, and biological examples exist.
• Phase 18 now has a cross-distributional-parameter correlation gate: residual rho12, constant fitted random-effect block correlations, predictor-dependent q=2 corpair() routes, and known sampling covariance V remain separate layers, while non-Gaussian, slope-level, shape, inflation, hurdle, one-inflation, and rho12 random-effect covariance surfaces stay outside Wave A until focused gates close.
• Phase 18 now has a coordinate spatial Gaussian mu one-slope smoke surface under inst/sim/, covering seeded data generation, drmTMB() fitting for spatial(1 + x | site, coords = coords), parameter summaries, aggregate output, manifests, and failure ledgers.
• Phase 18 now has a non-zero-inflated Poisson mu random-effect smoke surface under inst/sim/, covering seeded data generation, drmTMB() fitting for (1 | id) + (0 + x | id), parameter summaries, aggregate output, manifests, and failure ledgers.
• Phase 18 non-zero-inflated Poisson mu random-effect smoke output now includes Wald interval rows and coverage summaries for fixed log-mean coefficients, while random-effect SD rows remain visible as missing-SE interval rows until profile producers are attached.
• Phase 18 non-zero-inflated Poisson mu random-effect smoke output now includes direct profile-likelihood interval rows and coverage summaries for the fitted random-intercept and independent random-slope SD targets.
• Phase 18 now has a non-zero-inflated NB2 mu random-effect smoke surface under inst/sim/, covering seeded data generation, drmTMB() fitting for (1 | id) + (0 + x | id) with fixed-effect sigma ~ z overdispersion, parameter summaries, aggregate output, manifests, and failure ledgers.
• Phase 18 non-zero-inflated NB2 mu random-effect smoke output now includes Wald interval rows and coverage summaries for fixed log-mean and log-overdispersion coefficients, while random-effect SD rows remain visible as missing-SE interval rows until profile producers are attached.
• Phase 18 non-zero-inflated NB2 mu random-effect smoke output now includes direct profile-likelihood interval rows and coverage summaries for the fitted random-intercept and independent random-slope SD targets.
• NB2 mu random effects now have a focused weak-SD boundary diagnostic test, exercising check_drm() lower-boundary reporting for a near-zero fitted random-intercept SD before larger Phase 18 grids vary the true SD.
• Phase 18 Poisson and NB2 mu random-effect condition helpers now build true crossed condition grids, including random-effect SDs and, for NB2, overdispersion settings.
• Phase 18 now has a paired count-family mu random-effect pilot helper that combines ready Poisson and NB2 surfaces into one optional aggregate, manifest, failure-ledger, Wald-coverage, and profile-coverage output.
• Phase 18 now has a plot-data helper for paired Poisson/NB2 mu random-effect pilot outputs, preparing aggregate, coverage, manifest, and failure tables for later figure-gallery work.
• Phase 18 now has a count-pilot figure-gallery report template for bias, RMSE, interval coverage, manifests, and warning/error ledgers from paired Poisson/NB2 mu random-effect pilots.
• Phase 18 now has count-pilot gallery helper plumbing that writes plot-ready CSV inputs and renders a checked local HTML gallery artifact from a paired Poisson/NB2 mu random-effect pilot object.
• Phase 18 now has an end-to-end count-gallery smoke runner that executes a tiny paired Poisson/NB2 mu random-effect pilot and renders the Florence-facing HTML gallery from the resulting tables.
• The Phase 18 count-pilot gallery now has a first Florence visual-polish pass, with horizontal estimand labels, shared palette/theme helpers, captions, and MCSE-aware coverage ranges when available.
• The figure gallery now has a Florence visual-repair pass: inference summaries render as Confidence Eye displays where finite interval bounds are available, tutorial plots share explicit palettes, discrete comparison and empirical marginal displays no longer fall back to default black styling, status-strip labels have better contrast, and simulation coverage/power examples show replicate-block proportions plus aggregate binomial MCSE intervals instead of treating simulation uncertainty as required Confidence Eyes.
• The Simulation & Comparison plot-grammar article now carries the raincloud lesson into the Phase 18 simulation lane: bias displays show replicate-level errors plus mean/MCSE intervals in fixed surface facets, while RMSE keeps a separate aggregate point/MCSE panel instead of being visually collapsed with signed bias or mistaken for a mean absolute-error cloud.
• The Phase 18 count-pilot gallery now follows the same accuracy-display contract: bias and RMSE use fixed family facets, show MCSE intervals when available, and explain that replicate-error clouds require replicate-level output rather than aggregate CSVs alone.
• The figure gallery now keeps correlation layers visually separate, with faceted residual rho12, ordinary group, and phylogenetic corpairs()-style rows plus status-strip boundaries for structured-effect layers that were not yet fitted at the time of 0.1.3.
• The figure gallery now includes a source-map table that maps each display to its fitted object or fixture, extractor or plotter, interval source, and current support boundary.
• The figure gallery now shows the supported fixed-effect univariate mu emmeans route, including factor-conditioned and interaction grids, an empirical marginal_parameters() summary, and unsupported boundaries for non-mu or blocked emmeans targets.
• The figure gallery now separates residual sigma, ordinary group-level SDs, conditional random-slope deviations, and fitted sd(site) surfaces, with unavailable random-effect SD intervals shown as an explicit plotting boundary.
• The visualization grammar now records the Florence closeout: plot_parameter_surface() and plot_corpairs() remain the exported helpers, most gallery-specific displays stay as ggplot2 recipes, and simulation or failure-ledger helpers wait for stable Phase 18 result schemas.
• The figure gallery now includes distributional-parameter panels for mu, sigma, Student-t nu, zero-inflation probability zi, and residual rho12, with explicit response-scale labels and interval provenance.
• The pkgdown site now includes a user-facing figure gallery for model interpretation plots, confidence bands, correlation displays, random-effect scale surfaces, and simulation operating-characteristic figures.
• The pkgdown site now includes an improving-convergence guide explaining when the default optimizer budget is enough, how to use drm_control(optimizer = ...), how to interpret check_drm() rows, and when to separate optimization from Wald uncertainty with se = FALSE. The guide now also separates residual rho12, phylogenetic mean-mean, and ordinary group-level correlation boundaries for bivariate structured fits, and warns that larger data sets help only when they add information that separates those covariance layers.
• Phase 18 now has a pre-simulation readiness matrix that separates fitted, smoke-tested, interval-ready, weak-boundary-tested, planned, and blocked surfaces before broad simulation reports are written.
• Phase 18 now has a meta_V(V = V) summary-smoke runner that executes vector and dense known-covariance pilot replicates and returns grouped bias, RMSE, and Monte Carlo error summaries.
• Phase 18 now has a Gaussian location-scale summary-smoke runner that executes two or more pilot replicates and returns grouped bias, RMSE, and Monte Carlo error summaries from the smoke output.
• Phase 18 now has Monte Carlo uncertainty helpers for simulation summaries, including MCSEs for mean error, RMSE, proportions, and explicit interval-coverage summaries when lower and upper interval columns are present.
• Phase 18 now has a synthetic interval-coverage smoke helper for testing coverage-table plumbing before Wald, profile, or bootstrap interval methods are attached.
• Phase 18 smoke reports can now display optional aggregate, manifest, and warning/error ledger CSVs while still rendering when those files are not supplied.
• Phase 18 smoke report rendering is now covered by a skip-aware test with tiny aggregate, manifest, and warning/error ledger CSV fixtures.
• Phase 18 interval work now has a producer contract for Wald, profile, and bootstrap interval tables, including reported scale, method, status, and correlation-scale rules.
• Phase 18 now has a Fisher-z back-transformed Wald interval helper for correlation summaries, complementing raw-correlation Wald intervals from the generic helper.
• Phase 18 now has a generic Wald interval-table helper for parameter summaries that already contain estimates and standard errors, recording interval method, reported scale, status, and failure messages.
• Phase 18 now has a small aggregation helper for simulation parameter summaries, reporting replicate counts, bias, RMSE, absolute error, empirical standard error, convergence rate, Hessian rate, warning rate, and elapsed time by explicit grouping columns.
• Phase 18 now has a first end-to-end meta_V(V = V) smoke runner under inst/sim/run/, covering vector and dense known sampling covariance cells through DGP, drmTMB() fit, saved RDS output, and combined parameter summaries.
• Phase 18 meta_V(V = V) pilot summaries now carry standard errors for estimated mu coefficients and response-scale fitted residual sigma, while keeping known sampling covariance V out of interval targets.
• Phase 18 meta_V(V = V) summary-smoke output now includes Wald interval rows and coverage summaries for estimated mu coefficients and fitted residual sigma.
• Phase 18 now has a first end-to-end Gaussian location-scale smoke runner under inst/sim/run/, wiring the cell registry, seeded DGP, drmTMB() fit, pilot summariser, saved RDS output, and combined parameter table for one small surface.
• Phase 18 Gaussian location-scale pilot summaries now carry fixed-effect standard errors when the fitted model exposes them, preparing that surface for real Wald interval coverage checks.
• Phase 18 Gaussian location-scale summary-smoke output now includes formula-coefficient Wald interval rows and coverage summaries.
• Phase 18 now has a Gaussian meta-analysis meta_V(V = V) simulation pilot covering vector and dense known sampling covariance, including smoke tests that keep known V out of interval targets.
• Phase 18 now has a Gaussian location-scale simulation pilot: optional inst/sim/ helpers generate mu ~ x, sigma ~ z data and summarise one fitted model into a truth/estimate/error table.
• Phase 18 now has an optional inst/sim/ skeleton with reproducible seed-table and cell-registry helpers plus a CRAN-safe smoke test for simulation-run bookkeeping.
• Phase 18 now has an ADEMP-style simulation-programme blueprint in docs/design/41-phase-18-simulation-programme.md, including first-wave surfaces, estimands, methods, performance measures with Monte Carlo uncertainty, and the next three implementation slices.
• The meta-analysis tutorial and design examples now use meta_V(V = V) as the preferred known sampling covariance spelling, with meta_known_V(V = V) described only as a compatibility alias.
• meta_V() interval safety is now tested for Gaussian meta-analysis fits: profile_targets() keeps estimated sigma, random-effect SD, and bivariate rho12 targets visible while never treating known sampling covariance V as an estimated confidence-interval target.
• meta_V() now gives a clearer reserved-boundary error for proportional sampling-variance arguments such as meta_V(w = w, scale = "proportional"), meta_V(w = w), or meta_V(V = V, scale = "exact"). Diagonal/vector meta_V(V = V) can still use ordinary likelihood weights, while full matrix-V fits reject non-unit weights until joint-block weighting has a separate design.
• meta_V(V = V) is now accepted as the preferred additive known sampling covariance marker for Gaussian meta-analysis, routing to the same likelihood path as meta_known_V(V = V). The proportional branch meta_V(w = w, scale = "proportional") remains deliberately unimplemented and errors before fitting.
• Slice 204 now records the meta_V() API decision: the preferred additive known-covariance spelling is meta_V(V = V), without a positional response/value argument, and meta_known_V(V = V) is retained as a compatibility alias rather than a separate likelihood path.
• Slice 203 now records the post-202 Phase 17 return block for meta-analysis hardening. The next targets are the preferred meta_V(V = V) spelling and compatibility story, additive vector/matrix known V, proportional-variance design boundaries, interval safety, and reader examples.
• The Slice 202 pre-simulation gate now keeps broad Phase 18 comprehensive simulation closed until the post-202 Phase 17 hardening block is complete. A narrow Poisson mu random-effect pilot simulation is allowed, but meta-analysis hardening around meta_V()/known V, interval safety, and reader-facing examples should come before broad simulation claims.
• Focused Poisson random-effect recovery added a factor-predictor random-intercept case and a weak-SD boundary case. Historical note, superseded by later NB2, ordinal, and exact row-specific structured/component gates: this entry establishes only the ordinary Poisson mu path and does not define the current class-wide boundary.
• The validation-debt register now contains a Slice 201 non-Gaussian failure ledger. It names the convergence, boundary, identifiability, interval, and runtime risks that Phase 18 should measure or exclude before broad simulation claims are made.
• animal() and relmat() are now exported and documented as planned structured-effect markers, giving the reference index the intended animal -> phylo -> spatial -> lower-level known-dependence path while keeping gr() as a reserved legacy marker. These markers are parsed for roadmap examples and rejected by drmTMB() until fitted likelihood, diagnostics, profile-target, and recovery-test evidence exists (#147).
• summary(conf.int = TRUE) now handles fitted non-Gaussian models whose summary has no parameter rows to receive intervals, including cumulative-logit ordinal models. Wald fixed-effect intervals are still reported where fixed effects exist, and empty coefficient or parameter tables keep explicit interval-status columns instead of erroring.
• The model map, family chooser, and structural-dependence article now teach the structural-dependence ladder in biological order: animal(), phylo(), spatial(), combined phylogenetic-plus-spatial models, and lower-level relmat() known-dependence matrices. The same pages keep the non-Gaussian random-effect boundary visible before the comprehensive simulation phase.
• Shape random effects now have a dedicated unsupported-boundary message and Student-t test gate. Current Student-t nu formulas remain fixed-effect tail-shape models outside the exact Q-Series nu ~ phylo(1 | id, tree = tree) local-fit gate; skew-normal fixed-effect nu models residual slant, while skew-normal and skew-t shape random effects need separate likelihood recovery before nu/tau random effects or ID-level skew(id) ~ x style models are added.
• Poisson mu now supports ordinary unlabelled random intercepts and independent numeric random slopes for non-zero-inflated Poisson models, such as (1 | group) and (0 + x | group) on the log-mean predictor. The fitted SDs appear in sdpars$mu, ranef()/random_effects$mu, and profile_targets() as direct log_sd_mu targets; correlated Poisson random-slope blocks, labelled covariance blocks, zero-inflated Poisson random effects, and cross-parameter non-Gaussian covariance remain planned.
• Non-Gaussian sigma random effects have family-specific gates. Ordinary NB2, lognormal, and Gamma admit one independent log-sigma random intercept; Student-t, beta, beta-binomial, truncated NB2, and hurdle NB2 remain fixed-effect only. Slopes, labelled or combined scale effects, and structured scale effects still require likelihood, recovery, extractor, interval, and documentation evidence before fitting.
• Bivariate random-slope boundary errors now distinguish the fitted slope-only mu1/mu2 target, the fitted matching q=2 sigma1/sigma2 scale-slope target, the fitted same-response q2 mu/sigma slope target, the smoke-artifact-routed q=4 and q=6 location-only targets, the diagnostic-artifact-routed first q8 all-endpoint target, and broader p8/q8 endpoint variants that remain closed.
• Gaussian location random-slope blocks now support ordinary unstructured numeric multi-slope mu terms such as (1 + x1 + x2 | id) and labelled variants. The first public path has q=3 recovery, sdpars$mu, corpars$re_cov, corpairs(), summary(), and profile_targets() coverage; larger q blocks are advanced fits whose sample-size cost remains a simulation target.
• Gaussian location-scale models now fit more than one independent matched labelled mu/sigma random-intercept covariance block, for example matching (1 | p | id) and (1 | q | site) terms in both formulas. Each block reports its own corpars$mu_sigma, corpairs(class = "mean-scale"), summary(), and profile_targets() row.
• check_drm() and profile-interval tests now cover two independent univariate mu/sigma random-intercept covariance blocks, reporting one diagnostic row per mean-scale block and confirming the second eta_cor_mu_sigma profile target.
• Gaussian residual-scale random slopes now have an explicit multiple-independent-term boundary: sigma ~ z + (1 | id) + (0 + w1 | id) + (0 + w2 | id) fits separate log-sigma random-effect SDs with correlations fixed at zero. Historical note, superseded by the current 0.6.0 guidance: unlabelled ordinary correlated blocks such as (1 + x | id) and multi-slope variants are now fitted; labelled residual-scale slope covariance and cross-formula mu-sigma slope covariance remain planned.
• The earlier phylogenetic random-slope boundary was superseded by the structured-slope parity slices: phylo(1 + x | species, tree = tree) now fits one univariate Gaussian mu slope as an independent intercept and slope field. Historical note, superseded again by the current 0.6.0 guidance: the exact q1 phylogenetic sigma one-slope route is fitted and inference-ready with caveats; additional multiple or labelled phylogenetic-slope layouts outside the exact fitted ledger cells and structured slope correlations remain planned.
• Random-effect planning now includes a one-slope-per-layer gate before the non-Gaussian revisit, separating fitted ordinary Gaussian mu, Gaussian sigma, univariate mean-scale, bivariate intercept, phylogenetic intercept, and coordinate-spatial one-slope surfaces from the remaining Gaussian double-hierarchical limits.
• Spatial one-slope coverage now confirms a profile-likelihood interval for the coordinate-spatial mu slope-field SD. Historical note, superseded by the current 0.6.0 guidance: the q1 spatial sigma one-slope route now has point-fit/extractor evidence, but its interval gate remains blocked; additional multiple or labelled spatial-slope layouts outside the exact fitted ledger cells and spatial slope correlations remain planned.
• The double-hierarchical endpoint map now reflects the current Gaussian boundary after the random-slope gate: q > 2 ordinary mu, independent sigma slopes, multiple univariate mean-scale intercept blocks, coordinate-spatial one-slope support, matching bivariate q=4/q=6 location smoke routes, the matching bivariate q=2 scale-slope route, the first same-response q2 mu/sigma slope route, and the first q8 all-endpoint route are fitted, while q8 coverage/power evidence and spatial q=4 covariance surfaces remain planned.
• Non-Gaussian random-effect planning now has a first mu random-intercept gate: Poisson and NB2-style count likelihoods are the first candidates, while continuous, bounded, ordinal, zero-inflation, hurdle, shape, and structured non-Gaussian random-effect paths keep explicit unsupported messages until their recovery tests exist.
• emmeans::emmeans() now supports the first narrow drmTMB path: fixed-effect univariate mu estimated marginal means with retained model frames and fixed-effect covariance available. Generic emmeans pairwise contrasts on that returned mu grid are covered by a small parity test.
• The public emmeans::emmeans() boundary tests now confirm that transformed-response formulas such as log(y) ~ x error before an emmGrid is returned, keeping the first bridge limited to untransformed response formulas and explicit transformed-scale prediction tables through prediction_grid().
• The public emmeans::emmeans() boundary tests now confirm that bivariate Gaussian fits error with the unsupported "biv_gaussian" model type before an emmGrid is returned, instead of falling through to a generic missing-mu message.
• The public emmeans::emmeans() boundary tests now confirm that cumulative-logit ordinal fits still error before an emmGrid is returned, with guidance toward prediction_grid() and predict_parameters() for explicit prediction tables.
• The public emmeans::emmeans() boundary tests now confirm that hurdle NB2 fits still error before an emmGrid is returned, with guidance toward prediction_grid() and predict_parameters() for explicit prediction tables.
• The public emmeans::emmeans() boundary tests now confirm that zero-inflated NB2 fits error with the unsupported "zi_nbinom2" model type before an emmGrid is returned, matching the existing zero-inflated Poisson boundary.
• The public emmeans::emmeans() boundary tests now confirm that zero-inflated Poisson fits still error before an emmGrid is returned, with guidance toward prediction_grid() and predict_parameters() for explicit prediction tables.
• The fixed-effect univariate mu emmeans::emmeans() path now explicitly checks the type = "response" argument path, so response-scale EMMs requested directly from emmeans() must match predict(dpar = "mu", type = "response").
• Fixed-effect prediction matrices and the first univariate mu emmeans::emmeans() bridge now preserve fitted ordered-factor coding when newdata or an emmeans reference grid supplies the same levels as an ordinary factor, so ordered polynomial columns still align with fitted coefficients.
• Fixed-effect prediction matrices now accept character newdata values for fitted factor levels, ignore unused factor columns, and reject unknown or missing factor levels with a clear predictor-specific error before model-matrix construction.
• Fixed-effect prediction matrices now validate that newdata supplies every predictor required by the requested distributional parameter and that required predictor values are complete before model-matrix construction.
• Fixed-effect prediction matrices now reject non-finite numeric values such as Inf in required predictors before model-matrix construction.
• Fixed-effect prediction matrices now reject newdata values that produce non-finite transformed-predictor columns, such as log(size) when size = 0, before returning a non-finite prediction.
• marginal_parameters() and predict_parameters() now have explicit coverage for fitted random-effect scale model names such as sd(id), returning the random-effect-sd-model component, preserving row labels in long prediction tables, and averaging supplied direct-SD rows in marginal summaries.
• prediction_grid() now has explicit integration coverage for fitted direct-SD predictors: grids over predictors such as w in sd(id) ~ w can feed predict_parameters(..., dpar = "sd(id)") and marginal_parameters(..., by = "w").
• Random-effect scale predictions now accept character newdata values for fitted sd(group) factor levels and reject unknown levels with a predictor-specific error before random-effect scale model-matrix construction.
• Random-effect scale predictions now have explicit boundary coverage for newdata containers: non-data-frame inputs error, while zero-row data frames return named length-zero numeric vectors on both link and response scales.
• Random-effect scale predictions now have explicit coverage that multi-row newdata returns one value per supplied row, preserves rownames(newdata), uses response scale by default, and matches exp(link) when compared with type = "link".
• Random-effect scale predictions now have explicit multiple-target coverage: when a fit includes formulas such as sd(id) ~ w_id and sd(site) ~ w_site, each requested dpar validates its own required predictors, ignores sibling-target extras, and names the missing target-specific predictor.
• Random-effect scale predictions now validate the raw predictors required by sd(group) formulas in supplied newdata, so missing columns, missing values, and non-finite numeric values error before random-effect scale model-matrix construction.
• Random-effect scale predictions such as predict(fit, dpar = "sd(id)", newdata = ...) now reject newdata values that produce non-finite transformed-predictor columns, such as log(w_pos) when w_pos = 0, before returning an infinite link- or response-scale SD prediction.
• The fixed-effect univariate mu emmeans::emmeans() path now explicitly checks factor-conditioned reference grids such as emmeans(fit, ~ habitat | season, at = list(x = 0.25)), so conditional EMM rows must match predict(dpar = "mu") with the same factor levels.
• The fixed-effect univariate mu emmeans::emmeans() path now explicitly checks interaction formulas such as habitat * x on an explicit reference grid, so conditional EMMs must match predict(dpar = "mu") at the same interaction design point.
• The fixed-effect univariate mu emmeans::emmeans() path now explicitly checks multiple numeric at values, so conditional grids such as emmeans(fit, ~ habitat | x, at = list(x = c(-0.25, 0.75))) must match row-wise predict(dpar = "mu") on the same grid.
• The fixed-effect univariate mu emmeans::emmeans() path now explicitly checks cov.reduce = FALSE, so EMMs that average over unreduced numeric covariate levels must match predict(dpar = "mu") averaged over the same reference grid.
• The fixed-effect univariate mu emmeans::emmeans() path now explicitly checks custom numeric covariate reduction, so cov.reduce = stats::median must move the reference grid to median(x) rather than the default mean.
• The fixed-effect univariate mu emmeans::emmeans() path now explicitly checks the default numeric covariate-reduction rule, so emmeans(fit, ~ habitat) must match predict(dpar = "mu") at the mean of the numeric covariate used by the reference grid.
• The fixed-effect univariate mu emmeans::emmeans() path now has explicit parity coverage for formulas with offset(log(exposure)), so exposure-adjusted count-rate grids must match predict(dpar = "mu") on both link and response scales.
• The fixed-effect univariate mu emmeans::emmeans() path now has explicit recover-data coverage for transformed predictors such as log(size), so reference grids supplied through at are checked against predict(dpar = "mu").
• Bivariate, zero-inflated, hurdle, ordinal expected-score, random-effect, structured-effect, non-mu, slope, and interval-specialized emmeans targets still error before an emmGrid is returned.
• The model-workflow article now shows the first optional emmeans::emmeans() example for fixed-effect univariate mu, keeping adjusted means separate from predict_parameters() tables and from unsupported sigma, random-effect, bivariate, zero-inflated, hurdle, ordinal, and slope workflows. Broader drmTMB-specific contrast helpers remain a separate future contract.
• The model-workflow article now shows how to build an explicit prediction_grid() for a fitted random-effect scale model such as sd(site) ~ reef_cover, then pass that grid through predict_parameters(..., dpar = "sd(site)") and marginal_parameters(..., by = "reef_cover") without treating random-effect SDs as residual sigma or raw responses.
• The model-map article now routes fitted random-effect SD surfaces through prediction_grid(), predict_parameters(..., dpar = "sd(group)"), and marginal_parameters(), with the random-effect-sd-model component kept separate from residual sigma.
• predict_parameters(conf.int = TRUE) now adds Wald fixed-effect confidence intervals for supplied newdata grids when the requested distributional parameter has an ordinary fixed-effect basis. The table fills std.error, conf.low, conf.high, conf.level, conf.status = "wald", and interval_source = "wald" for supported rows, while fitted-row requests and direct random-effect scale models keep explicit unavailable interval status.
• Installation docs now point tagged-preview users to pak::pak("itchyshin/drmTMB@v0.1.2").
• docs/design/39-visualization-grammar.md now records the Phase 17 visualization and marginal-effects research contract. The note uses ggplot2, tidybayes, ggdist, emmeans, ggeffects, marginaleffects, diagnostic plotting packages, and figure-composition tools as design sources while keeping drmTMB data-first and dependency-light. The model-workflow article now states that predict_parameters() and marginal_parameters() are data tables that plotting helpers can consume, not plotters themselves.
• plot_corpairs() now provides the first optional ggplot2 display for explicit corpairs() tables. It draws one point per fitted correlation row, adds interval segments only when finite conf.low and conf.high bounds have supported interval provenance, can facet by a supplied table column such as level, and keeps correlation level, class, display interval status, and interval source attached to the plotted data.
• plot_parameter_surface() now provides the first optional ggplot2 plotting helper for long tables returned by predict_parameters(). It plots existing point estimates, keeps interval provenance columns attached to the data, and leaves EMMs, contrasts, and slope plots for later tested helpers.
• plot_parameter_surface() now draws confidence bands for continuous x-values and interval bars for discrete x-values when the supplied table already contains finite conf.low and conf.high bounds with real conf.status and interval_source provenance. It still does not compute confidence intervals, and rows with interval_source = "not_available" remain visibly interval-free.
• The model-workflow and model-map articles now show prediction-surface confidence bands as a table-first workflow: prediction_grid() -> predict_parameters(conf.int = TRUE) -> plot_parameter_surface(), with conf.status, conf.level, and interval_source left visible.
• plot_parameter_surface() now labels single-parameter panels with the fitted distributional parameter and prediction scale, such as sigma estimate (response scale), while keeping the generic Estimate label when multiple parameters are plotted together.
• prediction_grid() now builds explicit newdata grids for predict_parameters() and marginal_parameters(). The first contract supports focal predictors, supplied values, conditioned nuisance predictors, mean-reference grids, and empirical counterfactual grids while recording the grid rule as metadata.
• The bivariate-coscale tutorial now shows a fitted corpairs() table flowing into plot_corpairs(..., facet = "level"), separating residual rho12 from group-level correlation rows in the displayed workflow.
• docs/design/39-visualization-grammar.md now records the pre-export contract that plot_corpairs() follows: consume explicit corpairs() tables, keep correlation levels/classes visible, draw intervals only from finite confidence bounds with supported provenance, and test residual, ordinary group-level, phylogenetic, derived-unavailable, empty-table, and missing-ggplot2 cases before export.
• The model-map article now includes a Phase 17 visualization decision table that routes raw responses, fitted parameter surfaces, empirical marginal summaries, correlations, interval tables, and diagnostics to the current data helpers before readers choose a plotting style.
• The Reference index now makes the post-fit path explicit: fitting, checking, summaries, predictions, uncertainty, and extractors are grouped under “Model fitting and post-fit tools”, while exported plotting helpers appear under “Visualization”. The current exported plotting helpers are plot_parameter_surface() and plot_corpairs().
• predict_parameters() and marginal_parameters() now include interval provenance columns. The first contract reports conf.status = "not_requested" and interval_source = "not_available" so downstream tables and future plots cannot imply confidence intervals that were not computed.
• confint(), summary(), and corpairs() originally rejected bootstrap interval methods while the simulate-refit contract was still untested. That boundary is now superseded for confint(..., method = "bootstrap") on selected direct fitted-object targets; method = "parametric_bootstrap" and bootstrap routing through summary() or corpairs() still fail explicitly.
• summary() now reports delta-method standard errors for direct response-scale parameter rows, including constant sigma, residual rho12, random-effect SDs, and random-effect correlations, when TMB::sdreport() succeeds. Descriptive fitted ranges and derived variance ratios keep missing standard errors. This release originally recommended profile likelihood broadly for fitted SD and correlation targets; the cell-specific 0.6.0 guidance above supersedes that blanket recommendation.
• summary() profile summaries now keep fixed-effect Wald 95% confidence intervals while adding profile-likelihood 95% confidence intervals for selected direct targets such as sigma. Printed parameter tables no longer show duplicated minimum and maximum columns for constant direct parameters where those values equal the estimate.
• The model-workflow article now includes a compact guide to reading ordinary summary() output. The guide maps coefficients, parameters, covariance, derived, and confint components to the interpretation task and then points readers to fixef(), sigma(), rho12(), ranef(), corpairs(), and profile_targets() when they need more detail.
• The model-workflow article now shows empirical marginalization with prediction_grid(..., margin = "empirical") and marginal_parameters(..., by = "temperature"), separating conditioned prediction rows from averages over the fitted-row covariate distribution.
• The model-workflow article now adds raw-data-plus-model display rules for Phase 17: show observed responses on the observed-response scale, draw fitted mu and sigma surfaces from an explicit predict_parameters() table, keep interval provenance visible, and do not place raw response points on sigma, sigma^2, rho12, random-effect SD, or correlation axes.
• The Gaussian location-scale tutorial now includes a response-scale interpretation ladder for fixed mean slopes, residual-SD and residual-variance ratios, random-slope SDs, residual-scale random-slope SDs, and sd(group) slopes. It also adds a trait-named parrot beak-length equation block defining mu, sigma, body mass, forest habitat, and scale coefficients before the worked growth example. The worked growth example shows profile_targets(fit_growth) and a compact fitted translation table so readers can report mean growth, predictability, and among-group variation on the correct scale.
• The phylogenetic-spatial tutorial is now framed as a structural-dependence route: first phylogeny, then coordinate spatial dependence, then the planned phylogeny-plus-spatial endpoint. The article defines the conceptual combined equation, names the phylogenetic and spatial SDs, and keeps simultaneous phylo() plus spatial() syntax marked as planned until multiple structural mu layers have implementation and identifiability checks.

drmTMB 0.1.1 (2026-05-10)

• docs/design/34-validation-debt-register.md now backs the stable-core matrix with an evidence and debt ledger. Each advertised surface is marked as covered, partial, opt-in, or blocked, with tests, diagnostics, interval status, docs, and explicit debt recorded before the project expands the surface.
• README and the “What can I fit today?” model-map article now include a stable-core feature matrix. The matrix separates stable fitted surfaces, first-slice implementations, opt-in large-data controls, and planned or rejected neighbouring syntax, with interval and diagnostic status attached to each row.
• check_drm() now reports full-matrix meta_known_V(V = V) fits as dense known-covariance notes with retained dimension, storage, density, size, rank, and conditioning, making clear that dense known covariance is a small-to-moderate path until sparse or block-sparse V storage has implementation and benchmark evidence.
• confint() profile rows now include lightweight profile diagnostics through profile.boundary and profile.message. Successful intervals currently report "ok" unless transformed SD intervals are close to zero or transformed correlation intervals are close to the correlation boundary; profile failure messages now explicitly point to boundary, one-sided, non-monotone, or failed-inner-optimization profiles as possible causes.
• confint() and summary() now use an explicit conf.status column for interval output. Successful confint() rows report conf.status = "wald" or "profile", while summary(conf.int = TRUE) marks parameter rows that need newdata, are derived-only, are ready but not selected in the current call, or are unavailable for Wald intervals.
• summary() and profile_targets() now expose the first derived variance-ratio summaries without claiming derived confidence intervals. Simple Gaussian random-intercept repeatability and univariate phylogenetic signal appear as point-estimate rows with target_type = "derived", transformation = "variance_ratio", and profile_note = "derived_target". When intervals are requested, these rows report derived_interval_unavailable; confint(..., method = "profile") fails before launching an unsupported derived profile.
• Direct profile-likelihood intervals for random-effect SDs and correlations now have focused Slice 55 coverage across the currently fitted ordinary, phylogenetic, and coordinate-spatial surfaces. The tests verify the first spatial sd:mu:spatial(1 | site) interval, ordinary and phylogenetic constant corpairs(conf.int = TRUE) rows, bivariate phylogenetic SD and mean-mean correlation targets in summary(conf.int = TRUE, method = "profile"), and the continued separation of derived covariance intervals from direct SD/correlation intervals.
• Row-specific profile-likelihood intervals now have focused coverage for response-scale sigma, sigma1, sigma2, residual rho12, and fitted q=2 ordinary or phylogenetic corpair() values supplied through newdata. The tests verify multi-row bivariate scale intervals, fitted latent-correlation intervals for both ordinary and phylogenetic q=2 corpair() routes, and early errors for ambiguous newdata requests such as multiple parm values, non-data-frame inputs, or empty grids.
• profile_targets() now treats its returned table as a tested namespace contract. Target rows use controlled target_type, profile_ready, profile_note, and transformation values, and memory-light fits created with drm_control(keep_tmb_object = FALSE) now keep listing direct target names while marking them with profile_note = "tmb_object_required" instead of implying that direct profile intervals can be run.
• confint() now wraps direct TMB::tmbprofile() calls with clearer target-specific errors. Users can still tune profile controls such as ystep, ytol, and parm.range, but drmTMB now blocks attempts to override the internal obj, name, lincomb, or trace arguments through ... and reports the profile_targets() name when profiling or profile-interval extraction fails.
• drm_control(se = FALSE) now skips TMB::sdreport() while keeping optimized fits usable for coefficients, fitted values, residuals, prediction, simulation, log-likelihood, and profile-likelihood routes that retain fit$obj. Fits also survive sdreport() failure with fit$uncertainty$status = "failed", while summary(), vcov(), and check_drm() report the skipped or failed uncertainty state explicitly.
• drm_control() now reserves future start, fixed-parameter map, fallback-optimizer, and multi-start control names so they cannot be silently passed to nlminb() through a plain optimizer list. Profile-likelihood calls also re-pin the TMB object to the selected opt$par before profiling, keeping mutable TMB state aligned with the chosen optimum.
• drm_control(aggregate_gaussian = TRUE) now fits the first sufficient-statistic aggregation path for univariate Gaussian fixed-effect models. Repeated rows are grouped by processed mu and sigma design state, TMB evaluates the Gaussian likelihood with n, sum(y), and sum(y^2) cells, and fitted-row predictions and residuals remain one value per original model row. Random effects, direct-SD formulas, structured effects, known sampling covariance, bivariate models, non-Gaussian families, non-unit likelihood weights, and combined sparse fixed-effect matrices remain planned.
• drm_control(keep_model_frame = FALSE) now also drops nested model-frame caches for direct random-effect SD models and fitted q=2 corpair() regression models after their model matrices and group metadata have been retained. This keeps the memory-light fitted-object path aligned with sd_phylo() and latent-correlation features.
• drm_control(sparse_fixed = TRUE) now fits the first sparse fixed-effect path for univariate Gaussian mu fixed effects with intercept-only sigma. The fitted object keeps the mu design as a sparse Matrix, prediction follows the fitted sparse flag, and check_drm() reports sparse fixed-effect design matrices; random effects, known covariance, phylogenetic or spatial terms, bivariate models, non-Gaussian models, and sparse sigma remain planned.
• bench/large-phylo-location.R now records aggregate_gaussian, requested and fitted aggregation-cell counts, aggregation compression ratio, and largest aggregation cell size, and can run a non-phylogenetic aggregation smoke benchmark with --structured none --aggregate-gaussian true --aggregation-cells 100.
• bench/large-phylo-location.R now records structured and sparse_fixed settings and can run a non-phylogenetic sparse fixed-effect smoke benchmark with --structured none --factor-heavy true --sparse-fixed true.
• check_drm() now includes the density of the largest retained fixed-effect design block in the fixed_effect_design_size row, making high-cardinality mostly-zero designs easier to distinguish from genuinely dense designs and confirming when a fitted object retains sparse fixed-effect matrices.
• The optional bench/large-phylo-location.R benchmark now records the largest retained fixed-effect design block, its column count, nonzero count, and density, and bench/summarize-results.R includes those fields when present.
• Univariate Gaussian phylogenetic location models now support the Family B direct-SD formula sd_phylo(species) ~ x_species. The implementation uses a non-centred unit phylogenetic base effect scaled at observed tips by species-level SD predictors, giving marginal tip covariance D_tip A_tip D_tip without assigning predictors to internal tree nodes. coef(), predict(), sdpars, and profile_targets() expose the fitted SD surface.
• Bivariate Gaussian phylogenetic location models now support Family B direct-SD formulas sd_phylo1(species) ~ x_species and sd_phylo2(species) ~ x_species for matching mu1 and mu2 phylogenetic location effects. The implementation keeps the latent phylogenetic mean-mean correlation constant, exposes response-specific species SD surfaces through coef(), predict(), and sdpars, and rejects mixtures with all-four q=4 phylogenetic location-scale blocks.
• biv_gaussian() now supports Family B direct location random-effect SD formulas for labelled bivariate location random intercepts: sd1(id) ~ x_group targets the mu1 random-effect SD and sd2(id) ~ x_group targets the mu2 random-effect SD. Predictors must be constant within the named group, and unsupported scale-random-effect SD names such as sd_sigma1() / sd_sigma2() plus same-group q=4 Family A mixtures are rejected to avoid mixing direct SD models with scale-formula random effects.
• biv_gaussian() now supports the first predictor-dependent latent random-effect correlation models for q=2 location blocks. Ordinary grouped blocks use corpair(id, level = "group", block = "p", from = "mu1", to = "mu2") ~ x_group; phylogenetic blocks use corpair(species, level = "phylogenetic", block = "p", from = "mu1", to = "mu2") ~ ecology beside matching phylo(1 | p | species, tree = tree) terms. The fitted link-scale coefficients appear in coef(), summary(), vcov(), and profile_targets(), while corpairs() reports the response-scale mean, range, and number of group/species correlation values. Predictors must be constant within group/species; location-scale, scale-scale, q=4, and spatial corpair() regressions remain planned.
• The fitted q=2 phylogenetic corpair() route uses two independent unit phylogenetic fields with species-specific loadings. This gives a positive-definite nonstationary covariance model, preserves the same-species local correlation interpretation, and reduces to the existing constant bivariate phylogenetic covariance when the correlation predictor is constant. A CRAN-safe broad-trend recovery test now checks that a positive species-level correlation predictor recovers the ordering of fitted phylogenetic correlations without hitting the correlation guard. The first implementation target is mu1-mu2; phylogenetic location-scale and scale-scale correlation regressions remain q=4 extensions.
• A new large-data workflow article documents current memory-light fit controls, practical post-fit output cautions, and the optional bench/large-phylo-location.R benchmark harness for Gaussian phylogenetic location models.
• check_drm() now reports optimizer evaluation counts, dense fixed-effect design size, finite fixed-effect standard errors, near-boundary random-effect standard deviations, univariate mu/sigma mean-scale covariance diagnostics, bivariate same-response mu/sigma diagnostics, bivariate mu1/mu2 and sigma1/sigma2 random-intercept covariance diagnostics, ordinary q=4 bivariate location-scale covariance diagnostics, bivariate phylogenetic mu1/mu2 covariance diagnostics, phylogenetic q=4 location-scale covariance diagnostics, coordinate-spatial mu diagnostics, and univariate or bivariate sd_phylo*() direct-SD surface diagnostics, helping users diagnose large, difficult, or weakly identified fits before interpreting estimates.
• biv_gaussian() now fits matching intercept-only phylo(1 | species, tree = tree) or labelled phylo(1 | p | species, tree = tree) terms in mu1 and mu2 as correlated phylogenetic location effects. It also fits the first matching labelled all-four phylogenetic q=4 block across mu1, mu2, sigma1, and sigma2, reporting four endpoint SDs and all six latent phylogenetic correlations while keeping residual rho12 separate. Partial, unlabelled, mismatched, and slope phylogenetic q=4 forms remain rejected.
• profile_targets() lists the fitted-model target names that can be passed to confint(), including whether each target is ready for direct profile-likelihood intervals. It distinguishes group-level covariance targets such as cor:mu_sigma:cor(mu:(Intercept),sigma:(Intercept) | p | id), cor:mu_sigma:cor(mu1:(Intercept),sigma1:(Intercept) | p | id), cor:mu:cor(mu1:(Intercept),mu2:(Intercept) | p | id), cor:sigma:cor(sigma1:(Intercept),sigma2:(Intercept) | p | id), and cor:phylo:cor(mu1:(Intercept),mu2:(Intercept) | phylo | species) from residual rho12. Ordinary q=4 theta_re_cov correlations and full phylogenetic q=4 theta_phylo correlations are listed as derived unstructured-correlation targets, while block-diagonal phylogenetic q=4 fallback fits expose direct constant block-correlation targets.
• confint() now returns Wald fixed-effect confidence intervals by default and can compute profile-likelihood intervals for explicit direct targets such as fixef:mu:x, constant sigma, sd:mu:(1 + x | id):(Intercept), sd:mu:phylo(1 | species), cor:mu:cor((Intercept),x | id), the first univariate mu/sigma, bivariate mu1/mu2, bivariate sigma1/sigma2, block-diagonal bivariate phylogenetic mu1/mu2 and sigma1/sigma2 random-effect correlations, and constant residual rho12. It also profiles row-specific response-scale sigma, sigma1, sigma2, rho12, and fitted q=2 ordinary or phylogenetic corpair() values when newdata is supplied. Full phylogenetic q=4 correlations are currently reported as derived targets, and direct fallback targets still need fit-specific profile diagnostics before being interpreted as usable intervals.
• corpairs() now accepts conf.int = TRUE for profile-likelihood correlation-pair intervals where the fitted target is profile-ready. Rows that are not interval-ready, such as predictor-dependent residual rho12 summaries that need newdata or derived q=4 unstructured-correlation rows, now carry an explicit conf.status instead of silently omitting bounds.
• biv_gaussian() now fits the first bivariate group-level covariance blocks: matching labelled random intercepts in mu1/mu2 and in sigma1/sigma2. The fitted group-level SDs appear in sdpars$mu or sdpars$sigma, the same-parameter random-intercept correlations appear in corpars$mu or corpars$sigma and corpairs(), and residual rho12 remains a separate within-observation correlation.
• biv_gaussian() now fits same-response cross-parameter random-intercept covariance blocks, such as matching (1 | p | id) terms in mu1 and sigma1, or a separate (1 | q | id) pair in mu2 and sigma2. The fitted mean-scale correlations appear in corpars$mu_sigma, corpairs(), and profile_targets().
• biv_gaussian() now fits an intercept-only ordinary q=4 location-scale covariance block when the same labelled (1 | p | id) term appears in mu1, mu2, sigma1, and sigma2. The block estimates four group-level SDs and all six latent random-effect correlations while keeping residual rho12 separate.
• corpairs() now accepts group and block filters so users can directly subset fitted group-level covariance rows while keeping residual rho12 rows separate. It also accepts location-class aliases such as class = "location-location" and class = "location-scale" for the existing mean-mean and mean-scale rows, matching the reserved corpair() formula terminology without renaming current output.
• drm_control() is now exported and provides the first large-data storage controls for drmTMB(): users can pass optimizer settings through optimizer = list(...), drop stored complete-case data with keep_data = FALSE, drop stored model frames after fitting with keep_model_frame = FALSE, and drop the retained TMB automatic-differentiation object with keep_tmb_object = FALSE.
• drm_formula() now reserves explicit coefficient-specific random-effect SD syntax such as sd(id, dpar = "mu", coef = "x1") ~ x_group for future random-slope scale models. drmTMB() rejects these formulas until the likelihood, covariance diagnostics, and simulation tests exist.
• drm_formula() now uses singular endpoint-specific corpair(group, level = "...", block = "...", from = "mu1", to = "mu2") ~ x syntax for predictor-dependent latent random-effect correlations. The first fitted paths are ordinary and phylogenetic q=2 mu1/mu2; spatial, location-scale, scale-scale, and q=4 variants remain parsed or documented as later targets. Use rho12 = ~ x for residual correlation and corpairs() to extract fitted latent correlations.
• drmTMB() now fits the first labelled cross-formula covariance block for univariate Gaussian location-scale models: matching y ~ x + (1 | p | id) and sigma ~ z + (1 | p | id) random intercepts. The fitted mean-scale correlation appears in corpars$mu_sigma and corpairs().
• Gaussian residual-scale random slopes are now implemented for univariate Gaussian sigma formulas as independent terms such as sigma ~ z + (0 + w | id). Historical note, superseded by the current 0.6.0 guidance: unlabelled ordinary correlated blocks such as (1 + x | id) are now fitted; labelled residual-scale slope covariance and cross-formula mu-sigma slope covariance remain planned.
• Installation docs now point tagged-preview users to pak::pak("itchyshin/drmTMB@v0.1.1").
• marginal_parameters() averages long-format distributional-parameter predictions over fitted rows or supplied newdata groups, providing the first simple marginalisation surface for mean, scale, shape, and residual-correlation summaries.
• nbinom2() and the zero-inflated, zero-truncated, and hurdle NB2 routes now share an internal count-kernel helper that avoids observed-count loops for large counts while preserving the small-overdispersion Poisson limit. Deterministic high-count tests compare the optimized objective against independent stats::dnbinom() calculations.
• Developer documentation now includes a C++ modularization source map that identifies safe header-only helper extraction, hidden model_type probe branches, required test gates, and the template pieces that should not move during the first refactor pass.
• Phase 6d stable-core validation is locally closed with a stable-core feature matrix, validation-debt register, failure-safe sdreport() controls, optimizer/start/map contract, dense covariance guardrails, count-kernel hardening, and a C++ modularization source map. GitHub Actions remains the PR-side gate.
• predict_parameters() returns long-format predictions for fitted distributional parameters such as mu, sigma, nu, and rho12, giving interpretation tables and future plotting or marginalisation helpers one shared data surface.
• summary() now reports a response-scale parameter table for fitted scale, shape, random-effect SD, and correlation quantities, with opt-in Wald or profile-likelihood confidence intervals through conf.int = TRUE, including direct profile intervals for the first fitted group-level covariance rows. It also includes a covariance component with fitted random-effect variance and covariance point summaries for currently fitted registry-backed covariance blocks and the first bivariate phylogenetic mu1/mu2 mean-mean row; derived covariance intervals remain unavailable until a nonlinear interval method is implemented, and the covariance table marks that interval status explicitly.

drmTMB 0.1.0 (2026-05-10)

• bf() now stores parsed formula entries for distributional parameters, including bivariate rho12, meta-analysis meta_known_V(V = V), and random-effect scale syntax.
• beta() now fits fixed-effect beta mean-scale models for strict continuous proportions in (0, 1), using logit(mu), log(sigma), and internal precision phi = 1 / sigma^2; fitted() returns mu and sigma(fit) returns the public scale parameter.
• beta_binomial() now fits fixed-effect beta-binomial mean-overdispersion models for cbind(successes, failures) responses, using logit(mu), log(sigma), and internal beta precision phi = 1 / sigma^2. fitted() returns the success probability mu, and sigma(fit) returns the public extra-binomial variation scale.
• biv_gaussian() now fits fixed-effect bivariate Gaussian location-scale-coscale models with separate mu1, mu2, sigma1, sigma2, and predictor-dependent rho12 formulas.
• check_drm() now provides a first-pass diagnostic table for optimizer convergence, fixed gradients, Hessian status, dropped rows, scale positivity, rho12 boundaries, Student-t nu boundary behaviour, known sampling covariance summaries, random-effect replication, and weak random-slope design checks.
• corpairs() now returns a long table of fitted correlation pairs that already exist in a drmTMB fit, currently residual bivariate rho12 summaries and ordinary group-level mu random-effect correlations.
• drm_formula() is now the primary formula constructor; bf() remains a short alias.
• drm_formula(mvbind(y1, y2) ~ x) is now implemented as shorthand for identical bivariate Gaussian location formulas, expanding internally to mu1 = y1 ~ x and mu2 = y2 ~ x.
• drmTMB() now fits Gaussian location-scale models with fixed effects, random intercepts, labelled random intercepts such as (1 | p | id), independent numeric random slopes, and ordinary labelled or unlabelled correlated random intercept-slope blocks in the mu formula, such as bf(y ~ x1 + (1 | id) + (0 + x1 | id), sigma ~ x1), bf(y ~ x1 + (1 + x1 | id), sigma ~ x1), and bf(y ~ x1 + (1 + x1 | p | id), sigma ~ x1).
• cumulative_logit() now fits fixed-effect univariate ordinal location models for ordered responses, using ordered cutpoints, Pr(y_i <= k) = logit^-1(theta_k - mu_i), and a fixed latent logistic scale. fitted() returns the expected ordered-category score, and ordinal scale/discrimination formulas remain planned.
• drmTMB() now fits fixed-effect Gamma mean-CV models for positive responses with family = Gamma(link = "log"), where mu is the response mean and sigma is the coefficient of variation.
• drmTMB() now fits fixed-effect Poisson mean models for count responses with family = poisson(link = "log"), including standard R exposure offsets in the mu formula such as count ~ habitat + offset(log(trap_nights)). The same family route supports fixed-effect zero-inflated Poisson models via zi ~ predictors, where mu is the conditional Poisson mean, zi is the structural-zero probability, and fitted() returns (1 - zi) * mu. Overdispersion, random effects, known sampling covariance, and bivariate count models remain planned for this route.
• drmTMB() now supports a top-level weights = argument for non-negative row log-likelihood multipliers, with weights(fit) returning the processed weights after model-row filtering. Known sampling variance or covariance remains separate and should use meta_V(V = V), with meta_known_V(V = V) as a compatibility alias.
• nbinom2() now fits fixed-effect negative-binomial 2 mean-dispersion models for overdispersed counts, with log(mu), log(sigma), and Var(y) = mu + sigma^2 * mu^2. The mu formula supports standard R exposure offsets such as count ~ habitat + offset(log(trap_nights)). Here sigma is an overdispersion scale, not a residual standard deviation or size parameter. The same family route supports fixed-effect zero-inflated NB2 models via zi ~ predictors.
• drmTMB() now accepts family = c(gaussian(), gaussian()) and family = list(gaussian(), gaussian()), routing both to the implemented bivariate Gaussian location-coscale likelihood. Mixed composed families such as c(gaussian(), poisson()) remain planned and currently error clearly.
• Gaussian residual-scale random intercepts are implemented in the sigma formula, for example bf(y ~ x1 + (1 | id), sigma ~ x1 + (1 | id)). These model residual-scale heterogeneity and are distinct from random-effect scale formulae such as sd(id) ~ x_group.
• Gaussian random-effect scale formulae are implemented for one or more distinct unlabelled mu random intercepts, for example bf(y ~ x1 + (1 | id) + (1 | site), sigma ~ x2, sd(id) ~ x_group, sd(site) ~ site_type). Each sd(group) predictor must be constant within the named group after missing-row filtering.
• Gaussian known-covariance meta-analysis with meta_known_V(V = V) is now covered by targeted validation when combined with random-effect scale formulae such as sd(id) ~ x_group, using an independent dense marginal-likelihood comparator.
• Gaussian mu random-effect correlations from correlated blocks are exposed as corpars$mu, keeping group-level labels such as p separate from residual bivariate rho12.
• Profile-likelihood confidence intervals are documented as a planned inference phase with an explicit target namespace, such as sd:mu:(1 | id), cor:mu:cor((Intercept),x | id), and fixef:rho12:(Intercept), plus boundary flags and nonlinear derived-quantity guidance.
• deviance(), df.residual(), and nobs() now work for drmTMB fits, making base-R model summaries and comparison helpers more complete.
• fitted() now returns family-specific response summaries: mu for implemented Gaussian-like, Gamma, beta, beta-binomial, Poisson, and NB2 mean models, the expected ordered-category score for cumulative-logit ordinal models, the arithmetic response mean for lognormal models, (1 - zi) * mu for zero-inflated Poisson and zero-inflated NB2 models, (1 - hu) * mu / (1 - Pr_NB2(0)) for hurdle NB2 models, and a two-column mu1/mu2 matrix for bivariate Gaussian models.
• fixef() now returns distributional fixed-effect coefficients and acts as a mixed-model-friendly alias for coef().
• lognormal() now fits fixed-effect univariate lognormal location-scale models for positive responses, with mu and sigma defined on the log-response scale, fitted() returning the arithmetic response mean, and simulation plus likelihood tests checked against stats::dlnorm().
• meta_known_V(V = V) now fits Gaussian meta-analysis with diagonal or dense full known sampling covariance using family = gaussian().
• meta_vcov_bivariate() now builds row-paired dense sampling covariance matrices for bivariate Gaussian meta-analysis with known within-study covariance, and meta_known_V(V = V) now fits complete-row bivariate Gaussian known-V models by adding that sampling covariance to the fitted residual covariance from sigma1, sigma2, and rho12.
• ranef() now returns fitted conditional random-effect blocks, including ordinary mu, residual-scale sigma, phylo_mu, and the first spatial_mu blocks when present.
• rho12() now returns response-scale residual correlations from bivariate Gaussian location-coscale fits, with type = "link" available for Fisher-z-like linear predictors using the guarded transform rho12 = 0.999999 * tanh(eta_rho12).
• student() now fits fixed-effect univariate Student-t location-scale-shape models with mu, sigma, and nu formulas. The nu parameter is modelled as nu = 2 + exp(eta_nu) for a stable finite-variance robust continuous family.
• truncated_nbinom2() now fits zero-truncated negative-binomial 2 models for positive counts, with ordinary mu random intercepts and independent numeric slopes allowed in non-hurdle models. mu and sigma describe the untruncated NB2 component, fitted() returns the conditional positive-count mean mu / (1 - Pr_NB2(0)), and sigma(fit) returns the NB2 overdispersion scale. Adding hu ~ predictors still fits the corresponding fixed-effect hurdle NB2 model; hurdle random effects outside the exact Q-Series hu ~ relmat(1 | id, Q = Q) local-fit gate, correlated zero-truncated random slopes, and sigma random effects remain planned.
• drmTMB() now fits phylogenetic random intercepts and one numeric phylogenetic random slope in the univariate Gaussian location formula with phylo(1 | species, tree = tree) and phylo(1 + x | species, tree = tree), using an ultrametric branch-length tree and the sparse augmented A-inverse path. It also fits coordinate-based spatial random intercepts and one numeric spatial mu slope in the univariate Gaussian location formula with spatial(1 | site, coords = coords) and spatial(1 + x | site, coords = coords). The slope paths estimate independent intercept and slope fields with the same fixed structured precision and separate SDs, labelled with terms such as phylo(1 | species), phylo(0 + x | species), spatial(1 | site), and spatial(0 + x | site).
• animal() and relmat() now fit one numeric univariate Gaussian mu slope beside their fitted random-intercept paths, for example animal(1 + x | id, pedigree = ped) and relmat(1 + x | id, K = K), using independent intercept and slope fields with separate SDs. Planned structured-effect markers outside the first fitted one-slope paths, such as standalone or partial phylogenetic scale terms, spatial terms in sigma, spatial(1 | site, mesh = mesh), multiple structured slopes, slope correlations, spatial q=4 blocks, and predictor-dependent spatial corpair() formulas, are parsed by drm_formula() and rejected by drmTMB() with planned-feature errors until their TMB likelihoods and recovery tests are implemented.
• The “Which scale are you modelling?” tutorial now includes a copy-run scale audit with fitted output and interpretations for sigma ~, weights =, preferred meta_V(V = V), sd(group) ~, and bivariate rho12 ~ syntax.
• Tutorial prose now clarifies that ordinary likelihood weights and known sampling covariance are separate concepts, and that dense full meta_V(V = V) paths currently reject non-unit weights; deprecated meta_known_V(V = V) remains only a compatibility alias.
• The Gaussian location-scale tutorial now includes a worked growth example with equations, fitted summary() output, response-scale sigma interpretation, and a table mapping mean growth and residual SD back to the biological question.
• The bivariate location-coscale tutorial now includes worked activity-boldness equations, output-reading guidance for rho12, and a response-scale residual-correlation curve along a disturbance gradient.
• The bivariate Gaussian coscale phase now has a closure audit in the roadmap: rho12(), corpairs(), bivariate known sampling covariance, row likelihood weights, mvbind() shorthand, residual diagnostics, and unsupported bivariate random-effect syntax are all documented as implemented or planned in one place.
• The response-family tutorial now starts with an at-a-glance table that maps common measurement processes to implemented families, distributional parameters, the meaning of sigma, and the main current limitation.
• The meta-analysis tutorial now includes a worked restoration example with equations, fitted summary() output, response-scale residual heterogeneity interpretation, and a clearer distinction between preferred meta_V(V = V) and ordinary likelihood weights =.
• The phylogenetic-spatial tutorial now includes a worked thermal-tolerance example for the implemented phylo(1 | species, tree = tree) path, with equations, fitted output, tree/species validation guidance, bivariate phylogenetic and coordinate-spatial reading guidance for corpairs(), summary(fit)$covariance, direct profile targets, the first labelled q=4 phylogenetic location-scale syntax, and clearer marking of implemented coordinate-spatial intercept, one-slope, and q=2 bivariate location paths versus planned mesh/SPDE and spatial q=4 paths.
• The getting-started article and pkgdown tutorial menu now provide a clearer learning path from biological or statistical questions to the matching tutorial and distributional parameter.
• Public documentation now pairs symbolic model equations with matching R syntax for the first Gaussian location-scale, random-effect scale, bivariate rho12, meta-analysis, and phylogenetic examples, and clarifies planned spatial coords versus mesh inputs.
• The likelihood design now includes a central TMB model_type routing table, aligned with the implemented source map, including the hidden phylogenetic prior parity branch used only by tests.
• residuals() now returns whitened Pearson residuals for bivariate Gaussian fits, and vcov() now uses coefficient-level row and column names.
• Initial project scaffold.