ビジネスの要因分析に効く回帰分析と統計モデリングの実践ガイド

序章 なぜ集計の先に回帰分析が必要なのか

単純比較が誤解を生む構造

「施策を実施した店舗は、していない店舗より売上が高い。だから施策は有効だ」。この推論の危うさは、多くの方が直感的に気づかれると思います。施策を実施したのが大型店ばかりだったなら、売上の差は施策の効果ではなく店舗規模の差かもしれません。データ分析の言葉で言えば、比較したい要因以外の条件がそろっていない、ということです。

回帰分析の本質的な価値は、この「条件をそろえる」操作を、データの上で近似的に実現できることにあります。店舗規模や立地といった他の要因をモデルに入れることで、「それらが同じだったとしたら、施策の有無で売上はどれだけ違うか」という問いに答えられるようになります。単純集計から重回帰へ進むことは、単なる手法の高度化ではなく、問いの立て方の進化だと考えています。

「説明」と「予測」を区別する

本ガイドで最も重視するのは、説明と予測の区別です。回帰分析には2つの目的があります。ひとつは説明、すなわち「どの要因が、どれだけ効いているのか」を知ること。もうひとつは予測、すなわち「次の値はいくつになるか」を当てることです。同じ手法を使っていても、この2つでは変数の選び方も、モデルの評価方法も、気にすべき落とし穴もまったく異なります。

予測が目的なら、係数の解釈性を多少犠牲にしても精度を追求でき、その作法は機械学習編で扱いました。一方、説明が目的なら、係数そのものが成果物です。係数が信用できなくなる多重共線性は深刻な問題になり、変数選択には統計指標だけでなく因果の見取り図が必要になります。本ガイドは後者、説明のための回帰を軸に据えつつ、両者の使い分けを随所で整理していきます。

AI時代における統計モデリングの位置づけ

機械学習や生成AIが強力になった現在でも、経営の現場で求められる説明責任は変わりません。「なぜこの施策なのか」「どの要因にどれだけ投資すべきか」という問いには、予測精度ではなく、要因と効果の関係を定量的に示すことが求められます。回帰分析と統計モデリングは、この説明責任に応えるための、いまなお中心的な技術です。マーケティング効果の測定（MMM）、価格設定、与信、医療・品質管理まで、係数の解釈が価値を持つ領域で広く使われ続けています。

本ガイドが繰り返し立ち返る軸は3つです。説明と予測を区別すること。モデルの前提を診断すること。そしてデータの生成のされ方に合わせてモデルを選ぶこと。12章を通じてこの3つを繰り返し確認していきます。

Anagraftでは、AIプロジェクトの構想・課題設計から、データ分析・機械学習モデルの開発、AI人材の育成まで一貫したご支援を行っています。ご相談は、以下よりお問い合わせください。

お問い合わせ

全体の地図

本ガイドの構成です。第1部から順に積み上がりますが、各章は参照用に独立して読める形にもしています。

第1章 単回帰から重回帰へ、複数の要因を同時に扱う

回帰分析とは何をする道具か

本ガイドは「要因を説明するための回帰分析」を軸にした実践ガイドです。統計学基礎編にあたる『データから正しく結論を導く統計学の実践ガイド』では記述統計や仮説検定を、機械学習編にあたる『機械学習の仕組みと使いどころがわかる実践ガイド』では予測を主目的とするモデリングを扱いました。本ガイドはその中間に位置し、「何が結果を動かしているのか」を数値で説明するための道具として回帰分析を扱います。

回帰分析とは、結果として観測される1つの数値（目的変数）を、それに影響を与えていると考えられる1つ以上の数値（説明変数）を使って説明・予測する統計的な枠組みです。目的変数は被説明変数、応答変数とも呼ばれ、説明変数は独立変数、予測変数とも呼ばれます。ビジネスの現場でよく出てくる「何が売上を動かしているのか」「どの施策が離職率に影響しているのか」「顧客単価を左右している要因は何か」といった問いは、いずれも目的変数（売上・離職率・顧客単価）を説明変数（広告費・上司との面談頻度・購入頻度など）で説明するという形に翻訳できます。回帰分析は、こうした問いに対して「どの要因が」「どの程度の大きさで」「どのくらいの確からしさで」影響しているのかを、数値として言語化するための最も基本的な統計的道具です。

回帰分析には大きく分けて2つの使い道があります。1つは、未知のデータに対して目的変数の値を当てにいく予測目的の利用、もう1つは、説明変数それぞれが目的変数にどの程度の影響を与えているかを解釈する説明目的の利用です。両者は同じ数式を使いますが、予測目的では当てはまりの良さ（予測精度）が最優先されモデルの内部構造が複雑でも構わないのに対し、説明目的では各係数の解釈可能性と統計的な妥当性が優先され、モデルが単純で読み解きやすいことに価値が置かれます。機械学習編で扱った予測モデルの多くは前者に軸足を置いていましたが、本ガイドは後者、すなわち「係数を読んで意思決定に使う」ための回帰分析に主眼を置きます。

単回帰の復習：最小二乗法・回帰係数・決定係数

重回帰を理解するための前提として、単回帰（説明変数が1つだけの回帰）の要点を振り返ります。詳しい導出や検定の考え方は統計学基礎編で扱っているため、ここでは重回帰を理解するために必要な最小限の要点だけを整理します。

単回帰は、目的変数\( y \)と説明変数\( x \)の関係を、\( y = b_0 + b_1 x + \varepsilon \)という直線の式で表すモデルです。\( b_0 \)は切片（\( x=0 \)のときの\( y \)の予測値）、\( b_1 \)は回帰係数（\( x \)が1単位増えたときに\( y \)がどれだけ変化するかを表す傾き）、\( \varepsilon \)は誤差項です。この\( b_0 \)と\( b_1 \)を、手元のデータに最もよく当てはまる形で決定する方法が最小二乗法です。「よく当てはまる」とは、実際の観測値\( y_i \)とモデルによる予測値\( \hat{y}_i \)の差（残差）を2乗して全データ分足し合わせた残差平方和が、最も小さくなるという意味です。最小二乗法は、この残差平方和を最小にする\( b_0, b_1 \)を数式的に一意に求める方法であり、後述する重回帰でも同じ考え方がそのまま拡張されます。

モデルの当てはまりの良さを表す指標が決定係数\( R^2 \)です。決定係数は、目的変数の全体のばらつき（分散）のうち、モデルによって説明できた割合を表し、切片を含む通常の回帰を当てはめたデータの上では0から1の値を取ります（切片を置かないモデルや、当てはめに使っていない別のデータで計算した場合には負の値になることもあります）。\( R^2 = 0.7 \)であれば、目的変数のばらつきの70%がその説明変数によって説明できている、という読み方になります。回帰係数がその変数の「効果の大きさと向き」を表すのに対し、決定係数はモデル全体の「説明力」を表すという役割の違いを押さえておくと、この先の議論が整理しやすくなります。

単回帰の限界：要因を1つずつ調べることの問題

現実のビジネス現象は、たいてい複数の要因が同時に絡み合って結果を作り出しています。売上は広告費だけでなく、価格設定や季節要因、競合の動きなど複数の変数の影響を受けて決まります。ここで陥りやすい誤りが、「広告費と売上の単回帰」「価格と売上の単回帰」「気温と売上の単回帰」というように、要因を1つずつ取り出して単回帰を繰り返し、それぞれの係数を並べて「効果の大きさ」として報告してしまうことです。

この進め方が問題になるのは、単回帰の係数には、本来は他の変数に帰属するはずの影響が混ざり込んでしまうためです。単回帰は「その1つの説明変数だけで目的変数のばらつきをすべて引き受ける」形でモデルを組んでいるため、実際には別の変数によって動いている部分まで、その変数の効果であるかのように計算されてしまいます。特に、複数の説明変数どうしが互いに相関している場合（広告費を増やす時期にちょうど値下げキャンペーンも重なりやすい、といった実務上よくある状況）、単回帰の係数は大きく歪んだ値になり得ます。この「他の変数の影響が入り込んでしまう」問題は交絡と呼ばれ、第3章で詳しく扱います。

この問題を実際のコードで確認します。以下は、広告費・価格・気温の3つの要因から売上が決まる架空のデータを生成し、まず広告費だけを説明変数とした単回帰を行うコードです。

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

np.random.seed(42)
n = 200

# 広告費(万円)・価格(円)・気温(度)から売上(万円)が決まる架空データ
ad_spend = np.random.uniform(10, 100, n)
price = np.random.normal(1000, 150, n)
temperature = np.random.uniform(5, 35, n)
noise = np.random.normal(0, 8, n)

sales = 30 + 0.5 * ad_spend - 0.02 * price + 1.2 * temperature + noise

df = pd.DataFrame({
    "sales": sales,
    "ad_spend": ad_spend,
    "price": price,
    "temperature": temperature,
})

# 広告費だけを使った単回帰
X_simple = sm.add_constant(df["ad_spend"])
model_simple = sm.OLS(df["sales"], X_simple).fit()
print(model_simple.params)
print(f"R-squared: {model_simple.rsquared:.3f}")

このデータは、売上が実際には「広告費に0.5万円、価格に-0.02万円、気温に1.2万円をそれぞれ掛けた値の合計」に近い形で決まるよう、あらかじめ数値を設定して生成しています。ところが広告費だけを使った単回帰を実行すると、広告費の係数は約0.447、決定係数は0.429という結果になります。真の効果である0.5からずれているうえ、決定係数が示す通り、売上のばらつきのうち広告費だけで説明できるのは全体の43%程度にとどまり、残りの6割弱は価格や気温といった、モデルに含まれていない要因によるばらつきとして扱われてしまっています。なお、この例では広告費・価格・気温を互いに独立に生成しているため、係数が0.5からずれているのは交絡による偏りではなく、他の要因によるばらつきをまるごと誤差として抱え込んだことによる推定の不確かさです。説明変数どうしが相関しているときに係数がどこまで歪むかは、第3章の店舗規模の例で具体的に確認します。次の節以降で見る重回帰では、この3つの要因を同時にモデルに組み込むことで、決定係数も係数の精度も大きく改善します。

重回帰モデルの形：複数の説明変数を同時に扱う

単回帰の限界を解消するのが重回帰です。重回帰は、説明変数を1つに限らず、\( p \)個の説明変数を同時にモデルへ組み込みます。数式で表すと次のようになります。

\( y = b_0 + b_1 x_1 + b_2 x_2 + \cdots + b_p x_p + \varepsilon \)

\( y \)は目的変数、\( x_1, x_2, \ldots, x_p \)は\( p \)個の説明変数、\( b_0 \)は切片、\( b_1, \ldots, b_p \)は各説明変数に対応する回帰係数です。単回帰との違いは、説明変数と回帰係数の組が複数存在する点だけであり、モデルの基本的な考え方は単回帰の自然な拡張になっています。

誤差項\( \varepsilon \)は、モデルに含めた\( p \)個の説明変数だけでは説明しきれなかった、目的変数のばらつきの残り部分を表します。現実には、モデルに含めきれなかった要因（競合の販促活動、天候の急変、突発的なニュースなど）や、そもそも観測が難しい要因、あるいは本質的にランダムな変動が必ず存在するため、誤差項が完全にゼロになることはありません。重回帰モデルを組むという作業は、この誤差項をできる限り小さく、かつ体系的な偏りのない状態に抑え込むために、どの説明変数を組み込むべきかを検討する作業であるとも言い換えられます。

偏回帰係数の意味：「他の変数を一定としたとき」の効果

重回帰を実務で使いこなすうえで、最も重要な概念が偏回帰係数の解釈です。重回帰モデルにおける各係数\( b_1, \ldots, b_p \)は、単回帰の係数とは意味が異なります。重回帰における\( b_1 \)は、「\( x_2, \ldots, x_p \)を一定の値に固定した（統制した）状態で、\( x_1 \)を1単位変化させたときに\( y \)が平均してどれだけ変化するか」を表しています。この「他の変数を一定とした」という条件付きの解釈が、単回帰の係数と重回帰の係数を分ける決定的な違いであり、本ガイド全体で繰り返し用いる基本的な考え方です。

先ほどの売上の例で言えば、重回帰における広告費の係数は「価格と気温が変化しなかったと仮定した場合に、広告費だけを1万円増やしたら売上が平均してどれだけ変わるか」を表します。単回帰における広告費の係数が「価格や気温も一緒に動いてしまう現実の中で、広告費と売上がどれだけ一緒に動いたか」という、いわば要因が混ざったままの相関に近い数値だったのに対し、重回帰の係数は他の要因の影響を切り離した、より純粋な「広告費固有の効果」に近い数値になります。この違いを理解しないまま重回帰の係数を「単純な相関の強さ」として読んでしまうと、施策の効果を見誤る危険があります。

重回帰における各係数は、「その変数単独の効果」ではなく「他の説明変数を一定に保った条件のもとでの効果」を意味します。この条件付きの解釈を抜きにして係数の大きさを比較すると、単回帰と同じ混同を重回帰の中に持ち込んでしまうため、常に「他の変数を一定として」という枕詞を付けて係数を読む習慣が欠かせません。

この「他を一定とする」という操作は、実験室で1つの条件だけを変えて他の条件を固定する対照実験に近い発想です。回帰分析は、実際には何も統制されていない観測データの中から、数式の力を借りて疑似的にこの比較を再現していると捉えることもできます。ただし、この解釈が成り立つためには説明変数どうしの相関が過度に強くないことが前提です。説明変数どうしが強く相関していると、「他の変数を一定として」という条件そのものが現実のデータにほとんど存在しない仮想的な状況になり、係数の推定が不安定になります。この問題（多重共線性）は第5章で詳しく扱います。

最小二乗法の多変数への拡張

重回帰における係数の推定も、単回帰と同じく最小二乗法によって行われます。すなわち、実際の売上\( y_i \)とモデルによる予測値\( \hat{y}_i \)の残差を2乗して全データ分足し合わせた残差平方和を最小にするように、\( b_0, b_1, \ldots, b_p \)のすべてを同時に決定します。単回帰では調整するパラメータが\( b_0, b_1 \)の2つだけでしたが、重回帰では変数の数だけパラメータが増え、それらすべてを同時に最適化する点が異なります。

変数の数が増えると手計算では対応できなくなるため、実務上はすべての説明変数をまとめて行列として扱い、\( \hat{b} = (X^{\top}X)^{-1}X^{\top}y \)という行列計算によって一度に解を求めます（\( X \)は説明変数の行列、\( y \)は目的変数のベクトルです）。この行列表記の詳細まで理解する必要はなく、「変数がいくつ増えても残差平方和を最小にするという発想自体は変わらない」という点さえ押さえておけば、重回帰の仕組みを実務で使ううえで十分です。実際の計算はstatsmodelsのようなライブラリが内部で処理してくれます。

自由度調整済み決定係数：変数を増やすほど決定係数が上がる問題への対処

重回帰を使ってモデルの説明力を評価する際に注意すべき性質があります。それは、説明変数を追加すればするほど、決定係数\( R^2 \)は下がることがなく、必ず上がるか、少なくとも変わらないという性質です。これは、追加した変数が目的変数と本当に関係があるかどうかに関わらず成り立ちます。極端な例では、目的変数と本質的には無関係な、単なる乱数の列を説明変数として加えるだけでも、決定係数はわずかながら必ず上昇します。手元のデータに偶然生じたわずかな相関を、最小二乗法が拾い上げてしまうためです。

この性質があるため、決定係数だけを見て「変数を増やしたらモデルが良くなった」と判断するのは危険です。そこで使われるのが自由度調整済み決定係数（Adjusted \( R^2 \)）です。自由度調整済み決定係数は、説明変数の数\( p \)とデータ数\( n \)を使って次のように計算され、変数を増やすことに対してペナルティを課します。

\( \bar{R}^2 = 1 – (1 – R^2)\dfrac{n – 1}{n – p – 1} \)

意味のある変数を追加すれば自由度調整済み決定係数も上昇しますが、目的変数への説明力がほとんどない変数を追加すると、決定係数のわずかな上昇分がペナルティに食われて自由度調整済み決定係数はむしろ低下します。複数のモデル案（採用する変数の組み合わせが異なるモデル）を比較する際には、決定係数ではなく自由度調整済み決定係数を基準にするのが安全です。以下は、意味のある3つの説明変数（広告費・価格・気温）に、目的変数とは無関係な乱数の変数を5個追加した場合に、決定係数と自由度調整済み決定係数がどう動くかを確認するコードです。

rng = np.random.default_rng(123)
for i in range(5):
    df[f"noise_var{i+1}"] = rng.normal(0, 1, n)

base_cols = ["ad_spend", "price", "temperature"]
noise_cols = [f"noise_var{i+1}" for i in range(5)]

model_base = sm.OLS(df["sales"], sm.add_constant(df[base_cols])).fit()
model_noisy = sm.OLS(df["sales"], sm.add_constant(df[base_cols + noise_cols])).fit()

print(f"3変数モデル: R2={model_base.rsquared:.4f}, 調整済R2={model_base.rsquared_adj:.4f}")
print(f"8変数モデル: R2={model_noisy.rsquared:.4f}, 調整済R2={model_noisy.rsquared_adj:.4f}")

実行すると、3変数モデルでは決定係数0.8029、自由度調整済み決定係数0.7999であるのに対し、無関係な乱数の変数を5個追加した8変数モデルでは、決定係数は0.8069とわずかに上昇する一方、自由度調整済み決定係数は0.7988とむしろ下がります。決定係数だけを見ていると「変数を増やしたのでモデルが改善した」と誤解しかねませんが、自由度調整済み決定係数を確認することで、追加した変数が実質的な説明力を持たない無駄な変数であったことを見抜くことができます。どの変数をモデルに残すべきかというモデル選択の考え方は第6章で体系的に扱います。

statsmodelsによる重回帰の実行と読み方

ここまでの内容を踏まえ、広告費・価格・気温の3つの説明変数をすべて使って売上を説明する重回帰を、statsmodelsのOLS（Ordinary Least Squares、最小二乗法）クラスで実行します。

X = sm.add_constant(df[["ad_spend", "price", "temperature"]])
model = sm.OLS(df["sales"], X).fit()
print(model.summary())

add_constant()は、切片\( b_0 \)を推定するための定数項の列（すべて1の列）を説明変数の行列に追加する関数です。これを忘れると、切片がゼロに固定された（原点を通ることを強制された）モデルになってしまうため、通常はほぼ必ず付けておく処理だと考えて構いません。fit()を呼び出すことで最小二乗法によるパラメータ推定が実行され、summary()メソッドで推定結果の全体像をまとめて確認できます。実際の出力（主要な部分を抜粋）は次のようになります。

OLS Regression Results
==============================================================================
Dep. Variable:                  sales   R-squared:                       0.803
Model:                            OLS   Adj. R-squared:                  0.800
Method:                 Least Squares   F-statistic:                     266.2
No. Observations:                 200   Prob (F-statistic):           7.57e-69
Df Residuals:                      196
Df Model:                            3
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          28.5158      4.409      6.468      0.000      19.821      37.211
ad_spend        0.4983      0.022     22.817      0.000       0.455       0.541
price          -0.0193      0.004     -4.868      0.000      -0.027      -0.011
temperature     1.2206      0.065     18.855      0.000       1.093       1.348
==============================================================================

この出力を、上から順に読み解いていきます。まず上段のブロックはモデル全体に関する情報です。No. Observationsはデータ数（この例では200件）、Df Model（モデルの自由度）は説明変数の数（3個）、Df Residuals（残差の自由度）はデータ数から説明変数の数と切片の分を引いた値（200-3-1=196）を表します。R-squaredは先述の決定係数で、この重回帰モデルでは0.803、すなわち売上のばらつきの80.3%が3つの説明変数によって説明できていることを示します。単回帰（広告費のみ）における決定係数0.429と比べると、価格と気温を加えたことで説明力が大幅に向上したことが分かります。Adj. R-squaredは自由度調整済み決定係数で、この例では0.800とR-squaredに近い値になっており、3つの説明変数がいずれも実質的な説明力を持っていることを裏付けています。

F-statisticとProb （F-statistic）は、「すべての説明変数の係数が同時にゼロである（=このモデルには目的変数を説明する力が全くない）」という仮説を検定した結果です。Prob （F-statistic）がここでの値である7.57e-69のように極めて小さい（0.05を大きく下回る）場合、このモデル全体には統計的に意味のある説明力があると判断できます。

次に、係数の表（下段のブロック）を1列ずつ見ていきます。

係数の解釈を具体的な数値で確認します。広告費（ad_spend）の係数0.4983は、価格と気温を一定に保った状態で広告費を1万円増やすと、売上が平均して約0.50万円（5,000円）増加することを意味します。価格（price）の係数-0.0193は、広告費と気温を一定に保った状態で価格を1円上げると、売上が平均して約0.0193万円（約19円）減少することを、気温（temperature）の係数1.2206は、広告費と価格を一定に保った状態で気温が1度上がると、売上が平均して約1.22万円増加することを、それぞれ表しています。今回はデータを生成した際に設定した真の係数（広告費0.5、価格-0.02、気温1.2）を知っているため答え合わせができますが、推定された係数がいずれも真の値に近い数値になっており、最小二乗法によって元のパラメータがおおむね正しく復元できていることが確認できます。

信頼区間についても補足します。広告費の95%信頼区間は[0.455, 0.541]であり、この範囲に真の係数0.5が収まっていることが分かります。信頼区間が0をまたいでいない（下限・上限ともに正、あるいはともに負）ことは、その係数が統計的に有意である（p値が小さい）こととほぼ対応した情報であり、係数の大きさとその不確実性の幅を同時に確認できる点で、p値だけを見るよりも実務的に有用な情報を含んでいます。

なお、今回の出力にはCond. No.（条件数）が大きいという注意書きが付随することがあります。これは、広告費・価格・気温がそれぞれ大きく異なるスケール（単位）を持っていることに起因するもので、この例では説明変数どうしが強く相関しているために生じる本質的な多重共線性の兆候ではありません。条件数は説明変数の単位に強く左右される数値のため、単独で良し悪しを判断する指標には向きません。すべての説明変数を標準化したうえで計算し、30を超えるようなら説明変数どうしの関係を疑う、という使い方が一般的です。説明変数どうしの相関が実際に問題となるケースの見分け方と対処法は、より直接的な指標であるVIFとあわせて第5章で扱います。

本章では、単回帰の考え方を出発点に、複数の要因を同時に扱う重回帰の基本形と、その中核にある偏回帰係数の解釈、そしてstatsmodelsによる実行結果の読み方までを一通り確認しました。次章では、この重回帰モデルを実務でより丁寧に運用するための作法として、カテゴリ変数の扱いや交互作用項、係数の大きさを比較しやすくする標準化といった、実務上頻出する係数解釈のトピックを扱います。

この章を深めたい方への参考書籍

『回帰分析（新装版）』（佐和隆光、朝倉書店）：数式による導出を丁寧に追いながら回帰分析の理論的な基盤を解説した定評ある一冊で、本章で扱った最小二乗法や偏回帰係数の意味をより厳密に理解したい読者に適しています。実務書というよりは理論書に近い構成のため、数式に抵抗のない読者が本章の内容を深掘りする際の副読本として位置づけられます。

第2章 回帰係数の読み方と伝え方

前章では、単回帰から重回帰へと視野を広げ、複数の説明変数を同時に扱うことで「他の条件を一定とした場合の関係」を捉えられるようになりました。しかし、重回帰分析の出力にある係数の並びは、正しい読み方を知らないまま眺めると、かえって誤った判断を招きます。係数がプラスだから効果がある、数字が大きいから重要、p値が0.05を下回っていないから意味がない、といった単純化は、実務でよく見かける誤読です。この章では、重回帰分析の推定結果を正しく解釈し、経営会議やレポートの場で誤解なく伝えるところまでを、実務作法として順番に扱います。

この章で扱う内容は次の通りです。

• 係数の単位と大きさ：「1単位増えるとyがどれだけ変わるか」という読み方と、単位が違えば大きさを比較できないという制約
• 標準化偏回帰係数：影響度を比較するための道具と、その限界
• p値と信頼区間の読み方：「係数がゼロ」という帰無仮説と、有意でないことの正しい意味
• カテゴリ変数（ダミー変数）の係数解釈：基準カテゴリとの差であること
• 交互作用項：効果が条件によって変わることを表現する方法
• 対数変換した場合の係数解釈：弾力性とパーセント変化の近似
• 係数解釈の落とし穴集：外挿、因果と相関の混同、変数の追加削除による係数の変動
• 経営会議への伝え方：係数表をそのまま見せず、1変数1メッセージに翻訳する

係数の単位と大きさをどう読むか

重回帰分析における偏回帰係数 \( b_1 \) の基本的な読み方は、「他の説明変数をすべて一定に保ったまま、\( x_1 \) が1単位増えたときに、\( y \) が平均的に \( b_1 \) だけ変化する」というものです。この「他を一定として」という条件が、単回帰の係数との大きな違いです。単回帰の係数は\( x \)と\(y\)の間に他の変数を挟まない素朴な関係を表しますが、重回帰の係数は、他の説明変数がもたらす影響を取り除いた上での、その変数固有の関係を表します。例えば、広告費と来店客数と売上を同時にモデルに入れた場合の広告費の係数は、来店客数が同じであれば広告費を増やしても売上にどれだけの上乗せがあるか、という限定つきの効果を示します。

ここで実務上つまずきやすいのが、係数の大きさをそのまま比較して「この変数の方が影響力が大きい」と結論づけてしまうことです。偏回帰係数の大きさは、その説明変数の単位に強く依存します。広告費（万円単位）の係数が0.8で、店舗面積（平方メートル単位）の係数が120だったとしても、これは店舗面積の方が150倍重要だという意味にはなりません。広告費を「円」単位で測り直せば係数は0.00008になりますし、店舗面積を「坪」単位に変えれば係数はまた別の値になります。単位を変えれば数字も変わるという事実を忘れると、係数表を見ただけで変数の重要度を順位づけるという誤りを犯しやすくなります。単位に依存しない形で影響度を比較したい場合には、次に説明する標準化偏回帰係数を使う必要があります。

標準化偏回帰係数：影響度を比較する道具とその限界

説明変数どうしの「効果の大きさ」を単位に左右されずに比較したい場合によく使われるのが、標準化偏回帰係数（standardized partial regression coefficient）です。考え方は標準化（z得点）と同じで、目的変数\(y\)とすべての説明変数をあらかじめ平均0・標準偏差1に変換してから回帰分析を行います。こうして得られる係数は、「説明変数が標準偏差1個分増えたときに、目的変数が標準偏差にして何個分動くか」を表す、単位に依存しない指標になります。非標準化係数から標準化係数を求める式は次の通りです。

\( b_j^{\ast} = b_j \times \dfrac{s_{x_j}}{s_y} \)

ここで \( b_j \) は非標準化係数、\( s_{x_j} \) は説明変数\( x_j \)の標本標準偏差、\( s_y \) は目的変数の標本標準偏差です。標準化係数どうしであれば、単位の異なる広告費（万円）と店舗面積（平方メートル）を、「それぞれが1標準偏差分動いたときに売上を何標準偏差分動かすか」という共通のものさしで比較できます。

ただし、標準化係数には見落とされがちな限界がいくつかあります。第一に、標準化係数はその変数がデータの中でどれだけばらついているか（標準偏差の大きさ）に依存します。店舗面積のばらつきが小さく広告費のばらつきが大きいデータであれば、それだけで広告費の標準化係数が相対的に大きく出ます。これは変数固有の重要性というより、手元のデータでどれだけ変化の幅が観測されたかを反映しているにすぎません。第二に、ダミー変数（0か1しか取らない変数）を標準化すると、その標準偏差は取っている値の比率に依存するため、連続変数の標準化係数と並べて比較する妥当性には議論があります。第三に、後の章で扱う多重共線性（説明変数どうしの強い相関）が存在すると、標準化の有無にかかわらず係数の推定自体が不安定になります。標準化係数は単位をそろえて比較するための便利な道具ではありますが、それが指し示す影響度は、あくまで手元のデータにおける相対的な目安であり、絶対的な重要度のランキングではないという前提を忘れないことが大切です。

p値と信頼区間の読み方

回帰分析の出力には、各係数についてp値と標準誤差、そして多くの場合95%信頼区間が併記されます。係数の検定の統計的な考え方をひとことで整理すると、係数の検定は「この説明変数の真の係数はゼロである（\(y\)に影響を与えていない）」という帰無仮説を出発点に、手元のデータがその帰無仮説のもとでどれだけ起こりにくい結果かをp値で評価する仕組みです。帰無仮説を棄却できるかどうかは、係数の点推定値そのものの大きさだけでなく、その推定値がどれだけ精度よく求まっているか（標準誤差の小ささ）にも左右されます。標準誤差はサンプルサイズや誤差のばらつき、説明変数のばらつきによって決まるため、同じ点推定値でもデータの条件が違えば有意になるかどうかは変わってきます。p値が慣例的なしきい値（0.05など）を下回れば、係数がゼロだとする帰無仮説を棄却し、「統計的に有意な関係がある」と判断します。

ここで実務上の誤解が生まれやすいのが、「有意でない（p値がしきい値を上回った）ならば、その変数には効果がない」という読み替えです。これは正確ではありません。p値が大きいという結果は、「真の効果がゼロである」ことを積極的に示すものではなく、「手元のデータ量やばらつきでは、効果がゼロでないと言い切るだけの証拠が集まらなかった」ことを意味するにすぎません。サンプルサイズが小さい、あるいは説明変数のばらつきが乏しいデータでは、実際には意味のある効果があっても有意という判定に届かないことがよく起こります。逆にサンプルサイズが非常に大きい場合には、実務上ほとんど意味のない小さな効果でも有意と判定されることがあります。p値は「効果の有無」ではなく「偶然のばらつきだけでは説明しにくいかどうか」を語る指標であり、効果の大きさや実務上の重要性を直接教えてくれるものではないという点を、意思決定の場面では意識しておく必要があります。

