論文一覧に戻る 🗺 概念マップ 統計データ分析コンペ 教育用再現集
2025年 統計データ分析コンペティション ｜ 審査員奨励賞 [高校生の部]

消滅可能性自治体を救うには

⏱️ 推定読了時間: 約30分
早稲田大学系属 早稲田実業学校
📝 3行で分かる要約

目次

  1. 研究概要と背景
  2. データ：SSDSE-A 521自治体
  3. 無相関検定
  4. バックワード変数選択による重回帰
  5. 事例分析：川島町
  6. 外れ値の箱ひげ図診断
  7. まとめ
  8. 📥 データの準備
  9. 💼 実社会での応用
  10. ⚠️ よくある誤解
  11. 📖 用語集
  12. 📐 手法ガイド
  13. 🚀 発展の可能性
  14. 🎯 自分でやってみよう
  15. 🤔 Q&A

🎯 この記事を読むと何ができるようになるか

📥 データの準備（再現コードを動かす前に）

このページの分析を自分で再現するには、以下の手順でデータを準備してください。コードの編集は不要です。

1
データをダウンロードする 統計センターの SSDSE 配布ページから、以下のファイルをダウンロードします。
SSDSE-A-2025.csv　← SSDSE-A（市区町村データ）📥 直接DL
⬇ SSDSEダウンロードページを開く
2
ファイルを所定のフォルダに配置する ダウンロードしたCSVを、プロジェクトの data/raw/ フォルダに入れます。
2026 統計・データ解析コンペ/ ├── code/ │ └── 2025_H5_3_shorei.py ← 実行するスクリプト └── data/ └── raw/ SSDSE-A-2025.csv ← ここに置く
3
スクリプトをそのまま実行する ターミナルでプロジェクトルートに移動し、以下を実行します。
python3 code/2025_H5_3_shorei.py
図は html/figures/ に自動保存されます。
概
研究概要と背景

2014年の「消滅可能性都市」リスト公表以降、地方自治体の人口維持は日本の最重要政策課題の一つとなっている。「消滅可能性」の基準は20〜39歳の若年女性人口の減少率であり、本研究はこれを目的変数として統計的分析を行った。

「消滅可能性自治体」とは 2010〜2040年の30年間で「20〜39歳の若年女性人口が50%以上減少する」と推計される自治体。消滅可能性が高いのは人口0.5〜2万人の小規模自治体に多く、本研究の分析対象。
分析の流れ
SSDSE-A
521自治体
抽出
→
無相関検定
（候補変数）
→
バックワード
変数選択
(全p<0.05)
→
川島町
事例分析

SSDSE-A 無相関検定 バックワード変数選択 事例分析

データ：SSDSE-A 521自治体

データ抽出の手順

SSDSE（社会・人口統計体系データセット）-A は市区町村レベルの統計データを収録する。全1728自治体から「人口0.5万〜2万人」の521自治体を抽出。

条件自治体数理由
全自治体1728SSDSE-A 全データ
人口0.5万〜2万人521消滅可能性が特に高い規模帯
欠損値除外後521（合成）本教育用コードでは合成データ使用

説明変数候補（人口1万人あたり）

カテゴリ変数想定効果
保育・教育保育所数正（子育て環境が良いと若年女性が残る）
幼稚園数正
雇用製造業従業者割合正（安定した雇用）
農業従業者割合負（高齢化・低賃金）
福祉生活保護率負（経済的困窮）
医療・福祉老人ホーム定員数、病院病床数？
住環境住宅地平均地価正（経済的活力の代理）
やってみよう■ 実データ読み込み（SSDSE-A-2025: 市区町村データ）
📝 コード
1
2
3
4
5
6
df_a_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-A-2025.csv'),
    encoding='cp932', header=1
)
df_a = df_a_raw.iloc[1:].copy()
df_a.columns = df_a_raw.iloc[0].values
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • pd.read_csv(...) でCSVを読み込みます。encoding='cp932' は日本語Windows由来の文字コード、header=1 は「2行目を列名として使う」。
💡 Python TIPS df['A'] / df['B'] — pandasの列同士の四則演算は要素ごと（element-wise）。forループ不要なのが強み。
やってみよう■ 実データ読み込み（SSDSE-A-2025: 市区町村データ） — 数値変換
📝 コード
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
# 数値変換
num_cols = [
    '総人口', '総人口（女）', '15歳未満人口（女）', '15～64歳人口（女）',
    '65歳以上人口（女）', '75歳以上人口（女）',
    '就業者数', '就業者数（女）',
    '第１次産業就業者数', '第２次産業就業者数', '第３次産業就業者数',
    '従業者数（民営）（製造業）',
    '保育所等数（基本票）', '幼稚園数', '小学校数', '中学校数', '高等学校数',
    '一般病院数', '一般診療所数', '医師数',
    '農家数（販売農家）', '農家数（自給的農家）',
]
for c in num_cols:
    if c in df_a.columns:
        df_a[c] = pd.to_numeric(df_a[c], errors='coerce')
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • このステップでは前のステップで作ったデータを加工しています。コードを上から順に読んでみてください。
💡 Python TIPS Seriesの .map() は「1対1の置き換え」、.apply() は「関数を当てる」。辞書なら .map()、ロジックなら .apply()。
やってみよう■ 実データ読み込み（SSDSE-A-2025: 市区町村データ） — 人口5,000〜20,000人の自治体に絞る（消滅可能性自治体の規模）
📝 コード
21
22
23
24
25
26
27
28
# 人口5,000〜20,000人の自治体に絞る（消滅可能性自治体の規模）
df = df_a[
    (df_a['総人口'] >= 5000) & (df_a['総人口'] <= 20000)
].copy().reset_index(drop=True)

print("=" * 60)
print(f"■ 対象自治体数: N={len(df)}（人口5,000〜20,000人）")
print("=" * 60)
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • このステップでは前のステップで作ったデータを加工しています。コードを上から順に読んでみてください。
💡 Python TIPS [式 for x in リスト] はリスト内包表記。forループでappendする代わりに1行でリストを作れます。
やってみよう■ 実データ読み込み（SSDSE-A-2025: 市区町村データ） — ── 特徴量エンジニアリング ─────────────────────────────────────
📝 コード
29
30
31
32
33
34
35
# ── 特徴量エンジニアリング ─────────────────────────────────────
pop = df['総人口'].values.clip(1)
pop_f = df['総人口（女）'].values.clip(1)
emp  = df['就業者数'].values.clip(1)

# 目的変数: 15〜64歳女性人口比率（%）— 消滅可能性に直結する指標
y_vals = df['15～64歳人口（女）'].values / pop_f * 100
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • このステップでは前のステップで作ったデータを加工しています。コードを上から順に読んでみてください。
💡 Python TIPS r, p = stats.pearsonr(...) — Pythonは複数戻り値を同時に受け取れる（タプルアンパック）。
やってみよう■ 実データ読み込み（SSDSE-A-2025: 市区町村データ） — 説明変数（人口1万人あたり等）
📝 コード
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
# 説明変数（人口1万人あたり等）
nursery       = df['保育所等数（基本票）'].values / pop * 10000    # 保育所数（万人対）
kindergarten  = df['幼稚園数'].values / pop * 10000               # 幼稚園数（万人対）
elementary    = df['小学校数'].values / pop * 10000               # 小学校数（万人対）
junior_high   = df['中学校数'].values / pop * 10000               # 中学校数（万人対）
mfg_ratio     = df['従業者数（民営）（製造業）'].values / emp * 100 # 製造業従業者割合（%）
agri_ratio    = df['第１次産業就業者数'].values / emp * 100        # 農業従業者割合（%）
female_emp_rate = df['就業者数（女）'].values / df['15～64歳人口（女）'].values.clip(1) * 100  # 女性就業率
hospital      = df['一般病院数'].values / pop * 10000             # 一般病院数（万人対）
clinic        = df['一般診療所数'].values / pop * 10000           # 一般診療所数（万人対）
farmer_rate   = df['農家数（販売農家）'].values / pop * 10000      # 販売農家数（万人対）

