地価に関する最適モデルの構築と手法提案

import os
import warnings
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 sklearn.linear_model import Ridge, Lasso, LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler

warnings.filterwarnings('ignore')

# ── パス設定 ──────────────────────────────────────────────────────────
DATA_DIR = 'data/raw'
FIG_DIR  = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family': 'Hiragino Sans',
    'axes.unicode_minus': False,
    'figure.dpi': 150,
    'axes.spines.top': False,
    'axes.spines.right': False,
})

# ── SSDSE-B-2026.csv 読み込み ──────────────────────────────────────
print("=== SSDSE-B 読み込み ===")
df_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932',
    header=0,
    skiprows=[1],   # 日本語ラベル行をスキップ（先頭行=コード名を使用）
)

# 都道府県レベルのみ（地域コード = R + 5桁数字）
df_raw = df_raw[df_raw['Code'].str.match(r'^R\d{5}$', na=False)].copy()

# 年度列の名前統一（先頭列が年度）
year_col = df_raw.columns[0]   # 'SSDSE-B-2026'
df_raw[year_col] = pd.to_numeric(df_raw[year_col], errors='coerce')

# 2022年度のみ使用
df = df_raw[df_raw[year_col] == 2022].copy()
print(f"2022年度データ: {len(df)}都道府県")

# ── 数値変換 ────────────────────────────────────────────────────────
NUM_COLS = [
    'C5401',    # 標準価格（住宅地, 千円/m²）
    'C5403',    # 標準価格（商業地, 千円/m²）
    'A1101',    # 総人口
    'A1302',    # 15〜64歳人口（生産年齢人口）
    'F3102',    # 有効求人数（新規）
    'F3103',    # 有効求職者数（新規）
    'G7101',    # 宿泊者数（延べ）
    'E6302',    # 大学学生数
    'I510120',  # 病院数
    'B4101',    # 年平均気温
    'A4103',    # 合計特殊出生率
    'A5101',    # 転入者数
]
for c in NUM_COLS:
    df[c] = pd.to_numeric(df[c], errors='coerce')

# ── 派生変数の計算（実データのみ） ─────────────────────────────────
# 生産年齢人口割合 (%)
df['working_age_pct'] = df['A1302'] / df['A1101'] * 100

# 有効求人倍率 = 有効求人数 / 有効求職者数
df['job_ratio'] = df['F3103'] / df['F3102']

# 宿泊者数 per capita
df['tourism_pc'] = df['G7101'] / df['A1101']

# 大学学生率 (%)
df['univ_rate'] = df['E6302'] / df['A1101'] * 100

# 病院密度（人口1万対）
df['hospital_density'] = df['I510120'] / df['A1101'] * 10000

# 転入率 (%)
df['inflow_rate'] = df['A5101'] / df['A1101'] * 100

# 目的変数・説明変数の列名
TARGET  = 'C5401'   # 住宅地標準価格（千円/m²）
PRED_COLS = [
    'A1101',           # 総人口（万人換算しない）
    'working_age_pct', # 生産年齢人口割合
    'job_ratio',       # 有効求人倍率
    'tourism_pc',      # 宿泊者数per capita
    'univ_rate',       # 大学学生率
    'hospital_density',# 病院密度
    'B4101',           # 年平均気温
    'A4103',           # 合計特殊出生率
    'inflow_rate',     # 転入率
]
PRED_LABELS = [
    '総人口',
    '生産年齢人口割合',
    '有効求人倍率',
    '宿泊者数per capita',
    '大学学生率',
    '病院密度',
    '年平均気温',
    '合計特殊出生率',
    '転入率',
]

# ── 分析用データ整理 ────────────────────────────────────────────────
df_model = df[['Prefecture', TARGET] + PRED_COLS].dropna().copy()
df_model = df_model.reset_index(drop=True)
print(f"欠損除去後サンプル数: {len(df_model)}")