この限界を補う情報が信頼区間です。95%信頼区間は、係数の点推定値（最も尤もらしい1つの値）だけでなく、「真の係数がどの範囲に収まっていそうか」という幅を示してくれます。例えば、ある施策の効果の係数が「点推定0.8、95%信頼区間は[-0.1, 1.7]」だった場合、点推定だけを見れば効果がありそうに見えますが、信頼区間がゼロをまたいでいる（下限がマイナス）ことから、効果がほぼゼロである可能性も、逆に相応に大きい可能性も、データからは排除しきれていないことが分かります。反対に「点推定0.8、95%信頼区間は[0.6, 1.0]」であれば、効果の大きさにある程度の幅はあるものの、ゼロに近い値である可能性はかなり低いと言えます。同じ点推定値でも、信頼区間の広さによって意思決定の確からしさはまったく違ってきます。p値による有意・非有意の二択の判定だけでなく、信頼区間の幅そのものを読む習慣を持つことが、係数を正しく解釈する上での基本作法です。

カテゴリ変数（ダミー変数）の係数解釈

説明変数が「地域」や「顧客セグメント」「曜日」のようなカテゴリ（質的変数）である場合、そのままでは回帰式に投入できないため、ダミー変数（0か1の値を取る変数）に変換して扱います。カテゴリがA・B・Cの3種類であれば、そのうち1つ（例えばA）を基準カテゴリとして選び、残りのB・Cそれぞれについて「そのカテゴリに該当すれば1、それ以外は0」というダミー変数を2つ作成し、回帰式に投入します。カテゴリの数が\(k\)個であれば、作成するダミー変数は\(k-1\)個になります。

ここで重要なのは、ダミー変数の係数が「そのカテゴリの絶対的な水準」ではなく「基準カテゴリと比べたときの差」を表しているという点です。地域ダミー（基準:東京）の「大阪」ダミーの係数が\(-50\)（万円）だったとすると、これは「他の説明変数を一定としたとき、大阪は東京と比べて売上が平均50万円少ない」という意味であり、大阪の売上水準そのものを表しているわけではありません。基準カテゴリを東京から大阪に変更すれば、今度は東京ダミーの係数が「大阪と比べた東京の差」として\(+50\)前後の値で出てきます。どのカテゴリを基準に選ぶかによって係数表に並ぶ数字の見え方（どのカテゴリの差が目立つか、符号がプラスかマイナスか）は変わりますが、これはモデルの当てはまりの良さ（決定係数やF検定の結果）を変えるものではなく、あくまで表現上の違いにすぎません。実務上は、比較の基準として最も自然なカテゴリ（全社平均に近いセグメント、施策前の状態など）を基準に選ぶと、係数表がそのまま報告に使いやすくなります。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(0)
n = 300

region = rng.choice(["東京", "大阪", "名古屋"], size=n, p=[0.5, 0.3, 0.2])
ad_spend = rng.normal(100, 20, n)

# 地域ごとの水準差を仕込んで生成
base = {"東京": 0, "大阪": -50, "名古屋": -20}
sales = (
    300
    + 2.0 * ad_spend
    + np.array([base[r] for r in region])
    + rng.normal(0, 30, n)
)
df = pd.DataFrame({"sales": sales, "ad_spend": ad_spend, "region": region})

# 東京を基準カテゴリにしたモデル
model_tokyo = smf.ols(
    "sales ~ ad_spend + C(region, Treatment(reference='東京'))", data=df
).fit()
print(model_tokyo.params)

# 大阪を基準カテゴリに変えたモデル
model_osaka = smf.ols(
    "sales ~ ad_spend + C(region, Treatment(reference='大阪'))", data=df
).fit()
print(model_osaka.params)
print(model_tokyo.rsquared, model_osaka.rsquared)

このコードでは、statsmodelsのformula API（smf.ols）を使い、C()関数でカテゴリ変数を指定しています。Treatment(reference='東京')のように基準カテゴリを明示的に指定できる点が実務上便利です。基準カテゴリを東京から大阪に変えても、決定係数（rsquared）は変わらず、係数の数字の表示のされ方だけが変わることをこのコードで確認できます。

交互作用項：効果が条件によって変わることを表現する

ここまでの重回帰モデルは、各説明変数の効果が他の説明変数の水準にかかわらず一定であることを前提にしていました。しかし実務では、「広告費の効果は季節によって変わる」「値引きの効果は顧客セグメントによって変わる」といった、条件によって効果の大きさそのものが変化する現象がよく見られます。この「効果が条件によって変わる」という関係を回帰モデルに組み込む道具が交互作用項（interaction term）です。

広告費（\(x_1\)）の効果が季節（\(x_2\)）によって変わる状況を式で表すと、次のようになります。

\( y = b_0 + b_1 x_1 + b_2 x_2 + b_3 (x_1 \times x_2) + \varepsilon \)

交互作用項を含むモデルでは、\(x_1\)の効果を表す偏微分（\(x_1\)が1単位増えたときの\(y\)の変化）は \( b_1 + b_3 x_2 \) となり、\(x_2\)の値によって変わることが式の上からも分かります。ここで注意したいのが、\(b_1\)（主効果の係数）の意味です。\(b_1\)は「\(x_2\)がゼロのときの\(x_1\)の効果」を表しており、すべての水準を通じた平均的な効果ではありません。季節がカテゴリ変数の場合、主効果\(b_1\)は「基準カテゴリの季節における広告費の効果」を意味し、他の季節ダミーとの交互作用項の係数は、その季節において効果が基準の季節と比べてどれだけ違うかを表します。連続変数どうしの交互作用では、\(x_2\)がゼロという値が現実的でないことも多いため、あらかじめ\(x_2\)を平均で中心化（データから平均を引く）しておくと、主効果の係数が「\(x_2\)が平均的な水準にあるときの\(x_1\)の効果」という、解釈しやすい意味を持つようになります。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(1)
n = 400

season = rng.choice(["春", "夏", "秋", "冬"], size=n)
ad_spend = rng.normal(100, 25, n)

# 広告費の効果が季節で変わるデータを生成
effect = {"春": 1.5, "夏": 3.0, "秋": 1.2, "冬": 0.8}
sales = (
    250
    + np.array([effect[s] for s in season]) * ad_spend
    + rng.normal(0, 40, n)
)
df = pd.DataFrame({"sales": sales, "ad_spend": ad_spend, "season": season})

# 交互作用項を含むモデル(春を基準カテゴリとする)
model_int = smf.ols(
    "sales ~ ad_spend * C(season, Treatment(reference='春'))", data=df
).fit()
print(model_int.params)
print(model_int.pvalues)

formula APIでは「ad_spend * C(season, ...)」と書くだけで、広告費の主効果・季節ダミーの主効果・両者の交互作用項を自動的に展開してくれます（「ad_spend + C(season, ...) + ad_spend:C(season, ...)」と同じ結果になります）。出力される交互作用項の係数（例えば「ad_spend:C(season)[T.夏]」）がプラスで有意であれば、基準カテゴリの春と比べて、夏は広告費1単位あたりの押し上げ効果が大きい、と解釈します。連続変数どうしの交互作用では、中心化した変数を使ってformulaを組み立てると主効果の係数が読みやすくなります。

交互作用項のある回帰式は、主効果の係数だけを取り出して「この変数の効果は\(b_1\)である」と単純に語ることができません。交互作用の相手となる変数がどの値・どのカテゴリを取っているかによって、実際の効果は変わります。交互作用を含むモデルの結果を報告する際は、係数表をそのまま渡すのではなく、「代表的な条件（季節、セグメントなど）ごとに効果がどう変わるか」を具体的な数値や図に翻訳して伝えることが欠かせません。

対数変換した場合の係数解釈

売上や所得のように右に裾を引く分布を扱う際、目的変数や説明変数を対数変換してから回帰分析を行うことがあります。対数変換を行うと、係数の意味は通常の線形回帰とは変わり、パーセント（割合）で解釈できるようになります。

目的変数と説明変数の両方を対数変換したモデル（log-logモデル、\(\ln y = b_0 + b_1 \ln x + \varepsilon\)）では、係数\(b_1\)は経済学でいう弾力性（elasticity）そのものであり、「\(x\)が1%増えたとき、\(y\)がおよそ\(b_1\)%変化する」とそのまま読むことができます。価格弾力性（価格が1%上がったときに需要が何%変わるか）を求める分析は、この形の代表例です。

一方、目的変数だけを対数変換し、説明変数は元のスケールのまま扱うモデル（log-linearモデル、\(\ln y = b_0 + b_1 x + \varepsilon\)）では、係数の読み方に注意が必要です。厳密には、\(x\)が1単位増えたときに\(y\)は \( (\exp(b_1) – 1) \times 100 \) %変化する、というのが正確な解釈です。\(b_1\)の値が小さい（目安として0.1以下程度）場合には、\(\exp(b_1) – 1\)は\(b_1\)にほぼ等しくなるため、「\(x\)が1単位増えるとyがおよそ\(b_1 \times 100\)%変化する」という近似が広く使われています。しかし\(b_1\)が大きい値になると近似の誤差は無視できません。例えば\(b_1=0.5\)の場合、近似では「50%の変化」となりますが、正確な計算では \( (\exp(0.5)-1)\times100 \approx 64.9 \) %であり、両者には10%以上の開きが生じます。近似で済ませてよいかは係数の大きさを見て判断します。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(2)
n = 500

price = rng.uniform(500, 3000, n)
loyalty_years = rng.uniform(0, 10, n)

# 価格弾力性-0.8、継続年数の効果+0.05を仕込んだ需要データ
demand = np.exp(
    5.0
    - 0.8 * np.log(price)
    + 0.05 * loyalty_years
    + rng.normal(0, 0.3, n)
)
df = pd.DataFrame({"demand": demand, "price": price, "loyalty_years": loyalty_years})

# log-logモデル(価格弾力性を推定)
model_loglog = smf.ols("np.log(demand) ~ np.log(price) + loyalty_years", data=df).fit()
print(model_loglog.params)

b_loyalty = model_loglog.params["loyalty_years"]
approx_pct = b_loyalty * 100
exact_pct = (np.exp(b_loyalty) - 1) * 100
print(f"近似: {approx_pct:.2f}%, 正確な値: {exact_pct:.2f}%")

formula APIの中でnp.log()をそのまま書けるのがsmf.olsの利点で、np.log(demand)、np.log(price)のように目的変数・説明変数を個別に変換できます。log(price)の係数は価格弾力性としてそのまま読み、loyalty_yearsのように元のスケールのままの説明変数の係数は、近似のパーセント表示と\(\exp(b)-1\)による正確なパーセント表示の両方を確認する習慣をつけておくと、報告時の誤差を避けられます。

係数解釈の落とし穴集

係数解釈で特に見落とされやすい3つの落とし穴を整理しておきます。

• 外挿の危険：回帰係数は、手元のデータが観測された範囲の中でしか、その関係の妥当性を保証していません。広告費が月50万円から200万円の範囲でしか観測されていないのに、係数を使って「広告費を1000万円に増やせば売上はこれだけ伸びる」と予測することは、データの裏付けのない外挿です。多くの実務上の関係は、ある範囲を超えると効果が頭打ちになったり、逆転したりすることが珍しくありません。
• 因果と相関の混同：偏回帰係数は「他の説明変数を一定としたときの、統計的な関連の強さ」を表しますが、これは自動的に「\(x\)を変化させれば\(y\)が変化する」という因果関係を意味するわけではありません。両者を混同すると、施策の効果を過大に見積もる危険があります。この論点は次の第3章で、交絡変数の調整という切り口から詳しく扱います。
• 変数の追加・削除による係数の変動：モデルに新しい説明変数を加えたり、既存の変数を取り除いたりすると、残った変数の係数の大きさや符号、有意性が変わることがあります。これは、説明変数どうしが互いに相関している（多重共線性がある）場合に典型的に起こる現象で、「なぜモデルを少し変えただけで結論が変わるのか」という疑問につながりやすい部分です。この現象の仕組みと対処法は、第5章の多重共線性で詳しく取り上げます。

経営会議への伝え方：係数表をそのまま見せない

回帰分析の結果を経営会議やレポートで共有する場面では、統計ソフトが出力する係数表（偏回帰係数・標準誤差・t値・p値の表）をそのままスライドに貼り付けることは避けた方がよい場合がほとんどです。係数表は分析者には必要な情報の一覧ですが、統計の専門用語に馴染みのない読み手には、数字の羅列以上の意味を持ちません。

実務で有効なのは、係数表の情報を「1つの変数につき1つの意味のあるメッセージ」に翻訳し直す作業です。例えば「広告費の係数は2.3、95%信頼区間は[0.8, 3.8]」という1行は、「広告費を100万円追加すると、売上はおよそ80万円から380万円の範囲で増え、平均的には230万円程度の増加が見込まれる」という文章に置き換えられます。信頼区間が狭ければ「確度の高い見積もり」、広ければ「方向性は見えるが幅のある見積もり」という確からしさの違いも一緒に伝わります。「有意/非有意」という専門用語でそのまま報告するよりも、信頼区間の幅を使って効果の大きさへの自信の度合いを語る方が、意思決定者にとって扱いやすい情報になります。

特に3行目のように、有意水準を満たさなかった係数を「効果がない」と言い切って報告すると、実際には一定の効果がある施策を早々に打ち切る判断につながりかねません。有意でない結果は「効果がないと確認された」のではなく「効果の有無を判断するにはデータが足りない、あるいはばらつきが大きい」という留保つきの結果である、という伝え方をする必要があります。分析者には統計的に正確な解釈と、意思決定者に伝わる平易な言葉づかいの両方が求められます。

この章を深めたい方への参考書籍

『入門はじめての多変量解析』（石村貞夫・石村光資郎、東京図書）：回帰分析をはじめとする多変量解析の考え方を、数式の展開を丁寧に追いながら解説した入門書です。偏回帰係数や標準化係数、ダミー変数の扱いといったこの章のテーマを、手を動かしながら基礎から確認したい場合の土台になります。

第3章 説明のための回帰、予測のための回帰と交絡の調整

前章までは、複数の説明変数を使った重回帰の基本的な組み方と、係数をどう読むかという実務作法を見てきました。本章では、回帰分析という道具の使い方を左右する、より根本的な論点を扱います。それは「その回帰分析を、何のために行うのか」という目的の問題です。技術的な手続きは同じでも、目的が変われば正しい進め方が変わる、というのが本章で伝えたい要点です。

同じ最小二乗法、同じ回帰係数、同じ決定係数を使っていても、分析の目的が「なぜその結果になったのかを説明したい」のか、「次に起きることをできるだけ正確に当てたい」のかによって、変数を選ぶ基準も、共線性への向き合い方も、モデルの評価方法もまったく別物になります。この違いを意識せずに回帰分析を組み立ててしまうと、説明力の高い変数を予測モデルから削ってしまったり、逆に予測精度だけを追い求めた変数選択の結果を経営判断の根拠にしてしまったりする失敗が起こります。本章で扱う「説明と予測の違い」は、本ガイド全体を貫く中心的なテーゼの1つです。

説明のための回帰と予測のための回帰は何が違うか

回帰分析を使う目的は、大きく2つに分けられます。1つは「ある結果に対して、どの要因がどの程度影響しているのかを明らかにしたい」という説明目的です。広告費と店舗改装のどちらが売上への影響が大きいか、離職率に効いているのは残業時間か上司との関係か、といった問いがこれにあたります。もう1つは「入力データから、まだ分かっていない結果をできるだけ正確に当てたい」という予測目的です。来月の需要量を当てる、退職しそうな社員を事前に検知する、といった問いがこれにあたります。

この2つの目的は、同じ回帰分析という技術を使うにもかかわらず、実務上の作法が大きく異なります。次の表に主な違いを整理します。

この表からも分かるとおり、多重共線性は説明目的では致命的な問題になりやすい一方、予測目的では致命傷にならないことがあります。複数の説明変数が互いに強く相関していても、それらをまとめて使った予測値そのものは安定していることが多いためです。ただし個々の係数を見て「この変数の影響は大きい/小さい」と解釈しようとする場面では、共線性が大きな問題になります。多重共線性そのものの診断方法（VIFなど）は第5章で詳しく扱います。また予測目的におけるモデル選択と交差検証の考え方は第6章で、機械学習の枠組みでの汎化性能の評価は本シリーズの別冊（機械学習の実践ガイド）でより深く扱います。

回帰分析に取り組む前に、まず「これは説明のための分析か、予測のための分析か」を自分自身に問い直すことが、変数選択からモデル評価まですべての工程の方針を決めます。目的を曖昧にしたまま作業を進めると、どちらの基準で判断すればよいか分からないまま迷走しやすくなります。

単純な集計が誤る理由と回帰による交絡の調整

ここからは、説明のための回帰がなぜ単純な平均比較やクロス集計より優れているのかを、具体的な数値例で確認します。鍵になるのは「交絡（confounding）」という概念です。交絡とは、注目している2つの変数（たとえば施策の有無と成果）の両方に影響を与える第三の変数が存在するために、見かけ上の関係が歪んでしまう現象を指します。

次のような架空の状況を考えます。ある小売チェーンが、業績支援のための販促施策を一部の店舗に導入しました。施策導入後の売上成長率を、施策を導入した店舗と導入していない店舗で単純に比較したところ、次のような結果になったとします。

• 施策を導入した店舗の平均売上成長率:6.2パーセント
• 施策を導入していない店舗の平均売上成長率:6.3パーセント

この数字だけを見ると、「施策には効果がない、むしろわずかにマイナスかもしれない」という結論になりそうです。しかしこのチェーンでは、業績が伸び悩んでいた大型店舗を中心に施策を投入していたという経緯がありました。大型店舗はすでに売上規模が大きく、成長余地そのものが小型店舗より小さいため、施策の有無にかかわらず成長率は低く出やすい傾向があります。つまり「店舗規模」が、施策の有無（大型店ほど施策が投入されやすい）と、売上成長率（大型店ほど成長率が低い）の両方に影響する交絡変数になっているのです。この構造を無視して施策あり/なしの単純平均だけを比べると、施策の効果が過小評価され、符号までひっくり返って見えてしまいます。層別クロス集計で見ると規模ごとの傾向が逆転するように見える現象は、俗にシンプソンのパラドックスと呼ばれるものの一種で、回帰分析における交絡問題はこれの連続変数版と考えると理解しやすくなります。

実際に店舗規模を大型/小型の2群に分けて、群ごとに施策あり/なしの平均を比較すると、次のように施策の効果がはっきり見えてきます。

店舗規模ごとに見れば、小型店・大型店のどちらでも施策ありの成長率が施策なしを2から3ポイント程度上回っています。ところが施策は大型店に偏って投入されていたため、規模を無視して全体を平均すると、成長率が低い大型店のデータが施策ありグループに多く混ざり、成長率が高い小型店のデータが施策なしグループに多く混ざる形になります。この結果、施策ありグループの平均が押し下げられ、施策なしグループの平均が押し上げられて、見かけ上は効果がないかのような数字になってしまうのです。

この層別クロス集計と同じ調整を、回帰分析では「店舗規模を説明変数としてモデルに含める」という形で自動的に行うことができます。次のコードは、この状況をシミュレーションで再現し、店舗規模を無視した回帰と、店舗規模を統制した回帰の係数を比較する例です。

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

rng = np.random.default_rng(42)
n = 300

# 店舗規模(1=大型店, 0=小型店)
store_size = rng.binomial(1, 0.5, n)
# 施策は大型店に手厚く投入される(業績支援の必要性が高いと判断されたため)
p_treat = np.where(store_size == 1, 0.75, 0.25)
treatment = rng.binomial(1, p_treat)

true_effect = 3.0  # 施策が売上成長率に与える真の効果(ポイント)
noise = rng.normal(0, 4, n)
# 大型店であること自体が成長率を6ポイント押し下げる、という設定
sales_growth = 8 - 6 * store_size + true_effect * treatment + noise

df = pd.DataFrame({
    'store_size': store_size,
    'treatment': treatment,
    'sales_growth': sales_growth,
})

# 店舗規模を無視した(交絡を統制しない)回帰
X_naive = sm.add_constant(df[['treatment']])
model_naive = sm.OLS(df['sales_growth'], X_naive).fit()

# 店舗規模を説明変数に加えた(交絡を統制した)回帰
X_adj = sm.add_constant(df[['treatment', 'store_size']])
model_adj = sm.OLS(df['sales_growth'], X_adj).fit()

print('店舗規模を無視した場合の施策係数:', round(model_naive.params['treatment'], 2))
print('店舗規模を統制した場合の施策係数:', round(model_adj.params['treatment'], 2))
print('店舗規模を統制した場合の店舗規模係数:', round(model_adj.params['store_size'], 2))

このコードを実行すると、店舗規模を無視した回帰では施策の係数がおよそマイナス0.1と、効果がほぼゼロかわずかにマイナスという結果になります。これは前述の単純平均比較と同じ結論です。一方、店舗規模を説明変数として加えた回帰では、施策の係数はおよそプラス2.9となり、シミュレーションで設定した真の効果（プラス3.0）にかなり近い値が回帰係数として復元されます。店舗規模の係数もおよそマイナス5.9と、設定したマイナス6.0に近い値です。このように、交絡変数を説明変数としてモデルに加えるだけで、単純集計では見えなかった効果を回帰分析が正しく取り出してくれるのです。これが「単純集計より重回帰が有効である」という主張の具体的な中身になります。

「他を一定として」の意味と限界

重回帰分析で複数の説明変数を同時にモデルへ入れることの本質は、「ある変数の係数を、他の変数の値を一定に保った状態での影響として推定する」という点にあります。前節の例で言えば、「店舗規模が同じ条件のもとで、施策の有無が売上成長率にどれだけ差をつけるか」を取り出す操作にあたります。この「他を一定として」という発想があるからこそ、回帰分析は単純な平均比較よりも交絡に強い分析手法になっています。

ただしこの力には限界があります。第一に、モデルに加えていない変数については、当然ながら「一定に保つ」ことができません。次節で扱う「観測されていない交絡」の問題がこれにあたります。第二に、「他を一定として」という操作自体が、現実には起こりえない、あるいは意味をなさない仮想的な比較になっている場合があります。たとえば「従業員数を一定として店舗の売場面積だけを増やした場合の売上効果」を回帰係数から読み取ろうとしても、実際の店舗運営では売場面積を広げれば従業員数も同時に増やすのが通常であり、「従業員数だけを固定して売場面積だけを動かす」という状況自体が現実離れしている可能性があります。回帰係数はあくまで数式上の操作であり、その操作が業務上の意思決定として意味を持つかどうかは、分析者が別途吟味する必要があります。

統制すべきでない変数：媒介変数とコライダー

「交絡変数は回帰に加えて統制するべき」という考え方を裏返して、「変数はとにかく多く入れておけば安全」と誤解してしまうと、今度は別の落とし穴にはまります。回帰分析に加えるべきではない変数が、少なくとも2種類存在するためです。

1つ目は媒介変数（mediator）です。媒介変数とは、注目している原因（施策など）が結果（成果）に影響を及ぼす際に、その効果が実際に伝わっていく経路上に位置する変数のことです。たとえば「研修プログラムの受講」が「従業員のスキルスコア向上」を経て「売上成長率の向上」につながっている、という構造を考えます。この場合、スキルスコアは研修と売上成長率の間の媒介変数にあたります。もし回帰モデルに研修の受講有無とスキルスコアの両方を説明変数として入れてしまうと、研修が売上成長率に与える効果のうち、スキルスコアを経由する部分がスキルスコアの係数に吸収されてしまい、研修そのものの係数は本来の効果より小さく、極端な場合はほぼゼロに見えてしまいます。効果の通り道そのものを統制してしまうと、効果が消えて見えるのです。

次のコードは、この現象を簡単なシミュレーションで確認する例です。研修の受講（training）がスキルスコア（skill_score）を経由してのみ売上成長率（sales_growth）に影響するという構造を設定し、媒介変数を統制する前と後で研修の係数がどう変わるかを見ます。

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

rng = np.random.default_rng(7)
n = 300

training = rng.binomial(1, 0.5, n)  # 研修プログラムの受講有無

a = 8.0  # 研修がスキルスコアに与える効果
skill_score = 50 + a * training + rng.normal(0, 5, n)

b = 0.6  # スキルスコアが売上成長率に与える効果
# 研修そのものには、スキルスコアを経由する以外の直接効果はない設定
sales_growth = 10 + b * skill_score + rng.normal(0, 3, n)

df = pd.DataFrame({
    'training': training,
    'skill_score': skill_score,
    'sales_growth': sales_growth,
})

# 媒介変数(スキルスコア)を統制しない回帰、総合効果を推定
X1 = sm.add_constant(df[['training']])
model_total = sm.OLS(df['sales_growth'], X1).fit()

# 媒介変数を統制した回帰
X2 = sm.add_constant(df[['training', 'skill_score']])
model_direct = sm.OLS(df['sales_growth'], X2).fit()

print('媒介変数を統制しない場合の研修係数:', round(model_total.params['training'], 2))
print('媒介変数を統制した場合の研修係数:', round(model_direct.params['training'], 2))

このシミュレーションでは、研修がスキルスコアを8ポイント押し上げ、スキルスコア1ポイントが売上成長率を0.6ポイント押し上げる設定にしているため、研修の総合的な効果は理論上8かける0.6でおよそ4.8になります。実際にスキルスコアを統制しない回帰を実行すると、研修の係数はおよそ5.1とこの理論値に近い値になります。ところがスキルスコアを説明変数として加えた回帰では、研修の係数はおよそ0.5まで縮小し、ほぼ効果がないかのような結果になります。これは研修の効果が消えたのではなく、効果の伝達経路であるスキルスコアの係数（およそ0.62）にほぼ吸収されてしまったために起きています。「研修に効果があるかどうか」を知りたいのであれば、媒介変数であるスキルスコアを説明変数に加えるべきではありません。

2つ目はコライダー（collider、合流点）と呼ばれる変数です。コライダーとは、2つの独立した原因が共通して影響を与えている結果側の変数のことを指します。たとえば「営業スキル」と「運の良さ」は本来互いに無関係だとしても、両方が「大口契約の獲得」という共通の結果に影響しているとします。この「大口契約の獲得」がコライダーです。ここで、大口契約を獲得した社員だけに絞り込んで営業スキルと運の良さの関係を調べると、奇妙な負の関係が浮かび上がることがあります。大口契約を獲得したという条件のもとでは、営業スキルが低い社員は運がよほど良くなければそこに含まれず、逆に運が悪かった社員は営業スキルがよほど高くなければそこに含まれないためです。結果として、営業スキルと運の良さが本来は無関係であっても、コライダーで層別した集団の中だけを見ると見かけ上の関連が生まれてしまいます。回帰分析でコライダーにあたる変数（あるいはその代理変数）をうっかり説明変数に含めてしまうと、これと同じ歪みがモデルの中に持ち込まれます。どの変数が交絡でどの変数が媒介やコライダーなのかは、データの相関関係だけからは判別できません。変数同士の因果関係についての仮説を、分析に着手する前に整理しておく必要があります。

回帰による因果推論の限界：観測されていない交絡

ここまでの議論で、回帰分析に交絡変数を加えることで単純集計よりも実態に近い効果を推定できることを見てきました。ただしここには重要な前提が隠れています。それは「交絡になっている変数が、データとして観測できている」という前提です。

観測データを使った回帰分析で調整できるのは、あくまでモデルに説明変数として含めた、観測された交絡だけです。実際の意思決定の現場では、施策の有無を決める判断にどのような要因が影響していたかを、分析者があとから完全に把握できるとは限りません。店長の経営手腕、地域の競合状況、施策導入前の一時的な業績変動など、データとして記録されていない要因が施策の割り当てと成果の両方に影響していれば、それは未観測交絡としてモデルの外に残り続けます。未観測交絡が存在する限り、回帰係数がどれだけ理論値に近い値を示していても、それが真の因果効果と一致している保証はありません。

この限界を正面から扱うのが、ランダム化比較試験（RCT）や差分の差分法、傾向スコア、操作変数法といった因果推論の手法群です。これらの手法は、観測データだけでは解決できない未観測交絡の問題に、実験計画や準実験的な状況設定によって別の角度から取り組みます。本ガイドは回帰分析と統計モデリングの技術的な扱い方を中心に扱う巻であるため、因果推論の枠組みそのものについては本シリーズの別冊（因果推論と効果検証の実践ガイド）で改めて体系的に扱う予定です。ここでは「回帰分析は魔法のように因果関係を取り出す道具ではなく、観測された交絡を調整する道具である」という限界を押さえておくことが重要です。

回帰分析で交絡を統制できるのは、あくまで説明変数としてモデルに含めた、観測済みの変数に限られます。未観測の交絡が存在する可能性を常に念頭に置き、「この回帰係数は本当に因果効果と呼べるのか、それとも観測できた範囲での調整済みの相関にとどまるのか」を区別する姿勢が、説明目的の回帰分析には欠かせません。

変数を追加すると係数が動く理由：省略変数バイアス

前節までの数値例でも見たとおり、モデルに変数を1つ加えるだけで、他の変数の係数が大きく変化することがあります。この現象は省略変数バイアス（omitted variable bias）と呼ばれ、直感的には次のように理解できます。

ある変数Xの回帰係数は、本来「Xの効果」と「モデルに入れていない変数の効果のうち、Xと相関している部分」が混ざり合った値になります。モデルに交絡変数Zを追加していない場合、Xの係数には、Zが持つ本来の効果のうち、XとZの相関を通じて漏れ込んでくる分がそのまま乗ってしまいます。店舗規模の例で言えば、施策の係数には本来の施策効果（プラス3程度）に加えて、店舗規模の効果（マイナス6程度）のうち施策と店舗規模の相関（大型店ほど施策を受けやすい）を通じて漏れ込んだ分が上乗せされ、結果としてほぼゼロないしわずかにマイナスの係数になっていました。Zをモデルに追加すると、この漏れ込みがZの係数として明示的に切り出されるため、Xの係数は本来の効果に近づいていきます。

この直感から実務上重要な示唆が2つ導かれます。1つは、係数がモデルに入れる変数の組み合わせによって変わるのは異常なことではなく、むしろ説明変数の間に相関がある限り自然に起こる現象だという点です。もう1つは、変数を追加した際に係数が大きく動いた場合、それは無視できないほど強い交絡がそこに存在していたことを示すシグナルとして受け止めるべきだという点です。逆に、変数を追加しても係数がほとんど動かないのであれば、その変数は少なくとも今回のモデルにおいて深刻な交絡源ではなかったと判断できます。

