New parameter spreadIntercept: a spread formula with an intercept

NEWS.md

@@ -1,5 +1,6 @@
 # fireSense_dataPrepFit (development version)
 
+- New parameter `spreadIntercept` (default `FALSE`): with `TRUE` the spread formula is `~ 1 + ...`, so `fireSense_spreadFit` fits a free intercept with centred covariates (`fireSenseUtils::spreadInterceptTxt`) and the coefficients describe variation about the centre. Covariates are rescaled to 0-1, so without an intercept the level of the linear predictor is set by the coefficients alone and the climate coefficient trades off against the fuel and non-forest ones (ELF 13.1: the other coefficients explain a median 73% of CMD's variance in a final population). `FALSE` keeps the formula string `~ 0 + ...` byte for byte. Needs the `fireSenseUtils` change in PredictiveEcology/fireSenseUtils#131 (floor to be set once it has a version).
 - New parameter `minCovariateProp` (default `0.05`): a land-cover class covering less than this share of the ELF's flammable pixels (most recent data year) gets no spread coefficient of its own. A rare non-forest class joins the non-forest group nearest its burn rate; with less treed wetland than this there is no `treedWetland_agb` and its tree AGB stays in `dom_agb_*`/`sec_agb_*`. Treed wetland is under 5% in 7 of the 10 fitted ELFs. Needs fireSenseUtils >= 0.2.3.9081 (PredictiveEcology/fireSenseUtils#129). Fits and cached `dataPrepBuild`/`prepSpreadFitData` results of ELFs whose covariates change are invalid.
 - Fixed: with one fitted ELF in the SpreadFit ledger, `Init()` now sets `nonForestedLCCGroups` and `missingLCCgroup` to that
   ELF's groups, so this run's ignition and escape fits use the same non-forest columns (`nfLCC_*`) that

R/spreadFormula.R

@@ -0,0 +1,12 @@
+#' The spread formula
+#'
+#' `~ 0 + RHS`, which has no intercept, or `~ 1 + RHS` with `intercept = TRUE`
+#' (`fireSenseUtils::spreadInterceptTxt` is its coefficient in `fireSense_spreadFit`). Without an
+#' intercept the string is what it always was, so cache keys that hold it do not move.
+#'
+#' @param RHS character; the covariate terms joined with `" + "`.
+#' @param intercept logical; `TRUE` for a formula with an intercept.
+#' @return a character string.
+spreadFormula <- function(RHS, intercept = FALSE) {
+  paste0("~ ", if (isTRUE(intercept)) "1" else "0", " + ", RHS)
+}

fireSense_dataPrepFit.R

@@ -113,6 +113,14 @@ defineModule(sim, list(
                     NA, NA, paste("Name of the ledger file in `spreadFitGoogleDriveFolder`: study area polygons with",
                                   "their fitted SpreadFit parameters. `\"latest\"` (the default) takes each polygon's",
                                   "fit from the most recent ledger file that has it (`fireSenseUtils::latestSpreadFits()`).")),
+    defineParameter("spreadIntercept", "logical", FALSE, NA, NA,
+                    paste("If `TRUE`, the spread formula has an intercept (`~ 1 + ...`) and `fireSense_spreadFit` fits",
+                          "it, with the covariates centred, so the coefficients describe variation and the intercept",
+                          "the level. Covariates are rescaled to 0-1 and so are all >= 0: without an intercept the",
+                          "level of the linear predictor is set by the coefficients alone, and a climate coefficient trades off",
+                          "against the fuel and non-forest ones (ELF 13.1). `FALSE` (default) keeps `~ 0 + ...`.",
+                          "Ignored if `sim$fireSense_spreadFormula` is supplied. Needs a `fireSense_spreadFit` and a",
+                          "`fireSense_spreadPredict` that know the intercept.")),
     defineParameter("targetFuelClasses", "numeric", 5, 1, 7,
                     "the target number of unique fuel classes when using semi-automated approach"),
     defineParameter("useRasterizedFireForSpread", "logical", FALSE, NA, NA,
@@ -939,7 +947,7 @@ prepare_SpreadFit <- function(sim) {
   sim$fireSense_nonAnnualSpreadFitCovariates <- nonAnnuals
 
   if (is.null(sim$fireSense_spreadFormula)) {
-    sim$fireSense_spreadFormula <- paste0("~ 0 + ", RHS)
+    sim$fireSense_spreadFormula <- spreadFormula(RHS, Par$spreadIntercept)
   }
 
   return(invisible(sim))

tests/testthat/test-metadata.R

@@ -97,7 +97,7 @@ test_that("parameters are the expected names", {
            "forestedLCC", "fuelClassCol", "fuelCovariates", "heldOutFold", "igAggFactor",
            "minBufferSize", "minCovariateProp", "nonflammableLCC",
            "nonForestCanBeYoungAge", "scanfiVersion", "sppEquivCol", "spreadFitFilename",
-           "spreadFitGoogleDriveFolder", "targetFuelClasses",
+           "spreadFitGoogleDriveFolder", "spreadIntercept", "targetFuelClasses",
            "useRasterizedFireForSpread", "whichModulesToPrepare"))
   )
 })

tests/testthat/test-spreadFormulaYoungAgeOnce.R

@@ -23,6 +23,7 @@ test_that("the spread formula lists a veg-table youngAge column once, not twice"
   env$sim <- new.env()
   env$sim$climateVariablesForFire <- list(spread = "CMD")
   env$sim$fireSense_spreadFormula <- NULL
+  env$Par <- list(spreadIntercept = FALSE) # the module's `Par`, as it is with the default
   ## the shape fireSenseVegData has after joinFireBuffersToVeg() and setnames(..., "buffer",
   ## "burned"): a youngAge column that is not all-zero, as in ELFs 14.3/14.4
   env$fireSenseVegData <- data.table::data.table(

tests/testthat/test-spreadIntercept.R

@@ -0,0 +1,35 @@
+## `spreadIntercept` (default FALSE): the spread formula has an intercept, `~ 1 + RHS`, which
+## fireSense_spreadFit fits with centred covariates. Off, the formula string is byte-identical to what it
+## has always been, so every cache key that holds it is unchanged.
+
+test_that("spreadFormula() gives '~ 0 + RHS' by default and '~ 1 + RHS' with an intercept", {
+  rhs <- "CMD + youngAge + fuelA"
+  expect_identical(spreadFormula(rhs), "~ 0 + CMD + youngAge + fuelA")
+  expect_identical(spreadFormula(rhs, intercept = FALSE), "~ 0 + CMD + youngAge + fuelA")
+  expect_identical(spreadFormula(rhs, intercept = TRUE), "~ 1 + CMD + youngAge + fuelA")
+  ## the default string is the one the module built before the parameter existed
+  expect_identical(spreadFormula(rhs), paste0("~ 0 + ", rhs))
+  ## and the intercept is what terms() sees
+  expect_equal(attr(terms(as.formula(spreadFormula(rhs, TRUE))), "intercept"), 1L)
+  expect_equal(attr(terms(as.formula(spreadFormula(rhs))), "intercept"), 0L)
+  expect_identical(fireSenseUtils::spreadDesignCols(spreadFormula(rhs, TRUE)),
+                   c(fireSenseUtils::spreadInterceptTxt, "CMD", "youngAge", "fuelA"))
+})
+
+test_that("the module defines spreadIntercept, a logical that is FALSE by default", {
+  def <- Filter(function(x) is.call(x) && identical(x[[1]], as.name("defineModule")),
+                parse(testthat::test_path("..", "..", "fireSense_dataPrepFit.R"), keep.source = FALSE))
+  params <- eval(as.list(def[[1]][[3]])[["parameters"]], envir = asNamespace("SpaDES.core"))
+  i <- match("spreadIntercept", params$paramName)
+  expect_false(is.na(i))
+  expect_identical(params$paramClass[[i]], "logical")
+  expect_identical(params$default[[i]], FALSE)
+})
+
+test_that("prepare_SpreadFit() builds the formula with spreadFormula() from the spreadIntercept parameter", {
+  src <- paste(deparse(body(prepare_SpreadFit)), collapse = "\n")
+  expect_match(src, "spreadFormula\\(RHS, Par\\$spreadIntercept\\)")
+  expect_no_match(src, '"~ 0 \\+ "', fixed = FALSE)
+  ## the formula is only built when none was supplied
+  expect_match(src, "is.null\\(sim\\$fireSense_spreadFormula\\)")
+})