diff --git a/src/TMB/ll_ltriangle.hpp b/src/TMB/ll_ltriangle.hpp index 4af48961..dc9fd0e5 100644 --- a/src/TMB/ll_ltriangle.hpp +++ b/src/TMB/ll_ltriangle.hpp @@ -85,8 +85,11 @@ Type ll_ltriangle(objective_function* obj) { Type absdev = CppAD::CondExpGe(dev, Type(0.0), dev, -dev); Type inside = scalelog - absdev; // floor the density argument to keep the likelihood finite (and act as a - // soft barrier) if a candidate support fails to cover a data point. - Type floored = CppAD::CondExpGt(inside, Type(1e-8), inside, Type(1e-8)); + // soft barrier) if a candidate support fails to cover a data point. The + // floor is scaled to `scalelog` (the peak density argument) so the soft + // barrier behaves consistently regardless of the magnitude of the data. + Type floor = Type(1e-8) * scalelog; + Type floored = CppAD::CondExpGt(inside, floor, inside, floor); nll -= weight(i) * (log(floored) - Type(2.0) * log(scalelog) - log(left(i))); }; if(left(i) < right(i)){ // censored values diff --git a/tests/testthat/test-ltriangle.R b/tests/testthat/test-ltriangle.R index 9387764d..680196ce 100644 --- a/tests/testthat/test-ltriangle.R +++ b/tests/testthat/test-ltriangle.R @@ -25,3 +25,34 @@ test_that("ltriangle", { expect_snapshot_value(ssd_rltriangle(2), style = "deparse") }) }) + +test_that("ltriangle fit is invariant to scaling the concentrations", { + data <- ssddata::ccme_boron + fit <- ssd_fit_dists(data, dists = "ltriangle") + + data_scaled <- data + data_scaled$Conc <- data_scaled$Conc * 1000 + fit_scaled <- ssd_fit_dists(data_scaled, dists = "ltriangle") + + est <- estimates(fit) + est_scaled <- estimates(fit_scaled) + + # scaling concentrations by a constant only shifts locationlog by its log and + # leaves scalelog unchanged, so hazard concentrations scale by the constant. + # The scale-relative density floor keeps this invariance across magnitudes. + expect_equal( + est_scaled$ltriangle.scalelog, + est$ltriangle.scalelog, + tolerance = 1e-5 + ) + expect_equal( + est_scaled$ltriangle.locationlog, + est$ltriangle.locationlog + log(1000), + tolerance = 1e-5 + ) + expect_equal( + ssd_hc(fit_scaled, proportion = 0.05)$est, + ssd_hc(fit, proportion = 0.05)$est * 1000, + tolerance = 1e-4 + ) +})