顎・足・枕が多様な観光客にどのような影響を及ぼすのか

df_b = pd.read_csv(DATA_PATH, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)]
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)

print("=" * 60)
print(f"■ SSDSE-B-2026 読み込み完了: {len(df)}都道府県（2022年度）")
print("=" * 60)

num_cols = [
    '総人口', '65歳以上人口',
    '延べ宿泊者数', '外国人延べ宿泊者数',
    '旅館営業施設数（ホテルを含む）', '旅館営業施設客室数（ホテルを含む）',
    '食料費（二人以上の世帯）', '消費支出（二人以上の世帯）',
    '交通・通信費（二人以上の世帯）',
    '年平均気温',
]
for c in num_cols:
    df[c] = pd.to_numeric(df[c], errors='coerce')

# 目的変数
df['宿泊者数_per_capita'] = df['延べ宿泊者数'] / df['総人口']          # 一人当たり宿泊需要
df['外国人比率']          = df['外国人延べ宿泊者数'] / df['延べ宿泊者数']  # 外国人シェア

# 説明変数
df['食料費割合']    = df['食料費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）']   # 顎（食）
df['交通費割合']    = df['交通・通信費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）']  # 足（交通）
df['旅館数密度']    = df['旅館営業施設数（ホテルを含む）'] / df['総人口'] * 10000          # 枕（宿泊インフラ）
df['客室数_per_施設'] = df['旅館営業施設客室数（ホテルを含む）'] / df['旅館営業施設数（ホテルを含む）']  # 施設規模
df['気温']          = df['年平均気温']
df['高齢化率']      = df['65歳以上人口'] / df['総人口']

# 分析列を揃えてNaN行を削除
FEATURES = ['食料費割合', '交通費割合', '旅館数密度', '客室数_per_施設', '気温', '高齢化率']
OUTCOMES = ['宿泊者数_per_capita', '外国人比率']
df_model = df[['都道府県'] + FEATURES + OUTCOMES].dropna().reset_index(drop=True)

print(f"■ 欠損除去後: {len(df_model)}都道府県\n")
print(df_model[FEATURES + OUTCOMES].describe().round(4))

ALL_VARS = FEATURES + OUTCOMES
ALL_LABELS = FEATURE_LABELS + ['宿泊者数\nper capita', '外国人比率']
corr_mat = df_model[ALL_VARS].corr().values
n_vars = len(ALL_VARS)

fig1, ax1 = plt.subplots(figsize=(10, 8))
im = ax1.imshow(corr_mat, cmap='RdBu_r', vmin=-1, vmax=1)
ax1.set_xticks(range(n_vars))
ax1.set_yticks(range(n_vars))
ax1.set_xticklabels(ALL_LABELS, rotation=30, ha='right', fontsize=9)
ax1.set_yticklabels(ALL_LABELS, fontsize=9)
ax1.set_title('顎（食料費）・足（交通費）・枕（旅館密度）と観光需要の相関行列\n（SSDSE-B 2022年度 都道府県別）',
              fontsize=11, fontweight='bold')
for i in range(n_vars):
    for j in range(n_vars):
        r = corr_mat[i, j]
        txt_col = 'white' if abs(r) > 0.55 else 'black'
        ax1.text(j, i, f'{r:.2f}', ha='center', va='center',
                 fontsize=9, color=txt_col, fontweight='bold')
plt.colorbar(im, ax=ax1, shrink=0.8, label='相関係数')
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2024_H5_6_fig1_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2024_H5_6_fig1_corr.png")

x2 = df_model['旅館数密度'].values
y2 = df_model['宿泊者数_per_capita'].values
prefs = df_model['都道府県'].tolist()

r2_val, p2_val = stats.pearsonr(x2, y2)

fig2, ax2 = plt.subplots(figsize=(10, 7))
ax2.scatter(x2, y2, color='#E65100', s=50, alpha=0.7, zorder=3)

# 都道府県ラベル（全都道府県）
HIGHLIGHT = {'山梨県': '#C62828', '長野県': '#C62828', '静岡県': '#C62828',
             '北海道': '#1565C0', '沖縄県': '#1565C0', '鳥取県': '#2E7D32'}
for i, pref in enumerate(prefs):
    color = HIGHLIGHT.get(pref, '#555555')
    fontsize = 8.5 if pref in HIGHLIGHT else 7
    fw = 'bold' if pref in HIGHLIGHT else 'normal'
    ax2.annotate(pref, (x2[i], y2[i]),
                 textcoords='offset points', xytext=(4, 3),
                 fontsize=fontsize, color=color, fontweight=fw)

# ハイライト点
for pref, col in HIGHLIGHT.items():
    idx = prefs.index(pref) if pref in prefs else None
    if idx is not None:
        ax2.scatter(x2[idx], y2[idx], color=col, s=100, zorder=5)

# 回帰直線
z2 = np.polyfit(x2, y2, 1)
xs2 = np.linspace(x2.min() - 0.1, x2.max() + 0.1, 200)
ax2.plot(xs2, np.poly1d(z2)(xs2), 'r-', linewidth=2, alpha=0.8,
         label=f'回帰直線  r={r2_val:.2f}, p={p2_val:.4f}')

