From 8fca0a97e155dbc98487067397b4d167500d6926 Mon Sep 17 00:00:00 2001 From: Eliot McIntire Date: Mon, 5 Oct 2026 09:49:44 -0700 Subject: [PATCH] Predict with the intercept and centres of a ledger row A row with an intercept and stored covariate centres is predicted with both; a row without them predicts as before. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01CwcjqqK59FmTJscyi7xUqv --- NEWS.md | 2 + fireSense_spreadPredict.R | 45 ++++++++++++---- tests/testthat/test-spreadIntercept.R | 78 +++++++++++++++++++++++++++ 3 files changed, 116 insertions(+), 9 deletions(-) create mode 100644 tests/testthat/test-spreadIntercept.R diff --git a/NEWS.md b/NEWS.md index b33e4ce..f6ef5cf 100644 --- a/NEWS.md +++ b/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. diff --git a/fireSense_spreadPredict.R b/fireSense_spreadPredict.R index 29b0330..83bb4e7 100644 --- a/fireSense_spreadPredict.R +++ b/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 diff --git a/tests/testthat/test-spreadIntercept.R b/tests/testthat/test-spreadIntercept.R new file mode 100644 index 0000000..bb287f3 --- /dev/null +++ b/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 +})