ソーシャルキャピタルと
地域コミュニティ活動の関係分析

region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東',
    '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄'
}
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12'
}

df['地域'] = df['都道府県'].map(region_map)

pop = df['総人口']  # 単位: 人

# 保育所密度（保育所等数 / 総人口 × 10,000）
df['保育所密度'] = df['保育所等数'] / pop * 10000

# 高齢化率（65歳以上 / 総人口）
df['高齢化率'] = df['65歳以上人口'] / pop * 100

# 年少人口率（15歳未満 / 総人口）
df['年少人口率'] = df['15歳未満人口'] / pop * 100

# 医療施設密度（一般病院数 + 診療所数 / 万人）
df['医療施設密度'] = (df['一般病院数'] + df['一般診療所数']) / pop * 10000

# リサイクル率（環境への意識 = ソーシャルキャピタルの代理）
df['リサイクル率'] = df['ごみのリサイクル率']

# 教育進学率（高校卒業→進学率）
df['進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100

# 婚姻率（人口千人あたり）
df['婚姻率'] = df['婚姻件数'] / pop * 1000

# 消費支出に占める教養娯楽費率（文化活動への関与）
df['教養娯楽費率'] = df['教養娯楽費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100

# 消費支出に占める教育費率
df['教育費率'] = df['教育費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100

# 旅行・観光活動（延べ宿泊者 / 万人）
df['観光密度'] = df['延べ宿泊者数'] / pop * 10000

# 保育所定員充足率（在所児 / 定員）
df['保育充足率'] = df['保育所等在所児数'] / df['保育所等定員数'] * 100

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 scipy import stats
from sklearn.preprocessing import StandardScaler
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import pdist

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

FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

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)

# 2022年度断面データ
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
assert len(df) == 47, f"47都道府県が揃っていません: {len(df)}件"

vars_all = [
    '合計特殊出生率', '保育所密度', '高齢化率', '年少人口率',
    '医療施設密度', 'リサイクル率', '進学率', '婚姻率',
    '教養娯楽費率', '教育費率', '観光密度', '保育充足率'
]
df_ana = df[['都道府県', '地域'] + vars_all].dropna().copy()
print(f"分析対象: {len(df_ana)}都道府県")

sc_vars = [
    '保育所密度', '高齢化率', 'リサイクル率', '進学率',
    '婚姻率', '教養娯楽費率', '教育費率', '医療施設密度'
]
corr_labels = [
    '保育所密度', '高齢化率', 'リサイクル率', '進学率',
    '婚姻率', '教養娯楽費率', '教育費率', '医療施設密度'
]
corr_matrix = df_ana[sc_vars].corr()

y_var = '合計特殊出生率'
x_vars = ['保育所密度', '高齢化率', 'リサイクル率', '進学率',
          '婚姻率', '教養娯楽費率', '医療施設密度']

y = df_ana[y_var]
X_raw = df_ana[x_vars]

# 標準化
scaler = StandardScaler()
X_std = scaler.fit_transform(X_raw)
X_std_df = pd.DataFrame(X_std, columns=x_vars, index=df_ana.index)
X_std_df = sm.add_constant(X_std_df)

ols_model = sm.OLS(y, X_std_df).fit()
print(ols_model.summary())

# 標準化偏回帰係数と95%CI
coef = ols_model.params[1:]       # 定数項を除く
ci   = ols_model.conf_int().iloc[1:]
pval = ols_model.pvalues[1:]

r2      = ols_model.rsquared
r2_adj  = ols_model.rsquared_adj
f_stat  = ols_model.fvalue
f_pval  = ols_model.f_pvalue

print(f"\nR² = {r2:.3f}, Adj R² = {r2_adj:.3f}")
print(f"F統計量 = {f_stat:.2f}, p値 = {f_pval:.4f}")

cluster_vars = ['保育所密度', '高齢化率', 'リサイクル率', '婚姻率', '医療施設密度', '進学率']
X_clust = df_ana[cluster_vars].values
X_clust_std = StandardScaler().fit_transform(X_clust)
Z = linkage(X_clust_std, method='ward')

# 個々の相関係数（fig1用）
r_birth_hoiku, p_birth_hoiku = stats.pearsonr(df_ana['保育所密度'], df_ana['合計特殊出生率'])
print(f"\n保育所密度 vs 合計特殊出生率: r={r_birth_hoiku:.3f}, p={p_birth_hoiku:.4f}")

fig, ax = plt.subplots(figsize=(14, 6))
pref_names = df_ana['都道府県'].tolist()

dn = dendrogram(Z, labels=pref_names, ax=ax,
                leaf_rotation=90, leaf_font_size=9,
                color_threshold=Z[-4, 2])

ax.set_title('図4: Ward法による47都道府県のコミュニティ特性クラスタリング（2022年度）',
             fontsize=13, pad=12)
ax.set_ylabel('クラスター間距離（Ward距離）', fontsize=11)
ax.set_xlabel('都道府県', fontsize=11)
ax.grid(True, axis='y', alpha=0.3)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_U5_5_fig4.png'), bbox_inches='tight')
plt.close(fig)
print("fig4 saved")

print("\n========== HTML用 統計サマリー ==========")
print(f"分析年度: 2022年度")
print(f"分析対象: {len(df_ana)}都道府県")
print(f"決定係数 R²: {r2:.3f}")
print(f"自由度調整済み R²: {r2_adj:.3f}")
print(f"F統計量: {f_stat:.2f}  p値: {f_pval:.4f}")
print(f"保育所密度 vs 合計特殊出生率: r = {r_birth_hoiku:.3f}, p = {p_birth_hoiku:.4f}")

print("\n-- 標準化偏回帰係数 --")
for v in x_vars:
    b = ols_model.params[v]
    p = ols_model.pvalues[v]
    ci_lo = ols_model.conf_int().loc[v, 0]
    ci_hi = ols_model.conf_int().loc[v, 1]
    sig = '*' if p < 0.05 else ''
    print(f"  {v}: β={b:.3f}  CI=[{ci_lo:.3f}, {ci_hi:.3f}]  p={p:.4f} {sig}")

print("\n-- 相関行列（ソーシャルキャピタル変数）--")
for i, v1 in enumerate(sc_vars):
    for j, v2 in enumerate(sc_vars):
        if j > i:
            r, p = stats.pearsonr(df_ana[v1], df_ana[v2])
            if abs(r) > 0.4:
                sig = '**' if p < 0.01 else ('*' if p < 0.05 else '')
                print(f"  {v1} ↔ {v2}: r={r:.3f} p={p:.4f} {sig}")

print("\n-- 都道府県別クラスター（閾値 = Z[-4,2]）--")
from scipy.cluster.hierarchy import fcluster
labels_clust = fcluster(Z, t=4, criterion='maxclust')
for i, pref in enumerate(pref_names):
    print(f"  {pref}: クラスター{labels_clust[i]}")

print("\n全図の保存先:", FIG_DIR)
print("スクリプト完了")