Add log-triangular (ltriangle) distribution

[joethorley]

Add a log-triangular distribution where log(concentration) follows a symmetric triangular distribution with mode locationlog and half-width scalelog. Wired in following the llogis pattern:

• R/ltriangle.R: ssd_pltriangle/qltriangle/rltriangle/eltriangle, the base symmetric-triangular helpers (p/q/rtriangle_ssd), the concentration-scale helpers (pltriangle_ssd/qltriangle_ssd) used in model averaging, and sltriangle starting values.
• src/TMB/ll_ltriangle.hpp: TMB negative log-likelihood (CondExp-based CDF, soft floor on the density) registered in ssdtools_TMBExports.cpp.
• dist_data: ltriangle row (bcanz=FALSE, tails=FALSE, npars=2, valid=TRUE, bound=FALSE) - bounded support, no tails, no bounded parameter.
• ssd_pmulti/qmulti/rmulti and params.R documentation extended.
• tests/testthat/test-ltriangle.R mirrors test-llogis.R (test_dist plus snapshots); test-dists.R lists updated; affected snapshots regenerated.
• distributions.Rmd gains a log-triangular section.

[claude]

Claude Code Review

Claude Code Review is paused for this repository. To reconnect it, an admin of this repository's GitHub organization (or the account owner, for personal repositories) who can also manage your Claude organization's Code Review settings needs to re-link GitHub in Code Review settings. This is a one-time step.

_{Tip: disable this comment in your organization's Code Review settings.}

[beckyfisher]

PR #169 review findings — ltriangle (log-triangular) distribution

Repo: poissonconsulting/ssdtools · PR: #169 · Author: joethorley
Base: dev · Head: add-triangle-distribution · 10 commits, 25 files, +638/−7
Reviewer: Rebecca Fisher · Reviewed at: 8b9ea754 ("test log-triangular")
Status: DRAFT — read-only review, nothing submitted to GitHub.

---

1. Phase 1 — mechanical checks

Check	Result	Attribution
devtools::document() sync	PASS (with caveat)	Local artifact, not a PR defect
devtools::check()	1 ERROR, 1 NOTE	Both pre-existing/environmental
devtools::test()	FAIL 3 / PASS 1358 / SKIP 39	Pre-existing on dev — not this PR
Snapshot stability	PASS	No ltriangle snapshot rewritten
air format --check .	PASS	Clean, exit 0 (air.toml present)

1.1 Docs in sync — PASS with caveat

document() rewrote only NAMESPACE, and only as a formatting change: the committed file uses
multi-line importFrom(pkg, a, b) blocks; my locally installed roxygen2 8.0.0 collapses them to
one-per-line importFrom(pkg,a). No export()/S3method() line changed, and all four ltriangle
exports are present and correct:

export(ssd_eltriangle) export(ssd_pltriangle) export(ssd_qltriangle) export(ssd_rltriangle)

Discarded with git checkout -- NAMESPACE man/. No action for the author — this is a roxygen2
version skew on my machine. (Worth noting DESCRIPTION carries no RoxygenNote field, so there is
nothing pinning the expected roxygen version.)

1.2 devtools::check() — 1 ERROR, 1 NOTE

• checking whether package 'ssdtools' can be installed ... [57s] OK — the TMB backend including
  ll_ltriangle.hpp compiles cleanly. No compiler warnings surfaced.
• checking examples ... [20s] OK, re-building of vignette outputs ... [20s/22s] OK
  (so the new distributions.Rmd section knits).
• ERROR — checking tests: the 3 test-censor.R snapshot failures below.
• NOTE — hidden files: Found the following hidden files and directories: .git. Artifact of
  checking from a git worktree. Not a PR issue.

1.3 devtools::test() — 3 failures, all pre-existing

FAILURE 'test-censor.R:47:3'  Snapshot of `path` has changed.
FAILURE 'test-censor.R:54:3'  Snapshot of `path` has changed.
FAILURE 'test-censor.R:61:3'  Snapshot of `path` has changed.
[ FAIL 3 | WARN 0 | SKIP 39 | PASS 1358 ]

Root cause: the committed snapshot header is
Chemical,Species,Conc,Group,Units,Medium,right; the regenerated one is
Chemical,Species,Conc,Group,Units,right. My installed ssddata 1.0.0 has
ccme_boron columns Chemical,Species,Conc,Group,Units — no Medium column.

Attribution confirmed: I built a second worktree at upstream/dev (ff47b9ce) and ran the
censor tests there — identical 3 failures, [ FAIL 3 | PASS 30 ]. This is an ssddata version
mismatch in my environment, not caused by PR #169. (CI installs ssddata from
open-AIMS/ssddata, which presumably still has Medium.) No action for the author.

