📌 補足セクション — 決定係数を SSDSE-B-2026 で確かめる

本セクションは「決定係数」を 47都道府県データ（SSDSE-B-2026）で具体的に確認するための追加教材です。 例として総人口で課税対象所得を説明したときの R²を扱います。

🎨 直感で掴む — 決定係数

📐 数式または定義

🧮 実値で計算してみる — 決定係数

# 最小コード（直書き）
df = pd.read_csv('data/raw/SSDSE-B-2026.csv')

import pandas as pd
import numpy as np

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', header=[0,1,2])
# 列名を 3 段ヘッダの最下段だけ採用（コード列： A1101, D3201 等）
df.columns = [c[-1] for c in df.columns]

# 2022 年の 47都道府県スナップショット
sub = df[df['年度コード'] == 2022].copy()
x = sub['A1101'].astype(float)   # 総人口
y = sub['D3201'].astype(float)   # 課税対象所得

# 決定係数の基礎統計
x_mean, y_mean = x.mean(), y.mean()
beta1 = ((x - x_mean) * (y - y_mean)).sum() / ((x - x_mean) ** 2).sum()
beta0 = y_mean - beta1 * x_mean

print(f'n = {len(x)}')             # 47
print(f'beta1 = {beta1:,.4f}')     # 傾き
print(f'beta0 = {beta0:,.4f}')     # 切片
print(f'相関係数 = {x.corr(y):.4f}')  # 0.95+ になる

# 残差・決定係数も計算
y_hat = beta0 + beta1 * x
resid = y - y_hat
ss_res = (resid ** 2).sum()
ss_tot = ((y - y_mean) ** 2).sum()
r2 = 1 - ss_res / ss_tot
print(f'R^2 = {r2:.4f}')

⚠️ 補足の落とし穴

• 都道府県データはサンプルサイズ 47 が固定：標本数を増やせないため、 統計的検定の検出力は中程度に留まります。 解釈時にこの限界を意識してください。
• スケールの違い：総人口（百万単位）と課税対象所得（百億単位）など、 桁が大きく異なる変数は標準化して扱うのが安全です。
• 東京・大阪等の影響点：少数の大都市が回帰や相関の結果を大きく動かすことがあります。 影響統計量（てこ比、 クックの距離）の確認を推奨します。

🔗 関連用語（補足リンク）

🔬 補強：R² と分散分解 — TSS = ESS + RSS を SSDSE-B-2026 で実演

R² は「目的変数の分散のうち、 モデルが説明できる割合」と定義されますが、 この分子・分母を実データで分解すると、 R² がなぜ 0〜1 を取るのか、 マイナスになり得るのか、 なぜ訓練データでは「変数を増やすと必ず R² が上がる」のかが明快になります。 SSDSE-B-2026 を使って実演します。

📐 三つの平方和

TSS = ESS + RSS
R² = ESS / TSS = 1 − RSS / TSS

• TSS（Total Sum of Squares）= Σ(yᵢ − ȳ)² … 目的変数全体の変動。 モデル無関係の「説明したい量」。
• ESS（Explained SS）= Σ(ŷᵢ − ȳ)² … モデルが説明した変動。
• RSS（Residual SS）= Σ(yᵢ − ŷᵢ)² … 残差。 モデルが説明できなかった変動。

🧮 SSDSE-B-2026 で実値計算

import pandas as pd, numpy as np
from sklearn.linear_model import LinearRegression

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='utf-8', skiprows=1)
X = df[['一般世帯数', '15歳未満人口']].dropna()
y = df.loc[X.index, '人口（総数）']

model = LinearRegression().fit(X, y)
y_hat = model.predict(X)
y_mean = y.mean()

TSS = ((y - y_mean) ** 2).sum()
ESS = ((y_hat - y_mean) ** 2).sum()
RSS = ((y - y_hat) ** 2).sum()

print(f'TSS = {TSS:.0f}')
print(f'ESS = {ESS:.0f}')
print(f'RSS = {RSS:.0f}')
print(f'TSS - (ESS+RSS) = {TSS-ESS-RSS:.2f}  (理論上 0)')
print(f'R² = ESS/TSS = {ESS/TSS:.4f}')
print(f'R² = 1-RSS/TSS = {1-RSS/TSS:.4f}')
# 出力例: R² ≈ 0.999（人口は世帯数でほぼ完全予測される）

📊 同じデータで「変数を 1 つずつ追加」して R² 上昇を観察

from sklearn.linear_model import LinearRegression

base_cols = ['一般世帯数']
add_cols = ['15歳未満人口', '65歳以上人口', '出生数', '死亡数']

