気温と脳卒中の発症リスクについて

import os
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Patch

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 = '原論文の報告値の可視化（再計算ではない）'

# ============================================================
# 原論文の報告値（PDF原本より転記。捏造なし）
# ============================================================

# ── 図1（原論文）2018年 日本人の死因割合（人口動態調査）──────────
cause_labels = ['悪性新生物\n（がん） 27%', '心疾患 15%',
                '脳血管疾患\n（脳卒中） 8%', 'その他 49%']
cause_values = [27, 15, 8, 49]
cause_colors = ['#4C72B0', '#2E5C9A', '#C62828', '#B0B0B0']
cause_explode = [0, 0, 0.10, 0]

# ── 表4（原論文）各モデル都市の相関係数（月間平均気温×脳卒中死亡率, 2019年）
#    北日本地域 = 札幌・仙台・新潟
corr_cities = ['札幌市', '仙台市', '新潟市', 'さいたま市',
               '大阪市', '広島市', '福岡市', '熊本市']
corr_values = [-0.80, -0.60, -0.43, -0.72, -0.66, -0.19, -0.63, -0.58]
corr_is_north = [True, True, True, False, False, False, False, False]

# ── 仮説3（原論文4.3節）相関・偏相関のまとめ ─────────────────
#   気温×脳卒中（床暖房の影響を除く前の相関）      : -0.33
#   気温×脳卒中（床暖房の影響を除いた後の偏相関）  : -0.40（仮説に反し強まった）
#   気温×床暖房普及率（相関）                       : -0.43
#   床暖房普及率×脳卒中（外気温の影響を除く偏相関） : -0.25（弱い負の相関, 1月）
pcorr_labels = ['気温×脳卒中\n（床暖房を除く前）',
                '気温×脳卒中\n（床暖房を除いた後・偏相関）',
                '気温×床暖房\n普及率',
                '床暖房普及率×脳卒中\n（気温を除く偏相関・1月）']
pcorr_values = [-0.33, -0.40, -0.43, -0.25]
pcorr_colors = ['#4C72B0', '#C62828', '#6A1B9A', '#2E7D32']

# ── 仮説の検証結果まとめ（原論文5節の解釈より）──────────────
hypotheses = [
    ('仮説1-1', '月間平均気温が低いほど脳卒中死亡率は高い', '支持',
     '8都市すべてで負の相関（表4）。ただし札幌市は例外的に死亡率が低い'),
    ('仮説1-2', '日較差（最高−最低気温差）が大きいほど死亡率は高い', '不支持',
     '札幌は負・仙台はほぼ無相関。1日の気温差は影響なし'),
    ('仮説2-1', '月間平均気温が低い年は死亡率が高い', '不支持',
     '各年1月で回帰直線がほぼ水平。年単位の絶対気温は無関係'),
    ('仮説2-2', '前年より気温が低い年は死亡率が高い', '不支持',
     '都市によりばらつき。一概には言えない'),
    ('仮説3-1', '床暖房普及率が高い都市ほど死亡率は低い', '支持（弱）',
     '外気温を除いた偏相関 −0.25。札幌の低さを説明'),
    ('仮説3-2', '床暖房を除くと気温と死亡率の相関は弱まる', '不支持',
     '偏相関は −0.33→−0.40 と逆に強化。室内要因を除いても相関は残る'),
]

# 図1: 2018年 日本人の死因割合（原論文 図1 の報告値）
# ============================================================
fig1, ax1 = plt.subplots(figsize=(8, 7))
wedges, texts, autotexts = ax1.pie(
    cause_values, labels=cause_labels, colors=cause_colors,
    explode=cause_explode, startangle=90, counterclock=False,
    autopct='', textprops={'fontsize': 13}, wedgeprops={'edgecolor': 'white', 'linewidth': 1.5}
)
ax1.set_title('図1：2018年 日本人の死因割合\n脳卒中（脳血管疾患）は全体の8%',
              fontsize=15, fontweight='bold', pad=18)
ax1.text(0, -1.35, f'出典: 人口動態調査（2018年）／{NOTE}',
         ha='center', va='center', fontsize=10, color='#666')
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_H2_fig1.png'), dpi=150, bbox_inches='tight')
plt.close()
print('図1（死因割合）保存完了')

# 図2: 各モデル都市の相関係数（原論文 表4 の報告値）
# ============================================================
fig2, ax2 = plt.subplots(figsize=(10, 6.5))
order = np.argsort(corr_values)          # 強い負相関を上に
cities_o = [corr_cities[i] for i in order]
vals_o   = [corr_values[i] for i in order]
north_o  = [corr_is_north[i] for i in order]
colors_o = ['#C62828' if n else '#4C72B0' for n in north_o]

y = np.arange(len(cities_o))
bars = ax2.barh(y, vals_o, color=colors_o, edgecolor='white', height=0.68)
for yi, v in zip(y, vals_o):
    ax2.text(v - 0.02, yi, f'{v:.2f}', va='center', ha='right',
             fontsize=12, color='white', fontweight='bold')
ax2.set_yticks(y)
ax2.set_yticklabels(cities_o, fontsize=12)
ax2.set_xlim(-0.9, 0.05)
ax2.axvline(0, color='#333', lw=1)
ax2.set_xlabel('相関係数（月間平均気温 × 脳卒中死亡率, 2019年）', fontsize=12)
ax2.set_title('図2：各モデル都市の相関係数（原論文 表4）\nすべての都市で負の相関＝気温が低いほど脳卒中死亡率は高い',
              fontsize=13, fontweight='bold')
