外国人人口と市区町村の特性との関係性
重回帰・変数選択・決定木・ランダムフォレストによる探索的分析

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.tree import DecisionTreeRegressor, plot_tree
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

FIG_DIR = 'html/figures'
DATA_A = 'data/raw/SSDSE-A-2025.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ===== ステップ1: SSDSE-A(市区町村)から原論文 表1 の変数を構成 =====
# 1行目=変数コード、2行目=年度、3行目=日本語名 → コード行を列名に使う
raw = pd.read_csv(DATA_A, encoding='cp932', header=0, skiprows=[1, 2])
raw = raw.rename(columns={raw.columns[0]: 'code',
                          'Prefecture': '都道府県',
                          'Municipality': '市区町村'})
num = raw.apply(pd.to_numeric, errors='coerce')
print(f"【SSDSE-A-2025】{len(raw)}市区町村")
print('  原論文は SSDSE-2019A(2015年国勢調査ベース)+e-Stat追加4変数の29変数。')
print('  SSDSE-A-2025 では 未婚者割合・課税対象所得・学歴割合(大卒)・')
print('  刑法犯認知件数 の4変数が未収録 → 25変数で部分再現する')

df = pd.DataFrame({'都道府県': raw['都道府県'], '市区町村': raw['市区町村']})
pop = num['A1101'].replace(0, np.nan)  # 総人口(2020)。人口0の町村は除外

# 目的変数: 外国人人口(10万人あたり) = log(外国人/総人口*100000 + 1)
df['外国人人口'] = np.log(num['A1700'] / pop * 100000 + 1)

# 説明変数(原論文 表1 の計算式に従う。単位は標準化するため影響しない)
df['総人口'] = np.log(pop)
df['自然人口増減率'] = (num['A4101'] - num['A4200']) / pop * 100
df['社会的人口増減率'] = (num['A5101'] - num['A5102']) / pop * 100
df['核家族世帯割合'] = num['A810102'] / num['A710101'] * 100
df['単独世帯割合'] = num['A810105'] / num['A710101'] * 100
df['人口密度'] = np.log(pop / num['B1103'])          # 可住地面積あたり
emp = {'農業林業': 'C220832', '建設業': 'C220836', '製造業': 'C220837',
       '情報通信業': 'C220839', '運輸郵便業': 'C220840',
       '卸売小売業': 'C220841', '学術研究専門技術': 'C220846',
       '宿泊飲食サービス業': 'C220847', '教育学習支援業': 'C220849',
       '医療福祉': 'C220850', 'サービス業他': 'C220852'}
for name, code in emp.items():   # 従業者割合(外国人労働者が多い11産業)
    df[f'従業者割合({name})'] = num[code] / num['C2208'] * 100
df['農家数'] = np.log(num['C310201'] / pop * 100000 + 1)  # 販売農家
df['経常収支比率'] = num['D2203']
df['土木費割合'] = num['D320308'] / num['D3203'] * 100
df['教育費割合'] = num['D320310'] / num['D3203'] * 100
df['完全失業率'] = num['F1107'] / (num['F1107'] + num['F1102']) * 100
df['図書館数'] = np.log(num['G1401'] / pop * 100000 + 1)
df['医師数'] = np.log(num['I6100'] / pop * 100000 + 1)
df['保育所等数'] = np.log(num['J250302'] / pop * 100000 + 1)

# 総人口ゼロ(2020年国勢調査の福島県双葉町)と欠損値のある行を除外
df = df.dropna().reset_index(drop=True)
X_cols = [c for c in df.columns if c not in ('都道府県', '市区町村', '外国人人口')]
print(f"\n分析対象: {len(df)}市区町村 / 説明変数 {len(X_cols)}個")

# 原論文と同じく全変数を平均0・分散1に標準化
Z = df[X_cols + ['外国人人口']].copy()
Z = (Z - Z.mean()) / Z.std(ddof=0)

# ===== ステップ2: 重回帰(全変数)と残差プロット → 図1 =====
X_full = sm.add_constant(Z[X_cols])
ols_full = sm.OLS(Z['外国人人口'], X_full).fit()
stud = ols_full.get_influence().resid_studentized_external
out_mask = np.abs(stud) > 4      # 外れ値: スチューデント化残差 |t|>4
print('\n=== 残差プロットによる外れ値の確認(原論文 4.1節の手順) ===')
print(f"外れ値と判定: {out_mask.sum()}市区町村(|スチューデント化残差|>4)")
print(df.loc[out_mask, ['都道府県', '市区町村']].to_string(index=False))
print('  ※ 原論文はSSDSE-2019Aの残差プロットで6箇所を外れ値として除外')

fig1, ax1 = plt.subplots(figsize=(8, 6))
ax1.scatter(ols_full.fittedvalues, ols_full.resid, s=12, alpha=0.5,
            color='#3F51B5', label='市区町村')