# 東京都フラグ
tokyo_mask = df_model['Prefecture'].str.contains('東京', na=False)
df_no_tokyo = df_model[~tokyo_mask].copy().reset_index(drop=True)
print(f"東京除外後サンプル数: {len(df_no_tokyo)}")

# ── OLS 重回帰（statsmodels, 東京除外）────────────────────────────
print("\n=== OLS 重回帰（statsmodels, 東京除外）===")
X_ols = sm.add_constant(df_no_tokyo[PRED_COLS].astype(float))
y_ols = df_no_tokyo[TARGET].astype(float)
res_ols = sm.OLS(y_ols, X_ols).fit()
print(res_ols.summary())
print(f"\nAIC = {res_ols.aic:.2f}")
print(f"BIC = {res_ols.bic:.2f}")
print(f"R²  = {res_ols.rsquared:.4f}")
print(f"Adj.R² = {res_ols.rsquared_adj:.4f}")

# ── 標準化係数の計算 ───────────────────────────────────────────────
scaler = StandardScaler()
X_scaled_nt = scaler.fit_transform(df_no_tokyo[PRED_COLS].astype(float))
y_nt        = df_no_tokyo[TARGET].astype(float).values

# 標準化 OLS
lr = LinearRegression().fit(X_scaled_nt, y_nt)
std_coefs = lr.coef_

# ── sklearn モデル比較（5分割 CV R²）── 東京除外 ──────────────────
print("\n=== モデル比較（5分割 CV R²）===")
cv = 5

# OLS（LinearRegression）
cv_ols   = cross_val_score(LinearRegression(), X_scaled_nt, y_nt, cv=cv, scoring='r2')
# Ridge（α=1.0）
cv_ridge = cross_val_score(Ridge(alpha=1.0), X_scaled_nt, y_nt, cv=cv, scoring='r2')
# Lasso（α=0.1）
cv_lasso = cross_val_score(Lasso(alpha=0.1, max_iter=10000), X_scaled_nt, y_nt, cv=cv, scoring='r2')

print(f"OLS   CV-R²: {cv_ols.mean():.4f}  ± {cv_ols.std():.4f}")
print(f"Ridge CV-R²: {cv_ridge.mean():.4f}  ± {cv_ridge.std():.4f}")
print(f"Lasso CV-R²: {cv_lasso.mean():.4f}  ± {cv_lasso.std():.4f}")

# Lasso 係数（変数選択確認）
lasso_model = Lasso(alpha=0.1, max_iter=10000).fit(X_scaled_nt, y_nt)
print("\nLasso 係数（ゼロ = 変数選択で除外）:")
for label, coef in zip(PRED_LABELS, lasso_model.coef_):
    mark = "  ← 除外" if abs(coef) < 1e-6 else ""
    print(f"  {label}: {coef:.4f}{mark}")

# ── Figure 1: 人口 vs 住宅地価格（東京外れ値）────────────────────
print("\nFigure 1: 人口 vs 住宅地価格（対数軸）...")
fig, ax = plt.subplots(figsize=(9, 6))

ax.scatter(df_no_tokyo['A1101'] / 1e6,
           df_no_tokyo[TARGET],
           color='#1565C0', alpha=0.75, s=60, zorder=3, label='各都道府県（東京除外）')

# 東京を別途プロット
if tokyo_mask.any():
    df_tokyo = df_model[tokyo_mask]
    ax.scatter(df_tokyo['A1101'] / 1e6,
               df_tokyo[TARGET],
               color='#C62828', s=120, zorder=4, marker='*',
               label='東京都（外れ値）')
    for _, row in df_tokyo.iterrows():
        ax.annotate('東京都',
                    xy=(row['A1101'] / 1e6, row[TARGET]),
                    xytext=(8, -14), textcoords='offset points',
                    fontsize=10, color='#C62828', fontweight='bold')