for i in range(len(add_cols) + 1):
    cols = base_cols + add_cols[:i]
    Xi = df[cols].dropna()
    yi = df.loc[Xi.index, '人口（総数）']
    r2 = LinearRegression().fit(Xi, yi).score(Xi, yi)
    print(f'説明変数 {len(cols)} 個: R² = {r2:.6f}')
# R² は単調非減少。 無関係な変数を入れても上がるのが落とし穴。

⚠️ R² がマイナスになる場合 — 何が起きているか

R² < 0 は、 数式上は「モデルが平均値で予測するより悪い」を意味します。 ホールドアウト検証や、 切片無しモデル、 過学習で大きく外れた検証データで起こりがちです。

from sklearn.model_selection import train_test_split

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, shuffle=False)
model_bad = LinearRegression(fit_intercept=False).fit(X_tr, y_tr)
print('test R² (切片なし):', model_bad.score(X_te, y_te))
# 切片無しで強引にフィットすると test R² がマイナスになり得る

📐 R² と相関係数の関係 — 単回帰のみ

単回帰では R² = r²（r は説明変数 X と目的変数 Y の相関係数の絶対値）。 ただし重回帰では R² ≠ Σrⱼ² なので、 「個別変数の貢献」と R² 全体は直接対応しません（→ 標準化β や偏相関で評価）。

💡 使い分け：モデル比較なら調整 R² や AIC、 報告書は R² 単独で OK、 予測誤差を生の単位で見たいなら RMSE、 と目的別に並走させましょう。

🌐 補強2：R² の派生・代替指標 — McFadden / Nagelkerke / Cox-Snell / 調整 R²

古典 R² は OLS 線形回帰のための指標ですが、 ロジスティック回帰やポアソン回帰など 非線形・非ガウス モデルでも「説明力」を測りたい場面は多い。 そこで考案されたのが「擬似 R²」（pseudo R²）と呼ばれる代替指標群です。 主要 4 つを比較整理します。

🔬 「対数尤度」に基づく擬似 R² の発想

非ガウスモデルでは「残差平方和」が定義しにくいため、 代わりに 対数尤度 を使います。 完全モデル（M）と帰無モデル（切片のみ、0）の対数尤度を比べ、 「モデルが帰無モデルからどれだけ改善したか」を尤度の比で測るのが擬似 R² の核心です。

🐍 SSDSE-B-2026 で擬似 R² を計算

import pandas as pd
import statsmodels.api as sm
import numpy as np

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='utf-8', skiprows=1)
df = df.dropna(subset=['人口（総数）', '一般世帯数', '65歳以上人口'])

# 二値ターゲット：「高齢化率が中央値以上か」
df['high_aging'] = ((df['65歳以上人口'] / df['人口（総数）']) > 0.30).astype(int)

X = sm.add_constant(df[['一般世帯数']])
y = df['high_aging']

logit = sm.Logit(y, X).fit(disp=0)
ll_M = logit.llf
ll_0 = logit.llnull
n = len(y)

mcfadden = 1 - ll_M / ll_0
cox_snell = 1 - np.exp(2 * (ll_0 - ll_M) / n)
max_cox = 1 - np.exp(2 * ll_0 / n)
nagelkerke = cox_snell / max_cox

print(f'McFadden R²  = {mcfadden:.4f}')
print(f'Cox-Snell R² = {cox_snell:.4f}')
print(f'Nagelkerke R² = {nagelkerke:.4f}')
# 値の絶対比較は要注意：McFadden 0.2-0.4 で「良いフィット」とされる

⚠️ 擬似 R² の「絶対値での比較」は危険

古典 R² が 0.5 と擬似 R² が 0.5 では意味がまったく違います。 McFadden の 0.2 ≈ 古典 R² の 0.5 と言われることもあり、 「擬似」と付くだけあって直接比較できません。 報告時は必ず「どの擬似 R² か」を明記しましょう。

📐 R² の信頼区間 — ブートストラップで推定

古典 R² の標本分布は複雑なので、 信頼区間はブートストラップで求めるのが実務的です。

from sklearn.linear_model import LinearRegression
import numpy as np

df2 = df.dropna(subset=['人口（総数）', '一般世帯数']).copy()
X2 = df2[['一般世帯数']].values
y2 = df2['人口（総数）'].values

r2_boot = []
rng = np.random.default_rng(42)
for _ in range(1000):
    idx = rng.integers(0, len(y2), len(y2))
    r2 = LinearRegression().fit(X2[idx], y2[idx]).score(X2[idx], y2[idx])
    r2_boot.append(r2)

print(f'R² 95% CI: [{np.percentile(r2_boot, 2.5):.4f}, {np.percentile(r2_boot, 97.5):.4f}]')
# 都道府県データのような小標本では区間が広くなることに注意

💡 論文掲載時のチェック：(1) R² 単独で結論しない、 (2) 調整 R² も併記、 (3) 非線形モデルなら擬似 R² と尤度比検定、 (4) 小標本ならブートストラップ CI、 を意識すると説得力が増します。