1.4 Snapshot stability — PASS

Running the suite rewrote no ltriangle-related snapshot. _snaps/ltriangle.md,
_snaps/estimates.md, _snaps/multi.md and _snaps/data/dist_data.csv were all stable — so the
snapshots committed in the PR are deterministic and current. The only churn was the 3 censor
.new.csv files (above) and deletion of .png snapshots, which happens because plot tests are
skip_on_os-skipped on Linux and testthat prunes them as unused. All discarded with
git checkout -- tests/testthat/_snaps.

1.5 Air — PASS

air format --check . exits 0 with no output. The PR is fully Air-formatted.

---

2. Phase 2A — code-reading findings

Overall the PR follows the llogis pattern closely and faithfully; the wiring is consistent,
alphabetically ordered, and complete. Divergences found are concentrated in the two places where
bounded support genuinely differs from llogis.

F1 — R/pqr.R:116-124 · quantile endpoints ignore bounded support (non-blocking, but new class of bug)

.qd() hard-codes the extreme quantiles before ever dispatching to the distribution's own function:

if (p == 0) { if (!.lgt) return(0); return(-Inf) }
if (p == 1) { return(Inf) }
do.call(fun, args = args)

This is correct for every pre-existing distribution (all have support unbounded above, and →0
below). ltriangle is the first distribution in the package with bounded support, so it is the
first to violate the assumption. Measured on locationlog=0, scalelog=sqrt(6):

	returned	correct (= its own limit)
ssd_qltriangle(0)	0	exp(-sqrt(6)) = 0.086338
ssd_qltriangle(1)	Inf	exp(+sqrt(6)) = 11.582435
ssd_qltriangle(1e-12)	0.08633793	— (converges correctly)

