男女別の自殺に至る要因について

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

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

FIG_DIR = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

# すべての図に添える注記（本ページの図はすべて報告値の可視化）
NOTE = '※原論文（河又杏香 2021）の報告値の可視化であり、本スクリプトで再計算した値ではない。'
print('セットアップ完了。図の出力先:', FIG_DIR)

# ── 原論文 表6〜10：男女別・要因別の自殺率との相関係数（報告値の転記）──
# 「男性有意」「女性有意」= 原論文で * が付されたもの（|r|>0.2 かつ p<0.05）
rows = [
    # (指標, 要因分類, 男性r, 女性r, 男性有意, 女性有意)
    ('降水日数',            '自然環境',  0.1892666,  -0.2720383, False, False),
    ('降水量',              '自然環境', -0.1264948,  -0.3876141, False, True),
    ('森林・人工林率',       '自然環境',  0.2891686,  -0.1450716, False, False),
    ('平均月収',            '社会的',   -0.3681953,   0.2598133, True,  False),
    ('最低賃金',            '社会的',   -0.4985958,   0.2394296, True,  False),
    ('人口密度',            '社会的',   -0.3484553,   0.1531512, True,  False),
    ('教師一人あたり生徒数','社会的',  -0.3917334,   0.2162066, True,  False),
    ('完全失業率',          '社会的',    0.08647537,  0.03220396, False, False),
    ('カルシウム',          '食生活',   -0.124765,    0.3667838, False, True),
    ('ビタミンD',           '食生活',   -0.1002074,   0.05206311,False, False),
    ('ビタミンA',           '食生活',   -0.22498,     0.2489419, False, False),
    ('酒類',               '食生活',   -0.3887164,   0.01988522,True,  False),
    ('離婚率',              '身辺事情', -0.09331021, -0.02504564,False, False),
    ('父子・母子世帯率',    '身辺事情',  0.2260748,  -0.2988605, False, True),
    ('入院患者率',          '身辺事情',  0.1125827,  -0.332586,  False, True),
    ('外来患者率',          '身辺事情', -0.0273847,  -0.2530422, False, False),
    ('学業',               '生活行動', -0.2988691,  -0.07702816,True,  False),
    ('仕事',               '生活行動',  0.2714662,  -0.08670267,False, False),
    ('介護関係',            '生活行動',  0.009313488,-0.1951387, False, False),
    ('趣味・娯楽の総数',    '生活行動', -0.5175879,   0.2080845, True,  False),
    ('スポーツ',            '生活行動', -0.560301,    0.1939965, True,  False),
    ('睡眠',               '生活行動',  0.6107483,  -0.06961662,True,  False),
    ('休養・くつろぎ',      '生活行動',  0.2736665,  -0.1304227, False, False),
    ('旅行・行楽',          '生活行動', -0.4192268,   0.2343289, True,  False),
]
rep = pd.DataFrame(rows, columns=['指標','要因分類','男性','女性','男性有意','女性有意'])
print('報告された相関係数（24指標）：')
print(rep[['要因分類','指標','男性','女性']].to_string(index=False))
print()
print('有意(p<0.05)な指標数  —  男性:', int(rep['男性有意'].sum()),
      ' 女性:', int(rep['女性有意'].sum()))

# ── 図1：24指標×男女別の相関係数（表6〜10の報告値）──
order = rep.iloc[::-1].reset_index(drop=True)   # 上から降水日数になるよう逆順
y = np.arange(len(order))
fig, axes = plt.subplots(1, 2, figsize=(13, 10), sharey=True)
for ax, sex in zip(axes, ['男性', '女性']):
    vals = order[sex].values
    sig  = order[sex + '有意'].values
    colors = []
    for v, s in zip(vals, sig):
        if s:
            colors.append('#C62828' if v >= 0 else '#1565C0')   # 濃色=有意
        else:
            colors.append('#EF9A9A' if v >= 0 else '#90CAF9')   # 淡色=非有意
    ax.barh(y, vals, color=colors, edgecolor='white')
    ax.axvline(0, color='#333', lw=1.1)
    ax.axvline(0.2, color='#bbb', ls=':', lw=0.8)
    ax.axvline(-0.2, color='#bbb', ls=':', lw=0.8)
    ax.set_xlim(-0.72, 0.72)
    ax.set_title('%s自殺率との相関' % sex, fontsize=13, fontweight='bold')
    ax.set_xlabel('相関係数 r（報告値）', fontsize=11)
    for yi, (v, s) in enumerate(zip(vals, sig)):
        if s:
            ax.text(v + (0.02 if v >= 0 else -0.02), yi, '*', va='center',
                    ha='left' if v >= 0 else 'right', fontsize=15, fontweight='bold', color='#222')
    ax.grid(True, axis='x', alpha=0.25)
