thermostrife.inference

Case-crossover null + conditional logit + stratified permutation, plus the descriptive H1 tests ported from the original analysis.py.

.. py:module:: thermostrife.inference

Case-crossover null model + H1 descriptive tests + 1σ rescaling.

The H1 tests (Wilcoxon signed-rank, sign, bootstrap CI) are ports of the original analysis.py scaffold kept under _legacy/; the H2 case-crossover engine is the headline contribution of this module. All tests run unconditionally on every event regardless of significance, following the project rule that decision logic belongs in interpretation, not in the pipeline.

.. py:function:: wilcoxon_signed_rank(anomaly: numpy.ndarray) -> dict :module: thermostrife.inference

One-sided Wilcoxon signed-rank that median(anomaly) > 0.

.. py:function:: sign_test(anomaly: numpy.ndarray) -> dict :module: thermostrife.inference

Binomial sign test on P(anomaly > 0) > 0.5.

.. py:function:: bootstrap_mean_ci(anomaly: np.ndarray, n_boot: int = 10000, alpha: float = 0.05, rng: np.random.Generator | None = None) -> dict :module: thermostrife.inference

Bootstrap percentile CI for the mean anomaly.

.. py:function:: daylight_hours(lat_deg: float, when: ~datetime.date) -> float :module: thermostrife.inference

Closed-form solar daylight hours for when at latitude lat_deg.

Uses the standard astronomical-twilight-free approximation (sunrise/sunset by solar declination); accurate to ~10 min outside the polar circles, which is well inside the noise of our daily-Tmax inference and avoids a hard dependency on astral.

:returns: Daylight in hours. 0.0 at polar night, 24.0 at polar day.

.. py:function:: build_case_crossover_frame(events: list[dict]) -> pandas.DataFrame :module: thermostrife.inference

Stack per-event case + control rows into one long DataFrame.

Each events element must contain:

  • event_id (str): stratum identifier.

  • lat (float), lon (float): for the daylight covariate.

  • when (date): event day.

  • tmax_event_c (float): event-day Tmax (the case row).

  • baseline (DataFrame): one row per control day, index = date, column tmax in °C.

Returns a long DataFrame with columns [event_id, day, is_case, tmax_c, daylight_h], one row per (event, day). Cases get is_case = 1; controls get 0.

.. py:function:: case_crossover_conditional_logit(frame: pandas.DataFrame, *, covariates: list[str] | None = None) -> dict :module: thermostrife.inference

Conditional logistic regression on matched case-control sets.

Fits

   logit P(case | event_id) = β · tmax_c + γ · covariates

using a per-event-id stratum (conditional likelihood); the stratum intercepts are integrated out. exp(β) is the odds ratio of uprising-on-this-day per +1 °C above the local same-month baseline.

:param frame: Output of :func:build_case_crossover_frame. :param covariates: Extra columns to include alongside tmax_c (default: ["daylight_h"]).

:returns: Dict with the headline estimates and metadata.

.. py:function:: stratified_permutation(frame: pd.DataFrame, n_perm: int = 10000, rng: np.random.Generator | None = None) -> dict :module: thermostrife.inference

Within-event label-shuffle test on the mean case-vs-control gap.

Under H0 the case-day Tmax is exchangeable with the baseline-day Tmax within the same event stratum. For each event we shuffle the case label across (case + controls), compute the (event-mean of case-Tmax) minus (event-mean of control-Tmax), and average over events.

Two-sided p = fraction of permutations whose absolute statistic is ≥ the observed statistic.

.. py:function:: stratify_case_crossover(events: list[dict], key_fn, *, min_events: int = 5, covariates: list[str] | None = None) -> dict :module: thermostrife.inference

Split events into strata by key_fn(event) -> str and run H2 per stratum.

Returns a dict mapping stratum label to the :func:case_crossover_conditional_logit result dict. Strata with fewer than min_events events are skipped with a {'skipped': True, 'reason': ...} marker so the caller can still see them in the report.

.. py:function:: event_anomaly_profile(events: list[dict], offsets: tuple[int, …] = (-7, -2, -1, 0, 1, 2, 7), fetch_fn=) -> pandas.DataFrame :module: thermostrife.inference

Fetch per-event Tmax at each offset and return as a long frame.

For each event we look up Tmax on day when + offset via the same source that resolved the event day, subtract the event’s baseline mean to get an anomaly, and stack the rows. Empty offsets (the fetcher returned None) are skipped silently — the superposed-epoch plotter handles the resulting ragged n per offset.

Returns columns: event_id, offset_days, anomaly_C.

.. py:function:: h3_within_event_contrast(events: list[dict], window_offsets: tuple[int, …] = (0, -1), surround_offsets: tuple[int, …] = (-7, 7), fetch_fn=) -> dict :module: thermostrife.inference

Per-event paired contrast: mean anomaly {t, t-1} minus {t-7, t+7}.

For each resolved event, fetch Tmax on day t + offset for each offset in window_offsets ∪ surround_offsets via the same source/station that resolved the event day, convert to anomaly by subtracting the existing baseline mean, then take

   diff = mean(anomaly @ window_offsets) − mean(anomaly @ surround_offsets)

A positive diff means the event day and its immediate neighbour sit higher above the local baseline than days a week away — the signature of a “hot day triggers riot” mechanism. A flat profile (diff ≈ 0) is consistent with the alternative “hot week happened to contain an event” explanation.

Tests the across-event distribution of diff with a one-sided Wilcoxon signed-rank against H₁: median > 0. Returns the n actually used (events where every offset resolved), the mean and median diff, a bootstrap 95 % CI on the mean, and the p-value.

.. py:function:: hsiang_sigma_rescaled(events: list[dict]) -> dict :module: thermostrife.inference

Per-event z-score (event - baseline_mean) / baseline_std, averaged.

Following Burke, Hsiang & Miguel (2015), expressing each event’s anomaly as a multiple of its own baseline-window standard deviation makes effects comparable across stations and climates. Their headline number for interpersonal violence is +2.4 % per 1 σ contemporaneous warming.

Returns the mean z across events, a bootstrap 95 % CI, and the fraction of events with z > 0.

.. py:function:: benjamini_hochberg(pvalues: dict[str, float], alpha: float = 0.05) -> dict :module: thermostrife.inference

Benjamini–Hochberg step-up FDR correction.

:param pvalues: Dict mapping test name → raw p-value. :param alpha: Target false-discovery rate (default 0.05).

:returns: Dict with per-test raw_p, bh_adjusted_p (the standard monotone-adjusted q-value), bonferroni_adjusted_p, and the boolean bh_reject / bonferroni_reject flags. Also exposes the family-level threshold pair so the report can quote them.