周辺環境からみた新型コロナウイルス感染症の要因分析

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()

# 数値列を変換
num_cols = [
    '総人口', '65歳以上人口', '延べ宿泊者数', '外国人延べ宿泊者数',
    '旅館営業施設客室数（ホテルを含む）', '旅館営業施設数（ホテルを含む）',
    '年平均気温', '交通・通信費（二人以上の世帯）', '消費支出（二人以上の世帯）',
    '月間有効求人数（一般）', '標準価格（平均価格）（住宅地）',
]
for c in num_cols:
    if c in df_b.columns:
        df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

print("=" * 60)
print(f"■ SSDSE-B 都道府県データ: {df_b.shape[0]}行 × {df_b.shape[1]}列")
print(f"  年度範囲: {df_b['年度'].min()} 〜 {df_b['年度'].max()}")
print(f"  都道府県数: {df_b['都道府県'].nunique()}")
print("=" * 60)

df_b['宿泊密度'] = df_b['延べ宿泊者数'] / df_b['総人口'].clip(1)
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'].clip(1) * 100  # %
df_b['外国人宿泊率'] = df_b['外国人延べ宿泊者数'] / df_b['延べ宿泊者数'].clip(1) * 100  # %

df_2019 = df_b[df_b['年度'] == 2019].copy()
df_2020 = df_b[df_b['年度'] == 2020].copy()

print(f"\n■ 2019年断面: N={len(df_2019)}")
print(f"  宿泊密度 平均={df_2019['宿泊密度'].mean():.2f}, 中央値={df_2019['宿泊密度'].median():.2f}")

ENV_VARS = {
    '高齢化率(%)':      '高齢化率',
    '宿泊施設客室数':   '旅館営業施設客室数（ホテルを含む）',
    '交通通信費(万円)': '交通・通信費（二人以上の世帯）',
    '消費支出(万円)':   '消費支出（二人以上の世帯）',
    '有効求人数(千)':   '月間有効求人数（一般）',
    '地価住宅地':       '標準価格（平均価格）（住宅地）',
}

df_2019 = df_2019.copy()
df_2019['高齢化率'] = df_2019['高齢化率']
df_2019['宿泊施設客室数']   = pd.to_numeric(df_2019['旅館営業施設客室数（ホテルを含む）'], errors='coerce')
df_2019['交通通信費(万円)'] = pd.to_numeric(df_2019['交通・通信費（二人以上の世帯）'], errors='coerce') / 10000
df_2019['消費支出(万円)']   = pd.to_numeric(df_2019['消費支出（二人以上の世帯）'], errors='coerce') / 10000
df_2019['有効求人数(千)']   = pd.to_numeric(df_2019['月間有効求人数（一般）'], errors='coerce') / 1000
df_2019['地価住宅地']       = pd.to_numeric(df_2019['標準価格（平均価格）（住宅地）'], errors='coerce')

XVARS = ['高齢化率', '宿泊施設客室数', '交通通信費(万円)', '消費支出(万円)', '有効求人数(千)', '地価住宅地']
XVAR_LABELS = {
    '高齢化率':         '高齢化率（%）',
    '宿泊施設客室数':   '宿泊施設客室数',
    '交通通信費(万円)': '交通・通信費（万円）',
    '消費支出(万円)':   '消費支出（万円）',
    '有効求人数(千)':   '月間有効求人数（千件）',
    '地価住宅地':       '住宅地地価（千円/m²）',
}

df_2019_clean = df_2019[['都道府県', '宿泊密度'] + XVARS].dropna()
N = len(df_2019_clean)
print(f"\n■ 2019年 分析用データ: N={N}（欠損除外後）")

df_merge = pd.merge(
    df_2019[['都道府県', '宿泊密度']].rename(columns={'宿泊密度': '宿泊密度_2019'}),
    df_2020[['都道府県', '宿泊密度']].rename(columns={'宿泊密度': '宿泊密度_2020'}),
    on='都道府県'
)
df_merge['変化率'] = (df_merge['宿泊密度_2020'] - df_merge['宿泊密度_2019']) / df_merge['宿泊密度_2019'].clip(0.001) * 100
df_merge = df_merge.sort_values('変化率').reset_index(drop=True)

