From 7646328c4f63aa698f6d3e946cb70e7fe36cf198 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sat, 9 May 2026 13:50:43 -0700 Subject: [PATCH 1/7] Initialize PWF baseline for #20 + #43 Combined planning baseline for cd_compare defaults (#20) and window-vs-window p-value (#43). Issues are explicitly coupled in #20's issue body ('Tightly coupled with #43; both should land together'), so single branch / single PR. Co-Authored-By: Claude Opus 4.7 (1M context) --- planning/active/findings.md | 115 +++++++++++++++++++++++++++++++++++ planning/active/progress.md | 11 ++++ planning/active/task_plan.md | 69 +++++++++++++++++++++ 3 files changed, 195 insertions(+) create mode 100644 planning/active/findings.md create mode 100644 planning/active/progress.md create mode 100644 planning/active/task_plan.md diff --git a/planning/active/findings.md b/planning/active/findings.md new file mode 100644 index 0000000..1242571 --- /dev/null +++ b/planning/active/findings.md @@ -0,0 +1,115 @@ +# Findings — cd_compare defaults + window-vs-window p-value (#20 + #43) + +## Issue #20 — context + +`cd_compare()` requires the user to specify both time windows +explicitly. There's no sensible default that helps a first-time +user understand what climate departure looks like for their AOI. + +**Status (post-v0.2.1):** vignettes have settled on a de-facto +framing — recent **2015-2025** vs WMO-style **1951-1980**, paired +with both 75-yr and 45-yr Theil-Sen trends. Both `peace-fwcp` and +`kootenay-lake` use this consistently. The methodology lit review +from #54 (archived at +`planning/archive/2026-05-issue-53-snow-lit-review/findings.md`) +backs the choice — Mann-Kendall + Theil-Sen with no prewhitening +per Yue & Wang 2002 is correct for our 76-yr series with strong +trends. + +**Remaining scope:** wire those windows in as `cd_compare()` +defaults so a first-time user gets a sensible answer without +having to supply both windows. Document the rationale (cumulative +impact framing — "2 °C warmer" punchline rather than per-year +rate) in roxygen + a brief note in the vignettes. + +Tightly coupled with #43. + +## Issue #43 — context + +`cd_compare()` returns a window-vs-window difference (mean_diff / +pct_change) but no measure of statistical confidence. Readers of +the vignette want to know "is this window-shift real or noise?" + +The Mann-Kendall p-value from `cd_trend()` is currently used as a +proxy in the FWCP Peace vignette but answers a different question +("is there a monotonic trend?") — not "do the two windows differ?" +Two windows can differ significantly without a monotonic trend +(step change, U-shape), and a non-significant trend doesn't always +mean the windows are statistically indistinguishable. + +**Proposal:** add `test = c("t", "wilcox", NULL)` to `cd_compare()`. +Welch t-test (default) or Mann-Whitney U. Output gains `p_value` +column (NA when `test = NULL`). Small samples (< 8 in either +window) → return NA with a warning. Document the independence +assumption. + +## Codebase exploration + +### `R/cd_compare.R` (current) + +```r +cd_compare <- function(x, window_a, window_b, method = "mean_diff") { + method <- match.arg(method, c("mean_diff", "pct_change")) + mean_a <- x |> + dplyr::filter(.data$year %in% window_a) |> + dplyr::summarise(mean_a = mean(.data$value, na.rm = TRUE), .by = c("variable", "period")) + mean_b <- x |> + dplyr::filter(.data$year %in% window_b) |> + dplyr::summarise(mean_b = mean(.data$value, na.rm = TRUE), .by = c("variable", "period")) + n_a <- length(unique(x$year[x$year %in% window_a])) + n_b <- length(unique(x$year[x$year %in% window_b])) + if (n_a < 2) warning("window_a has fewer than 2 years of data", call. = FALSE) + if (n_b < 2) warning("window_b has fewer than 2 years of data", call. = FALSE) + out <- dplyr::left_join(mean_a, mean_b, by = c("variable", "period")) + out$difference <- switch(method, + mean_diff = out$mean_a - out$mean_b, + pct_change = (out$mean_a - out$mean_b) / abs(out$mean_b) * 100 + ) + out$method <- method + out +} +``` + +Returns 6 cols: `variable, period, mean_a, mean_b, difference, method`. + +### Existing tests (`tests/testthat/test-cd_compare.R`) + +4 tests covering mean_diff, pct_change, small-window warning, and +multi-variable handling. All use `1956:1960 / 1951:1955` explicit +windows on a 1951:1960 toy series — independent of the new +defaults. + +### Vignette callsites + +- `peace-fwcp.Rmd:534` — `cmp <- cd::cd_compare(ts, 2015:2025, 1951:1980)` +- `peace-fwcp.Rmd:656` — same in per-ecoregion loop +- `kootenay-lake.Rmd:561` — same +- `kootenay-lake.Rmd:688` — same in per-ecoregion loop +- `R/cd_plot_comparison.R:18` — example; same windows + +The `cd_plot_comparison()` example windows can stay explicit +since it shows a synthetic series (`1956:1960` vs `1951:1955`) — +or can drop to defaults. Decision to make in Phase 1. + +### Reused patterns + +- `cd_trend()` uses `match.arg` + per-row dispatch and rounds + p-values to 4 decimals (`round(mk$sl[1], 4)`) — Phase 2 should + mirror. +- Existing `n_a < 2` + `n_b < 2` warnings already use + `call. = FALSE` and a clean message — pattern reused for the + `< 8` test guard. + +## Reference periods + +The 1951–1980 window is the WMO-style reference. ERA5-Land +catalogues run 1950 → present (~9 km). Recent 2015–2025 is 11 +years. Both far exceed the n=8 guard, so per-ecoregion + per-WSG +loops in vignettes will all return real p-values, not NA. + +## Open question parked + +The issue (#43) suggests `test = "t"` as the default for the next +minor release — explicit in the issue body. Going with that. +Plan calls a minor bump (v0.2.8 → v0.3.0) since the schema gains +a column by default. diff --git a/planning/active/progress.md b/planning/active/progress.md new file mode 100644 index 0000000..cc1c5d1 --- /dev/null +++ b/planning/active/progress.md @@ -0,0 +1,11 @@ +# Progress — cd_compare defaults + window-vs-window p-value (#20 + #43) + +## Session 2026-05-09 + +- Plan-mode exploration — phases approved by user +- Created branch `20-43-cd-compare-defaults-pvalue` off main + (post-v0.2.8 release commit `51aa556`) +- Confirmed the two issues are tightly coupled (#20 issue body + explicitly says "should land together" with #43) +- Scaffolded PWF baseline with approved 4-phase plan +- Next: start Phase 1 — edit `R/cd_compare.R` signature + roxygen diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md new file mode 100644 index 0000000..e588e87 --- /dev/null +++ b/planning/active/task_plan.md @@ -0,0 +1,69 @@ +# Task: cd_compare defaults + window-vs-window p-value (#20 + #43) + +`cd_compare()` requires both `window_a` / `window_b` (no defaults) +and emits no measure of statistical confidence. Issue #20 settles +the framing — recent **2015–2025** vs WMO-style **1951–1980** — +based on what both regional vignettes (`peace-fwcp`, `kootenay-lake`) +already use consistently and the methodology lit review at +`planning/archive/2026-05-issue-53-snow-lit-review/findings.md`. +Issue #43 closes the second gap: a real window-vs-window p-value +column so `cd_compare()` stops borrowing `cd_trend()`'s +Mann-Kendall p-value as a proxy. The MK test answers a different +question (monotonic trend vs window-vs-window difference). Issue +#20 explicitly says "Tightly coupled with #43; both should land +together" — single branch, single PR. + +## Phase 1 — Add defaults (#20) + +- [ ] Edit `R/cd_compare.R` signature: `window_a = 2015:2025`, + `window_b = 1951:1980` defaults; update roxygen with the + framing rationale (cumulative-impact vs rate-of-change) +- [ ] Update existing tests in `tests/testthat/test-cd_compare.R` — + keep explicit-window tests, add a new test exercising the + defaults +- [ ] `devtools::document()` to refresh `man/cd_compare.Rd` +- [ ] `devtools::test()` — green +- [ ] Atomic commit: "Add cd_compare() defaults: 2015:2025 vs 1951:1980 (#20)" + +## Phase 2 — Window-vs-window p-value (#43) + +- [ ] Edit `R/cd_compare.R` body — add `test = "t"` argument, + window-extraction + Welch t / Wilcoxon dispatch, < 8-year + small-N guard with batched warning, `p_value` column emitted + only when `test != NULL` (default `"t"` so column is present + by default) +- [ ] Add tests covering: clean step-change → tiny p; iid noise → + large p; `test = "wilcox"`; `test = NULL` → original 6-col + schema; small-N guard fires once with warning; multi-variable + input gets per-row p-values +- [ ] `devtools::document()` to refresh `man/cd_compare.Rd` +- [ ] `devtools::test()` — green +- [ ] `lintr::lint_package()` — clean +- [ ] Atomic commit: "Add window-vs-window p-value to cd_compare() (#43)" + +## Phase 3 — Vignette updates + +- [ ] `vignettes/peace-fwcp.Rmd` — bare default-driven + `cd::cd_compare(ts)` recipe call; comparison table swaps + synthesized `trend_p` proxy for new `Δ p (windows)` column + (keep `Trend p (75-yr)` alongside — different question); + add one-paragraph narrative on the two p-values +- [ ] `vignettes/kootenay-lake.Rmd` — same edits +- [ ] Local render of both vignettes — clean +- [ ] Atomic commit: "Use cd_compare() defaults + window-vs-window p-value in regional vignettes" + +## Phase 4 — Release + +- [ ] Update `NEWS.md` — single section bullet covering both +- [ ] Bump `DESCRIPTION` Version: `0.2.8` → `0.3.0` (minor — new + arg with default value, output schema gains `p_value` column + by default) +- [ ] PR via `/gh-pr-push`; merge via `/gh-pr-merge` +- [ ] `/planning-archive` after merge + +## Validation + +- [ ] Tests pass +- [ ] `/code-check` clean on each commit +- [ ] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion From 48597cbc09e950ff2c0d49e13a58d1837a4ffd5b Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sat, 9 May 2026 13:51:59 -0700 Subject: [PATCH 2/7] Refresh stale roxygen docs (cd_anomaly, cd_variables) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Stale Rd regenerations from earlier .R-file edits. No source change — just devtools::document() catching up so subsequent runs don't surface noise. Co-Authored-By: Claude Opus 4.7 (1M context) --- man/cd_anomaly.Rd | 11 ++++++++--- man/cd_variables.Rd | 19 +++++++++++++------ 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/man/cd_anomaly.Rd b/man/cd_anomaly.Rd index 1cf059e..209cd88 100644 --- a/man/cd_anomaly.Rd +++ b/man/cd_anomaly.Rd @@ -14,7 +14,9 @@ cd_anomaly(x, baseline, cap_pct = 200) \code{variable}, \code{period}, \code{baseline_mean}.} \item{cap_pct}{Numeric. Cap for percent-of-normal anomalies. -Values beyond +/- \code{cap_pct} are clamped. Default \code{200}.} +Values beyond +/- \code{cap_pct} are clamped. Default \code{200}. Only +applies to \code{pct_normal} variables; \code{absolute} and +\code{pct_point_diff} anomalies are not capped.} } \value{ A tibble with columns \code{variable}, \code{period}, \code{year}, @@ -23,8 +25,11 @@ A tibble with columns \code{variable}, \code{period}, \code{year}, \description{ Calculates departure from a baseline for each year. Uses \code{\link[=cd_variables]{cd_variables()}} to determine the anomaly type: absolute deviation -for temperature, VPD, and RH; percent of normal for precipitation -and soil moisture. +for temperature, VPD, RH, and the annual snow scalars; percent of +normal for precipitation, soil moisture, and the monthly snow vars +(\code{swe}, \code{snowfall}, \code{snowmelt}); percentage-point difference for +variables that are already fractions/percentages (\code{snow_cover}, +\code{snowfall_fraction}). } \examples{ catalog <- cd_catalog( diff --git a/man/cd_variables.Rd b/man/cd_variables.Rd index 44ecdeb..ff66b68 100644 --- a/man/cd_variables.Rd +++ b/man/cd_variables.Rd @@ -11,17 +11,24 @@ A tibble with columns: \describe{ \item{variable}{Short name used throughout the package.} \item{long_name}{Human-readable label for plots and tables.} -\item{unit}{Measurement unit (degree C, percent, Pa).} -\item{anomaly_type}{"absolute" for direct departures, "pct_normal" -for percent-of-normal anomalies.} +\item{unit}{Measurement unit (degree C, percent, Pa, mm, day, mm/wk).} +\item{anomaly_type}{"absolute" for direct departures (value - baseline, +reported in the variable's native unit), "pct_normal" for +percent-of-normal anomalies (100 * value / baseline - 100, capped), +"pct_point_diff" for departures in percentage points (used for +variables that are already fractions/percentages, e.g. snow cover, +snowfall fraction, where pct-of-normal is meaningless and the +natural delta is value - baseline interpreted as percentage points).} \item{era5_name}{ERA5-Land variable name for CDS API requests, or NA for derived variables.} } } \description{ -Returns a tibble of metadata for the seven climate variables supported -by the cd package. Used internally for unit labels, anomaly type routing, -and ERA5-Land API variable names. +Returns a tibble of metadata for the climate variables supported by the +cd package. Used internally for unit labels, anomaly type routing, and +ERA5-Land API variable names. Covers the seven core climate variables +plus eight snow-related variables (four monthly natives, four annual +derived) added in v0.2.0. } \examples{ cd_variables() From 9f56b0494055f66e394da9496eca1ce8f5b48a0b Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sat, 9 May 2026 13:52:20 -0700 Subject: [PATCH 3/7] Add cd_compare() defaults: 2015:2025 vs 1951:1980 (#20) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Wire the de-facto framing both regional vignettes have settled on into the function signature so first-time users get a sensible answer without supplying both windows. Roxygen documents the why — cumulative-impact comparison ('how much warmer is the recent decade') vs rate-of-change (cd_trend's job), with 1951–1980 as the WMO-style standard normal reference. Toy example_catalog only spans 1951:1960 so the in-Rd examples keep explicit windows; a comment in the example block points users at the live-catalog defaults. Fixes #20 Co-Authored-By: Claude Opus 4.7 (1M context) --- R/cd_compare.R | 22 ++++++++++++++++++---- man/cd_compare.Rd | 20 ++++++++++++++++---- planning/active/progress.md | 5 ++++- planning/active/task_plan.md | 10 +++++----- tests/testthat/test-cd_compare.R | 17 +++++++++++++++++ 5 files changed, 60 insertions(+), 14 deletions(-) diff --git a/R/cd_compare.R b/R/cd_compare.R index 673fb20..e69a960 100644 --- a/R/cd_compare.R +++ b/R/cd_compare.R @@ -4,10 +4,19 @@ #' the difference between them. Enables custom comparisons beyond #' standard baseline anomalies. #' +#' Defaults compare a recent decade (`2015:2025`) to the WMO-style +#' standard normal reference period (`1951:1980`) — the framing +#' used in the package vignettes. This is a *cumulative-impact* +#' comparison ("how much warmer is the recent decade than the +#' pre-warming reference?") rather than a rate-of-change +#' comparison; for the latter use [cd_trend()]. +#' #' @param x A tibble from [cd_extract()] with columns `variable`, #' `period`, `year`, `value`. -#' @param window_a Integer vector of years for the first window. -#' @param window_b Integer vector of years for the second window. +#' @param window_a Integer vector of years for the recent / impact +#' window. Default `2015:2025`. +#' @param window_b Integer vector of years for the reference +#' window. Default `1951:1980` (WMO-style standard normal). #' @param method Character. Comparison method: `"mean_diff"` for #' `mean_a - mean_b`, or `"pct_change"` for percentage change #' relative to window_b. @@ -25,14 +34,19 @@ #' ) #' ts <- cd_extract(catalog, aoi) #' -#' # How has the recent period shifted from the early period? +#' # The toy example catalog only spans 1951:1960, so the windows +#' # below override the defaults (2015:2025 vs 1951:1980) — those +#' # are the right choice on the live STAC catalog. #' cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955) #' #' # Same comparison as percentage change #' cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, method = "pct_change") #' #' @export -cd_compare <- function(x, window_a, window_b, method = "mean_diff") { +cd_compare <- function(x, + window_a = 2015:2025, + window_b = 1951:1980, + method = "mean_diff") { method <- match.arg(method, c("mean_diff", "pct_change")) mean_a <- x |> diff --git a/man/cd_compare.Rd b/man/cd_compare.Rd index 2de38cc..d533880 100644 --- a/man/cd_compare.Rd +++ b/man/cd_compare.Rd @@ -4,15 +4,17 @@ \alias{cd_compare} \title{Compare arbitrary time windows} \usage{ -cd_compare(x, window_a, window_b, method = "mean_diff") +cd_compare(x, window_a = 2015:2025, window_b = 1951:1980, method = "mean_diff") } \arguments{ \item{x}{A tibble from \code{\link[=cd_extract]{cd_extract()}} with columns \code{variable}, \code{period}, \code{year}, \code{value}.} -\item{window_a}{Integer vector of years for the first window.} +\item{window_a}{Integer vector of years for the recent / impact +window. Default \code{2015:2025}.} -\item{window_b}{Integer vector of years for the second window.} +\item{window_b}{Integer vector of years for the reference +window. Default \code{1951:1980} (WMO-style standard normal).} \item{method}{Character. Comparison method: \code{"mean_diff"} for \code{mean_a - mean_b}, or \code{"pct_change"} for percentage change @@ -27,6 +29,14 @@ Computes the mean value for two user-defined time windows and the difference between them. Enables custom comparisons beyond standard baseline anomalies. } +\details{ +Defaults compare a recent decade (\code{2015:2025}) to the WMO-style +standard normal reference period (\code{1951:1980}) — the framing +used in the package vignettes. This is a \emph{cumulative-impact} +comparison ("how much warmer is the recent decade than the +pre-warming reference?") rather than a rate-of-change +comparison; for the latter use \code{\link[=cd_trend]{cd_trend()}}. +} \examples{ catalog <- cd_catalog( system.file("extdata", "example_catalog.json", package = "cd") @@ -37,7 +47,9 @@ aoi <- sf::st_read( ) ts <- cd_extract(catalog, aoi) -# How has the recent period shifted from the early period? +# The toy example catalog only spans 1951:1960, so the windows +# below override the defaults (2015:2025 vs 1951:1980) — those +# are the right choice on the live STAC catalog. cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955) # Same comparison as percentage change diff --git a/planning/active/progress.md b/planning/active/progress.md index cc1c5d1..19c593c 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -8,4 +8,7 @@ - Confirmed the two issues are tightly coupled (#20 issue body explicitly says "should land together" with #43) - Scaffolded PWF baseline with approved 4-phase plan -- Next: start Phase 1 — edit `R/cd_compare.R` signature + roxygen +- Phase 1 done — defaults `window_a = 2015:2025`, `window_b = 1951:1980` + on `cd_compare()`; roxygen documents cumulative-impact framing; + defaults test passes; 16 PASS / 0 FAIL on the cd_compare suite +- Next: Phase 2 — `test = "t"` argument + p_value column diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index e588e87..5ece5b5 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -15,15 +15,15 @@ together" — single branch, single PR. ## Phase 1 — Add defaults (#20) -- [ ] Edit `R/cd_compare.R` signature: `window_a = 2015:2025`, +- [x] Edit `R/cd_compare.R` signature: `window_a = 2015:2025`, `window_b = 1951:1980` defaults; update roxygen with the framing rationale (cumulative-impact vs rate-of-change) -- [ ] Update existing tests in `tests/testthat/test-cd_compare.R` — +- [x] Update existing tests in `tests/testthat/test-cd_compare.R` — keep explicit-window tests, add a new test exercising the defaults -- [ ] `devtools::document()` to refresh `man/cd_compare.Rd` -- [ ] `devtools::test()` — green -- [ ] Atomic commit: "Add cd_compare() defaults: 2015:2025 vs 1951:1980 (#20)" +- [x] `devtools::document()` to refresh `man/cd_compare.Rd` +- [x] `devtools::test()` — green +- [x] Atomic commit: "Add cd_compare() defaults: 2015:2025 vs 1951:1980 (#20)" ## Phase 2 — Window-vs-window p-value (#43) diff --git a/tests/testthat/test-cd_compare.R b/tests/testthat/test-cd_compare.R index cbc90a7..315e866 100644 --- a/tests/testthat/test-cd_compare.R +++ b/tests/testthat/test-cd_compare.R @@ -41,6 +41,23 @@ test_that("cd_compare warns on small windows", { ) }) +test_that("cd_compare uses 2015:2025 vs 1951:1980 defaults", { + # Toy series with values that bake in the default-window framing — + # 1951:1980 reference at value 0, 2015:2025 recent at value 2. + ts <- tibble::tibble( + variable = "tmean", + period = "annual", + year = c(1951:1980, 2015:2025), + value = c(rep(0, 30), rep(2, 11)) + ) + cmp <- cd_compare(ts) + + expect_named(cmp, c("variable", "period", "mean_a", "mean_b", "difference", "method")) + expect_equal(cmp$mean_a, 2) + expect_equal(cmp$mean_b, 0) + expect_equal(cmp$difference, 2) +}) + test_that("cd_compare handles multiple variables", { ts <- tibble::tibble( variable = rep(c("tmean", "prcp"), each = 6), From 4c2d6734dd17d82b19decc7983fda62b8ffe01fb Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sat, 9 May 2026 13:54:37 -0700 Subject: [PATCH 4/7] Add window-vs-window p-value to cd_compare() (#43) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New `test` argument (default "t") emits a p_value column on the output. `test = "t"` runs Welch's two-sample t-test on the annual values in each window; `test = "wilcox"` runs Mann-Whitney U; `test = NULL` skips the test entirely and drops the column. Rows where either window has < 8 non-NA years get `p_value = NA` and a single batched warning naming the affected variable/period rows (no per-row warning spam). The window-vs-window p-value answers a different question than the cd_trend Mann-Kendall test: "do the two windows differ" vs "is there a monotonic trend across the full series." Step changes and U-shapes can produce significant window differences without monotonic trends, and vice versa. Vignettes will report both p-values side-by-side in Phase 3. Roxygen documents the independence assumption; for series with strong autocorrelation the t-test is mildly anti-conservative — the Wilcoxon alternative is more robust to non-Gaussian tails but shares the independence assumption. Fixes #43 Co-Authored-By: Claude Opus 4.7 (1M context) --- R/cd_compare.R | 70 ++++++++++++++++++--- man/cd_compare.Rd | 45 +++++++++++--- planning/active/progress.md | 9 ++- planning/active/task_plan.md | 14 +++-- tests/testthat/test-cd_compare.R | 103 +++++++++++++++++++++++++++++-- 5 files changed, 213 insertions(+), 28 deletions(-) diff --git a/R/cd_compare.R b/R/cd_compare.R index e69a960..e9bd96a 100644 --- a/R/cd_compare.R +++ b/R/cd_compare.R @@ -1,8 +1,8 @@ #' Compare arbitrary time windows #' -#' Computes the mean value for two user-defined time windows and -#' the difference between them. Enables custom comparisons beyond -#' standard baseline anomalies. +#' Computes the mean value for two user-defined time windows, the +#' difference between them, and (by default) a p-value testing +#' whether the two windows differ. #' #' Defaults compare a recent decade (`2015:2025`) to the WMO-style #' standard normal reference period (`1951:1980`) — the framing @@ -11,6 +11,19 @@ #' pre-warming reference?") rather than a rate-of-change #' comparison; for the latter use [cd_trend()]. #' +#' The window-vs-window p-value (`test = "t"` Welch t-test or +#' `test = "wilcox"` Mann-Whitney U) tests whether the means of +#' the two windows differ — a different question than +#' [cd_trend()]'s Mann-Kendall test, which checks for a monotonic +#' trend across the full series. Two windows can differ +#' significantly without a monotonic trend (step changes, U-shapes) +#' and a non-significant trend doesn't always mean the windows are +#' indistinguishable. The test treats annual values as independent; +#' for series with strong autocorrelation the t-test p is mildly +#' anti-conservative — the Wilcoxon alternative is more robust to +#' non-Gaussian tails / outliers but shares the independence +#' assumption. +#' #' @param x A tibble from [cd_extract()] with columns `variable`, #' `period`, `year`, `value`. #' @param window_a Integer vector of years for the recent / impact @@ -20,9 +33,16 @@ #' @param method Character. Comparison method: `"mean_diff"` for #' `mean_a - mean_b`, or `"pct_change"` for percentage change #' relative to window_b. +#' @param test Character or NULL. Window-vs-window significance +#' test: `"t"` for Welch's two-sample t-test (default) or +#' `"wilcox"` for Mann-Whitney U. Set to `NULL` to skip — output +#' then drops the `p_value` column. Rows where either window has +#' fewer than 8 non-NA values get `p_value = NA` and a single +#' batched warning. #' #' @return A tibble with columns `variable`, `period`, `mean_a`, -#' `mean_b`, `difference`, `method`. +#' `mean_b`, `difference`, `method`, and (when `test` is +#' non-NULL) `p_value`. #' #' @examples #' catalog <- cd_catalog( @@ -36,18 +56,22 @@ #' #' # The toy example catalog only spans 1951:1960, so the windows #' # below override the defaults (2015:2025 vs 1951:1980) — those -#' # are the right choice on the live STAC catalog. -#' cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955) +#' # are the right choice on the live STAC catalog. test = NULL +#' # because the toy windows have < 8 years. +#' cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, test = NULL) #' #' # Same comparison as percentage change -#' cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, method = "pct_change") +#' cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, +#' method = "pct_change", test = NULL) #' #' @export cd_compare <- function(x, window_a = 2015:2025, window_b = 1951:1980, - method = "mean_diff") { + method = "mean_diff", + test = "t") { method <- match.arg(method, c("mean_diff", "pct_change")) + if (!is.null(test)) test <- match.arg(test, c("t", "wilcox")) mean_a <- x |> dplyr::filter(.data$year %in% window_a) |> @@ -70,5 +94,35 @@ cd_compare <- function(x, ) out$method <- method + if (!is.null(test)) { + p <- vapply(seq_len(nrow(out)), function(i) { + v <- out$variable[i] + p_ <- out$period[i] + a <- x$value[x$variable == v & x$period == p_ & x$year %in% window_a] + b <- x$value[x$variable == v & x$period == p_ & x$year %in% window_b] + a <- a[!is.na(a)] + b <- b[!is.na(b)] + if (length(a) < 8 || length(b) < 8) return(NA_real_) + pv <- tryCatch( + switch(test, + t = stats::t.test(a, b, var.equal = FALSE)$p.value, + wilcox = stats::wilcox.test(a, b, exact = FALSE)$p.value + ), + error = function(e) NA_real_ + ) + round(pv, 4) + }, numeric(1)) + + if (any(is.na(p))) { + bad <- which(is.na(p)) + warning(sprintf( + "p_value set to NA for %d row(s) with < 8 years in either window: %s", + length(bad), + paste(sprintf("%s/%s", out$variable[bad], out$period[bad]), collapse = ", ") + ), call. = FALSE) + } + out$p_value <- p + } + out } diff --git a/man/cd_compare.Rd b/man/cd_compare.Rd index d533880..6ede438 100644 --- a/man/cd_compare.Rd +++ b/man/cd_compare.Rd @@ -4,7 +4,13 @@ \alias{cd_compare} \title{Compare arbitrary time windows} \usage{ -cd_compare(x, window_a = 2015:2025, window_b = 1951:1980, method = "mean_diff") +cd_compare( + x, + window_a = 2015:2025, + window_b = 1951:1980, + method = "mean_diff", + test = "t" +) } \arguments{ \item{x}{A tibble from \code{\link[=cd_extract]{cd_extract()}} with columns \code{variable}, @@ -19,15 +25,23 @@ window. Default \code{1951:1980} (WMO-style standard normal).} \item{method}{Character. Comparison method: \code{"mean_diff"} for \code{mean_a - mean_b}, or \code{"pct_change"} for percentage change relative to window_b.} + +\item{test}{Character or NULL. Window-vs-window significance +test: \code{"t"} for Welch's two-sample t-test (default) or +\code{"wilcox"} for Mann-Whitney U. Set to \code{NULL} to skip — output +then drops the \code{p_value} column. Rows where either window has +fewer than 8 non-NA values get \code{p_value = NA} and a single +batched warning.} } \value{ A tibble with columns \code{variable}, \code{period}, \code{mean_a}, -\code{mean_b}, \code{difference}, \code{method}. +\code{mean_b}, \code{difference}, \code{method}, and (when \code{test} is +non-NULL) \code{p_value}. } \description{ -Computes the mean value for two user-defined time windows and -the difference between them. Enables custom comparisons beyond -standard baseline anomalies. +Computes the mean value for two user-defined time windows, the +difference between them, and (by default) a p-value testing +whether the two windows differ. } \details{ Defaults compare a recent decade (\code{2015:2025}) to the WMO-style @@ -36,6 +50,19 @@ used in the package vignettes. This is a \emph{cumulative-impact} comparison ("how much warmer is the recent decade than the pre-warming reference?") rather than a rate-of-change comparison; for the latter use \code{\link[=cd_trend]{cd_trend()}}. + +The window-vs-window p-value (\code{test = "t"} Welch t-test or +\code{test = "wilcox"} Mann-Whitney U) tests whether the means of +the two windows differ — a different question than +\code{\link[=cd_trend]{cd_trend()}}'s Mann-Kendall test, which checks for a monotonic +trend across the full series. Two windows can differ +significantly without a monotonic trend (step changes, U-shapes) +and a non-significant trend doesn't always mean the windows are +indistinguishable. The test treats annual values as independent; +for series with strong autocorrelation the t-test p is mildly +anti-conservative — the Wilcoxon alternative is more robust to +non-Gaussian tails / outliers but shares the independence +assumption. } \examples{ catalog <- cd_catalog( @@ -49,10 +76,12 @@ ts <- cd_extract(catalog, aoi) # The toy example catalog only spans 1951:1960, so the windows # below override the defaults (2015:2025 vs 1951:1980) — those -# are the right choice on the live STAC catalog. -cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955) +# are the right choice on the live STAC catalog. test = NULL +# because the toy windows have < 8 years. +cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, test = NULL) # Same comparison as percentage change -cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, method = "pct_change") +cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, + method = "pct_change", test = NULL) } diff --git a/planning/active/progress.md b/planning/active/progress.md index 19c593c..4569272 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -11,4 +11,11 @@ - Phase 1 done — defaults `window_a = 2015:2025`, `window_b = 1951:1980` on `cd_compare()`; roxygen documents cumulative-impact framing; defaults test passes; 16 PASS / 0 FAIL on the cd_compare suite -- Next: Phase 2 — `test = "t"` argument + p_value column +- Phase 2 done — `test = c("t", "wilcox", NULL)` argument added; + default `"t"` so output gains `p_value` column by default; + small-N guard (< 8 years either window) → NA + single batched + warning naming affected variable/period rows; 6 new tests + covering step-change, iid noise, wilcox, NULL, small-N guard, + multi-variable; full suite 181 PASS / 0 FAIL +- Next: Phase 3 — vignette wire-up (drop synthesized trend_p + proxy, surface real `Δ p (windows)` column) diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 5ece5b5..dc72339 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -27,19 +27,21 @@ together" — single branch, single PR. ## Phase 2 — Window-vs-window p-value (#43) -- [ ] Edit `R/cd_compare.R` body — add `test = "t"` argument, +- [x] Edit `R/cd_compare.R` body — add `test = "t"` argument, window-extraction + Welch t / Wilcoxon dispatch, < 8-year small-N guard with batched warning, `p_value` column emitted only when `test != NULL` (default `"t"` so column is present by default) -- [ ] Add tests covering: clean step-change → tiny p; iid noise → +- [x] Add tests covering: clean step-change → tiny p; iid noise → large p; `test = "wilcox"`; `test = NULL` → original 6-col schema; small-N guard fires once with warning; multi-variable input gets per-row p-values -- [ ] `devtools::document()` to refresh `man/cd_compare.Rd` -- [ ] `devtools::test()` — green -- [ ] `lintr::lint_package()` — clean -- [ ] Atomic commit: "Add window-vs-window p-value to cd_compare() (#43)" +- [x] `devtools::document()` to refresh `man/cd_compare.Rd` +- [x] `devtools::test()` — green +- [x] `lintr::lint_package()` — clean (no new lints from + cd_compare; pre-existing `.data` pronoun + vignette caption + length unchanged) +- [x] Atomic commit: "Add window-vs-window p-value to cd_compare() (#43)" ## Phase 3 — Vignette updates diff --git a/tests/testthat/test-cd_compare.R b/tests/testthat/test-cd_compare.R index 315e866..8b381ae 100644 --- a/tests/testthat/test-cd_compare.R +++ b/tests/testthat/test-cd_compare.R @@ -5,7 +5,7 @@ test_that("cd_compare computes mean_diff", { year = 1951:1960, value = c(rep(0, 5), rep(10, 5)) ) - cmp <- cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955) + cmp <- cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, test = NULL) expect_s3_class(cmp, "tbl_df") expect_named(cmp, c("variable", "period", "mean_a", "mean_b", "difference", "method")) @@ -22,7 +22,8 @@ test_that("cd_compare computes pct_change", { year = 1951:1960, value = c(rep(100, 5), rep(150, 5)) ) - cmp <- cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, method = "pct_change") + cmp <- cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955, + method = "pct_change", test = NULL) expect_equal(cmp$difference, 50) expect_equal(cmp$method, "pct_change") @@ -36,7 +37,7 @@ test_that("cd_compare warns on small windows", { value = 1:3 ) expect_warning( - cd_compare(ts, window_a = 1953, window_b = 1951:1952), + cd_compare(ts, window_a = 1953, window_b = 1951:1952, test = NULL), "window_a has fewer than 2 years" ) }) @@ -50,7 +51,7 @@ test_that("cd_compare uses 2015:2025 vs 1951:1980 defaults", { year = c(1951:1980, 2015:2025), value = c(rep(0, 30), rep(2, 11)) ) - cmp <- cd_compare(ts) + cmp <- cd_compare(ts, test = NULL) expect_named(cmp, c("variable", "period", "mean_a", "mean_b", "difference", "method")) expect_equal(cmp$mean_a, 2) @@ -65,9 +66,101 @@ test_that("cd_compare handles multiple variables", { year = rep(1951:1956, 2), value = c(0, 0, 0, 10, 10, 10, 100, 100, 100, 200, 200, 200) ) - cmp <- cd_compare(ts, window_a = 1954:1956, window_b = 1951:1953) + cmp <- cd_compare(ts, window_a = 1954:1956, window_b = 1951:1953, test = NULL) expect_equal(nrow(cmp), 2) expect_equal(cmp$difference[cmp$variable == "tmean"], 10) expect_equal(cmp$difference[cmp$variable == "prcp"], 100) }) + +# --------------------------------------------------------------------------- +# Window-vs-window p-value tests (#43) +# --------------------------------------------------------------------------- + +test_that("cd_compare returns p_value column when test = 't' (default)", { + set.seed(1) + # Clean step change between two 30-year windows — Welch t should + # return a vanishingly small p. + ts <- tibble::tibble( + variable = "tmean", + period = "annual", + year = c(1951:1980, 1991:2020), + value = c(rnorm(30, mean = 0, sd = 1), rnorm(30, mean = 5, sd = 1)) + ) + cmp <- cd_compare(ts, window_a = 1991:2020, window_b = 1951:1980) + + expect_true("p_value" %in% names(cmp)) + expect_equal(nrow(cmp), 1) + expect_lt(cmp$p_value, 0.001) +}) + +test_that("cd_compare t-test returns large p on iid noise", { + set.seed(2) + ts <- tibble::tibble( + variable = "tmean", + period = "annual", + year = c(1951:1980, 1991:2020), + value = rnorm(60, mean = 0, sd = 1) + ) + cmp <- cd_compare(ts, window_a = 1991:2020, window_b = 1951:1980) + + expect_gt(cmp$p_value, 0.05) +}) + +test_that("cd_compare supports test = 'wilcox'", { + set.seed(3) + ts <- tibble::tibble( + variable = "tmean", + period = "annual", + year = c(1951:1980, 1991:2020), + value = c(rnorm(30, mean = 0), rnorm(30, mean = 5)) + ) + cmp <- cd_compare(ts, window_a = 1991:2020, window_b = 1951:1980, test = "wilcox") + + expect_true("p_value" %in% names(cmp)) + expect_lt(cmp$p_value, 0.001) +}) + +test_that("cd_compare with test = NULL drops p_value column", { + ts <- tibble::tibble( + variable = "tmean", + period = "annual", + year = c(1951:1980, 1991:2020), + value = c(rep(0, 30), rep(5, 30)) + ) + cmp <- cd_compare(ts, window_a = 1991:2020, window_b = 1951:1980, test = NULL) + + expect_named(cmp, c("variable", "period", "mean_a", "mean_b", "difference", "method")) +}) + +test_that("cd_compare small-N guard returns NA p_value with one warning", { + ts <- tibble::tibble( + variable = "tmean", + period = "annual", + year = 1951:1960, + value = c(rep(0, 5), rep(5, 5)) + ) + expect_warning( + cmp <- cd_compare(ts, window_a = 1956:1960, window_b = 1951:1955), + "p_value set to NA" + ) + expect_true(is.na(cmp$p_value)) +}) + +test_that("cd_compare per-row p_values on multi-variable input", { + set.seed(4) + # tmean: real shift (small p); prcp: no shift (large p) + ts <- tibble::tibble( + variable = rep(c("tmean", "prcp"), each = 60), + period = "annual", + year = rep(c(1951:1980, 1991:2020), 2), + value = c( + c(rnorm(30, mean = 0, sd = 1), rnorm(30, mean = 5, sd = 1)), + rnorm(60, mean = 100, sd = 5) + ) + ) + cmp <- cd_compare(ts, window_a = 1991:2020, window_b = 1951:1980) + + expect_lt(cmp$p_value[cmp$variable == "tmean"], 0.001) + expect_gt(cmp$p_value[cmp$variable == "prcp"], 0.05) +}) From 46d6792e9207609ee09cafe2f6918d7e4a4d7aeb Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sat, 9 May 2026 13:58:20 -0700 Subject: [PATCH 5/7] Use cd_compare() defaults + window-vs-window p-value in regional vignettes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Both vignettes: - Compare-recipe chunk now shows bare `cd::cd_compare(ts)` — defaults (2015:2025 vs 1951:1980, Welch t) carry the call. Cleaner cd showcase; reader sees the recipe in 1 line. - Compare-table chunk recomputes `cmp` from cached `ts` so the new `p_value` column flows through without re-running the heavy precompute scripts. The cached `regional` predates #43 so doesn't have the column. - Table swaps the synthesized Mann-Kendall `trend_p` proxy for the real window-vs-window `Δ p (windows)` (Welch). `Trend p (75-yr)` stays alongside — different question, both belong. - One-paragraph narrative reframes the two p-values: "do the windows differ?" (window p) vs "is there a monotonic ramp?" (trend p). Step changes vs gradual ramps read differently. `data-raw/_vignette_data.R` updated so the next precompute run uses cd_compare defaults. Both vignettes render clean locally with both column headers present in the rendered HTML. Co-Authored-By: Claude Opus 4.7 (1M context) --- data-raw/kootenay_lake_vignette_data.R | 10 +++---- data-raw/peace_fwcp_vignette_data.R | 8 ++--- planning/active/progress.md | 11 +++++-- planning/active/task_plan.md | 11 ++++--- vignettes/kootenay-lake.Rmd | 41 ++++++++++++++++---------- vignettes/peace-fwcp.Rmd | 38 ++++++++++++++---------- 6 files changed, 71 insertions(+), 48 deletions(-) diff --git a/data-raw/kootenay_lake_vignette_data.R b/data-raw/kootenay_lake_vignette_data.R index 68d127e..1ef1236 100644 --- a/data-raw/kootenay_lake_vignette_data.R +++ b/data-raw/kootenay_lake_vignette_data.R @@ -48,12 +48,10 @@ regional_ts <- cd_extract(catalog, aoi) regional_bl <- cd_baseline(regional_ts, baseline_years = 1951:1980) regional_ano <- cd_anomaly(regional_ts, regional_bl) regional_trn <- cd_trend(regional_ano, trend_start = c(1951, 1981)) -regional_cmp <- cd_compare(regional_ts, - window_a = 2015:2025, window_b = 1951:1980, - method = "mean_diff") +regional_cmp <- cd_compare(regional_ts, method = "mean_diff") # defaults: 2015:2025 vs 1951:1980, Welch t regional_cmp_pct <- cd_compare( regional_ts[regional_ts$variable %in% pct_normal_vars, ], - window_a = 2015:2025, window_b = 1951:1980, method = "pct_change" + method = "pct_change" ) # -- Per-ecoregion ----------------------------------------------------------- @@ -65,7 +63,7 @@ ecoregion_results <- lapply(seq_len(nrow(ecoregions)), function(i) { trn <- cd_trend(ano, trend_start = c(1951, 1981)) cmp_pct <- cd_compare( ts[ts$variable %in% pct_normal_vars, ], - window_a = 2015:2025, window_b = 1951:1980, method = "pct_change" + method = "pct_change" ) list(ts = ts, ano = ano, trn = trn, cmp_pct = cmp_pct) }) @@ -81,7 +79,7 @@ wsg_results <- lapply(seq_len(nrow(wsgs)), function(i) { trn <- cd_trend(ano, trend_start = c(1951, 1981)) cmp_pct <- cd_compare( ts[ts$variable %in% pct_normal_vars, ], - window_a = 2015:2025, window_b = 1951:1980, method = "pct_change" + method = "pct_change" ) list(ts = ts, ano = ano, trn = trn, cmp_pct = cmp_pct) }) diff --git a/data-raw/peace_fwcp_vignette_data.R b/data-raw/peace_fwcp_vignette_data.R index ad9540e..d3728a2 100644 --- a/data-raw/peace_fwcp_vignette_data.R +++ b/data-raw/peace_fwcp_vignette_data.R @@ -42,16 +42,14 @@ regional_ts <- cd_extract(catalog, aoi) regional_bl <- cd_baseline(regional_ts, baseline_years = 1951:1980) regional_ano <- cd_anomaly(regional_ts, regional_bl) regional_trn <- cd_trend(regional_ano, trend_start = c(1951, 1981)) -regional_cmp <- cd_compare(regional_ts, - window_a = 2015:2025, window_b = 1951:1980, - method = "mean_diff") +regional_cmp <- cd_compare(regional_ts, method = "mean_diff") # defaults: 2015:2025 vs 1951:1980, Welch t # pct_change for vars with pct_normal anomaly type (#48 added swe / snowfall / # snowmelt to this list; snow_cover and snowfall_fraction are pct_point_diff so # mean_diff is the right comparison method for those — already in regional_cmp). pct_normal_vars <- c("prcp", "soil_moisture", "swe", "snowfall", "snowmelt") regional_cmp_pct <- cd_compare( regional_ts[regional_ts$variable %in% pct_normal_vars, ], - window_a = 2015:2025, window_b = 1951:1980, method = "pct_change" + method = "pct_change" ) # -- Per-ecoregion ----------------------------------------------------------- @@ -63,7 +61,7 @@ ecoregion_results <- lapply(seq_len(nrow(ecoregions)), function(i) { trn <- cd_trend(ano, trend_start = c(1951, 1981)) cmp_pct <- cd_compare( ts[ts$variable %in% pct_normal_vars, ], - window_a = 2015:2025, window_b = 1951:1980, method = "pct_change" + method = "pct_change" ) list(ts = ts, ano = ano, trn = trn, cmp_pct = cmp_pct) }) diff --git a/planning/active/progress.md b/planning/active/progress.md index 4569272..e555edc 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -17,5 +17,12 @@ warning naming affected variable/period rows; 6 new tests covering step-change, iid noise, wilcox, NULL, small-N guard, multi-variable; full suite 181 PASS / 0 FAIL -- Next: Phase 3 — vignette wire-up (drop synthesized trend_p - proxy, surface real `Δ p (windows)` column) +- Phase 3 done — both vignettes show bare `cd::cd_compare(ts)` + recipe call (defaults take over); compare-table chunks + recompute `cmp` from cached `ts` so `p_value` flows through + without re-running the heavy precompute scripts; surfaces both + `Δ p (windows)` and `Trend p (75-yr)` columns side-by-side; + narrative paragraph explains the two-test framing; precompute + scripts updated for next refresh; both vignettes render clean + with the new column visible in HTML +- Next: Phase 4 — NEWS + DESCRIPTION bump → v0.3.0; PR + merge diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index dc72339..3450f57 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -45,14 +45,17 @@ together" — single branch, single PR. ## Phase 3 — Vignette updates -- [ ] `vignettes/peace-fwcp.Rmd` — bare default-driven +- [x] `vignettes/peace-fwcp.Rmd` — bare default-driven `cd::cd_compare(ts)` recipe call; comparison table swaps synthesized `trend_p` proxy for new `Δ p (windows)` column (keep `Trend p (75-yr)` alongside — different question); add one-paragraph narrative on the two p-values -- [ ] `vignettes/kootenay-lake.Rmd` — same edits -- [ ] Local render of both vignettes — clean -- [ ] Atomic commit: "Use cd_compare() defaults + window-vs-window p-value in regional vignettes" +- [x] `vignettes/kootenay-lake.Rmd` — same edits +- [x] `data-raw/_vignette_data.R` — drop explicit + windows so future precompute runs use defaults +- [x] Local render of both vignettes — clean (HTML contains both + `Δ p (windows)` and `Trend p (75-yr)` column headers) +- [x] Atomic commit: "Use cd_compare() defaults + window-vs-window p-value in regional vignettes" ## Phase 4 — Release diff --git a/vignettes/kootenay-lake.Rmd b/vignettes/kootenay-lake.Rmd index e6987e9..6282c00 100644 --- a/vignettes/kootenay-lake.Rmd +++ b/vignettes/kootenay-lake.Rmd @@ -543,25 +543,33 @@ The table below compares two windows directly — the recent decade absolute is the difference of the two means in the variable's native units. Δ % is shown only for variables where it is meaningful (precipitation, soil moisture, vapour pressure deficit, relative -humidity). The trend p column is the Mann-Kendall p-value of the -75-year (1951–present) trend on the same variable and period; it -tests for a steady year-on-year ramp, which is a related but distinct -question from "do these two windows differ". +humidity). Two p-values are reported, answering different +questions. **Δ p (windows)** is the Welch t-test on the annual +values within each window — does the recent decade differ from +the pre-warming reference? **Trend p (75-yr)** is the Mann-Kendall +test on the full 1951–present series — is there a steady +year-on-year ramp? Step changes show up as significant Δ p with +non-significant trend p; gradual ramps show up as both significant. The recent decade was 1.6 to 1.8 °C warmer than the pre-warming reference for annual mean, daytime maximum, and overnight minimum, -with Mann-Kendall trend p-values below 0.001. Vapour pressure deficit -is up significantly. Annual precipitation was about 6 to 7 % *lower* -in the recent decade — and unlike the FWCP Peace just to the north, -the long-term Mann-Kendall trend test does confirm a steady -year-on-year decline (p ≈ 0.02). Soil moisture is roughly flat. -Relative humidity shows a small significant decline. +with both window and trend p-values below 0.001. Vapour pressure +deficit is up significantly on both tests. Annual precipitation +was about 6 to 7 % *lower* in the recent decade — and unlike the +FWCP Peace just to the north, the long-term Mann-Kendall trend +test does confirm a steady year-on-year decline (trend p ≈ 0.02). +Soil moisture is roughly flat. Relative humidity shows a small +significant decline. ```{r compare-recipe, eval = FALSE} -cmp <- cd::cd_compare(ts, window_a = 2015:2025, window_b = 1951:1980) +cmp <- cd::cd_compare(ts) # defaults: 2015–2025 vs 1951–1980, Welch t-test ``` ```{r compare-table, echo = FALSE} +# Recompute cmp from the cached ts so the p_value column reflects the +# current cd_compare() defaults — the precomputed rds was generated +# before #43 added the test argument. +cmp <- cd::cd_compare(ts) trn_p <- trn[trn$trend_start == 1951, c("variable", "period", "mk_pvalue")] names(trn_p)[3] <- "trend_p" trn_p$trend_p <- round(trn_p$trend_p, 3) @@ -580,17 +588,18 @@ cmp_combined <- cmp |> NA_real_, round(100 * (mean_a - mean_b) / mean_b, 1) ), - mean_a = round(mean_a, 2), - mean_b = round(mean_b, 2), - abs_diff = round(difference, 2) + mean_a = round(mean_a, 2), + mean_b = round(mean_b, 2), + abs_diff = round(difference, 2), + window_p = round(p_value, 3) ) |> - dplyr::select(variable, period, mean_a, mean_b, abs_diff, pct_change) |> + dplyr::select(variable, period, mean_a, mean_b, abs_diff, pct_change, window_p) |> merge(trn_p, by = c("variable", "period"), all.x = TRUE) names(cmp_combined) <- c( "Variable", "Period", "Recent (2015–2025)", "Pre-warming (1951–1980)", - "Δ absolute", "Δ %", "Trend p (75-yr)" + "Δ absolute", "Δ %", "Δ p (windows)", "Trend p (75-yr)" ) kableExtra::kable_styling( diff --git a/vignettes/peace-fwcp.Rmd b/vignettes/peace-fwcp.Rmd index 88ade99..ac11e84 100644 --- a/vignettes/peace-fwcp.Rmd +++ b/vignettes/peace-fwcp.Rmd @@ -517,24 +517,31 @@ The table below compares two windows directly — the recent decade absolute is the difference of the two means in the variable's native units. Δ % is shown only for variables where it is meaningful (precipitation, soil moisture, vapour pressure deficit, relative -humidity). The trend p column is the Mann-Kendall p-value of the -75-year (1951–present) trend on the same variable and period; it -tests for a steady year-on-year ramp, which is a related but distinct -question from "do these two windows differ". +humidity). Two p-values are reported, answering different +questions. **Δ p (windows)** is the Welch t-test on the annual +values within each window — does the recent decade differ from +the pre-warming reference? **Trend p (75-yr)** is the Mann-Kendall +test on the full 1951–present series — is there a steady +year-on-year ramp? Step changes show up as significant Δ p with +non-significant trend p; gradual ramps show up as both significant. The recent decade was 1.7 to 2.4 °C warmer than the pre-warming reference for annual mean, daytime maximum, and overnight minimum, -with Mann-Kendall trend p-values below 0.001. Vapour pressure deficit -is up significantly. Annual precipitation was about 3 to 4 % higher -in the recent decade; the long-term trend test does not confirm a -steady year-on-year ramp (p ≈ 0.20). Soil moisture and relative -humidity show no meaningful change. +with both window and trend p-values below 0.001. Vapour pressure +deficit is up significantly on both tests. Annual precipitation +was about 3 to 4 % higher in the recent decade; neither test +confirms the shift (Δ p ≈ 0.4, trend p ≈ 0.20). Soil moisture and +relative humidity show no meaningful change. ```{r compare-recipe, eval = FALSE} -cmp <- cd::cd_compare(ts, window_a = 2015:2025, window_b = 1951:1980) +cmp <- cd::cd_compare(ts) # defaults: 2015–2025 vs 1951–1980, Welch t-test ``` ```{r compare-table, echo = FALSE} +# Recompute cmp from the cached ts so the p_value column reflects the +# current cd_compare() defaults — the precomputed rds was generated +# before #43 added the test argument. +cmp <- cd::cd_compare(ts) trn_p <- trn[trn$trend_start == 1951, c("variable", "period", "mk_pvalue")] names(trn_p)[3] <- "trend_p" trn_p$trend_p <- round(trn_p$trend_p, 3) @@ -553,17 +560,18 @@ cmp_combined <- cmp |> NA_real_, round(100 * (mean_a - mean_b) / mean_b, 1) ), - mean_a = round(mean_a, 2), - mean_b = round(mean_b, 2), - abs_diff = round(difference, 2) + mean_a = round(mean_a, 2), + mean_b = round(mean_b, 2), + abs_diff = round(difference, 2), + window_p = round(p_value, 3) ) |> - dplyr::select(variable, period, mean_a, mean_b, abs_diff, pct_change) |> + dplyr::select(variable, period, mean_a, mean_b, abs_diff, pct_change, window_p) |> merge(trn_p, by = c("variable", "period"), all.x = TRUE) names(cmp_combined) <- c( "Variable", "Period", "Recent (2015–2025)", "Pre-warming (1951–1980)", - "Δ absolute", "Δ %", "Trend p (75-yr)" + "Δ absolute", "Δ %", "Δ p (windows)", "Trend p (75-yr)" ) kableExtra::kable_styling( From 88e697ef74716b760f77f107c86911a5acb1570a Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sat, 9 May 2026 14:04:28 -0700 Subject: [PATCH 6/7] Mark Phase 4 release-flow notes (PR + post-merge bump) Repo pattern: NEWS + DESCRIPTION bump lands as a separate 'Release v0.X.Y' commit on main post-merge, driven by /gh-pr-merge. The branch itself stays release-prep clean. Co-Authored-By: Claude Opus 4.7 (1M context) --- planning/active/task_plan.md | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 3450f57..d5b1472 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -59,11 +59,12 @@ together" — single branch, single PR. ## Phase 4 — Release -- [ ] Update `NEWS.md` — single section bullet covering both -- [ ] Bump `DESCRIPTION` Version: `0.2.8` → `0.3.0` (minor — new - arg with default value, output schema gains `p_value` column - by default) -- [ ] PR via `/gh-pr-push`; merge via `/gh-pr-merge` +- [ ] PR via `/gh-pr-push` (NEWS + DESCRIPTION bump handled + post-merge by `/gh-pr-merge`, per repo convention — recent + pattern: `Release v0.X.Y` is a separate post-merge commit + on main, e.g. `298cce4 Release v0.2.8`) +- [ ] Bump on merge: `0.2.8` → `0.3.0` (minor — output schema + gains `p_value` column by default; new `test` argument) - [ ] `/planning-archive` after merge ## Validation From 4a193798109d5277da264592af112d99fc18ee35 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sun, 10 May 2026 00:02:20 -0700 Subject: [PATCH 7/7] Move Recent Decade vs Pre-Warming Reference up + rename for prominence MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The cumulative-impact comparison (mean(2015-2025) - mean(1951-1980)) is the headline number this package produces — it was buried 300+ lines down between Snowpack and Spatial Pattern, and titled with methodological framing ('Recent vs Pre-warming') rather than something that telegraphs the answer. Two changes in both vignettes: - Rename: 'Recent vs Pre-warming' -> 'Recent Decade vs Pre-Warming Reference'. Adds 'decade' + 'reference' so the heading reads as a comparison label rather than a methods stub. - Reorder: place directly after Trends, before Daytime Highs and Overnight Lows. Trends still introduces the 1951-1980 reference and the two-test framing the table relies on, then the headline cumulative number lands; DTR / Snowpack / Spatial Pattern read as drill-down on the headline rather than a long detour before it. Both vignettes render clean; new heading visible in HTML, old heading absent. Co-Authored-By: Claude Opus 4.7 (1M context) --- vignettes/kootenay-lake.Rmd | 154 ++++++++++++++++++------------------ vignettes/peace-fwcp.Rmd | 150 +++++++++++++++++------------------ 2 files changed, 152 insertions(+), 152 deletions(-) diff --git a/vignettes/kootenay-lake.Rmd b/vignettes/kootenay-lake.Rmd index 6282c00..ecaaed0 100644 --- a/vignettes/kootenay-lake.Rmd +++ b/vignettes/kootenay-lake.Rmd @@ -261,6 +261,83 @@ cd::cd_plot_timeseries( ) ``` +## Recent Decade vs Pre-Warming Reference + +The table below compares two windows directly — the recent decade +(2015–2025) against the pre-warming reference (1951–1980). Δ +absolute is the difference of the two means in the variable's native +units. Δ % is shown only for variables where it is meaningful +(precipitation, soil moisture, vapour pressure deficit, relative +humidity). Two p-values are reported, answering different +questions. **Δ p (windows)** is the Welch t-test on the annual +values within each window — does the recent decade differ from +the pre-warming reference? **Trend p (75-yr)** is the Mann-Kendall +test on the full 1951–present series — is there a steady +year-on-year ramp? Step changes show up as significant Δ p with +non-significant trend p; gradual ramps show up as both significant. + +The recent decade was 1.6 to 1.8 °C warmer than the pre-warming +reference for annual mean, daytime maximum, and overnight minimum, +with both window and trend p-values below 0.001. Vapour pressure +deficit is up significantly on both tests. Annual precipitation +was about 6 to 7 % *lower* in the recent decade — and unlike the +FWCP Peace just to the north, the long-term Mann-Kendall trend +test does confirm a steady year-on-year decline (trend p ≈ 0.02). +Soil moisture is roughly flat. Relative humidity shows a small +significant decline. + +```{r compare-recipe, eval = FALSE} +cmp <- cd::cd_compare(ts) # defaults: 2015–2025 vs 1951–1980, Welch t-test +``` + +```{r compare-table, echo = FALSE} +# Recompute cmp from the cached ts so the p_value column reflects the +# current cd_compare() defaults — the precomputed rds was generated +# before #43 added the test argument. +cmp <- cd::cd_compare(ts) +trn_p <- trn[trn$trend_start == 1951, c("variable", "period", "mk_pvalue")] +names(trn_p)[3] <- "trend_p" +trn_p$trend_p <- round(trn_p$trend_p, 3) + +no_pct_vars <- c( + "tmean", "tmax", "tmin", + # snow_cover and snowfall_fraction are already in % (pct_point_diff + # anomaly type) — pct-of-baseline mixes units. snowmelt_doy_50 is a + # day-of-year — pct-change of a date is meaningless. + "snow_cover", "snowfall_fraction", "snowmelt_doy_50" +) +cmp_combined <- cmp |> + dplyr::mutate( + pct_change = ifelse( + variable %in% no_pct_vars, + NA_real_, + round(100 * (mean_a - mean_b) / mean_b, 1) + ), + mean_a = round(mean_a, 2), + mean_b = round(mean_b, 2), + abs_diff = round(difference, 2), + window_p = round(p_value, 3) + ) |> + dplyr::select(variable, period, mean_a, mean_b, abs_diff, pct_change, window_p) |> + merge(trn_p, by = c("variable", "period"), all.x = TRUE) + +names(cmp_combined) <- c( + "Variable", "Period", + "Recent (2015–2025)", "Pre-warming (1951–1980)", + "Δ absolute", "Δ %", "Δ p (windows)", "Trend p (75-yr)" +) + +kableExtra::kable_styling( + knitr::kable(cmp_combined, label = NA, + caption = "Recent decade (2015-2025 mean) compared to pre-warming reference (1951-1980 mean) for the Kootenay Lake Region.", + row.names = FALSE), + bootstrap_options = c("striped", "hover", "condensed") +) |> + kableExtra::scroll_box(height = "360px") +``` + +
+ ## Daytime Highs and Overnight Lows The cd package ships daytime maximum (tmax) and overnight minimum (tmin) @@ -536,83 +613,6 @@ decade. For aquatic ecosystems downstream, this is a loss of late-season cold-water input to streams during the warmest, most thermally stressful weeks of the year. -## Recent vs Pre-warming - -The table below compares two windows directly — the recent decade -(2015–2025) against the pre-warming reference (1951–1980). Δ -absolute is the difference of the two means in the variable's native -units. Δ % is shown only for variables where it is meaningful -(precipitation, soil moisture, vapour pressure deficit, relative -humidity). Two p-values are reported, answering different -questions. **Δ p (windows)** is the Welch t-test on the annual -values within each window — does the recent decade differ from -the pre-warming reference? **Trend p (75-yr)** is the Mann-Kendall -test on the full 1951–present series — is there a steady -year-on-year ramp? Step changes show up as significant Δ p with -non-significant trend p; gradual ramps show up as both significant. - -The recent decade was 1.6 to 1.8 °C warmer than the pre-warming -reference for annual mean, daytime maximum, and overnight minimum, -with both window and trend p-values below 0.001. Vapour pressure -deficit is up significantly on both tests. Annual precipitation -was about 6 to 7 % *lower* in the recent decade — and unlike the -FWCP Peace just to the north, the long-term Mann-Kendall trend -test does confirm a steady year-on-year decline (trend p ≈ 0.02). -Soil moisture is roughly flat. Relative humidity shows a small -significant decline. - -```{r compare-recipe, eval = FALSE} -cmp <- cd::cd_compare(ts) # defaults: 2015–2025 vs 1951–1980, Welch t-test -``` - -```{r compare-table, echo = FALSE} -# Recompute cmp from the cached ts so the p_value column reflects the -# current cd_compare() defaults — the precomputed rds was generated -# before #43 added the test argument. -cmp <- cd::cd_compare(ts) -trn_p <- trn[trn$trend_start == 1951, c("variable", "period", "mk_pvalue")] -names(trn_p)[3] <- "trend_p" -trn_p$trend_p <- round(trn_p$trend_p, 3) - -no_pct_vars <- c( - "tmean", "tmax", "tmin", - # snow_cover and snowfall_fraction are already in % (pct_point_diff - # anomaly type) — pct-of-baseline mixes units. snowmelt_doy_50 is a - # day-of-year — pct-change of a date is meaningless. - "snow_cover", "snowfall_fraction", "snowmelt_doy_50" -) -cmp_combined <- cmp |> - dplyr::mutate( - pct_change = ifelse( - variable %in% no_pct_vars, - NA_real_, - round(100 * (mean_a - mean_b) / mean_b, 1) - ), - mean_a = round(mean_a, 2), - mean_b = round(mean_b, 2), - abs_diff = round(difference, 2), - window_p = round(p_value, 3) - ) |> - dplyr::select(variable, period, mean_a, mean_b, abs_diff, pct_change, window_p) |> - merge(trn_p, by = c("variable", "period"), all.x = TRUE) - -names(cmp_combined) <- c( - "Variable", "Period", - "Recent (2015–2025)", "Pre-warming (1951–1980)", - "Δ absolute", "Δ %", "Δ p (windows)", "Trend p (75-yr)" -) - -kableExtra::kable_styling( - knitr::kable(cmp_combined, label = NA, - caption = "Recent decade (2015-2025 mean) compared to pre-warming reference (1951-1980 mean) for the Kootenay Lake Region.", - row.names = FALSE), - bootstrap_options = c("striped", "hover", "condensed") -) |> - kableExtra::scroll_box(height = "360px") -``` - -
- ## Spatial Pattern The zonal mean reduces a region this size to a single number; the spatial diff --git a/vignettes/peace-fwcp.Rmd b/vignettes/peace-fwcp.Rmd index ac11e84..6bb5381 100644 --- a/vignettes/peace-fwcp.Rmd +++ b/vignettes/peace-fwcp.Rmd @@ -253,6 +253,81 @@ cd::cd_plot_timeseries( ) ``` +## Recent Decade vs Pre-Warming Reference + +The table below compares two windows directly — the recent decade +(2015–2025) against the pre-warming reference (1951–1980). Δ +absolute is the difference of the two means in the variable's native +units. Δ % is shown only for variables where it is meaningful +(precipitation, soil moisture, vapour pressure deficit, relative +humidity). Two p-values are reported, answering different +questions. **Δ p (windows)** is the Welch t-test on the annual +values within each window — does the recent decade differ from +the pre-warming reference? **Trend p (75-yr)** is the Mann-Kendall +test on the full 1951–present series — is there a steady +year-on-year ramp? Step changes show up as significant Δ p with +non-significant trend p; gradual ramps show up as both significant. + +The recent decade was 1.7 to 2.4 °C warmer than the pre-warming +reference for annual mean, daytime maximum, and overnight minimum, +with both window and trend p-values below 0.001. Vapour pressure +deficit is up significantly on both tests. Annual precipitation +was about 3 to 4 % higher in the recent decade; neither test +confirms the shift (Δ p ≈ 0.4, trend p ≈ 0.20). Soil moisture and +relative humidity show no meaningful change. + +```{r compare-recipe, eval = FALSE} +cmp <- cd::cd_compare(ts) # defaults: 2015–2025 vs 1951–1980, Welch t-test +``` + +```{r compare-table, echo = FALSE} +# Recompute cmp from the cached ts so the p_value column reflects the +# current cd_compare() defaults — the precomputed rds was generated +# before #43 added the test argument. +cmp <- cd::cd_compare(ts) +trn_p <- trn[trn$trend_start == 1951, c("variable", "period", "mk_pvalue")] +names(trn_p)[3] <- "trend_p" +trn_p$trend_p <- round(trn_p$trend_p, 3) + +no_pct_vars <- c( + "tmean", "tmax", "tmin", + # snow_cover and snowfall_fraction are already in % (pct_point_diff + # anomaly type) — pct-of-baseline mixes units. snowmelt_doy_50 is a + # day-of-year — pct-change of a date is meaningless. + "snow_cover", "snowfall_fraction", "snowmelt_doy_50" +) +cmp_combined <- cmp |> + dplyr::mutate( + pct_change = ifelse( + variable %in% no_pct_vars, + NA_real_, + round(100 * (mean_a - mean_b) / mean_b, 1) + ), + mean_a = round(mean_a, 2), + mean_b = round(mean_b, 2), + abs_diff = round(difference, 2), + window_p = round(p_value, 3) + ) |> + dplyr::select(variable, period, mean_a, mean_b, abs_diff, pct_change, window_p) |> + merge(trn_p, by = c("variable", "period"), all.x = TRUE) + +names(cmp_combined) <- c( + "Variable", "Period", + "Recent (2015–2025)", "Pre-warming (1951–1980)", + "Δ absolute", "Δ %", "Δ p (windows)", "Trend p (75-yr)" +) + +kableExtra::kable_styling( + knitr::kable(cmp_combined, label = NA, + caption = "Recent decade (2015-2025 mean) compared to pre-warming reference (1951-1980 mean) for the FWCP Peace Region.", + row.names = FALSE), + bootstrap_options = c("striped", "hover", "condensed") +) |> + kableExtra::scroll_box(height = "360px") +``` + +
+ ## Daytime Highs and Overnight Lows The cd package ships daytime maximum (tmax) and overnight minimum (tmin) @@ -510,81 +585,6 @@ decade. For aquatic ecosystems downstream, this is a loss of late- season cold-water input to streams during the warmest part of the year. -## Recent vs Pre-warming - -The table below compares two windows directly — the recent decade -(2015–2025) against the pre-warming reference (1951–1980). Δ -absolute is the difference of the two means in the variable's native -units. Δ % is shown only for variables where it is meaningful -(precipitation, soil moisture, vapour pressure deficit, relative -humidity). Two p-values are reported, answering different -questions. **Δ p (windows)** is the Welch t-test on the annual -values within each window — does the recent decade differ from -the pre-warming reference? **Trend p (75-yr)** is the Mann-Kendall -test on the full 1951–present series — is there a steady -year-on-year ramp? Step changes show up as significant Δ p with -non-significant trend p; gradual ramps show up as both significant. - -The recent decade was 1.7 to 2.4 °C warmer than the pre-warming -reference for annual mean, daytime maximum, and overnight minimum, -with both window and trend p-values below 0.001. Vapour pressure -deficit is up significantly on both tests. Annual precipitation -was about 3 to 4 % higher in the recent decade; neither test -confirms the shift (Δ p ≈ 0.4, trend p ≈ 0.20). Soil moisture and -relative humidity show no meaningful change. - -```{r compare-recipe, eval = FALSE} -cmp <- cd::cd_compare(ts) # defaults: 2015–2025 vs 1951–1980, Welch t-test -``` - -```{r compare-table, echo = FALSE} -# Recompute cmp from the cached ts so the p_value column reflects the -# current cd_compare() defaults — the precomputed rds was generated -# before #43 added the test argument. -cmp <- cd::cd_compare(ts) -trn_p <- trn[trn$trend_start == 1951, c("variable", "period", "mk_pvalue")] -names(trn_p)[3] <- "trend_p" -trn_p$trend_p <- round(trn_p$trend_p, 3) - -no_pct_vars <- c( - "tmean", "tmax", "tmin", - # snow_cover and snowfall_fraction are already in % (pct_point_diff - # anomaly type) — pct-of-baseline mixes units. snowmelt_doy_50 is a - # day-of-year — pct-change of a date is meaningless. - "snow_cover", "snowfall_fraction", "snowmelt_doy_50" -) -cmp_combined <- cmp |> - dplyr::mutate( - pct_change = ifelse( - variable %in% no_pct_vars, - NA_real_, - round(100 * (mean_a - mean_b) / mean_b, 1) - ), - mean_a = round(mean_a, 2), - mean_b = round(mean_b, 2), - abs_diff = round(difference, 2), - window_p = round(p_value, 3) - ) |> - dplyr::select(variable, period, mean_a, mean_b, abs_diff, pct_change, window_p) |> - merge(trn_p, by = c("variable", "period"), all.x = TRUE) - -names(cmp_combined) <- c( - "Variable", "Period", - "Recent (2015–2025)", "Pre-warming (1951–1980)", - "Δ absolute", "Δ %", "Δ p (windows)", "Trend p (75-yr)" -) - -kableExtra::kable_styling( - knitr::kable(cmp_combined, label = NA, - caption = "Recent decade (2015-2025 mean) compared to pre-warming reference (1951-1980 mean) for the FWCP Peace Region.", - row.names = FALSE), - bootstrap_options = c("striped", "hover", "condensed") -) |> - kableExtra::scroll_box(height = "360px") -``` - -
- ## Spatial Pattern The zonal mean reduces a region this size to a single number; the spatial