実務ガイドライン：目的から逆算して変数を選ぶ

本章の内容を実務でどう使うかを整理します。まず分析に着手する前に、その回帰分析が説明目的なのか予測目的なのかを明確にします。

説明目的の場合は、変数を機械的に総当たりで試すのではなく、まず「何が何に影響していると考えられるか」という因果関係の見取り図を簡単にでも描いてから、説明変数を選ぶことが望ましいやり方です。見取り図を描く際には、注目している原因と結果の両方に影響を与えている交絡変数を洗い出して説明変数に含める一方、原因と結果の間をつなぐ媒介変数や、結果側で合流するコライダーにあたる変数は、統制すべきではない変数として除外を検討します。この見取り図はごく簡単な箇条書きや矢印の図で構いません。重要なのは、モデルに入れる変数と入れない変数を、相関の強さだけでなく、変数同士の関係についての仮説に基づいて選ぶという手順を踏むことです。

• その回帰分析は説明目的か予測目的か、着手前に明確にする
• 説明目的なら、注目する原因と結果の両方に影響する交絡変数を洗い出してから説明変数を選ぶ
• 媒介変数やコライダーにあたる変数は、統制すべきではない変数として除外を検討する
• 変数を追加して係数が大きく動いたら、そこに交絡があったサインとして受け止める
• 観測データの回帰で調整できるのは観測された交絡だけであることを常に意識する

一方、予測目的の場合は、説明目的ほど因果関係の仮説にこだわる必要はありません。極端に言えば、予測に有効であれば理論的な意味づけが曖昧な変数を使っても構わない場面すらあります。その代わりに重視すべきは、手元のデータだけに過剰に適合していないか、まだ見ていないデータに対してもきちんと当てはまるかという汎化性能の確認です。この確認には交差検証（cross-validation）と呼ばれる手法が広く使われます。交差検証の具体的な手順とモデル選択の考え方は第6章で、機械学習的な観点からの汎化性能の評価はさらに本シリーズの別冊で扱います。説明目的の作法をそのまま予測目的の作法として持ち込んでしまう、あるいはその逆をしてしまうという取り違えが、実務における回帰分析の混乱の多くを引き起こしています。本章で整理した目的の違いを、次章以降の技術的な内容を読み進める際の軸として持っておいてください。

この章を深めたい方への参考書籍

『調査観察データの統計科学』（星野崇宏、岩波書店）：ランダム化されていない観測データから因果効果を推定する際の考え方を、傾向スコアなどの手法とあわせて体系的に解説した書籍です。本章で触れた「観測された交絡しか調整できない」という限界を、より厳密な統計的枠組みで理解したい読者に適しています。

第4章 回帰の前提と残差診断

前章までで、単回帰から重回帰へと話を広げ、係数の解釈の作法や、交絡を調整して「説明のための回帰」を組み立てる考え方を見てきました。回帰分析は、目的変数と説明変数の列さえそろっていれば、statsmodelsであれscikit-learnであれ、数行のコードで係数とp値、決定係数を返してくれます。この手軽さは回帰分析の大きな利点です。ただし、コードが結果を返したこと自体は、その結果が信頼できることを保証しません。

最小二乗法によるOLS（Ordinary Least Squares）は、データがどのような性質を持っていても、残差の二乗和を最小にする係数を計算して返します。説明変数と目的変数の関係が実は曲線的であっても、データの中に極端な値が紛れ込んでいても、モデルの前提が大きく崩れていても、model.fit()は何かしらの係数とp値を出力し、エラーは出ません。つまり、出力された数字が信頼できるかどうかを判断する責任は、そのまま分析者の側に残ります。計算が通ったことと、結果を経営判断の根拠に使えることは別の問題だと考えています。

本章では、回帰モデルを当てはめた後に確認すべき「前提」と、その前提が満たされているかを確認するための「残差診断」の手順を整理します。前提が崩れたときに何が起きるのかを正しく理解しておくことは、モデルの出力をそのまま鵜呑みにせず、経営判断の材料として使ってよいかどうかを見極めるうえで欠かせない作業です。

線形回帰の主な前提を整理する

OLSによる線形回帰が「良い性質を持つ推定量である」と言えるのは、いくつかの前提が満たされている場合に限られます。代表的な前提は次の4つです。

• 線形性：目的変数と説明変数の関係が、係数に関して線形の式で近似できること
• 誤差の独立性：ある観測の誤差が、他の観測の誤差と無関係であること
• 等分散性：誤差のばらつきの大きさが、説明変数の値によらず一定であること
• 正規性：誤差が正規分布に従うこと（検定や区間推定といった推論の場面で関係する前提）

ここで強調しておきたいのは、これらの前提が崩れたときに「何が信用できなくなるか」は、前提の種類によって異なるという点です。すべてをひとくくりに「前提違反だからモデルは無効」と考えるのではなく、どの部分が壊れているのかを切り分けて考える必要があります。

まず線形性が崩れている場合です。本来は曲線的な関係にあるデータに直線をあてはめてしまうと、係数の推定値そのものが体系的にゆがみます（バイアスを持ちます）。これは点推定の段階で信頼性を失う、最も深刻な前提違反だといえます。

次に誤差の独立性です。時系列データを回帰する際などによく問題になりますが、ある時点の誤差が前後の時点の誤差と相関している（自己相関がある）場合、係数の点推定自体はゆがみませんが、標準誤差は過小に見積もられがちです。標準誤差が実際より小さく出るということは、p値も実際より小さく出やすいということであり、本当は有意でない関係を「有意である」と誤って判断してしまうリスクが高まります。

等分散性が崩れている場合（不均一分散、英語ではheteroscedasticityと呼ばれます）も、係数の点推定自体には影響しませんが、標準誤差の計算が歪みます。通常のOLSの標準誤差は、誤差の分散がどの観測でも一定であるという前提のもとで導かれているため、この前提が崩れると、標準誤差が過大にも過小にもなり得ます。結果として、p値や信頼区間の解釈を誤るリスクが生まれます。

最後に正規性です。ここは誤解されやすい点なので、はっきりさせておきます。誤差が正規分布に従うという前提は、OLSによる係数の点推定そのものには必須ではありません。ガウス・マルコフの定理により、説明変数と誤差項が無相関で（この条件を外生性と呼びます）、誤差の分散が有限で均一であり独立である限り、OLS推定量は正規分布を仮定しなくても、線形不偏推定量の中で最小分散を持つことが保証されています。逆に言えば、外生性が崩れているとき、たとえば第3章で扱った省略変数バイアスがあるときや、目的変数と説明変数が互いに影響し合っているとき、説明変数に測定誤差が含まれているときには、係数の点推定そのものがバイアスを持ちます。この外生性は残差プロットからは直接見えないため、データがどう作られたかという理解で担保するほかない前提です。正規性が主に効いてくるのは、係数のt検定やF検定、信頼区間の構築といった「推論」の段階、それも特にサンプルサイズが小さい場合です。サンプルサイズが十分に大きければ、中心極限定理の働きにより、正規性から多少外れていても検定統計量の分布は近似的に正規分布に従うため、実務上の影響は小さくなる傾向があります。

前提が崩れたときに失われるものを整理すると、線形性の違反は係数の点推定そのものをゆがめる最も重い問題であり、独立性と等分散性の違反は、点推定は保ったまま標準誤差とp値、信頼区間の信頼性を損ないます。正規性の違反は、主に小標本での検定・区間推定の精度に影響するものであり、OLSの点推定そのものに必須の前提ではなく、大標本ではその影響は相対的に小さくなります。「前提が崩れている」とひとくくりにせず、どの前提が崩れ、結果としてどの数字が信用できなくなるのかを切り分けて考えることが、診断の出発点になります。

残差プロットの読み方：残差 対 予測値

前提が満たされているかどうかを確認する最も基本的な道具が「残差プロット」です。横軸に予測値（モデルが計算したあてはめ値）、縦軸に残差（実測値から予測値を引いた差）をとった散布図で、モデルを当てはめたら真っ先に確認すべきグラフだといえます。

前提がおおむね満たされているとき、残差プロットは0の水平線を中心にして、特定のパターンを持たず、ランダムに散らばった雲のような形になります。逆に、この雲の形に何らかの規則性が見えたときは、前提のどこかが崩れているサインです。代表的なパターンを整理します。

扇形のパターン（不均一分散）は、売上や所得のように値の大きさに応じて絶対的なばらつきも大きくなるデータでよく見られます。例えば店舗の規模が大きくなるほど、売上予測の誤差も絶対額として大きくなりやすい、といった状況です。曲線的なパターン（非線形性）は、本来は逓減的な効果（投入量を増やすほど追加の効果が小さくなる関係など）を、直線のモデルで無理に近似しようとしたときによく現れます。うねりのパターン（自己相関）は、時系列データを扱う回帰で特に注意が必要で、この点は第12章で時系列データ特有の回帰の落とし穴としてあらためて扱います。

残差プロットは、時系列順や、モデルに含めていない別の変数（地域、店舗規模、担当者など）を横軸にとっても有用です。予測値に対しては問題がなさそうに見えても、除外した変数に対して残差プロットを描くと明確なパターンが浮かび上がることがあり、これはその変数をモデルに追加すべきという重要な手がかりになります。

Q-Qプロットで正規性を視覚的に確認する

正規性の確認によく使われるのが、Q-Qプロット（Quantile-Quantile Plot）です。これは、理論的な正規分布の分位点（パーセンタイル）を横軸に、残差の同じ分位点を縦軸にとってプロットしたもので、残差が正規分布に近ければ、点はおおむね45度の直線に沿って並びます。

Q-Qプロットの読み方としては、次のような崩れ方に注目します。

• 両端（裾）が直線から上下に大きく逸脱している場合、正規分布よりも裾が重い（極端な値が出やすい）分布であることを示唆します
• 点全体が緩やかなS字を描いている場合、分布が左右非対称（歪んでいる）であることを示唆します
• 中央付近はほぼ直線に乗っているが両端だけずれている場合、大部分のデータは正規性からそれほど外れていないものの、一部に外れ値的な観測が含まれている可能性があります

ここで繰り返し確認しておきたいのは、Q-Qプロットで多少のずれが見つかったからといって、直ちにモデルが使い物にならないと判断する必要はない、という点です。正規性はあくまで検定や信頼区間の精度に関わる前提であり、サンプルサイズが数百、数千という規模であれば、中心極限定理の働きによって検定統計量の分布への影響は限定的になることが多いためです。一方で、サンプルサイズが数十件程度と小さい場合や、裾の重さが極端な場合には、通常のt検定やF検定の前提が崩れている可能性を踏まえ、後述する頑健な手法での確認を検討する価値があります。

外れ値・てこ比・影響点：クックの距離で見分ける

残差診断でしばしば混同されるのが、「外れ値」と「影響点」という2つの概念です。この2つは似ているようで、意味するところがまったく異なります。

• 外れ値：モデルによる予測値から大きく外れた実測値を持つ観測。残差が大きい点を指します
• てこ比（leverage）：説明変数の値そのものが、他の観測から離れた極端な位置にある観測。目的変数の値とは無関係に、説明変数の分布の中で「端」に位置しているかどうかで決まります
• 影響点（influential point）：その観測を取り除くと、回帰係数の推定値が大きく変わってしまう観測を指します

ここで重要なのは、残差が大きい（外れ値である）ことと、モデルに大きな影響を与える（影響点である）ことは、必ずしも一致しないという点です。ある観測が説明変数の分布の中心付近にあり、かつ残差が大きいとしても、その1点を取り除いて回帰をやり直しても係数はほとんど変わらない、ということがあります。逆に、説明変数の値が極端（てこ比が高い）で、なおかつその点がおおむね回帰直線の延長線上に乗っている場合、残差自体は小さいにもかかわらず、その1点の有無で係数の傾きが大きく変わってしまう、ということも起こり得ます。つまり影響点になりやすいのは「てこ比が高く、なおかつ残差もある程度大きい」観測であり、この2つの要素を同時に考慮する指標が必要になります。

この2つの要素（残差の大きさとてこ比の高さ）を1つの数値にまとめた指標が「クックの距離（Cook’s Distance）」です。クックの距離は、ある観測を除いてモデルを当て直したときに、全体の予測値がどれだけ変化するかを要約した指標で、値が大きいほど、その観測が係数の推定に強い影響を与えていることを意味します。経験的な目安として\( 4/n \)（nはサンプルサイズ）を超える観測、あるいは他の観測に比べて突出して大きい値を持つ観測を、まず確認すべき候補として扱うことが多いです。

クックの距離が大きい観測を見つけたとき、機械的にその行を削除してモデルを当て直す、という対応は避けるべきです。まず確認すべきは、その観測がなぜ極端な値になっているのかという原因です。入力ミスやシステム上のデータ不整合であれば修正または除外が妥当ですが、実際に起きた特殊な事象（大型商談、突発的な需要変動、災害などの外部要因）を反映した正当なデータである場合、それを取り除くことはむしろ実態を歪め、都合の良い結果だけを残す操作になりかねません。影響点への対応は、統計的な処理である以前に、データの生成過程を理解する調査から始める必要があります。

前提が崩れているときの対処法を整理する

残差診断によって何らかの前提の崩れが確認できた場合、いくつかの対処法があります。問題の種類によって、有効な手当ては異なります。

変数変換：非線形なパターンや不均一分散が見られる場合、目的変数や説明変数に対数変換や平方根変換を施すことで、関係を線形に近づけたり、分散を安定させたりできることがあります。売上や所得、来店客数のように、値が大きくなるほど絶対的なばらつきも大きくなる（相対的なばらつきはむしろ一定に近い）データでは、対数変換によって不均一分散が緩和されるケースがよく見られます。ただし変換後は係数の解釈も変わるため（対数変換した目的変数の係数は、元のスケールでは相対的な変化率に近い意味を持つようになるなど）、変換の意味を正しく理解したうえで使う必要があります。

頑健標準誤差（HC系）：不均一分散が疑われる一方で、係数の点推定自体は妥当だと考えられる場合、係数はそのままに、標準誤差の計算方法だけを不均一分散に頑健なものに切り替えるという方法があります。これが「頑健標準誤差」で、統計ソフトウェアの表記ではHC0からHC3といった記号で呼ばれます。HCはheteroscedasticity-consistent、つまり不均一分散があっても整合的に計算できる、という意味の略語です。HC3は小標本でも保守的な（広めの）標準誤差を返す設計になっており、実務ではまずHC3から試すことが多い選択肢です。モデルの構造そのものを変えずに、p値や信頼区間の信頼性だけを底上げできる、手軽で有効な対処法だといえます。

頑健回帰：外れ値の影響を受けにくい形で係数そのものを推定し直したい場合には、Huber損失など、通常の二乗誤差よりも外れ値への感度を抑えた損失関数を使う「頑健回帰」という選択肢があります。通常のOLSは誤差を2乗するため極端な値の影響を強く受けますが、Huber損失は一定以上の誤差については2乗ではなく絶対値に近い扱いに切り替えることで、外れ値の影響を抑えます。外れ値を都合よく取り除くのではなく、外れ値の影響を統計的に緩和したモデルで係数を求め直す、という考え方です。

モデルの作り直し：変数変換や頑健な手法でも対応しきれないほど根本的に前提が崩れている場合、線形回帰という枠組みそのものを見直す必要があります。目的変数が0か1の二値である、件数のようにカウントデータである、値が必ず正であるなど、目的変数の性質そのものが正規分布を前提とする通常の線形回帰になじまない場合、一般化線形モデル（GLM）という、より広い枠組みで目的変数の性質に応じた分布とリンク関数を選び直すアプローチが有効です。GLMの考え方については第7章であらためて詳しく扱います。

診断の実務フロー

回帰モデルの前提確認は、一度きりの作業ではなく、繰り返しのプロセスとして組み込むのが実務的です。おおまかな流れは次のようになります。

1. モデルを当てはめる。まずは素直にOLSで回帰を実行し、係数・p値・決定係数を確認します
2. 残差を見る。残差vs予測値プロット、Q-Qプロット、クックの距離を確認し、前提の崩れや影響点の有無を点検します
3. 手当てする。見つかった問題の種類に応じて、変数変換、頑健標準誤差、頑健回帰、モデルの作り直しのいずれかを適用します
4. 再度見る。手当てを施した後のモデルについて、再び残差プロットとQ-Qプロットを確認し、問題が解消されたか、あるいは新たな問題が生じていないかを点検します

このサイクルは、1回で完結するとは限りません。変数変換によって非線形性は解消されても、別の説明変数との関係で新たな不均一分散が見えてくることもあります。診断と手当てを反復しながら、モデルの妥当性を少しずつ高めていくという姿勢が求められます。

回帰モデルを当てはめた後に確認すべき項目を、チェックリストとして整理します。

• 残差vs予測値プロットに扇形やカーブ、うねりといったパターンが見られないか
• Q-Qプロットで残差が45度線から大きく外れていないか（特に裾の部分）
• クックの距離が突出して大きい観測がないか、あればその原因を調査したか
• てこ比の高い観測が、係数の推定に不釣り合いな影響を与えていないか
• 時系列データの場合、残差に自己相関のパターンが出ていないか
• 前提の崩れが見つかった場合、その種類に応じた適切な手当てを施したか
• 手当て後に、再度残差プロットとQ-Qプロットで確認したか

残差診断で前提の崩れを見つけ、手当てを施しても、まだ検討すべき論点は残ります。とりわけ説明変数どうしの相関が強い場合、係数の推定はそれ自体が別の形で不安定になります。次章では、この多重共線性の問題を取り上げます。

Pythonで確認する

ここまでの内容を、statsmodelsを使って一通り確認します。まずは回帰モデルを当てはめたうえで、残差vs予測値プロットとQ-Qプロットを描きます。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# df: 実務データを想定した例(説明変数x1, x2、目的変数y)。手元のデータに置き換えて使うテンプレートです
X = sm.add_constant(df[["x1", "x2"]])
model = sm.OLS(df["y"], X).fit()
print(model.summary())

fitted = model.fittedvalues
resid = model.resid

# 残差 vs 予測値プロット(線形性・等分散性の確認)
plt.figure(figsize=(6, 4))
plt.scatter(fitted, resid, alpha=0.6)
plt.axhline(0, color="red", linestyle="--")
plt.xlabel("予測値")
plt.ylabel("残差")
plt.title("残差 vs 予測値プロット")
plt.show()

# Q-Qプロット(正規性の確認)
sm.qqplot(resid, line="45", fit=True)
plt.title("残差のQ-Qプロット")
plt.show()

残差vs予測値プロットで点が0の水平線を中心に特定のパターンなくばらついていれば、線形性・等分散性についてはひとまず問題が小さいと判断できます。Q-Qプロットで点がおおむね45度線に沿っていれば、正規性についても大きな逸脱はないと考えられます。

続いて、クックの距離を使って影響点の候補を確認します。

from statsmodels.stats.outliers_influence import OLSInfluence
from statsmodels.graphics.regressionplots import influence_plot

influence = OLSInfluence(model)
cooks_d, _ = influence.cooks_distance

n = len(df)
threshold = 4 / n  # 目安となる閾値(4/nを超えると要確認とされることが多い)
suspects = np.where(cooks_d > threshold)[0]
print(f"クックの距離が閾値{threshold:.4f}を超えるサンプル数: {len(suspects)}")

check = df.iloc[suspects].assign(cooks_d=cooks_d[suspects])
print(check.sort_values("cooks_d", ascending=False).head())

# 残差・てこ比・クックの距離を1枚で俯瞰する
fig, ax = plt.subplots(figsize=(7, 5))
influence_plot(model, ax=ax, criterion="cooks")
plt.show()

クックの距離が閾値を超える観測が見つかった場合は、該当する行の元データを確認し、なぜその値になっているのかを調べることが次のステップになります。influence_plotは、横軸にてこ比、縦軸に標準化残差をとり、点の大きさでクックの距離の大小を表現するため、外れ値とてこ比、影響の大きさを1枚のグラフで俯瞰するのに便利です。

最後に、不均一分散が疑われる場合の頑健標準誤差の使い方を確認します。

model_hc3 = sm.OLS(df["y"], X).fit(cov_type="HC3")
print(model_hc3.summary())

# 通常の標準誤差とHC3の標準誤差を比較する
comparison = pd.DataFrame({
    "coef": model.params,
    "se_default": model.bse,
    "se_HC3": model_hc3.bse,
})
comparison["se_ratio"] = comparison["se_HC3"] / comparison["se_default"]
print(comparison)

cov_type引数に”HC3″を指定するだけで、係数の推定値はそのままに、標準誤差だけが不均一分散に頑健な計算方式に切り替わります。通常の標準誤差とHC3の標準誤差を比較し、両者に大きな差がある場合は、通常のOLSのp値や信頼区間をそのまま信用せず、HC3側の結果を優先して判断するのが安全です。

この章を深めたい方への参考書籍

『統計モデルと推測』（松井秀俊・小泉和之、講談社）：回帰モデルにおける推測の考え方を、最尤法や検定の枠組みと合わせて丁寧に解説した一冊です。本章で整理した「点推定は保たれても、標準誤差や検定の精度は前提に依存する」という論点を、より数理的な視点から深く理解したい読者に適しています。

第5章 多重共線性とその対処

前章では、回帰分析が成立するための前提条件を確認し、残差診断によってモデルの粗さを見抜く方法を扱いました。前提が崩れていないかを確かめる作業は、モデルの妥当性を確認するための基本的な点検です。本章で扱う多重共線性は、この点検だけでは見つけにくい問題です。

多重共線性の扱いが難しいのは、モデルの当てはまり自体は悪くならないことが多い点にあります。決定係数は高いまま、残差もきれいに見えるのに、個々の係数の値だけが信用できなくなるという現象が起きます。第3章で「説明のための回帰と予測のための回帰は目的が異なる」という整理をしましたが、多重共線性はこの区別が最も重要になる論点だと考えています。本章では、多重共線性とは何か、何が起きるのか、どう見つけるのか、なぜ起きるのか、そしてどう対処するのかという順で整理します。最後に「交差検証さえしていれば共線性は気にしなくてよい」という誤解についても説明します。

多重共線性とは何か：説明変数同士が似ている状態

重回帰分析は、複数の説明変数を使って目的変数を説明するモデルです。この枠組みが前提としているのは、それぞれの説明変数が目的変数に対して独立した情報を持っているということです。ところが実務のデータでは、説明変数同士が強く相関し合っている場面が頻繁に発生します。この状態を多重共線性（Multicollinearity）と呼びます。

多重共線性には程度の差があります。1つは完全な共線性で、ある説明変数が他の説明変数の完全な線形結合で表せてしまう状態です。たとえば「税抜価格」「消費税額」「税込価格」の3つを同時に説明変数へ入れてしまうと、税込価格は残り2つの単純な合計であるため、最小二乗法の計算そのものが数学的に破綻します。統計ソフトはエラーを返すか、いずれかの変数を自動的に落とすかたちで対応します。

もう1つは強い共線性で、変数同士が完全な線形関係にはないものの、相関が非常に強い状態です。こちらは計算自体は問題なく実行できてしまうため、モデルを作成した本人が異常に気づかないまま結果を報告してしまうことが少なくありません。実務で問題になるのは、ほとんどの場合この強い共線性の方です。本章で扱う「対処が必要な多重共線性」も、以降はこちらを指します。

多重共線性で起きること：標準誤差の膨張と符号の反転

強い共線性がある状態で回帰モデルを推定すると、主に次の3つの現象が起こります。

1つ目は、係数の標準誤差が大きく膨らむことです。標準誤差はその係数の推定値がどれくらいの不確実性を伴っているかを表す指標で、これが大きくなると、係数のp値も大きくなり、本来は目的変数に影響を与えている変数であっても「統計的に有意ではない」という結果になりやすくなります。相関の強い変数同士が目的変数への影響を奪い合ってしまい、それぞれの変数に固有の効果を切り分けられなくなるためです。

2つ目は、係数の符号が理論や常識と逆転することです。たとえば住宅価格を予測するモデルで、床面積と部屋数がともに説明変数に含まれているとします。両者は強く相関しているため、床面積の係数がプラス（広いほど高い）になる一方で、部屋数の係数がマイナス（部屋数が多いほど安い）になるといった、直感に反する結果が出ることがあります。これは部屋数が価格を下げているのではなく、床面積という共通の要因を2つの変数で分け合った結果、計算上の帳尻合わせとして符号が歪んでいるだけです。

3つ目は、データをわずかに追加・変更しただけで係数が大きく動くことです。共線性が強い状態では、各変数への影響の配分が不安定な均衡の上に成り立っており、少数のデータ点が変わるだけで配分がまったく別の組み合わせに移ってしまいます。同じ母集団から取った標本のはずなのに、月によって係数の大きさや符号が入れ替わるようであれば、多重共線性を疑う十分な根拠になります。

ここで重要なのは、これら3つの現象が起きている一方で、モデル全体の予測精度、すなわち決定係数やテストデータに対するRMSE（予測値と実測値の差を二乗して平均し、平方根を取った予測誤差の指標。値が小さいほど予測が当たっている）はそれほど悪化しないことが多いという事実です。目的変数を予測するという観点では、床面積と部屋数のどちらにどれだけの重みを割り振っても、両者を合わせた予測値自体はさほど変わらないためです。この「説明は壊れているが予測はそこそこ機能する」という非対称性こそが、多重共線性という問題の本質です。第3章で整理したとおり、モデルの目的が「Aを1単位増やすとBがどれだけ変わるか」を語る説明のための回帰であれば深刻な欠陥ですが、目的変数の値そのものを当てにいく予測のための回帰であれば、致命傷にならないことも珍しくありません。

見つけ方：相関行列の限界とVIF

多重共線性を疑ったとき、まず思いつくのが説明変数同士の相関行列を確認することです。これは有効な第一歩ですが、それだけでは不十分です。相関行列が捉えられるのは2変数間の関係だけであり、3つ以上の変数が組み合わさって初めて強い共線性が生まれるケースを見逃してしまうためです。たとえば変数A、B、Cのどの2つを取っても相関係数がさほど高くないのに、Aが「BとCの合計に近い」という関係を持っている場合、2変数間の相関だけを見ていては共線性の存在に気づけません。

そこで使われるのが分散拡大要因（Variance Inflation Factor、VIF）です。VIFは、ある説明変数\( X_j \)を目的変数の代わりに置き、残りのすべての説明変数で\( X_j \)自身を説明する補助的な回帰を行い、そこから得られる決定係数\( R_j^2 \)を使って次のように定義されます。

\( \text{VIF}_j = \dfrac{1}{1 – R_j^2} \)

\( R_j^2 \)は「\( X_j \)が他の説明変数によってどれだけ説明されてしまうか」を表します。他の変数から強く予測できてしまう変数ほど\( R_j^2 \)は1に近づき、分母である\( 1 – R_j^2 \)は0に近づくため、VIFは急激に大きくなります。逆に他の変数と無関係な変数は\( R_j^2 \)がほぼ0となり、VIFは1に近い値をとります。VIFは説明変数ごとに1つずつ計算され、目的変数がいくつあっても目的変数自体は計算に関与しない点に注意してください。あくまで説明変数同士の関係だけを見る指標です。

VIFの目安としてよく引用されるのが10という値で、VIFが10を超える変数は強い共線性の疑いが濃いとされます。より厳しい基準を採用する分野では5を境界に使うこともあります。ただしこれらの数値はあくまで経験的な目安であり、絶対的な合否ラインではありません。VIFが8だから安全、11だから危険と機械的に判断するのではなく、VIFの大きさと、係数の符号や有意性の不自然さ、データを変えたときの安定性といった他の兆候をあわせて総合的に判断する姿勢が重要です。分野によっては、サンプルサイズが十分に大きい場合はVIFが多少高くても係数の標準誤差への影響が相対的に小さくなることもあり、目安の数値だけを一人歩きさせないよう注意が必要です。

import pandas as pd
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm

# X_train は説明変数のみを含むDataFrame(目的変数は含めない)。手元のデータに置き換えて使うテンプレートです
X_const = sm.add_constant(X_train)

vif_data = pd.DataFrame()
# 定数項(const)の列はVIFの対象外なので1列目を除いて計算する
vif_data["変数"] = X_const.columns[1:]
vif_data["VIF"] = [
    variance_inflation_factor(X_const.values, i)
    for i in range(1, X_const.shape[1])
]
print(vif_data.sort_values("VIF", ascending=False))

ここでの注意点は、VIFを計算する際に定数項（切片）を加えたうえで計算し、定数項自体のVIFは無視するという手順です。statsmodelsのvariance_inflation_factor関数は、切片を含めずに計算すると値が不自然に歪むことが知られているため、add_constantで定数項を追加した行列を渡し、range(1, ...)として定数項の列（0列目）を計算対象から外すのが正しい使い方です。この手順を誤ると、VIFの値そのものが信頼できなくなるため注意してください。

実務データに共線性が生まれる構造的な理由

多重共線性は分析者の不注意だけで生まれるわけではありません。多くの場合、実務データが持つ構造そのものに原因があります。典型的なのは、複数の変数が同じ背景要因によって同時に動いてしまうケースです。

マーケティング分析では、広告費とキャンペーンの実施回数がよく似た動きをします。予算が潤沢な月ほど広告費も増え、同時にキャンペーンの実施回数も増える傾向があるため、両者は強く相関しやすくなります。小売の需要予測では、気温とアイスクリームの出荷量が並行して動くのと同様に、気温、湿度、来店客数といった季節性に紐づく変数群がまとめて同じ方向に振れます。企業の財務分析では、売上高、従業員数、資産規模のように「企業の規模」という共通の背景因子に紐づく変数が軒並み相関することも珍しくありません。

これらに共通するのは、変数同士が直接の因果関係を持っているのではなく、観測されていない共通の要因（予算規模、季節、企業規模など）によって同時に動かされているという構造です。データを収集する側からすれば自然に発生する相関であり、分析者が意図的に似た変数を選んでしまったというより、実務データの性質そのものが共線性を生み出しやすいと理解しておく必要があります。

共線性のあるデータで係数が暴れる様子を確認する