axes[0].set_yticks(y)
axes[0].set_yticklabels(order['指標'].tolist(), fontsize=10)
fig.suptitle('図1　５要因２４指標と自殺率の相関係数（原論文 表6〜10 の報告値の可視化）',
             fontsize=14, fontweight='bold')
fig.text(0.5, 0.005, NOTE + '  濃色=有意(p<0.05, *)、淡色=非有意。点線は|r|=0.2。',
         ha='center', fontsize=9, color='#B71C1C')
fig.tight_layout(rect=[0, 0.035, 1, 0.965])
fig.savefig(os.path.join(FIG_DIR, '2021_H5_1_fig1.png'), bbox_inches='tight')
plt.close(fig)
male_sig = rep.loc[rep['男性有意'], '指標'].tolist()
female_sig = rep.loc[rep['女性有意'], '指標'].tolist()
print('fig1 saved')
print('男性で有意:', '、'.join(male_sig))
print('女性で有意:', '、'.join(female_sig))

# ── 図2：男性自殺率の単純相関 vs 偏相関（原論文 表7・11・12）──
# 最低賃金または人口密度の影響を除去した偏相関で疑似相関を判定する
mp = pd.DataFrame([
    ('森林・人工林率',       '人口密度',  0.2891686,  0.09068578),
    ('教師一人あたり生徒数','人口密度', -0.3917334, -0.2357294),
    ('最低賃金',            '人口密度', -0.4985958, -0.3938469),
    ('平均月収',            '最低賃金', -0.3681953,  0.09068578),
    ('人口密度',            '最低賃金', -0.3484553, -0.1100749),
    ('酒類',               '最低賃金', -0.3887164, -0.1681852),
    ('学業',               '最低賃金', -0.2988691, -0.2637391),
    ('趣味・娯楽の総数',    '最低賃金', -0.5175879, -0.2218496),
    ('スポーツ',            '最低賃金', -0.560301,  -0.3359372),
    ('睡眠',               '最低賃金',  0.6107483,  0.4310617),
    ('旅行・行楽',          '最低賃金', -0.4192268, -0.1382047),
], columns=['指標','除去変数','単純相関','偏相関'])
mp['真の相関'] = mp['偏相関'].abs() > 0.2   # |偏相関|>0.2 は疑似相関でない

fig, ax = plt.subplots(figsize=(11, 8.5))
yy = np.arange(len(mp))[::-1]
h = 0.38
ax.barh(yy + h/2, mp['単純相関'], height=h, color='#B0BEC5', label='単純相関 r')
pcol = ['#2E7D32' if t else '#C62828' for t in mp['真の相関']]
ax.barh(yy - h/2, mp['偏相関'], height=h, color=pcol)
ax.axvline(0, color='#333', lw=1.1)
ax.axvline(0.2, color='#bbb', ls=':', lw=0.8); ax.axvline(-0.2, color='#bbb', ls=':', lw=0.8)
ax.set_yticks(yy)
ax.set_yticklabels(['%s\n（除去: %s）' % (i, r) for i, r in zip(mp['指標'], mp['除去変数'])], fontsize=9)
ax.set_xlabel('相関係数 / 偏相関係数（報告値）', fontsize=11)
ax.set_xlim(-0.72, 0.72)
handles = [mpatches.Patch(color='#B0BEC5', label='単純相関 r'),
           mpatches.Patch(color='#2E7D32', label='偏相関（|r|>0.2：疑似相関でない）'),
           mpatches.Patch(color='#C62828', label='偏相関（|r|≦0.2：疑似相関）')]
ax.legend(handles=handles, loc='lower left', fontsize=9, framealpha=0.9)
ax.set_title('図2　男性自殺率：単純相関と偏相関の比較（原論文 表7・11・12 の報告値）',
             fontsize=12.5, fontweight='bold')