ax1.scatter(ols_full.fittedvalues[out_mask], ols_full.resid[out_mask],
            s=140, facecolors='none', edgecolors='#C62828', linewidths=1.8,
            label=f'外れ値として除外({out_mask.sum()}件)')
ax1.axhline(0, color='gray', linewidth=1)
ax1.set_xlabel('予測値(標準化後)')
ax1.set_ylabel('残差')
ax1.set_title('図1: 重回帰(25変数)の残差プロット\n'
              '【SSDSE-A-2025 による実データ部分再現。原論文 図1 は'
              'SSDSE-2019A・29変数】', fontsize=11)
ax1.legend()
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2019_U5_2_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print('図1 保存完了(実データ部分再現)')

# 外れ値を除いて以降の分析に使う(原論文と同じく全手法で共通のデータ)
dfc = df.loc[~out_mask].reset_index(drop=True)
Zc = Z.loc[~out_mask].reset_index(drop=True)
y = Zc['外国人人口']
X = Zc[X_cols]

# ===== ステップ3: 重回帰の推定(実データ部分再現) =====
ols = sm.OLS(y, sm.add_constant(X)).fit()
coef = pd.DataFrame({'係数': ols.params, 'P値': ols.pvalues}).drop('const')
coef_sorted = coef.reindex(coef['係数'].abs().sort_values(ascending=False).index)
print('\n=== 重回帰の推定結果(実データ部分再現・係数の絶対値順) ===')
print(coef_sorted.round(4).to_string())
print(f"決定係数 R2 = {ols.rsquared:.4f} / N = {len(y)}")
print('  ※ 原論文(29変数)の報告値では 人口密度 0.375 が最大で、')
print('    製造業 0.287・課税対象所得 0.236・未婚者割合 0.196 が続く')

# ===== ステップ4: 図2 原論文 表2 の報告値の可視化(再計算ではない) =====
# 原論文 表2「重回帰分析の推定結果」の29変数の係数とP値。
# ※印は SSDSE-A-2025 未収録(e-Stat追加)の4変数。
# 「従業者割合(卸売業、小売業)」は原論文表2では行ラベルが
# 「宿泊、飲食サービス」と表記されているが、表1・表3の変数順との
# 対応から卸売業、小売業の誤記とみられる(値はそのまま用いる)。
paper_tab2 = [  # (変数名, 係数, P値)
    ('人口密度', 0.375057, 2.44e-16),
    ('従業者割合(製造業)', 0.287207, 2.69e-09),
    ('課税対象所得※', 0.236293, 7.95e-12),
    ('未婚者割合※', 0.196142, 7.82e-11),
    ('単独世帯割合', 0.156925, 1.93e-09),
    ('農家数', 0.141056, 6.04e-06),
    ('従業者割合(宿泊業、飲食サービス業)', 0.137586, 8.76e-06),
    ('完全失業率', -0.123417, 2.76e-08),
    ('従業者割合(農業、林業)', 0.120916, 1.18e-06),
    ('自然人口増減率', -0.10529, 0.000101),
    ('核家族世帯割合', -0.08347, 0.000763),
    ('従業者割合(医療、福祉)', -0.067113, 0.034773),
    ('刑法犯認知件数※', 0.059249, 0.008059),
    ('総人口', 0.050128, 0.101313),
    ('歳出における教育費割合', 0.047338, 0.007494),
    ('従業者割合(情報通信業)', -0.047043, 0.036935),
    ('図書館数', 0.043784, 0.009019),
    ('従業者割合(運輸業、郵便業)', 0.04116, 0.053491),
    ('従業者割合(建設業)', -0.038802, 0.178217),
    ('従業者割合(教育、学習支援業)', -0.038387, 0.08657),
    ('経常収支比率', -0.032771, 0.103782),
    ('社会的人口増減率', -0.03223, 0.105888),
    ('学歴割合(大卒)※', -0.026888, 0.497624),
    ('保育所等数', 0.025764, 0.150709),
    ('従業者割合(サービス業他)', -0.011944, 0.573534),
    ('従業者割合(卸売業、小売業)', 0.010875, 0.696603),
    ('従業者割合(学術研究、専門・技術)', -0.005469, 0.78124),
    ('医師数', 0.005705, 0.792476),
    ('歳出における土木費割合', 0.002082, 0.902164),
]
labels = [t[0] for t in paper_tab2][::-1]
betas = [t[1] for t in paper_tab2][::-1]
pvals = [t[2] for t in paper_tab2][::-1]
colors = ['#C62828' if (p < 0.05 and b > 0) else
          '#1565C0' if (p < 0.05 and b < 0) else '#B0BEC5'
          for b, p in zip(betas, pvals)]
