Predict with the intercept and centres of a ledger row

NEWS.md

@@ -1,5 +1,7 @@
 # fireSense_spreadPredict (development version)
 
+- A ledger row from a fit with an intercept (a parameter named `fireSenseUtils::spreadInterceptTxt`, from `fireSense_dataPrepFit`'s `spreadIntercept`) is predicted with it: `spreadProbOneELF()` adds a column of 1s for it after rescaling and centring, and subtracts the row's stored covariate centres (`fireSenseUtils::spreadFitCovCentreTxt`, `ledgerCovCentre()`) from the rescaled covariates, as the fit did. With several ELFs each row uses its own. A row without them, such as every row written so far, predicts exactly as before. Needs the `fireSenseUtils` change in PredictiveEcology/fireSenseUtils#131 (floor to be set once it has a version).
+
 - Predictions no longer average parameter sets (`spreadProbOneELF()` took `rowMeans()` of one map per row, and `yearSpreadSDOf()` the mean `yearSpreadSD`, which broke the pairing of each set's `yearSpreadSD` with its own coefficients). New parameter `.rep` (integer, default 1): replicate `.rep` uses the one whole parameter set `((rep - 1) %% number of sets) + 1` of each ELF, and `fireSense_SpreadSD` is that set's `yearSpreadSD`. `.rep` is a SpaDES-aware parameter: `SpaDES.project::setupProject()` sets `.globals$.rep` from the experiment's `.rep` (PredictiveEcology/SpaDES.project#190), so nothing needs wiring by hand.
 
 - reqdPkgs now lists `data.table`, which the module calls (`copy`, `data.table`, `setDT`) but did not list. Version 1.1.2.

fireSense_spreadPredict.R

@@ -151,7 +151,8 @@ spreadPredictRun <- function(sim) {
     params1 <- paramSetForRep(sa$params[[1]], P(sim)$.rep)
     pred <- spreadProbOneELF(covs, params = params1, covMinMax = sim$covMinMax_spread,
                              formula = sim$fireSense_spreadFormula, yr = time(sim),
-                             maxFireSpread = P(sim)$maxFireSpread, lowerSpreadProb = P(sim)$lowerSpreadProb)
+                             maxFireSpread = P(sim)$maxFireSpread, lowerSpreadProb = P(sim)$lowerSpreadProb,
+                             covCentre = ledgerCovCentre(sa, 1L))
     sim$fireSense_SpreadSD <- yearSpreadSDOf(params1)
   } else {
     ids <- as.character(sa[[fireSenseUtils::polygonIDTxt]])
@@ -171,7 +172,8 @@ spreadPredictRun <- function(sim) {
       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)
+                            lowerSpreadProb = P(sim)$lowerSpreadProb, byParams = TRUE,
+                            covCentre = ledgerCovCentre(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)
@@ -240,9 +242,13 @@ ELFblendWeights <- function(elfRas, template, pixelID, ids, width) {
 #'   `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) {
+                             byParams = FALSE, covCentre = NULL) {
   moduleName <- "fireSense_spreadPredict"
   covs <- copy(covs)
   ## the per-year random effect is not a covariate coefficient: fireSense_burn applies it (fireSense_SpreadSD)
@@ -265,7 +271,9 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
 
   ## every fitted coefficient needs a covariate: without one the logistic gets no covariates and
   ## fails later in rowMeans(), with no hint of the cause
-  fitted <- setdiff(names(params), unlist(fireSenseUtils::logisticParamNames))
+  ## the intercept, if the fit has one, is a coefficient with no covariate: it multiplies a column of 1s
+  hasIntercept <- fireSenseUtils::spreadInterceptTxt %in% names(params)
+  fitted <- setdiff(names(params), c(unlist(fireSenseUtils::logisticParamNames), fireSenseUtils::spreadInterceptTxt))
   noCov <- setdiff(fitted, names(covs))
   if (length(noCov)) {
     stop(
@@ -278,7 +286,7 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
 
   ## the covariates this ELF was fitted with
   needed <- if (isTRUE(byParams)) {
-    setdiff(names(params), unlist(fireSenseUtils::logisticParamNames))
+    fitted
   } else {
     terms <- delete.response(terms.formula(as.formula(formula)))
     all.vars(reformulate(attr(terms, "term.labels"), intercept = attr(terms, "intercept")))
@@ -297,6 +305,8 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
   # integers x 1000, the form `spreadProbFromIntegerCovs` expects
   shortAnnDTx1000 <- toX1000(list(covs))[[1]] |> setDT()
   colsToUse <- setdiff(names(covs), "pixelID")
+  ## what the coefficients multiply: the intercept first, as in the fit (fireSenseUtils::spreadDesignCols())
+  designCols <- c(if (hasIntercept) fireSenseUtils::spreadInterceptTxt, colsToUse)
 
   ## youngAge is mutually exclusive with every other non-climate covariate, exactly as in the fit
   ## (fireSense_spreadFit::spreadFitPrep()): wherever youngAge is non-zero, fuel biomass, nfLCC_*
@@ -314,16 +324,17 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
                               yr = yr,
                               covMinMax = covMinMax,
                               mutuallyExclusive = mutuallyExclusive,
-                              colsToUse = colsToUse,
+                              colsToUse = designCols,
                               doAssertions = FALSE,
                               logisticPars = params,
-                              maxFireSpread = maxFireSpread
+                              maxFireSpread = maxFireSpread,
+                              covCentre = covCentre
                               )
 
-  mat <- as.matrix(shortAnnDT[, ..colsToUse])
+  mat <- as.matrix(shortAnnDT[, ..designCols])
 
   par <- params[1L, ] |> as.vector() |> unlist()
-  covPars <- par[intersect(names(par), colsToUse)]
+  covPars <- par[intersect(names(par), designCols)]
   logisticPars <- par[setdiff(names(par), names(covPars))]
   # Make sure the order is correct in the matrix
   matching <- intersect(names(covPars), colnames(mat))
@@ -339,6 +350,22 @@ spreadProbOneELF <- function(covs, params, covMinMax, formula, yr, maxFireSpread
 
 yearSpreadSDTxt <- "yearSpreadSD"
 
+#' The covariate centres of one ledger row
+#'
+#' A fit with an intercept centres its covariates and stores the centres in the ledger row
+#' (`fireSenseUtils::spreadFitCovCentreTxt`). Prediction subtracts the same values.
+#'
+#' @param sa `sim$studyAreaWithSpreadParams`, one row per ELF.
+#' @param i integer; the row.
+#' @return A named list, or `NULL` for a row without centres: a fit without an intercept, or a row written
+#'   before they were stored. Such a row predicts exactly as it always did.
+ledgerCovCentre <- function(sa, i) {
+  col <- fireSenseUtils::spreadFitCovCentreTxt
+  if (!col %in% names(sa)) return(NULL)
+  centre <- sa[[col]][[i]]
+  if (is.list(centre) && length(centre)) centre else NULL
+}
+
 #' The parameter set a replicate uses
 #'
 #' Replicate `rep` uses row `((rep - 1) %% NROW(params)) + 1` of an ELF's fitted parameter sets, whole: its

tests/testthat/test-spreadIntercept.R

@@ -0,0 +1,78 @@
+## A ledger row from a fit with an intercept holds a coefficient named fireSenseUtils::spreadInterceptTxt
+## and, in the column fireSenseUtils::spreadFitCovCentreTxt, the covariate means the fit subtracted from the
+## rescaled covariates. Prediction applies both. A row without them predicts as it always did.
+##
+## x = intercept + sum(beta * (rescaled covariate - centre)); the toy covariates are those of
+## test-spreadProb-values.R: MDC / 200, youngAge, fuelA / 8.
+
+icptName <- fireSenseUtils::spreadInterceptTxt
+centreCol <- fireSenseUtils::spreadFitCovCentreTxt
+centreI <- list(MDC = 0.45, youngAge = 0.2, fuelA = 0.3)
+
+paramsI <- function(intercept = 0.4) {
+  data.frame(maxAsymptote = 0.25, hillSlope1 = 2, inflectionPoint1 = 1, `(Intercept)` = intercept,
+             MDC = 1.5, youngAge = -0.5, fuelA = 1, check.names = FALSE)
+}
+inputsI <- function(intercept = TRUE, centre = TRUE) {
+  ins <- toyInputs(params = if (intercept) paramsI() else toyParams(),
+                   formula = if (intercept) "~ 1 + MDC + youngAge + fuelA" else "~ MDC + youngAge + fuelA - 1")
+  if (centre) ins$studyAreaWithSpreadParams[[centreCol]] <- list(centreI)
+  ins
+}
+
+## the linear predictor at each toy cell, by hand
+xHand <- function(intercept = 0.4, centre = list(MDC = 0, youngAge = 0, fuelA = 0)) {
+  cv <- toyCovariates()
+  o <- order(cv$pixelID)
+  mdc <- cv$MDC[o] / 200; ya <- cv$youngAge[o]; fuel <- cv$fuelA[o] / 8
+  ## in pixelID order, which is cell order for the cells that have one
+  intercept + 1.5 * (mdc - centre$MDC) - 0.5 * (ya - centre$youngAge) + 1 * (fuel - centre$fuelA)
+}
+cells <- c(1, 2, 3, 4, 6, 7, 8)
+
+test_that("a row with an intercept and centres predicts the logistic of intercept + centred covariates", {
+  v <- predVals(toyRun(inputsI()))
+  expect_equal(v[cells], handLogistic3(xHand(0.4, centreI)), tolerance = 1e-7)
+  expect_true(all(is.na(v[c(5, 9)])))
+})
+
+test_that("an intercept without centres, and centres without an intercept, each apply on their own", {
+  v1 <- predVals(toyRun(inputsI(centre = FALSE)))
+  expect_equal(v1[cells], handLogistic3(xHand(0.4)), tolerance = 1e-7)
+  v2 <- predVals(toyRun(inputsI(intercept = FALSE)))
+  expect_equal(v2[cells], handLogistic3(xHand(0, centreI)), tolerance = 1e-7)
+})
+
+test_that("an intercept with a centre is the same model as the intercept i - sum(beta * centre) with none", {
+  shift <- 0.4 - (1.5 * centreI$MDC - 0.5 * centreI$youngAge + 1 * centreI$fuelA)
+  a <- predVals(toyRun(inputsI()))
+  ins <- toyInputs(params = paramsI(shift), formula = "~ 1 + MDC + youngAge + fuelA")
+  expect_equal(predVals(toyRun(ins)), a, tolerance = 1e-9)
+})
+
+test_that("a row without an intercept or centres predicts exactly as before", {
+  v <- predVals(toyRun(inputsI(intercept = FALSE, centre = FALSE)))
+  expect_equal(v[cells], handLogistic3(xHand(0)), tolerance = 1e-7)
+  expect_identical(v, predVals(toyRun()))
+  expect_null(ledgerCovCentre(toyInputs()$studyAreaWithSpreadParams, 1L))
+})
+
+test_that("an old ledger row's missing centre is NULL, also beside a row that has one", {
+  sa <- data.frame(polygonID = c("A", "B"))
+  sa[[centreCol]] <- list(centreI, NULL)         # as rbindlist(fill = TRUE) leaves an older row
+  expect_identical(ledgerCovCentre(sa, 1L), centreI)
+  expect_null(ledgerCovCentre(sa, 2L))
+  expect_null(ledgerCovCentre(data.frame(polygonID = "A"), 1L))
+})
+
+test_that("with several ELFs each uses its own intercept and centre; an older ELF's row predicts as before", {
+  ins <- multiInputs()
+  ## ELF A: intercept and centre (MDC, fuelA); ELF B: none. All coefficients 0 but the intercept.
+  ins$studyAreaWithSpreadParams$params[[1]][[icptName]] <- 0.5
+  ins$studyAreaWithSpreadParams$params[[1]] <- ins$studyAreaWithSpreadParams$params[[1]][c(1:3, 6, 4:5)]
+  ins$studyAreaWithSpreadParams[[centreCol]] <- list(list(MDC = 0.1, fuelA = 0.2), NULL)
+  v <- predVals(toyRun(ins, params = list(ELFblendWidth = 6000)))
+  ## coefficients are 0, so the centre does not matter here: A is lower + (0.25 - 0.13) * plogis(2 * 0.5)
+  expect_equal(v[1], 0.13 + 0.12 * stats::plogis(2 * 0.5), tolerance = 1e-9)
+  expect_equal(v[10], 0.20, tolerance = 1e-9)    # B: no intercept, as before
+})