🔬 補強3：R² の挙動をシミュレーション実験で理解する

R² が「サンプルサイズ」「ノイズの大きさ」「変数の数」でどう動くかを、 SSDSE-B-2026 を素材にした擬似実験で確認します。 教科書の数式だけでは見えない「現場感覚」を養うのが目的です。

🧪 実験 1：ノイズを増やすと R² はどう減るか

import pandas as pd, numpy as np
from sklearn.linear_model import LinearRegression

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='utf-8', skiprows=1)
df = df.dropna(subset=['人口（総数）', '一般世帯数'])

X = df[['一般世帯数']].values
y_clean = df['人口（総数）'].values

print('ノイズ標準偏差倍率 → R²')
y_std = y_clean.std()
rng = np.random.default_rng(0)
for k in [0, 0.1, 0.3, 0.5, 1.0, 2.0, 5.0]:
    noise = rng.normal(0, k * y_std, size=len(y_clean))
    y_noisy = y_clean + noise
    r2 = LinearRegression().fit(X, y_noisy).score(X, y_noisy)
    print(f'  {k:.1f}σ : R² = {r2:.4f}')
# ノイズ k=1σ で R² が半分、 k=5σ でほぼ 0 に

🧪 実験 2：「無関係な変数」を増やすと訓練 R² と汎化 R² はどう乖離するか

import pandas as pd, numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='utf-8', skiprows=1)
df = df.dropna(subset=['人口（総数）', '一般世帯数'])

X = df[['一般世帯数']].copy().reset_index(drop=True)
y = df['人口（総数）'].reset_index(drop=True)

rng = np.random.default_rng(42)
print('追加無関係変数の数 → 訓練R² / CV R²')
for n_extra in [0, 5, 10, 20, 30, 40]:
    Xx = X.copy()
    for j in range(n_extra):
        Xx[f'noise_{j}'] = rng.normal(0, 1, size=len(Xx))
    train_r2 = LinearRegression().fit(Xx, y).score(Xx, y)
    cv_r2s = []
    for tr, te in KFold(n_splits=5, shuffle=True, random_state=1).split(Xx):
        m = LinearRegression().fit(Xx.iloc[tr], y.iloc[tr])
        cv_r2s.append(m.score(Xx.iloc[te], y.iloc[te]))
    print(f'  +{n_extra} 個: 訓練 R² = {train_r2:.4f}, CV R² = {np.mean(cv_r2s):.4f}')
# 訓練 R² は単調増加、 CV R² は途中から低下 = 過学習の証拠

🧪 実験 3：標本サイズが小さいと R² は過大評価される

import pandas as pd, numpy as np
from sklearn.linear_model import LinearRegression

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='utf-8', skiprows=1)
df = df.dropna(subset=['人口（総数）', '一般世帯数'])

print('標本サイズ → R² のブートストラップ平均と分散')
rng = np.random.default_rng(0)
for n in [5, 10, 20, 47]:
    r2s = []
    for _ in range(500):
        idx = rng.choice(len(df), n, replace=True)
        m = LinearRegression().fit(df[['一般世帯数']].iloc[idx],
                                    df['人口（総数）'].iloc[idx])
        r2s.append(m.score(df[['一般世帯数']].iloc[idx],
                           df['人口（総数）'].iloc[idx]))
    print(f'  n={n:3d}: mean R² = {np.mean(r2s):.4f}, std = {np.std(r2s):.4f}')
# 小標本では R² が高めに出るが分散も大きい

📊 実験結果のまとめ表

🌐 R² と関連指標の使い分け早見表

💡 シミュレーション思考の効用：教科書を読むだけでなく、 自分でデータを動かしてみることで、 R² が「単なる数値」でなく「データ構造の鏡」であることが体感できます。 実験コードを写経して、 値を変えてみることをおすすめします。

📚 補強4：R² を扱う実務 TIPS と参考文献

R² は誰でも知っている指標ですが、 「正しく報告する」ためのチェックリストは意外と整理されていません。 ここでは論文掲載・実務報告で役立つ実践 TIPS を整理します。

📋 R² 報告チェックリスト

🐍 報告に使える「R² ダッシュボード」関数

import pandas as pd, numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.metrics import mean_squared_error, mean_absolute_error