print("\n■ 2019→2020 宿泊密度変化率 上位・下位5件:")
print(df_merge[['都道府県', '変化率']].tail(5).to_string())
print("...")
print(df_merge[['都道府県', '変化率']].head(5).to_string())

fig1, ax1 = plt.subplots(figsize=(10, 5))

years  = ts['年度'].values
values = ts['宿泊密度_平均'].values

# COVID期（2020〜2021）シェード
ax1.axvspan(2019.5, 2021.5, alpha=0.18, color='#C62828', label='COVID-19 流行期（2020〜2021）')

ax1.plot(years, values, 'o-', color='#1565C0', linewidth=2.2, markersize=6,
         label='延べ宿泊者数（人口一人あたり）全国平均', zorder=3)

# 2019→2020 矢印注釈
y19 = ts[ts['年度'] == 2019]['宿泊密度_平均'].values
y20 = ts[ts['年度'] == 2020]['宿泊密度_平均'].values
if len(y19) and len(y20):
    ax1.annotate('', xy=(2020, y20[0]), xytext=(2019, y19[0]),
                 arrowprops=dict(arrowstyle='->', color='#C62828', lw=2.0))
    ax1.text(2019.5, (y19[0] + y20[0]) / 2 + 0.05, '▼COVID衝撃', color='#C62828',
             fontsize=9, ha='center')

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('延べ宿泊者数（人口一人あたり：泊/人）', fontsize=12)
ax1.set_title('延べ宿泊者数（人口あたり）の推移 2012〜2023\n（都道府県平均・SSDSE-B実データ）',
              fontsize=13, fontweight='bold')
ax1.set_xticks(years)
ax1.set_xticklabels([str(y) for y in years], rotation=45, ha='right')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_U5_6_fig1_ts.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2022_U5_6_fig1_ts.png")

fig2, ax2 = plt.subplots(figsize=(8, 7))

# 列ラベルを短くする
label_map = {
    '宿泊密度':         '宿泊密度\n(一人あたり)',
    '高齢化率':         '高齢化率\n(%)',
    '宿泊施設客室数':   '宿泊施設\n客室数',
    '交通通信費(万円)': '交通・\n通信費',
    '消費支出(万円)':   '消費支出',
    '有効求人数(千)':   '有効求人数',
    '地価住宅地':       '住宅地地価',
}
cm_data = df_2019_clean[corr_cols].rename(columns=label_map)
cm = cm_data.corr()

im = ax2.imshow(cm.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax2, fraction=0.046, pad=0.04, label='相関係数')

labels = list(cm.columns)
ax2.set_xticks(range(len(labels)))
ax2.set_yticks(range(len(labels)))
ax2.set_xticklabels(labels, fontsize=9)
ax2.set_yticklabels(labels, fontsize=9)

for i in range(len(labels)):
    for j in range(len(labels)):
        v = cm.values[i, j]
        tc = 'white' if abs(v) > 0.6 else 'black'
        ax2.text(j, i, f'{v:.2f}', ha='center', va='center', fontsize=9, color=tc)

ax2.set_title('相関ヒートマップ：宿泊密度と環境変数\n（2019年断面, 都道府県 N=47）',
              fontsize=12, fontweight='bold')
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2022_U5_6_fig2_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2022_U5_6_fig2_corr.png")

fig3, ax3 = plt.subplots(figsize=(8, 5))

betas = model_std.params[1:]
ses   = model_std.bse[1:]
pvals = model_std.pvalues[1:]

xlabels = [XVAR_LABELS.get(v, v) for v in XVARS]
colors3 = ['#C62828' if p < 0.05 else '#90A4AE' for p in pvals]
y_pos   = np.arange(len(XVARS))

ax3.barh(y_pos, betas, color=colors3, alpha=0.78, edgecolor='white', height=0.55)
ax3.errorbar(betas, y_pos, xerr=1.96 * ses, fmt='none', color='#333',
             capsize=4, linewidth=1.3, zorder=5)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos)