# 主要都市をラベル
for _, row in df_no_tokyo.iterrows():
    if row['A1101'] > 5e6 or row[TARGET] > 90000:
        ax.annotate(row['Prefecture'],
                    xy=(row['A1101'] / 1e6, row[TARGET]),
                    xytext=(4, 4), textcoords='offset points',
                    fontsize=8, color='#333333')

ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('総人口（百万人, 対数軸）', fontsize=12)
ax.set_ylabel('住宅地標準価格（千円/m², 対数軸）', fontsize=12)
ax.set_title('総人口と住宅地標準価格の関係\n（2022年度、47都道府県）',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2023_H3_fig1_scatter.png'), bbox_inches='tight')
plt.close()

print("Figure 1 saved.")

# ── Figure 2: OLS 標準化係数棒グラフ ─────────────────────────────
print("Figure 2: OLS 標準化係数棒グラフ...")

# |β| でソート
sort_idx = np.argsort(np.abs(std_coefs))[::-1]
sorted_labels = [PRED_LABELS[i] for i in sort_idx]
sorted_coefs  = std_coefs[sort_idx]

# p値（元のOLSから対応する順番）
pvals_ols = res_ols.pvalues[1:].values  # 定数項を除く
sorted_pvals = pvals_ols[sort_idx]

colors_bar = ['#C62828' if c < 0 else '#1565C0' for c in sorted_coefs]
sig_marks  = ['***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))
              for p in sorted_pvals]

fig, ax = plt.subplots(figsize=(10, 6))
y_pos = np.arange(len(sorted_labels))
bars  = ax.barh(y_pos, sorted_coefs, color=colors_bar, alpha=0.85, height=0.65)

for bar, mark in zip(bars, sig_marks):
    x = bar.get_width()
    offset = 0.012 if x >= 0 else -0.012
    ax.text(x + offset, bar.get_y() + bar.get_height() / 2,
            mark, va='center', ha='left' if x >= 0 else 'right',
            fontsize=11, color='#333333')

ax.axvline(0, color='black', lw=1)
ax.set_yticks(y_pos)
ax.set_yticklabels([f'{l}' for l in sorted_labels], fontsize=11)
ax.set_xlabel('標準化回帰係数（β）', fontsize=12)
ax.set_title('住宅地価格モデル：OLS 標準化係数\n（東京除外、|β|降順、*** p<0.001, ** p<0.01, * p<0.05）',
             fontsize=12, fontweight='bold')

neg_p = mpatches.Patch(color='#C62828', alpha=0.85, label='負の効果')
pos_p = mpatches.Patch(color='#1565C0', alpha=0.85, label='正の効果')
ax.legend(handles=[pos_p, neg_p], fontsize=10, loc='lower right')
ax.grid(axis='x', alpha=0.3)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2023_H3_fig2_coef.png'), bbox_inches='tight')
plt.close()

print("Figure 2 saved.")

# ── Figure 3: OLS vs Ridge vs Lasso CV-R² 比較 ───────────────────
print("Figure 3: モデル比較（CV-R²）...")
model_names  = ['OLS\n(最小二乗法)', 'Ridge\n(α=1.0)', 'Lasso\n(α=0.1)']
cv_means     = [cv_ols.mean(),   cv_ridge.mean(),   cv_lasso.mean()]
cv_stds      = [cv_ols.std(),    cv_ridge.std(),     cv_lasso.std()]
bar_colors_m = ['#1565C0', '#2E7D32', '#E65100']

fig, ax = plt.subplots(figsize=(9, 6))
x_pos = np.arange(len(model_names))

# プロット：バーは 0 基準、正・負どちらも表示
bars_m = ax.bar(x_pos, cv_means, yerr=cv_stds, capsize=8,
                color=bar_colors_m, alpha=0.85, width=0.5,
                error_kw={'elinewidth': 2, 'ecolor': '#333333'})

for i, (mean_val, std_val) in enumerate(zip(cv_means, cv_stds)):
    # ラベルを棒の外側に配置
    y_label = mean_val + std_val + 0.05 if mean_val >= 0 else mean_val - std_val - 0.12
    va = 'bottom' if mean_val >= 0 else 'top'
    ax.text(i, y_label,
            f'R²={mean_val:.3f}',
            ha='center', va=va, fontsize=11, fontweight='bold')

ax.set_xticks(x_pos)
ax.set_xticklabels(model_names, fontsize=12)
ax.set_ylabel('5分割交差検証 R²（平均 ± 標準偏差）', fontsize=11)

# y軸範囲：負のR²も見えるように
y_min = min(cv_means) - max(cv_stds) - 0.3
y_max = max(cv_means) + max(cv_stds) + 0.3
ax.set_ylim(max(-2.5, y_min), min(1.15, y_max))

ax.set_title('住宅地価格モデル：OLS vs Ridge vs Lasso\n（5分割 CV R²、東京除外、n=46）',
             fontsize=12, fontweight='bold')
ax.axhline(0, color='black', lw=1.5, ls='--', label='R²=0（ベースライン）')
ax.axhline(res_ols.rsquared, color='#888888', lw=1.2, ls=':',
           label=f'訓練データ R²={res_ols.rsquared:.3f}')
ax.legend(fontsize=10)
ax.grid(axis='y', alpha=0.3)

# 注釈：過学習の説明
note = ('負のCV-R²はモデルの過学習を示す。\n'
        'n=46に対して9変数は多すぎる。\n'
        'Ridge/Lassoは正則化でやや改善。')
ax.text(0.97, 0.03, note, transform=ax.transAxes,
        fontsize=9, color='#555555', ha='right', va='bottom',
        bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.85))

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2023_H3_fig3_cv.png'), bbox_inches='tight')
plt.close()