def r2_dashboard(X, y, model=None):
    # R² 報告に必要な指標を一覧で返す
    if model is None:
        model = LinearRegression()
    n, p = X.shape
    model.fit(X, y)
    y_hat = model.predict(X)
    r2 = model.score(X, y)
    adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)
    cv_r2 = cross_val_score(model, X, y, cv=5, scoring='r2').mean()
    rmse = np.sqrt(mean_squared_error(y, y_hat))
    mae = mean_absolute_error(y, y_hat)
    return {
        'n':       n,
        'p':       p,
        'R²':      round(r2, 4),
        '調整 R²': round(adj_r2, 4),
        'CV R²':   round(cv_r2, 4),
        'RMSE':    round(rmse, 2),
        'MAE':     round(mae, 2),
    }

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='utf-8', skiprows=1)
df = df.dropna(subset=['人口（総数）', '一般世帯数', '15歳未満人口'])
X = df[['一般世帯数', '15歳未満人口']]
y = df['人口（総数）']
dashboard = r2_dashboard(X, y)
for k, v in dashboard.items():
    print(f'  {k}: {v}')

📐 R² の業界別「目安」

「R² が低いから悪いモデル」と短絡しないこと。 分野によっては R²=0.3 でも価値ある発見になり得ます。 大切なのは「同分野のベースラインと比べてどうか」です。

📚 主要参考文献

• Cohen, J. (1988) Statistical Power Analysis for the Behavioral Sciences — R² と効果量の関係。
• Kvalseth, T.O. (1985) "Cautionary Note about R²" — R² の落とし穴を 7 つ整理した古典。
• Nakagawa & Schielzeth (2013) — 混合モデルの R²（限界 R² と条件付き R²）。
• Anderson-Sprecher (1994) — 同じ R² でも意味が違う 7 ケース。

💡 結論：R² は「便利だが万能ではない」指標。 報告するときは必ず「どんな R² か」「何と比べたか」を明示し、 他の指標と併走させることで、 説得力のある分析になります。

🌐 補強5：R² と ANOVA — F 統計量と分散分析表で読み解く

R² は ANOVA（分散分析）と直結します。 OLS の標準的な「F 検定」は実は「R² の有意性検定」と数学的に等価。 SSDSE-B-2026 で並べて確認します。

📐 R² と F の関係

F = [R² / p] / [(1 − R²) / (n − p − 1)]

つまり R² が大きいほど F も大きく、 帰無仮説「すべての係数 = 0」が棄却されやすくなります。

🐍 SSDSE-B-2026 で F と R² を比較

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

df = pd.read_csv("data/raw/SSDSE-B-2026.csv", encoding="utf-8", skiprows=1)
df.columns = [c.replace("（", "_").replace("）", "_") for c in df.columns]
df = df.dropna(subset=["人口_総数_", "一般世帯数", "15歳未満人口"])

model = smf.ols("人口_総数_ ~ 一般世帯数 + Q("15歳未満人口")", data=df).fit()
print(model.summary())
# R² と F 統計量が表示される。 両者は同一情報の別表現。
anova = sm.stats.anova_lm(model)
print(anova)
# 分散分析表：説明変数ごとの SS（平方和）と F 値

📊 分散分析表の読み方

💡 R² ＝ ANOVA の使い分け

• R²：説明力の「割合」を一目で見たいとき
• F 検定：「モデル全体が有意か」を検定したいとき
• 個別 t 検定：「各係数が有意か」を判定したいとき
• 分散分析表：各変数の「相対寄与」を比較したいとき

💡 使い分けのコツ：R² は「全体の説明力」、 ANOVA は「変数ごとの相対寄与」、 t 検定は「個別の有意性」。 3 つを併走させると、 モデルの全体像が立体的に見えてきます。

📚 補強6：R² FAQ — 7 つの本音質問

Q1. R² が 0.95 なのに予測が外れます。 なぜ？

訓練データの R² は内側のフィット。 ホールドアウトテストデータの R² を確認してください。 過学習の可能性大。

Q2. R² がマイナスになりました。 どうすれば？

テストデータでは普通に起こります。 切片なしモデル、 訓練と分布が違うテストデータ、 過学習などが原因。

Q3. R² 0.3 はモデルが悪い？

分野次第。 心理学・社会科学なら良好、 物理学なら不十分。 同分野の典型 R² と比較しましょう。

Q4. R² と相関係数 r の二乗は同じ？

単回帰では r² = R²。 重回帰では別物。 重回帰で「個別変数の貢献率」を測りたいなら偏 R² または標準化β。

Q5. 説明変数を増やすと必ず R² が上がるのですか？

はい、 訓練データの R² は単調非減少。 だから「変数を増やすたびに R² が上がった」だけでは進歩の証拠になりません。 調整 R² や CV R² を使いましょう。

Q6. ロジスティック回帰で R² を計算したい

McFadden、 Cox-Snell、 Nagelkerke のいずれか。 「擬似 R²」と総称されますが、 古典 R² と数値を直接比較できない点に注意。

Q7. R² が高いから因果関係がある？

違います。 R² は「説明力」であって「因果」ではない。 高い R² + 強い理論的根拠 + 内生性対策があってはじめて因果に近づきます。