自治体間の人口流動性を考慮した
潜在的な人手不足の可視化

positioning = {
    '農業・林業': {
        'risk': ['北海道', '岩手県', '青森県', '長野県', '鹿児島県'],
        'grow': ['宮崎県', '秋田県', '島根県'],
    },
    '卸売業・小売業': {
        'risk': ['北海道', '鹿児島県', '熊本県', '和歌山県', '長崎県', '秋田県'],
        'grow': ['都市部', '沖縄県'],
    },
    '運輸業・郵便業': {
        'risk': ['北海道', '大阪府'],
        'grow': ['埼玉県', '千葉県'],
    },
    '情報通信業': {
        'risk': [],
        'grow': ['東京都が極端に上昇（一極集中）'],
    },
    '不動産業・物品賃貸業': {
        'risk': ['関西圏', '北海道', '沖縄県'],
        'grow': ['東京都が突出（一極集中）'],
    },
}
print('\n=== [2] 報告値の可視化：ポジショニングマップの第二象限＝人手不足リスク地域（原論文 図1）===')
print('※ 産業別従業者数(SSDSE-A)はSSDSE-Bに無いため再計算不可。以下はPDF原本からの転記。')
for ind, v in positioning.items():
    risk = '・'.join(v['risk']) if v['risk'] else '（該当なし）'
    print(f'  【{ind}】人手不足リスク: {risk}')
    print(f'            伸び・期待  : {"・".join(v["grow"])}')

factors = {
    '第2因子': {
        'ind': '情報通信業・金融業・不動産業/物品賃貸業（重要度が高く増加数も多い）',
        'pref': ['東京都'],
        'note': '都市部（大阪・愛知・宮城・広島・福岡・関東圏）も弱く反応＝東京一極集中の要因',
    },
    '第3因子': {
        'ind': '農業・林業（重要度が高いのに増加数が低い）',
        'pref': ['北海道', '青森県', '岩手県'],
        'note': '東北で農林業の潜在的人手不足リスク。近畿以西はリスクが低い',
    },
    '第10因子': {
        'ind': '卸売業・小売業（重要度が本研究で最大なのに増加数が低い）',
        'pref': ['北海道', '青森県', '岩手県', '福島県', '大阪府', '京都府', '高知県', '千葉県', '埼玉県'],
        'note': '東北＋西日本＋東京近郊（千葉・埼玉）で卸小売の潜在的人手不足リスク',
    },
}
print('\n=== [3] 報告値の可視化：NMF主要因子と反応が強い都道府県（原論文 表2・図2）===')
print('※ NMF（10因子）は産業別従業者数が入力。SSDSE-Bで再計算不可のためPDF原本から転記。')
for fac, v in factors.items():
    print(f'  【{fac}】{v["ind"]}')
    print(f'            反応が強い地域: {"・".join(v["pref"])}')
    print(f'            解釈          : {v["note"]}')

df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}', na=False)].copy()
df_b['年度'] = df_b['年度'].astype(int)
d = df_b[df_b['年度'] == 2023].copy()

d['転入超過数'] = d['転入者数（日本人移動者）'] - d['転出者数（日本人移動者）']
d['転入超過率'] = d['転入超過数'] / d['総人口'] * 1000          # 人口千人あたり
d['高齢化率']  = d['65歳以上人口'] / d['総人口'] * 100
d['地域区分']  = d['都道府県'].map(region_map)
d = d.sort_values('転入超過率', ascending=False).reset_index(drop=True)

net_in = d[d['転入超過数'] > 0]
print('=== [1] 実再現：転入超過（社会増減）の集計 ===')
print(f'対象: SSDSE-B 2023年 / {len(d)}都道府県（転入者数−転出者数, 日本人移動者）')
print(f'転入超過（社会増）の都道府県: {len(net_in)} / {len(d)}')
print('  →', '・'.join(net_in['都道府県'].tolist()))
print(f'転出超過（社会減）の都道府県: {len(d) - len(net_in)} / {len(d)}')
print('原論文の動機（2018年）: 東京圏＋愛知・大阪・福岡が転入超過。'
      '2023年でも東京圏＋大阪・福岡が転入超過で構図は共通（愛知は2023年は転出超過）。')

print('\n=== [4] 実再現：転入超過率ランキング（SSDSE-B 2023, 人口千人あたり）===')
print('［上位＝転入超過］')
for _, r in d.head(6).iterrows():
    print(f'  {r["都道府県"]:<5} 転入超過 {int(r["転入超過数"]):>+7,} 人 / 千人あたり {r["転入超過率"]:+.2f}')
print('［下位＝転出超過］')
for _, r in d.tail(6).iterrows():
    print(f'  {r["都道府県"]:<5} 転入超過 {int(r["転入超過数"]):>+7,} 人 / 千人あたり {r["転入超過率"]:+.2f}')

r_val, p_val = stats.pearsonr(d['転入超過率'], d['高齢化率'])
print('\n=== [5] 実再現：転入超過率 × 高齢化率 の相関 ===')
print(f'ピアソン相関 r = {r_val:.3f}  (p = {p_val:.2e}, N={len(d)})')

a = pd.read_csv('data/raw/SSDSE-A-2025.csv', encoding='cp932', header=2)
a = a[a['地域コード'].str.match(r'^R\d+', na=False)].copy()
a = a[pd.to_numeric(a['総人口'], errors='coerce') > 0]
num = lambda c: pd.to_numeric(a[c], errors='coerce')

d2 = pd.DataFrame({'市区町村': (a['都道府県'] + a['市区町村']).values})
d2['転出超過率'] = ((num('転出者数（日本人移動者）') - num('転入者数（日本人移動者）')) / num('総人口') * 100).values
d2['生産年齢割合'] = (num('15～64歳人口') / num('総人口') * 100).values
d2['第3次産業割合'] = (num('第３次産業就業者数') / num('就業者数') * 100).values
d2 = d2.dropna()
print(f'市区町村数: {len(d2)}')
print('転出超過率 上位5（人口流出＝潜在的人手不足の供給源）:')
for _, r in d2.nlargest(5, '転出超過率').iterrows():
    print(f"  {r['市区町村']}: 転出超過 {r['転出超過率']:.2f}%/年")
r1, p1 = stats.pearsonr(d2['転出超過率'], d2['生産年齢割合'])
print(f'転出超過率 × 生産年齢割合: r = {r1:+.3f} (p={p1:.1e})')