実際に共線性を持つ合成データを作り、係数がどれだけ不安定になるかを確認します。真の関係式を\( y = 3 x_1 + 2 x_2 + \varepsilon \)と設定したうえで、説明変数\( x_1 \)と\( x_2 \)の相関を強くしたケースと、ほぼ無相関にしたケースを比較します。

import numpy as np
import statsmodels.api as sm

rng = np.random.default_rng(42)
n = 200

def make_data(corr, n=n, rng=rng):
    # x1とx2の相関をcorrで指定した2変量正規分布から生成
    mean = [0, 0]
    cov = [[1, corr], [corr, 1]]
    x1, x2 = rng.multivariate_normal(mean, cov, size=n).T
    # 真の関係式: y = 3*x1 + 2*x2 + ノイズ
    y = 3 * x1 + 2 * x2 + rng.normal(0, 1, n)
    X = sm.add_constant(np.column_stack([x1, x2]))
    return X, y

for corr in [0.0, 0.7, 0.95, 0.99]:
    coefs = []
    # 同じ相関構造から30回サンプリングし、係数の暴れ幅を確認
    for trial in range(30):
        X, y = make_data(corr, rng=np.random.default_rng(trial))
        model = sm.OLS(y, X).fit()
        coefs.append(model.params[1:])  # x1, x2の係数のみ
    coefs = np.array(coefs)
    print(f"相関={corr}: x1係数の標準偏差={coefs[:,0].std():.3f}, "
          f"x2係数の標準偏差={coefs[:,1].std():.3f}")

このコードを実行すると、x1とx2の相関を0から0.99まで引き上げるにつれて、係数（x1、x2それぞれ）の標本間でのばらつき（標準偏差）が急激に増大する様子が確認できます。相関が0の場合、30回の試行を通じてx1の係数は2.85から3.15、x2の係数は1.84から2.12という狭い範囲に収まり、いずれも真の値（3と2）のすぐ近くで安定しています。ところが相関を0.99まで上げると、同じ生成過程からサンプリングしているにもかかわらず、x1の係数は2.18から4.12、x2の係数は0.90から2.87まで散らばり、標準偏差は6倍に膨らみます。x2にいたっては、真の値2の半分以下に沈む試行と、1.4倍を超える試行が同居することになります。サンプルサイズをさらに小さくすれば、係数の符号が反転する試行すら現れます。真の関係式は一切変えていないにもかかわらず、係数の推定値だけがこれほど不安定になるという事実は、多重共線性が「モデルの誤り」ではなく「情報を切り分けられないことによる推定の不確実性」であることを端的に示しています。

対処法の一覧：目的に応じた選び分け

多重共線性への対処法は1つに絞られるものではなく、分析の目的やデータの性質に応じて選び分けるべきものです。代表的な対処法を次の表に整理します。

中心化についてはやや説明を補います。交互作用項（\( x_1 \times x_2 \)）や2乗項（\( x_1^2 \)）を説明変数に加えると、これらの項は元の変数と強く相関しやすく、共線性を生み出します。ところが、この共線性は変数の値そのものの原点（0の位置）をどこに置くかという計算上の理由で生じている部分が大きく、各変数から平均値を引いて中心化してから交互作用項・多項式項を作ることで、この種の共線性はかなり緩和されます。ただし中心化が効くのはこうした構造的な共線性に限られ、広告費とキャンペーン回数のように、変数そのものが実質的に同じ情報を運んでいる場合には、中心化をしても相関の強さ自体は変わらないため効果は限定的です。

多重共線性への対処は「症状を消すこと」が目的ではありません。VIFの数値を基準値以下に下げること自体を目標にしてしまうと、目的変数を説明するうえで本来必要な変数まで機械的に削ってしまう恐れがあります。まず「このモデルは説明のためか、予測のためか」を明確にし、その目的に照らして、どの対処法がどんな代償を伴うかを理解したうえで選ぶという順序が欠かせません。

VIFを計算し、対処前後で係数の安定性を比較する

変数を絞る対処法を実際に試し、VIFと係数の変化を確認します。ここでは、共線性のある\( x_1 \)、\( x_2 \)に加えて、これらと無関係な\( x_3 \)を含むデータセットを用意します。

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

rng = np.random.default_rng(0)
n = 300

# x1とx2は強い相関(0.97)、x3は独立
mean = [0, 0]
cov = [[1, 0.97], [0.97, 1]]
x1, x2 = rng.multivariate_normal(mean, cov, size=n).T
x3 = rng.normal(0, 1, n)
y = 3 * x1 + 2 * x2 + 1.5 * x3 + rng.normal(0, 1, n)

df = pd.DataFrame({"x1": x1, "x2": x2, "x3": x3})

def show_vif(X):
    Xc = sm.add_constant(X)
    vif = pd.DataFrame({
        "変数": Xc.columns[1:],
        "VIF": [variance_inflation_factor(Xc.values, i)
                for i in range(1, Xc.shape[1])]
    })
    return vif

# 対処前: x1, x2, x3すべてを投入
X_full = df[["x1", "x2", "x3"]]
model_full = sm.OLS(y, sm.add_constant(X_full)).fit()
print("対処前のVIF")
print(show_vif(X_full))
print(model_full.params)

# 対処後: 強く相関するx2を除外し、x1を代表変数として残す
X_reduced = df[["x1", "x3"]]
model_reduced = sm.OLS(y, sm.add_constant(X_reduced)).fit()
print("対処後(x2を除外)のVIF")
print(show_vif(X_reduced))
print(model_reduced.params)

対処前の結果では、x1とx2のVIFがいずれも10を大きく上回る一方、無関係なx3のVIFは1に近い値にとどまります。x1とx2の係数の標準誤差も大きく、符号や大きさが真の値（3と2）から離れた推定になりやすい状態です。x2を除外した対処後のモデルでは、残ったx1の係数はx1とx2の効果を合わせて引き受けるかたちに変化し、標準誤差も小さくなって安定します。ただしこの場合、x1の係数はもはや「x2を固定したうえでのx1固有の効果」ではなく、「x1とx2が共有する成分全体の効果」に近い意味へと変わっている点に注意が必要です。変数を削るという操作は、単に不要な変数を取り除く作業ではなく、残った変数の意味そのものを変えてしまう操作であることを理解しておく必要があります。

「交差検証さえしていれば共線性は無視してよい」という誤解

ここまでの内容を踏まえると、「多重共線性があっても交差検証で予測精度を確認しておけば問題ない」という考え方が誤解を含んでいることが分かります。交差検証は、モデルが未知のデータに対してどれだけ正確に予測できるかを検証する手続きであり、第4章や本章で見てきたとおり、多重共線性があっても予測精度自体はそれほど悪化しないことが多いため、交差検証のスコアだけを見れば「問題ない」という結果が出ることは十分にあり得ます。

しかし交差検証は、個々の係数の値が安定しているか、符号が理論と整合しているかを検証する手続きではありません。予測というゴールに対しては交差検証で十分な保証が得られても、「x1が1単位増えるとyがどれだけ変わるか」という係数の解釈というゴールに対しては、何の保証にもなりません。交差検証のスコアが良好だからといって係数の符号反転や標準誤差の膨張が解消されているわけではなく、両者は別の問題を検証する別の手続きだと理解しておく必要があります。

したがって、多重共線性への向き合い方は、モデルの目的によって次のように整理できます。予測を主目的とし、係数を説明資料や意思決定の根拠として使う予定がないのであれば、交差検証によって予測性能さえ確認できていれば、多重共線性を過度に気にする必要は薄くなります。一方で、施策担当者に「この変数を増やせば効果がどれだけ見込めるか」を説明する、あるいは係数そのものを経営判断の根拠として提示するのであれば、交差検証のスコアがどれだけ良くても、VIFを確認し、必要な対処を講じないまま係数を報告することは避けるべきです。目的を明確にしないまま「交差検証を通ったから大丈夫」という判断だけで進めることが、この誤解の核心にある落とし穴です。

この章を深めたい方への参考書籍

『多変量解析入門 線形から非線形へ』（小西貞則、岩波書店）は、重回帰分析から一般化線形モデル、さらに非線形回帰までを一貫した視点で扱っており、本章で扱ったVIFや多重共線性の数理的な背景を、より厳密な行列表現に基づいて理解し直したい読者に向いています。実務での対処法だけでなく、なぜそうした現象が数式のうえで起きるのかを掘り下げたい場合の副読本としておすすめできます。

第6章 情報量規準と交差検証によるモデル選択

前章では、説明変数同士が強く相関する多重共線性が係数の推定を不安定にする仕組みと、VIFによる検知や変数整理といった対処法を扱いました。多重共線性への対応も含め、重回帰モデルを組み立てる過程では、常に「どの説明変数を入れるか、入れないか」という選択がついて回ります。年齢、性別、地域、購買履歴、Web行動ログなど候補となる変数が10個あれば、組み合わせだけで1000通りを超えるモデルが作れてしまいます。本章では、この膨大な候補から1つ、あるいは少数のモデルに絞り込んでいくモデル選択の考え方を整理します。

モデル選択の基本：当てはまりとシンプルさのバランス

説明変数の候補がp個あるとき、定数項だけのモデルから全変数を含むモデルまで、理論上は\( 2^p \)通りのモデルが候補になり得ます。変数が10個であれば1,024通り、20個であれば100万通りを超えます。実務でこれらすべてを目視で比較することは現実的ではなく、何らかの基準に沿って機械的に絞り込む仕組みが必要になります。

モデル選択の根底にあるのは、当てはまりの良さと、モデルのシンプルさという、しばしば相反する2つの価値観のバランスをどう取るかという問題です。手元のデータへの当てはまりだけを追求するなら、説明変数を増やせば増やすほど有利になります。しかし変数を増やしすぎたモデルは、手元のデータの持つ偶然のノイズまで拾い上げてしまい、新しいデータへの予測力がかえって落ちることがあります。この現象は過学習と呼ばれます。過学習は機械学習の文脈で語られることが多いですが、統計モデリングでも本質的に同じ問題として生じます。手元のデータへの当てはまりと、まだ見ていないデータへの予測力は別物であり、両者を混同しないことがモデル選択の出発点になります。

そこで統計学では古くから、必要以上に複雑なモデルを避け、目的を達成できる範囲でできるだけシンプルなモデルを選ぶべきだ、という考え方が重視されてきました。これをパーシモニーの原則（倹約の原則）と呼びます。同程度の説明力を持つモデルが複数あるなら、変数の数が少なく解釈しやすいモデルを優先する、という発想です。本章で扱う情報量規準や交差検証は、この「当てはまりの良さ」と「シンプルさ」という2つの軸を1つの数値に落とし込んで比較するための道具です。以降では、その代表であるAICとBIC、そして実際にデータを分けて予測性能を測る交差検証を順に取り上げ、最後に説明を目的とした分析でのモデル選択の考え方を整理します。

決定係数で選んではいけない理由

モデル選択の基準として、決定係数\( R^2 \)を使いたくなる場面は少なくありません。しかし第1章で確認したとおり、決定係数には見過ごせない性質があります。説明変数をどれだけ意味のないものであっても追加すると、決定係数は必ず同じか、わずかでも上昇します。これは決定係数が「残差の二乗和をどれだけ減らせたか」という当てはまりの良さだけを測る指標であり、変数を増やすこと自体へのペナルティを一切持たないためです。乱数で生成しただけの、理論的に何の意味もない変数を追加しても、決定係数はほぼ確実に上昇します。したがって決定係数の大小だけを根拠にモデルを比較すると、常に変数を最大限詰め込んだモデルが選ばれてしまい、モデル選択としての役割を果たしません。

この欠点を補うために考案されたのが自由度調整済み決定係数です。変数の数（自由度）を考慮した調整を加えることで、意味のない変数を追加すると値が下がるように設計されています。決定係数がモデルへの当てはまりだけを見る指標であるのに対し、こちらは変数を増やすことへの一定のペナルティを織り込んでいる、と理解しておくとよいでしょう。

ただし自由度調整済み決定係数にも限界があります。ペナルティの強さが控えめであるため、ほとんど説明力を持たない変数を追加しても値がわずかに上昇してしまうことが珍しくありません。そのため実務では、これを「モデルの当てはまりを大まかに把握するための補助指標」と位置づけ、この後で扱う情報量規準や交差検証と組み合わせて使うのが安全な扱い方になります。決定係数だけでモデルの優劣を決め切る使い方は避けるべきです。

AIC（赤池情報量規準）：当てはまりとペナルティのバランス

モデル選択の実務でもっとも広く使われる指標がAIC（赤池情報量規準、Akaike Information Criterion）です。AICは次の式で定義されます。

\( AIC = -2\log L + 2k \)

\( \log L \)はモデルの最大対数尤度、つまりモデルが手元のデータをどれだけうまく説明できているかを表す量で、値が大きいほど当てはまりが良いことを意味します。\( -2\log L \)はその符号を反転させたものなので、値が小さいほど当てはまりが良い、と読み替えられます。\( k \)はモデルに含まれるパラメータの数（説明変数の数に切片や誤差分散などを加えたもの）です。つまりAICは「当てはまりの悪さを表す\( -2\log L \)に、パラメータの数に応じたペナルティ\( 2k \)を上乗せした数値」であり、AICが小さいモデルほど、当てはまりの良さとシンプルさのバランスが優れていると評価されます。

AICの理論的な位置づけは、正確に理解しておく価値があります。AICは「どのモデルが真実か」を判定するための指標ではありません。AICが推定しようとしているのは、そのモデルを未知の新しいデータに適用したときに生じる予測誤差の期待値、より正確には、モデルが表す確率分布と、データを生み出している真の（未知の）確率分布との間の隔たり（カルバック・ライブラー情報量）です。AICは「真のモデルを言い当てる」ことではなく、「候補の中で、未知データへの予測性能がもっとも良いと期待できるモデルはどれか」を推定するための指標だと理解しておくことが重要です。

AICを実務で使ううえで欠かせないのが、比較できる範囲の理解です。AICの絶対値そのものには意味がなく、複数のモデルのAICを比較して相対的な優劣を見る指標です。ただし、この比較が成立するのは、同じデータ、同じ目的変数を対象にしたモデル同士に限られます。目的変数を対数変換したモデルと、変換前のままのモデルでは尤度の計算対象そのものが変わるため、AICを単純に比較できません。また、欠損値の扱いによってモデルごとに実質的に使われているデータの行数が異なる場合も、比較の前提が崩れます。AICで複数モデルを比較する際は、目的変数の定義とデータの行数がすべてのモデルで揃っていることを確認してください。

AICの差の大きさをどう読むかについても、目安が知られています。もっともAICが小さいモデルを基準にしたとき、他のモデルとの差が2未満であれば両者はほぼ同等の支持を受けていると考えられます。差が4から7程度であれば劣っている側への支持はある程度弱いとみなされ、差が10を超えると劣っている側はほとんど支持されないと判断されるのが一般的な目安です。あくまで経験的な指針であり、厳密な検定のような明確な合否の基準ではない点には注意してください。

BICとの違い：ペナルティの強さと目的

AICとしばしば並べて語られる指標にBIC（ベイズ情報量規準、Bayesian Information Criterion）があります。BICは次の式で定義されます。

\( BIC = -2\log L + k\log(n) \)

AICとの違いは、パラメータ数\( k \)にかけるペナルティの大きさにあります。AICのペナルティは常に\( 2k \)で一定なのに対し、BICのペナルティは\( k\log(n) \)、つまりサンプルサイズ\( n \)の対数に比例して大きくなります。サンプルサイズが8を超えると\( \log(n) \)は2を上回るため、実務データの多くでBICはAICよりも複雑なモデルに厳しいペナルティを課します。結果として、AICとBICで異なるモデルを選んだ場合、BICの方がより変数の少ないシンプルなモデルを選ぶ傾向があります。

この違いの背景には、両者が目指しているものの違いがあります。BICは、候補の中に真の（データを生成した）モデルが含まれているという前提のもとで、サンプルサイズを大きくしていくと確率1でその真のモデルを選び出すという性質（一致性）を持つように設計されています。BICは「真のモデルを探し当てる」という思想に基づいた指標です。これに対しAICは、真のモデルが候補に含まれているかどうかにかかわらず、未知データに対する予測誤差の期待値をできるだけ小さくするモデルを選ぶという、予測の良さを重視する思想に基づいています。この思想の違いを押さえておくと、両者の結果が食い違ったときにどちらを重視すべきかの判断がしやすくなります。実務での目的が将来のデータに対する予測精度であれば、AICを主軸に据えるのが自然な選択です。

交差検証によるモデル選択

情報量規準がモデルの尤度とパラメータ数から計算する理論的な指標であるのに対し、交差検証は実際にデータを分割し、モデルの予測性能を直接測定するアプローチです。手元のデータを学習用と検証用に分割し、学習用データだけでモデルを推定したうえで、検証用データに対する予測誤差（回帰であれば二乗誤差の平方根、RMSEなど）を計算します。この分割を複数回、位置を変えながら繰り返し、誤差の平均を取る手法がk分割交差検証で、機械学習の分野で汎化性能を測る標準的な方法として広く使われています。

統計モデリングの文脈における交差検証の位置づけは、機械学習の予測モデル構築とやや異なります。統計モデリングでは係数の解釈や統計的な有意性が重視される場面が多い一方、モデルが実データに対してどれだけの予測力を持つかという観点は情報量規準だけでは直接確認できません。交差検証は、情報量規準が理論的に推定している「未知データへの当てはまりの良さ」を、実際にデータを分けて確認する手段だと位置づけると理解しやすくなります。サンプルサイズが小さい場合や、モデルの分布の仮定に自信が持てない場合には、理論的な近似に頼る情報量規準よりも実測に基づく交差検証の結果を優先したい場面もあります。

AICと交差検証は、まったく無関係の指標ではありません。線形回帰モデルにおいては、1つのデータ点だけを検証用に残し、残り全部を学習に使うことを全データ点について繰り返すLOOCV（1つ抜き交差検証）による評価結果と、AICによる評価結果が、サンプルサイズを大きくしていく極限で漸近的に近い挙動を示すことが知られています。この関係は、AICが「未知データへの予測誤差の期待値を推定する」という交差検証と同じ目的を、計算コストの軽い理論式によって近似的に実現しようとしている指標であることを裏づけています。実務上は、サンプルサイズが十分大きく計算資源に制約がある場合はAICを優先し、サンプルサイズが小さい、あるいは非線形なモデルで尤度の計算が難しい場合には交差検証を使う、という使い分けが現実的です。両方を並べて確認し、結論が一致することを確かめる進め方がもっとも手堅い方法です。

ステップワイズ法の功罪

説明変数の候補が多いとき、AICなどの指標を基準に、変数の追加や除去を自動的に繰り返しながらモデルを絞り込んでいく手法をステップワイズ法と呼びます。何もない状態から変数を1つずつ追加していく前進選択法、全変数を含む状態から不要な変数を1つずつ取り除いていく後退消去法、そして両者を組み合わせ、追加と除去を交互に試しながら進める変数増減法（狭義のステップワイズ法）が代表的な方式です。

ステップワイズ法の便利さは明らかです。数十個の候補変数があっても、人手による総当たりの比較を行うことなく、機械的にAICの改善するモデルへ向かって探索を進めてくれます。総当たりで\( 2^p \)通りを比較するのが現実的でない場面では、探索の範囲を絞り込む道具として一定の価値があります。

一方でステップワイズ法には、実務でしばしば見過ごされがちな危うさがいくつも指摘されています。第一に、多重検定の問題です。ステップワイズ法は変数を追加するかどうかの判断のたびに統計的な有意性の判定を繰り返し行っています。候補変数が多いほど、この繰り返しの中で偶然有意な結果が出てしまう変数が紛れ込む確率が高まり、本来は関係のない変数がモデルに残ってしまうことがあります。第二に、係数のバイアスです。選ばれた変数の係数は、選択の過程自体がデータに合わせて有利な変数を拾い上げているため、真の効果より大きめに推定されがちであることが知られています。第三に、理論なき選択という問題です。ステップワイズ法は変数同士の因果関係や、業務上その変数を含めるべき理由をまったく考慮せず、純粋に統計的な当てはまりの改善だけを基準に判断を進めます。結果として、業務知識の観点からは不自然な変数の組み合わせが選ばれることがあります。

• 候補変数への判定を繰り返すことで生じる多重検定の問題
• 選ばれた変数の係数が真の値より過大に推定されやすいという係数バイアス
• 因果関係や業務上の理由を考慮しない、理論に基づかない変数選択になりやすい
• データが少し変わるだけで、選ばれる変数の組み合わせが大きく変わることがある不安定さ

これらの問題を踏まえると、ステップワイズ法を実務で使う際の位置づけは、最終的な変数の確定手段としてではなく、候補となる変数を大まかに絞り込むための一次スクリーニングの道具として使う、というものが妥当です。提示された候補モデルをそのまま採用するのではなく、選ばれた変数、選ばれなかった変数の両方について、業務知識や理論的な妥当性に照らして人間が最終判断を下す工程を必ず挟むことが、この手法の危うさを抑える実務上の作法になります。

説明目的のモデル選択における統計指標の位置づけ

ここまで扱ってきたAIC、BIC、交差検証は、いずれも予測の良さやモデルの当てはまりを測る指標です。しかし第3章で扱ったとおり、回帰分析には予測を目的とする使い方と、ある要因が結果にどう影響するかを説明することを目的とする使い方の2種類があり、この目的の違いによって、モデル選択で重視すべき基準もまったく異なってきます。

説明を目的としたモデルにおいて、統計指標だけを基準にモデルを選ぶことには根本的な限界があります。ここでは第3章で扱った交絡変数、すなわち注目している要因と結果の両方に影響を与えている第三の変数の議論が関わってきます。ある変数が結果に与える影響を正しく取り出すためには、交絡を生んでいる変数を、たとえその変数自体がAICを大きく改善しなくても、モデルに含めておく必要があります。交絡変数を除いてAICがわずかに改善したとしても、それは「予測にとって不要な変数を削った」結果であって、「注目している効果を正しく推定できるようになった」結果ではありません。逆に、統計的な説明力が乏しい変数でも、業務上あるいは理論上、結果に影響を与えることが分かっている要因であれば、モデルに残しておくべき場面があります。

したがって説明を目的とした回帰分析では、統計指標は変数選択を最終的に決定する主役ではなく、あくまで補助的な位置づけとして扱うべきです。主役となるのは、第3章で整理した因果の見取り図（どの変数がどの変数に影響し、どこに交絡が生じているかという構造の理解）と、対象領域の業務知識です。AICや交差検証が示す「予測に効くかどうか」という情報は変数候補を検討する材料の1つとして参照しつつ、最終的な採否は因果構造と業務知識に照らして判断する、という順序を守ることが誠実なモデル選択の姿勢になります。

予測を目的としたモデル選択では、AIC・BIC・交差検証といった指標に沿って当てはまりと複雑さのバランスを取ることが主軸になります。これに対して説明を目的としたモデル選択では、統計指標はあくまで参考情報にとどめ、因果の見取り図と業務知識に基づいて、交絡変数を含めるべきかどうかを判断することが主軸になります。目的が異なれば、同じAICという数値の使い方も変わってくる、という点を常に意識しておく必要があります。

モデル選択の実務フロー

ここまでの内容を踏まえ、実務でモデル選択に取り組む際の大まかな流れを整理します。指標の計算に入る前に、目的をはっきりさせておくことが、後の判断のぶれを防ぐ鍵になります。

モデル選択の実務フローは、おおむね次の4つの段階を経て進めます。

1. 目的の明確化: このモデルは予測のためか、それとも特定の要因の効果を説明するためかを最初に確認する。目的によって、この後の判断で重視すべき基準が変わる
2. 候補モデルの設計: 業務知識や第3章の因果の見取り図に基づき、含めるべき変数、除いてよい変数の当たりをつけたうえで、比較したい候補モデルをいくつか用意する。ステップワイズ法を使う場合も、ここでの一次スクリーニングにとどめる
3. 指標での比較: AIC・BIC・自由度調整済み決定係数・交差検証によるRMSEなど、複数の指標を並べて候補モデルを比較する。1つの指標だけに頼らず、複数の指標が同じ方向を示しているかを確認する
4. 頑健性チェック: データを時期や属性で分割し、それぞれの部分データでも同じ変数、同じ結論が支持されるかを確認する。特定の期間やセグメントだけに依存した、たまたまの結果でないかを検証する

この4段階のうち、実務で省略されがちなのが最後の頑健性チェックです。統計指標の数値上はもっとも優れて見えるモデルであっても、データを2つに分けて別々に推定したときに、選ばれる変数や係数の符号が大きく変わるようであれば、その結果はデータ全体のたまたまの特徴を拾ってしまっている可能性があります。特に説明を目的としたモデルでは、「この要因は結果に効いている」という結論を業務判断に使うことになるため、データを分けても結論が変わらないかという確認は、最終判断を下す前の重要な一手間になります。

ここからは、Pythonでモデル選択の指標を計算する例を見ていきます。まず、statsmodelsで複数の候補モデルのAICを比較する流れを確認します。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

# 住宅価格を床面積・築年数・駅距離・部屋数から予測する想定データ
rng = np.random.default_rng(0)
n = 300
df = pd.DataFrame({
    "floor_area": rng.normal(70, 20, n),
    "age": rng.integers(0, 40, n),
    "dist_station": rng.normal(10, 5, n).clip(min=0.5),
    "rooms": rng.integers(1, 6, n),
})
noise = rng.normal(0, 30, n)
df["price"] = (
    50 + 3.2 * df["floor_area"] - 1.1 * df["age"]
    - 0.8 * df["dist_station"] + noise
)

# 説明変数の組み合わせが異なる3つの候補モデルを用意する
model_a = smf.ols("price ~ floor_area", data=df).fit()
model_b = smf.ols("price ~ floor_area + age + dist_station", data=df).fit()
model_c = smf.ols(
    "price ~ floor_area + age + dist_station + rooms", data=df
).fit()

for name, m in [("床面積のみ", model_a), ("床面積+築年数+距離", model_b),
                ("全変数", model_c)]:
    print(f"{name}: AIC={m.aic:.1f}  BIC={m.bic:.1f}  "
          f"調整済みR2={m.rsquared_adj:.3f}")

このコードを実行すると、床面積だけの単純なモデルよりも、築年数と駅からの距離を加えたモデルの方がAICもBICも小さくなり（値が小さいほど良い）、当てはまりと予測の両面で優れていることが確認できます。一方で、価格の生成に本来関与していない部屋数まで加えた全変数モデルでは、AICやBICが3変数モデルよりも改善しないことがあります。これは、部屋数という変数が追加のパラメータ数に見合うだけの説明力を持っていないことを、ペナルティが検出している例です。決定係数だけを見ていると、変数を増やすほど数値が上がるため、この判断はできません。

続いて、AIC・BIC・調整済み決定係数に交差検証のRMSEを加えた比較表を作る例を示します。理論的な評価と、実際にデータを分けて測定した予測誤差を同じ表で見比べます。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.model_selection import KFold
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

formulas = {
    "床面積のみ": "price ~ floor_area",
    "床面積+築年数+距離": "price ~ floor_area + age + dist_station",
    "全変数": "price ~ floor_area + age + dist_station + rooms",
}

kf = KFold(n_splits=5, shuffle=True, random_state=0)
results = []

for name, formula in formulas.items():
    fitted = smf.ols(formula, data=df).fit()

    x_cols = [c.strip() for c in formula.split("~")[1].split("+")]
    X = df[x_cols].values
    y = df["price"].values

    fold_rmse = []
    for train_idx, test_idx in kf.split(X):
        lr = LinearRegression().fit(X[train_idx], y[train_idx])
        pred = lr.predict(X[test_idx])
        fold_rmse.append(
            mean_squared_error(y[test_idx], pred) ** 0.5
        )

    results.append({
        "モデル": name,
        "AIC": round(fitted.aic, 1),
        "BIC": round(fitted.bic, 1),
        "調整済みR2": round(fitted.rsquared_adj, 3),
        "CV RMSE(5分割平均)": round(np.mean(fold_rmse), 1),
    })

comparison_table = pd.DataFrame(results)
print(comparison_table.to_string(index=False))

この比較表を眺める際は、AICとBICは値が小さいほど、調整済み決定係数は値が大きいほど、CVによるRMSEは値が小さいほど良いモデルであることを念頭に置いてください。理論的な指標であるAIC・BICと、実測に基づく交差検証のRMSEが同じモデルを支持していれば、その結果への確信を強めてよい根拠になります。逆に両者の結論が食い違う場合は、サンプルサイズの大きさや、モデルの分布の仮定が妥当かどうかを見直す手がかりになります。いずれの場合も、この表の数値だけでモデルを確定させるのではなく、前節の実務フローに沿って、業務知識との整合性や頑健性のチェックを経たうえで最終的な採否を判断することが欠かせません。

この章を深めたい方への参考書籍

『情報量規準』（小西貞則・北川源四郎、朝倉書店）：AICの提唱者の1人である赤池弘次の理論的背景から、BICやその他の情報量規準との数理的な関係までを体系的に扱った専門書です。本章で直感的な説明にとどめたAICの導出過程や、カルバック・ライブラー情報量との関係を、数式のレベルで深く理解したい読者に向いています。統計モデリングを専門的に扱う立場であれば、一度は目を通しておきたい定番書です。

第7章 線形回帰を拡張する一般化線形モデル（GLM）

通常の線形回帰が適さない場面：0/1、件数、割合、金額のデータ

ここまでの章では、目的変数が連続的な値を取り、その誤差がおおむね正規分布に従うという前提のもとで、重回帰モデルを扱ってきました。売上金額や気温、テストの点数のように、理論上は正負どちらの方向にも連続的に動き得る量を予測する場面では、この前提は概ね妥当です。しかし実務で扱いたい目的変数は、必ずしもこの形をしているとは限りません。

典型的には次のようなデータが挙げられます。顧客が解約したかどうか、与信審査に通ったかどうかといった0か1かの二値データ。工場のある工程で1日に発生する不良品の件数や、Webサイトへの問い合わせ件数といった、0以上の整数しか取らないカウントデータ。ある商品カテゴリの中で特定商品が占める割合や、生産ラインにおける不良率といった、0から1の範囲に収まる比率データ。そして、購入金額や保険金の支払額のように、理論上は0以上の値しか取らず、しばしば右に裾を引く分布をする非負の金額データです。なお金額データにガンマ分布を使う場合は、ガンマ分布が0そのものを取らない点に注意が必要です。購入しなかった顧客の0円が混ざるようなデータでは、0かどうかを分ける段階と金額を説明する段階に分けた二部モデルや、0を許容するTweedie分布を検討することになります。