print("Figure 3 saved.")

# ── Figure 4: 予測値 vs 実測値（OLS, 東京除外）─────────────────
print("Figure 4: 予測値 vs 実測値...")
y_pred_ols = res_ols.fittedvalues.values
y_actual   = y_ols.values

fig, ax = plt.subplots(figsize=(8, 7))

ax.scatter(y_actual, y_pred_ols, color='#1565C0', alpha=0.75, s=70, zorder=3)

# 45度線
lims = [min(y_actual.min(), y_pred_ols.min()) * 0.9,
        max(y_actual.max(), y_pred_ols.max()) * 1.1]
ax.plot(lims, lims, 'k--', lw=1.5, label='完全予測（y=ŷ）', alpha=0.6)
ax.set_xlim(lims)
ax.set_ylim(lims)

# 主要都市ラベル
for i, row in df_no_tokyo.iterrows():
    if row[TARGET] > 80000 or y_pred_ols[i] > 80000 or \
       abs(y_actual[i] - y_pred_ols[i]) > 20000:
        ax.annotate(row['Prefecture'],
                    xy=(y_actual[i], y_pred_ols[i]),
                    xytext=(5, 5), textcoords='offset points',
                    fontsize=8, color='#555555')

# RMSE / R²
rmse = np.sqrt(np.mean((y_actual - y_pred_ols) ** 2))
r2   = res_ols.rsquared
ax.text(0.05, 0.92,
        f'R² = {r2:.4f}\nRMSE = {rmse:,.0f} 千円/m²\nAIC = {res_ols.aic:.1f}',
        transform=ax.transAxes, fontsize=11,
        verticalalignment='top',
        bbox=dict(boxstyle='round', facecolor='#EFF3FF', alpha=0.85))