ax2.legend(handles=[Patch(facecolor='#C62828', label='北日本地域（札幌・仙台・新潟）'),
                    Patch(facecolor='#4C72B0', label='その他の地域')],
           loc='lower left', fontsize=10, framealpha=0.9)
ax2.grid(True, axis='x', alpha=0.3)
ax2.text(0.99, 0.02, NOTE, transform=ax2.transAxes, ha='right', va='bottom',
         fontsize=9, color='#888',
         bbox=dict(boxstyle='round', facecolor='#F5F5F5', alpha=0.8))
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_H2_fig2.png'), dpi=150, bbox_inches='tight')
plt.close()
print('図2（相関係数）保存完了')

# 図3: 偏相関分析の結果（原論文 4.3節 仮説3-1/3-2 の報告値）
# ============================================================
fig3, ax3 = plt.subplots(figsize=(10, 6.5))
x = np.arange(len(pcorr_labels))
bars = ax3.bar(x, pcorr_values, color=pcorr_colors, edgecolor='white', width=0.6)
for xi, v in zip(x, pcorr_values):
    ax3.text(xi, v - 0.02, f'{v:.2f}', ha='center', va='top',
             fontsize=13, color='white', fontweight='bold')
ax3.set_xticks(x)
ax3.set_xticklabels(pcorr_labels, fontsize=10.5)
ax3.axhline(0, color='#333', lw=1)
ax3.set_ylim(-0.55, 0.08)
ax3.set_ylabel('（偏）相関係数', fontsize=12)
ax3.set_title('図3：床暖房普及率を第3変数とした偏相関分析（原論文 4.3節）\n'
              '床暖房の影響を除くと相関は −0.33 → −0.40 とむしろ強まった（仮説に反する結果）',
              fontsize=12.5, fontweight='bold')
# −0.33 → −0.40 の強化を示す矢印
ax3.annotate('', xy=(0.72, -0.40), xytext=(0.28, -0.33),
             arrowprops=dict(arrowstyle='->', color='#C62828', lw=2.2,
                             connectionstyle='arc3,rad=-0.35'))
ax3.text(2.5, -0.50, '室内要因（床暖房）を除いても、外気温と脳卒中死亡率の相関は残る',
         ha='center', va='center', fontsize=11, color='#C62828', fontweight='bold')
ax3.grid(True, axis='y', alpha=0.3)
ax3.text(0.99, 0.02, NOTE, transform=ax3.transAxes, ha='right', va='bottom',
         fontsize=9, color='#888',
         bbox=dict(boxstyle='round', facecolor='#F5F5F5', alpha=0.8))
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_H2_fig3.png'), dpi=150, bbox_inches='tight')
plt.close()
print('図3（偏相関分析）保存完了')

# 図4: 6つの仮説と検証結果のまとめ（原論文5節の解釈より）
# ============================================================
fig4, ax4 = plt.subplots(figsize=(11, 6.5))
ax4.axis('off')
verdict_color = {'支持': '#2E7D32', '支持（弱）': '#F9A825', '不支持': '#C62828'}
n = len(hypotheses)
row_h = 1.0
for i, (hid, desc, verdict, note) in enumerate(hypotheses):
    yy = (n - i - 1) * row_h
    ax4.add_patch(plt.Rectangle((0, yy + 0.06), 11, row_h - 0.12,
                                facecolor='#F8F9FA' if i % 2 else '#FFFFFF',
                                edgecolor='#DEE2E6', lw=0.8))
    ax4.text(0.15, yy + row_h/2, hid, fontsize=12, fontweight='bold',
             va='center', ha='left', color='#333')
    ax4.text(1.35, yy + row_h/2, desc, fontsize=11, va='center', ha='left')
    ax4.add_patch(plt.Rectangle((7.05, yy + 0.22), 1.55, row_h - 0.44,
                                facecolor=verdict_color[verdict], alpha=0.9))
    ax4.text(7.82, yy + row_h/2, verdict, fontsize=11, fontweight='bold',
             va='center', ha='center', color='white')
    ax4.text(8.85, yy + row_h/2, note, fontsize=8.5, va='center', ha='left', color='#555')
ax4.set_xlim(0, 11)
ax4.set_ylim(-0.2, n * row_h + 0.6)
ax4.text(0, n * row_h + 0.15, '図4：6つの仮説と検証結果のまとめ',
         fontsize=14, fontweight='bold', ha='left', color='#C62828')
ax4.text(11, -0.15, NOTE, fontsize=9, color='#888', ha='right', va='top')
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_H2_fig4.png'), dpi=150, bbox_inches='tight')
plt.close()
print('図4（仮説まとめ）保存完了')

# ============================================================
# コンソール出力（HTML埋め込み確認用）
# ============================================================
print('\n=== 原論文の報告値（転記）===')
print('■ 2018年 死因割合: 悪性新生物27% / 心疾患15% / 脳卒中8% / その他49%')
print('■ 表4 相関係数（月間平均気温×脳卒中死亡率, 2019年）')
for c, v, nth in zip(corr_cities, corr_values, corr_is_north):
    print(f'   {c:<8} r = {v:+.2f}  {"[北日本]" if nth else ""}')
print('■ 仮説3の（偏）相関')
for l, v in zip(['気温×脳卒中(除く前)', '気温×脳卒中(除いた後・偏相関)',
                 '気温×床暖房普及率', '床暖房×脳卒中(気温除く偏相関・1月)'], pcorr_values):
    print(f'   {l:<28} = {v:+.2f}')
print('\n全図の生成が完了しました。')
print(f'出力先: {os.path.abspath(FIG_DIR)}')