こうしたデータに、通常の線形回帰（最小二乗法による重回帰）をそのまま当てはめると、いくつかの不具合が表面化します。1つ目は、予測値が現実にはあり得ない範囲にはみ出してしまう問題です。解約の有無（0/1）を目的変数にして線形回帰を当てはめると、説明変数の値によっては予測値が1.3や-0.2といった、確率としては意味を持たない値になってしまいます。件数を目的変数にした場合も同様で、マイナスの件数という現実には存在しない予測が平気で出力されます。

2つ目は、分散の構造が目的変数の実態と噛み合わないという問題です。通常の線形回帰は、目的変数の平均がどこにあっても誤差のばらつき（分散）は一定であるという前提（等分散性）を置いています。しかしカウントデータでは、平均的な発生件数が多い条件ほど、件数のばらつきそのものも大きくなるという性質が自然に見られます。1日あたりの問い合わせ件数が平均2件の窓口と平均200件の窓口とでは、日によるブレの大きさもまったく違うはずです。この分散が平均水準によって変化するという実態を、線形回帰の等分散性の前提は正しく表現できません。

3つ目は、係数の解釈が不自然になりやすいという問題です。解約の有無に線形回帰を当てはめる手法は線形確率モデルと呼ばれ、係数を「説明変数が1単位増えると解約確率が何ポイント変化するか」という形でそのまま解釈できる手軽さから実務で使われることもあります。しかし前述の範囲外予測の問題に加えて、確率が0や1に近づくにつれて実際の増減の余地が小さくなっていくはずだ、という確率特有の非線形な振る舞いを表現できない点は、看過できない弱点として残ります。

データの生成のされ方に合わせてモデルを選ぶという発想

この不具合を解消する鍵は、「目的変数がどのようなメカニズムで生成されているデータなのか」という視点に立ち返ることにあります。二値の結果は成功か失敗かの試行を積み重ねたベルヌーイ的な出来事の結果であり、件数データは一定期間内に離散的な出来事が発生する回数であり、金額データはしばしば多数の小さな要因が掛け合わさって生じる非負の連続量です。それぞれの生成メカニズムには、その性質を素直に表現できる確率分布が対応しています。

この発想は、確率分布そのものを扱う統計学の基礎的な単元と直接つながっています。二値の結果は二項分布、単位期間あたりの発生回数はポアソン分布、連続的な誤差の蓄積は正規分布というように、目的変数の性質に応じて自然に対応する確率分布が存在するという考え方です。回帰分析とは、本来「目的変数の分布の平均（期待値）を、説明変数を使ってどう表現するか」という営みであり、その平均を表現する際に、目的変数の実態に合った確率分布を出発点に据えるという発想に立てば、線形回帰は数ある選択肢の中の1つの特殊なケースに過ぎないと分かります。

この「データの生成のされ方に合わせてモデルの前提となる確率分布を選ぶ」という統計モデリングの基本姿勢を、線形予測子とリンク関数という2つの仕組みで体系的に扱えるように拡張した枠組みが、本章で扱う一般化線形モデル（Generalized Linear Model、以下GLM）です。GLMは1970年代にネルダーとウェダバーンによって整理された考え方であり、通常の線形回帰、ロジスティック回帰、ポアソン回帰といった一見バラバラに見える手法を、共通の1つの理論的な枠組みの中に位置づけ直します。

GLMを構成する3つの部品

GLMは、次の3つの部品の組み合わせとして定義されます。この3つを目的変数の性質に応じて使い分けることで、線形回帰では扱えなかった目的変数にも同じ枠組みで対応できるようになります。

• 確率分布（ランダム成分）: 目的変数のばらつきの型を表す部品です。目的変数\( y \)が、どのような確率分布に従っていると仮定するかを指定します。正規分布、二項分布、ポアソン分布、ガンマ分布などが、目的変数の性質に応じて選ばれます。この分布族は「指数型分布族」と呼ばれるグループに属している必要があります。指数型分布族とは、確率の計算式を指数関数を使った共通の書式に書き直せる分布の集まりで、正規分布・二項分布・ポアソン分布・ガンマ分布はいずれもこれに含まれます。書式がそろっているおかげで、分布が違っても同じ推定の手続きでパラメータを求められる、というのがGLMがこの条件を置いている理由です。
• 線形予測子（システマティック成分）: 説明変数の重み付き和のことで、通常の重回帰式\( \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p \)とまったく同じ形をしています。この\( \eta \)（イータ）を線形予測子と呼びます。GLMにおいても、説明変数を組み合わせる部分は線形回帰と変わらず線形の式のままである、という点が重要です。
• リンク関数: 線形予測子\( \eta \)と、目的変数の期待値\( \mu \)（ミュー、目的変数の平均）とをつなぐ変換式です。\( g(\mu) = \eta \)という形で表され、線形予測子はどのような値も取り得る一方、目的変数の期待値には（確率なら0から1、件数なら0以上というように）取り得る範囲の制約があるため、この値域の違いを吸収する変換としてリンク関数を用います。

この3つの部品のうち、線形予測子は常に説明変数の重み付き和のままであり、変化するのは確率分布とリンク関数の組み合わせだけだという点が、GLMという枠組みの見通しの良さを支えています。目的変数がどのような性質を持っていても、「説明変数を線形に組み合わせる」という重回帰の基本発想そのものは崩さずに済み、変えるべきは分布とリンク関数の選び方だけである、と捉え直すことができます。

正規線形モデルはGLMの特殊ケース

この枠組みに立つと、これまでの章で扱ってきた通常の重回帰（正規線形モデル）が、実はGLMの中でも最も単純な特殊ケースであることが分かります。通常の重回帰は、確率分布に正規分布を、リンク関数に恒等関数（何も変換しない、\( g(\mu) = \mu \)という関数）を選んだ場合のGLMそのものです。線形予測子\( \eta \)がそのまま目的変数の期待値\( \mu \)になるため、見慣れた\( y = \beta_0 + \beta_1 x_1 + \cdots + \varepsilon \)という式の形に一致します。

この見方の意義は、単なる理論上の整理にとどまりません。「通常の線形回帰とロジスティック回帰、ポアソン回帰はまったく別の手法である」という捉え方から、「同じGLMという枠組みの中で、確率分布とリンク関数の指定だけが異なる、地続きの手法群である」という捉え方に変わることで、目的変数の性質が変わるたびに手法をゼロから学び直す必要がなくなります。線形予測子の作り方、係数の考え方、多重共線性への配慮、モデル選択の発想といった、これまでの章で積み上げてきた知識の大部分は、GLM全体に共通する基礎知識としてそのまま活かせます。

GLMは「確率分布」「線形予測子」「リンク関数」という3つの部品の組み合わせであり、通常の線形回帰は「正規分布+恒等リンク」という特殊ケースに過ぎません。目的変数の性質（二値・件数・比率・非負の金額など）に応じて分布とリンク関数を選び直すことで、線形回帰では対応できなかった多様なデータに、同じ考え方の延長線上でモデルを当てはめられるようになります。

主なGLMの対応表：目的変数のタイプ別に分布とリンク関数を選ぶ

実務でGLMを使う際にまず押さえておきたいのは、目的変数がどのようなタイプのデータかによって、選ぶべき確率分布とリンク関数がある程度定型的に決まるという点です。次の対応表は、代表的な組み合わせと、それぞれがビジネスの現場でどのような場面に対応するかをまとめたものです。

この表からも分かる通り、目的変数が二値かどうか、整数の件数かどうか、負の値を取り得るかどうか、といったごく基本的なデータの性質を確認するだけで、選ぶべき確率分布はかなりの程度絞り込めます。次章以降で詳しく扱うロジスティック回帰（二項分布+ロジットリンク）とポアソン回帰（ポアソン分布+対数リンク）は、この表の中でも実務での登場頻度が特に高い代表例にあたります。

最尤法：最小二乗法に代わるパラメータ推定の考え方

通常の線形回帰では、残差の二乗和を最小にする係数を選ぶ最小二乗法によってパラメータを推定してきました。しかしGLMでは正規分布以外の確率分布も扱うため、二乗和という発想がそのままでは意味をなさない場面が出てきます。そこでGLM全体を通じて使われるパラメータ推定の考え方が、最尤法（さいゆうほう、Maximum Likelihood Estimation）です。

最尤法の直感は、「手元に観測されたデータが得られる確率（尤度）が最も高くなるように、パラメータを選ぶ」というものです。ある係数の組み合わせを仮に置いたとき、その仮定のもとで実際に観測されたデータが生じる確率、あるいは確率密度を計算できます。この確率を尤度と呼び、係数の値を色々と動かしながらこの尤度が最大になる点を探す、というのが最尤法の基本発想です。「このデータが実際に観測された以上、そのデータが最も出やすくなるようなパラメータの値が、最も尤（もっと）もらしい」という考え方だと捉えると分かりやすいと思います。

ここで重要なのは、目的変数が正規分布に従うと仮定した場合、最尤法によって導かれる係数の推定値は、最小二乗法によって導かれる係数の推定値と完全に一致するという事実です。つまり、これまでの章で扱ってきた最小二乗法は、実は「正規分布を仮定した最尤法」という、より一般的な推定原理の1つの表れに過ぎなかったことになります。GLMでは、正規分布以外の分布（二項分布やポアソン分布など）に対しても、それぞれの分布の尤度を最大にするという同じ原理でパラメータを推定します。ほとんどの統計ソフトウェアは、この尤度を数値的に最大化するアルゴリズム（反復重み付け最小二乗法など）を内部で自動的に実行してくれるため、利用者が最適化の計算そのものを手で行う必要はありません。

GLMの評価：逸脱度（deviance）とAICの適用

通常の線形回帰では、モデルの当てはまりの良さを決定係数\( R^2 \)で、モデル選択をAICで評価してきました。GLMでは、正規分布以外の分布も扱う都合上、決定係数に代わる指標として逸脱度（deviance）という概念が使われます。

逸脱度は、当てはめたモデルの尤度と、データに完全に適合する理想的な（飽和）モデルの尤度とを比較し、その差を2倍したもので定義されます。ざっくり言えば、モデルが理想的な当てはまりからどれだけ乖離しているかを、尤度の言葉で測った指標であり、値が小さいほど当てはまりが良いことを意味します。正規分布を仮定したGLM（つまり通常の線形回帰）における逸脱度は、実は残差二乗和とほぼ同じ役割を果たす量になります。この対応関係からも、逸脱度が最小二乗法における残差二乗和の考え方を、より一般の分布に拡張した指標であることが分かります。

AIC（赤池情報量規準）についても、GLMの枠組みでそのまま拡張して適用できます。AICは尤度を基礎にした指標であるため、「最尤法によって尤度を最大化する」というGLM共通の推定原理と自然に整合しており、確率分布やリンク関数の異なる複数のモデルを比較する際にも、前章で扱ったのと同じ考え方（AICが小さいモデルほど、説明力と複雑さのバランスが良い）で利用できます。ただし、異なる確率分布を仮定したモデル同士のAICを比較する場合には、尤度の計算の前提が異なる点に注意が必要です。同じ確率分布族の中で、説明変数の組み合わせを変えて比較する使い方が最も安全で、分布そのものが異なるモデル間の比較は慎重に扱う必要があります。

逸脱度を使った当てはまりの確認では、逸脱度をその自由度で割った値（尺度化逸脱度）がおよそ1に近いかどうかも、実務上の目安としてよく確認されます。この値が1から大きく離れて大きい場合、想定した確率分布ではデータのばらつきをうまく説明できていない（過分散などが起きている）可能性を示唆しており、次章以降で扱うポアソン回帰の過分散問題とも関わってきます。

モデル選択の実務：どの分布を選ぶかを決める手順

実務でGLMを使う際にまず悩むのは、目的変数に対してどの確率分布を選べばよいかという判断です。次のような手順で確認していくと、判断の抜け漏れを減らせます。

1. 目的変数が連続値か離散値かを確認します。連続値であれば正規分布かガンマ分布、離散値であれば二項分布やポアソン分布、負の二項分布が候補になります。
2. 目的変数が取り得る値の範囲を確認します。負の値も取り得るなら正規分布、0以上の整数しか取らないなら二項分布かポアソン系、0か1の2値しか取らないなら二項分布、0から1の間の連続的な比率なら二項分布かベータ分布が候補になります。
3. 目的変数の平均とばらつき（分散）の関係を確認します。平均が大きくなるほど分散も大きくなる傾向が強い件数データでは、ポアソン分布（平均=分散という制約を持つ分布）がまず候補になりますが、分散が平均を大きく上回る過分散が見られる場合は、負の二項分布への切り替えを検討します。この点は第9章で詳しく扱います。
4. 候補となる分布とリンク関数の組み合わせでGLMを当てはめ、逸脱度やAIC、残差の診断プロットを確認しながら、当てはまりの悪いモデルを候補から外していきます。
5. 最終的に残った候補の中から、当てはまりの良さだけでなく、係数の解釈のしやすさや、意思決定の現場での説明のしやすさも踏まえて、採用するモデルを決定します。

この手順を踏まえると、GLMの実務は「とりあえず正規分布で回帰式を当ててみる」という発想から、「目的変数がどのように生成されたデータなのかをまず観察し、それに見合った分布を選んでからモデルを組む」という発想への転換であることが分かります。この転換こそが、統計モデリングという考え方の核心にあたります。

Pythonコード例：statsmodelsによるGLMの基本形

Pythonでは、statsmodelsライブラリがsm.GLMクラス、およびR言語ライクな数式記法を使えるsmf.glm関数によって、GLMを直接扱えるインターフェイスを提供しています。まずは、目的変数と説明変数の行列を直接渡すsm.GLMの基本形を確認します。ここでは、ある工程の稼働条件から不良品の発生件数を予測する、ポアソン分布を仮定したGLMを例にします。

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

# ある工程における、稼働時間・シフト・不良件数のダミーデータ
rng = np.random.default_rng(0)
n = 200
operating_hours = rng.uniform(4, 12, size=n)
night_shift = rng.integers(0, 2, size=n)

# 稼働時間が長く、夜勤であるほど不良件数が増える生成過程を仮定
lam = np.exp(-1.2 + 0.18 * operating_hours + 0.35 * night_shift)
defect_count = rng.poisson(lam)

X = pd.DataFrame({
    "operating_hours": operating_hours,
    "night_shift": night_shift,
})
X = sm.add_constant(X)  # 切片項を追加
y = defect_count

# 確率分布にポアソン分布、リンク関数に対数リンクを指定
model = sm.GLM(y, X, family=sm.families.Poisson())
result = model.fit()

print(result.summary())
print(f"逸脱度(deviance): {result.deviance:.2f}")
print(f"AIC: {result.aic:.2f}")

ここでfamily=sm.families.Poisson()という指定が、GLMの3部品のうち「確率分布」を選ぶ部分にあたります。statsmodelsの各family（sm.families.Gaussian()、sm.families.Binomial()、sm.families.Poisson()、sm.families.Gamma()、sm.families.NegativeBinomial()など）は、それぞれ標準的なリンク関数をあらかじめ内蔵しているため、多くの場合はfamilyを指定するだけでリンク関数まで自動的に決まります。上記の例ではポアソン分布の標準リンクである対数リンクが自動的に使われています。

係数の推定結果であるresult.paramsやresult.summary()の出力の見方は、これまでの章で扱ってきた重回帰の出力と基本的に共通しており、各係数、標準誤差、p値、信頼区間が並びます。ただし対数リンクを経由しているため、係数の解釈にはひと工夫が必要になります。この解釈の詳細は、次章以降でロジスティック回帰・ポアソン回帰それぞれについて丁寧に扱います。

数式記法を使いたい場合は、Rのモデル式に近い書き方ができるsmf.glmを使うと、データフレームの列名をそのまま数式の中で扱えて便利です。同じデータに対して、リンク関数を明示的に指定する書き方も含めて確認します。

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

df = X.drop(columns="const").copy()
df["defect_count"] = y

# 数式記法でGLMを指定。リンク関数を明示的に対数リンクとして渡す例
formula = "defect_count ~ operating_hours + night_shift"
model2 = smf.glm(
    formula=formula,
    data=df,
    family=sm.families.Poisson(link=sm.families.links.Log()),
)
result2 = model2.fit()

print(result2.summary())

family=sm.families.Poisson(link=sm.families.links.Log())という書き方をすると、確率分布とリンク関数の組み合わせを明示的に指定できます。二項分布であればsm.families.Binomial(link=sm.families.links.Logit())、ガンマ分布であればsm.families.Gamma(link=sm.families.links.Log())というように、familyの引数にlinkを渡すことで、対応表に挙げた組み合わせを自由に切り替えられます。この共通のインターフェイスこそが、GLMという1つの枠組みの中で多様な目的変数を扱えることの実装面での現れだと言えます。

本章のまとめと次章予告

本章では、目的変数が0/1、件数、割合、非負の金額といった、通常の線形回帰の前提が崩れる場面において、確率分布・線形予測子・リンク関数という3つの部品の組み合わせでモデルを柔軟に組み立て直せる、GLMという枠組みを見てきました。正規線形モデルがGLMの特殊ケースに過ぎないという見方に立てば、これまでの章で積み上げてきた回帰分析の知識の多くは、そのままGLM全体の基礎知識として活かせます。パラメータの推定原理が最小二乗法から最尤法へと一般化され、モデルの当てはまりの評価も決定係数から逸脱度・AICへと一般化される、という対応関係も押さえておくと、次章以降の理解がスムーズになります。

次章では、GLMの中でも実務での登場頻度が最も高いロジスティック回帰を、二項分布とロジットリンクという組み合わせに焦点を当てて詳しく扱います。オッズ比による係数の解釈や、分類精度の評価指標といった、ロジスティック回帰ならではの実務的な論点を掘り下げます。続く第9章では、件数データを扱うポアソン回帰と、その前提が崩れたときに直面する過分散という問題、そしてその対処法である負の二項回帰について扱います。

この章を深めたい方への参考書籍

『データ解析のための統計モデリング入門』（久保拓弥、岩波書店）：GLMを、確率分布・線形予測子・リンク関数という3要素から一般化線形モデルとして丁寧に解説しており、本章で扱った内容をより数式的な裏付けとともに深く理解したい読者に適した定番の1冊です。植物の種子数データなどを題材にポアソン回帰を導入する構成は、本章の対応表の考え方をより実感を伴って掴み直す助けになると思います。

第8章 要因分析の道具としてのロジスティック回帰

本シリーズの機械学習編第6章では、ロジスティック回帰を「解約するかどうか」「与信が通るかどうか」といった二値の結果をあらかじめ予測するための分類器として扱いました。本章で扱うのは同じロジスティック回帰ですが、目的が異なります。予測目的の使い方は機械学習編で扱ったため、本章は要因分析の道具として、すなわち「どの要因が結果にどれだけ影響しているか」を語るための使い方に絞って解説します。第7章で整理した一般化線形モデル（GLM）の枠組み、確率分布とリンク関数と線形予測子の3点セットに沿って考えると、ロジスティック回帰はその最初の実践例にあたります。二項分布とロジットというリンク関数の組み合わせが、なぜ二値の結果を扱うのに適しているのかというところから話を始めます。

二値の結果を確率で説明するGLM

これまでの章で扱ってきた重回帰分析は、売上や所得のように連続した値を取る目的変数を前提にしていました。しかし実務でよく直面するのは、「解約したかしなかったか」「成約したか失注したか」「不良品が出たか出なかったか」のように、結果が2つの値しか取らない場面です。このような目的変数を通常の重回帰でそのまま扱うと、第4章で確認した残差の等分散性や正規性の前提が崩れ、係数の標準誤差や信頼区間が正しく計算できなくなります。二値の結果が従う分布はベルヌーイ分布（何回か繰り返せば二項分布）であり、連続値を前提にした重回帰の枠組みとは相性がよくありません。

ここで役立つのが第7章のGLMの考え方です。GLMは、目的変数が従う確率分布と、線形予測子（説明変数の重み付き和）を確率分布のパラメータに変換するリンク関数を指定することで、重回帰の技術を様々な種類のデータへ拡張します。ロジスティック回帰は、確率分布に二項分布を、リンク関数にロジット（対数オッズ）を採用したGLMの特殊ケースです。数式で書くと次のようになります。