ax.grid(True, axis='x', alpha=0.25)
fig.text(0.5, 0.005, NOTE, ha='center', fontsize=9, color='#B71C1C')
fig.tight_layout(rect=[0, 0.03, 1, 1])
fig.savefig(os.path.join(FIG_DIR, '2021_H5_1_fig2.png'), bbox_inches='tight')
plt.close(fig)
survive = mp.loc[mp['真の相関'], '指標'].tolist()
pseudo  = mp.loc[~mp['真の相関'], '指標'].tolist()
print('fig2 saved')
print('偏相関でも |r|>0.2（疑似相関でない）:', '、'.join(survive))
print('偏相関で |r|≦0.2（疑似相関と判定）:', '、'.join(pseudo))
print('人口密度と最低賃金の相関 r = 0.8108688（強い正の相関・報告値）')

# ── 図3：女性自殺率の単純相関 vs 偏相関（病院数の影響を除去・原論文 表9・15）──
byoin_r, byoin_p = -0.3854302, 0.00746   # 人口千人あたり病院数×女性自殺率（追加指標）
fp = pd.DataFrame([
    ('降水量',           -0.3876141, -0.2364077),
    ('カルシウム',        0.3667838,  0.1793106),
    ('母子・父子世帯率', -0.2988605, -0.1149966),
    ('入院患者率',       -0.332586,   0.06085137),
], columns=['指標','単純相関','偏相関'])
fp['真の相関'] = fp['偏相関'].abs() > 0.2

fig, ax = plt.subplots(figsize=(11, 6))
yy = np.arange(len(fp))[::-1]
h = 0.38
ax.barh(yy + h/2, fp['単純相関'], height=h, color='#B0BEC5', label='単純相関 r')
pcol = ['#2E7D32' if t else '#C62828' for t in fp['真の相関']]
ax.barh(yy - h/2, fp['偏相関'], height=h, color=pcol)
# 追加指標：人口千人あたり病院数（除去に使用した変数自身）
ax.barh([len(fp) + 0.2], [byoin_r], height=0.5, color='#6A1B9A')
ax.text(byoin_r - 0.02, len(fp) + 0.2, 'r=%.3f\np=%.5f' % (byoin_r, byoin_p),
        va='center', ha='right', fontsize=9, color='#4A148C')
ax.axvline(0, color='#333', lw=1.1)
ax.axvline(0.2, color='#bbb', ls=':', lw=0.8); ax.axvline(-0.2, color='#bbb', ls=':', lw=0.8)
labels = fp['指標'].tolist()
ax.set_yticks(list(yy) + [len(fp) + 0.2])
ax.set_yticklabels(labels + ['人口千人あたり病院数（追加・除去に使用）'], fontsize=9)
ax.set_xlabel('相関係数 / 偏相関係数（報告値）', fontsize=11)
ax.set_xlim(-0.62, 0.62)
handles = [mpatches.Patch(color='#B0BEC5', label='単純相関 r'),
           mpatches.Patch(color='#2E7D32', label='偏相関（|r|>0.2）'),
           mpatches.Patch(color='#C62828', label='偏相関（|r|≦0.2：疑似相関）'),
           mpatches.Patch(color='#6A1B9A', label='追加指標：病院数×女性自殺率')]
ax.legend(handles=handles, loc='lower left', fontsize=9, framealpha=0.9)
ax.set_title('図3　女性自殺率：単純相関と偏相関（病院数の影響を除去）（原論文 表9・15 の報告値）',
             fontsize=11.5, fontweight='bold')
ax.grid(True, axis='x', alpha=0.25)
fig.text(0.5, 0.005, NOTE, ha='center', fontsize=9, color='#B71C1C')
fig.tight_layout(rect=[0, 0.04, 1, 1])
fig.savefig(os.path.join(FIG_DIR, '2021_H5_1_fig3.png'), bbox_inches='tight')
plt.close(fig)
print('fig3 saved')
print('人口千人あたり病院数×女性自殺率: r=%.4f, p=%.5f（有意）' % (byoin_r, byoin_p))
print('病院数除去後に |r|>0.2 が残ったのは:', '、'.join(fp.loc[fp['真の相関'], '指標'].tolist()))

