from __future__ import annotations import hashlib import json import math from collections.abc import Iterable from dataclasses import dataclass from pathlib import Path from typing import Any import numpy as np import pandas as pd from ..evaluation_trace import ( EvaluationGateResult, EvaluationMetric, EvaluationRunRecord, ) from ..models import utc_now from ..sandbox import sha256_file from ..weather import WeatherStudy, load_weather_research_registry from .weather_conditions import build_game_weather_summaries ANALYSIS_VERSION = "1.2.0" START_SEASON = 2018 BOOTSTRAP_REPETITIONS = 4_000 MEANINGFUL_WIND_MPS = 6.7056 # 15 mph CALM_WIND_MPS = 4.4704 # 10 mph HEAVY_RAIN_MM_PER_HOUR = 2.5 PBP_COLUMNS = { "game_id", "qb_dropback", "qb_kneel", "qb_spike", "qb_scramble", "rush_attempt", "pass_attempt", "complete_pass", "air_yards", "epa", "wp", "score_differential", "down", "game_seconds_remaining", "aborted_play", "play_deleted", "field_goal_attempt", "field_goal_result", "kick_distance", "extra_point_attempt", "extra_point_result", "interception", "fumble_lost", } OUTCOMES: dict[str, tuple[str, str]] = { "neutral_dropback_rate": ("Neutral early-down dropback rate", "share"), "rushing_rate": ("Neutral early-down rushing rate", "share"), "completion_rate": ("Completion rate", "share"), "pass_epa_per_dropback": ("Pass EPA per dropback", "epa-per-dropback"), "average_depth_of_target": ("Average target depth", "yards"), "deep_target_rate": ("Deep-target rate", "share"), "deep_completion_rate": ("Deep completion rate", "share"), "deep_target_epa": ("EPA per deep target", "epa-per-target"), "turnover_rate": ("Turnover rate", "turnovers-per-play"), "field_goal_attempts": ("Combined field-goal attempts", "attempts-per-game"), "field_goal_accuracy": ("Field-goal accuracy", "share"), "kicker_fantasy_points": ("Combined kicker fantasy points", "points-per-game"), "total_points": ("Combined points", "points-per-game"), } @dataclass(frozen=True) class StudyAnalysisDesign: outcomes: tuple[str, ...] weather_controls: tuple[str, ...] = () fixed_effects: tuple[str, ...] = () @dataclass(frozen=True) class StudyContrast: contrast_id: str treatment: pd.Series control: pd.Series description: str condition_label: str control_label: str STUDY_DESIGNS = { "weather-kickers-roof": StudyAnalysisDesign( outcomes=("field_goal_attempts", "field_goal_accuracy", "kicker_fantasy_points") ), "weather-rain-playcalling": StudyAnalysisDesign( outcomes=( "neutral_dropback_rate", "average_depth_of_target", "completion_rate", "turnover_rate", ), weather_controls=("wind_speed_mps", "air_temp_c"), ), "weather-deep-targets-rain": StudyAnalysisDesign( outcomes=("deep_target_rate", "deep_completion_rate", "deep_target_epa"), weather_controls=("wind_speed_mps", "air_temp_c"), ), "weather-rain-wind-interaction": StudyAnalysisDesign( outcomes=( "neutral_dropback_rate", "pass_epa_per_dropback", "field_goal_accuracy", "total_points", ), weather_controls=("wind_speed_mps", "air_temp_c"), ), "weather-snow-offense": StudyAnalysisDesign( outcomes=( "neutral_dropback_rate", "rushing_rate", "pass_epa_per_dropback", "field_goal_accuracy", ), weather_controls=("wind_speed_mps", "air_temp_c"), ), "weather-snow-wind-interaction": StudyAnalysisDesign( outcomes=( "neutral_dropback_rate", "average_depth_of_target", "pass_epa_per_dropback", "total_points", ), weather_controls=("wind_speed_mps", "air_temp_c"), ), "weather-kickers-wind": StudyAnalysisDesign( outcomes=("field_goal_attempts", "field_goal_accuracy", "kicker_fantasy_points"), # Wind defines the treatment. Adjusting for it again would erase the estimand. weather_controls=("air_temp_c",), ), "weather-deep-targets-wind": StudyAnalysisDesign( outcomes=("deep_target_rate", "deep_completion_rate", "deep_target_epa"), weather_controls=("air_temp_c",), ), "weather-domes-offense": StudyAnalysisDesign( outcomes=( "total_points", "neutral_dropback_rate", "pass_epa_per_dropback", "deep_target_rate", ), fixed_effects=("home_team", "away_team"), ), "weather-open-closed-roof": StudyAnalysisDesign( outcomes=( "total_points", "pass_epa_per_dropback", "field_goal_accuracy", "kicker_fantasy_points", ), fixed_effects=("stadium_id", "home_team", "away_team"), ), } STUDY_OUTCOMES = {study_id: design.outcomes for study_id, design in STUDY_DESIGNS.items()} STUDY_CONTRAST_REGISTRATIONS = { "weather-kickers-roof": ["dome/closed vs outdoor/open"], "weather-rain-playcalling": ["observed rain vs explicitly rain-free"], "weather-deep-targets-rain": ["2.5+ mm maximum hourly rain vs rain-free; snow excluded"], "weather-rain-wind-interaction": ["rain plus 15+ mph wind vs dry below 15 mph"], "weather-snow-offense": ["observed snow vs explicitly snow-free"], "weather-snow-wind-interaction": ["snow plus 15+ mph wind vs no snow below 15 mph"], "weather-kickers-wind": [ "10+ mph wind vs below 10 mph", "12+ mph wind vs below 10 mph", "15+ mph wind vs below 10 mph", ], "weather-deep-targets-wind": ["15+ mph wind vs below 10 mph"], "weather-domes-offense": ["dome/closed vs outdoor/open"], "weather-open-closed-roof": [ "open vs closed within retractable-roof stadiums observed in both states" ], } def _read_required(path: Path, columns: set[str]) -> pd.DataFrame: available = set(pd.read_csv(path, nrows=0).columns) missing = sorted(columns - available) if missing: raise ValueError(f"{path.name} is missing required columns: {missing}") return pd.read_csv(path, usecols=sorted(columns), low_memory=False) def _aggregate_pbp(path: Path, game_ids: set[str]) -> pd.DataFrame: plays = _read_required(path, PBP_COLUMNS) plays = plays[plays["game_id"].isin(game_ids)].copy() if plays.empty: return pd.DataFrame() valid = ~plays["play_deleted"].fillna(0).eq(1) & ~plays["aborted_play"].fillna(0).eq(1) dropback = ( valid & plays["qb_dropback"].fillna(0).eq(1) & ~plays["qb_kneel"].fillna(0).eq(1) & ~plays["qb_spike"].fillna(0).eq(1) ) designed_rush = ( valid & plays["rush_attempt"].fillna(0).eq(1) & ~plays["qb_scramble"].fillna(0).eq(1) & ~plays["qb_kneel"].fillna(0).eq(1) ) neutral = ( (dropback | designed_rush) & plays["wp"].between(0.2, 0.8) & plays["score_differential"].between(-14, 14) & plays["down"].isin([1, 2]) & plays["game_seconds_remaining"].gt(120) ) attempt = dropback & plays["pass_attempt"].fillna(0).eq(1) deep = attempt & plays["air_yards"].ge(20) field_goal = valid & plays["field_goal_attempt"].fillna(0).eq(1) field_goal_made = field_goal & plays["field_goal_result"].eq("made") extra_point_made = ( valid & plays["extra_point_attempt"].fillna(0).eq(1) & plays["extra_point_result"].eq("good") ) distance = pd.to_numeric(plays["kick_distance"], errors="coerce") field_goal_points = np.select( [field_goal_made & distance.ge(50), field_goal_made & distance.ge(40), field_goal_made], [5, 4, 3], default=0, ) turnover = valid & ( plays["interception"].fillna(0).eq(1) | plays["fumble_lost"].fillna(0).eq(1) ) working = plays.assign( is_play=(dropback | designed_rush).astype(int), is_neutral=neutral.astype(int), neutral_dropback=(neutral & dropback).astype(int), is_dropback=dropback.astype(int), is_attempt=attempt.astype(int), is_complete=(attempt & plays["complete_pass"].fillna(0).eq(1)).astype(int), target_depth=plays["air_yards"].where(attempt), pass_epa=plays["epa"].where(dropback), is_deep=deep.astype(int), deep_complete=(deep & plays["complete_pass"].fillna(0).eq(1)).astype(int), deep_epa=plays["epa"].where(deep), is_turnover=turnover.astype(int), is_field_goal=field_goal.astype(int), field_goal_made=field_goal_made.astype(int), field_goal_distance=distance.where(field_goal), kicker_points=field_goal_points + extra_point_made.astype(int), ) return working.groupby("game_id", as_index=False).agg( scrimmage_plays=("is_play", "sum"), neutral_plays=("is_neutral", "sum"), neutral_dropbacks=("neutral_dropback", "sum"), dropbacks=("is_dropback", "sum"), pass_attempts=("is_attempt", "sum"), completions=("is_complete", "sum"), air_yards_total=("target_depth", "sum"), pass_epa_total=("pass_epa", "sum"), deep_targets=("is_deep", "sum"), deep_completions=("deep_complete", "sum"), deep_epa_total=("deep_epa", "sum"), turnovers=("is_turnover", "sum"), field_goal_attempts=("is_field_goal", "sum"), field_goals_made=("field_goal_made", "sum"), average_field_goal_distance=("field_goal_distance", "mean"), kicker_fantasy_points=("kicker_points", "sum"), ) def build_weather_game_metrics( *, schedule_path: Path, pbp_paths: Iterable[Path], weather_primary: pd.DataFrame, latest_completed_season: int, ) -> pd.DataFrame: schedule_columns = { "game_id", "season", "game_type", "week", "roof", "stadium_id", "home_team", "away_team", "home_score", "away_score", } schedule = _read_required(schedule_path, schedule_columns) schedule = schedule[ schedule["season"].between(START_SEASON, latest_completed_season) & schedule["game_type"].eq("REG") ].copy() if schedule["game_id"].duplicated().any(): raise ValueError("The schedule repeats a game identifier") parts = [ aggregated for path in pbp_paths if not (aggregated := _aggregate_pbp(path, set(schedule["game_id"]))).empty ] if not parts: raise ValueError("No play-by-play rows joined to the completed schedule") plays = pd.concat(parts, ignore_index=True) if plays["game_id"].duplicated().any(): raise ValueError("Play-by-play inputs repeat games across files") metrics = schedule.merge(plays, on="game_id", how="inner", validate="one_to_one") metrics["total_points"] = metrics["home_score"] + metrics["away_score"] metrics["neutral_dropback_rate"] = metrics["neutral_dropbacks"] / metrics["neutral_plays"] metrics["rushing_rate"] = 1 - metrics["neutral_dropback_rate"] metrics["completion_rate"] = metrics["completions"] / metrics["pass_attempts"] metrics["pass_epa_per_dropback"] = metrics["pass_epa_total"] / metrics["dropbacks"] metrics["average_depth_of_target"] = metrics["air_yards_total"] / metrics["pass_attempts"] metrics["deep_target_rate"] = metrics["deep_targets"] / metrics["pass_attempts"] metrics["deep_completion_rate"] = metrics["deep_completions"] / metrics["deep_targets"] metrics["deep_target_epa"] = metrics["deep_epa_total"] / metrics["deep_targets"] metrics["turnover_rate"] = metrics["turnovers"] / metrics["scrimmage_plays"] metrics["field_goal_accuracy"] = metrics["field_goals_made"] / metrics["field_goal_attempts"] weather_columns = [ "game_id", "air_temp_c", "wind_speed_mps", "precip_game_window_mm", "precip_max_1h_mm", "rain_observed", "snow_observed", "nearest_report_minutes", ] metrics = metrics.merge( weather_primary[weather_columns], on="game_id", how="left", validate="one_to_one", ) metrics["roof_group"] = np.where( metrics["roof"].isin(["dome", "closed"]), "indoor", np.where(metrics["roof"].isin(["outdoors", "open"]), "outdoor", "unknown"), ) return metrics.sort_values("game_id").reset_index(drop=True) def _rng(label: str) -> np.random.Generator: return np.random.default_rng(int.from_bytes(hashlib.sha256(label.encode()).digest()[:8], "big")) def _bootstrap_difference( treatment: np.ndarray, control: np.ndarray, *, label: str ) -> tuple[float, float, float]: rng = _rng(label) differences: list[np.ndarray] = [] for count in range(0, BOOTSTRAP_REPETITIONS, 200): repetitions = min(200, BOOTSTRAP_REPETITIONS - count) treatment_mean = rng.choice( treatment, size=(repetitions, len(treatment)), replace=True ).mean(axis=1) control_mean = rng.choice(control, size=(repetitions, len(control)), replace=True).mean( axis=1 ) differences.append(treatment_mean - control_mean) samples = np.concatenate(differences) lower, upper = np.quantile(samples, [0.025, 0.975]) p_value = min( 1.0, 2 * min(float((samples <= 0).mean()), float((samples >= 0).mean())), ) return float(lower), float(upper), p_value def _condition_masks(metrics: pd.DataFrame, study_id: str) -> tuple[pd.Series, pd.Series, str]: """Return the strictest registered contrast used by a study's sample gate.""" if study_id == "weather-kickers-roof": return ( metrics["roof_group"].eq("indoor"), metrics["roof_group"].eq("outdoor"), "fixed/retractable roof recorded closed or dome vs outdoor/open", ) rain = metrics["rain_observed"].eq(True) dry = metrics["rain_observed"].eq(False) snow = metrics["snow_observed"].eq(True) no_snow = metrics["snow_observed"].eq(False) wind = metrics["wind_speed_mps"].ge(MEANINGFUL_WIND_MPS) calm = metrics["wind_speed_mps"].lt(MEANINGFUL_WIND_MPS) calm_below_10 = metrics["wind_speed_mps"].lt(CALM_WIND_MPS) if study_id in {"weather-kickers-wind", "weather-deep-targets-wind"}: return wind, calm_below_10, "observed 15+ mph wind vs below 10 mph" if study_id == "weather-domes-offense": return ( metrics["roof_group"].eq("indoor"), metrics["roof_group"].eq("outdoor"), "fixed/retractable roof recorded closed or dome vs outdoor/open", ) if study_id == "weather-open-closed-roof": eligible_stadiums = set( metrics.loc[metrics["roof"].isin(["open", "closed"])] .groupby("stadium_id")["roof"] .agg(lambda values: frozenset(values.dropna())) .loc[lambda values: values.map(lambda states: {"open", "closed"}.issubset(states))] .index ) in_scope = metrics["stadium_id"].isin(eligible_stadiums) return ( in_scope & metrics["roof"].eq("open"), in_scope & metrics["roof"].eq("closed"), "recorded open vs closed games within retractable-roof stadiums observed in both states", ) if study_id == "weather-rain-playcalling": return rain, dry, "observed rain vs explicitly rain-free" if study_id == "weather-deep-targets-rain": return ( metrics["precip_max_1h_mm"].ge(HEAVY_RAIN_MM_PER_HOUR) & rain & ~snow, dry & ~snow, f"maximum hourly rain at least {HEAVY_RAIN_MM_PER_HOUR:g} mm vs rain-free", ) if study_id == "weather-rain-wind-interaction": return rain & wind, dry & calm, "observed rain plus 15+ mph wind vs dry and below 15 mph" if study_id == "weather-snow-offense": return snow, no_snow, "observed snow vs explicitly snow-free" if study_id == "weather-snow-wind-interaction": return snow & wind, no_snow & calm, "observed snow plus 15+ mph wind vs no snow and calm" raise ValueError(f"No deterministic weather condition is implemented for {study_id}") def _study_contrasts(metrics: pd.DataFrame, study_id: str) -> list[StudyContrast]: if study_id == "weather-kickers-wind": control = metrics["wind_speed_mps"].lt(CALM_WIND_MPS) return [ StudyContrast( contrast_id=f"wind-{threshold_mph}-mph", treatment=metrics["wind_speed_mps"].ge(threshold_mph * 0.44704), control=control, description=f"observed {threshold_mph}+ mph wind vs below 10 mph", condition_label=f"{threshold_mph}+ mph wind", control_label="Below 10 mph", ) for threshold_mph in (10, 12, 15) ] treatment, control, description = _condition_masks(metrics, study_id) labels = { "weather-kickers-roof": ("Dome or closed roof", "Outdoor or open roof"), "weather-rain-playcalling": ("Observed rain", "Explicitly rain-free"), "weather-deep-targets-rain": ("Heavy observed rain", "Rain-free"), "weather-rain-wind-interaction": ("Rain plus 15+ mph wind", "Dry below 15 mph"), "weather-snow-offense": ("Observed snow", "Explicitly snow-free"), "weather-snow-wind-interaction": ("Snow plus 15+ mph wind", "No snow below 15 mph"), "weather-deep-targets-wind": ("15+ mph wind", "Below 10 mph"), "weather-domes-offense": ("Dome or closed roof", "Outdoor or open roof"), "weather-open-closed-roof": ("Open retractable roof", "Closed retractable roof"), }[study_id] return [ StudyContrast( contrast_id="primary", treatment=treatment, control=control, description=description, condition_label=labels[0], control_label=labels[1], ) ] def _holm(rows: list[dict[str, Any]]) -> None: ordered = sorted(enumerate(rows), key=lambda pair: pair[1]["p_value"]) running = 0.0 count = len(rows) adjusted = [1.0] * count for rank, (index, row) in enumerate(ordered): running = max(running, min(1.0, row["p_value"] * (count - rank))) adjusted[index] = running for row, value in zip(rows, adjusted, strict=True): row["holm_adjusted_p"] = value def _adjusted_effect( frame: pd.DataFrame, *, outcome: str, treatment_column: str, weather_controls: tuple[str, ...], fixed_effects: tuple[str, ...], ) -> tuple[float, float, float, float, str, float]: controls = ["week", "season", *weather_controls, *fixed_effects] if outcome == "field_goal_accuracy": controls.append("average_field_goal_distance") selected = frame.dropna(subset=[outcome, treatment_column, *controls]).copy() season = pd.get_dummies(selected["season"].astype(str), drop_first=True, dtype=float).to_numpy() columns = [ np.ones(len(selected)), selected[treatment_column].to_numpy(dtype=float), (selected["week"].to_numpy(dtype=float) - 9) / 10, ] labels = ["week and season fixed effects"] for weather_control in weather_controls: center = 15 if weather_control == "air_temp_c" else 0 scale = 10 if weather_control == "air_temp_c" else 5 columns.append((selected[weather_control].to_numpy(dtype=float) - center) / scale) if weather_controls: labels.append(" and ".join(name.replace("_", " ") for name in weather_controls)) if outcome == "field_goal_accuracy": columns.append((selected["average_field_goal_distance"].to_numpy(dtype=float) - 40) / 10) labels.append("average field-goal distance") columns.append(season) for fixed_effect in fixed_effects: columns.append( pd.get_dummies( selected[fixed_effect].astype(str), drop_first=True, dtype=float ).to_numpy() ) if fixed_effects: labels.append( " and ".join(name.replace("_", " ") for name in fixed_effects) + " fixed effects" ) design = np.column_stack(columns) values = selected[outcome].to_numpy(dtype=float) inverse = np.linalg.pinv(design.T @ design) coefficients = inverse @ design.T @ values residuals = values - design @ coefficients leverage = np.sum(design * (design @ inverse), axis=1) adjusted_residuals = residuals / np.clip(1 - leverage, 1e-8, None) covariance = inverse @ (design.T @ ((adjusted_residuals[:, None] ** 2) * design)) @ inverse standard_error = float(np.sqrt(max(0, covariance[1, 1]))) estimate = float(coefficients[1]) if standard_error == 0: robust_p_value = 0.0 if estimate != 0 else 1.0 else: # Two-sided normal approximation for the same HC3 estimate shown in the article. robust_p_value = math.erfc(abs(estimate / standard_error) / math.sqrt(2)) return ( estimate, standard_error, estimate - 1.96 * standard_error, estimate + 1.96 * standard_error, "; ".join(labels), robust_p_value, ) def analyze_conditional_weather( metrics: pd.DataFrame, studies: Iterable[WeatherStudy], *, data_through: int, ) -> dict[str, Any]: decisions: list[dict[str, Any]] = [] effects: list[dict[str, Any]] = [] observed_seasons = metrics.loc[metrics["nearest_report_minutes"].notna(), "season"] observed_weather_data_through = ( int(observed_seasons.max()) if not observed_seasons.empty else None ) for study in studies: if study.study_id not in STUDY_OUTCOMES: continue design = STUDY_DESIGNS[study.study_id] treatment_mask, control_mask, contrast = _condition_masks(metrics, study.study_id) treatment_games = int(treatment_mask.sum()) control_games = int(control_mask.sum()) complete_games = treatment_games + control_games eligible = ( complete_games >= study.minimum_complete_games and treatment_games >= study.minimum_condition_games ) decision = "eligible" if eligible else study.underpowered_policy reasons = [] if complete_games < study.minimum_complete_games: reasons.append( f"{complete_games} comparable games; {study.minimum_complete_games} required" ) if treatment_games < study.minimum_condition_games: reasons.append( f"{treatment_games} condition games; {study.minimum_condition_games} required" ) decisions.append( { "study_id": study.study_id, "title": study.title, "contrast": contrast, "condition_games": treatment_games, "control_games": control_games, "comparable_games": complete_games, "minimum_complete_games": study.minimum_complete_games, "minimum_condition_games": study.minimum_condition_games, "data_through": int(metrics.loc[treatment_mask | control_mask, "season"].max()), "decision": decision, "reasons": reasons, } ) study_rows: list[dict[str, Any]] = [] for registered_contrast in _study_contrasts(metrics, study.study_id): comparison_mask = registered_contrast.treatment | registered_contrast.control comparison = metrics[comparison_mask].copy() comparison["study_treatment"] = registered_contrast.treatment[comparison_mask].astype( int ) for outcome in design.outcomes: treatment = ( metrics.loc[registered_contrast.treatment, outcome] .dropna() .to_numpy(dtype=float) ) control = ( metrics.loc[registered_contrast.control, outcome].dropna().to_numpy(dtype=float) ) if len(treatment) < 2 or len(control) < 2: continue lower, upper, p_value = _bootstrap_difference( treatment, control, label=f"{study.study_id}:{registered_contrast.contrast_id}:{outcome}", ) adjusted = _adjusted_effect( comparison, outcome=outcome, treatment_column="study_treatment", weather_controls=design.weather_controls, fixed_effects=design.fixed_effects, ) label, units = OUTCOMES[outcome] study_rows.append( { "study_id": study.study_id, "publication_decision": decision, "contrast_id": registered_contrast.contrast_id, "contrast": registered_contrast.description, "condition_label": registered_contrast.condition_label, "control_label": registered_contrast.control_label, "contrast_comparable_games": int(comparison_mask.sum()), "outcome": outcome, "label": label, "units": units, "condition_games": len(treatment), "control_games": len(control), "condition_estimate": float(treatment.mean()), "control_estimate": float(control.mean()), "unadjusted_difference": float(treatment.mean() - control.mean()), "bootstrap_ci_95_lower": lower, "bootstrap_ci_95_upper": upper, "bootstrap_p_value": p_value, "p_value": adjusted[5], "adjusted_difference": adjusted[0], "robust_standard_error": adjusted[1], "adjusted_ci_95_lower": adjusted[2], "adjusted_ci_95_upper": adjusted[3], "controls": adjusted[4], } ) _holm(study_rows) effects.extend(study_rows) results = { "schema_version": "1.0", "analysis_version": ANALYSIS_VERSION, "data_through": ( str(observed_weather_data_through) if observed_weather_data_through else None ), "football_data_through": str(data_through), "observed_weather_data_through": ( str(observed_weather_data_through) if observed_weather_data_through else None ), "source_lag": ( { "source_key": "noaa-global-hourly", "requested_through": str(data_through), "available_through": str(observed_weather_data_through), "handling": ( "Weather-dependent studies stop at the latest available season; " "the next refresh retries automatically." ), } if observed_weather_data_through is not None and observed_weather_data_through < data_through else None ), "estimand": "Observed associations, not causal effects or pregame forecast value", "games_with_play_by_play": int(metrics["game_id"].nunique()), "games_with_observed_weather": int(metrics["nearest_report_minutes"].notna().sum()), "study_decisions": decisions, "result_sha256": "", } results["result_sha256"] = _records_hash([results | {"result_sha256": None}, *effects]) return {"results": results, "effects": effects} def _records_hash(value: Any) -> str: def canonical(item: Any) -> Any: if isinstance(item, dict): return {key: canonical(entry) for key, entry in item.items()} if isinstance(item, list): return [canonical(entry) for entry in item] if isinstance(item, (np.integer,)): return int(item) if isinstance(item, (float, np.floating)): return round(float(item), 12) if np.isfinite(item) else None return item encoded = json.dumps(canonical(value), sort_keys=True, separators=(",", ":"), allow_nan=False) return hashlib.sha256(encoded.encode()).hexdigest() def _write_json(path: Path, value: Any) -> None: path.write_text(json.dumps(value, indent=2, sort_keys=True, allow_nan=False), encoding="utf-8") def run_conditional_weather_studies( *, schedule_path: Path, pbp_paths: Iterable[Path], observation_dir: Path, output_dir: Path, artifact_root: str, latest_completed_season: int, code_sha: str, run_id: str, started_at: str, completed_at: str | None = None, ) -> EvaluationRunRecord: paths = tuple(pbp_paths) output_dir.mkdir(parents=True, exist_ok=False) registry = load_weather_research_registry() studies = [study for study in registry.studies if study.study_id in STUDY_OUTCOMES] registration = { "schema_version": "1.0", "program_id": registry.program_id, "analysis_version": ANALYSIS_VERSION, "code_sha": code_sha, "locked_at": started_at, "locked_before_results": True, "population": "NFL regular-season games from 2018 through the latest completed season", "observed_weather_window": "one hour before kickoff through four hours after kickoff", "meaningful_wind_mps": MEANINGFUL_WIND_MPS, "heavy_rain_mm_per_hour": HEAVY_RAIN_MM_PER_HOUR, "uncertainty": "4,000 game-level bootstrap draws and HC3 robust OLS intervals", "multiplicity": "Holm adjustment within each preregistered study family", "studies": [ { "study_id": study.study_id, "outcomes": list(STUDY_OUTCOMES[study.study_id]), "contrasts": STUDY_CONTRAST_REGISTRATIONS[study.study_id], "weather_controls": list(STUDY_DESIGNS[study.study_id].weather_controls), "fixed_effects": list(STUDY_DESIGNS[study.study_id].fixed_effects), "minimum_complete_games": study.minimum_complete_games, "minimum_condition_games": study.minimum_condition_games, "underpowered_policy": study.underpowered_policy, } for study in studies ], "interpretation": "association, not causal effect or pregame forecast skill", } _write_json(output_dir / "registration.json", registration) weather = build_game_weather_summaries( schedule_path=schedule_path, acquisition_dir=observation_dir, latest_completed_season=latest_completed_season, ) primary = pd.DataFrame(weather["primary"]) metrics = build_weather_game_metrics( schedule_path=schedule_path, pbp_paths=paths, weather_primary=primary, latest_completed_season=latest_completed_season, ) first = analyze_conditional_weather(metrics, studies, data_through=latest_completed_season) second = analyze_conditional_weather(metrics, studies, data_through=latest_completed_season) reproducible = first["results"]["result_sha256"] == second["results"]["result_sha256"] source_manifest = json.loads( (observation_dir / "source-manifest.json").read_text(encoding="utf-8") ) by_filename: dict[str, dict[str, Any]] = {} for record in source_manifest["sources"]: filename = record.get("filename") source_name = str(record.get("source_name", "")) if not filename and source_name.startswith("NOAA Global Hourly "): _, _, _, station_id, season = source_name.split() filename = f"{station_id}-{season}.csv" if filename: by_filename[str(filename)] = record raw_observation_files = sorted((observation_dir / "global-hourly").glob("*.csv")) provenance_valid = all( path.name in by_filename and by_filename[path.name].get("rights_status") == "approved" and by_filename[path.name].get("sha256") == sha256_file(path) for path in raw_observation_files ) coverage = pd.DataFrame(weather["coverage"]) completed_coverage = coverage[coverage["usable_time_games"] > 0] historical_coverage = float( completed_coverage["usable_time_games"].sum() / completed_coverage["games"].sum() ) decisions = first["results"]["study_decisions"] rain_ready = next( ( row["decision"] == "eligible" for row in decisions if row["study_id"] == "weather-rain-playcalling" ), False, ) rain_agreement = weather["station_agreement"]["rain_classification_agreement"] wind_comparable = weather["station_agreement"]["wind_comparable_games"] wind_pair_mae = weather["station_agreement"]["wind_mean_absolute_difference_mps"] gates = [ EvaluationGateResult( gate_id="rights-and-provenance", status="passed" if provenance_valid else "failed", details="Every used NOAA observation matches an approved source receipt and SHA-256.", evidence_hrefs=[f"{artifact_root}/source-manifest.json"], ), EvaluationGateResult( gate_id="observation-coverage", status="passed" if historical_coverage >= 0.9 else "failed", details=f"Released seasons have {historical_coverage:.1%} usable game-window coverage.", evidence_hrefs=[f"{artifact_root}/coverage.csv"], ), EvaluationGateResult( gate_id="paired-station-sensitivity", status="passed" if rain_agreement is not None and rain_agreement >= 0.9 else "failed", details="Primary and secondary stations meet the locked rain agreement floor.", evidence_hrefs=[f"{artifact_root}/station-agreement.json"], ), EvaluationGateResult( gate_id="paired-station-wind-sensitivity", status=( "passed" if wind_comparable is not None and wind_comparable >= 1_000 and wind_pair_mae is not None and wind_pair_mae <= 2.5 else "failed" ), details=( "At least 1,000 paired games are required with primary-secondary wind mean " "absolute difference no greater than 2.5 m/s." ), evidence_hrefs=[f"{artifact_root}/station-agreement.json"], ), EvaluationGateResult( gate_id="rain-study-sample", status="passed" if rain_ready else "failed", details="The main rain study has enough comparable and exposed games for review.", evidence_hrefs=[f"{artifact_root}/study-decisions.json"], ), EvaluationGateResult( gate_id="reproducibility", status="passed" if reproducible else "failed", details="Two independent executions produced the same canonical result hash.", evidence_hrefs=[f"{artifact_root}/reproducibility.json"], ), EvaluationGateResult( gate_id="temporal-separation", status="passed", details="Observed conditions are marked research-only and cannot enter predictions.", evidence_hrefs=[f"{artifact_root}/registration.json"], ), ] status = "passed" if all(gate.status == "passed" for gate in gates) else "withheld" primary.replace({np.nan: None}).to_csv(output_dir / "game-weather.csv", index=False) metrics.replace({np.nan: None}).to_csv(output_dir / "game-metrics.csv", index=False) pd.DataFrame(first["effects"]).to_csv(output_dir / "study-effects.csv", index=False) coverage.to_csv(output_dir / "coverage.csv", index=False) _write_json(output_dir / "station-agreement.json", weather["station_agreement"]) _write_json(output_dir / "study-decisions.json", {"studies": decisions}) _write_json(output_dir / "results.json", first["results"]) _write_json( output_dir / "reproducibility.json", { "first_sha256": first["results"]["result_sha256"], "second_sha256": second["results"]["result_sha256"], "matched": reproducible, }, ) _write_json(output_dir / "source-manifest.json", source_manifest) _write_json( output_dir / "gates.json", { "passed": status == "passed", "gates": [gate.model_dump(mode="json") for gate in gates], }, ) evaluation_metrics = [ EvaluationMetric( metric_id="released-weather-coverage", label="Released-season observed weather coverage", value=historical_coverage, units="share", population="US outdoor/open regular-season games through the latest NOAA release", ) ] if rain_agreement is not None: evaluation_metrics.append( EvaluationMetric( metric_id="rain-station-agreement", label="Rain classification station agreement", value=rain_agreement, units="share", population="Games with comparable primary and secondary station reports", ) ) if wind_pair_mae is not None: evaluation_metrics.append( EvaluationMetric( metric_id="wind-station-mean-absolute-difference", label="Primary-secondary wind mean absolute difference", value=wind_pair_mae, units="m/s", population="Games with comparable primary and secondary station wind readings", ) ) evaluation = EvaluationRunRecord( run_id=run_id, workload_id="weather-observed-conditions", run_kind="research", status=status, started_at=started_at, completed_at=completed_at or utc_now(), code_sha=code_sha, source_sha256={ "schedule": sha256_file(schedule_path), "observations": sha256_file(observation_dir / "source-manifest.json"), **{ f"pbp-{path.name.removeprefix('play_by_play_').removesuffix('.csv.gz')}": ( sha256_file(path) ) for path in paths }, }, registration_href=f"{artifact_root}/registration.json", artifact_hrefs=[ f"{artifact_root}/{name}" for name in ( "results.json", "game-weather.csv", "game-metrics.csv", "study-effects.csv", "study-decisions.json", "coverage.csv", "station-agreement.json", "source-manifest.json", "reproducibility.json", ) ], metrics=evaluation_metrics, gates=gates, total_cost_usd=0, ) _write_json(output_dir / "evaluation-run.json", evaluation.model_dump(mode="json")) return evaluation