\( \log\left(\dfrac{p}{1-p}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p \)

左辺の \( p \) は「結果が1（解約、成約、不良など、関心のある事象）である確率」を表します。右辺は重回帰とまったく同じ形の線形結合です。つまりロジスティック回帰は、確率\(p\)そのものではなく、確率を対数オッズという形に変換した量を、重回帰の技術でそのまま予測している仕組みだと理解できます。この「右辺は重回帰と同じ、左辺の変換だけが違う」という見方を持っておくと、次に登場するポアソン回帰（第9章）でも同じ発想が繰り返されることが見えてきます。

オッズとオッズ比の理解

ロジスティック回帰の係数を正しく読むには、確率とオッズの違いを押さえておく必要があります。確率\(p\)は「その事象が起きる割合」で、0から1の範囲に収まります。これに対してオッズは、事象が起きる確率と起きない確率の比、\(p/(1-p)\)のことです。解約確率が0.2であれば、オッズは\(0.2/0.8=0.25\)となり、「解約は継続の4分の1の起きやすさ」と読みます。確率が0.8であればオッズは4になり、「解約は継続の4倍起きやすい」という意味になります。

ロジスティック回帰が確率そのものではなくオッズの対数を線形予測子として使うのには理由があります。確率\(p\)は0から1という範囲に制約されているため、説明変数の重み付き和（理論上マイナス無限大からプラス無限大まで取り得る）と直接結びつけると、値の取り得る範囲が食い違います。オッズ\(p/(1-p)\)は0から無限大までの値を取れるようになりますが、まだ下限がゼロという制約が残ります。さらに対数を取ることで、対数オッズは実数全体（マイナス無限大からプラス無限大）を取り得るようになり、線形予測子と値の範囲が一致します。この「対数オッズが線形結合と一致する」という性質を、ロジットの線形性と呼びます。この変換によって、確率という扱いにくい量を、重回帰と同じ加法的な形で表現できるようになる点が、ロジスティック回帰の要になっています。

オッズ比は、あるグループのオッズと別のグループのオッズの比です。たとえば問い合わせ履歴のある顧客のオッズが0.5、ない顧客のオッズが0.2であれば、オッズ比は\(0.5/0.2=2.5\)となり、「問い合わせ履歴があると解約のオッズがおよそ2.5倍になる」と表現します。オッズ比が1より大きければその要因は事象を起きやすくする方向に、1より小さければ起きにくくする方向に働き、1であれば影響がないと解釈します。

ここで実務上もっとも誤解が生じやすいのが、オッズ比を「確率の比（リスク比）」と同じものだと思い込んでしまうことです。オッズ比とリスク比は、事象の発生率（ベースラインの確率）が低いときにはよく似た値になりますが、発生率が高くなるにつれて両者は大きく乖離していきます。具体的な数値で確認します。

発生率が5%程度の稀な事象であれば、オッズ比2.0は「確率がおよそ2倍になる」という近似的な読み方をしてもさほど誤差は出ません。しかし発生率が40%まで上がると、オッズ比2.0は「確率が1.43倍になる」ことにしか対応せず、「2倍になる」という説明は事実よりかなり誇張された印象を与えてしまいます。解約分析のように解約率が数十パーセントに達することもある領域では、この乖離は無視できない大きさになります。

それでもオッズ比が統計の実務で標準的な指標であり続けているのには理由があります。オッズ比は、症例対照研究（先に結果が起きたグループと起きなかったグループを集めてから要因を遡って調べる研究デザイン）のように、発生率そのものを直接推定できないデータからでも計算できるという性質を持っています。加えて、ロジスティック回帰の係数がそのまま対数オッズ比に対応するという数学的な扱いやすさも、オッズ比が定着した理由の一つです。実務上大切なのは、オッズ比という指標そのものを避けることではなく、発生率が高い分析ではオッズ比とリスク比が乖離することを理解した上で、必要に応じて後述する限界効果や、モデルが出力する実際の確率も併記するという姿勢です。

オッズ比は「確率が何倍になるか」ではなく「オッズが何倍になるか」を表す数字です。ベースラインの発生率が低い場合は確率の比に近づきますが、発生率が高い場合はオッズ比の方が実際の確率の変化より大きく見えてしまいます。解約率や成約率が高い分析でオッズ比を報告するときは、あわせて実際の確率の変化（ベースラインの確率とモデルが示す確率を並べて示す）も添えることが、誤解を避ける実務上の作法です。

係数の解釈：exp（係数）の読み方

ロジスティック回帰の推定結果に並ぶ係数\(\beta_j\)は、対数オッズの世界での「1単位あたりの変化量」です。これをオッズの世界に戻すには、\(\exp(\beta_j)\)を計算します。対数の世界での足し算は元の世界では掛け算に対応するため、\(\exp(\beta_j)\)がそのままオッズ比になります。この変換さえ押さえておけば、重回帰で係数を読むときと同じ感覚で、ロジスティック回帰の係数も読み進めることができます。

連続変数の場合、\(\exp(\beta_j)\)は「他の説明変数を一定に保ったまま、その変数が1単位増えたときにオッズが何倍になるか」を表します。たとえば契約期間（月数）の係数が\(-0.03\)であれば、\(\exp(-0.03)\approx0.97\)となり、「契約期間が1か月延びるごとに解約のオッズがおよそ3%下がる」と読みます。カテゴリ変数の場合は、重回帰のダミー変数と同じ考え方で、基準カテゴリとの比較になります。プラン（基準:プレミアム）のうち「ベーシック」ダミーの係数が0.4であれば、\(\exp(0.4)\approx1.49\)となり、「ベーシックプランはプレミアムプランと比べて解約のオッズがおよそ1.49倍」と解釈します。基準カテゴリの選び方によって数字の見え方が変わる点も、第2章で扱った重回帰のダミー変数の解釈とまったく同じです。

点推定値だけでなく、信頼区間つきで語ることも欠かせません。オッズ比の95%信頼区間は、係数の信頼区間の両端をそれぞれ\(\exp\)で変換して求めます。信頼区間が1をまたいでいれば、そのオッズ比は統計的に有意な差があるとは言い切れず、逆に区間が1を上回る範囲（あるいは下回る範囲）に収まっていれば、方向性を伴った関係があると判断できます。第2章で扱ったp値と信頼区間の関係と同じ論理が、ここでもそのまま成り立ちます。statsmodelsのconf_int()が返す信頼区間は、既定では係数の標準誤差を正規分布近似で使うワルド型の信頼区間です。サンプル数が少ない場合やサンプル数に対してイベント数が乏しい場合は、この近似の精度が落ちることがあるため、後述するサンプル数の目安とあわせて確認しておくと安心です。オッズ比の点推定だけを取り出して報告するのではなく、区間の幅を添えて語ることが、係数解釈における基本作法だと考えています。

実務例：解約要因の分析

ここまでの内容を、月額制のオンライン学習サービスを想定した解約要因分析の合成データで確認します。契約期間（月数）、週あたりの利用頻度、過去3か月の問い合わせ件数、契約プランの4つの説明変数から、解約の有無を説明するモデルを組み立てます。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(0)
n = 3000

tenure_months = rng.integers(1, 37, n)      # 契約期間(月)
usage_freq = rng.poisson(3, n)              # 週あたりの利用回数
num_inquiries = rng.poisson(0.5, n)         # 過去3か月の問い合わせ件数
plan = rng.choice(["ベーシック", "プレミアム"], size=n, p=[0.7, 0.3])

# 解約に至る対数オッズを仕込む(契約期間・利用頻度が上がると解約しにくく、
# 問い合わせが多い、ベーシックプランだと解約しやすい)
z = (
    -0.03 * tenure_months
    - 0.30 * usage_freq
    + 0.45 * num_inquiries
    + 0.40 * (plan == "ベーシック")
    - 0.20
)
prob_churn = 1 / (1 + np.exp(-z))
churn = rng.binomial(1, prob_churn)

df = pd.DataFrame({
    "tenure_months": tenure_months,
    "usage_freq": usage_freq,
    "num_inquiries": num_inquiries,
    "plan": plan,
    "churn": churn,
})

model = smf.logit(
    "churn ~ tenure_months + usage_freq + num_inquiries + "
    "C(plan, Treatment(reference='プレミアム'))",
    data=df,
).fit()
print(model.summary())

smf.logitはformula APIでロジスティック回帰を推定する関数で、C()によるカテゴリ変数の扱いや基準カテゴリの指定は、第2章で確認したsmf.olsとまったく同じ書き方ができます。model.summary()の出力には、各係数と標準誤差、z値、p値が並びます。ここから係数をオッズ比に変換し、信頼区間とあわせて一覧にします。

conf = model.conf_int()
conf.columns = ["ci_low", "ci_high"]

odds_table = np.exp(pd.concat([model.params, conf], axis=1))
odds_table.columns = ["odds_ratio", "ci_low", "ci_high"]
odds_table["p_value"] = model.pvalues
print(odds_table.round(3))

この結果、tenure_monthsとusage_freqのオッズ比が1を下回り、その95%信頼区間も1をまたいでいなければ、「契約期間が長い顧客ほど、また利用頻度が高い顧客ほど、解約のオッズが統計的に有意に下がる」と読めます。num_inquiriesのオッズ比が1を上回っていれば、「問い合わせが多い顧客ほど解約しやすい」ことになり、これは製品への不満やトラブルが解約の予兆になっている可能性を示唆します。C(plan)[T.ベーシック]のオッズ比が1を上回っていれば、同じ契約期間・利用頻度・問い合わせ件数の顧客同士を比べても、ベーシックプランの方が解約しやすい、という関係が見えてきます。

ここから導ける施策の示唆はいくつか考えられます。利用頻度の低下は解約の先行指標になり得るため、一定期間ログインが減った顧客に対してオンボーディングやリマインドの施策を打つこと。問い合わせ件数が多い顧客は不満を抱えている可能性が高いため、問い合わせ対応の質を高めるか、問い合わせに至る前の製品側の改善を優先すること。ベーシックプランの解約オッズが高いのであれば、プレミアムプランへのアップグレード提案や、ベーシックプラン利用者向けのサポート強化を検討すること、などです。これらはあくまでモデルが示す関連であり、施策の効果を確定させるには、本シリーズの因果推論と効果検証編で扱う効果検証（A/Bテストなど）の手続きが別途必要になる点は付け加えておきます。

限界効果：オッズ比が伝わらない相手への説明

オッズ比は「何倍になるか」という掛け算の言葉であり、統計に馴染みのない相手にはかえって伝わりにくいことがあります。「利用頻度が週1回増えると解約のオッズが0.76倍になる」と言われても、それが実際にどれくらい解約を減らすことになるのか、直感的にはつかみにくいという声はよく聞かれます。こうした相手には、オッズ比の代わりに限界効果（marginal effect）を使うと伝わりやすくなります。限界効果とは、ある説明変数が1単位増えたときに、確率そのものが何ポイント変化するかを表す指標です。

ロジスティック回帰では、シグモイド型の曲線を使って確率を計算しているため、説明変数が1単位増えたときの確率の変化幅は、他の説明変数の水準によって変わります。契約直後の顧客と契約3年目の顧客とでは、利用頻度が1回増えたときの解約確率の下がり方が同じとは限りません。この性質があるため、限界効果は代表的な条件（全観測点の平均、あるいは各観測点ごとに計算した値の平均）を定めた上で算出するのが一般的です。観測点ごとに計算した限界効果を平均する方法は平均限界効果（average marginal effects、AME）と呼ばれ、statsmodelsのget_margeff()メソッドで求められます。

margeff = model.get_margeff(at="overall", method="dydx")
print(margeff.summary())

この出力は、各説明変数について「その変数が1単位増えたときに、解約確率が平均して何ポイント変化するか」を示します。たとえばusage_freqの平均限界効果が\(-0.04\)であれば、「利用頻度が週1回増えると、解約確率はおよそ4ポイント下がる」とそのまま説明できます。オッズ比よりも直感的に理解しやすく、経営層への説明や施策の優先順位づけの材料として使いやすい表現です。

get_margeff()のat引数には、限界効果をどの条件で計算するかを指定します。at=”overall”を指定すると、観測点ごとに限界効果を計算してから平均する平均限界効果（AME）が得られます。一方、at=”mean”を指定すると、すべての説明変数をあらかじめ平均値に固定した上で、その代表的な1点における限界効果を計算します（これは平均的な限界効果、MEMと呼ばれ、AMEとは計算の順序が異なります）。カテゴリ変数が含まれる場合、平均値に固定するという操作の意味が分かりにくくなるため、実務ではAMEの方が解釈しやすく、報告にも使いやすい場面が多くなります。いずれの方法でも、限界効果はあくまである条件のもとでの近似的な変化量であり、契約期間や利用頻度が極端な値の顧客には同じ変化幅がそのまま当てはまらない点は留意しておく必要があります。

分離（separation）の問題

ロジスティック回帰を実データで動かしていると、まれに係数の推定値が異常に大きな値になったり、標準誤差が発散したように見えたりすることがあります。多くの場合、これは分離（separation）と呼ばれる現象が原因です。分離とは、ある説明変数（や説明変数の組み合わせ）の値によって、目的変数の1と0を完全に、あるいはほぼ完全に分け切れてしまう状態を指します。

たとえば、ある期間より前に契約した顧客が1人も解約しておらず、それ以降に契約した顧客だけが解約している、というデータがあると、契約時期という変数だけで解約の有無を完全に説明できてしまいます。この場合、最尤推定は「その変数の係数をどこまでも大きくすればするほど、尤度（データの当てはまりの良さ）が上がり続ける」という状況に陥り、反復計算が収束しないか、非常に大きな（実務的には無意味な）係数値に落ち着いてしまいます。分離が起きやすいのは、サンプル数が少ない場合、カテゴリ変数のある水準に該当するサンプルが極端に少ない場合、説明変数の数に対してイベント数（解約や成約など関心のある事象の件数）が少なすぎる場合などです。

分離が疑われるときの対処法はいくつかあります。

• まずクロス集計や分割表で、疑わしい説明変数と目的変数の組み合わせを確認し、完全に、あるいはほぼ完全に分かれているカテゴリや区間がないかを点検する
• 該当するカテゴリのサンプルを他のカテゴリと統合する、あるいは問題を起こしている変数を一旦モデルから外して係数の挙動を確認する
• 第10章で扱う正則化回帰（L2正則化を加えたロジスティック回帰、いわゆるRidge型のペナルティ）を使い、係数が無限に発散するのを防ぐ
• Firth法（バイアス補正した最尤推定、罰則付き尤度を用いることで分離が起きている状況でも有限の係数と標準誤差が得られる手法）の利用を検討する。Firth法はstatsmodelsの標準機能には含まれていませんが、専用の外部パッケージや、統計解析向けの他のソフトウェアで実装が提供されています

係数の絶対値が極端に大きい、標準誤差が異常に大きい、あるいは推定が収束しなかった旨の警告が出た場合は、まず分離を疑うことが大切です。データを確認せずにそのまま結果を報告すると、実態のない極端な係数を「非常に強い効果がある」と誤読してしまう危険があります。

サンプル数と稀な事象

ロジスティック回帰の係数の安定性は、単純なサンプル数ではなく、目的とする事象（解約、成約、不良など）がどれだけ発生しているか、というイベント数によって左右されます。全体のサンプル数が数万件あっても、解約という事象が数十件しか含まれていなければ、係数の推定は不安定になりやすく、標準誤差も大きくなります。

目安としてよく参照されるのが、1変数あたりのイベント数（events per variable、EPV）という考え方です。説明変数の数に対してイベント数が少なすぎると、係数の推定値にバイアスが生じたり、標準誤差が過小・過大に評価されたりしやすいことが、疫学分野を中心とした研究で報告されています。厳密な閾値があるわけではありませんが、EPVがおおむね10を下回ると係数の推定が不安定になりやすく、10から20程度を確保することが望ましいという目安が広く紹介されています。説明変数が5個のモデルであれば、少なくとも50件、できれば100件程度のイベント数が欲しいという計算になります。

解約率が1%未満のような非常に稀な事象を扱う場合、全体のサンプル数を増やすだけでは対応しきれないこともあります。こうした場面では、稀な事象に該当するサンプルを意図的に多めに抽出し（ケースコントロール型のサンプリング）、後から抽出比率に応じて補正を行うといった工夫が使われます。いずれの場合も、モデルを組む前に「関心のある事象が何件あるか」をまず確認する習慣が、後になって係数の不安定さに悩まされることを防ぎます。

適合度の確認

ロジスティック回帰の当てはまりの良さを確認する際、重回帰の決定係数のような単純な指標はそのままでは使えません。代わりによく参照されるのが逸脱度（deviance）です。逸脱度は、モデルの対数尤度を基にした指標で、値が小さいほど当てはまりが良いことを示します。切片だけのモデル（説明変数を一切使わない、最も単純なモデル）の逸脱度をnull deviance、実際に説明変数を投入したモデルの逸脱度をresidual devianceと呼び、両者の差が大きいほど、投入した説明変数がモデルの当てはまりを改善していると解釈できます。

もう一つよく登場するのが疑似R2（pseudo R-squared、代表的なものにMcFaddenの疑似R2があります）です。これはnull devianceとresidual devianceの比から計算され、statsmodelsのロジットモデルではmodel.prsquaredとして直接取得できます。ただし疑似R2は重回帰の決定係数とは尺度がまったく異なり、同じ0.3という値でも意味合いが違います。McFaddenの疑似R2は0.2から0.4程度でもかなり良い当てはまりとされることが多く、重回帰の決定係数の感覚でそのまま高低を判断すると誤解を招きます。疑似R2はモデル間の相対的な比較や大まかな目安として参考程度に使うにとどめ、それだけでモデルの良し悪しを断定しないことが実務上の姿勢として妥当です。

逸脱度や疑似R2に加えて、キャリブレーション（較正）的な見方も欠かせません。モデルが「解約確率60%」と判定した顧客を多数集めたとき、実際の解約割合がおよそ60%前後になっているかどうかを確認する視点です。逸脱度や疑似R2はモデル全体の当てはまりを1つの数字に要約しますが、キャリブレーションは確率の水準ごとに実際の発生割合と突き合わせることで、モデルの出力をそのまま業務判断に使ってよいかを確認する手立てになります。この確認は、次章のポアソン回帰や、不均衡データを扱う場面でも同じ考え方が繰り返し登場します。

この章を深めたい方への参考書籍

『ロジスティック回帰分析』（丹後俊郎・山岡和枝・高木晴良、朝倉書店）：ロジスティック回帰の数理的な背景から、係数の解釈、モデルの適合度評価、サンプルサイズの設計まで体系的にまとめた専門書です。本章で扱ったオッズ比の信頼区間や分離の問題、イベント数と係数の安定性といった論点を、より厳密な統計理論の立場から掘り下げたい読者に向いています。

第9章 件数データを扱うポアソン回帰と過分散への対処

第7章ではGLM（一般化線形モデル）という枠組みそのものを、第8章では二項分布とロジットリンクを組み合わせたロジスティック回帰を取り上げました。本章ではGLMの枠組みを、来店客数や故障件数といった「件数データ」に適用します。件数データは経営やオペレーションの現場で日常的に扱う指標であるにもかかわらず、通常の線形回帰をそのまま当てはめると不自然な結果を招きやすいという特徴があります。

件数データにGLMを適用する代表的な手法がポアソン回帰です。本章では、ポアソン回帰の考え方と係数の解釈、規模の違いを補正するオフセットの使い方、そして実務でしばしば直面する過分散という問題への対処法を、合成データを使った具体例とともに整理します。

カウントデータの性質

件数データとは、店舗の来客数、設備の故障件数、コールセンターへの着信数、製品の不良品数のように「ある期間・ある単位の中で何回起きたか」を数えた変数のことです。件数データには、通常の連続量とは異なる3つの性質があります。

• 必ず0以上の整数しか取らない（マイナスの来客数や1.5件の故障は存在しない）
• 実務データでは0の観測が多く含まれやすい（閑散日、無故障日、着信ゼロの時間帯など）
• 平均が大きくなるほど、ばらつき（分散）も大きくなる傾向がある（来客数が多い店舗ほど日々の変動幅も大きい）

この性質を無視して通常の線形回帰（最小二乗法）を件数データに当てはめると、いくつかの不具合が生じます。まず、線形回帰の予測値は理論上マイナス無限大からプラス無限大までの値を取り得るため、説明変数の組み合わせによっては「来客数がマイナス3人」といった、現実にはあり得ない予測値を出してしまうことがあります。また、線形回帰は誤差の分散が説明変数の値によらず一定であることを前提にしていますが、件数データでは平均が大きい観測ほど分散も大きくなるため、この前提が崩れています。この結果、係数の標準誤差やp値が歪み、統計的な検定の信頼性が損なわれます。

こうした不具合を避けるために、件数データにはあらかじめ「0以上の整数で、平均とともに分散が変化する」という性質を織り込んだ確率分布を使う必要があります。その代表がポアソン分布であり、これを目的変数の確率分布として採用したGLMがポアソン回帰です。

ポアソン回帰の仕組み：対数リンクによる件数のモデル化

ポアソン分布は、単位時間や単位区域あたりに稀な事象が何回起きるかを表す確率分布で、平均発生回数を表すパラメータ \( \lambda \) （ラムダ）を1つだけ持ちます。ポアソン分布の大きな特徴は、この \( \lambda \) が分布の期待値であると同時に分散でもある、つまり「平均と分散が一致する」という点です。この性質は本章の後半で扱う過分散の議論の出発点になります。

ポアソン回帰は、GLMの枠組みに沿って、目的変数の確率分布にポアソン分布を、リンク関数に対数（log）を採用したモデルです。説明変数 \( x_1, x_2, \ldots, x_p \) から求めた線形予測子を、対数リンクを介して \( \lambda \) に変換します。

\( \log(\lambda_i) = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \cdots + \beta_p x_{pi} \)

対数リンクを使う理由は、線形回帰のように線形予測子をそのまま \( \lambda \) として扱うと、係数の値によっては \( \lambda \) がマイナスになり得るためです。件数の平均発生回数がマイナスになることはあり得ませんから、これは都合が悪い性質です。対数リンクを挟むことで、線形予測子がどんな実数値を取っても、\( \lambda = \exp(\beta_0+\beta_1x_1+\cdots) \) は必ずプラスの値になります。指数関数の値域がプラスの実数全体である、という性質をそのまま利用した工夫です。

この対数リンクは、係数の解釈のしかたにも直接影響します。線形回帰であれば「係数の分だけ目的変数が足し算で変化する」と読みますが、ポアソン回帰では対数の世界で足し算になっている分、元の件数の世界では掛け算になります。具体的には、説明変数 \( x_j \) が1単位増えたとき、期待件数は \( \exp(\beta_j) \) 倍になります。この \( \exp(\beta_j) \) を発生率比（rate ratio）と呼び、ロジスティック回帰のオッズ比と同様に、係数を実務担当者に説明できる言葉に変換する役割を持ちます。\( \exp(\beta_j) \) が1より大きければその要因は件数を増やす方向に、1より小さければ減らす方向に働き、1であれば影響がないと解釈します。

実務例：店舗来客数の要因分析

ここからは、店舗の1日あたり来客数を題材にした合成データで、ポアソン回帰の一連の流れを確認します。週末か平日か、広告出稿額、近隣競合店舗数、気温という4つの説明変数が来客数にどう影響しているかを分析するという想定です。

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

rng = np.random.default_rng(0)
n = 300

# 週末フラグ、広告出稿額(万円)、近隣競合店舗数、気温(度)を用意する
weekend = rng.integers(0, 2, n)
ad_spend = rng.normal(10, 3, n).clip(0)
competitor = rng.integers(0, 4, n)
temperature = rng.normal(20, 8, n)

# 対数リンクの逆変換(exp)でlambdaを作り、ポアソン分布から来客数を生成する
log_lambda = 2.5 + 0.35 * weekend + 0.04 * ad_spend - 0.15 * competitor + 0.01 * temperature
lam = np.exp(log_lambda)
visitors = rng.poisson(lam)

df_store = pd.DataFrame({
    "visitors": visitors,
    "weekend": weekend,
    "ad_spend": ad_spend,
    "competitor": competitor,
    "temperature": temperature,
})

# family=sm.families.Poisson()を指定するとポアソン回帰になる
poisson_store = smf.glm(
    "visitors ~ weekend + ad_spend + competitor + temperature",
    data=df_store,
    family=sm.families.Poisson(),
).fit()
print(poisson_store.summary())

coef_table = pd.DataFrame({
    "coef": poisson_store.params,
    "exp(coef)": np.exp(poisson_store.params),
})
print(coef_table)

推定結果の見方はロジスティック回帰と共通しています。summary()に出力される係数（coef）そのものは対数の世界の値なので、そのままでは「週末になると来客数がどれだけ増えるか」を直感的に語れません。そこで各係数に \( \exp \) を取った発生率比を確認します。

この表の数値は生成に使ったパラメータに対応する目安であり、実際に乱数で生成したデータをフィッティングすると誤差の分だけ値がずれます。重要なのは、係数を対数のまま報告するのではなく、必ず \( \exp \) を取った発生率比に変換したうえで「何倍になるか」「何パーセント変化するか」という言葉に翻訳して業務側に伝えるという手順です。ロジスティック回帰のオッズ比と同じく、この変換を挟むことで初めて係数が経営判断に使える情報になります。

オフセット：規模あたりの発生率の比較

前節の例では、来客数を左右する要因として広告出稿額や競合店舗数を取り上げましたが、実務ではもう1つ見落としやすい要因があります。それは「観測対象の規模や機会の大きさ」です。営業時間が長い店舗ほど、稼働時間が長い設備ほど、単純に観測される件数は多くなります。1日の営業時間が15時間の店舗と6時間の店舗を、来客数の絶対数だけで比較するのは公平ではありません。

この「規模」を明示的にモデルへ組み込む仕組みがオフセットです。営業時間や稼働時間などの規模を表す変数をエクスポージャー（exposure）と呼び、その対数を線形予測子にそのまま（係数を1に固定した状態で）加えます。

\( \log(\lambda_i) = \log(\text{exposure}_i) + \beta_0 + \beta_1 x_{1i} + \cdots \)

この式は両辺から \( \log(\text{exposure}_i) \) を移項すると、\( \log(\lambda_i / \text{exposure}_i) = \beta_0+\beta_1x_1+\cdots \) という形になります。つまりオフセットを加えたポアソン回帰は、実質的に「件数そのもの」ではなく「規模あたりの発生率（単位時間あたりの故障件数など）」を目的変数として説明変数で説明していることになります。オフセットの係数をあらかじめ1に固定するのは、規模が2倍になれば件数もちょうど2倍になるはず、という比例関係を仮定しているためです。この仮定が成り立たない場合は、\( \log(\text{exposure}_i) \) を通常の説明変数としてモデルに加え、係数を推定させる方法もあります。

設備の稼働時間あたりの故障件数を例に、オフセットを使ったポアソン回帰を確認します。設備の経過年数と定期メンテナンスの有無が故障率に与える影響を、稼働時間の違いを補正したうえで推定します。

rng = np.random.default_rng(1)
n = 200

op_hours = rng.uniform(50, 500, n)       # 設備ごとの稼働時間
machine_age = rng.uniform(0, 15, n)      # 設備の経過年数
maintenance = rng.integers(0, 2, n)      # 定期メンテナンス実施の有無

# 単位時間あたりの故障率(rate)を経過年数とメンテナンスで決め、
# 稼働時間を掛けた期待故障件数からポアソン乱数で故障件数を生成する。
# あわせて、同じ経過年数でも壊れやすい個体とそうでない個体があるという
# 設備ごとの個体差を、平均1のガンマ分布で掛け合わせる(過分散の原因になる)
log_rate = -4.0 + 0.12 * machine_age - 0.5 * maintenance
rate = np.exp(log_rate)
expected_failures = rate * op_hours
frailty = rng.gamma(shape=1 / 0.5, scale=0.5, size=n)  # 平均1、分散0.5の個体差
failures = rng.poisson(expected_failures * frailty)

df_machine = pd.DataFrame({
    "failures": failures,
    "machine_age": machine_age,
    "maintenance": maintenance,
    "op_hours": op_hours,
})

# offsetにはあらかじめlogを取った値を渡す(生の稼働時間をそのまま渡してはいけない)
poisson_machine = smf.glm(
    "failures ~ machine_age + maintenance",
    data=df_machine,
    family=sm.families.Poisson(),
    offset=np.log(df_machine["op_hours"]),
).fit()
print(poisson_machine.summary())
print("経過年数の発生率比:", np.exp(poisson_machine.params["machine_age"]))
print("メンテナンスの発生率比:", np.exp(poisson_machine.params["maintenance"]))

offset引数には、稼働時間そのものではなく、あらかじめnp.log()で対数を取った値を渡す点が実務でつまずきやすいポイントです。生の値をそのまま渡すと、線形予測子に稼働時間そのものが足し合わされてしまい、意図した「規模あたりの発生率」のモデルになりません。オフセットを正しく組み込むと、machine_ageやmaintenanceの係数は「稼働時間の違いを取り除いたうえでの、故障のしやすさそのものへの影響」を表すようになり、稼働時間が長い設備ほど故障件数が多く見えるという見かけ上の影響を排除できます。

過分散：ポアソン分布の前提が崩れる問題

ポアソン回帰には、先に触れたとおり「平均と分散が一致する」という強い前提があります。しかし実務データでは、この前提が成り立たず、分散が平均を上回る現象がしばしば観測されます。これを過分散（overdispersion）と呼びます。

過分散が生じる典型的な原因は、モデルに含めていない要因によるばらつきです。たとえば店舗の来客数であれば、近隣でのイベント開催、突発的な悪天候、SNSでの話題化といった、説明変数に含めていない要因が日によって来客数を大きく上下させます。設備の故障件数であれば、個体差（同じ経過年数でも壊れやすい個体とそうでない個体がある）や、観測期間内でのメンテナンス品質のばらつきなどが該当します。こうした「モデルが捉えきれていない不均一さ」が積み重なると、観測される分散はポアソン分布が仮定する水準よりも大きくなります。

過分散を見逃したままポアソン回帰の結果を使うと、標準誤差が実際よりも小さく計算され、本来は誤差の範囲内に収まるはずの係数まで「統計的に有意」と判定されてしまいます。この見せかけの有意性をもとに、効果のない施策を効果ありと誤認したり、影響の小さい要因を過大評価して意思決定してしまったりする実害につながります。過分散のチェックは、ポアソン回帰を実務で使ううえで省略できない工程です。

過分散の診断：逸脱度とピアソンカイ二乗の目安

過分散が起きているかどうかを簡便に確認する方法として、モデルの逸脱度（deviance）またはピアソンカイ二乗統計量を、残差の自由度（df_resid）で割った値を見る方法があります。ポアソン分布の前提どおりに平均と分散が一致していれば、この比率はおおむね1に近い値になります。目安として、比率が1.5から2を超えてくると軽度の過分散、2を大きく上回るようであれば無視できない過分散が生じている、と判断されることが一般的です。

もう1つの確認方法は、残差プロットです。予測された \( \lambda \) の値（fittedvalues）に対して逸脱度残差やピアソン残差を散布図に描き、予測値が大きくなるにつれて残差の散らばりが扇状に広がっていないかを確認します。ばらつきが予測値とともに拡大しているようであれば、過分散が疑われます。

対処法：準ポアソンと負の二項回帰

過分散が確認された場合の対処法は大きく2つあります。1つ目は、準ポアソン（quasi-Poisson）と呼ばれる考え方です。これは、係数の点推定値自体はそのままに、標準誤差だけを推定された分散拡大係数（dispersion）で調整するという方法です。statsmodelsでは、GLMのfit()メソッドにscale=”X2″を指定することで、ピアソンカイ二乗に基づく分散拡大係数を用いた標準誤差の調整を行えます。

2つ目は、負の二項回帰（negative binomial regression）です。こちらはポアソン分布そのものを、分散が平均よりも大きくなることをあらかじめ許容する負の二項分布に置き換えるアプローチです。負の二項分布は、平均 \( \mu \) に加えてばらつきパラメータ \( \alpha \) を持ち、分散は \( \mu + \alpha \mu^2 \) という形で表されます。\( \alpha \) が0に近づくほどポアソン分布に近い挙動になり、\( \alpha \) が大きいほど平均に対して分散が大きくなることを許容します。準ポアソンが標準誤差だけを事後的に補正するのに対し、負の二項回帰は分布の形そのものを変えて \( \alpha \) を他の係数と同時に推定する点が異なります。

先ほどの設備の故障件数データを使って、過分散の診断、準ポアソンによる標準誤差の調整、負の二項回帰による再推定を順に行います。

# 逸脱度とピアソンカイ二乗を残差の自由度で割り、過分散の程度を確認する
deviance_ratio = poisson_machine.deviance / poisson_machine.df_resid
pearson_ratio = poisson_machine.pearson_chi2 / poisson_machine.df_resid
print("逸脱度/自由度:", deviance_ratio)
print("ピアソンカイ二乗/自由度:", pearson_ratio)

# 準ポアソン的な調整(点推定は変えず、標準誤差だけをX2で補正する)
poisson_quasi = smf.glm(
    "failures ~ machine_age + maintenance",
    data=df_machine,
    family=sm.families.Poisson(),
    offset=np.log(df_machine["op_hours"]),
).fit(scale="X2")
print(poisson_quasi.summary())

# 負の二項回帰(ばらつきパラメータalphaを他の係数と同時に最尤推定する)
nb_machine = smf.negativebinomial(
    "failures ~ machine_age + maintenance",
    data=df_machine,
    offset=np.log(df_machine["op_hours"]),
).fit()
print(nb_machine.summary())
print("推定されたばらつきパラメータalpha:", nb_machine.params["alpha"])

逸脱度/自由度やピアソンカイ二乗/自由度が1を大きく上回っていれば過分散が疑われ、標準のポアソン回帰の標準誤差をそのまま信用すべきではない、というサインになります。準ポアソンの結果では係数の点推定値そのものは通常のポアソン回帰と一致しますが、標準誤差（std err）とp値が過分散の分だけ大きくなり、より保守的な（誤って有意と判定しにくい）検定結果が得られます。負の二項回帰の結果では、alphaが0から統計的に離れているほど過分散が実際に存在することの裏付けになり、machine_ageやmaintenanceの係数についても、ポアソン回帰より広めの信頼区間が得られるのが通常です。どちらの方法も、過分散を無視した場合に比べて偽陽性（見せかけの有意）を防ぐ方向に働きます。

ゼロ過剰：ゼロが多すぎるデータへの対応

過分散とは別に、件数データではゼロの観測が想定以上に多く含まれる場合があります。たとえば、ある製品の月間購入個数を顧客ごとに集計すると、そもそも購入しない顧客が大半を占め、ゼロが極端に多い分布になることがあります。このようなデータにポアソン回帰や負の二項回帰をそのまま当てはめても、ゼロの多さをうまく説明しきれず、当てはまりが悪くなることがあります。

この問題への対応として、ゼロ過剰ポアソン回帰（zero-inflated Poisson regression）とハードルモデル（hurdle model）という2つのアプローチが知られています。ゼロ過剰ポアソン回帰は、観測されるゼロを「そもそも発生する見込みがない構造的なゼロ」と「発生する見込みはあるがたまたま0回だったゼロ」の2種類に分けて考え、前者をロジスティック回帰的な二値モデルで、後者を通常のポアソン分布でモデル化し、両者を混合させる枠組みです。ハードルモデルは、まず「0回か、1回以上か」を二値分類モデルで説明し、1回以上だった場合の具体的な件数は、0を除いた切断ポアソン分布や切断負の二項分布で説明するという、2段階の枠組みを取ります。statsmodelsにはZeroInflatedPoissonやZeroInflatedNegativeBinomialPといったクラスが用意されており、通常のポアソン回帰や負の二項回帰と同様にformulaベースで指定して推定できます。どちらのモデルも、ゼロの構造まで踏み込んだ分析が必要な場面で検討する選択肢として押さえておくとよいでしょう。

件数データを分析する際の実務フローの目安は次のとおりです。

1. 目的変数が0以上の整数（件数）であることを確認する
2. まずポアソン回帰（必要に応じてオフセットを追加）で基本モデルを当てはめる
3. 逸脱度/自由度やピアソンカイ二乗/自由度で過分散の有無を確認する
4. 過分散が見られる場合は、準ポアソンによる標準誤差の調整、または負の二項回帰へ切り替える
5. ゼロの割合が想定以上に多く、負の二項回帰でも当てはまりが悪い場合は、ゼロ過剰ポアソンやハードルモデルの適用を検討する

次章では、説明変数の数が多く多重共線性が疑われる場面や、予測精度を重視する場面で使われる正則化回帰（Ridge・Lasso）を取り上げます。ポアソン回帰や負の二項回帰も、リンク関数を通した線形予測子の枠組みは共通しているため、正則化の考え方はGLM全般に応用が利きます。

この章を深めたい方への参考書籍

『一般化線形モデル入門（原著第2版）』（Annette J. Dobson、共立出版）：ポアソン分布や負の二項分布を含む指数型分布族とリンク関数の理論的な位置づけを、数式の導出とともに体系的に解説しています。本章で扱った過分散の診断や負の二項回帰への拡張についても、なぜその手法が理論的に妥当なのかという背景から理解したい読者に向いた1冊です。

第10章 正則化回帰による係数の安定化

第5章では、説明変数どうしが強く相関していると係数の推定が不安定になる「多重共線性」という問題を取り上げました。VIFで共線性の度合いを診断し、変数を間引いたり主成分に置き換えたりする対処法を紹介しましたが、変数を削る以外にもう1つの有力な選択肢があります。それが本章のテーマである「正則化」です。

正則化という言葉自体は、機械学習編の第4章ですでに登場しています。そこでは主に「過学習を防ぐ手段」としてRidge回帰とLasso回帰を紹介しました。本章では同じ手法を回帰分析の文脈から見直します。過学習対策という側面に加えて、共線性下での係数の安定化、説明変数の数がサンプル数に迫る、あるいは超えるような高次元データでの推定、そして自動的な変数選択という、統計モデリングの実務に直結する3つの役割に焦点を当てます。あわせて、正則化された係数をどう解釈すべきかという、実務でしばしば誤解される論点も整理します。

正則化が回帰分析にもたらす3つの効果

通常の最小二乗法（OLS）は、残差平方和を最小にする係数を何の制約もなく探します。この「制約のなさ」が、データの性質によっては推定を不安定にする原因になります。正則化は、係数の大きさに何らかのペナルティを課すことで、この不安定さを抑え込む手法です。回帰分析の実務で正則化が役立つ場面は、大きく3つに整理できます。

• 共線性下での係数の安定化：第5章で見たとおり、説明変数どうしの相関が強いと、OLSの係数は少しのデータの変化で大きく振れてしまいます。正則化は係数の大きさにペナルティを課すことで、この振れを抑え、符号が入れ替わったり極端に大きな値になったりする現象を緩和します。
• サンプル数を上回る変数の扱い：説明変数の数\( p \)がサンプル数\( n \)に迫る、あるいは\( p > n \)となる高次元のデータでは、OLSはそもそも一意に解を求められません。正則化はこの状況でも安定した解を与えてくれる、実務上ほぼ唯一の現実的な選択肢です。
• 自動的な変数選択：後述するLasso回帰は、ペナルティのかけ方の性質上、寄与の小さい変数の係数をちょうどゼロに縮めます。人手による変数の取捨選択を、学習の過程に組み込んでしまう発想です。

3つ目の変数選択の効果はLassoに固有のものですが、1つ目と2つ目の効果はRidge・Lassoの両方に共通しています。まずはこの2つの代表的な正則化手法の仕組みを順に見ていきます。

Ridge回帰：L2正則化による係数の縮小

Ridge回帰は、通常の残差平方和に「係数の2乗和」をペナルティとして加えた式を最小化します。

\( \text{Ridge}: \quad \sum_i (y_i – \hat{y}_i)^2 + \alpha \sum_{j} \beta_j^{2} \)

\( \beta_j \)は各説明変数の係数、\( \alpha \)はペナルティの強さを決めるハイパーパラメータです。\( \alpha \)を大きくするほど係数は全体的に0へ近づいていきますが、この2乗ペナルティ（L2正則化）には、係数をちょうどゼロにするほどの力はほとんど働きません。したがってRidgeは「係数を縮める」効果はあっても「変数を捨てる」効果は基本的に持たない、という性質を押さえておく必要があります。

共線性への強さは、Ridgeの実務上の価値の中心です。直感的に説明すると、互いに強く相関する2つの変数がある場合、OLSはどちらか一方に極端に大きな係数を、もう一方に極端に小さい（場合によっては符号が逆の）係数を割り振ることで、たまたま手元のデータに最もよく適合させようとします。この割り振り方は本質的に不安定で、データが少し変わるだけで係数の大小関係が入れ替わってしまいます。Ridgeのペナルティは係数の2乗和を抑えようとするため、極端に偏った割り振りよりも、似た変数どうしに係数を均等に分配する解のほうがペナルティが小さくなり、選ばれやすくなります。結果として、相関の強い変数群には近い大きさの係数が分配され、推定全体が安定します。

共線性のある人工データで、OLSとRidgeの係数の挙動を比較してみます。

import numpy as np
from sklearn.linear_model import LinearRegression, Ridge

rng = np.random.default_rng(0)
n = 100

# x1とx2を強く相関させる(共通因子+わずかなノイズ)
common = rng.normal(size=n)
x1 = common + rng.normal(scale=0.05, size=n)
x2 = common + rng.normal(scale=0.05, size=n)
x3 = rng.normal(size=n)  # 無関係な変数
X = np.column_stack([x1, x2, x3])

# 真の関係はx1とx2の合計に依存し、係数の"真の分け方"は一意ではない
y = 2.0 * x1 + 2.0 * x2 + 0.0 * x3 + rng.normal(scale=0.3, size=n)

ols = LinearRegression().fit(X, y)
ridge = Ridge(alpha=5.0).fit(X, y)

print("OLS係数  :", np.round(ols.coef_, 2))
print("Ridge係数:", np.round(ridge.coef_, 2))

このコードを何度か乱数シードを変えて実行すると、OLSの係数はx1とx2の間で大きく偏ります。50通りのシードで試すと、x1の係数は0.95から2.81、x2の係数は1.18から3.05の範囲に散らばり、あるシードでは0.95対3.05、別のシードでは2.81対1.18というように、どちらに重みが寄るかが実行のたびに入れ替わります。真の分け方はどちらも2.0であり、合計の4.0はほぼ再現できているにもかかわらず、その内訳だけが定まらないのです。一方Ridgeの係数はx1が1.87から2.00、x2が1.89から2.01と、いずれもおおむね1.9前後の狭い範囲に収まり、乱数シードを変えても配分が大きく崩れません。x3の係数はどちらの手法でも0に近い値になります。共線性がある場面でRidgeが選ばれる理由は、まさにこの「配分の安定性」にあります。

Lasso回帰：L1正則化による変数選択

Lasso回帰は、ペナルティを係数の2乗和ではなく絶対値の和に変える手法です。

\( \text{Lasso}: \quad \sum_i (y_i – \hat{y}_i)^2 + \alpha \sum_{j} |\beta_j| \)

式の見た目はRidgeとよく似ていますが、絶対値を使うL1正則化には、Ridgeにはない重要な性質があります。ペナルティ関数の形状の違いにより、最適化の過程で寄与の小さい変数の係数がちょうどゼロに押しつぶされやすいのです。幾何学的には、L1ペナルティが作る制約領域が角のある形をしており、誤差を最小にする点がその角（つまりある係数がゼロになる点）に一致しやすいため、と説明されます。厳密な導出は専門書に譲りますが、実務上は「Lassoは係数をゼロにしてくれる」という性質だけ押さえておけば十分です。

この性質により、Lassoは推定と同時に変数選択を行っているとみなせます。説明変数が数十〜数百と多く、そのうち実際に効いているのはごく一部だろうと想定される場面（スパースな解が妥当だと考えられる場面）で威力を発揮します。手作業で1つずつ変数を吟味する代わりに、\( \alpha \)を1つ調整するだけで、寄与の小さい変数を機械的にふるい落とせるという実務上の利点があります。

Lassoと交差検証を組み合わせて、多数の候補変数の中から実際に効いている変数を絞り込む例を見てみます。標準化は交差検証の各フォールド内で行う必要があるため、Pipelineにまとめてから探索します。あわせて、\( \alpha \)を大きくしていくと各係数がどのようにゼロへ収束していくかを示す「正則化パス」も確認します。

import numpy as np
from sklearn.linear_model import Lasso, lasso_path
from sklearn.model_selection import GridSearchCV
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(1)
n_samples, n_features = 100, 30

X = rng.normal(size=(n_samples, n_features))
# 効いているのは最初の5個のみ、残り25個は無関係なノイズ
true_coef = np.zeros(n_features)
true_coef[:5] = [3.0, -2.5, 2.0, -1.5, 1.0]
y = X @ true_coef + rng.normal(scale=0.5, size=n_samples)

# 標準化とLassoをPipelineにまとめ、フォールドごとに標準化する
pipe = make_pipeline(StandardScaler(), Lasso(max_iter=10000))
param_grid = {"lasso__alpha": np.logspace(-2, 1, 30)}
search = GridSearchCV(
    pipe, param_grid, cv=5, scoring="neg_root_mean_squared_error"
)
search.fit(X, y)

best_lasso = search.best_estimator_.named_steps["lasso"]
selected = np.where(np.abs(best_lasso.coef_) > 1e-6)[0]

print("交差検証で選ばれたalpha:", round(search.best_params_["lasso__alpha"], 4))
print("非ゼロ係数が残った変数番号:", selected)
print("Lasso係数(抜粋):", np.round(best_lasso.coef_[:8], 2))

# 正則化パス:alphaを動かしたときの係数の推移(可視化用に全体を標準化)
X_scaled = StandardScaler().fit_transform(X)
alphas, coefs, _ = lasso_path(X_scaled, y, alphas=np.logspace(-2, 1, 50)[::-1])
print("パス計算に使ったalphaの個数:", len(alphas))
print("alphaが最大のときの非ゼロ係数の数:", np.sum(np.abs(coefs[:, 0]) > 1e-6))
print("alphaが最小のときの非ゼロ係数の数:", np.sum(np.abs(coefs[:, -1]) > 1e-6))

実行すると、この探索はおおむね最初の5個の変数（真に効いている変数）の係数を非ゼロのまま残し、残り25個の係数のほとんどをゼロに落とすことが確認できます。正則化パスの結果は、\( \alpha \)が大きいうちはすべての係数がゼロに押しつぶされ、\( \alpha \)を小さくしていくにつれて重要な変数から順に係数が非ゼロになっていく様子を示します。この推移を横軸\( \alpha \)（または\( \log \alpha \)）、縦軸に各変数の係数をとって描いた図が「正則化パス図」で、どの変数がどの強さのペナルティまで生き残るかを一目で確認できる、Lassoの分析でよく使われる図です。

Elastic Net：RidgeとLassoの折衷

Lassoの変数選択は便利な一方、弱点もあります。互いに強く相関する変数の集団があると、Lassoはその集団の中からどれか1つだけを選び、残りの係数をゼロにしてしまう傾向があります。しかもデータがわずかに変わるだけで、どの変数が選ばれるかが入れ替わることがあり、選ばれた変数の解釈に迷いが生じやすくなります。「広告費」と「販促費」のように意味的に近く相関の高い変数の片方だけが残り、もう片方がゼロになるといった状況です。

Elastic Netは、L1ペナルティとL2ペナルティを一定の比率で組み合わせることで、この弱点を緩和する手法です。

\( \text{ElasticNet}: \quad \sum_i (y_i – \hat{y}_i)^2 + \alpha \left( \rho \sum_{j} |\beta_j| + \frac{1-\rho}{2} \sum_{j} \beta_j^{2} \right) \)

\( \rho \)（scikit-learnではl1_ratioで指定します）がL1成分の混合比率です。\( \rho = 1 \)ならLassoそのもの、\( \rho = 0 \)ならRidgeそのものになります。L2成分を混ぜておくことで、相関の強い変数群に対してはRidgeと同様「まとめて残す」方向に働きつつ、L1成分の効果で寄与の小さい変数はゼロに落とすという、両者の折衷的な挙動になります。相関の強い変数の集団をまとめてモデルに残したいが、変数選択の効果も欲しいという場面での定番の選択肢です。

正則化の強さの選び方：交差検証と1標準誤差ルール

Ridge・Lasso・Elastic Netのいずれも、ペナルティの強さ\( \alpha \)をどう決めるかが実務上の中心的な作業になります。基本方針は前章までに扱った交差検証です。訓練データを複数のフォールドに分割し、\( \alpha \)の候補ごとに検証誤差を計算して、平均的に誤差が最も小さくなる\( \alpha \)を採用します。scikit-learnにはRidgeCV、LassoCV、ElasticNetCVが用意されており、候補のリストを渡すだけでこの探索を自動化できます。

単純に検証誤差が最小になる\( \alpha \)を選ぶ方法に加えて、「1標準誤差ルール（one-standard-error rule）」と呼ばれる考え方も知られています。これは、各\( \alpha \)における検証誤差にはフォールド間のばらつき（標準誤差）があることに着目し、「検証誤差が最小値から標準誤差1つ分の範囲に収まる\( \alpha \)の中で、最も正則化の強い（モデルが最も簡潔になる）ものを選ぶ」という方法です。誤差がほとんど変わらないのであれば、より単純で説明しやすいモデルを優先しようという発想で、Rのglmnetパッケージで広く使われている考え方です。scikit-learnの標準機能には1標準誤差ルールの自動実装はありませんが、LassoCVなどが持つmse_path_属性からフォールドごとの検証誤差を取り出せるため、平均と標準誤差を自分で計算すれば同様の選択が行えます。精度が僅差であれば、変数の少ない簡潔なモデルを選ぶという発想自体は、\( \alpha \)の自動選択に頼り切らない実務判断として持っておく価値があります。

標準化の必須性と切片の扱い（再確認）

正則化を適用するうえで欠かせない前提が、説明変数の標準化（平均0・分散1へのスケール変換）です。ペナルティは係数\( \beta_j \)の値そのものにかかるため、「年齢（20〜60程度）」と「年収（数百万〜数千万円）」のようにスケールの異なる変数が混在していると、スケールの小さい変数の係数が相対的に大きくなりやすく、同じ強さのペナルティでも変数によって実質的な影響度が変わってしまいます。標準化でスケールを揃えて初めて、\( \alpha \)による縮小を全変数に公平にかけられます。

切片（定数項）の扱いにも注意が必要です。切片は目的変数の全体的な水準を表す項であり、通常は正則化の対象に含めません。scikit-learnのRidge、Lasso、ElasticNetはいずれもfit_intercept=Trueがデフォルトで、切片は自動的にペナルティの対象外として扱われます。標準化と正則化付き回帰は必ずPipelineにまとめ、交差検証の各フォールドで標準化のパラメータを訓練部分だけから計算するようにします。検証データを含めて標準化してしまうと、検証データの情報が学習側に漏れ出す「データリーク」を招きます。

正則化回帰の係数の解釈

ここまで見てきたRidge・Lasso・Elastic Netは、いずれも実務で頻繁に使われる強力な手法ですが、係数の解釈に関しては、通常のOLSやGLMとは異なる注意が必要です。この点を誤解すると、正則化で得られた係数をそのまま「各変数の効果の大きさ」として報告してしまい、誤った意思決定につながりかねません。

正則化によって縮小された係数は、各変数が目的変数に与える効果の大きさを偏りなく推定したものではありません。ペナルティによって意図的に0へ引き寄せられた値であるため、真の効果を過小に見積もっている可能性が常にあります。また、Ridge・Lasso・Elastic Netの係数には、通常のOLSのような標準誤差やp値が原理的に定義しにくく、多くのソフトウェアでも出力されません。正則化回帰は「効果の大きさを説明する道具」ではなく、「予測精度を確保しつつ変数を絞り込むための道具」であると割り切って使うことが実務上の基本姿勢になります。

この特性を踏まえると、変数の効果を経営層に説明する目的でモデルを使う場合には、次のような2段構えの手順が実務でよく用いられます。まずLassoやElastic Netを使って多数の候補変数から寄与の大きそうな変数を絞り込み、その次に、選ばれた変数だけを使って通常のOLS（あるいはロジスティック回帰やポアソン回帰などのGLM）を改めて当てはめ直すという流れです。2段目のOLS・GLMであれば、これまでの章で扱ってきたとおり標準誤差やp値、信頼区間を通常どおり計算でき、係数の大きさや有意性を根拠をもって説明できるようになります。

1. 1段目（スクリーニング）：候補変数が多いデータに対してLasso（またはElastic Net）を適用し、交差検証で選んだ\( \alpha \)のもとで係数が非ゼロになった変数を抽出します。
2. 2段目（効果の推定）：1段目で選ばれた変数だけを使い、通常のOLSやGLMを当てはめ直します。ここで得られる係数・標準誤差・p値をもとに、効果の大きさや有意性を議論します。

この2段構えの手順には、統計学的に「post-selection inference（選択後推論）」と呼ばれる難しさが伴う点にも触れておきます。1段目で変数を選ぶ作業自体がデータを見て行われているため、2段目のOLSが出すp値や信頼区間は、あらかじめ変数を決めていた場合の理論どおりには機能しない、つまり本来より過信しやすい方向に歪む可能性があることが知られています。この歪みを補正する統計的な手法（選択的推論などの研究分野）も存在しますが、専門的な内容のため本ガイドでは立ち入りません。実務での現実的な対応としては、2段目で得られたp値や信頼区間を額面どおりに厳密なものとして扱わず、あくまで参考値として扱ったうえで、可能であれば変数選択に使ったデータとは別のデータ（ホールドアウトした検証データ）で効果を再確認する、あるいは業務知識に基づく妥当性の吟味を並行して行う、といった慎重さが求められます。

説明目的での正則化：変数のスクリーニング

これまでの内容を踏まえると、「説明のための回帰」（第3章）という観点から見た正則化の位置づけが見えてきます。候補となりうる説明変数が数十から数百に及び、どれが本当に効いているのか見当がつかない、という初期段階で、正則化は強力なスクリーニングの道具になります。業務知識だけを頼りに1つずつ変数を吟味していくのは時間がかかりますし、見落としも生じやすい作業です。Lassoによって機械的に候補を絞り込んでおき、そのうえで残った変数について、業務知識と照らし合わせながら前段の2段構えの手順でOLS・GLMによる詳細な効果検証に進む、という流れは実務上効率的です。

一方で、正則化はあくまでスクリーニングの道具であり、それ自体が「説明のための最終的なモデル」にはなりにくいという点は繰り返し強調しておく価値があります。交絡変数の調整（第3章）のように、理論的・業務的な根拠から明示的に含めるべき変数がある場合、その変数がLassoの選択でたまたま除外されてしまうと、交絡を調整しきれないモデルになってしまうおそれがあります。業務上どうしても含めるべき変数は、正則化の選択結果によらず強制的にモデルに残す、という判断も実務では珍しくありません。

ビジネスへの含意

正則化回帰は、共線性の緩和・高次元データへの対応・自動的な変数選択という3つの実務的な効果を備えた、応用範囲の広い手法です。特にLassoによる変数選択は、候補変数が多く見通しの立たない初期の分析段階で威力を発揮します。ただし、正則化された係数を効果の大きさとしてそのまま報告してしまうと、意思決定を誤らせるおそれがあります。予測やスクリーニングの道具としての正則化と、効果を説明するための道具としてのOLS・GLMを役割分担させ、選ばれた変数で改めて効果を推定し直すという2段構えの運用を徹底することが、実務での正しい使い方だといえます。

次章では、線形の直線・平面では捉えきれない非線形な関係を、解釈のしやすさを保ちながら扱う一般化加法モデル（GAM）を取り上げます。多項式回帰よりも柔軟に、かつブラックボックス化を避けながら非線形性を表現する考え方を見ていきます。

この章を深めたい方への参考書籍

『スパース回帰分析とパターン認識』（梅津佑太・西井龍映・上田勇祐、講談社）：Lasso・Elastic Netをはじめとするスパース推定の考え方を、数理的な背景から丁寧に解説した一冊です。本章で直感的に説明した「なぜLassoが係数をゼロにするのか」という仕組みを、より厳密な理論とともに深掘りしたい読者に適しています。

第11章 直線では捉えられない関係を扱うGAMと非線形回帰

前章では、説明変数が多い場面や係数が不安定になりやすい場面において、係数の大きさに罰則を課すRidge回帰・Lasso回帰・Elastic Netを取り上げ、正則化という発想がモデルを複雑にしすぎないための有効な手立てになることを確認しました。本章で扱うのは、正則化とは別の角度からモデルの複雑さに切り込む論点です。すなわち、説明変数と目的変数の関係そのものが、そもそも直線では表せない形をしている場合に、どう向き合うかという問題になります。

これまでの章で扱ってきた回帰モデルは、基本的に「説明変数が1単位増えると、目的変数はいつも同じ量だけ増える、あるいは減る」という直線的な関係を前提にしてきました。しかし実務のデータには、この前提がそもそも成り立たない関係が数多く存在します。本章では、直線を無理に当てはめたときに何を見落とすのかを確認したうえで、変数変換や多項式回帰といった素朴な対処法、そしてスプラインと一般化加法モデル（Generalized Additive Model、GAM）という、関係の形そのものを柔軟に推定する枠組みを順に取り上げます。

現実の関係は直線ばかりではない

ビジネスで扱うデータの多くは、説明変数を増やせば目的変数も比例して増える、という単純な関係にはなっていません。典型的な例を3つ挙げます。

• 広告費と売上の関係：広告費を投じ始めた当初は売上への効果が大きいものの、一定額を超えると効果の伸びは徐々に鈍化していきます。この逓減する効果を、広告費と売上が常に一定の比率で比例するという直線の式で表そうとすると、少額の広告費の効果を過小評価するか、高額の広告費の効果を過大評価するかのどちらかに偏ってしまいます
• 気温と売上の関係：飲料や空調機器の売上は、気温が低すぎても高すぎても伸び悩み、ある快適な気温帯で最も伸びるという山型の関係を持つことがあります。この関係を直線で捉えようとすると、気温が高いほど売上が増える（あるいは減る）という一方向の関係しか表現できず、山の頂点付近の実態を取りこぼしてしまいます
• 年齢と購買行動の関係：年齢と購買単価や購買頻度の関係は、若年層と高齢層でそれぞれ異なる山や谷を持つ、単調とは言えない複雑な形になることが少なくありません。年齢を1つの直線的な係数だけで説明しようとすると、年代ごとに異なる購買パターンの起伏がすべて1本の傾きに押し込められてしまいます

これらの関係に共通するのは、目的変数の変化率が説明変数の水準によって変わるという点です。線形回帰は、この変化率（係数）を全区間で一定の値と仮定するモデルであるため、変化率そのものが変化する関係には原理的に向いていません。それでも直線を無理に当てはめてしまうと、残差にはある水準の説明変数の範囲でまとまって正、別の範囲でまとまって負になるという系統的な偏りが現れます。第4章で扱った残差診断の観点からいえば、これはモデルの前提が崩れている典型的な兆候であり、見た目の決定係数がそこそこ高くても、経営判断の根拠として使うべきではない状態だといえます。加えて、こうした偏りを見逃したまま「広告費を増やせば売上は一定の比率で伸び続ける」といった説明を続けてしまうと、実際には効果が頭打ちになっている領域への追加投資を正当化してしまう恐れがあります。

素朴な対処：変数変換と多項式項

直線では捉えきれない関係に対する最も手軽な対処が、説明変数そのものを変換してから回帰に投入する方法です。代表的なのが対数変換や平方根変換で、広告費のように「増えるほど効果が鈍化する」逓減型の関係には、広告費の対数（あるいは平方根）を説明変数として使うことで、鈍化する曲線をそのまま直線の回帰式で近似できます。対数変換は値が大きくなるほど変化を圧縮する性質を持つため、逓減効果を持つ変数との相性が良く、実務でも広く使われている手法です。ただしこの方法は、単調に増加（または減少）しながら鈍化する関係にしか対応できず、気温と売上のような山型の関係には使えません。

山型のように、増加してから減少に転じる関係には、多項式項を追加する方法が使われます。説明変数の2乗項を加えれば、放物線のような山型・谷型の関係を表現でき、3乗項以上を加えればさらに複雑な波形の関係も近似できます。多項式回帰は、既存の線形回帰の枠組みに変数を追加するだけで実装できるため、導入のハードルが低いという利点があります。

しかし多項式項には無視できない限界が2つあります。1つは、次数を上げるほどデータの端（説明変数の値が最小・最大に近い領域）で予測値が不自然に暴れやすくなることです。高次の多項式は、データが密に存在する中央付近の当てはまりを良くしようとするあまり、データが疎になる端の領域で振れ幅の大きな曲線を描いてしまうことがあり、この現象はルンゲ現象として知られています。とりわけ、学習データの範囲を超えて外挿する場面では、この暴れが極端な予測値を生み出す危険があります。もう1つは、次数の選び方そのものが難しいという問題です。次数を上げれば学習データへの当てはまりは必ず良くなりますが、それが未知のデータへの予測精度の向上を意味するとは限らず、多くの場合は交差検証で次数ごとの汎化性能を比較しながら決める必要があります。しかも次数を1つ変えるだけで曲線の形が全区間にわたって変わってしまうため、ある区間だけを局所的に微調整するといった融通が利きにくい点も、多項式回帰特有の使いにくさです。

スプラインとノット：区分的に曲線をつなぐ

多項式回帰の弱点である「1本の数式で全区間を表そうとする無理」を解消する発想が、スプラインです。スプラインは、説明変数の取りうる範囲をいくつかの区間に区切り、区間ごとに低い次数（3次が使われることが多い）の滑らかな曲線を当てはめたうえで、区間の境界で曲線の値と傾きが滑らかにつながるように制約をかける手法です。この区間の境界点を、ノット（knot、節点）と呼びます。区間ごとに別々の曲線を使うため、ある区間で急激に曲がる関係があっても、その暴れが遠く離れた区間の曲線の形にまで波及しにくいという利点があります。

スプラインを使ううえで避けて通れないのが、ノットをいくつ、どこに置くかという設計上の判断です。かつてはこの判断を分析者が経験則で決めることが一般的でしたが、ノットの数が少なすぎれば曲線は単純になりすぎ、多すぎれば学習データの細かいノイズまで追いかけてしまいます。そこで現在主流になっているのが、あらかじめノットを多めに配置しておき、曲線の曲がり具合（2階微分の大きさなど）に罰則を課すことで、ノットの数そのものではなく罰則の強さによって曲線の滑らかさを調整するという考え方です。この罰則の強さは、交差検証や一般化交差検証（GCV）、あるいは後述するGAMの推定手順の中で、データから自動的に決められます。ノットの位置を人手で試行錯誤するのではなく、罰則付きの推定によって滑らかさの調整を自動化するという考え方が、次に説明するGAMの基礎になっています。

一般化加法モデル（GAM）という枠組み

スプラインによる罰則付きの曲線推定を、重回帰のように複数の説明変数を同時に扱えるモデルへと拡張した枠組みが、一般化加法モデル（GAM）です。GAMは、目的変数\( y \)を次の式で表します。

\( y = b_0 + f_1(x_1) + f_2(x_2) + \cdots + f_p(x_p) \)

線形回帰との違いは明確です。線形回帰では各説明変数に係数\( b_j \)を1つずつ掛け合わせていたのに対し、GAMでは各説明変数に、データから推定される滑らかな関数\( f_j \)を1つずつ割り当てます。この\( f_j \)がスプラインによって表現され、直線にも、山型にも、逓減する曲線にも、データが示す形に応じて自在に変化します。

ここで重要なのは、複数の関数を「足し合わせる」という加法の構造そのものは、線形回帰から変えていない点です。この構造を保っているおかげで、それぞれの変数の効果を、他の変数の効果と切り離して個別に眺めることができます。線形回帰が持つ「1つの変数の効果を単独で解釈できる」という利点を保ちながら、その効果の形を直線に限定せず、データが示す曲線の形をそのまま推定できるようにしたモデルがGAMです。

GAMの本質は、線形回帰が持つ解釈のしやすさと、非線形の関係を捉える柔軟さの間に位置を取ることにあります。各変数の効果を独立した関数として推定するため、決定木の集合体のような複雑なモデルに比べて「どの変数が、どんな形で目的変数に効いているか」を直接示しやすく、かつ直線という制約からは解放されているという、両者の良いところ取りを狙った枠組みです。

GAMの読み方：部分効果プロット

GAMを実務で使ううえで最も重要な出力物が、部分効果プロット（partial dependence plot）です。これは、ある1つの説明変数\( x_j \)の値を横軸に、その変数がもたらす効果\( f_j(x_j) \)を縦軸にとって描いた曲線で、多くの実装では推定の不確実性を示す信頼区間の帯も合わせて描画されます。線形回帰であれば「係数がいくつで、符号がプラスかマイナスか」という数値でしか変数の効果を語れませんが、GAMの部分効果プロットは、その変数が動くと目的変数がどう動くかを曲線の形でそのまま示してくれます。

この読み方が実務で強みを発揮するのは、経営層への説明の場面です。「広告費の係数は0.42で、p値は0.01未満でした」という説明よりも、「広告費が一定額を超えたあたりから効果の伸びが鈍っている、この曲線を見てください」という部分効果プロットのほうが、直感的に理解しやすく、次のアクション（広告費の追加投資をどこまで続けるべきか）にも直結しやすいという特徴があります。山型の関係や、ある閾値を境に効果が変わる関係も、統計的な専門用語を使わずに1枚の図で伝えられる点は、GAMが持つ大きな利点です。

平滑化の調整：過剰と過小のバランス

GAMの各関数\( f_j \)は、罰則付きの推定によって滑らかさが調整されます。この罰則の強さの調整には、注意すべき2つの失敗の方向があります。1つは、罰則が弱すぎて曲線が学習データの細かいノイズまで追いかけてしまう状態です。この状態になると、部分効果プロットの曲線はギザギザとした不自然な形になり、新しいデータに対する予測性能は悪化します。これは典型的な過学習の症状です。もう1つは、罰則が強すぎて曲線がほとんど直線に戻ってしまう状態です。この場合、せっかくデータに山型や逓減型の構造があっても、それをうまく拾えず、線形回帰とほとんど変わらない結果しか得られなくなります。適切な罰則の強さは、この2つの失敗の間のどこかにあり、多くの実装では交差検証や、モデルの当てはまりの良さと複雑さの釣り合いを測る指標（一般化交差検証や周辺尤度など）を使って、データから自動的に探索されます。

この調整の結果を要約する数値が、有効自由度（Effective Degrees of Freedom、EDF）です。EDFは、その変数の効果を表現するためにモデルがどれだけの「曲がる自由度」を使ったかを示す指標で、直線の回帰係数1つ分に相当する自由度が1だとすると、EDFが1に近ければその変数の効果はほぼ直線的であることを意味し、EDFが4や8のように大きくなるほど、複雑に曲がりくねった関係を表現していることを意味します。EDFは多くのGAMの実装が推定結果のサマリーとして出力するため、「この変数は本当に非線形の効果を持っているのか、それとも単に直線で十分だったのか」を確認する簡便な診断材料として使えます。

GLMの拡張としてのGAM

GAMの位置づけを正確に理解するには、第7章で扱った一般化線形モデル（GLM）の枠組みを思い出すと分かりやすくなります。GLMは、目的変数の分布（正規分布、二項分布、ポアソン分布など）を指定したうえで、線形予測子\( \eta = b_0 + b_1 x_1 + b_2 x_2 + \cdots \)を、リンク関数を介して分布の平均に結びつけるという枠組みでした。GAMは、この線形予測子\( \eta \)の中身だけを、線形の項の和から滑らかな関数の和\( \eta = b_0 + f_1(x_1) + f_2(x_2) + \cdots \)に置き換えたものです。分布の指定とリンク関数という骨組みはそのまま引き継いでいます。

この位置づけがあるため、第8章で扱ったロジスティック回帰や第9章で扱ったポアソン回帰も、そのままGAMへ拡張できます。目的変数が2値であれば二項分布とロジットリンクを組み合わせたロジスティックGAMになり、目的変数が件数であればポアソン分布と対数リンクを組み合わせたポアソンGAMになります。GAMは、GLMという枠組みを崩す新しいモデルではなく、GLMの表現力を「直線の和」から「曲線の和」に広げた拡張だと理解しておくと、これまでの章との接続がすっきりします。

使い分けの整理：GAM・変数変換・GBDT

非線形の関係に向き合う手段は、GAMだけではありません。ここまで見てきた変数変換や多項式回帰、そして機械学習編第8章で扱った勾配ブースティング決定木（GBDT）も、いずれも非線形の関係を捉える手段です。分析の目的とデータの性質に応じて、どの手段を選ぶべきかを次の表に整理させました。

この整理から導かれる実務上の使い分けは、次のように言い換えられます。関係の形があらかじめ想定でき、説明変数も少数であれば、対数変換のような単純な変数変換で十分なことが多く、モデルをわざわざ複雑にする必要はありません。関係の形は分からないが、少数の説明変数それぞれについて、どんな形で目的変数に効いているのかを把握し、経営層や現場に説明したいという目的であれば、GAMが最も適した選択肢になります。一方、説明変数の数が多く、変数同士の交互作用も豊富に存在し、個々の変数の効果よりも予測精度そのものを優先したいという場面では、GBDTのような機械学習モデルが適しています。GAMは「精度を犠牲にしてでも形を見たい」という要求に応えるモデルであり、精度の最大化そのものを目的とする場面では、機械学習編で扱う手法群と比較検討する価値があります。

説明変数が少なく、それぞれの効果の形を経営判断に使いたいならGAM、関係の形がすでに見当がついている単純な非線形なら変数変換、変数の数や交互作用が多く精度を最優先するならGBDTという3択で、目的とデータの性質から逆算して手段を選ぶという順序が、非線形の回帰分析における基本的な判断軸になります。

Pythonで確認する多項式回帰の次数比較

まず、山型の関係を持つ合成データに対して多項式回帰の次数を変えながら当てはめ、次数を上げすぎたときに何が起きるかを確認します。

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

rng = np.random.default_rng(0)
n = 200

# 気温と売上の山型(逆U字)の関係を模した合成データ
temperature = rng.uniform(0, 35, n)
sales = -0.15 * (temperature - 22) ** 2 + 500 + rng.normal(0, 15, n)

X = temperature.reshape(-1, 1)
y = sales

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0
)