VAR_NAMES = [
    '保育所数', '幼稚園数', '小学校数', '中学校数',
    '製造業従業者割合', '農業従業者割合', '女性就業率',
    '一般病院数', '一般診療所数', '農家比率'
]

X_raw = np.column_stack([
    nursery, kindergarten, elementary, junior_high,
    mfg_ratio, agri_ratio, female_emp_rate,
    hospital, clinic, farmer_rate
])
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • このステップでは前のステップで作ったデータを加工しています。コードを上から順に読んでみてください。
💡 Python TIPS x if cond else y は三項演算子。リスト内包表記と組み合わせると、forとifを1行で書けます。
やってみよう■ 実データ読み込み（SSDSE-A-2025: 市区町村データ） — 欠損を含む行を除外
📝 コード
59
60
61
62
63
64
65
66
# 欠損を含む行を除外
valid = np.all(np.isfinite(X_raw), axis=1) & np.isfinite(y_vals)
X_raw = X_raw[valid]
y_vals = y_vals[valid]
df_valid = df[valid].reset_index(drop=True)

N = len(y_vals)
print(f"有効自治体数（欠損除外後）: N={N}")
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • このステップでは前のステップで作ったデータを加工しています。コードを上から順に読んでみてください。
💡 Python TIPS df[col]（1列）と df[[col1, col2]]（複数列）でカッコの数が違います。リストを渡していると覚えるとミスを減らせます。
やってみよう■ 実データ読み込み（SSDSE-A-2025: 市区町村データ） — ── 川島町（埼玉県）のインデックスを特定 ─────────────────────────────
📝 コード
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
# ── 川島町（埼玉県）のインデックスを特定 ─────────────────────────────
kawashima_mask = df_valid['市区町村'].str.contains('川島', na=False)
kawashima_idx = df_valid[kawashima_mask].index.tolist()
has_kawashima = len(kawashima_idx) > 0
if has_kawashima:
    kaw_idx = kawashima_idx[0]
    print(f"川島町: index={kaw_idx}, 15-64歳女性比率={y_vals[kaw_idx]:.2f}%")
    print(f"  保育所数={nursery[valid][kaw_idx]:.2f}, 製造業割合={mfg_ratio[valid][kaw_idx]:.2f}%")
else:
    print("川島町は対象外（人口範囲外または欠損）")

df_main = pd.DataFrame(X_raw, columns=VAR_NAMES)
df_main['若年女性比率'] = y_vals

print()
print(df_main[['若年女性比率'] + VAR_NAMES].describe().round(2))
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • .describe() — 件数・平均・標準偏差・四分位・最大/最小を一括計算。データの素性チェックに必須。
💡 Python TIPS s[:-n]「末尾n文字を除く」／s[n:]「先頭n文字を除く」。スライス [start:stop:step] はリスト・タプル・文字列共通の基本ワザです。
やってみよう■ 実データ読み込み（SSDSE-A-2025: 市区町村データ） — IQR外れ値チェック
📝 コード
83
84
85
86
87
88
89
90
91
# IQR外れ値チェック
Q1 = df_main['若年女性比率'].quantile(0.25)
Q3 = df_main['若年女性比率'].quantile(0.75)
IQR = Q3 - Q1
lower_iqr = Q1 - 1.5 * IQR
upper_iqr = Q3 + 1.5 * IQR
mask_iqr = (df_main['若年女性比率'] >= lower_iqr) & (df_main['若年女性比率'] <= upper_iqr)
n_out = (~mask_iqr).sum()
print(f"\n  IQR法外れ値: {n_out}件（IQR={IQR:.2f}, 範囲=[{lower_iqr:.2f}, {upper_iqr:.2f}]）")
▼ 実行結果
============================================================
■ 対象自治体数: N=523（人口5,000〜20,000人）
============================================================
有効自治体数（欠損除外後）: N=523
川島町: index=186, 15-64歳女性比率=52.46%
  保育所数=1.03, 製造業割合=53.17%

       若年女性比率    保育所数    幼稚園数    小学校数  ...   女性就業率   一般病院数  一般診療所数     農家比率
count  523.00  523.00  523.00  523.00  ...  523.00  523.00  523.00   523.00
mean    48.08    3.21    0.72    3.30  ...   90.17    0.78    6.95   344.71
std      4.91    1.60    1.21    2.18  ...    9.89    0.82    2.81   245.41
min     35.27    0.00    0.00    0.00  ...   54.73    0.00    1.38     0.00
25%     44.96    1.93    0.00    1.76  ...   83.56    0.00    4.97   144.79
50%     48.05    2.99    0.00    2.87  ...   90.83    0.69    6.72   290.98
75%     51.59    4.21    1.17    4.13  ...   96.79    1.29    8.30   517.94
max     63.22    9.74   13.03   14.69  ...  130.03    7.71   18.16  1321.37

[8 rows x 11 columns]

  IQR法外れ値: 1件（IQR=6.63, 範囲=[35.01, 61.53]）
💡 解説
  • このステップでは前のステップで作ったデータを加工しています。コードを上から順に読んでみてください。
💡 Python TIPS np.cumsum(arr) は累積和、np.linspace(a, b, n) は「aからbを等間隔でn個」。NumPyの定石です。
3. 無相関検定
1
無相関検定

まず若年女性人口増加と関係しそうな指標を二変量レベルで洗い出すことが有効だと考えられる。 その理由は保育・雇用・福祉など候補変数が多く、いきなり重回帰を組むと過適合のリスクがあるからである。 ここでは521自治体という大標本の利点に着目し、ピアソン相関と無相関検定を用いる。 保育所数・製造業従業者割合といった社会基盤指標が有意な正相関を示す結果が期待される。

各説明変数候補と若年女性人口増加率（目的変数）の相関係数を算出し、帰無仮説「相関なし」を検定する。有意な相関が認められた変数を回帰分析の候補とする。

相関係数棒グラフ
図2：各変数と若年女性人口増加率の相関係数。赤は有意（p<0.05）な変数。
📌 この回帰係数プロットの読み方
このグラフは
重回帰分析の各説明変数の係数（影響の強さと向き）をバーや点で表したグラフ。
読み方
右（プラス方向）に伸びるバーは「この変数が増えると目的変数も増える」正の影響。左（マイナス方向）は逆。
なぜそう解釈できるか
エラーバー（誤差棒）が0をまたいでいない変数が統計的に有意（p < 0.05）。バーが長いほど影響が大きい。
有意な相関が見られた変数
  • 保育所数（正）：子育て支援が若年女性の定着に貢献
  • 製造業従業者割合（正）：安定した雇用機会の重要性
  • 生活保護率（負）：貧困環境からの離脱傾向

DS LEARNING POINT 1

N=521での無相関検定の検出力

検出力（power）とはサンプルサイズN が大きいほど小さな相関でも統計的に有意になる。N=521では|r|≥0.09程度でもp<0.05になることがある。このため「統計的有意」と「実質的意義（effect size）」は別途確認が必要。

from scipy import stats # N=521での検出力 N = 521 r, p = stats.pearsonr(X[:, 0], y) # Cohen(1988)の基準: r=0.1(小), 0.3(中), 0.5(大) if abs(r) >= 0.5: effect = "大（Large）" elif abs(r) >= 0.3: effect = "中（Medium）" elif abs(r) >= 0.1: effect = "小（Small）" else: effect = "微小" print(f"r={r:.4f}, p={p:.4f}, 効果量: {effect}")
やってみよう■ 図の生成（4枚）
📝 コード
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
fig1, ax1 = plt.subplots(figsize=(9, 5))
corrs  = [stats.pearsonr(X_raw[:, i], y_vals)[0] for i in range(len(VAR_NAMES))]
pvals2 = [stats.pearsonr(X_raw[:, i], y_vals)[1] for i in range(len(VAR_NAMES))]
cols2  = ['#C62828' if p < 0.05 else '#9E9E9E' for p in pvals2]
y_pos2 = np.arange(len(VAR_NAMES))

ax1.barh(y_pos2, corrs, color=cols2, alpha=0.75, edgecolor='white', height=0.6)
ax1.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax1.set_yticks(y_pos2)
ax1.set_yticklabels(VAR_NAMES, fontsize=10)
ax1.set_xlabel('ピアソン相関係数', fontsize=11)
ax1.set_title('各変数と若年女性比率の相関係数\n（無相関検定: 赤=p<0.05）',
              fontsize=12, fontweight='bold')