ax2.set_xlabel('旅館数密度（10,000人当たり旅館数）【枕の代理指標】', fontsize=11)
ax2.set_ylabel('一人当たり延べ宿泊者数（泊/人）', fontsize=11)
ax2.set_title('旅館インフラ密度（枕）と観光宿泊需要の関係\n（SSDSE-B 2022年度 都道府県別）',
              fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2024_H5_6_fig2_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2024_H5_6_fig2_scatter.png")

fig3, ax3 = plt.subplots(figsize=(10, 5.5))

x_pos3 = np.arange(len(FEATURES))
width = 0.38
bars1 = ax3.bar(x_pos3 - width/2, beta1, width,
                color='#E65100', alpha=0.78, edgecolor='white', label='一人当たり延べ宿泊者数')
bars2 = ax3.bar(x_pos3 + width/2, beta2, width,
                color='#1565C0', alpha=0.78, edgecolor='white', label='外国人延べ宿泊者数比率')

for bars, betas in [(bars1, beta1), (bars2, beta2)]:
    for bar, b in zip(bars, betas):
        va = 'bottom' if b >= 0 else 'top'
        offset = 0.01 if b >= 0 else -0.01
        ax3.text(bar.get_x() + bar.get_width() / 2,
                 bar.get_height() + offset,
                 f'β={b:.2f}', ha='center', va=va, fontsize=8.5, fontweight='bold')

ax3.axhline(0, color='gray', linestyle='--', linewidth=0.8)
ax3.set_xticks(x_pos3)
ax3.set_xticklabels(FEATURE_LABELS, fontsize=10)
ax3.set_ylabel('標準化回帰係数 β', fontsize=11)
ax3.set_title('顎・足・枕の標準化回帰係数\n一人当たり宿泊者数モデル vs 外国人比率モデル',
              fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(axis='y', alpha=0.3)

ymax = max(abs(beta1).max(), abs(beta2).max())
ax3.set_ylim(-ymax * 1.35, ymax * 1.45)

# R²テキスト
ax3.text(0.02, 0.97,
         f'R²（宿泊者数モデル）= {model1.rsquared:.3f}\nR²（外国人比率モデル）= {model2.rsquared:.3f}',
         transform=ax3.transAxes, fontsize=9, va='top',
         bbox=dict(facecolor='#F5F5F5', edgecolor='#BDBDBD', boxstyle='round,pad=0.4'))

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

import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import statsmodels.api as sm
from sklearn.preprocessing import StandardScaler
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'
os.makedirs(FIG_DIR, exist_ok=True)

DATA_PATH = 'data/raw/SSDSE-B-2026.csv'

def run_ols(y_col, df_m, feature_cols, label):
    X = sm.add_constant(df_m[feature_cols].values)
    y = df_m[y_col].values
    model = sm.OLS(y, X).fit()
    print(f"\n{'='*60}")
    print(f"■ OLS: 目的変数 = {label} ({y_col})")
    print(model.summary(xname=['const'] + feature_cols))
    return model

model1 = run_ols('宿泊者数_per_capita', df_model, FEATURES, '一人当たり延べ宿泊者数')
model2 = run_ols('外国人比率',          df_model, FEATURES, '外国人延べ宿泊者数比率')

# 標準化係数
scaler = StandardScaler()
X_std = scaler.fit_transform(df_model[FEATURES].values)

def std_coefs(y_col, X_s, df_m):
    y = df_m[y_col].values
    y_s = (y - y.mean()) / y.std()
    m = sm.OLS(y_s, sm.add_constant(X_s)).fit()
    return m.params[1:]

beta1 = std_coefs('宿泊者数_per_capita', X_std, df_model)
beta2 = std_coefs('外国人比率',          X_std, df_model)

FEATURE_LABELS = ['食料費割合\n（顎）', '交通費割合\n（足）', '旅館数密度\n（枕）',
                  '客室数/施設', '気温', '高齢化率']

print("\n■ 標準化回帰係数（一人当たり宿泊者数モデル）:", beta1.round(4))
print("■ 標準化回帰係数（外国人比率モデル）:",        beta2.round(4))

top10 = df_model.nlargest(10, '宿泊者数_per_capita').reset_index(drop=True)
top10_prefs  = top10['都道府県'].tolist()
top10_vals   = top10['宿泊者数_per_capita'].values
top10_fratio = top10['外国人比率'].values

fig4, ax4a = plt.subplots(figsize=(10, 6))
ax4b = ax4a.twinx()

colors4 = ['#C62828' if p == '鳥取県' else '#E65100' for p in top10_prefs]
bars4 = ax4a.bar(range(10), top10_vals, color=colors4, alpha=0.78, edgecolor='white',
                 label='一人当たり延べ宿泊者数（左軸）')
ax4b.plot(range(10), top10_fratio * 100, 'o-',
          color='#1565C0', linewidth=2, markersize=7, label='外国人比率 %（右軸）')

for i, (v, f) in enumerate(zip(top10_vals, top10_fratio)):
    ax4a.text(i, v + 0.02, f'{v:.2f}', ha='center', va='bottom', fontsize=8.5, fontweight='bold')

ax4a.set_xticks(range(10))
ax4a.set_xticklabels(top10_prefs, rotation=25, ha='right', fontsize=10)
ax4a.set_ylabel('一人当たり延べ宿泊者数（泊/人）', fontsize=11)
ax4b.set_ylabel('外国人延べ宿泊者数比率（%）', fontsize=11, color='#1565C0')
ax4b.tick_params(axis='y', colors='#1565C0')
ax4a.set_title('一人当たり延べ宿泊者数 上位10都道府県\n（SSDSE-B 2022年度／外国人比率も併記）',
               fontsize=12, fontweight='bold')

# 凡例を統合
lines1, labs1 = ax4a.get_legend_handles_labels()
lines2, labs2 = ax4b.get_legend_handles_labels()
ax4a.legend(lines1 + lines2, labs1 + labs2, fontsize=9, loc='upper right')
ax4a.grid(axis='y', alpha=0.3)

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

print("\n全図の生成完了（4枚）")