学習意欲を高めるために
～日本の教育を見直す～

# ── SSDSE-D（生活行動・生活時間、都道府県別）の読み込み ─────────────────
# cp932・2行目に日本語見出し（header=1）。男女計（0_総数）かつ47都道府県だけを残す。
d = pd.read_csv(DATA_D, encoding='cp932', header=1)
d = d[d['男女の別'].astype(str).str.contains('総数')]
d = d[d['都道府県'] != '全国'].copy()
YCOL = '学習・自己啓発・訓練の総数'   # ＝学習意欲の代理指標（行動者率 %）
d[YCOL] = pd.to_numeric(d[YCOL], errors='coerce')
print('SSDSE-D 読み込み: N =', len(d), '都道府県  （調査年: 2021年）')
print(f'学習・自己啓発・訓練の行動者率: 平均 {d[YCOL].mean():.1f}%  '

# ── 原論文の報告値（2016年・PDF本文より転記） ───────────────────────────
# 学習・自己啓発・訓練の行動者率 と 各指標 の相関係数ρ（符号・数値は原論文どおり）
REPORTED = {
    '安定した生活': {
        'エンゲル係数':   0.2457423,   # p=0.09591（無相関）
        'スポーツ時間':   0.5453607,
        '睡眠時間':      -0.7728923,
    },
    '教育': {
        '学業時間':      -0.5725335,
        '教育費':         0.6677363,
        '教養娯楽費':     0.6740558,
    },
    '趣味・娯楽': {
        '読書':           0.8688196,
        '美術鑑賞':       0.7690495,
        '演芸・演劇・舞踊': 0.8742667,
        '音楽鑑賞':       0.8826958,
        '楽器の演奏':     0.8826888,
        'テレビ等':      -0.7092702,
        '遊園地等の見物':  0.8179552,
        'ゲーム':         0.7908306,
        '旅行・行楽':     0.6752246,
        '図書館数':      -0.3794661,   # p=0.008523（100万人当たり）
    },
    '自然環境': {
        '森林面積割合':   -0.5211408,
        '田の地積割合':   -0.1162438,   # p=0.4365（無相関）
        '畑の地積割合':   -0.1106442,   # p=0.4591（無相関）
    },
}
# 無相関性を棄却できなかった（原論文）指標
NS_ITEMS = {'エンゲル係数', '田の地積割合', '畑の地積割合'}

print()
print('=== 原論文の報告値（学習意欲との相関係数ρ・2016年） ===')
for theme, items in REPORTED.items():
    print(f'[{theme}]')
    for name, r in items.items():
        mark = '（無相関）' if name in NS_ITEMS else ''

# SSDSE-D の列名 → 原論文の指標名 と 報告値
D_COLS = {
    '趣味としての読書（マンガを除く）':                 ('読書',           0.8688196),
    '美術鑑賞(テレビ・スマートフォン・パソコンなどは除く)': ('美術鑑賞',       0.7690495),
    '演芸・演劇・舞踊鑑賞(テレビ・スマートフォン・パソコンなどは除く)': ('演芸・演劇・舞踊', 0.8742667),
    'ＣＤ・スマートフォンなどによる音楽鑑賞':            ('音楽鑑賞',       0.8826958),
    '楽器の演奏':                                     ('楽器の演奏',     0.8826888),
    '遊園地、動植物園、水族館などの見物':               ('遊園地等の見物',  0.8179552),
    '旅行・行楽の総数':                               ('旅行・行楽',     0.6752246),
}
print()
print('=== 実再現：SSDSE-D(2021年調査)で相関を再計算し報告値と比較 ===')
print(f"{'指標':<12}{'報告値(2016)':>14}{'再計算(2021)':>14}{'p値':>10}")
recompute = {}
for col, (name, rep) in D_COLS.items():
    x = pd.to_numeric(d[col], errors='coerce')
    m = x.notna() & d[YCOL].notna()
    r, p = stats.pearsonr(x[m], d[YCOL][m])
    recompute[name] = r
    print(f'{name:<12}{rep:>14.4f}{r:>14.4f}{p:>10.4f}')
print('→ 趣味・娯楽の行動者率は調査年が違っても強い正の相関がおおむね再現される。')

# 図2【実再現】読書 × 学習意欲（SSDSE-D 2021）
# ══════════════════════════════════════════════════════════════
xcol = '趣味としての読書（マンガを除く）'
x = pd.to_numeric(d[xcol], errors='coerce')
m = x.notna() & d[YCOL].notna()
r2, p2 = stats.pearsonr(x[m], d[YCOL][m])
slope, intercept, _, _, _ = stats.linregress(x[m], d[YCOL][m])
print()
print('=== 図2 読書 散布図：再計算 r ===')
print(f'読書行動者率 × 学習・自己啓発・訓練の行動者率: r = {r2:.4f} (p = {p2:.2e}, N = {m.sum()})')

# 図3【報告値の可視化】学級規模：総人口の交絡と偏相関
# ══════════════════════════════════════════════════════════════
# 原論文の報告値：
#   corr(総人口, 教員一人当たり児童・生徒数) … 交絡の大きさ（全て正）
#   partial corr(教員当たり児童・生徒数, 学習意欲 | 総人口) … 偏相関（弱い正）
levels = ['小学校', '中学校', '高等学校']
confound = [0.7324618, 0.7046322, 0.7825836]   # 総人口との相関
partial  = [0.4840228, 0.5989611, 0.5718934]   # 偏相関（p小学=3.238e-09）
print()
print('=== 学級規模：総人口の交絡と偏相関（原論文の報告値） ===')
print(f"{'学校種':<8}{'総人口との相関':>16}{'偏相関(学習意欲)':>18}")
for lv, c, pr in zip(levels, confound, partial):
    print(f'{lv:<8}{c:>16.7f}{pr:>18.7f}')

fig3, ax3 = plt.subplots(figsize=(8.5, 5.2))
xpos = np.arange(len(levels)); w = 0.38
b1 = ax3.bar(xpos - w/2, confound, w, color='#b0bec5', alpha=0.95,
             edgecolor='white', label='総人口との相関（交絡）')
b2 = ax3.bar(xpos + w/2, partial,  w, color='#5cb85c', alpha=0.95,
             edgecolor='white', label='偏相関：学習意欲との関係（総人口を除去）')
for b, v in list(zip(b1, confound)) + list(zip(b2, partial)):
    ax3.text(b.get_x() + b.get_width()/2, v + 0.01, f'{v:.3f}',
             ha='center', va='bottom', fontsize=9.5)
ax3.set_xticks(xpos); ax3.set_xticklabels(levels, fontsize=11)
ax3.set_ylim(0, 1)
ax3.set_ylabel('相関係数', fontsize=11)
ax3.set_title('図3: 教員一人当たり児童・生徒数と学習意欲【報告値の可視化】\n'
              '※ 再計算ではなく原論文の報告値。総人口の影響を除くと弱い正の偏相関',
              fontsize=12, fontweight='bold', pad=10)
ax3.legend(loc='upper right', fontsize=9.5)
ax3.grid(axis='y', alpha=0.3, linestyle='--')
ax3.set_facecolor('#FAFAFA')
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2021_H5_3_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print('図3 保存完了（報告値の可視化）')

# 図4【実再現】学習・自己啓発・訓練の行動者率 都道府県ランキング
# （原論文の図1（地図ヒートマップ）を、実データのランキングとして再表現）
# ══════════════════════════════════════════════════════════════
rank = d[['都道府県', YCOL]].dropna().sort_values(YCOL, ascending=True)
nat = d[YCOL].mean()
print()
print('=== 図4 学習・自己啓発・訓練 行動者率 都道府県ランキング（実再現・2021） ===')
top = rank.tail(5).iloc[::-1]
bot = rank.head(5)
print('▼ 上位5:')
for _, r in top.iterrows():
    print(f'  {r["都道府県"]:<6} {r[YCOL]:.1f}%')
print('▼ 下位5:')
for _, r in bot.iterrows():
    print(f'  {r["都道府県"]:<6} {r[YCOL]:.1f}%')
print(f'全国平均（47都道府県の平均）: {nat:.1f}%')

fig4, ax4 = plt.subplots(figsize=(8.5, 11))
cmap = plt.cm.YlOrRd
norm = plt.Normalize(rank[YCOL].min(), rank[YCOL].max())
bar_colors = cmap(norm(rank[YCOL].values))
ax4.barh(range(len(rank)), rank[YCOL], color=bar_colors,
         edgecolor='white', linewidth=0.4)
ax4.set_yticks(range(len(rank)))