ax1.invert_yaxis()
ax1.grid(axis='x', alpha=0.3)
for i, (r, p) in enumerate(zip(corrs, pvals2)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    ax1.text(r + (0.005 if r >= 0 else -0.005), i, f' {r:.3f}{sig}',
             va='center', ha='left' if r >= 0 else 'right', fontsize=8.5)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_H5_3_fig1_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2025_H5_3_fig1_corr.png")
▼ 実行結果
図1保存: 2025_H5_3_fig1_corr.png
💡 解説
  • fig, ax = plt.subplots(...) — 図全体（fig）と軸（ax）を作る定番。以降は ax.bar(...) 等で操作。
  • ax.axhline / ax.axvline — 水平／垂直の点線。平均線や基準線として定番。
  • stats.pearsonr(x, y) — Pearson相関係数 r と p値を同時に返します。
💡 Python TIPS x if cond else y は三項演算子。リスト内包表記と組み合わせると、forとifを1行で書けます。
2
バックワード変数選択による重回帰

前節の保育所数・製造業従業者割合・生活保護率が有意な相関を示した結果を踏まえると、 これら要因が互いに絡みつつ若年女性の定着を左右すると考えられる。 これを検証する必要があるが、その手法としてp値基準のバックワード変数選択を伴う重回帰に着目した。 交絡を統制した上でも保育・雇用の効果が頑健に残り、政策的処方箋を導ける結果が期待される。

有意な相関が確認された変数を全て投入したOLS重回帰から出発し、最も有意でない変数（p値が最大）を順次除外する。全変数のp<0.05になるまで繰り返す。

y_i = β₀ + Σ βⱼxᵢⱼ + εᵢ

除外基準: argmax{p値} の変数を除外（p≥0.05の間繰り返す）
バックワード選択後の回帰係数
図3：バックワード変数選択後の回帰係数（全p<0.05）。保育所数・製造業従業者割合が最も強い効果を示す。
📌 この回帰係数プロットの読み方
このグラフは
重回帰分析の各説明変数の係数（影響の強さと向き）をバーや点で表したグラフ。
読み方
右（プラス方向）に伸びるバーは「この変数が増えると目的変数も増える」正の影響。左（マイナス方向）は逆。
なぜそう解釈できるか
エラーバー（誤差棒）が0をまたいでいない変数が統計的に有意（p < 0.05）。バーが長いほど影響が大きい。
バックワード選択の結果 全変数のp<0.05となるモデルへの収束後、以下の変数が残存：
  • 保育所数（正）：保育サービスの充実が最重要
  • 製造業従業者割合（正）：雇用の安定性
  • 生活保護率（負）：貧困環境は定着を妨げる

DS LEARNING POINT 2

バックワード変数選択の実装

p値基準のバックワード選択は解釈が簡単だが、変数間の相関が高い場合に不安定になる欠点がある。AIC基準（前セクション参照）と組み合わせると より頑健。

import statsmodels.api as sm current_vars = list(range(n_vars)) while True: X_cur = sm.add_constant(X[:, current_vars]) model = sm.OLS(y, X_cur).fit() pvals = model.pvalues[1:] # 定数項を除く max_p = pvals.max() if max_p < 0.05: break # 全変数が有意 → 終了 # 最も有意でない変数を削除 worst_idx = pvals.argmax() removed_var = current_vars.pop(worst_idx) print(f"除外: {VAR_NAMES[removed_var]}（p={max_p:.4f}）") print(f"最終変数: {[VAR_NAMES[i] for i in current_vars]}")
やってみよう図図2 (fig2_backward): バックワード選択後の回帰係数
📝 コード
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
fig2, ax2 = plt.subplots(figsize=(9, max(4, len(current_vars) * 0.8)))
n_sel  = len(current_vars)
coefs3 = np.asarray(final_model.params[1:])
ses3   = np.asarray(final_model.bse[1:])
pvals3 = np.asarray(final_model.pvalues[1:])
cols3  = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E' for p in pvals3]
y_pos3 = np.arange(n_sel)

ax2.barh(y_pos3, coefs3, color=cols3, alpha=0.75, edgecolor='white', height=0.6)
ax2.errorbar(coefs3, y_pos3, xerr=1.96*ses3, fmt='none', color='#333',
             capsize=4, linewidth=1.2)
ax2.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax2.set_yticks(y_pos3)
ax2.set_yticklabels(selected_names, fontsize=10)
ax2.set_xlabel('回帰係数（±1.96SE）', fontsize=11)
ax2.set_title(f'バックワード変数選択後の重回帰係数\n'
              f'R²={final_model.rsquared:.3f}（全p<0.05）, N={N}自治体',
              fontsize=12, fontweight='bold')
ax2.invert_yaxis()
ax2.grid(axis='x', alpha=0.3)
from matplotlib.patches import Patch
ax2.legend(handles=[Patch(color='#C62828', alpha=0.75, label='p<0.01'),
                    Patch(color='#FF8F00', alpha=0.75, label='p<0.05')], fontsize=9)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_H5_3_fig2_backward.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2025_H5_3_fig2_backward.png")
▼ 実行結果
図2保存: 2025_H5_3_fig2_backward.png
💡 解説
  • import pandas as pd など — 必要なライブラリをまとめて呼び出します。as pd は短い別名（alias）。
  • fig, ax = plt.subplots(...) — 図全体（fig）と軸（ax）を作る定番。以降は ax.bar(...) 等で操作。
  • ax.axhline / ax.axvline — 水平／垂直の点線。平均線や基準線として定番。
💡 Python TIPS df[col]（1列）と df[[col1, col2]]（複数列）でカッコの数が違います。リストを渡していると覚えるとミスを減らせます。
3
事例分析：川島町

前節の保育所数・製造業従業者割合が増加に正、生活保護率が負の効果を持つという結果を踏まえると、 具体的な成功自治体ではこれらが揃って高水準であると考えられる。 これを検証する必要があるが、その手法として代表事例（川島町）の値を全国分布上にプロットする事例分析に着目した。 モデルの示す処方箋と実在する自治体の数値が一致する結果が期待される。

分析の最終段階として、統計モデルで特定された要因が実際の成功事例に当てはまるかを川島町（埼玉県）を例に確認する。

保育所数・製造業と若年女性増加率の散布図
図4：保育所数（左）・製造業従業者割合（右）と若年女性人口増加率の関係。川島町（赤星）は両変数で高い値を示す。
指標川島町（参考）平均分析結果との整合
保育所数（人口1万対）高い3.0程度一致（保育充実）
製造業従業者割合高い（~45%）22%一致（雇用安定）
若年女性人口増加率正（増加）負（減少傾向）一致
川島町が示す成功要因 川島町は首都圏へのアクセスが良く、製造業の雇用が豊富で、保育環境も整備されている。統計モデルが指摘する「保育・雇用の充実」という処方箋を体現した事例として解釈できる。
やってみよう図図3 (fig3_box): 若年女性比率の分布（箱ひげ図 + ヒストグラム）
📝 コード
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
fig3, axes3 = plt.subplots(1, 2, figsize=(12, 5))
fig3.suptitle(f'若年女性比率（15〜64歳女性/総女性人口）の分布\n'
              f'人口5,000〜20,000人 {N}自治体（SSDSE-A 実データ）',
              fontsize=12, fontweight='bold')

ax3a = axes3[0]
bp = ax3a.boxplot(df_main['若年女性比率'], vert=True, patch_artist=True,
                   boxprops=dict(facecolor='#BBDEFB', color='#1565C0'),
                   medianprops=dict(color='#C62828', linewidth=2),
                   flierprops=dict(marker='o', markerfacecolor='#FF8F00', markersize=4))
ax3a.axhline(upper_iqr, color='#C62828', linestyle='--', linewidth=1.2, alpha=0.7,
             label=f'IQR上限={upper_iqr:.1f}%')
ax3a.axhline(lower_iqr, color='#C62828', linestyle='--', linewidth=1.2, alpha=0.7,
             label=f'IQR下限={lower_iqr:.1f}%')
ax3a.set_ylabel('若年女性比率（%）', fontsize=11)
ax3a.set_title('箱ひげ図（IQR×1.5外れ値）', fontsize=11, fontweight='bold')
ax3a.legend(fontsize=9)
ax3a.grid(axis='y', alpha=0.3)

if has_kawashima:
    kaw_y = y_vals[kaw_idx]
    ax3a.scatter([1], [kaw_y], color='#2E7D32', s=80, zorder=5)
    ax3a.text(1.15, kaw_y, '←川島町', va='center', fontsize=9, color='#2E7D32')

ax3b = axes3[1]
ax3b.hist(df_main['若年女性比率'], bins=35, color='#1565C0', alpha=0.7, edgecolor='white')
mean_y = df_main['若年女性比率'].mean()
ax3b.axvline(mean_y, color='#C62828', linestyle='--',
             linewidth=1.5, label=f'平均={mean_y:.2f}%')
ax3b.set_xlabel('若年女性比率（%）', fontsize=11)
ax3b.set_ylabel('自治体数', fontsize=11)
ax3b.set_title('ヒストグラム', fontsize=11, fontweight='bold')
ax3b.legend(fontsize=9)
ax3b.grid(axis='y', alpha=0.3)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_H5_3_fig3_box.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2025_H5_3_fig3_box.png")
▼ 実行結果
図3保存: 2025_H5_3_fig3_box.png
💡 解説
  • fig, ax = plt.subplots(...) — 図全体（fig）と軸（ax）を作る定番。以降は ax.bar(...) 等で操作。
  • ax.axhline / ax.axvline — 水平／垂直の点線。平均線や基準線として定番。
💡 Python TIPS s[:-n]「末尾n文字を除く」／s[n:]「先頭n文字を除く」。スライス [start:stop:step] はリスト・タプル・文字列共通の基本ワザです。
6. 外れ値
4
外れ値の箱ひげ図診断

521自治体の若年女性人口増加率には、急激な人口増加・減少を示す外れ値が含まれる可能性がある。IQR×1.5法で外れ値を診断する。

若年女性人口増加率の分布
図1：521自治体の若年女性人口増加率の分布。箱ひげ図（左）とヒストグラム（右）。
📌 この箱ひげ図の読み方
このグラフは
データの分布（中央値・四分位範囲・外れ値）を箱と線で表したグラフ。
読み方
箱の中の線が中央値。箱の上下端が25%・75%点（四分位範囲）。箱の外の点が外れ値。
なぜそう解釈できるか
箱が高い位置にあるほど値が大きいグループ。箱の大きさがばらつきの大きさ。グループ間で箱が重なっていなければ有意差の証拠になりやすい。
外れ値の解釈 人口増加率が極端に高い自治体は、新興住宅地開発や工場立地など特殊要因による「一時的な人口流入」の可能性がある。極端に低い自治体は「既に消滅段階」にある場合もある。外れ値を機械的に除外するのではなく、その背景を確認することが重要。

DS LEARNING POINT 3

IQR法とその限界

IQR（四分位範囲）法は正規分布を仮定しない頑健な外れ値検出法だが、自治体データのように右裾が長い分布では多くの自治体が「外れ値」として検出される可能性がある。

Q1 = np.percentile(y, 25) Q3 = np.percentile(y, 75) IQR = Q3 - Q1 # 通常の箱ひげ図基準（×1.5） lower_15 = Q1 - 1.5 * IQR upper_15 = Q3 + 1.5 * IQR n_out_15 = ((y < lower_15) | (y > upper_15)).sum() # 論文では×1.5を検討 # 外れ値が多すぎる場合は ×2.0 や×3.0 を検討 print(f"IQR×1.5: {n_out_15}件の外れ値") print(f"外れ値除外後 N = {N - n_out_15}") # 外れ値の確認：どの自治体か？ outlier_mask = (y < lower_15) | (y > upper_15) print("外れ値自治体の特徴:", df[outlier_mask][['地域', 'y']].head())
やってみよう共通設定
📝 コード
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
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'
os.makedirs(FIG_DIR, exist_ok=True)
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • import pandas as pd など — 必要なライブラリをまとめて呼び出します。as pd は短い別名（alias）。
  • matplotlib.use('Agg') — グラフを画面表示せずファイルに保存するためのおまじない。
  • plt.rcParams['font.family'] — グラフの日本語表示用フォント指定（Macは Hiragino Sans、Windowsなら Yu Gothic 等）。
  • os.makedirs('html/figures', exist_ok=True) — 図の保存先フォルダを作る（既にあってもOK）。
💡 Python TIPS f"...{x}..." はf-string。文字列の中に {変数} と書くだけで埋め込めて、{x:.2f} のように書式も指定できます。
やってみよう■ Step1. 無相関検定
📝 コード
203
204
205
206
207
208
209
210
211
print("\n" + "=" * 60)
print("■ Step1. 無相関検定（ピアソン相関）")
print("=" * 60)
print(f"\n  {'変数':<18} {'r':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 46)
for name, col in zip(VAR_NAMES, X_raw.T):
    r, p = stats.pearsonr(col, y_vals)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {name:<18} {r:>8.4f} {p:>10.4f} {sig:>6}")
▼ 実行結果
============================================================
■ Step1. 無相関検定（ピアソン相関）
============================================================

  変数                        r         p値     有意
  ----------------------------------------------
  保育所数                -0.3134     0.0000    ***
  幼稚園数                 0.0214     0.6256   n.s.
  小学校数                -0.3552     0.0000    ***
  中学校数                -0.3681     0.0000    ***
  製造業従業者割合             0.3253     0.0000    ***
  農業従業者割合             -0.3035     0.0000    ***
  女性就業率               -0.6485     0.0000    ***
  一般病院数               -0.3384     0.0000    ***
  一般診療所数              -0.3153     0.0000    ***
  農家比率                -0.2870     0.0000    ***
💡 解説
  • stats.pearsonr(x, y) — Pearson相関係数 r と p値を同時に返します。
💡 Python TIPS Seriesの .map() は「1対1の置き換え」、.apply() は「関数を当てる」。辞書なら .map()、ロジックなら .apply()。
やってみよう■ Step2. バックワード変数選択（p<0.05）
📝 コード
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
print("\n" + "=" * 60)
print("■ Step2. バックワード変数選択（閾値 p < 0.05）")
print("=" * 60)

current_vars = list(range(len(VAR_NAMES)))
step = 0
while True:
    X_cur = sm.add_constant(X_raw[:, current_vars])
    model = sm.OLS(y_vals, X_cur).fit()
    pvals = np.asarray(model.pvalues[1:])  # 定数項を除く（numpy配列に変換）
    max_p_idx = int(np.argmax(pvals))
    max_p = pvals[max_p_idx]
    if max_p < 0.05:
        break
    removed_var_idx = current_vars[max_p_idx]
    step += 1
    print(f"  Step{step}: {VAR_NAMES[removed_var_idx]} を除外（p={max_p:.4f}）")
    current_vars.pop(max_p_idx)

selected_names = [VAR_NAMES[i] for i in current_vars]
final_model = sm.OLS(y_vals, sm.add_constant(X_raw[:, current_vars])).fit()
print(f"\n  最終選択変数: {selected_names}")
print(f"  R² = {final_model.rsquared:.4f}")
print(final_model.summary2())
▼ 実行結果
============================================================
■ Step2. バックワード変数選択（閾値 p < 0.05）
============================================================
  Step1: 幼稚園数 を除外（p=0.8944）
  Step2: 小学校数 を除外（p=0.8019）
  Step3: 農家比率 を除外（p=0.4418）
  Step4: 保育所数 を除外（p=0.2521）

  最終選択変数: ['中学校数', '製造業従業者割合', '農業従業者割合', '女性就業率', '一般病院数', '一般診療所数']
  R² = 0.5694
                 Results: Ordinary least squares
==================================================================
Model:              OLS              Adj. R-squared:     0.564    
Dependent Variable: y                AIC:                2720.1934
Date:               2026-05-18 11:24 BIC:                2750.0104
No. Observations:   523              Log-Likelihood:     -1353.1  
Df Model:           6                F-statistic:        113.7    
Df Residuals:       516              Prob (F-statistic): 4.38e-91 
R-squared:          0.569            Scale:              10.485   
--------------------------------------------------------------------
           Coef.    Std.Err.      t       P>|t|     [0.025    0.975]
--------------------------------------------------------------------
const     77.5655     1.5353    50.5229   0.0000   74.5494   80.5816
x1        -0.6234     0.1326    -4.7018   0.0000   -0.8839   -0.3629
x2         0.0376     0.0079     4.7723   0.0000    0.0221    0.0530
x3         0.0887     0.0218     4.0617   0.0001    0.0458    0.1316
x4        -0.3019     0.0190   -15.9170   0.0000   -0.3392   -0.2647
x5        -1.0513     0.1812    -5.8010   0.0000   -1.4074   -0.6953
x6        -0.3320     0.0539    -6.1560   0.0000   -0.4380   -0.2261
------------------------------------------------------------------
Omnibus:              6.322         Durbin-Watson:           1.485
Prob(Omnibus):        0.042         Jarque-Bera (JB):        6.150
Skew:                 -0.244        Prob(JB):                0.046
Kurtosis:             3.208         Condition No.:           1020 
==================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the
errors is correctly specified.
[2] The condition number is large, 1.02e+03. This might indicate
that there are strong multicollinearity or other numerical
problems.
💡 解説
  • sm.add_constant(X) — 切片項（定数1の列）を先頭に追加。statsmodelsで必須。
  • sm.OLS(y, X).fit() — 最小二乗法でモデルを推定。model.params, model.pvalues, model.conf_int() で結果取得。
💡 Python TIPS [式 for x in リスト] はリスト内包表記。forループでappendする代わりに1行でリストを作れます。
やってみよう■ Step3. 外れ値処理後の再推定
📝 コード
236
237
238
239
240
241
242
243
print("\n" + "=" * 60)
print("■ Step3. 外れ値除外（IQR×1.5）後の再推定")
print("=" * 60)

X_clean = X_raw[mask_iqr.values][:, current_vars]
y_clean = y_vals[mask_iqr.values]
clean_model = sm.OLS(y_clean, sm.add_constant(X_clean)).fit()
print(f"  除外後 N={mask_iqr.sum()}, R² = {clean_model.rsquared:.4f}")
▼ 実行結果
============================================================
■ Step3. 外れ値除外（IQR×1.5）後の再推定
============================================================
  除外後 N=522, R² = 0.5677
💡 解説
  • sm.add_constant(X) — 切片項（定数1の列）を先頭に追加。statsmodelsで必須。
  • sm.OLS(y, X).fit() — 最小二乗法でモデルを推定。model.params, model.pvalues, model.conf_int() で結果取得。
💡 Python TIPS r, p = stats.pearsonr(...) — Pythonは複数戻り値を同時に受け取れる（タプルアンパック）。
やってみよう図図4 (fig4_case): 事例分析 — 保育所数・製造業割合 vs 若年女性比率
📝 コード
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('保育・雇用環境と若年女性比率の関係（実データ）', fontsize=13, fontweight='bold')

# (a) 保育所数 vs 若年女性比率
ax4a = axes4[0]
x_nursery = df_main['保育所数'].values
ax4a.scatter(x_nursery, y_vals, color='#1565C0', s=20, alpha=0.4, label='各自治体', zorder=2)
if has_kawashima:
    ax4a.scatter([x_nursery[kaw_idx]], [y_vals[kaw_idx]],
                 color='#C62828', s=150, zorder=5, label='川島町', marker='*')
z4 = np.polyfit(x_nursery, y_vals, 1)
xs4 = np.linspace(x_nursery.min(), x_nursery.max(), 100)
ax4a.plot(xs4, np.poly1d(z4)(xs4), 'r-', linewidth=1.5, alpha=0.7)
r4, p4 = stats.pearsonr(x_nursery, y_vals)
ax4a.set_xlabel('保育所数（人口1万人あたり）', fontsize=11)
ax4a.set_ylabel('若年女性比率（%）', fontsize=11)
ax4a.set_title(f'保育所数との関係 r={r4:.3f} (p={p4:.4f})', fontsize=11, fontweight='bold')
ax4a.legend(fontsize=9)
ax4a.grid(True, alpha=0.3)
▼ 実行結果
このステップは print はしません。データや図が裏で更新されただけ。次のステップへ進みましょう。
💡 解説
  • fig, ax = plt.subplots(...) — 図全体（fig）と軸（ax）を作る定番。以降は ax.bar(...) 等で操作。
  • stats.pearsonr(x, y) — Pearson相関係数 r と p値を同時に返します。
💡 Python TIPS np.cumsum(arr) は累積和、np.linspace(a, b, n) は「aからbを等間隔でn個」。NumPyの定石です。
やってみよう図図4 (fig4_case): 事例分析 — 保育所数・製造業割合 vs 若年女性比率 — (b) 製造業従業者割合 vs 若年女性比率
📝 コード
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
# (b) 製造業従業者割合 vs 若年女性比率
ax4b = axes4[1]
x_mfg = df_main['製造業従業者割合'].values
ax4b.scatter(x_mfg, y_vals, color='#1565C0', s=20, alpha=0.4, label='各自治体', zorder=2)
if has_kawashima:
    ax4b.scatter([x_mfg[kaw_idx]], [y_vals[kaw_idx]],
                 color='#C62828', s=150, zorder=5, label='川島町', marker='*')
z4b = np.polyfit(x_mfg, y_vals, 1)
xs4b = np.linspace(x_mfg.min(), x_mfg.max(), 100)
ax4b.plot(xs4b, np.poly1d(z4b)(xs4b), 'r-', linewidth=1.5, alpha=0.7)
r4b, p4b = stats.pearsonr(x_mfg, y_vals)
ax4b.set_xlabel('製造業従業者割合（%）', fontsize=11)
ax4b.set_ylabel('若年女性比率（%）', fontsize=11)
ax4b.set_title(f'製造業従業者割合との関係 r={r4b:.3f} (p={p4b:.4f})',
               fontsize=11, fontweight='bold')
ax4b.legend(fontsize=9)
ax4b.grid(True, alpha=0.3)

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

print("\n全図の生成完了（4枚）")
print(f"  fig1_corr.png  : 無相関検定（相関係数棒グラフ）")
print(f"  fig2_backward.png: バックワード変数選択後の回帰係数")
print(f"  fig3_box.png   : 箱ひげ図・ヒストグラム")
print(f"  fig4_case.png  : 事例分析（川島町・散布図）")
▼ 実行結果
図4保存: 2025_H5_3_fig4_case.png

全図の生成完了（4枚）
  fig1_corr.png  : 無相関検定（相関係数棒グラフ）
  fig2_backward.png: バックワード変数選択後の回帰係数
  fig3_box.png   : 箱ひげ図・ヒストグラム
  fig4_case.png  : 事例分析（川島町・散布図）
💡 解説
  • stats.pearsonr(x, y) — Pearson相関係数 r と p値を同時に返します。
💡 Python TIPS f-stringの書式 {値:.2f}（小数2桁）、{値:,}（3桁区切り）、{値:>10}（右寄せ10桁）など、覚えると出力が一気に整います。

まとめ

主要な発見

SSDSE-A の521自治体（人口0.5〜2万人）のデータを用いたバックワード重回帰分析の結果：

  1. 保育環境の充実（正）：保育所数が多い自治体ほど若年女性が定着しやすい。子育て支援は消滅可能性回避の最重要施策。
  2. 製造業雇用（正）：安定した雇用が若年女性の経済的自立を支え、定住を促す。
  3. 貧困指標（負）：生活保護率が高い自治体では若年女性が流出する傾向。
  4. 川島町事例：統計モデルが示す「保育・製造業雇用の充実」を体現した成功事例。
政策への示唆 「消滅可能性自治体を救う」ためには、子育て支援（保育所の整備）と雇用創出（製造業の誘致・維持）の組み合わせが統計的に支持される。
教育的価値（この分析から学べること）
  • 消滅可能性自治体：若年女性人口の減少率で定義される指標。地方創生の重要キーワード。
  • 人口減少の構造：出生率・転出・高齢化が複合的に作用する。単一原因では説明できない。
  • 政策の効果検証：移住促進・子育て支援などの政策効果を地域比較で評価する考え方。

データ・コードのダウンロード

分析スクリプト（2025_H5_3_shorei.py）
データ出典
SSDSE-A 市区町村データ統計数理研究所 SSDSE（社会・人口統計体系）
消滅可能性都市リスト日本創成会議（増田レポート、2014年）
市区町村別保育所数厚生労働省 保育所等関連状況取りまとめ

本教育用コードは合成データを使用（np.random.seed(2027)）。実際の分析はSSDSE-Aの実データによる。

教育用再現コード ｜ 2025年 統計データ分析コンペティション 審査員奨励賞 [高校生の部] ｜ 早稲田大学系属 早稲田実業学校

⚠️ よくある誤解と注意点

統計分析の解釈で初心者がやりがちな勘違いをまとめます。特に「相関と因果の混同」「p値の過信」は研究現場でもよく起きる落とし穴です。本文を読む前にも、読んだ後にも、目を通してみてください。

❌ 「相関がある＝因果関係がある」ではない
疑似相関（spurious correlation）とは、見かけ上は関係があるように見えるが、実際は無関係、または第三の変数（交絡変数）が両方に影響しているだけの現象です。

古典例： アイスクリームの売上 と 水難事故件数 は強く相関するが、片方が他方を引き起こしているわけではない。両者とも「夏の暑さ」という第三の変数に引きずられているだけ。

論文を読むときの心構え： 「○○と△△に強い相関が見られた」だけで終わっている主張は、本当に因果関係があるのか、それとも第三の変数（人口・所得・地理など）が共通要因として効いているだけではないかを必ず疑ってください。
❌ 「p値が小さい＝重要な発見」ではない
p値が小さい（例えば p < 0.001）ことは「統計的に偶然とは考えにくい」という意味であって、「実用的に大きな効果がある」という意味ではありません。

例： 巨大なサンプルサイズ（n=100,000）では、相関係数 r=0.02 でも p < 0.001 になります。しかし r=0.02 は実用上ほぼ無視できる関係です。

正しい読み方： p値と効果量（係数の大きさ、相関係数の値）の両方をセットで判断してください。p値だけで「重要な発見」と結論づけるのは誤りです。
❌ 「回帰係数が大きい＝重要な変数」ではない
回帰係数の絶対値は、説明変数の単位に強く依存します。「年収（万円）」と「失業率（%）」の係数を直接比較しても意味がありません。

正しい比較方法： (1) 標準化係数（各変数を平均0・分散1に変換した上での係数）を使う、(2) 限界効果（変数を1標準偏差動かしたときのyの変化）で比較する。

また、係数の大きさが「因果関係の強さ」を意味するわけでもありません。あくまで「相関的な関連の強さ」です。
❌ 「外れ値を除外すれば正しい結果」ではない
外れ値（極端な値）を「目障りだから」「結果が綺麗にならないから」という理由で除外するのは分析の改ざんに近い行為です。

外れ値が示すもの： 本当に重要な情報（東京の超高密度、北海道の超低密度など）であることが多い。外れ値を取り除くと「日本全体の傾向」を見誤る原因になります。

正しい対処： (1) 外れ値の出現要因を調査する（なぜ東京だけ突出するのか）、(2) ノンパラメトリック手法（Spearman相関・Kruskal-Wallis）を使う、(3) 外れ値を含む結果と除外した結果の両方を提示し、解釈を読者に委ねる。
❌ 「サンプルサイズが大きい＝信頼できる」ではない
サンプルサイズ（n）が大きいと統計的検定の検出力は上がりますが、それは「偶然による誤差を減らす効果」にすぎません。

nが大きくても解消されない問題：
・選択バイアス（標本が偏っている）
・測定誤差（変数の定義が曖昧）
・欠損値のパターン（欠損がランダムでない）
・交絡変数の見落とし

例： 1万人にWeb調査して「ネット利用と幸福度は強く相関」と言っても、そもそも回答者がネットユーザー寄りに偏っているため、母集団全体の結論にはなりません。
❌ 「複雑なモデル＝より良い分析」ではない
ランダムフォレスト・ニューラルネット・複雑な階層モデルなど、高度な手法を使えば「良い分析」と感じがちですが、必ずしもそうではありません。

過学習（overfitting）の罠： モデルが複雑すぎると、訓練データの偶然のパターンまで学習してしまい、新しいデータでは予測精度が落ちます。

シンプルさの価値： 重回帰分析や相関分析は「結果が解釈しやすい」「再現性が高い」という大きな利点があります。複雑な手法はシンプルな手法で答えが出ない時の最後の手段です。
❌ 「多重共線性は気にしなくていい」ではない
多重共線性とは、説明変数同士の相関が極めて強い状態のこと。これを放置すると、回帰係数の符号や大きさが入れ替わる異常事態が起こります。

典型例： 「総人口」と「労働力人口」を同時に投入すると、両者の相関が r=0.99 になり、係数推定が極端に不安定になります。「総人口は正だが、労働力人口は負」のような解釈不能な結果になりがちです。

診断と対処：
・VIF（分散拡大係数）を計算し、VIF > 10 の変数を確認
・相関行列で |r| > 0.8 のペアをチェック
・対処法：一方を除外、合成変数（PCA）に変換、Ridge回帰で安定化
❌ 「R²が高い＝良いモデル」ではない
決定係数 R² はモデルの「当てはまりの良さ」を示しますが、R² が高くてもモデルが正しいとは限りません。

R² が高くなる罠：
・説明変数を増やせば R² は自動的に上がる（無関係な変数を追加してもR²は下がらない）
・時系列データでは、共通のトレンド（時間とともに増加）があるだけで R² が 0.9 を超える
・サンプルサイズが小さいとR²が過大評価される

代替指標： 調整済み R²（変数の数でペナルティ）、AIC・BIC（モデル選択基準）を併用してください。予測力の真の評価には交差検証（cross-validation）でテストデータの R² を見ること。
❌ 「ステップワイズで選んだ変数は重要」ではない
ステップワイズ法（バックワード・フォワード選択）は便利ですが、p値ベースの変数選択は再現性に問題があると批判されています。

問題点：
・同じデータでも実行順序によって最終モデルが変わる
・p値を繰り返し見ることで「偶然に有意な変数」を拾ってしまう（p-hacking）
・係数の標準誤差が過小評価され、信頼区間が嘘っぽくなる

より良い方法：
・事前に変数を理論で絞る（先行研究から候補を選ぶ）
・LASSO回帰（自動かつ統計的に正当化された変数選択）を使う
・交差検証で AIC/BIC 最小モデルを選ぶ
❌ 「線形回帰なら線形関係を前提にすべき」
重回帰分析は線形関係を前提とします。実際の関係が非線形なのに線形モデルで分析すると、本当の関係を見逃します。

非線形の例：
・U字型関係： 失業率と物価上昇率（フィリップス曲線）
・逓減効果： 所得と幸福度（年収 800万円までは強い正の効果、それ以上は飽和）
・閾値効果： 高齢化率と医療費（ある水準を超えると急激に上がる）

診断と対処：
・残差プロットで残差が0周辺に均等に分布しているか確認
・変数の対数変換・二乗項追加で非線形性を取り込む
・どうしても線形では捉えられないなら、機械学習（RF・GBM）を併用する
❌ 「データに当てはまった＝予測に使える」ではない
「過去のデータでフィットしたから将来も予測できる」と思うのは危険です。

過学習（overfitting）の例： 47都道府県のデータに10個の説明変数を投入すれば、ほぼ完璧にフィットします（自由度がほぼゼロ）。でもそのモデルを新しい年度に適用すると、予測精度はほぼランダム並みに落ちることがあります。

正しい予測力の評価：
・データを訓練用 70%とテスト用 30%に分割し、テスト用での予測精度を見る
・k分割交差検証（k-fold CV）で予測の安定性を確認
・「説明変数の数 ≪ サンプルサイズ」のバランスを意識（目安：n > 10 × 変数数）

📖 用語集（この記事に出てくる統計用語）

統計の基本用語を初心者向けに解説します。本文中で見慣れない言葉が出てきたら、ここに戻って確認してください。

p値
「効果がない」と仮定したときに、観察されたデータ（またはより極端なデータ）が得られる確率。0〜1の値で、慣例的に 0.05（5%）未満を「有意」と判断する。
有意水準
「偶然」と「意味のある違い」を分ける基準。通常 α=0.05（5%）を使う。p値 < α なら「有意」と判定。
信頼区間
「真の値はこの範囲にあるだろう」という幅。95%信頼区間 = 同じ実験を100回繰り返したら95回はこの範囲に真の値が入る。
サンプルサイズ
分析に使ったデータ点の数（n）。一般にnが大きいほど推定が安定し、わずかな差も検出できるようになる。
標準誤差
推定値（係数など）のばらつきの目安。標準誤差が小さいほど推定値が安定している。
正規分布
釣鐘型の左右対称な分布。多くのパラメトリック検定（t検定・F検定など）は「データが正規分布に従う」ことを仮定する。
因果と相関
「相関がある」と「原因と結果の関係（因果）」は別物。アイスクリームの売上と水難事故は相関するが、原因は両者とも「夏の暑さ」。
外れ値
他のデータから極端に離れた値。分析結果を歪める原因になるため、検出して除外するか別途扱う必要がある。
欠損値
データが取得できなかった部分（NaN・空白）。除外するか補完（平均代入・回帰代入など）するかが分析上の重要な判断点。
VIF
Variance Inflation Factor（分散拡大係数）。多重共線性の強さを示す指標。VIF > 10 で「強い多重共線性あり」と判断。
係数（回帰係数）
「説明変数 x が1単位増えたとき、目的変数 y が平均でどれだけ変化するか」を示す数値。正の値は正の影響、負の値は負の影響。
多重共線性
説明変数同士の相関が強すぎる状態。係数推定が不安定になり、解釈を誤る原因になる。VIF > 10 が警告サイン。
標準化係数
変数の単位の影響を取り除いた係数。複数の変数の影響の大きさを単位に依存せず比較するために使う。
決定係数 R²
回帰モデルが目的変数のばらつきの何%を説明できるかを示す指標。0〜1の値で、1に近いほどモデルの説明力が高い。

📐 使っている手法をわかりやすく解説

統計手法について「何のためか」「結果をどう読むか」を初心者向けに解説します。

◆ 統計の基本概念（どの論文にも共通）

🔍 p値（有意確率）とは
何？
「もし本当に効果がなかったとしたら、今回の結果（またはもっと極端な結果）が偶然起きる確率」のこと。
なぜ必要？
帰無仮説（「効果なし」の仮定）のもとで検定統計量の分布から計算する。
何がわかる？
「この関係は偶然ではなく、統計的に意味がある」と主張するための客観的な根拠になる。
読み方
p < 0.05（5%未満）を「統計的に有意」と判断するのが慣例。ただし「p値が小さい＝効果が大きい」ではない。効果量（係数の大きさ）とセットで判断する。
🗂️ ノンパラメトリック検定とは（なぜ使うのか）
何？
「データが正規分布に従う」という仮定を置かない検定手法の総称。Kruskal-Wallis検定・Mann-Whitney U検定などが代表例。
なぜ必要？
データの値ではなく「順位」に変換して検定統計量を計算する。外れ値や偏った分布に対しても安定して機能する。
何がわかる？
サンプルサイズが小さい・データが歪んでいる・外れ値がある場合でも、グループ差の有無を検定できる。
読み方
「なぜノンパラメトリックを選ぶのか」の理由を示すには、正規性検定（Shapiro-Wilk）の結果を添えるのが望ましい。結果の解釈は対応するパラメトリック検定と同様（p < 0.05 で有意差あり）。

◆ この論文で使われている手法

📈 重回帰分析
何？
複数の説明変数（原因候補）が1つの目的変数（結果）にどれだけ影響するかを同時に推定する手法。
どう使う？
目的変数 y を複数の説明変数 x₁, x₂, … で予測する式（y = a₁x₁ + a₂x₂ + … + b）を最小二乗法でフィットさせる。
何がわかる？
複数の要因が混在するなかで「どれが一番効いているか」を一度に検証できる。交絡変数を統制できる。
結果の読み方
係数（a₁, a₂…）のプラスは正の影響、マイナスは負の影響。p < 0.05 で統計的に有意。R²が1に近いほどモデルの説明力が高い。
⚠️ 注意点
(1) 多重共線性を必ずVIFで確認（VIF>10で警告）。(2) 線形性の仮定—関係が曲線なら対数変換や二乗項を追加。(3) 残差プロットで正規性・等分散性を確認。(4) サンプル数は最低でも「説明変数数×10」が目安。(5) 外れ値1つで係数が大きく変わるのでCook距離で確認。
🔗 相関分析
何？
2つの変数の「一緒に増減する傾向の強さと向き」を −1〜+1 の相関係数 r で数値化する手法。
どう使う？
散布図を描き、Pearson（連続データ）または Spearman（順序データ・外れ値に強い）の相関係数を計算する。
何がわかる？
「気温が高い県ほど熱中症指標が高い」などの傾向を素早く確認できる。変数選択の第一歩として使われることも多い。
結果の読み方
r > +0.7 は強い正の相関、r < −0.7 は強い負の相関、|r| < 0.3 はほぼ無相関。相関は因果関係を示すものではない点に注意。
⚠️ 注意点
(1) 多重共線性を必ずVIFで確認（VIF>10で警告）。(2) 線形性の仮定—関係が曲線なら対数変換や二乗項を追加。(3) 残差プロットで正規性・等分散性を確認。(4) サンプル数は最低でも「説明変数数×10」が目安。(5) 外れ値1つで係数が大きく変わるのでCook距離で確認。
🌿 Ward法クラスタリング
何？
データをグループ（クラスター）に自動分類する手法。グループ内のばらつきが最小になるよう統合していく。
どう使う？
統合後の「ばらつき増加」が最小になるペアを繰り返し合体させ、デンドログラム（樹形図）で可視化する。
何がわかる？
都道府県を「都市型」「農村型」などのグループに自動分類し、グループ間の特徴比較ができる。
結果の読み方
デンドログラムの切り位置でクラスター数を決める。各クラスターの変数平均を見てグループを命名・解釈する。
⚠️ 注意点
(1) 多重共線性を必ずVIFで確認（VIF>10で警告）。(2) 線形性の仮定—関係が曲線なら対数変換や二乗項を追加。(3) 残差プロットで正規性・等分散性を確認。(4) サンプル数は最低でも「説明変数数×10」が目安。(5) 外れ値1つで係数が大きく変わるのでCook距離で確認。
🎛️ AIC基準によるステップワイズ変数選択
何？
多数の候補変数からモデルの「精度」と「複雑さ」のバランスが最良な変数の組み合わせを自動選択する手法。
どう使う？
バックワード（全変数から除去）またはフォワード（空から追加）で、AIC最小を目指して変数を探索する。
何がわかる？
「30変数中で最も説明力が高い5変数はどれか」を客観基準で決められる。恣意的な変数選択を回避できる。
結果の読み方
AICは小さいほど良い。最終的に残った変数がモデルに「有効」と判断された変数。
⚠️ 注意点
(1) 多重共線性を必ずVIFで確認（VIF>10で警告）。(2) 線形性の仮定—関係が曲線なら対数変換や二乗項を追加。(3) 残差プロットで正規性・等分散性を確認。(4) サンプル数は最低でも「説明変数数×10」が目安。(5) 外れ値1つで係数が大きく変わるのでCook距離で確認。
↔️ VAR（ベクトル自己回帰）/ Granger因果検定
何？
複数の時系列変数が互いに影響し合う関係を分析する手法（VAR）と、「AがBの予測に役立つか」を検定する手法（Granger因果）。
どう使う？
VARは全変数を互いに説明変数として同時回帰。Granger因果はF検定でAのラグ変数がBの予測精度を向上させるかを確認する。
何がわかる？
「女性就業率と出生率はどちらが先に動くか」「リード・ラグ関係」を特定できる。
結果の読み方
Granger因果 p < 0.05 → 「Aの過去値はBの予測に役立つ」（ただし真の因果とは限らない）。
⚠️ 注意点
(1) 多重共線性を必ずVIFで確認（VIF>10で警告）。(2) 線形性の仮定—関係が曲線なら対数変換や二乗項を追加。(3) 残差プロットで正規性・等分散性を確認。(4) サンプル数は最低でも「説明変数数×10」が目安。(5) 外れ値1つで係数が大きく変わるのでCook距離で確認。

🚀 発展の可能性（結果 X → 新仮説 Y → 課題 Z）

この研究をさらに発展させるための3つの方向性を示します。「今回わかったこと（X）」から「次に検証すべき仮説（Y）」を立て、「具体的に何をするか（Z）」まで考えてみましょう。

① データ・時間的拡張
結果 X
本論文は特定の年度・地域の断面データ（または限られた時系列）で分析を行った。
新仮説 Y
より新しい年度のデータや市区町村レベルの細粒度データを使えば、知見の時間的頑健性や地域内格差を検証できる。
課題 Z
（1）統計センターから最新の SSDSE をダウンロードし、同じ分析を再実行する。（2）結果が変わった場合、その要因（コロナ・政策変化など）を考察する。（3）市区町村データ（SSDSE-A/C/F）で分析単位を細かくした場合の結果と比較する。
② 手法の発展：重回帰分析 の次のステップ
結果 X
本論文は 重回帰分析 を用いた推定を行った。
新仮説 Y
パネルデータ固定効果モデル（FE）による都道府県固有の差の統制 により、本分析では統制できていない問題を解消できる可能性がある。
課題 Z
（1）パネルデータ固定効果モデル（FE）による都道府県固有の差の統制 を実装し、本論文の係数推定と比較する。（2）操作変数法（IV）による内生性の解消 も試し、結果の頑健性を確認する。（3）推定結果の変化から、元の分析の仮定のどれが重要だったかを考察する。
③ 政策提言・実践への応用
結果 X
本論文は分析結果から特定の変数が目的変数に影響することを示した。
新仮説 Y
分析対象を日本全国から特定地域に絞ること、または逆に国際比較に拡張することで、政策の移転可能性と文脈依存性を検証できる。
課題 Z
（1）有意な変数を「政策で変えられるもの」と「変えにくいもの」に分類する。（2）政策で変えられる変数について、係数の大きさから「どれだけ変えればどれだけ効果があるか」を試算する。（3）自治体・政策立案者への提言として、実現可能なアクションプランを1枚にまとめる。

🎯 自分でやってみよう（5つのチャレンジ）

学んだだけでは身につきません。実際に手を動かすのが最強の学習方法です。本論文のスクリプトをベースに、以下のチャレンジに挑戦してみてください。難易度別に5つ用意しました。

★☆☆☆☆ 入門
CH1. 同じデータで分析を再現する
まずは付属の Python スクリプトをそのまま実行し、論文と同じ図を再現してみてください。
ポイント： 各図がどのコード行から生成されているか辿る。エラーが出たら原因を考える。
★★☆☆☆ 初級
CH2. 説明変数を1つ追加・除外して結果を比較
本論文の分析モデルから説明変数を1つ抜いて再実行、あるいは1つ追加して再実行してください。
ポイント： 係数・p値・R² がどう変わったか観察する。多重共線性が原因で結果が変わる例を見つけられたら理想的。
★★★☆☆ 中級
CH3. 別の年度・別の都道府県で同じ分析を試す
SSDSE の別の年度（例：2015年度・2020年度）または特定都道府県のみのデータで同じ分析を実行してください。
ポイント： 時代や地域によって結論が変わるか？ 変わるならその理由を考察する。
★★★★☆ 上級
CH4. 別の手法を組み合わせる
本論文の手法 + 1つの追加手法（例：重回帰 + LASSO、相関分析 + 主成分分析）で結果を比較してください。
ポイント： 手法の違いで結論が変わるか？ どちらが妥当かを「なぜ」とともに説明できるように。
★★★★★ 発展
CH5. オリジナルの問いを立てて分析する
本論文の手法を借りて、あなた自身の問いを立てて分析してください。 例：「カフェの数と幸福度に関連はあるか」「教育費の高い県は出生率も高いか」など。
ポイント： 問い・データ・手法・結論を1ページのレポートにまとめる。これがデータサイエンスの「実践」。
💡 ヒント： 詰まったら本サイトの他の論文（同じ手法を使っている）のスクリプトをコピーして組み合わせるのが効率的です。手法ガイド・用語集も参考に。

💼 この手法は実社会でこう使われている

本論文で学んだ手法は、研究の世界だけでなく、行政・企業・NPO の現場でも様々に活用されています。具体的なシーンを紹介します。

🏛️
行政の政策立案
都道府県・市区町村の政策担当者は、本論文と同様のデータ分析を用いて「どこに予算を投じれば効果が出るか」を検討します。 例えば医療費削減策、移住促進策、子育て支援策などの効果予測・効果検証に直結します。
🏢
企業のマーケティング・出店戦略
小売チェーン・サービス業の出店戦略では、地域特性（人口構成、所得、ライフスタイル）と売上の関係を本論文と同じ手法で分析します。 ECサイトでも顧客セグメント分析・購買要因分析に類似手法が使われます。
🏥
医療・公衆衛生
感染症の流行予測、医療資源配分の最適化、健康格差の地域要因分析などで、本論文の統計手法は標準的に使われています。 WHO・厚労省レベルの政策評価でも同じ手法が活躍しています。
📊
メディア・ジャーナリズム
新聞・テレビの社会調査記事、選挙予測、世論調査の分析でも、本論文と同じ手法（回帰分析・クラスタリングなど）が使われています。 データジャーナリズムの記事はこの種の分析が中核です。
🎓
学術研究（隣接分野）
経済学・社会学・公衆衛生学・教育学・地理学などの実証研究では、本論文と同じ手法が日常的に使われます。 専門誌に掲載される論文の8割以上が、こうした統計手法に基づいて結論を出しています。
💰
金融・保険業界
与信判断（融資審査）、保険料の地域別設定、不動産価格予測などで、本論文と同様のモデリング手法が広く活用されています。 統計分析の能力は金融業界の必須スキルになっています。

🤔 よくある質問（読者からの想定Q&A）

この論文を読んで初心者が抱きやすい疑問に、教育的観点から答えます。

Q1. この分析、自分でもできますか？
はい、できます。SSDSE データは無料で公開されており、Python の pandas, scikit-learn, statsmodels を使えば全く同じ手順で再現可能です。本ページ下部のスクリプトを実行するだけで結果が得られます。
Q2. 使われている手法は他の分野にも応用できますか？
十分応用可能です。本論文の[手法]は、医療・教育・経済・環境など他のドメインでも標準的に使われる手法です。データの中身（変数）を入れ替えるだけで、別の問いにも適用できます。
Q3. 結論は本当に「因果関係」を示していますか？
本論文は「観察データ」を使った分析であり、厳密な意味での「因果関係」を完全に証明したわけではありません。あくまで「強い関連が見られた」という事実を提示しているにとどまります。真の因果を示すには、無作為化比較試験（RCT）か、自然実験を活用したIV・DiD 等の手法が必要です。
Q4. データの最新版を使うとどうなりますか？
SSDSE は毎年更新されているため、最新版を使えば近年のトレンド（特にコロナ禍以降の変化）も含めて分析できます。ただし、結論が変わる可能性もあります。それ自体が新しい発見につながります。
Q5. もっと深く学ぶには何を読めばいいですか？
「計量経済学」「データサイエンス入門」「統計的因果推論」などのテキストが入門に向いています。Python の場合は『Python ではじめる機械学習』（オライリー）、R の場合は『R で学ぶ統計学』が定番です。本サイトの他の論文も読み比べてみてください。
← 前の論文
自己肯定感を高める要因となる共通特徴の探求
次の論文 →
不登校の原因究明と対策