d_raw = pd.read_csv('data/raw/SSDSE-D-2023.csv', encoding='cp932', header=1)
d_raw = d_raw[(d_raw['男女の別'] == '0_総数') & (d_raw['地域コード'] != 'R00000')]
d = d_raw.set_index('都道府県')
dv = lambda c: pd.to_numeric(d[c], errors='coerce')
b_raw = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', header=1)
b = b_raw[(b_raw['年度'] == 2021) & b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].set_index('都道府県')
bv = lambda c: pd.to_numeric(b[c], errors='coerce')

sleep = dv('睡眠')
hobby = dv('趣味・娯楽の総数')
divorce = bv('離婚件数') / bv('総人口') * 1000
r1, p1 = stats.pearsonr(sleep, divorce.reindex(sleep.index))
r2, p2 = stats.pearsonr(hobby, divorce.reindex(hobby.index))
print(f'睡眠時間 × 人口千対離婚件数 : r = {r1:+.3f} (p={p1:.4f})')
print(f'趣味・娯楽行動者率 × 離婚件数: r = {r2:+.3f} (p={p2:.4f})')
print('趣味・娯楽行動者率 下位5県:', '、'.join(hobby.nsmallest(5).index))

# ── 図4：男女で異なる自殺の関連要因（有意な相関のまとめ・報告値）──
theme = {
    '平均月収':'経済・労働', '最低賃金':'経済・労働', '人口密度':'経済・労働', '教師一人あたり生徒数':'経済・労働',
    '学業':'生活行動', '趣味・娯楽の総数':'生活行動', 'スポーツ':'生活行動', '睡眠':'生活行動', '旅行・行楽':'生活行動',
    '酒類':'食生活', 'カルシウム':'食生活',
    '降水量':'自然環境',
    '入院患者率':'健康・医療', '父子・母子世帯率':'健康・医療', '人口千人あたり病院数':'健康・医療',
}
tcol = {'経済・労働':'#1565C0', '生活行動':'#2E7D32', '食生活':'#E65100', '自然環境':'#607D8B', '健康・医療':'#C62828'}

m = rep[rep['男性有意']][['指標','男性']].rename(columns={'男性':'r'}).copy()
w = rep[rep['女性有意']][['指標','女性']].rename(columns={'女性':'r'}).copy()
w = pd.concat([w, pd.DataFrame([{'指標':'人口千人あたり病院数','r':-0.3854302}])], ignore_index=True)
for d in (m, w):
    d.sort_values('r', inplace=True)
    d.reset_index(drop=True, inplace=True)

fig, axes = plt.subplots(1, 2, figsize=(13, 6))
for ax, d, ttl in zip(axes, (m, w), ('男性で有意な相関', '女性で有意な相関')):
    yy = np.arange(len(d))
    cols = [tcol.get(theme.get(i, '自然環境'), '#999') for i in d['指標']]
    ax.barh(yy, d['r'], color=cols, edgecolor='white')
    ax.axvline(0, color='#333', lw=1.1)
    ax.set_yticks(yy); ax.set_yticklabels(d['指標'].tolist(), fontsize=10)
    ax.set_xlim(-0.65, 0.65)
    ax.set_xlabel('相関係数 r（報告値）', fontsize=11)
    ax.set_title(ttl, fontsize=13, fontweight='bold')
    ax.grid(True, axis='x', alpha=0.25)
handles = [mpatches.Patch(color=c, label=k) for k, c in tcol.items()]
fig.legend(handles=handles, loc='lower center', ncol=5, fontsize=9, framealpha=0.9, bbox_to_anchor=(0.5, -0.02))
fig.suptitle('図4　男女で異なる自殺の関連要因（有意な相関のまとめ・報告値）',
             fontsize=14, fontweight='bold')
fig.text(0.5, -0.06, NOTE, ha='center', fontsize=9, color='#B71C1C')
fig.tight_layout(rect=[0, 0.02, 1, 0.96])
fig.savefig(os.path.join(FIG_DIR, '2021_H5_1_fig4.png'), bbox_inches='tight')
plt.close(fig)
print('fig4 saved')
print('男性側：経済・労働系と生活行動系が中心')
print('女性側：健康・医療系（病院数・入院患者率）が中心')