"""
16_country_group_analysis.py
国别/区域分组分析模块：8大区域组的描述统计与组间比较
Country Group Analysis Module: descriptive statistics and between-group comparison
across the 8 OECD-38 region groups

用途 / Purpose:
    对 08 号脚本产出的 panel_re_scores.csv，按 oecd38_re_index_panel 表
    schema 中的 region_key 枚举（anglosphere / western_europe /
    latin_america / cee / nordic / southern_europe / other / east_asia）
    做分组统计：均值、标准差、2000-2023 变化量、组内排名。
    并做单因素 ANOVA 检验区域间 RE 均值差异是否显著。

依赖 / Dependencies: pandas, numpy, scipy
"""

import logging
from pathlib import Path

import numpy as np
import pandas as pd
from scipy.stats import f_oneway

logging.basicConfig(level=logging.INFO, format="%(asctime)s  %(message)s")
log = logging.getLogger("16_country_group")

DATA_DIR = Path("./data")

# ISO3 -> region_key 映射，须与 oecd38_re_index_panel 表 schema 的枚举一致
REGION_MAP = {
    "USA": "anglosphere", "GBR": "anglosphere", "AUS": "anglosphere",
    "CAN": "anglosphere", "NZL": "anglosphere", "IRL": "anglosphere",
    "DEU": "western_europe", "FRA": "western_europe", "NLD": "western_europe",
    "BEL": "western_europe", "AUT": "western_europe", "CHE": "western_europe",
    "LUX": "western_europe",
    "MEX": "latin_america", "CHL": "latin_america", "COL": "latin_america",
    "CRI": "latin_america",
    "POL": "cee", "CZE": "cee", "HUN": "cee", "SVK": "cee", "SVN": "cee",
    "EST": "cee", "LVA": "cee", "LTU": "cee",
    "SWE": "nordic", "NOR": "nordic", "DNK": "nordic", "FIN": "nordic",
    "ISL": "nordic",
    "ITA": "southern_europe", "ESP": "southern_europe", "PRT": "southern_europe",
    "GRC": "southern_europe",
    "TUR": "other", "ISR": "other", "COL2": "other",
    "JPN": "east_asia", "KOR": "east_asia",
}


def main():
    scores = pd.read_csv(DATA_DIR / "panel_re_scores.csv")
    scores["region_key"] = scores["iso3"].map(REGION_MAP)
    unmapped = scores[scores["region_key"].isna()]["iso3"].unique()
    if len(unmapped):
        log.warning("Countries without region mapping (check REGION_MAP): %s", list(unmapped))

    group_stats = (
        scores.groupby("region_key")["re_score"]
        .agg(mean="mean", std="std", min="min", max="max")
        .round(2)
        .reset_index()
    )

    # 2000 vs 2023 变化量，按区域汇总
    y0 = scores[scores["year"] == scores["year"].min()].groupby("region_key")["re_score"].mean()
    y1 = scores[scores["year"] == scores["year"].max()].groupby("region_key")["re_score"].mean()
    change = (y1 - y0).round(2).rename("re_score_change_2000_2023")
    group_stats = group_stats.merge(change, on="region_key", how="left")

    group_stats.to_csv(DATA_DIR / "region_group_statistics.csv", index=False)
    log.info("Region group statistics saved -> data/region_group_statistics.csv")

    # 单因素 ANOVA：检验区域间 RE 均值差异
    groups = [g["re_score"].values for _, g in scores.dropna(subset=["region_key"]).groupby("region_key")]
    if len(groups) > 1:
        f_stat, p_value = f_oneway(*groups)
        log.info("One-way ANOVA across region groups: F=%.3f, p=%.4g", f_stat, p_value)
        with open(DATA_DIR / "region_anova_result.txt", "w", encoding="utf-8") as f:
            f.write(f"One-way ANOVA (region_key groups on re_score)\nF-statistic: {f_stat:.4f}\np-value: {p_value:.6g}\n")


if __name__ == "__main__":
    main()