fig2, ax2 = plt.subplots(figsize=(9, 10))
ax2.barh(range(len(labels)), betas, color=colors)
ax2.set_yticks(range(len(labels)))
ax2.set_yticklabels(labels, fontsize=9)
ax2.axvline(0, color='gray', linewidth=1)
ax2.set_xlabel('標準化データでの回帰係数(原論文 表2 の報告値)')
ax2.set_title('図2: 重回帰分析(29変数)の推定係数\n'
              '【原論文 表2 の報告値の可視化(再計算ではない)。'
              '赤=正で有意・青=負で有意(P<0.05)・灰=非有意】\n'
              '※印は SSDSE-A-2025 未収録の e-Stat 追加変数', fontsize=11)
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2019_U5_2_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print('\n図2 保存完了(原論文 表2 の報告値の可視化)')

# ===== ステップ5: Best Subset Selection の結果(原論文の報告) =====
print('\n=== Best Subset Selection(原論文 4.2節の報告) ===')
print('  原論文はRのleapsで全29変数からCp基準の変数選択を行い、')
print('  23変数のモデルが最良(Cp最小)。除外されたのは 歳出の土木費割合・')
print('  医師数・大卒者割合・従業者割合(卸売、学術研究、サービス) の6つ。')
print('  Cpプロット(原論文 図2)は本スクリプトでは再現しない(原論文参照)')

# ===== ステップ6: 決定木回帰 → 図3(実データ部分再現) =====
tree = DecisionTreeRegressor(max_leaf_nodes=9, random_state=0).fit(X, y)
fig3, ax3 = plt.subplots(figsize=(13, 7))
plot_tree(tree, feature_names=X_cols, filled=True, rounded=True,
          fontsize=8, impurity=False, ax=ax3)
ax3.set_title('図3: 決定木回帰(25変数・葉9個)\n'
              '【SSDSE-A-2025 による実データ部分再現。原論文 図3 の木'
              '(最初の分割は課税対象所得)とは変数構成が異なる】',
              fontsize=11)
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2019_U5_2_fig3.png'), bbox_inches='tight')
plt.close(fig3)
used = [X_cols[i] for i in sorted(set(tree.tree_.feature)) if i >= 0]
print('\n=== 決定木回帰(実データ部分再現) ===')
print(f"分割に使われた変数: {'、'.join(used)}")
print('  ※ 原論文の木は 課税対象所得 → 大卒割合・未婚割合 の順に分割し、')
print('    ほかに製造業・医療福祉・農家数・完全失業率・核家族世帯割合を使用。')
print('    課税対象所得・大卒割合・未婚割合は SSDSE-A-2025 未収録')
print('図3 保存完了(実データ部分再現)')

# ===== ステップ7: ランダムフォレスト → 図4(実データ部分再現) =====
rf = RandomForestRegressor(n_estimators=500, random_state=0,
                           n_jobs=-1).fit(X, y)
imp = pd.Series(rf.feature_importances_, index=X_cols).sort_values()
fig4, ax4 = plt.subplots(figsize=(9, 8))
ax4.barh(range(len(imp)), imp.values, color='#2E7D32')
ax4.set_yticks(range(len(imp)))
ax4.set_yticklabels(imp.index, fontsize=9)
ax4.set_xlabel('変数重要度(不純度減少に基づく)')
ax4.set_title('図4: ランダムフォレストの変数重要度(25変数)\n'
              '【SSDSE-A-2025 による実データ部分再現。原論文 図4 は'
              'R の RandomForest・29変数で、課税対象所得が最重要】',
              fontsize=11)
fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2019_U5_2_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print('\n=== ランダムフォレストの変数重要度(実データ部分再現・上位8) ===')
print(imp.sort_values(ascending=False).head(8).round(4).to_string())
print('  ※ 原論文(29変数)では 課税対象所得 が最重要、未婚割合、')
print('    従業者割合(製造業、医療・福祉)、犯罪認知件数などが続く')
print('図4 保存完了(実データ部分再現)')

# ===== ステップ8: 10-fold クロスバリデーションで手法を比較 =====
cv = KFold(n_splits=10, shuffle=True, random_state=1)
models = {'重回帰分析': LinearRegression(),
          '決定木分析': DecisionTreeRegressor(max_leaf_nodes=9,
                                              random_state=0),
          'ランダムフォレスト': RandomForestRegressor(
              n_estimators=500, random_state=0, n_jobs=-1)}
print('\n=== 10-fold クロスバリデーション: 平均二乗誤差の平均値 ===')
print('【SSDSE-A-2025・25変数による実データ部分再現】')
for name, model in models.items():
    mse = -cross_val_score(model, X, y, cv=cv,
                           scoring='neg_mean_squared_error').mean()
    print(f"  {name:<12s}: {mse:.7f}")
print('【原論文 表4 の報告値(SSDSE-2019A・29変数、R による)】')
print('  重回帰分析: 0.4686977 / Best Subset Selection: 0.4641824')
print('  決定木分析: 0.527091  / ランダムフォレスト: 0.4305073')
print('  → 原論文でも本再現でもランダムフォレストが最小(=最も妥当)。')
print('    Best Subset Selection は本再現では実行していない(原論文参照)')