Predict with the fitted terms, not the spread formula's

NEWS.md

@@ -1,5 +1,7 @@
 # fireSense_spreadPredict (development version)
 
+- 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 1.2.0
 
 This release lets fire spread be forecast across study areas that contain several fire regions, each with its own fitted model, with predictions blended smoothly across the boundaries between them. Each simulation replicate now uses one complete set of fitted values instead of an average of all of them, so the spread of results across replicates reflects the uncertainty in the fit. The fitted year-to-year variation in fire spread is passed on to the fire simulation.

fireSense_spreadPredict.R

@@ -152,7 +152,7 @@ spreadPredictRun <- function(sim) {
   if (NROW(sa) == 1L) {
     params1 <- paramSetForRep(sa$params[[1]], P(sim)$.rep)
     pred <- spreadProbOneELF(covs, params = params1, covMinMax = sim$covMinMax_spread,
-                             formula = sim$fireSense_spreadFormula, yr = time(sim),
+                             yr = time(sim),
                              maxFireSpread = P(sim)$maxFireSpread, lowerSpreadProb = P(sim)$lowerSpreadProb,
                              covCentre = ledgerCovCentre(sa, 1L))
     sim$fireSense_SpreadSD <- yearSpreadSDOf(params1)
@@ -173,8 +173,8 @@ spreadPredictRun <- function(sim) {
       if (!length(these)) next
       parsI <- paramSetForRep(sa$params[[i]], P(sim)$.rep)
       p <- spreadProbOneELF(covs[these], params = parsI, covMinMax = sa$covMinMax_spread[[i]],
-                            formula = NULL, yr = time(sim), maxFireSpread = P(sim)$maxFireSpread,
-                            lowerSpreadProb = P(sim)$lowerSpreadProb, byParams = TRUE,
+                            yr = time(sim), maxFireSpread = P(sim)$maxFireSpread,
+                            lowerSpreadProb = P(sim)$lowerSpreadProb,
                             covCentre = ledgerCovCentre(sa, i))
       rows <- these[match(p$pixelID, covs$pixelID[these])]
       acc[rows] <- acc[rows] + w[rows, i] * p$spreadProb
@@ -236,21 +236,19 @@ ELFblendWeights <- function(elfRas, template, pixelID, ids, width) {
 #' Spread probability for the pixels of one ELF
 #'
 #' @param covs `data.table` of this ELF's pixels: `pixelID` and covariates (as `fireSense_dataPrepPredict`
-#'   makes them; fuel biomass logged).
+#'   makes them; fuel biomass logged). It can hold more than the fit used (every ELF's covariates, or a term the
+#'   spread formula has but the fit does not); only the covariates named in `params` are used.
 #' @param params one-row `data.frame`: the ELF's parameter set for this replicate (`paramSetForRep()`).
 #' @param covMinMax the ELF's `covMinMax_spread`.
-#' @param formula the spread formula, to check the covariates are all there.
-#' @param byParams logical; `TRUE` (several ELFs) ignores `formula`: the covariates are those named in
-#'   `params`, and the table is cut to them, since it holds every ELF's covariates.
 #' @param yr,maxFireSpread,lowerSpreadProb as for `fireSenseUtils::spreadProbFromIntegerCovs()` and
 #'   `fireSenseUtils::logisticAll()`.
 #' @param covCentre the ELF's `covCentre_spread` (`ledgerCovCentre()`): the covariate means the fit subtracted
 #'   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.
 #' @return `data.table` with `pixelID` and `spreadProb`, from the one parameter set in `params`.
-spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread, lowerSpreadProb,
-                             byParams = FALSE, covCentre = NULL) {
+spreadProbOneELF <- function(covs, params, covMinMax, yr, maxFireSpread, lowerSpreadProb,
+                             covCentre = NULL) {
   moduleName <- "fireSense_spreadPredict"
   covs <- copy(covs)
   ## the per-year random effect is not a covariate coefficient: fireSense_burn applies it (fireSense_SpreadSD)
@@ -286,23 +284,9 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
     )
   }
 
-  ## the covariates this ELF was fitted with
-  needed <- if (isTRUE(byParams)) {
-    fitted
-  } else {
-    terms <- delete.response(terms.formula(as.formula(formula)))
-    all.vars(reformulate(attr(terms, "term.labels"), intercept = attr(terms, "intercept")))
-  }
-  missing <- !needed %in% names(covs)
-  if (s <- sum(missing)) {
-    stop(
-      moduleName, "> '", needed[missing][1L], "'",
-      if (s > 1) paste0(" (and ", s - 1L, " other", if (s > 2) "s", ")"),
-      " not found in data objects."
-    )
-  }
-  ## with several ELFs the table holds every ELF's covariates; use this ELF's only
-  if (isTRUE(byParams)) covs <- covs[, c("pixelID", needed), with = FALSE]
+  ## the covariates this ELF was fitted with, not the spread formula's terms: they can differ (a hindcast
+  ## rebuilds the formula from another year's land cover), and `covs` can hold more than this ELF's
+  covs <- covs[, c("pixelID", fitted), with = FALSE]
 
   # integers x 1000, the form `spreadProbFromIntegerCovs` expects
   shortAnnDTx1000 <- toX1000(list(covs))[[1]] |> setDT()
@@ -340,11 +324,6 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
   logisticPars <- par[setdiff(names(par), names(covPars))]
   # Make sure the order is correct in the matrix
   matching <- intersect(names(covPars), colnames(mat))
-  missingCovs <- setdiff(colnames(mat), names(covPars))
-  if (length(missingCovs))
-    warning("There are covariates in the sim$fireSense_SpreadCovariates: \n",
-            paste0(missingCovs, collapse = ", "),
-            "\n...that are not in the sim$studyAreaWithSpreadParams")
 
   data.table(pixelID = shortAnnDT$pixelID,
              spreadProb = logisticAll(logisticPars, mat[, matching, drop = FALSE], covPars, lowerSpreadProb))

fireSense_spreadPredict.Rmd

@@ -53,7 +53,7 @@ Each year, predicts a raster of fire spread probabilities from the parameters fi
 
 ### Module inputs and parameters
 
-Two objects from *fireSense_spreadFit* are read from the `simList` though they are not declared as inputs: `studyAreaWithSpreadParams` (the fitted parameters) and `fireSense_spreadFormula` (every term must be a column of `fireSense_SpreadCovariates`).
+`studyAreaWithSpreadParams` (the fitted parameters, from *fireSense_spreadFit*) is read from the `simList` though it is not declared as an input. Prediction uses the covariates the fit has coefficients for, not the terms of `fireSense_spreadFormula`; each must be a column of `fireSense_SpreadCovariates`.
 The module stops if `studyAreaWithSpreadParams` has no parameters.
 `maxFireSpread` must have the same value in every module that defines it.
 

tests/testthat/test-failures.R

@@ -1,17 +1,5 @@
 ## Failure paths of the run event, with the message each one gives.
 
-test_that("a formula term that is not a covariate column stops, naming the term", {
-  ins <- toyInputs(formula = "~ MDC + youngAge + fuelA + fuelZ - 1")
-  expect_error(toyRun(ins), "'fuelZ' not found in data objects")
-})
-
-test_that("several missing terms are counted in the message", {
-  ins <- toyInputs(formula = "~ MDC + fuelY + fuelZ + fuelW - 1")
-  expect_error(toyRun(ins), "'fuelY' \\(and 2 others\\) not found in data objects")
-  ins <- toyInputs(formula = "~ MDC + fuelY + fuelZ - 1")
-  expect_error(toyRun(ins), "'fuelY' \\(and 1 other\\) not found in data objects")
-})
-
 test_that("a response on the left of the formula is ignored", {
   ins <- toyInputs(formula = "fires ~ MDC + youngAge + fuelA - 1")
   expect_equal(predVals(toyRun(ins))[1], 0.19, tolerance = 1e-7)
@@ -28,12 +16,6 @@ test_that("no fitted parameters stops with an explanation that names the run", {
   expect_error(toyRun(ins), "holds no fitted spread parameters for this run \\(toyRunName\\)")
 })
 
-test_that("a missing spread formula stops", {
-  ins <- toyInputs()
-  ins$fireSense_spreadFormula <- NULL
-  expect_error(toyRun(ins), "argument is not a valid model")
-})
-
 test_that("a fitted coefficient with no covariate stops at once, naming it and the covariates available", {
   ## a fit made with non-forest groups nfLCC_100_60 and nfLCC_40_50_80; the covariates were built
   ## with a single group `nf`, so neither coefficient has a column

tests/testthat/test-fittedTerms.R

@@ -0,0 +1,20 @@
+## A hindcast built its spread formula with a term (treedWetland_agb) the ledger fit does not have; the
+## covariates match the fit. Prediction must use the fitted terms, not the formula's.
+
+test_that("a formula term the fit does not have is ignored: the fitted terms predict", {
+  ref <- predVals(toyRun(toyInputs()))
+  ins <- toyInputs(formula = "~ MDC + youngAge + fuelA + fuelZ - 1")
+  expect_identical(predVals(toyRun(ins)), ref)
+})
+
+test_that("a formula missing a fitted term still predicts with every fitted term", {
+  ref <- predVals(toyRun(toyInputs()))
+  ins <- toyInputs(formula = "~ MDC - 1")
+  expect_identical(predVals(toyRun(ins)), ref)
+})
+
+test_that("covariate columns beyond the fitted terms are not used", {
+  covs <- toyCovariates()
+  data.table::set(covs, j = "treedWetland_agb", value = 1)
+  expect_identical(predVals(toyRun(toyInputs(covs = covs))), predVals(toyRun(toyInputs())))
+})

tests/testthat/test-spreadProb-values.R

@@ -140,12 +140,12 @@ test_that("maxFireSpread does not change the prediction", {
   expect_equal(a[4], 0.2491969, tolerance = 1e-6) # above 0.20: not capped
 })
 
-test_that("a covariate with no fitted coefficient is dropped with a warning", {
+test_that("a covariate with no fitted coefficient is not used", {
   covs <- toyCovariates()
   data.table::set(covs, NULL, "fuelB", 5)
   ins <- toyInputs(covs = covs)
   ins$covMinMax_spread$fuelB <- c(0, 10)
-  expect_warning(sim <- toyRun(ins), "fuelB\\s+\\.\\.\\.that are not in the sim\\$studyAreaWithSpreadParams")
+  expect_no_warning(sim <- toyRun(ins))
   ## fuelB contributes nothing
   expect_equal(predVals(sim)[c(1, 3)], c(0.19, 0.2408970), tolerance = 1e-6)
 })