First CRAN-targeted release candidate. It was first submitted to CRAN on 24 August 2026; this resubmission candidate follows CRAN's example feedback and is not yet accepted or available on CRAN. drmTMB fits distributional regression models -- location, scale, shape, zero inflation, and residual correlation -- for one or two responses, using Template Model Builder.
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.
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.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.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() 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.
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.
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.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.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.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.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.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.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.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().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.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.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.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.drmTMB() help previously claimed
offset(log(exposure)) support for zero-truncated negative-binomial mu
formulas, which the code rejected.drmTMB(..., engine = "julia") against committed
expected.toml numbers (drmTMB 0.6.0 pin). Measured locally: 11/11
pass when JuliaCall + DRM.jl are available.phylo() term; large-p phylo routes for the
count / non-Gaussian set remain. Beta-binomial stays gated.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.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.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.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.confint(fit)
(method = "wald"). It is the fastest fitted-object route when
TMB::sdreport() succeeded.confint(fit, parm = ..., method = "profile") after profile_targets(fit)
lists the row as profile-ready.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).See also the vignette First-week intervals: fit, profile, and boundary.
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.?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.drmTMB() fit or reader-facing reporting permission.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/.
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.
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.
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() 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.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.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.
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.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.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.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.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.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.(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.||.|| 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.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.(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.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.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.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.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.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.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.pdHess = TRUE, and all structured-SD, structured-correlation, and RMSE
gates passed.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.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.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.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.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.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.sigma structure, joint mu/sigma,
simultaneous structured providers, spatial/animal, bivariate responses, REML, intervals, coverage,
and supported claims remain outside this arc.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.theta_hat, random-effect/structured residual adequacy, and bivariate
joint (non-marginal) outputs remain separately authorized future work.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.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.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.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.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.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 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.
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.
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).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().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.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.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.
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.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.
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.
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.
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.
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.
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).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.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.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.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.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.mu random slopes, independent Gaussian sigma random slopes, and coordinate-spatial Gaussian mu slopes.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.rho12 grid writers, keeping their Wald/profile/bootstrap interval artifacts separate from the baseline first-wave runner.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.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).sigma ~ phylo(...) boundary is now superseded by the fitted intercept-only structured sigma route described above.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.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.mu random intercepts and independent slopes as the current fitted mixed-count route.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.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.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.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.animal(), and relmat() lanes, while keeping planned or blocked lanes in the failure ledger instead of fitted grids.phase18_dgp_gaussian_ls() helper to aims, DGP conditions, estimands, methods, performance measures, and Williams-style reporting checks before larger grids run.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.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.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.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.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.sigma1/sigma2 scale-slope route is fitted separately.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.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.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.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.relmat() model classes before broad simulation grids admit them.inst/sim/ that captures warnings, errors, elapsed time, session metadata, and optional RDS output for pilot simulation cells.plot_corpairs() now has a label argument so publication figures can use short row labels while keeping full correlation metadata in the source table.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.inst/sim/reports/, giving pilot surfaces a reader-facing structure for purpose, aggregate summaries, reader checks, and interpretation boundaries.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.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.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.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.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.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.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.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.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.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.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.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.mu random-effect condition helpers now build true crossed condition grids, including random-effect SDs and, for NB2, overdispersion settings.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.mu random-effect pilot outputs, preparing aggregate, coverage, manifest, and failure tables for later figure-gallery work.mu random-effect pilots.mu random-effect pilot object.mu random-effect pilot and renders the Florence-facing HTML gallery from the resulting tables.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.mu emmeans route, including factor-conditioned and interaction grids, an empirical marginal_parameters() summary, and unsupported boundaries for non-mu or blocked emmeans targets.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.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.mu, sigma, Student-t nu, zero-inflation probability zi, and residual rho12, with explicit response-scale labels and interval provenance.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.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.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.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.meta_V(V = V) summary-smoke output now includes Wald interval rows and coverage summaries for estimated mu coefficients and fitted residual sigma.inst/sim/run/, wiring the cell registry, seeded DGP, drmTMB() fit, pilot summariser, saved RDS output, and combined parameter table for one small surface.meta_V(V = V) simulation pilot covering vector and dense known sampling covariance, including smoke tests that keep known V out of interval targets.inst/sim/ helpers generate mu ~ x, sigma ~ z data and summarise one fitted model into a truth/estimate/error table.inst/sim/ skeleton with reproducible seed-table and cell-registry helpers plus a CRAN-safe smoke test for simulation-run bookkeeping.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.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.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.meta_V(V = V) spelling and compatibility story, additive vector/matrix known V, proportional-variance design boundaries, interval safety, and reader examples.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.mu path and does not define the current class-wide boundary.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.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.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.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.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.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.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.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.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.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.mu, Gaussian sigma, univariate mean-scale, bivariate intercept, phylogenetic intercept, and coordinate-spatial one-slope surfaces from the remaining Gaussian double-hierarchical limits.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.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.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.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().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.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.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.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.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.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").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.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.newdata supplies every predictor required by the requested distributional parameter and that required predictor values are complete before model-matrix construction.Inf in required predictors before model-matrix construction.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").newdata values for fitted sd(group) factor levels and reject unknown levels with a predictor-specific error before random-effect scale model-matrix construction.newdata containers: non-data-frame inputs error, while zero-row data frames return named length-zero numeric vectors on both link and response scales.newdata returns one value per supplied row, preserves rownames(newdata), uses response scale by default, and matches exp(link) when compared with type = "link".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.sd(group) formulas in supplied newdata, so missing columns, missing values, and non-finite numeric values error before random-effect scale model-matrix construction.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.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.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.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.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.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.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.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.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").mu, slope, and interval-specialized emmeans
targets still error before an emmGrid is returned.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.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.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.pak::pak("itchyshin/[email protected]").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.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.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.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.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.prediction_grid(..., margin = "empirical") and
marginal_parameters(..., by = "temperature"), separating conditioned
prediction rows from averages over the fitted-row covariate distribution.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.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.phylo() plus spatial() syntax marked as planned until multiple structural mu layers have implementation and identifiability checks.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.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.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.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.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.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.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.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.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().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.pak::pak("itchyshin/[email protected]").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.model_type probe branches, required test gates, and the template pieces that should not move during the first refactor pass.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.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.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.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.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.mu random-effect correlations from correlated blocks are exposed as corpars$mu, keeping group-level labels such as p separate from residual bivariate rho12.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.sigma ~, weights =, preferred meta_V(V = V), sd(group) ~, and bivariate rho12 ~ syntax.meta_V(V = V) paths currently reject non-unit weights; deprecated meta_known_V(V = V) remains only a compatibility alias.summary() output, response-scale sigma interpretation, and a table mapping mean growth and residual SD back to the biological question.rho12, and a response-scale residual-correlation curve along a disturbance gradient.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.sigma, and the main current limitation.summary() output, response-scale residual heterogeneity interpretation, and a clearer distinction between preferred meta_V(V = V) and ordinary likelihood weights =.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.rho12, meta-analysis, and phylogenetic examples, and clarifies planned spatial coords versus mesh inputs.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.