So the function is discontinuous at its own endpoints: it converges to 0.0863 as p→0 then
jumps to 0 at exactly p=0. Note this is a defect in shared infrastructure, not in
R/ltriangle.R — flagging it as something the PR exposes rather than introduces. Practical
blast radius is small (ssd_hc/ssd_hp don't evaluate at p=0/1), which is why nothing else caught
it. Fix would be to let qdist defer to the distribution when it has finite support.

F2 — src/TMB/ll_ltriangle.hpp:86-93 · the density floor is an attractor, not a barrier (BLOCKING)

Detailed in §3 with numerical evidence. Summary: the "soft floor" makes the negative log-likelihood
diverge to −∞ as scalelog → 0 once observations fall outside the candidate support, so it
rewards exactly the configuration it was meant to penalise.

F3 — src/TMB/ll_ltriangle.hpp:95-102 · censored branch has no floor at all (blocking-adjacent)

The uncensored branch is floored; the censored branch is not:

nll -= weight(i)*log(pright-pleft);   // no guard

If a censored interval lies entirely outside the candidate support, ptri returns the same value at
both ends, pright-pleft == 0, and log(0) → -Inf → non-finite objective. Reproduced: fitting
interval-censored data with an interval far below the cluster fails with
L-BFGS-B needs finite values of 'fn'. This asymmetry is inherited structurally from ll_llogis.hpp
(which also has no guard) — but for llogis the CDF is strictly increasing everywhere, so the
difference is never exactly 0. Bounded support makes it reachable.

F4 — dist_data flags — all correct ✅

ltriangle, bcanz=FALSE, tails=FALSE, npars=2L, valid=TRUE, bound=FALSE.

• bound=FALSE is genuinely correct, and the prompt's suspicion is resolved: R/data.R:52
  documents bound as "Whether one or more parameters have boundaries", consumed via
  is_bounds() (R/helpers.R:133) to apply L-BFGS-B box constraints. ltriangle's parameters are
  locationlog (unconstrained) and log_scalelog (unconstrained, exp-transformed). Parameter
  bounds ≠ support bounds. Correct as written.
• tails=FALSE is defensible (matches invpareto, the other non-tailed distribution) and has
  no computational effect — tails is used only as a filter in ssd_dists() (R/dists.R:52-53).
  Minor wording tension: R/data.R:49 glosses tails as "has both tails", and a symmetric
  triangular does have two (finite) tails. Reading it as "has both infinite tails" is what makes
  FALSE right. Not worth changing, possibly worth a doc clarification.

F5 — Model-averaging wiring — correct ✅

ssd_pmulti/ssd_qmulti/ssd_rmulti each gain ltriangle.weight=0, ltriangle.locationlog=0,
ltriangle.scalelog=3, alphabetically placed between lnorm_lnorm and weibull, matching the
existing convention exactly. R/params.R:199-201 documents all three. pltriangle_ssd /
qltriangle_ssd mirror pllogis_ssd / qllogis_ssd precisely.
Verified end-to-end: averaging lnorm + llogis + ltriangle on ccme_boron works, and ltriangle
takes weight 0.706 (best AIC, log_lik = -116, delta = 0).

F6 — R/ltriangle.R:116-145 · base helpers — correct, minor divergence in the good direction ✅

• ptriangle_ssd CDF branches verified against closed form.
• qtriangle_ssd: z = sqrt(2p) - 1 for p<=0.5, 1 - sqrt(2(1-p)) above — exactly van
  Straalen's inverse.
• rtriangle_ssd uses runif(n) + runif(n) - 1: the sum of two U(0,1) is triangular on [0,2] with
  mode 1, shifted to [-1,1] mode 0. Variance 2 × 1/12 = 1/6, so Var = scale²/6 — matches the
  SD = scale/√6 identity. Clean and correct.
• Minor divergence: ptriangle_ssd/qtriangle_ssd return rep(NaN, length(q)) on scale <= 0
  where plogis_ssd/qlogis_ssd return scalar NaN. The ltriangle version is more correct
  (length-preserving). No action; noting only because it is a deliberate departure from the template.

F7 — R/ltriangle.R:100-106 · starting value — sound ✅

halfwidth <- max(abs(logx - location)) guarantees the initial support covers the data, with a
1.1× cushion; the <= 0 fallback to 0.1 handles zero-spread data and is covered by a test. This is
a real improvement over blindly mirroring sllogis. (It does not, however, prevent the optimiser
from later leaving the covering region — see §3.)

F8 — TMB gradient at the mode — acceptable, worth a comment (nit)

absdev = CondExpGe(dev, 0, dev, -dev) gives |dev|, so the density has a slope sign-flip at the
mode and the log-likelihood is piecewise-smooth with kinks at each observation. This is inherent
to the triangular family, not a defect — L-BFGS-B tolerates it and every fit I ran converged. The
kink set has measure zero in the parameter space. Worth one comment in the header acknowledging the
non-smoothness so a future reader doesn't assume smooth-optimiser guarantees.

F9 — Tests — stronger than the llogis template ✅

test-ltriangle.R (178 lines) does more than mirror structure. Genuinely valuable additions:

• test_dist("ltriangle") + the standard p/q/r snapshots (template parity);
• scale-invariance test — Conc * 1000 shifts locationlog by log(1000) and leaves
  scalelog unchanged (this is the test that justifies the relative floor);
• monotonic finite HC across 1:99/100, and finite predict() across the full range;
• zero-log-spread starting-value edge case;
• two external-reference fits: endosulfan (fitdistrplus, Hose & Van den Brink 2004) and a
  US EPA SSD Toolbox permethrin reproduction (EPA/600/R-18/116 Table 2) with provenance and a public
  domain note. Reproducing the EPA Toolbox's a/b support endpoints in log10 is a nice check.

Gap: no test covers the failure modes in §3 — no test fits data where the optimum wants to leave
the covering region, and none covers censored data (the left < right TMB branch is untested for
this distribution).

F10 — Vignette — accurate, missing provenance (nit)

distributions.Rmd pdf (s-|y-μ|)/s² and the two-branch cdf are both correct and match the
implementation. Framing is honest — explicitly says "experimental", "not part of the default BCANZ
set", "only fitted when explicitly requested".
Flagged as requested: the section cites only Wikipedia. No origin reference for the
distribution in an SSD context — van Straalen (2002) and Stephan et al. (1985) are the
natural citations, and the US EPA SSD Toolbox is already cited in the test comments but not in the
user-facing docs. Recommend adding them.

---

3. Phase 2B — analytic verification battery

3.1 Result: 14/14 PASS

API introspected first: ssd_qltriangle(p, locationlog=0, scalelog=3, lower.tail, log.p);
probe ssd_qltriangle(0.5, locationlog=0, scalelog=sqrt(6)) returned exactly 1.00000000, so these
functions work on the concentration scale and comparisons were taken through log().

Check	Actual	Expected	Result
q(0.50) = mode	0.000000	0.000000	PASS
q(0.05) [van Straalen headline]	−1.674893	−1.674893	PASS
q(0.025)	−1.901767	−1.901767	PASS
q(0.10)	−1.354045	−1.354045	PASS
q(0.25)	−0.717439	−0.717439	PASS
symmetry q(0.95) = −q(0.05)	1.674893	1.674893	PASS
F(mode)	0.500000	0.5	PASS
F(q0.05)	0.050000	0.05	PASS
F(lower support)	0.000000	0	PASS
F(upper support)	1.000000	1	PASS
E[log X] (n=2e5)	−0.000423	0	PASS
Var[log X] (n=2e5)	0.994099	1	PASS
RNG support respected	min −2.4407, max 2.4397	within ±2.44949	PASS
ordering: ltri most conservative	−1.6753 < −1.6450 < −1.6212	ltri < lnorm < llogis	PASS

The reviewer-relevant claim holds: at matched log-scale mean/SD the bounded triangular is more
conservative in the lower tail than the infinite-tailed families
(ltri −1.6753 vs lnorm −1.6450 vs llogis −1.6212; analytic −1.6749 / −1.6449 / −1.6232).

The distribution math (p/q/r) is correct. ssd_qltriangle reproduces van Straalen's −1.675 to
six decimal places. No FAIL to interpret.

Reminder, as instructed: van Straalen's fitted h/q/HC values were deliberately not
used as targets. He fits by nonlinear least-squares on empirical plotting positions; ssdtools fits
by maximum likelihood. Only Appendix A's estimator-independent closed-form properties were
asserted, and every check above calls the distribution functions directly with no fitting, so
§4.1 isolates the p/q/r implementation from the estimation machinery. Differences between his
published fits and an ssdtools refit would be an expected estimator difference, not a bug.

3.2 The estimation machinery is where the problem is

Because §3.1 isolates p/q/r, I ran a separate battery against the fitting path. Normal use is
healthy — on ccme_boron, ltriangle fits (locationlog=2.5096, scalelog=2.8316, support
[0.725, 208.8] covering data [1, 70.7]), gives HC1/HC5/HC10/HC20 = 1.082 / 1.774 / 2.571 / 4.344,
bootstrap CIs work (se=0.605, lcl=1.33, ucl=3.87, 100 boot), and model averaging works.
ssd_hp() returning 0 below the support bound is correct for a bounded family, not an artefact.

But the floor breaks down under outliers. Replicating the TMB objective exactly
(ll_ltriangle.hpp:84-93), for an observation outside the candidate support the contribution is:

log(1e-8 · s) − 2·log(s) − log(y)  =  log(1e-8) − log(s) − log(y)

As s → 0, −log(s) → +∞. So the log-likelihood increases without bound and the nll diverges
to −∞, when for a true triangular it should be +∞ (density is exactly 0 out there).
Measured, with locationlog parked far from all data so every point is floored:

log_scalelog=-2    nll=   563.78
log_scalelog=-5    nll=   470.78
log_scalelog=-20   nll=     5.78
log_scalelog=-50   nll=  -924.22      <- lower nll = "better" fit
log_scalelog=-200  nll= -5574.22

The floor is an attractor, not a barrier. Two further consequences, both measured:

1. Zero gradient in locationlog. In the floored regime the objective is independent of
   locationlog — nll is identically 604.293499 at locationlog = 40, 50, 60, 70. Once a point
   is outside the support, nothing pulls the mode back toward it.
2. The relative floor removes its own guard at the limit. floor = 1e-8 * scalelog underflows
   to exactly 0 when scalelog underflows, so log(0) → -Inf/NaN:
   log_scalelog=-745 → nll=Inf; -746 → nll=NaN. This is the direct source of the observed
   L-BFGS-B needs finite values of 'fn' fit failures. (An absolute floor does not have this
   particular underflow, but diverges faster at −2·log(s), so it is not the fix either — commit
   76496b21 made the divergence rate slower, not absent.)

Observable behaviour, tight cluster of 30 points at ~10 plus one low outlier:

data	fitted support	outcome
ccme_boron (well-behaved)	[0.725, 208.8]	0/28 excluded ✅
30 tight + outlier ×1e−2	[6.36, 15.27]	1/31 excluded
30 tight + outlier ×1e−4	[6.36, 15.27]	1/31 excluded
30 tight + outlier ×1e−6	[6.36, 15.27]	1/31 excluded
30 tight + outlier ×1e−10	—	fit fails

The fitted support is identical regardless of how extreme the outlier is — the signature of the
zero gradient in (1). And the transition is n-dependent, exactly as the algebra predicts
(excluding a point costs −log(1e-8) ≈ 18.42; shrinking scalelog by factor k gains 2n·log k,
so exclusion becomes profitable around n ≈ 14):

n	fitted support	excluded
5	[4.6e−09, 5.2e+07]	0/6 ✅
10	[1.9e−08, 1.5e+08]	0/11 ✅
20	[8.35, 13.19]	2/21
40	[7.62, 13.79]	1/41
80	—	fit fails

Guideline-relevant consequence. On a 31-point dataset with one sensitive outlier, silently and
with no warning:

ltriangle fitted support = [7.494, 15];  data min = 1e-07;  points below support = 1
ltriangle HC5 = 8.36301
lnorm     HC5 = 0.0270088
ssd_hp() at the excluded observation (1e-07) = 0 %

HC5 is ~310× less protective, because the fit discarded the most sensitive species and then
reported that species as having zero probability of occurring. ssd_fit_dists() emits no warning
that the fitted support fails to cover the data. For an SSD package this is the failure mode that
matters most, and it is precisely the one that mirroring llogis could never surface.

Mitigating context, and why I'd still call this fixable rather than fatal: bcanz=FALSE keeps this
out of the default guideline path, the vignette labels it experimental, and every well-behaved
dataset I tried (including the PR's own endosulfan and permethrin fixtures) fits correctly.

---

4. Categorisation

Blocking

• B1 — The density floor rewards leaving the support (ll_ltriangle.hpp:86-93). The nll diverges
  to −∞ as scalelog → 0 once points are floored, so the "soft barrier" is an attractor. Observable
  as (a) silent fits whose support excludes data, with HC5 ~310× less protective than lnorm on the
  same data, and (b) hard L-BFGS-B needs finite values of 'fn' failures. Suggested directions:
  make the out-of-support penalty increase with distance (e.g. penalise absdev - scalelog
  quadratically) instead of flattening; and/or constrain scalelog >= max|log(y) - locationlog| so
  the support is covering by construction — which for a triangular is where the true MLE lives
  anyway. At minimum, warn when the fitted support does not cover the data.
• B2 — Censored branch is unguarded (ll_ltriangle.hpp:95-102). pright - pleft can be exactly
  0 under bounded support → log(0) → non-finite objective → fit failure. Needs the same treatment
  as the uncensored branch. Also currently untested for this distribution.

Non-blocking

• N1 — qdist extreme quantiles ignore bounded support (R/pqr.R:116-124). ssd_qltriangle(0)
  returns 0 and (1) returns Inf instead of exp(m∓s), making the function discontinuous at its
  own endpoints. Shared infrastructure — exposed by this PR, not introduced by it. Could reasonably
  be a follow-up issue rather than a blocker on Add log-triangular (ltriangle) distribution #169.
• N2 — Test coverage for the failure modes above. Add a fit test with an extreme outlier
  asserting the fitted support covers the data (or that a warning is raised), and a censored-data fit
  test exercising the left < right branch.
• N3 — Missing provenance in distributions.Rmd. Add van Straalen (2002) and Stephan et al.
  (1985); consider surfacing the US EPA SSD Toolbox reference already present in the test comments.

Nits

• n1 — One header comment in ll_ltriangle.hpp noting the log-likelihood is piecewise-smooth
  with kinks at each observation (inherent to the triangular; not a defect).
• n2 — R/data.R:49 glosses tails as "has both tails"; a symmetric triangular does have two
  finite tails. Consider "both unbounded tails" so tails=FALSE reads unambiguously.
• n3 — DESCRIPTION has no RoxygenNote, so nothing pins the roxygen2 version that produces the
  committed NAMESPACE formatting (see §2.1). Pre-existing, repo-wide.

Explicitly verified as correct — no action

bound=FALSE (parameter bounds, not support bounds — correct); model-averaging wiring; the p/q/r
math (14/14 analytic checks, van Straalen's −1.675 to 6 dp); RNG construction and variance identity;
starting-value hardening; scale invariance; Air formatting; TMB compilation.

[beckyfisher]

Both blocking defects are resolved.

B1 is fixed at the mechanism. Writing the log density in terms of the dimensionless
t = 1 - |log(y) - locationlog| / scalelog leaves the in-support likelihood algebraically
unchanged, and makes the out-of-support contribution to nll grow without bound as scalelog
shrinks rather than fall without bound. That is the sign reversal the previous floor lacked.

B2 is fixed, and the follow-up found the case that mattered more than the one originally
reported: softening log() on the censored interval mass left a flat plateau, so a left-censored
non-detect below the data was still silently excluded. Non-detects at the sensitive end are the
normal case in SSD data. The censored branch now applies the same scale-free barrier as the
uncensored one, on the distance from the mode to the nearest point of the censoring interval.

The correction on the direction of the censored divergence is accepted: nll went to +Inf
there, not -Inf, so only its finiteness was broken.

N1 is resolved upstream in #195 by removing the hard-coded p == 0 and p == 1 branches from
.qd(), which is the general fix rather than a special case for ltriangle. N2, N3 and both nits
are done, and the four new regression tests cover both original failure modes and the non-detect
case.

The algebra and branch conditions were checked by reading the header at dbd4b674; the figures
in the replies above are the author's measurements and were not independently reproduced.

Two follow-ups, neither blocking this PR:

• The coverage assertion. The reason given for omitting it, that the barrier makes non-coverage
  unreachable rather than unlikely, held for one day before the censored non-detect case disproved
  it. A post-fit check that the fitted support covers every observation is independent of the
  likelihood implementation and would have caught that case without requiring the insight that
  found it. Worth an issue.
• Bootstrap fits. ssd_hc(ci = TRUE) refits many resampled datasets against a much stiffer
  objective, and the new tests are all single fits. A bootstrap test on a dataset with an extreme
  low observation would cover that interaction.

[beckyfisher]

The log-triangular has been looked at before in the ANZ/Canadian work and the position has been consistently negative. Fox et al. (2021, ETC) describe it as a curious inclusion in the USEPA SSD Toolbox, on the basis that its tail characteristics aren't ones encountered in practice, and note that its only formal regulatory role is the Stephan et al. (1985) four-point fit to the most sensitive genera rather than a whole-of-dataset SSD. Fox et al. (2022) §4.4.3 then names it explicitly as the example of a distribution that is unrealistic for ecotoxicology because it has no left or right tails, and excludes it from the candidate set on that basis before any simulation work. The 2024 final report leaves the default set unchanged.

None of that is an argument against having it available in the package. But it does mean anything that makes it reachable through ssd_fit_bcanz() or the default candidate set would run against a standing published recommendation, and under the consultation rules that's a MAJOR-version Technical Committee matter. Happy to help frame it that way in the docs if useful.