for degree in [1, 2, 4, 10]:
    poly = PolynomialFeatures(degree=degree, include_bias=False)
    X_train_poly = poly.fit_transform(X_train)
    X_test_poly = poly.transform(X_test)

    model = LinearRegression().fit(X_train_poly, y_train)
    train_rmse = np.sqrt(mean_squared_error(y_train, model.predict(X_train_poly)))
    test_rmse = np.sqrt(mean_squared_error(y_test, model.predict(X_test_poly)))

    print(f"次数{degree}: 学習RMSE={train_rmse:.2f}, テストRMSE={test_rmse:.2f}")

このコードを実行すると、次数1（直線）では山型の構造を捉えきれずに学習・テストの両方のRMSEが大きいままである一方、次数2に上げると学習・テストのRMSEがそろって大きく改善します。ところが次数を4、10と上げていくと、学習RMSEはわずかに下がり続けるのに対し、テストRMSEはむしろ悪化に転じます。次数10のモデルで気温の範囲の端（0度付近や35度付近）における予測値を確認すると、実データの傾向からかけ離れた極端な値へ振れることも多く、多項式の次数を上げすぎることの危うさを具体的な数値で確認できます。

pyGAMでGAMを当てはめる

次に、GAMの当てはめをpyGAMというライブラリで実践します。pyGAMはPythonでGAMを扱うためのライブラリで、次のコマンドでインストールできます。