ax.set_xlabel('実測値（住宅地標準価格, 千円/m²）', fontsize=12)
ax.set_ylabel('OLS 予測値（千円/m²）', fontsize=12)
ax.set_title('OLS モデルの予測精度\n（東京除外、2022年度 46都道府県）',
             fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2023_H3_fig4_pred.png'), bbox_inches='tight')
plt.close()

print("Figure 4 saved.")

# ── 結果サマリー出力 ──────────────────────────────────────────────
print("\n" + "=" * 60)
print("分析完了：地価最適モデルの構築と手法提案")
print("=" * 60)
print(f"  サンプル数: {len(df_no_tokyo)}都道府県（東京除外）")
print(f"  目的変数: 住宅地標準価格（C5401, 千円/m²）")
print()
print(f"【OLS モデル性能】")
print(f"  R² = {res_ols.rsquared:.4f}, Adj.R² = {res_ols.rsquared_adj:.4f}")
print(f"  AIC = {res_ols.aic:.2f}, BIC = {res_ols.bic:.2f}")
print()
print(f"【5分割 CV R²】")
print(f"  OLS:   {cv_ols.mean():.4f} ± {cv_ols.std():.4f}")
print(f"  Ridge: {cv_ridge.mean():.4f} ± {cv_ridge.std():.4f}")
print(f"  Lasso: {cv_lasso.mean():.4f} ± {cv_lasso.std():.4f}")
print()
print("【出力ファイル】")
print("  html/figures/2023_H3_fig1_scatter.png - 人口 vs 住宅地価格")
print("  html/figures/2023_H3_fig2_coef.png    - OLS 標準化係数")
print("  html/figures/2023_H3_fig3_cv.png      - モデル比較（CV-R²）")
print("  html/figures/2023_H3_fig4_pred.png    - 予測値 vs 実測値")

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.decomposition import FactorAnalysis

b_raw = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', header=1)
b = b_raw[(b_raw['年度'] == 2022)
          & b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].set_index('都道府県')
e_raw = pd.read_csv('data/raw/SSDSE-E-2026.csv', encoding='cp932', header=2)
e = e_raw[e_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].set_index('都道府県')

# 目的変数: 住宅地の標準価格（円/m2）
price = pd.to_numeric(b['標準価格（平均価格）（住宅地）'], errors='coerce')

# 因子分析にかける観測変数（都市の性格を表す 8 変数）
X = pd.DataFrame(index=b.index)
X['人口密度'] = b['総人口'] / pd.to_numeric(e['可住地面積'], errors='coerce')
X['転入率'] = b['転入者数（日本人移動者）'] / b['総人口'] * 100
X['大学学生数割合'] = b['大学学生数'] / b['総人口'] * 100
X['求人倍率'] = pd.to_numeric(b['月間有効求人数（一般）'], errors='coerce') / pd.to_numeric(b['月間有効求職者数（一般）'], errors='coerce')
X['宿泊者数per千人'] = pd.to_numeric(b['延べ宿泊者数'], errors='coerce') / b['総人口'] * 1000
X['高齢化率'] = b['65歳以上人口'] / b['総人口'] * 100
X['出生率'] = pd.to_numeric(b['合計特殊出生率'], errors='coerce')
X['診療所密度'] = b['一般診療所数'] / b['総人口'] * 10000

Xz = (X - X.mean()) / X.std()  # 標準化（因子分析の前提）

# 因子分析（2因子・バリマックス回転）
fa = FactorAnalysis(n_components=2, rotation='varimax', random_state=0)
scores = fa.fit_transform(Xz.values)
loadings = pd.DataFrame(fa.components_.T, index=X.columns, columns=['因子1', '因子2'])
print('=== 因子負荷量（バリマックス回転後） ===')
print(loadings.round(3))

# 因子得点で log(地価) を回帰（原論文のフロー: 因子分析 → 重回帰）
sc = pd.DataFrame(scores, index=X.index, columns=['因子1', '因子2'])
model = sm.OLS(np.log(price), sm.add_constant(sc)).fit()
print()
print('=== log(住宅地価格) を因子得点で回帰 ===')
print(f'R2 = {model.rsquared:.3f}')
for c in ['因子1', '因子2']:
    print(f'  {c}: 係数 {model.params[c]:+.3f} (p={model.pvalues[c]:.2e})')
print()
print('因子1の得点 上位5県:', '、'.join(sc['因子1'].nlargest(5).index))
print('因子2の得点 上位5県:', '、'.join(sc['因子2'].nlargest(5).index))