Predict with climate outside the maxAsymptote ceiling when the ledger row asks for it

NEWS.md

@@ -1,5 +1,7 @@
 # fireSense_spreadPredict (development version)
 
+- A ledger row from a fit made with climate outside the `maxAsymptote` ceiling (`fireSense_spreadFit`'s new `spreadLink = "climateOutside"`) is predicted with that link: `spreadProbOneELF()` hands the row's stored climate reference (`fireSenseUtils::spreadFitClimateRefTxt`, `ledgerClimateRef()`) to `fireSenseUtils::spreadProbFromCovs()`, the one function the fit computes its spread probability with, so the climate term is added on the logit after the fuel ceiling exactly as in the fit; with several ELFs each row uses its own. `fireSense_burn` adds the year effect on the same logit, so nothing changes there. A row without a reference, such as every row written so far, predicts exactly as before (`spreadProbFromCovs()` with no reference is `logisticAll()`). Needs fireSenseUtils >= 0.2.3.9106.
+
 - A single fitted ELF now rescales its covariates with the `covMinMax_spread` of its ledger row, as several ELFs already did, not with `sim$covMinMax_spread`, which `fireSense_spreadFit` rebuilds in every run from that run's tables (in a hindcast, from a formula rebuilt for another year, which can lack a fitted term). A row without ranges (an older ledger) still uses `sim$covMinMax_spread`. A fitted covariate with no range now stops, naming it; before, it silently stayed on its raw scale.
 
 - A single fitted ELF now predicts with the covariates its fitted parameters name, as several ELFs already did, not the terms of `sim$fireSense_spreadFormula`. A hindcast whose formula had a term the fit did not (e.g. `treedWetland_agb`) stopped with "not found in data objects". Predictions where the formula matched the fit do not change. A fitted coefficient with no covariate still stops, naming it.

fireSense_spreadPredict.R

@@ -17,7 +17,7 @@ defineModule(sim, list(
   documentation = list("README.txt", "fireSense_spreadPredict.Rmd"),
   reqdPkgs = list("data.table", "magrittr", "Matrix", "methods", "terra", "SpaDES.core (>=3.0.4)", "stats",
                   "ggplot2", "viridis",
-                  "PredictiveEcology/fireSenseUtils@development (>= 0.2.3.9084)"),
+                  "PredictiveEcology/fireSenseUtils@development (>= 0.2.3.9106)"),
   parameters = bindrows(
     defineParameter(name = "lowerSpreadProb", class = "numeric", default = fireSenseUtils::spreadProbFloor,
                     desc = paste("Lower asymptote of the 2- and 3-parameter logistic. Default",
@@ -154,7 +154,7 @@ spreadPredictRun <- function(sim) {
     pred <- spreadProbOneELF(covs, params = params1, covMinMax = ledgerCovMinMax(sa, 1L, sim$covMinMax_spread),
                              yr = time(sim),
                              maxFireSpread = P(sim)$maxFireSpread, lowerSpreadProb = P(sim)$lowerSpreadProb,
-                             covCentre = ledgerCovCentre(sa, 1L))
+                             covCentre = ledgerCovCentre(sa, 1L), climateRef = ledgerClimateRef(sa, 1L))
     sim$fireSense_SpreadSD <- yearSpreadSDOf(params1)
   } else {
     ids <- as.character(sa[[fireSenseUtils::polygonIDTxt]])
@@ -175,7 +175,7 @@ spreadPredictRun <- function(sim) {
       p <- spreadProbOneELF(covs[these], params = parsI, covMinMax = sa$covMinMax_spread[[i]],
                             yr = time(sim), maxFireSpread = P(sim)$maxFireSpread,
                             lowerSpreadProb = P(sim)$lowerSpreadProb,
-                            covCentre = ledgerCovCentre(sa, i))
+                            covCentre = ledgerCovCentre(sa, i), climateRef = ledgerClimateRef(sa, i))
       rows <- these[match(p$pixelID, covs$pixelID[these])]
       acc[rows] <- acc[rows] + w[rows, i] * p$spreadProb
       accSD[rows] <- accSD[rows] + w[rows, i] * yearSpreadSDOf(parsI)
@@ -247,9 +247,13 @@ ELFblendWeights <- function(elfRas, template, pixelID, ids, width) {
 #'   from the rescaled covariates; `NULL` (a fit without an intercept, or an older ledger row) centres nothing.
 #'   A parameter set with `fireSenseUtils::spreadInterceptTxt` has an intercept, added as a column of 1s after
 #'   the covariates are rescaled and centred, as in the fit.
+#' @param climateRef the ELF's `climateRef_spread` (`ledgerClimateRef()`): with a fit made with climate outside the
+#'   `maxAsymptote` ceiling (`fireSense_spreadFit`'s `spreadLink = "climateOutside"`), the reference value of each
+#'   climate covariate; `fireSenseUtils::spreadProbFromCovs()` then applies that link, exactly as the fit did.
+#'   `NULL` (a fit with climate under the ceiling, or an older ledger row) is the link as it always was.
 #' @return `data.table` with `pixelID` and `spreadProb`, from the one parameter set in `params`.
 spreadProbOneELF <- function(covs, params, covMinMax, yr, maxFireSpread, lowerSpreadProb,
-                             covCentre = NULL) {
+                             covCentre = NULL, climateRef = NULL) {
   moduleName <- "fireSense_spreadPredict"
   covs <- copy(covs)
   ## the per-year random effect is not a covariate coefficient: fireSense_burn applies it (fireSense_SpreadSD)
@@ -333,8 +337,10 @@ spreadProbOneELF <- function(covs, params, covMinMax, yr, maxFireSpread, lowerSp
   # Make sure the order is correct in the matrix
   matching <- intersect(names(covPars), colnames(mat))
 
+  ## the one function the fit computes its spread probability with, so the prediction cannot differ from it
   data.table(pixelID = shortAnnDT$pixelID,
-             spreadProb = logisticAll(logisticPars, mat[, matching, drop = FALSE], covPars, lowerSpreadProb))
+             spreadProb = spreadProbFromCovs(mat[, matching, drop = FALSE], logisticPars, covPars, lowerSpreadProb,
+                                             climateRef = climateRef))
 }
 
 yearSpreadSDTxt <- "yearSpreadSD"
@@ -357,6 +363,24 @@ ledgerCovMinMax <- function(sa, i, default = NULL) {
   cmm
 }
 
+#' The climate reference of one ledger row
+#'
+#' A fit made with climate outside the `maxAsymptote` ceiling (`fireSense_spreadFit`'s `spreadLink = "climateOutside"`)
+#' stores the reference value of each climate covariate in the ledger row (`fireSenseUtils::spreadFitClimateRefTxt`),
+#' with `spreadLink` beside it. Prediction applies the same link with the same reference
+#' (`fireSenseUtils::spreadProbFromCovs()`).
+#'
+#' @param sa `sim$studyAreaWithSpreadParams`, one row per ELF.
+#' @param i integer; the row.
+#' @return A named list, or `NULL` for a row without a reference: a fit with climate under the ceiling, or a row
+#'   written before the link existed. Such a row predicts exactly as it always did.
+ledgerClimateRef <- function(sa, i) {
+  col <- fireSenseUtils::spreadFitClimateRefTxt
+  if (!col %in% names(sa)) return(NULL)
+  ref <- sa[[col]][[i]]
+  if (is.list(ref) && length(ref)) ref else NULL
+}
+
 #' The covariate centres of one ledger row
 #'
 #' A fit with an intercept centres its covariates and stores the centres in the ledger row

tests/testthat/test-climateOutside.R

@@ -0,0 +1,76 @@
+## A ledger row from a fit made with climate outside the maxAsymptote ceiling (fireSense_spreadFit's
+## spreadLink = "climateOutside") carries, in fireSenseUtils::spreadFitClimateRefTxt, the fit's reference value of
+## each climate covariate (on the rescaled scale), and "climateOutside" in fireSenseUtils::spreadLinkTxt. Prediction
+## applies the same link through fireSenseUtils::spreadProbFromCovs(), the function the fit used:
+## p0 = logistic of the other covariates, under the ceiling; p = plogis(qlogis(p0) + beta_clim * (climate - ref)),
+## uncapped. A row without a reference predicts as it always did.
+##
+## The toy covariates are those of test-spreadProb-values.R: MDC / 200 is the climate covariate, youngAge and
+## fuelA / 8 the rest; the toy parameters have hillSlope1 = 2, so p0 = 0.13 + 0.12 * plogis(2 * x_fuel), and the
+## climate shift is beta_clim * (MDC / 200 - ref) with no slope factor (the fit fixes hillSlope1 at 1).
+
+refCol <- fireSenseUtils::spreadFitClimateRefTxt
+linkCol <- fireSenseUtils::spreadLinkTxt
+cellsCO <- c(1, 2, 3, 4, 6, 7, 8)
+
+inputsCO <- function(ref = list(MDC = 0.4), link = "climateOutside") {
+  ins <- toyInputs()
+  if (!is.null(link)) ins$studyAreaWithSpreadParams[[linkCol]] <- link
+  ins$studyAreaWithSpreadParams[[refCol]] <- list(ref)
+  ins
+}
+
+## by hand, in cell order
+pHandCO <- function(ref = 0.4) {
+  cv <- toyCovariates(); o <- order(cv$pixelID)
+  mdc <- cv$MDC[o] / 200; ya <- cv$youngAge[o]; fuel <- cv$fuelA[o] / 8
+  p0 <- handLogistic3(-0.5 * ya + 1 * fuel)
+  plogis(qlogis(p0) + 1.5 * (mdc - ref))
+}
+
+test_that("a row with a climate reference predicts climate on the logit after the fuel ceiling, uncapped", {
+  v <- predVals(toyRun(inputsCO()))
+  expect_equal(v[cellsCO], pHandCO(0.4), tolerance = 1e-7)
+  expect_true(all(is.na(v[c(5, 9)])))
+  ## cell 4 (MDC 200, fuelA 8): p0 = 0.13 + 0.12 * plogis(2) = 0.2358, then + 1.5 * (1 - 0.4) on the logit
+  expect_equal(v[4], plogis(qlogis(handLogistic3(1)) + 0.9), tolerance = 1e-7)
+  expect_gt(v[4], 0.25)                      # above maxAsymptote: the ceiling no longer binds for climate
+  ## cell 1 (MDC 0, no fuel): p0 = 0.19, then 1.5 * (0 - 0.4) = -0.6 on the logit: below p0
+  expect_equal(v[1], plogis(qlogis(0.19) - 0.6), tolerance = 1e-7)
+  ## at climate = the reference the new link is p0: cell 3 has MDC 100 -> 0.5, so with ref 0.5 it is the
+  ## fuel-only logistic (fuelA 4 -> 0.5, x_fuel = 0.5)
+  v2 <- predVals(toyRun(inputsCO(ref = list(MDC = 0.5))))
+  expect_equal(v2[3], handLogistic3(0.5), tolerance = 1e-7)
+  expect_equal(v2[cellsCO], pHandCO(0.5), tolerance = 1e-7)
+})
+
+test_that("a row without a reference predicts exactly as before", {
+  base <- predVals(toyRun())
+  ## the column is there but empty: a fit with climate under the ceiling, written beside one without
+  expect_identical(predVals(toyRun(inputsCO(ref = list(), link = "ceiling"))), base)
+  expect_identical(predVals(toyRun(inputsCO(ref = NULL, link = NULL))), base)
+  ## and the reference changes the prediction only where climate differs from it
+  withRef <- predVals(toyRun(inputsCO(ref = list(MDC = 0.5))))
+  expect_false(isTRUE(all.equal(withRef, base)))
+  sa <- toyInputs()$studyAreaWithSpreadParams
+  expect_null(ledgerClimateRef(sa, 1L))
+  sa[[refCol]] <- list(list())
+  expect_null(ledgerClimateRef(sa, 1L))
+  sa[[refCol]] <- list(list(MDC = 0.4))
+  expect_identical(ledgerClimateRef(sa, 1L), list(MDC = 0.4))
+})
+
+test_that("with several ELFs each row applies its own reference, blended like the probabilities", {
+  ins <- multiInputs()
+  ## a climate effect for both ELFs; MDC is 100 everywhere, i.e. 0.5 rescaled
+  ins$studyAreaWithSpreadParams$params <- lapply(ins$studyAreaWithSpreadParams$params, function(p) { p$MDC <- 1; p })
+  ins$studyAreaWithSpreadParams[[linkCol]] <- "climateOutside"
+  ins$studyAreaWithSpreadParams[[refCol]] <- list(list(MDC = 0.5), list(MDC = 0.25))
+  v <- predVals(toyRun(ins, params = list(ELFblendWidth = 20000)))
+  pA <- 0.19                                     # at its reference: the fuel-only value (all other coefficients 0)
+  pB <- plogis(qlogis(0.20) + 1 * (0.5 - 0.25))  # 0.25 above its reference on the logit
+  rA <- c(1, 1, 1, 1, 1, 0.75, 0.5, 0.25, 0, 0); rB <- rev(rA)
+  expect_equal(v, (rA * pA + rB * pB) / (rA + rB), tolerance = 1e-9)
+  expect_equal(v[1], pA, tolerance = 1e-9)
+  expect_equal(v[10], pB, tolerance = 1e-9)
+})