pip install pygam

広告費（逓減効果）と気温（山型）という、形の異なる2つの非線形関係を持つ合成データを用意し、それぞれの変数に平滑化項を割り当てたGAMを当てはめます。

import numpy as np
from pygam import LinearGAM, s

rng = np.random.default_rng(0)
n = 300

# 広告費(逓減効果)と気温(山型)から売上を生成した合成データ
ad_spend = rng.uniform(0, 100, n)
temperature = rng.uniform(0, 35, n)
sales = (
    40 * np.log1p(ad_spend)
    - 0.2 * (temperature - 22) ** 2
    + 500
    + rng.normal(0, 10, n)
)

X = np.column_stack([ad_spend, temperature])
y = sales

# 各変数に平滑化項s()を割り当て、加法モデルとして推定する
gam = LinearGAM(s(0, n_splines=20) + s(1, n_splines=20))
gam.gridsearch(X, y)  # 平滑化パラメータをGCVなどの規準に基づいて自動探索する

gam.summary()

gam.gridsearchは、平滑化の強さを左右するパラメータ（ラムダ）の候補をいくつも試し、データへの当てはまりと複雑さの釣り合いをもとに最適な値を選び出します。summaryの出力には、広告費と気温それぞれの平滑化項について、前節で説明したEDFが表示されます。広告費の項のEDFが1に近ければ対数変換に近い緩やかな効果、気温の項のEDFが3から5程度あれば、山型のようにある程度曲がった効果を推定できていると読み取れます。

部分効果を可視化する

推定したGAMから、各変数の部分効果プロットを描画します。

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
titles = ["広告費の部分効果", "気温の部分効果"]

for i, ax in enumerate(axes):
    XX = gam.generate_X_grid(term=i)
    pdep, confi = gam.partial_dependence(term=i, X=XX, width=0.95)
    ax.plot(XX[:, i], pdep)
    ax.plot(XX[:, i], confi, c="gray", ls="--")
    ax.set_title(titles[i])

plt.tight_layout()
plt.show()

generate_X_gridは、指定した変数（termで指定するインデックス）の値を等間隔に並べたグリッドを作る関数で、partial_dependenceはそのグリッド上での効果の推定値（pdep）と、指定した信頼水準（width）に対応する信頼区間（confi）を返します。このコードを実行すると、広告費の部分効果は右肩上がりでありながら傾きが徐々に緩やかになる曲線として、気温の部分効果は22度付近を頂点とする山型の曲線として、それぞれ独立に描画されます。1つの回帰式の中に混在していた2つの異なる形の非線形関係が、変数ごとの図として分離して確認できる点が、GAMを使う実務上の価値です。なお、変数の数がさらに多い場合や、目的変数が2値・件数データである場合には、pyGAMのLogisticGAMやPoissonGAM、またはstatsmodelsが提供するGLMGamクラスを使っても同様の枠組みでGAMを当てはめられます。

この章を深めたい方への参考書籍

『統計的学習の基礎 データマイニング・推論・予測』（Trevor Hastie ほか、共立出版）：一般化加法モデルを含む統計的学習の理論的な背景を、数理的な観点から体系立てて解説した書籍です。本章で扱ったスプラインの罰則付き推定やEDFの考え方を、より厳密な定式化に基づいて理解し直したい読者に向いています。

第12章 時系列データに回帰を使うときの落とし穴

前章では、直線では捉えきれない非線形な関係をGAM（一般化加法モデル）によって滑らかに表現する方法を扱いました。ここまでの11章で、単回帰から重回帰、説明のための回帰、前提診断、多重共線性、モデル選択、GLM、正則化、GAMと、回帰分析という道具の使い方を一通り整理してきました。最終章となる本章では、視点を少し変えます。実務で扱うデータの多くは、時間の順序を持つ時系列データだからです。

月次の売上、日次のアクセス数やコンバージョン率、四半期ごとの業績など、経営層やDX推進担当者が日常的に目にする数字の大半は、時間軸に沿って並んでいます。ところが、これまでの11章で扱ってきた回帰分析の枠組みは、本来「データの並び順に意味がない」ことを前提にしています。時間順に並んだデータに、この前提のままの回帰を素朴に当てはめると、見た目にはもっともらしい、しかし実態とは異なる結論に導かれることがあります。本章では、時系列データに回帰を使うときに陥りやすい問題を4つに整理し、最後に「回帰で押し切れる場面」と「専用の時系列手法に切り替えるべき場面」の見極め方を示したうえで、本コラム全体の総括を行います。

見せかけの回帰：無関係な系列が高い決定係数を示す

時系列データを扱ううえで最初に押さえておくべき現象が、見せかけの回帰（spurious regression）です。これは、本来まったく無関係な2つの時系列を回帰分析にかけただけなのに、決定係数が高く、係数も統計的に有意という結果が出てしまう現象を指します。単なる偶然の一致ではなく、トレンドを持つ時系列同士を回帰すると構造的に起こりやすい現象であることが知られています。

この現象は実際に手を動かして確かめられます。無関係な2つのランダムウォーク（乱数の累積和で作られる系列）を用意し、回帰にかけてみます。ランダムウォークは、各時点の値が「1つ前の時点の値に、平均0のランダムな増分を加えたもの」として定義される系列で、次の増分がプラスに動くかマイナスに動くかは、過去の動きとは無関係に決まります。つまり、xとyという2本のランダムウォークを完全に独立な乱数から作れば、両者の間に本来の因果関係も相関関係も存在しません。

import numpy as np
import statsmodels.api as sm

def random_walk(n, rng):
    # 乱数の累積和でランダムウォークを生成する
    return np.cumsum(rng.normal(0, 1, n))

n = 200

# x, yは互いに無関係な乱数列から作った、まったくの無関係なランダムウォーク
x = random_walk(n, np.random.default_rng(56))
y = random_walk(n, np.random.default_rng(123))

X = sm.add_constant(x)
model = sm.OLS(y, X).fit()
print(f"係数: {model.params[1]:.3f}")
print(f"p値: {model.pvalues[1]:.5f}")
print(f"決定係数(R2): {model.rsquared:.3f}")

# 乱数シードを変えて500回繰り返し、偶然の一致でないことを確認する
count_significant = 0
r2_list = []
for trial in range(500):
    xt = random_walk(n, np.random.default_rng(trial))
    yt = random_walk(n, np.random.default_rng(trial + 10000))
    m = sm.OLS(yt, sm.add_constant(xt)).fit()
    r2_list.append(m.rsquared)
    if m.pvalues[1] < 0.05:
        count_significant += 1

print(f"500回中、p値0.05未満だった割合: {count_significant / 500:.2f}")
print(f"R2の平均値: {np.mean(r2_list):.3f}")

このコードを実行すると、単発の実行でも決定係数が0.5を超えるような結果が珍しくなく、乱数シードを変えて500回繰り返した場合でも、p値が0.05を下回る（統計的に有意と判定される）割合が本来期待される5%程度どころか、大幅に高い水準になります。互いに何の関係もない2本の乱数列を回帰にかけているにもかかわらず、です。実務の現場でこの現象を知らずに月次データの回帰分析を行うと、「広告費と全く無関係な指標」の間に強い相関を見出し、誤った因果関係を報告してしまう危険があります。

見せかけの回帰が起きるのは、統計ソフトの計算が間違っているからではありません。トレンドを持つ2つの系列は、それぞれが独自の方向に緩やかに動き続けるため、たまたま似た期間だけを切り取ると、片方が上がるときにもう片方も上がって見える局面が生まれやすくなります。この「たまたま同じ方向に動いて見える期間」を、通常の回帰分析は本物の関係と区別する仕組みを持っていません。

この現象が起きる直感的な理由は、通常の回帰分析（そして決定係数やp値の計算）が、各時点の観測値を互いに独立な情報として扱っていることにあります。ランダムウォークのような系列では、ある時点の値はそれ以前の値をすべて足し合わせた結果であるため、隣り合う時点の値は強く似通っており、実質的な独立な情報量はデータの点数よりずっと少なくなります。ところが回帰分析の標準誤差やp値の計算式は、データの点数がそのまま独立な情報の数であるという前提で作られています。この前提と実態のずれが、本来存在しない関係を統計的に有意な関係として誤認させる原因です。この現象の根っこにあるのが、次に説明する定常性という概念です。

定常性：平均・分散・自己相関が時間で変わらないこと

時系列分析における定常性（stationarity）とは、ある系列の統計的な性質、具体的には平均、分散、そして時点間の自己相関（ある時点の値が、その前後の時点の値とどれだけ似た動きをするか）の構造が、観測する時期によらず一定であることを指します。定常な系列であれば、今年の1月から6月までのデータから推定した平均や変動幅は、来年の1月から6月にもおおむね当てはまると期待できます。

これに対して、実務で扱うデータの多くは非定常です。月次売上は事業の成長にともなって右肩上がりのトレンドを持ちますし、日次のアクセス数は曜日や季節によって周期的に変動する季節性を持ちます。トレンドがある系列は、時期によって平均そのものが変わり続けるため定常性の定義を満たしません。前節で見たランダムウォークも、時間が経つほど値の散らばり（分散）が大きくなっていく性質を持つため、非定常な系列の代表例です。ランダムウォークのように、差分を1回取れば定常になる型の非定常性を、その系列は単位根を持つと表現します。単位根があるかどうかを統計的に判定するADF検定などの手続きも用意されており、時系列データを回帰にかける前の点検として使われます。

定常性が重要なのは、通常の回帰分析やその背後にある統計的な検定の理論の多くが、変数が定常であることを暗黙の前提としているためです。非定常な系列同士をそのまま回帰にかけると、前節で見た見せかけの回帰のように、本来存在しない関係を検出してしまうリスクが高まります。時系列データを回帰にかける前には、まず「この系列は定常か」を確認する姿勢が欠かせません。

対処の基本：差分・成長率・トレンド項・季節ダミー

非定常な時系列にそのまま回帰を当てはめる代わりに、実務でよく使われる対処法を整理します。それぞれ意味と限界が異なるため、目的に応じて選び分ける必要があります。

差分を取るという対処が、見せかけの回帰にどれだけ効くかを実際に確認します。先ほどの2本の無関係なランダムウォークについて、水準そのものを回帰した場合と、差分系列（隣り合う時点の差）に変換してから回帰した場合を比較します。

import numpy as np
import statsmodels.api as sm

def random_walk(n, rng):
    return np.cumsum(rng.normal(0, 1, n))

n = 200
x = random_walk(n, np.random.default_rng(56))
y = random_walk(n, np.random.default_rng(123))

# 水準そのものを回帰(見せかけの回帰が起こりやすいケース)
model_level = sm.OLS(y, sm.add_constant(x)).fit()

# 差分系列(前時点との差)に変換してから回帰
dx = np.diff(x)
dy = np.diff(y)
model_diff = sm.OLS(dy, sm.add_constant(dx)).fit()

print(f"水準の回帰: R2={model_level.rsquared:.3f}, p値={model_level.pvalues[1]:.3f}")
print(f"差分の回帰: R2={model_diff.rsquared:.3f}, p値={model_diff.pvalues[1]:.3f}")

実行すると、水準そのものを回帰した場合はR2が高く、p値も小さい（統計的に有意）結果が出やすいのに対し、差分系列に変換して回帰した場合はR2がほぼ0に近づき、p値も有意水準を上回ることが確認できます。ランダムウォークの増分はそれぞれ独立な乱数そのものであるため、差分を取った時点で「無関係な系列」という実態が正しく数値に反映されるようになります。ただし、差分や成長率への変換は万能ではありません。トレンドの影響は取り除けても、次に説明する残差の自己相関という別の問題は、差分を取るだけでは解消しない場合がある点に注意が必要です。

残差の自己相関：独立性の前提が崩れる

第4章では、回帰分析が成立するための前提の1つとして、誤差項（残差）が互いに独立であるべきという条件を確認しました。時系列データでは、この前提が崩れやすいという事情があります。ある月の予測誤差が大きくプラスに振れたとき、翌月の誤差も同じ方向にプラスへ振れやすい、という具合に、時間的に近い残差同士が似た値を取る現象を残差の自己相関と呼びます。

残差の自己相関を検出するための代表的な指標が、ダービン・ワトソン統計量（Durbin-Watson statistic、以下DW統計量）です。DW統計量は0から4の範囲を取り、2に近い値であれば自己相関がほぼないことを示し、0に近いほど強い正の自己相関（隣り合う残差が同じ方向に振れる）、4に近いほど強い負の自己相関（隣り合う残差が交互に振れる）を示します。

残差に自己相関があると、実害として何が起きるのかを確認します。通常の最小二乗法は、誤差が互いに独立であるという前提のもとで標準誤差を計算します。この前提が崩れて正の自己相関がある状態でこの計算式をそのまま使うと、実際の不確実性よりも標準誤差が小さく算出されてしまいます。標準誤差が過小評価されると、本来は統計的に有意でない係数まで「有意である」と誤って判定してしまう確率が上がります。第5章で扱った多重共線性が係数の値そのものを不安定にする問題だったのに対し、残差の自己相関は係数の値は正しくても、その信頼性の評価（標準誤差やp値）だけが歪むという、性質の異なる問題です。

この歪みへの実務的な手当てが、HAC標準誤差（Heteroskedasticity and Autocorrelation Consistent standard errors、不均一分散・自己相関に頑健な標準誤差）です。代表的な計算方法にニューウェイ・ウェスト（Newey-West）法があり、残差に自己相関や分散の不均一性があっても、標準誤差をより実態に近い値に補正してくれます。係数の推定値自体は通常の最小二乗法と変わりませんが、標準誤差とp値だけがより保守的な（自己相関を考慮した）値に置き換わります。

import numpy as np
import statsmodels.api as sm
from statsmodels.stats.stattools import durbin_watson

rng = np.random.default_rng(0)
n = 150

# AR(1)的な自己相関を持つノイズを作る(前時点の80%を引き継ぐ)
noise = np.zeros(n)
for t in range(1, n):
    noise[t] = 0.8 * noise[t - 1] + rng.normal(0, 1)

trend = np.arange(n) * 0.05
y = 10 + trend + noise
x = np.arange(n)
X = sm.add_constant(x)

# 通常の最小二乗法
model_ols = sm.OLS(y, X).fit()
dw = durbin_watson(model_ols.resid)
print(f"通常のOLS: 標準誤差={model_ols.bse[1]:.4f}, DW統計量={dw:.3f}")

# HAC(Newey-West)標準誤差で補正
model_hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": 4})
print(f"HAC補正後: 標準誤差={model_hac.bse[1]:.4f}")

このコードでは、残差にAR（1）過程（1つ前の残差の80%を引き継ぐ）を仕込んでいるため、DW統計量は2から大きく離れた低い値（強い正の自己相関を示す値）になります。この状態で通常の最小二乗法とHAC標準誤差を比べると、HAC補正後の標準誤差の方が大きくなる、つまり通常の最小二乗法は標準誤差を過小評価していたことが確認できます。DW統計量が2から大きくずれている場合は、係数のp値を鵜呑みにする前に、HAC標準誤差での再計算を検討する価値があります。

ラグの見落とし：効果は当月だけに現れるとは限らない

広告費が売上に与える影響を分析する場面を考えます。素朴なモデルは「当月の広告費」を説明変数に、「当月の売上」を目的変数にした単純な回帰です。しかし広告の効果は、必ずしも支出した当月のうちに使い切られるわけではありません。認知の獲得や検討期間を経て、翌月、翌々月にも効果が持ち越されることが実務では珍しくありません。

この持ち越し効果を表現するのが、ラグ変数（過去の時点の値を説明変数として使う変数）という発想です。当月の広告費に加えて、1か月前・2か月前の広告費も説明変数に加えるモデルを分布ラグモデル（distributed lag model）と呼びます。分布ラグモデルを使うことで、「広告効果は当月に何割、翌月に何割、翌々月に何割という形で分散して現れる」という構造をデータから推定できます。当月の値だけを説明変数に使うモデルは、こうした遅れて現れる効果をすべて見落としてしまい、広告の真の効果を過小評価する結果になりがちです。

もう1つ注意が必要なのが、目的変数自身の過去の値を説明変数として使うケースです。たとえば「先月の売上」を今月の売上を説明する変数の1つとして回帰式に加えることがあります。これは季節性や慣性的な動きを捉えるうえで有効な手法ですが、この種のモデルでは残差の自己相関が診断されにくくなる、あるいは係数の推定に偏りが生じやすくなるという技術的な注意点があります。目的変数の過去の値を説明変数に含めるモデルを組む際は、通常の最小二乗法をそのまま適用してよいかどうかを慎重に検討し、必要であれば専用の推定手法（本シリーズで予定している時系列分析編で扱う自己回帰モデルなど）への切り替えを検討するべき局面です。

リークの防止：時系列の交差検証は時間順を守る

本シリーズの機械学習編（「機械学習の仕組みと使いどころがわかる実践ガイド」）では、モデルの検証設計におけるリーク（データ漏洩）の問題を扱い、時系列データに通常のランダム分割による交差検証を使うと、未来のデータで過去を予測するという実運用ではあり得ない評価になってしまうことを説明しました。これは回帰分析にもそのまま当てはまる注意点です。

第6章では、モデル選択の手段としてAICと並んで交差検証を紹介しました。売上予測や需要予測のように時間の順序を持つデータに回帰モデルを適用する場合、この交差検証の分割方法にも同じ注意が必要です。ランダムに分割する交差検証（シャッフルを有効にしたKFoldなど）を時系列データに適用すると、未来の期間のデータで学習したモデルを、過去の期間のデータで評価するという組み合わせが生まれてしまいます。実運用のモデルは常に過去のデータから未来を予測することしかできないため、この評価は実態よりも高い精度を示してしまいます。

この問題への対処として、scikit-learnのTimeSeriesSplitのように、常に時間的に前のデータを訓練用、それより後のデータをテスト用に割り当てる分割手法を使うことが原則になります。回帰モデルを需要予測や売上予測といった時系列の文脈で構築し、その性能を交差検証で評価する場合は、通常のk-Fold CVではなく時間順を守った分割を使っているかを必ず確認する必要があります。

• 目的変数・説明変数のトレンドや季節性を確認したか。非定常な系列をそのまま回帰にかけていないか
• 差分・成長率・トレンド項・季節ダミーのいずれかで、非定常性への対処を検討したか
• 残差のDW統計量を確認したか。2から大きく外れている場合、HAC標準誤差での再検証を行ったか
• 広告費や販促といった施策変数について、当月だけでなくラグ変数（過去数か月分）を検討したか
• 交差検証の分割は時間順を守っているか（TimeSeriesSplitなど）。未来のデータで過去を予測する分割になっていないか

ここから先は時系列分析の専門領域

本章で扱った対処法（差分、成長率、トレンド項、季節ダミー、HAC標準誤差、ラグ変数）は、あくまで回帰分析という枠組みの中でできる範囲の工夫です。時系列データそのものの構造をより本格的に扱うには、ARIMA（自己回帰和分移動平均モデル）や状態空間モデルといった、時系列分析に特化した手法の世界に踏み込む必要があります。ARIMAは、系列自身の過去の値（自己回帰）と過去の誤差（移動平均）、そして差分による定常化を組み合わせて将来を予測する枠組みで、状態空間モデルは、トレンドや季節性、時間とともに変化する構造そのものを確率的なモデルとして表現する、より柔軟な枠組みです。

どちらを選ぶべきかの見極めは、目的によって変わります。ある施策変数が目的変数にどれだけ影響を与えるかを説明したい場合や、説明変数の情報を活かして予測したい場合は、本章で扱った工夫を加えた回帰分析で押し切れる場面が多くあります。一方で、目的変数自身の過去の動きのパターン（周期性、トレンドの変化、変動の大きさの変化など）を精緻に捉えて将来を予測したい場合や、外部の説明変数がほとんど使えず系列自身の情報だけで予測する必要がある場合は、ARIMAや状態空間モデルといった専用の手法に切り替えることを検討すべき局面です。この領域は、本シリーズの続刊である時系列分析編で、あらためて体系的に扱う予定です。

Anagraftでは、AIプロジェクトの構想・課題設計から、データ分析・機械学習モデルの開発、AI人材の育成まで一貫したご支援を行っています。ご相談は、以下よりお問い合わせください。

お問い合わせ

本編のまとめ：回帰分析12章を貫く3つの視点

ここまで12章にわたって、単回帰から始まり、重回帰の解釈、交絡の調整、前提の診断、多重共線性、モデル選択、GLM、ロジスティック回帰、ポアソン回帰、正則化回帰、GAM、そして本章の時系列の落とし穴まで、回帰分析という道具を多角的に見てきました。個々の章で扱った手法はそれぞれ異なりますが、全体を通じて繰り返し立ち返ってきた論点は、大きく3つに整理できます。

1つ目は、説明のための回帰と予測のための回帰を区別するという視点です。第3章で導入し、多重共線性（第5章）やモデル選択（第6章）など、随所でこの区別に立ち返ってきました。係数の値そのものに意味を持たせたいのか、目的変数の値を当てることを優先するのかによって、許容できる手法や気にすべき論点が変わります。2つ目は、回帰分析が成立するための前提を診断するという姿勢です。第4章で扱った残差診断から始まり、多重共線性、そして本章の定常性や残差の自己相関まで、モデルの当てはまりの良さ（決定係数）だけを見て安心せず、前提が崩れていないかを確かめる作業が重要だと考えています。この確認を怠ると、統計的に有意に見える結果が実態とかけ離れているという事態を招きかねません。

3つ目は、データがどのような過程で生まれたかに合わせてモデルを選ぶという発想です。目的変数が件数であればポアソン回帰（第9章）、確率であればロジスティック回帰（第8章）、非線形な関係が疑われればGAM（第11章）、時間の順序に意味があれば本章で扱った工夫や専用の時系列手法というように、データの生成過程を無視して機械的に最小二乗法を当てはめるのではなく、データの性質に合わせて手法を選び分けるという考え方です。回帰分析は単一のアルゴリズムではなく、この3つの視点を持ったうえで選択し、組み合わせていく、実務のための思考の枠組みだと捉えていただければと思います。

この章を深めたい方への参考書籍

『経済・ファイナンスデータの計量時系列分析』（沖本竜義、朝倉書店）は、定常性、単位根、見せかけの回帰、ARIMA、状態空間モデルといった時系列分析の基礎を、経済・金融データを題材に体系立てて解説した一冊です。本章で駆け足に触れた内容を、より厳密な数理的背景とともに学び直したい読者、また続刊で予定している時系列分析編の予習をしておきたい読者に向いています。

終章 回帰分析を意思決定に活かすために

12章で積み上げてきたこと

本ガイドでは、単回帰から重回帰への拡張と偏回帰係数の意味に始まり、係数を正しく読み伝える作法、説明と予測という2つの目的の区別、モデルの前提と診断、多重共線性、モデル選択、そして一般化線形モデル（ロジスティック回帰・ポアソン回帰）、正則化、GAM、時系列データの落とし穴までを一続きに見てきました。

序章で掲げた3つの軸を、あらためて確認しておきます。第一に、説明と予測を区別すること。何を成果物とするかで、変数の選び方も評価の仕方も変わります。第二に、モデルの前提を診断すること。回帰は走らせれば必ず数字が出ますが、その数字が信用できるかは残差と診断が教えてくれます。第三に、データの生成のされ方に合わせてモデルを選ぶこと。0か1か、件数か、金額か。目的変数の性質に合わせて分布とリンク関数を選ぶGLMの発想は、統計モデリングの核心です。

要因分析の実務チェックリスト

回帰分析による要因分析を実務で回すときの確認事項を、本ガイドの総復習を兼ねてまとめます。

1. 問いの確認：説明が目的か、予測が目的か。説明なら、どの要因の効果を知りたいのか、交絡しうる変数は何かの見取り図を描いたか（第3章）
2. モデルの選択：目的変数の性質（連続・二値・カウント）に合った分布とリンク関数を選んだか（第7〜9章）。関係の形は直線でよいか（第11章）
3. データの確認：時系列性はないか。あるなら定常性・自己相関の手当てをしたか（第12章）
4. 推定と診断：残差プロットと影響点を確認したか（第4章）。VIFで共線性を確認したか（第5章）
5. モデルの比較：候補モデルをAICや交差検証で比較し、結論がモデルの選び方に対して頑健かを確認したか（第6章）
6. 解釈と報告：係数を単位・信頼区間つきで読み、1変数1メッセージに翻訳したか。相関を因果と言い切っていないか（第2〜3章）

回帰分析の価値は、複雑な現実を「他を一定として、この要因はこれだけ効く」という言葉に翻訳できることにあります。その翻訳が信用に足るかどうかは、モデルを走らせた後の診断と、走らせる前の設計で決まります。

次に学ぶこと

本ガイドの先には、次の領域が広がっています。

• 因果推論と効果検証編（本シリーズ続刊）：回帰による交絡調整の先へ。A/Bテスト、差の差分法、傾向スコアなど、因果を主張するための専用設計を扱います
• 時系列分析編（本シリーズ続刊）：第12章で入口だけ扱った時系列の世界。ARIMA、状態空間モデルなどの専用手法を体系的に解説します
• ベイズ統計編（計画中）：本ガイドの統計モデリングをベイズの枠組みで捉え直し、階層モデルなどへ拡張します
• 既刊「機械学習の仕組みと使いどころがわかる実践ガイド」:予測を目的とする場合のモデリング作法はこちらをご覧ください
• 既刊「データから正しく結論を導く統計学の実践ガイド」：本ガイドの前提となる検定・区間推定の基礎はこちらで扱っています

Anagraftでは、AIプロジェクトの構想・課題設計から、データ分析・機械学習モデルの開発、AI人材の育成まで一貫したご支援を行っています。ご相談は、以下よりお問い合わせください。

お問い合わせ

本コラムで紹介した書籍一覧

各章末でご紹介した参考書籍の一覧です。学習の段階に合わせてお選びください。

• 『回帰分析（新装版）』（佐和隆光、朝倉書店）：回帰分析の理論を体系的に学べる古典的名著。時代を越えて読み継がれる定番の一冊です
• 『入門はじめての多変量解析』（石村貞夫・石村光資郎、東京図書）：重回帰を含む多変量解析の手順を、数式少なめに具体例で追える入門書です
• 『調査観察データの統計科学』（星野崇宏、岩波書店）：交絡の調整と因果推論を本格的に学びたい方への定番書。第3章の先の世界です
• 『統計モデルと推測』（松井秀俊・小泉和之、講談社）：推定・検定から回帰・モデル診断までをデータサイエンス視点で整理した現代的なテキストです
• 『多変量解析入門 線形から非線形へ』（小西貞則、岩波書店）：重回帰から正則化・非線形モデルまでを橋渡しする、数理をしっかり押さえた中級書です
• 『情報量規準』（小西貞則・北川源四郎、朝倉書店）：AICを中心とするモデル選択理論の本格的な解説書。第6章の理論的背景を深められます
• 『データ解析のための統計モデリング入門』（久保拓弥、岩波書店）：GLMから階層ベイズまでの統計モデリングの考え方を丁寧に導く、広く親しまれてきた一冊です
• 『ロジスティック回帰分析』（丹後俊郎・山岡和枝・高木晴良、朝倉書店）：ロジスティック回帰の実務適用を深く掘り下げた専門書です
• 『一般化線形モデル入門（原著第2版）』（Annette J. Dobson、共立出版）：GLMの標準的教科書。分布とリンク関数の理論を体系的に固められます
• 『スパース回帰分析とパターン認識』（梅津佑太・西井龍映・上田勇祐、講談社）：Lassoを中心とするスパース推定の理論と実装を日本語で学べる貴重な一冊です
• 『統計的学習の基礎 データマイニング・推論・予測』（Trevor Hastie ほか、共立出版）：GAMやスプラインを含む統計的学習の百科事典的名著。腰を据えて取り組む価値があります
• 『経済・ファイナンスデータの計量時系列分析』（沖本竜義、朝倉書店）：見せかけの回帰や単位根を含む計量時系列分析の定番教科書。第12章の先はこの一冊です

参考情報・出典

• statsmodels公式ドキュメント https://www.statsmodels.org/
• scikit-learn公式ドキュメント https://scikit-learn.org/stable/
• pyGAM公式ドキュメント https://pygam.readthedocs.io/