ax3.set_yticklabels(xlabels, fontsize=10)
ax3.set_xlabel('標準化回帰係数 Beta（±1.96SE）', fontsize=11)
ax3.set_title(f'OLS標準化係数（2019年断面, N={N}都道府県）\n'
              f'R²={model_std.rsquared:.3f}  目的変数：延べ宿泊者数（人口あたり）',
              fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

sig_patch   = mpatches.Patch(color='#C62828', alpha=0.78, label='p < 0.05')
insig_patch = mpatches.Patch(color='#90A4AE', alpha=0.78, label='n.s.')
ax3.legend(handles=[sig_patch, insig_patch], fontsize=9, loc='lower right')

for i, (b, p) in enumerate(zip(betas, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    ax3.text(b + (0.01 if b >= 0 else -0.01), i,
             f' {b:.3f}{sig}',
             va='center', ha='left' if b >= 0 else 'right', fontsize=8.5)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2022_U5_6_fig3_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2022_U5_6_fig3_coef.png")

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import statsmodels.api as sm
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

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

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

ts = df_b.groupby('年度')['宿泊密度'].mean().reset_index()
ts.columns = ['年度', '宿泊密度_平均']

print("\n■ 宿泊密度の年次推移（全国平均）:")
print(ts.to_string(index=False))

corr_cols = ['宿泊密度'] + XVARS
corr_mat = df_2019_clean[corr_cols].corr()
print("\n■ 相関行列（2019年断面）:")
print(corr_mat.round(3))

y_reg = df_2019_clean['宿泊密度'].values
X_reg = df_2019_clean[XVARS].values

# 標準化（Betaコefのため）
y_std = (y_reg - y_reg.mean()) / y_reg.std()
X_std = (X_reg - X_reg.mean(axis=0)) / X_reg.std(axis=0)

model_std = sm.OLS(y_std, sm.add_constant(X_std)).fit()
model_raw = sm.OLS(y_reg, sm.add_constant(X_reg)).fit()

print("\n■ OLS重回帰（標準化係数、2019年断面）:")
print(f"  R² = {model_std.rsquared:.4f}")
print(f"  {'変数':<20} {'Beta':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 50)
for i, var in enumerate(XVARS):
    beta = model_std.params[i + 1]
    p    = model_std.pvalues[i + 1]
    sig  = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {var:<20} {beta:>8.4f} {p:>10.4f} {sig:>6}")

fig4, ax4 = plt.subplots(figsize=(9, 13))

colors4 = ['#C62828' if v < 0 else '#1565C0' for v in df_merge['変化率']]
y_pos4  = np.arange(len(df_merge))

ax4.barh(y_pos4, df_merge['変化率'], color=colors4, alpha=0.75,
         edgecolor='white', height=0.75)
ax4.axvline(0, color='black', linewidth=0.8)
ax4.set_yticks(y_pos4)
ax4.set_yticklabels(df_merge['都道府県'], fontsize=9)
ax4.set_xlabel('変化率（%）', fontsize=11)
ax4.set_title('延べ宿泊者数（人口あたり）の変化率：2019年→2020年\n（都道府県別, SSDSE-B実データ）',
              fontsize=12, fontweight='bold')
ax4.grid(axis='x', alpha=0.3)

dec_patch = mpatches.Patch(color='#C62828', alpha=0.75, label='減少（赤）')
inc_patch = mpatches.Patch(color='#1565C0', alpha=0.75, label='増加（青）')
ax4.legend(handles=[dec_patch, inc_patch], fontsize=9, loc='lower right')

# 国全体平均線
nat_mean = df_merge['変化率'].mean()
ax4.axvline(nat_mean, color='#FF8F00', linestyle='--', linewidth=1.5,
            label=f'全国平均 {nat_mean:.1f}%')
ax4.legend(handles=[dec_patch, inc_patch,
                    mpatches.Patch(color='#FF8F00', label=f'全国平均 {nat_mean:.1f}%')],
           fontsize=9, loc='lower right')

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2022_U5_6_fig4_change.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2022_U5_6_fig4_change.png")

print("\n全図の生成完了（4枚）")
print("  2022_U5_6_fig1_ts.png     : 宿泊密度の時系列（COVID期シェード）")
print("  2022_U5_6_fig2_corr.png   : 相関ヒートマップ")
print("  2022_U5_6_fig3_coef.png   : OLS標準化係数プロット（2019年断面）")
print("  2022_U5_6_fig4_change.png : 都道府県別変化率 2019→2020（横棒）")