statsmodelsによる線形回帰入門――OLSの推定・検定・診断・予測

線形回帰モデル

本記事はstatsmodelsの線形回帰のサンプル(Linear Regression)を翻訳し、加筆したものだ。サンプルとして

日本語	statsmodels
最小二乗法	OLS
加重最小二乗法	WLS
一般化最小二乗法	GLS
再帰的最小二乗法	Recursive LS

などがある。ここではOLSのみを扱う。

　statsmodelsはPythonというプログラミング言語上で動く統計解析ソフトである。statsmodelsのサンプルを動かすにはPCにPythonがインストールされている必要がある。まだインストールされていない方はJupyter notebookのインストールを参照。Jupyter notebookはstatsmodelsを動かすのに大変便利だ。

• 線形回帰モデル

　statsmodelsでは線形回帰モデル
$y_i=\beta_0+\beta_1 x_{1i} +\cdots+\beta_{ki}x_{ki}+ \varepsilon_i $
または行列表示では
$\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta} +\boldsymbol{\varepsilon}$
の切片($\beta_0$)と回帰係数($\beta_1,\cdots,\beta_k$)を最小二乗法により推定する。$x$は説明変数、$\varepsilon_i$は誤差である。$y$の条件付平均を、説明変数の線形結合としてモデル化している。 最小二乗推定量の性質や統計的推測を考えるため、条件の役割はつぎのように整理される。

• $E[\varepsilon_i\mid X]=0$:不偏性に関係する
• 等分散・無相関：OLSが線形な推定量のうち、偏りがなく、ばらつきが最も小さい推定量となるための条件
• 正規性：有限標本での厳密な$t$検定や$F$検定に用いる
• 完全多重共線性がない：係数を一意に推定するために必要

最小二乗法による推定

人工的に生成したデータでモデルを確認

　$x_i$を説明変数、$y_i$を被説明変数、$\varepsilon_i$を誤差項として、つぎの単回帰モデルを考える。
$
y_i=\beta_0+\beta_1x_i+\varepsilon_i
$
行列を用いると、このモデルは
$
\boldsymbol{y}
=\boldsymbol{X}\boldsymbol{\beta}
+
\boldsymbol{\varepsilon}
$
と表される。ここで、$\boldsymbol{X}$は定数項を含む説明変数行列、$\boldsymbol{\beta}$は回帰係数ベクトル、$\boldsymbol{\varepsilon}$は誤差ベクトルである。

今回は、$\beta_0=1,\qquad \beta_1=0.1$として人工的なデータを生成する。そのデータに最小二乗法による回帰モデルを当てはめ、設定した$\beta_0$と$\beta_1$をどの程度正確に推定できるかを確認する。このように、あらかじめデータを生成する仕組みを決め、そのデータにモデルを当てはめる方法は、新しいモデルを学ぶ際に有効である。真の係数が分かっているため、推定結果がどの程度妥当であるかを直接確認できるからである。

　回帰分析そのものを学ぶのが初めてではない多くの人にとって、実際の分析において、推定されたモデルがどの程度役に立つのかという疑問が常にある。その問題を明らかにしていこう。

本記事では、新たな乱数標本を生成する分析例ごとに、np.random.default_rng(シード値)によって乱数生成器を初期化する。

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

　つぎに、人工データを生成する。ここでは、その出発点として、誤差項が標準正規分布
$$
\varepsilon_i\sim N(0,1)
$$
に従う単純な場合を考える。

nsample = 100
x = np.linspace(0, 10, nsample)
X = sm.add_constant(x)
beta = np.array([1, 0.1])
y_true = np.dot(X, beta)
rng = np.random.default_rng(1001)
e = rng.normal(loc=0, scale=1, size=nsample)
y = y_true + e

　np.linspace(0, 10, nsample)によって、0から10までの範囲に100個の値を等間隔に配置する。ここで作成した$x$は確率変数ではなく、あらかじめ定められた説明変数である。

　sm.OLS()は自動的には定数項を追加しない。そのため、ここではあらかじめsm.add_constant()によって、説明変数行列$\boldsymbol{X}$に定数項を加えている。
sm.add_constant(x)は$x$に定数項の列を加える。したがって、説明変数行列$\boldsymbol{X}$は、

\boldsymbol{X}
=\begin{pmatrix}
1 & x_1\\
1 & x_2\\
\vdots & \vdots\\
1 & x_{100}
\end{pmatrix}

となる。回帰係数は、

\boldsymbol{\beta}
=\begin{pmatrix}
1\\
0.1
\end{pmatrix}

と設定している。したがって、誤差を含まない値y_trueは、
$
\boldsymbol{y}_{\mathrm{true}}
=\boldsymbol{X}\boldsymbol{\beta}
$
によって計算される。これは、説明変数$\boldsymbol{X}$を所与としたときの被説明変数の条件付き平均、
$
E[\boldsymbol{Y}\mid\boldsymbol{X}]
=\boldsymbol{X}\boldsymbol{\beta}
$
に相当する。

　つぎに、標準正規分布から誤差を生成する。
$
\varepsilon_i\sim N(0,1)
$
この誤差をy_trueに加えることで、回帰モデルに従うデータを作成する。
$
\boldsymbol{y}
=\boldsymbol{y}_{\mathrm{true}}
+
\boldsymbol{e}
=\boldsymbol{X}\boldsymbol{\beta}
+
\boldsymbol{e}
$
理論上の$\boldsymbol{Y}$は確率変数であり、コード中のyは、生成された誤差に対応する一つの実現値である。

　作成したデータに最小二乗法を適用する。

model = sm.OLS(y, X)
res = model.fit()
print(res.summary())

model.fit()によって、設定した真の係数、

$$
\beta_0=1,\qquad \beta_1=0.1
$$

を標本データから推定する。

　サマリーレポートの最上部から最初の二重線までには、使用したデータと推定したモデルに関する基本情報が表示される。ここには、被説明変数、推定方法、観測数、決定係数、$F$統計量、対数尤度、AIC、BICなどが含まれる。

　二番目の区画には、切片と回帰係数の推定結果が表示される。各係数について、推定値、標準誤差、$t$統計量、$p$値、信頼区間を確認できる。

　三番目の区画には、主として残差とモデルの診断に関する統計量が表示される。ここで用いられる残差は、
$
e_i=y_i-\hat y_i
$
である。$\hat y_i$は、推定された回帰係数を用いて計算した当てはめ値であり、
$
\hat y_i=\hat\beta_0+\hat\beta_1x_i
$
と表される。これは、$X=x_i$が与えられたときの$Y_i$の条件付き平均$E[Y_i\mid X=x_i]$の推定値である。
　理論上の誤差項は、
$
\varepsilon_i=y_i-E[Y_i\mid X=x_i]
$
であるが、母集団の回帰係数は未知であるため、$\varepsilon_i$を直接観測することはできない。そこで、実際の分析では残差$ e_i$を用いて、誤差項の分布や自己相関などについて診断する。
　サマリーレポートに表示される結果は、個別の属性を利用して別々に出力することもできる。

print('Parameters: ', res.params)
print('R2: ', res.rsquared)

fig, ax = plt.subplots(figsize=(4,2.5))
ax.plot(x, y, 'o', label="data")
ax.plot(x, y_true, 'b-', label="True")
ax.plot(x, res.fittedvalues, 'r--.', label="OLS")
ax.legend(loc='best');

Parameters:  [1.10464712 0.10143395]
R2:  0.10598781888096265

　青い直線は、人工データを生成するときに設定した真の回帰直線であり、赤い破線は標本データから最小二乗法によって推定した回帰直線である。推定された切片と傾きはそれぞれ約$1.105$と$0.101$であり、真の値$1$と$0.1$からずれている。また、決定係数は$0.106$と小さい。このずれは、人工データに加えた誤差によって生じた標本変動によるものである。

　このモデルでは条件付き平均
$
E[Y_i\mid x_i]=\beta_0+\beta_1x_i
$
が$x_i$の線形関数として表される。一方、決定係数は$R^2=0.029$と低く、説明変数$x$によって説明される$y$の変動は小さい。これは、傾きによる変化に比べて誤差のばらつきが大きいためである。

　つぎに、被説明変数$y$と説明変数$x$の分布をヒストグラムで確認する。説明変数$x$は0から10までの値を等間隔に設定しているため、各区間の度数はほぼ等しくなる。

　なお、plt.hist(X)とすると、説明変数行列$\boldsymbol{X}$の2列が同時に描かれる。第1列は定数項として追加されたすべて1の列であり、第2列が説明変数$x$である。第1列は切片そのものではなく、切片を推定するための定数項の列である。

plt.figure(figsize=(4, 2.5))
plt.hist(y)
plt.xlabel("y")
plt.ylabel("Frequency")
plt.show()

plt.figure(figsize=(4, 2.5))
plt.hist(x)
plt.xlabel("x")
plt.ylabel("Frequency")
plt.show()

---

　ここで注意が必要なのは、回帰分析において正規分布を仮定しているのは、被説明変数$Y$そのものではなく、説明変数を所与とした誤差項であるという点である。このモデルでは、

$$
Y_i\mid X=x_i\sim N(1+0.1x_i,1)
$$

である。それぞれの$Y_i$は正規分布に従うが、その平均は$x_i$によって$1$から$2$まで変化する。したがって、すべての$y_i$をまとめて作成したヒストグラムは、一つの正規分布を表しているのではなく、平均の異なる複数の正規分布を混ぜ合わせた分布を表している。

　$x$を近似的に$0$から$10$までの一様分布とみなすと、$1+0.1X$は$1$から$2$までの一様分布になる。したがって、$Y$の分布は一様分布と標準正規分布の畳み込みであり、厳密には正規分布ではない。ただし、回帰直線による平均の変化に比べて誤差のばらつきが大きいため、そのヒストグラムは正規分布に近い形に見える。

　誤差項の分布を確認するためには、$y$のヒストグラムよりも、当てはめ値を差し引いた残差

$$
e_i=y_i-\hat y_i
$$

の分布を調べる方が適切である。残差は直接観測できない誤差項$\varepsilon_i$の代理であり、残差のヒストグラムやQ–Qプロットを用いることで、誤差項の正規性からの乖離を確認できる。

---

誤差と残差の確認

人工データでは、データ生成時に使用した真の誤差$\varepsilon$を観察できる。一方、実際の回帰分析で得られるのは残差res.residだ。

fig, ax = plt.subplots(figsize=(4, 2.5))

ax.hist(e, alpha=0.5, label="Error")
ax.hist(res.resid, alpha=0.5, label="Residual")
ax.set_xlabel("Value")
ax.set_ylabel("Frequency")
ax.legend()

plt.show()

　つぎに、誤差と残差を確認する。人工データでは、データ生成時に加えた誤差$\varepsilon_i$が分かっている。一方、回帰分析によって得られる残差は、
$
e_i=y_i-\hat y_i
$
である。一般の実データでは真の誤差$\varepsilon_i$を観測できないため、残差$e_i$を用いて誤差の性質を診断する。

線形の意味について

パラメータについて線形

　これまでは、説明変数 $x$ と被説明変数 $y$ の関係が直線で表されるモデルを考えてきた。ここでは、両者の関係が曲線で表される場合を考える。

　被説明変数のノイズを除いた部分を
$
\mu(x)=\beta_0+\beta_1x+\beta_2x^2
$
とし、観測値を
$
y=\mu(x)+\varepsilon
$
とする。ここで、$\beta_0,\beta_1,\beta_2$ は係数であり、$\varepsilon$ は平均0のノイズである。$\beta_0=1,\ \beta_1=0.1,\ \beta_2=0.01$ として人工的にデータを生成し、最小二乗法によって各係数を推定する。このモデルでは、$x^2$ が含まれているため、$x$ と $y$ の関係は直線ではない。

　しかし、モデルは
$
y=\beta_0\cdot1+\beta_1\cdot x+\beta_2\cdot x^2+\varepsilon
$
と表され、各係数 $\beta_0,\beta_1,\beta_2$ はすべて1次の形で現れている。$\varepsilon$は平均0、分散0.25の正規分布にしたがうとする。

　したがって、このモデルは説明変数 $x$ については非線形であるが、パラメータ $\beta$ については線形である。このようなモデルも線形回帰モデルに含まれ、通常の最小二乗法を用いて推定される。

import numpy as np
import statsmodels.api as sm
nsample = 100
x = np.linspace(0, 10, nsample)
# 説明変数 x と x^2 を並べる
X = np.column_stack((x, x**2))
# 定数項の列を加える
X = sm.add_constant(X)
# 真の係数
beta = np.array([1.0, 0.1, 0.01])
# 平均0、標準偏差0.5の正規ノイズ
rng = np.random.default_rng(1002)
e = rng.normal(loc=0.0, scale=0.5, size=nsample)
# 人工データを生成
y = X @ beta + e
# 最小二乗法による推定
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())

推定結果を見ると、$\beta_0$ と $\beta_2$ は統計的に有意である一方、$\beta_1$ の $p$ 値は大きく、真の値が0であるという帰無仮説を棄却できない。条件数 $\mathrm{Cond.No.}$ は144となっており、やや大きな値である。これは、説明変数として $x$ と $x^2$ を同時に用いているため、両者の間に強い相関があることや、変数の尺度が異なることによる。係数の推定値や $p$ 値は、生成された乱数標本によって変化する。ただし、ここでは乱数シードを1002に固定しているため、同じコードを実行すれば同じ結果が再現される。

　つぎに、ノイズを正規分布ではなく一様分布から生成し、ノイズの分布形状が係数の推定結果にどのような違いをもたらすかを確認する。ここでは、
$e\sim U\left(-\frac{\sqrt3}{2},\frac{\sqrt3}{2}\right)$とする。このノイズの平均は0である。

nsample = 100
x = np.linspace(0, 10, nsample)
X = np.column_stack((x, x**2))
X = sm.add_constant(X)
beta = np.array([1, 0.1, 0.01])
# 平均0、分散0.25の一様乱数
rng = np.random.default_rng(1003)
e = rng.uniform(low=-np.sqrt(3) / 2,high=np.sqrt(3) / 2,size=nsample)
y_true = np.dot(X, beta)
y = y_true + e
# Fit and summary
model = sm.OLS(y, X)
res = model.fit()
print(res.summary())

係数の推定値は、正規分布の場合と同様に、真の値に比較的近い値となっている。一方、個々の係数の $p$ 値や有意性の判定は異なっている。推定結果を見ると、$\beta_0$と $\beta_2$ は統計的に有意である一方、$\beta_1$ の $p$ 値は大きく、真の値が0であるという帰無仮説を棄却できない。ただし、これは一様分布を用いたためとは限らず、$\beta_1=0.1$と$\beta_2=0.01$ と小さく設定されていることや、今回生成された乱数標本の影響も受けている。

この1回の結果だけから、正規分布の場合と一様分布の場合に差がないと結論することはできない。しかし、誤差項について$E(\varepsilon\mid X)=0$が成立していれば、誤差が正規分布にしたがわなくても、最小二乗推定量は不偏推定量となる。したがって、誤差の正規性は最小二乗推定量の不偏性に必要な条件ではない。ただし、標本の大きさが小さい場合に、通常の $t$ 検定や $F$ 検定を厳密に行うためには、誤差の正規性が必要となる。標本の大きさが大きい場合には、推定量の漸近正規性を利用できる。

• 正規誤差の例では、$\beta_1$ の真の値は0ではないものの、$x$ と $x^2$ の強い相関と標本誤差のため、$\hat{\beta}_1$ の標準誤差が相対的に大きくなり、$\beta_1=0$ を棄却できなかった。一方、$\beta_2$ は係数の数値自体は小さいが、標準誤差に対して推定値が十分に大きかったため、$\beta_2=0$ が棄却された。個々の係数の有意性は、係数の大小ではなく、推定値と標準誤差の比によって決まる。

つぎに、真の回帰曲線、観測値、および最小二乗法によって推定された回帰曲線を同じグラフに表す。観測値にはばらつきがあるものの、推定された回帰曲線が真の回帰曲線をどの程度捉えているかを視覚的に確認できる。

fig, ax = plt.subplots(figsize=(4,2.5))
ax.plot(x, y, 'o', label="data")
ax.plot(x, y_true, 'b-', label="True")
ax.plot(x, res.fittedvalues, 'r--.', label="OLS")
ax.legend(loc='best');

plt.subplots(figsize=(4,2.5))
plt.hist(y)
plt.show()

よく「原因と結果の関係が直線で表されなければ、線形回帰は使えない」と説明される。しかし、この説明だけでは、線形回帰の意味を誤解しやすい。なお、回帰分析における説明変数と被説明変数は、必ずしも原因と結果を意味しない。

　例えば、
$
y=\beta_0+\beta_1x^2+\varepsilon
$
は、$x^2$を含むため、$x$と$y$の関係は直線ではない。したがって、$x$については非線形である。一方、パラメータ$\beta_0,\beta_1$は1次の形で現れているため、パラメータについては線形であり、通常の最小二乗法で推定できる。

　statsmodelsのリファレンスにあるつぎの例は、この違いを直感的に示している。
$
\mu(x)
=\beta_0x+\beta_1\sin(x)+\beta_2(x-5)^2+\beta_3
$
$
y=\mu(x)+\varepsilon
$
ここで、
$
\beta_0=0.5,\qquad
\beta_1=0.5,\qquad
\beta_2=-0.02,\qquad
\beta_3=5
$
とし、誤差$\varepsilon$は標準正規分布にしたがうものとする。このモデルには$\sin(x)$や$(x-5)^2$が含まれているため、説明変数$x$と被説明変数$y$の関係は明らかに非線形である。それでも、
$
\mu(x)
=\beta_0f_0(x)+\beta_1f_1(x)+\beta_2f_2(x)+\beta_3f_3(x)
$
というパラメータの線形結合になっているため、OLSで推定できる。線形回帰モデルの「線形」とは、回帰曲線が直線であることではなく、パラメータが線形に現れることを意味する。

nsample = 50
sig = 0.5
x = np.linspace(0, 20, nsample)
X = np.column_stack((x, np.sin(x), (x-5)**2, np.ones(nsample)))
beta = [0.5, 0.5, -0.02, 5.]
y_true = np.dot(X, beta)
rng = np.random.default_rng(1004)
y = y_true + sig * rng.normal(loc=0, scale=1, size=nsample)
# Fit and summary:
res = sm.OLS(y, X).fit()
print(res.summary())

推定結果はまずまずであり、各回帰係数は、データ生成時に設定した値に近い値となっている。この例は、線形回帰における「線形」の意味を明確に示している。ここで線形なのは、説明変数 $x$ と被説明変数 $y$ の関係ではなく、回帰関数が
$
\mu(x)
=\beta_0x+\beta_1\sin(x)+\beta_2(x-5)^2+\beta_3
$
のように、既知の関数 $x$、$\sin(x)$、$(x-5)^2$、$1$ の線形結合として表されていることである。$\sin(x)$ や $(x-5)^2$ が含まれているため、$x$ と $y$ の関係は明らかに非線形である。それでも、各パラメータは1次の形で現れているため、通常のOLSによって推定できる。

　つぎに、推定された回帰曲線を95％予測区間とともに表示する。まず、推定結果から予測値と予測区間を数値として取り出す。

print('Parameters: ', res.params)
print('Standard errors: ', res.bse)
print('Predicted values: ', res.predict())

Parameters:  [ 0.51268542  0.63430551 -0.02195314  5.01910585]
Standard errors:  [0.02143679 0.08427054 0.00188217 0.13899696]
Predicted values:  [ 4.47027735  5.0172552   5.51555309  5.93060194  6.24030845  6.43868501
  6.53683341  6.56112088  6.548848    6.54212035  6.58093064  6.69658807
  6.90657429  7.21167119  7.59583251  8.02882013  8.47117211  8.88068622
  9.21935361  9.45960356  9.5888333   9.61147856  9.54828484  9.43289894
  9.30634095  9.21026467  9.1801141   9.23930005  9.3953538   9.63868976
  9.94418032 10.27528429 10.59005048 10.84801017 11.01682601 11.07760496
 11.02800235 10.8826071  10.67054483 10.43069257 10.2052914  10.03300633
  9.94257378  9.94808025 10.04664796 10.21890763 10.43218085 10.64584899
 10.81802679 10.91244393]

sf = res.get_prediction().summary_frame(alpha=0.05)
fig, ax = plt.subplots(figsize=(4,2.5))
ax.plot(x, y, 'o', label="Data")
ax.plot(x, y_true, 'b-', label="True")
ax.plot(x, res.fittedvalues, 'r--.', label="Predicted")
ax.plot(x, sf["obs_ci_upper"] , 'r--')
ax.plot(x, sf["obs_ci_lower"] , 'r--')
legend = ax.legend(loc="best")

$y$は複雑に変動しているが、それでも最小二乗法で推定できる。回帰関数は、$x$、$\sin(x)$、$(x-5)^2$、定数項を係数付きで加えた線形結合だからである。線形回帰は、直線だけを足し合わせる方法ではない。

最後に、$y$のヒストグラムを確認する。

plt.subplots(figsize=(4,2.5))
plt.hist(y)
plt.show()

この $y$ のヒストグラムも、回帰関数とノイズを足し合わせた結果である。回帰関数が複雑に変動すれば、$y$ の分布にもその形が反映される。したがって、$y$ のヒストグラムだけからノイズの分布を判断することはできない。

つぎに、金融データにみられる極端な変動を意識して、ノイズをコーシー分布から生成してみよう。

rng = np.random.default_rng(1005)

# コーシー分布にしたがうノイズ
e = sig * rng.standard_cauchy(size=nsample)

y = y_true + e

res_cauchy = sm.OLS(y, X).fit()
print(res_cauchy.summary())

コーシー分布は正規分布よりも裾が非常に厚く、ときどき極端に大きな値を発生させる。平均と分散が存在しないため、誤差の平均が0、分散が一定で有限であるという通常の線形回帰の仮定は成立しない。

OLS自体は計算できるが、少数の極端な観測値に強く左右され、回帰係数や標準誤差、$t$ 検定の結果は不安定になる。金融収益率も裾の厚い分布を示すことが多いが、コーシー分布はその性質を強調した極端な例である。

データの特性の把握とモデル選択

ダミー変数

性別、地域、曜日、商品の種類などは、数値の大小によって表される変数ではなく、いくつかのカテゴリーに分類される変数である。このようなカテゴリー変数を回帰分析に取り入れるには、それぞれのカテゴリーに該当するかどうかを、0と1で表す。この0と1からなる変数をダミー変数という。

例えば、観測値がA、B、C、D、Eの5つのカテゴリーに分類されている場合を考える。まず、各カテゴリーを5回ずつ並べたデータを作る。

import numpy as np
import pandas as pd
import statsmodels.api as sm
category = ["A", "B", "C", "D", "E"]
category *= 5
category[:10]

この段階では、AからEまでが繰り返されている。カテゴリーごとに並べるために、データをソートする。

category = np.asarray(sorted(category))
category

array(['A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'C', 'C', 'C',
       'C', 'C', 'D', 'D', 'D', 'D', 'D', 'E', 'E', 'E', 'E', 'E'],
      dtype='<U1')

これをダミー変数に変換する。

design = pd.get_dummies(category, dtype=float)
design.drop_duplicates().head(5)

get_dummies()によって、AからEまでのカテゴリーが、それぞれ一つの列に変換される。カテゴリーAに属する観測ではA列が1となり、その他の列は0となる。カテゴリーBに属する観測ではB列が1となり、その他の列は0となる。このように、文字で表されたカテゴリーが、0と1からなる複数の数値変数に置き換えられる。

　つぎに、カテゴリーAからEに、それぞれ1から5までの値を対応させる。

y = np.repeat([1, 2, 3, 4, 5], 5) 
y

array([1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5,
       5, 5, 5])

となる。ここで、ダミー変数だけを説明変数として回帰を行う。

X = design
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())

この例は $y$ をカテゴリーから決定論的に作った完全適合の人工例であり、$R^2=1$、標準誤差はほぼゼロ（表示上の $10^{-16}$ は数値誤差）となる。」

このモデルには定数項を加えていない。そのため、AからEまでのすべてのダミー変数を説明変数として使用できる。

回帰モデルは、

$$
y_i
=\beta_A D_{Ai}
+\beta_B D_{Bi}
+\beta_C D_{Ci}
+\beta_D D_{Di}
+\beta_E D_{Ei}
+\varepsilon_i
$$

と表される。

$D_{Ai}$は、観測値$i$がカテゴリーAに属するときに1、そうでないときに0となる。同じように、$D_{Bi}$から$D_{Ei}$までが、それぞれのカテゴリーへの所属を表す。

推定される係数は、

$$
\hat\beta_A=1,\quad
\hat\beta_B=2,\quad
\hat\beta_C=3,\quad
\hat\beta_D=4,\quad
\hat\beta_E=5
$$

となる。

例えば、カテゴリーCに属する観測では、

$$
D_A=0,\quad D_B=0,\quad D_C=1,\quad D_D=0,\quad D_E=0
$$

であるから、予測値は、

$$
\hat y
=1\times0
+2\times0
+3\times1
+4\times0
+5\times0
=3
$$

となる。

この例では、各ダミー変数の係数が、それぞれのカテゴリーに対応する$y$の平均を表している。

ダミー変数を使うことによって、A、B、Cのような数値ではないカテゴリーも、回帰モデルの中で表現できることが分かる。

ここではダミー変数だけを説明変数として使用した。実際の回帰分析では、連続変数とダミー変数を同時に使用することが多い。つぎに、statsmodelsの公式例を用いて、その方法を確認する。

plt.subplots(figsize=(4,2.5))
plt.plot(y)
plt.plot(results.fittedvalues)

statsmodelsのサンプル

このモデルで使用している説明変数行列は、列の順に
$$
X=(x,\ d_1,\ d_2,\ 1)
$$
であり、係数ベクトルは
$$
\beta=(1,\ 3,\ -3,\ 10)^\top
$$
である。したがって、ノイズを加える前の値は
$$
y_{\mathrm{true}}
=1x+3d_1-3d_2+10
$$
となる。
ここで、グループ0を基準グループとしているため、グループ0に対応するダミー変数は説明変数行列から除かれている。各グループの回帰式はつぎのようになる。

　グループ0では$d_1=0,\ d_2=0$なので、
$$
y_{\mathrm{true}}=x+10
$$
となる。

　グループ1では$d_1=1,\ d_2=0$なので、
$$
y_{\mathrm{true}}=x+13
$$
となる。ダミー変数$d_1$によって、グループ0より切片が3だけ高くなる。

　グループ2では$d_1=0,\ d_2=1$なので、
$$
y_{\mathrm{true}}=x+7
$$
となる。ダミー変数$d_2$によって、グループ0より切片が3だけ低くなる。

したがって、このモデルは傾きがすべて1で共通し、切片だけが
$$
10,\quad 13,\quad 7
$$
と異なる3本の直線から構成されている。ダミー変数は、グループごとに回帰直線を上下に移動させる役割を果たしている。
推定結果では、係数はおおむねデータ生成時に設定した
$$
(1,\ 3,\ -3,\ 10)
$$
に近い値となっている。ただし、標本に標準正規乱数を加えているため、推定値は設定値と完全には一致しない。

　この例から、ダミー変数を使うことで、数値では直接表しにくいグループの違いを、基準グループからの切片の差として回帰モデルに組み込めることが分かる。
　さらにダミー変数の役割が直感的に理解できるように、statsmodelsのサンプルを見てみよう。

　50個の標本を生成する。groupsに0,1,2を設定する。それをもとにダミー変数を生成する。
説明変数xを0から20までの数値を50個に分割して生成する。それに３つのグループの内の最初のgroupsを取り除き加え、かつ定数(切片)を付加する。

　$\beta$を1,3,-3,10としてy_trueを生成する。ノイズは標準正規乱数とする。

nsample = 50
groups = np.zeros(nsample, int) # data typeを整数(int)に指定
groups[20:40] = 1
groups[40:] = 2
dummy = (groups[:,None] == np.unique(groups)).astype(float)
#dummy = sm.categorical(groups, drop=True)
x = np.linspace(0, 20, nsample)
# drop reference category
X = np.column_stack((x, dummy[:,1:])) # 2次元配列に変換
X = sm.add_constant(X, prepend=False) #定数は最初の列(True)ではなく、最後の列（False)に加える
beta = [1., 3, -3, 10]
y_true = np.dot(X, beta)
rng = np.random.default_rng(1006)
e = rng.normal(loc=0, scale=1, size=nsample)
y = y_true + e

# Inspect the data:最初の５つのデータを確認する。
print(X[:5,:])
print(y_true[:5])
print(y[:5])
print(groups)
print(dummy[:5,:])

[[0.         0.         0.         1.        ]
 [0.40816327 0.         0.         1.        ]
 [0.81632653 0.         0.         1.        ]
 [1.2244898  0.         0.         1.        ]
 [1.63265306 0.         0.         1.        ]]
[10.         10.40816327 10.81632653 11.2244898  11.63265306]
[11.16897556  8.41026794 11.69060646 13.02396705 12.27673746]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 2 2 2 2 2 2 2 2 2 2]
[[1. 0. 0.]
 [1. 0. 0.]
 [1. 0. 0.]
 [1. 0. 0.]
 [1. 0. 0.]]

２つのダミー変数の影響のある部分と無い部分の3つのグループから構成されるy_trueを可視化する。

fig, ax = plt.subplots(figsize=(8,2.5))
plt.plot(X[:,0],label='X')
plt.plot(X[:,1],label='dummy1')
plt.plot(X[:,2],label='dummy2')
plt.plot(X[:,3],label='const')
plt.plot(y_true,label='y_true')
plt.legend()
plt.show()

$X$と$\beta$の関係が見て取れる。

# Fit and summary:
res2 = sm.OLS(y, X).fit()
print(res2.summary())

結果は良好であるのでグラフに描いてみよう。ダミー変数の利用により、３つのグループの変数が区別されている。　このグラフでは、青線がノイズを加える前の真の値$y_{\mathrm{true}}$、赤の破線がOLSによる推定値である。両者がおおむね一致していることから、回帰モデルがデータ生成構造を適切に捉えていることが分かる。

　上下の赤い破線は予測区間である。これは平均的な回帰直線の不確実性だけでなく、個々の観測値に含まれるノイズも考慮しているため、推定された回帰直線の周囲に一定の幅を持つ。

sf = res2.get_prediction().summary_frame(alpha=0.05)

fig, ax = plt.subplots(figsize=(4,2.5))
ax.plot(x, y, 'o', label="Data")
ax.plot(x, y_true, 'b-', label="True")
ax.plot(x, res2.fittedvalues, 'r--.', label="Predicted")
ax.plot(x, sf["obs_ci_upper"] , 'r--')
ax.plot(x, sf["obs_ci_lower"] , 'r--')
legend = ax.legend(loc="best")

結合仮説の検定

　回帰分析では、説明変数や観測値をどのように選び、どのようなモデルを構成するかが重要である。モデルを比較する基準の一つに、サマリーレポートにも表示される自由度調整済み決定係数がある。一般には、自由度調整済み決定係数が大きいほど、説明力とモデルの複雑さのバランスがよいと考えられる。ただし、モデル選択には自由度調整済み決定係数以外にも、F検定、AIC、BICなど、さまざまな基準がある。ここでは、複数の回帰係数を同時に検定する結合仮説の検定について見ていく。

　統計的検定の基本については、つぎの記事も参照されたい。

統計的検定について

　統計的検定を理解するには、ある程度の時間と反復が必要である。焦らず、一つずつ仮説の意味を確認していくことが重要である。

F検定

　線形回帰のサマリーレポートに表示されるF検定は、通常、切片を除くすべての回帰係数がゼロであるという帰無仮説を検定している。今回のモデルは、
$
y=\beta_1x+\beta_2d_1+\beta_3d_2+\beta_4+\varepsilon
$
である。ここで、$d_1$と$d_2$はグループを表すダミー変数である。ここでは、2つのダミー変数の係数が同時にゼロであるという帰無仮説を検定する。
$
H_0:\beta_2=0,\quad \beta_3=0
$
対立仮説は、
$
H_1:\beta_2,\beta_3の少なくとも一方は0ではない
$
である。この帰無仮説は、行列を使うと、
$
R\beta=0
$
と表せる。係数ベクトルを
$
\beta=
\begin{pmatrix}
\beta_1\
\beta_2\
\beta_3\
\beta_4
\end{pmatrix}
$
とすると、制約行列$R$は、

R=
\begin{pmatrix}
0&1&0&0\\
0&0&1&0
\end{pmatrix}

となる。第1行の$[0,1,0,0]$は$\beta_2=0$という制約を表し、第2行の$[0,0,1,0]$は$\beta_3=0$という制約を表している。

R = [[0, 1, 0, 0],
     [0, 0, 1, 0]]

print(np.array(R))
print(res2.f_test(R))

実行結果は次のとおりである。

[[0 1 0 0]
 [0 0 1 0]]
<F test: F=122.7344666459456, p=3.612500226461874e-19, df_denom=46, df_num=2>

$p$値はきわめて小さいため、
$
H_0:\beta_2=\beta_3=0
$
は棄却される。したがって、$x$の値を一定としたとき、グループ間の切片の差を表す2つのダミー変数は、全体として回帰モデルに有意な情報を加えていると判断される。

　これは、3つのグループの条件付き平均がすべて同じであるという仮説が棄却されたことを意味する。より正確には、基準グループと比較した2つのグループ差の少なくとも一方がゼロではないと判断される。

　同じ仮説は、数式に近い形式でも指定できる。

print(res2.f_test("x2 = 0, x3 = 0"))

<F test: F=122.7344666459456, p=3.612500226461874e-19, df_denom=46, df_num=2>

　また、statsmodelsでは次のように記述することもできる。

print(res2.f_test("x2 = x3 = 0"))

<F test: F=122.7344666459457, p=3.612500226461817e-19, df_denom=46, df_num=2>

　結果はいずれもほぼ同じである。

　個々の$t$検定は、それぞれのダミー変数の係数を別々に検定する。一方、F検定では、複数の係数が同時にゼロであるかを検定できる。このため、ダミー変数全体がモデルに必要であるかを判断する場合には、結合仮説のF検定が適している。

グループ効果が小さい場合

　statsmodelsのリファレンスにあるこの例も非常に良い例である。Small group effects とはグループを表すダミー変数の係数が小さいという意味だ。よって、あるグループのダミー変数の係数を小さくしてみよう。この問題は仮説検定と密接な関係にある。ここでは前節のF検定と回帰係数のp値を用いて考えてみよう。

　実際のデータを生成する際に$\beta_2=0.3$を小さくし、$\beta_3=0.0$としてみる。回帰係数を見てみると、$\beta_2$は0.3なので帰無仮説のゼロに近くなる。$\beta_3$はゼロであるので帰無仮説と同じである。$\beta_2=0.3$はゼロではないため、理想的には帰無仮説$H_0:\beta_2=0$が棄却され、$\beta_3=0$については棄却されない結果を期待したくなる。しかし、1回の標本から必ずそのような結果が得られるわけではない。

beta = [1., 0.3, 0, 10]
y_true = np.dot(X, beta)
rng = np.random.default_rng(1006)
e = rng.normal(loc=0, scale=1, size=nsample)
y = y_true + e 
res3 = sm.OLS(y, X).fit()
print(res3.summary())

X2が$\beta_2$、X3が$\beta_3$であるので、どちらもp値は0.1よりも大きい。これでは帰無仮説を棄却できない。$\beta_2=0.3$はゼロではないため、帰無仮説$H_0:\beta_2=0$は実際には正しくない。しかし、係数0.3は誤差の大きさに比べて小さいため、この標本からは帰無仮説を棄却できなかった。一方、$\beta_3=0$についても帰無仮説は棄却されなかった。ただし、これらは1回の標本についての結果であり、真の係数だけからp値の大小が一意に決まるわけではない。

　つぎにF検定でみてみよう。サマリーレポートのF検定のp値は3.25e-38と非常に小さい。したがって、3つの回帰係数がゼロであるという帰無仮説は棄却されている。全体のF検定が棄却された主な理由は、$x$の係数が1であり、明確にゼロではないからである。したがって、この結果からダミー変数によるグループ効果が存在するとは結論できない。これをチェックしてみよう。

R = [[1, 0, 0, 0],[0, 1, 0, 0], [0, 0, 1, 0]]
print(res3.f_test(R))

<F test: F=688.5799991188836, p=3.2540012604738667e-38, df_denom=46, df_num=3>

p値は一致し、帰無仮説を棄却するという結論は同じである。つぎに1つ1つ見ていこう。

R = [[0, 1, 0, 0]]
print(res3.f_test(R))

<F test: F=0.21696237096650858, p=0.6435622766314202, df_denom=46, df_num=1>

R = [[0, 0, 1, 0]]
print(res3.f_test(R))

<F test: F=0.5166886961263107, p=0.47589330883667147, df_denom=46, df_num=1>

どちらも回帰係数のｐ値と同じという結果になった。1つの回帰係数について検定する場合、通常の線形回帰ではF統計量とt統計量の間に$F=t^2$という関係があるため、両者のp値は一致する。

　つぎに、$\beta_2$と$\beta_3$が同時にゼロであるという帰無仮説をF検定によって調べてみよう。

R = [[0, 1, 0, 0], [0, 0, 1, 0]]
print(res3.f_test(R))

<F test: F=0.314084681830452, p=0.7320112031900305, df_denom=46, df_num=2>

そうすると先ほどのF検定の結果とは異なりp値は0.732である程度の大きさをもっている。5％または10％の有意水準を上回っている。したがって、$\beta_2=\beta_3=0$という帰無仮説は棄却できない。

　仮説検定は、標本から帰無仮説の真偽を断定する手続きではない。あらかじめ第1種の過誤の確率を有意水準以下に抑えながら、帰無仮説に反する証拠が十分に得られたかどうかを判断する手続きである。そのため、帰無仮説を棄却できなかった場合には、通常、「帰無仮説を採択する」とは表現せず、「帰無仮説を棄却するだけの十分な証拠が得られなかった」と表現する。

　ここでの問題の原因は仮説検定の本質的な特性にあるのだ。それはつぎのように説明される。帰無仮説を平均値$\mu=50$としよう。図に書くと

というイメージだ。グレーの部分が棄却域である。実際の分布の$\mu=48$としてみよう。

そうすると棄却域のグレーの分は少し増える。しかし、帰無仮説を棄却したければ、$\mu$をもっとズラせばよいのである。そこで$\mu=56$にしてみよう。

詳しくはstatsmodelsによる統計的仮説検定 入門
を参照してほしい。

多重共線性

　Longley（1967）のデータは、説明変数間に強い多重共線性をもつことで知られている。多重共線性とは、ある説明変数が、ほかの説明変数の線形結合によってかなり正確に表される状態である。

　多重共線性が存在しても、通常の仮定のもとでは最小二乗推定量に偏りが生じるわけではない。しかし、回帰係数の標準誤差が大きくなり、係数の符号や大きさが、観測値やモデルのわずかな変更によって大きく変化することがある。そのため、個々の回帰係数を安定的に解釈することが難しくなる。

Longleyデータの概要はつぎのとおりである。

観測値の数 ― 16
説明変数の数 ― 6

変数名と内容：

    TOTEMP  - 総雇用者数（目的変数）
    GNPDEFL - GNPデフレーター
    GNP     - GNP
    UNEMP   - 失業者数
    ARMED   - 軍人の数
    POP     - 人口
    YEAR    - 年（1947年から1962年）

from statsmodels.datasets.longley import load_pandas

data = load_pandas()
y = data.endog
X = data.exog
X = sm.add_constant(X)

ols_model = sm.OLS(y, X)
ols_results = ols_model.fit()

print(ols_results.summary())

　メッセージは条件数から多重共線性を警告している。
サマリーレポートには、条件数が非常に大きいという警告が表示される。

　条件数が大きいことは、説明変数間に強い多重共線性が存在する可能性、または変数の尺度の違いなどによる数値計算上の問題が存在する可能性を示している。ただし、条件数が大きいという結果だけで、多重共線性の原因や程度を完全に判断することはできない。

条件数（Condition number）

　多重共線性を調べる方法の1つに条件数がある。statsmodelsのリファレンスでは、Greeneを参照し、条件数が20を超える場合には多重共線性に注意する必要があるとしている。ただし、この基準は変数の正規化方法や条件数の定義にも依存するため、絶対的な判定基準ではない。

　条件数は、数値計算上の問題の解きやすさを表す尺度である。条件数が小さい問題は良条件（well-conditioned）、条件数が大きい問題は悪条件（ill-conditioned）と呼ばれる。

　回帰分析における条件数は、説明変数行列$X$の最大特異値と最小特異値の比として表される。
$
\kappa(X)
=\frac{\sigma_{\max}(X)}
{\sigma_{\min}(X)}
$
　これは、$X^\mathsf{T}X$の最大固有値と最小固有値を用いると、
$
\kappa(X)
=\sqrt{
\frac{\lambda_{\max}(X^\mathsf{T}X)}
{\lambda_{\min}(X^\mathsf{T}X)}
}
$
と表すことができる。

　Longleyデータについて、定数項以外の説明変数の長さを1に正規化して条件数を計算する。

norm_x = X.to_numpy(dtype=float, copy=True)

for i, name in enumerate(X.columns):
    if name == "const":
        continue
    norm_x[:, i] = norm_x[:, i] / np.linalg.norm(norm_x[:, i])

norm_xtx = norm_x.T @ norm_x

eigs = np.linalg.eigvalsh(norm_xtx)
condition_number = np.sqrt(eigs.max() / eigs.min())

print(condition_number)

56240.870732127354

　条件数は約56,241であり、非常に大きい。この結果は、説明変数行列が線形従属に近く、回帰係数の推定が不安定になる可能性が高いことを示している。

　多重共線性が強い場合、データやモデルをわずかに変更しただけでも、回帰係数が大きく変化することがある。これを確認するため、16個の観測値のうち最後の2個を除き、14個の観測値でモデルを再推定する。

# 最後の2つの観測値を除いて再推定
ols_results2 = sm.OLS(y.iloc[:14], X.iloc[:14]).fit()

percentage_change = (
    (ols_results2.params - ols_results.params)
    / ols_results.params
    * 100
)

for name, change in percentage_change.items():
    print(f"{name}: {change:.2f}%")

const: 4.55%
GNPDEFL: -105.20%
GNP: -3.43%
UNEMP: 2.92%
ARMED: 3.32%
POP: 97.06%
YEAR: 4.64%

　観測値を2つ除いただけで、GNPデフレーターの係数は約105％、人口の係数は約97％変化している。この結果は、多重共線性が強いモデルでは、標本のわずかな変更によって個々の回帰係数が大きく変化することを示している。

　ただし、多重共線性が存在するからといって、直ちに説明変数を削除すべきであるとは限らない。予測が目的である場合には、個々の係数が不安定でもモデル全体として一定の予測力をもつ可能性がある。一方、個々の係数の解釈が目的である場合には、類似する説明変数の統合、変数の選択、追加データの取得などを検討する必要がある。

影響力のある観測値

DFBETAS

　先ほどは最後の2つの観測値を同時に除いて、回帰係数の変化を調べた。DFBETASでは、観測値を1つずつ除いてモデルを再推定し、それぞれの回帰係数がどの程度変化するかを調べる。

　観測値$i$を除いたときの係数を$\hat{\beta}_{j(i)}$とすると、DFBETASは概念的には、

\mathrm{DFBETAS}_{ij}
=\frac{
\hat{\beta}_j-\hat{\beta}_{j(i)}
}{
\hat{\beta}_jの標準誤差
}

と表される。分母の正確な定義には、観測値$i$を除いて推定した誤差分散が用いられる。

　したがって、DFBETASの絶対値が大きいほど、その観測値を除いたときに回帰係数が大きく変化する。符号が正であれば、全標本による係数が、その観測値を除いた場合の係数より大きいことを示す。符号が負であれば、その逆である。一般的な目安として、
$
|\mathrm{DFBETAS}_{ij}|>\frac{2}{\sqrt{N}}
$
を満たす観測値は、係数$j$に対して強い影響をもつ可能性があると考える。

　Longleyデータでは観測値数が$N=16$であるため、
$
\frac{2}{\sqrt{16}}=0.5
$
が基準となる。

infl = ols_results.get_influence()

dfbetas = infl.summary_frame().filter(regex="^dfb_")

threshold = 2 / np.sqrt(ols_results.nobs)

print("基準値:", threshold)
print(dfbetas)

基準値: 0.5
    dfb_const  dfb_GNPDEFL   dfb_GNP  dfb_UNEMP  dfb_ARMED   dfb_POP  dfb_YEAR
0   -0.016406    -0.234566 -0.045095  -0.121513  -0.149026  0.211057  0.013388
1   -0.020608    -0.289091  0.124453   0.156964   0.287700 -0.161890  0.025958
2   -0.008382     0.007161 -0.016799   0.009575   0.002227  0.014871  0.008103
3    0.018093     0.907968 -0.500022  -0.495996   0.089996  0.711142 -0.040056
4    1.871260    -0.219351  1.611418   1.561520   1.169337 -1.081513 -1.864186
5   -0.321373    -0.077045 -0.198129  -0.192961  -0.430626  0.079916  0.323275
6    0.315945    -0.241983  0.438146   0.471797  -0.019546 -0.448515 -0.307517
7    0.015816    -0.002742  0.018591   0.005064  -0.031320 -0.015823 -0.015583
8   -0.004019    -0.045687  0.023708   0.018125   0.013683 -0.034770  0.005116
9   -1.018242    -0.282131 -0.412621  -0.663904  -0.715020 -0.229501  1.035723
10   0.030947    -0.024781  0.029480   0.035361   0.034508 -0.014194 -0.030805
11   0.005987    -0.079727  0.030276  -0.008883  -0.006854 -0.010693 -0.005323
12  -0.135883     0.092325 -0.253027  -0.211465   0.094720  0.331351  0.129120
13   0.032736    -0.024249  0.017510   0.033242   0.090655  0.007634 -0.033114
14   0.305868     0.148070  0.001428   0.169314   0.253431  0.342982 -0.318031
15  -0.538323     0.432004 -0.261262  -0.143444  -0.360890 -0.467296  0.552421

　DFBETASの絶対値が0.5を超える観測値を抽出してみよう。

flags = dfbetas.abs() > threshold

years = load_pandas().exog["YEAR"].astype(int)

influential = pd.DataFrame({
    "YEAR": years[flags.any(axis=1)],
    "影響を受ける係数": flags[flags.any(axis=1)].apply(
        lambda row: ", ".join(
            row.index[row].str.replace("dfb_", "", regex=False)
        ),
        axis=1
    )
})

print(influential)

    YEAR                       影響を受ける係数
3   1950                 GNPDEFL, GNP, POP
4   1951    const, GNP, UNEMP, ARMED, POP, YEAR
9   1956           const, UNEMP, ARMED, YEAR
15  1962                       const, YEAR

1950年、1951年、1956年、1962年の観測値は、少なくとも1つの回帰係数についてDFBETASの絶対値が0.5を超えている。なかでも1951年の観測値は、多くの回帰係数に強い影響を与えている。1951年の主なDFBETASはつぎのとおりである。

\begin{aligned}
\mathrm{DFBETAS}_{\mathrm{const}}&=1.871,\\
\mathrm{DFBETAS}_{\mathrm{GNP}}&=1.611,\\
\mathrm{DFBETAS}_{\mathrm{UNEMP}}&=1.562,\\
\mathrm{DFBETAS}_{\mathrm{ARMED}}&=1.169,\\
\mathrm{DFBETAS}_{\mathrm{POP}}&=-1.082,\\
\mathrm{DFBETAS}_{\mathrm{YEAR}}&=-1.864.
\end{aligned}

つまり、1951年の観測値を1つ除くだけで、定数項、GNP、失業者数、軍人の数、人口およびYEARの係数が、それぞれの標準誤差と比較して大きく変化する。1956年の観測値も、定数項、失業者数、軍人の数およびYEARの係数に強い影響を与えている。また、1962年の観測値は、特に定数項とYEARの係数に影響を与えている。このように、Longleyデータでは特定の観測値が1つの係数だけでなく、複数の係数に同時に影響を与えている。これは、説明変数が互いに強く関連しているため、1つの観測値の変化が複数の係数に配分されるからである。ただし、DFBETASが基準値を超えたからといって、その観測値を機械的に削除してよいわけではない。DFBETASは、データ入力の誤り、特殊な事情、モデルの設定などを詳しく調べる必要がある観測値を示す診断指標である。

　Longleyデータでは、それぞれの観測値が1947年から1962年までの各年に対応している。そのため、観測値を除くことは、単に異常値を取り除くことではなく、特定の年の経済情報を分析から除外することを意味する。

Cookの距離

　DFBETASは、観測値が個々の回帰係数に与える影響を測定する。これに対してCookの距離は、1つの観測値が回帰モデル全体に与える影響をまとめて測定する指標である。
summary_frame()にはDFBETASだけでなく、Cookの距離、標準化残差、レバレッジ、DFFITSなども含まれている。

influence_frame = infl.summary_frame()

result = pd.DataFrame({
    "YEAR": years,
    "Cook_distance": influence_frame["cooks_d"]
})

print(result)

　Cookの距離については、経験的な目安として、
$
D_i>\frac{4}{N}
$
を用いることがある。

　Longleyデータでは、
$
\frac{4}{16}=0.25
$
である。

cook_threshold = 4 / ols_results.nobs

print("基準値:", cook_threshold)

print(
    result[
        result["Cook_distance"] > cook_threshold
    ]
)

    YEAR  Cook_distance
4   1951       0.6139
15  1962       0.4667

　この基準では、1951年と1962年の観測値が、モデル全体に強い影響を与える可能性がある。

　1950年と1956年のCookの距離もそれぞれ約0.244と0.235であり、基準値0.25にはわずかに届かないものの、注意すべき観測値である。Cookの距離の基準も絶対的なものではなく、観測値を詳しく調べるための目安として利用する。

　DFBETASとCookの距離の結果をまとめると、1951年の観測値は個々の係数とモデル全体の両方に強い影響を与えている。1962年の観測値もモデル全体への影響が大きい。一方、1950年と1956年は、モデル全体よりも特定の回帰係数に対する影響が大きい。

　Longleyデータは観測値が16個しかない一方で、定数項を含めて7個の係数を推定している。そのため、残差の自由度は9しかなく、1つの観測値が推定結果に与える影響が相対的に大きい。

　さらに、多重共線性によって個々の係数が不安定になっている。この例から分かるのは、非常に高い決定係数が得られていても、個々の回帰係数が安定しているとは限らないということである。

$$
R^2が高い
\quad\not\Rightarrow\quad
個々の回帰係数が安定している
$$

　したがって、回帰モデルを評価する際には、決定係数やp値だけでなく、条件数、DFBETAS、Cookの距離なども合わせて確認する必要がある。

回帰係数の推定と予測――標本内の当てはめと標本外予測

　現実のデータについて母集団回帰関数を直接知ることはできないため、観測された標本を用いて標本回帰関数を推定する。

　線形回帰モデルにおいて、最小二乗推定量が最良線形不偏推定量、すなわちBLUEとなるためには、主として次の条件が必要になる。

1. 回帰モデルがパラメータについて線形である
2. 説明変数行列が完全な列ランクをもち、完全な多重共線性が存在しない
3. 説明変数を所与とした誤差項の条件付き期待値がゼロである
4. 誤差項の条件付き分散が一定である
5. 異なる観測値の誤差項の間に相関がない

なお、誤差分散の推定や統計的検定を行うため、観測値の数が推定するパラメータの数より多いことが必要となる。

　数式で表すと、特に重要なのは
$$
E(\varepsilon\mid X)=0
$$
および
$$
\operatorname{Var}(\varepsilon\mid X)=\sigma^2\boldsymbol{I}_n
$$
である。後者は、誤差項の分散が一定であり、異なる観測値の誤差項の共分散がゼロであることを表している。

　これらに加えて、誤差項が正規分布に従うと仮定すれば、回帰係数に対する $t$ 検定や複数の係数に対する $F$ 検定を、有限標本において正確に行うことができる。ただし、誤差項の正規性は、最小二乗推定量がBLUEとなるための条件ではない。

　一方、モデルを将来の未知のデータの予測に利用する場合には、標本内での当てはまりだけでモデルを評価することはできない。母集団における変数間の関係を推定し、その回帰係数を解釈することと、未知の観測値を正確に予測することとは、目的が異なるからである。さらに、回帰係数に因果的な意味を与えるためには、通常の回帰分析とは別に、交絡などを考慮した因果識別の条件が必要になる。

　したがって、モデルを評価するときには、標本内における当てはまりと、標本外における予測性能を区別する必要がある。説明変数を追加したりモデルを複雑にしたりすれば、一般に標本内の当てはまりは改善する。しかし、その改善が標本外でも維持されるとは限らず、標本内の当てはまりがよいことは、標本外の予測性能がよいことを保証しない。

Prediction (in sample/out of sample)

plt.rc("figure", figsize=(4,2.5))
plt.rc("font", size=14)

# 人工データの生成
nsample = 50
sig = 0.25
x1 = np.linspace(0, 20, nsample)
X = np.column_stack((x1, np.sin(x1), (x1 - 5) ** 2))
X = sm.add_constant(X)
beta = [5.0, 0.5, 0.5, -0.02]
y_true = np.dot(X, beta)
rng = np.random.default_rng(1007)
y = y_true + sig * rng.normal(loc=0, scale=1, size=nsample)
plt.plot(x1,y,"o")

# 回帰係数の推定
olsmod = sm.OLS(y, X)
olsres = olsmod.fit()
print(olsres.summary())

標本内では、olsres.predict(X) によって得られる値と olsres.fittedvalues は同一である。

# 標本内の当てはめ値
ypred = olsres.predict(X)
print(ypred)
plt.plot(x1, y, "o", label="Data")
plt.plot(x1, ypred, "r", label="OLS prediction")
plt.plot(x1,olsres.fittedvalues,"+", label="fittedvalues")
plt.legend()

　つぎに、推定済みの回帰モデルに新しい説明変数を与え、標本外予測値を計算する。ここで新しく設定する $x_1$ の範囲は $20.5$ から $25$ であり、推定に使用した標本の範囲である $0$ から $20$ の外側にある。したがって、これは標本外予測であると同時に、外挿予測でもある。

　ここでの目的は過剰適合を検証することではなく、標本内の当てはめと標本外予測の計算上の違いを、statsmodelsのコードを通じて確認することである。

#　説明変数Xnewの新しい標本の生成と外挿予測
x1n = np.linspace(20.5, 25, 10)
Xnew = np.column_stack((x1n, np.sin(x1n), (x1n - 5) ** 2))
Xnew = sm.add_constant(Xnew)
ynewpred = olsres.predict(Xnew)  # predict out of sample
print(ynewpred)

[10.94392279 10.81555373 10.58805595 10.30042264 10.0039826   9.74983311
  9.57632954  9.49969465  9.51004665  9.57381863]

#　結果の描画
ynew_true = np.dot(Xnew, beta)

fig, ax = plt.subplots()
ax.plot(x1, y, "o", label="Data")
ax.plot(
    np.hstack((x1, x1n)),
    np.hstack((y_true, ynew_true)),
    "b-",
    label="True"
)
ax.plot(
    np.hstack((x1, x1n)),
    np.hstack((ypred, ynewpred)),
    "r",
    label="OLS prediction"
)
ax.legend(loc="best")

#　formula APIを用いた新しい説明変数に対する予測値
from statsmodels.formula.api import ols

data = {"x1": x1, "y": y}

res = ols("y ~ x1 + np.sin(x1) + I((x1-5)**2)", data=data).fit()

res.params

Intercept           4.997656
x1                  0.497429
np.sin(x1)          0.436318
I((x1 - 5) ** 2)   -0.019505
dtype: float64

res.predict(exog=dict(x1=x1n))

0    10.943923
1    10.815554
2    10.588056
3    10.300423
4    10.003983
5     9.749833
6     9.576330
7     9.499695
8     9.510047
9     9.573819
dtype: float64

formula_pred = res.predict(exog={"x1": x1n}).to_numpy()

np.allclose(ynewpred, formula_pred)

True

用語の使い分け：

fit()：回帰係数の推定
fittedvalues：標本内の当てはめ値
predict(X)：標本内では当てはめ値と同じ
predict(Xnew)：標本外予測値
今回の Xnew：外挿予測

発展：線形回帰の周辺モデル

　statsmodelsでは線形回帰モデル
$y_i=\beta_0+\beta_1 x_{1i} +\cdots+\beta_{ki}x_{ki}+ \varepsilon_i $
または行列表示では
$\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta} +\boldsymbol{\varepsilon}$
の切片($\beta_0$)と回帰係数($\beta_1,\cdots,\beta_k$)を

- 最小二乗法|OLS|
- 加重最小二乗法|WLS|
- 一般化最小二乗法|GLS|
- 再帰的最小二乗法|Recursive LS|

という4つの方法により推定する。

　GLSは、誤差の分散不均一性や自己相関を含む誤差共分散構造を考慮して推定する。WLSはそのうち、観測ごとに誤差分散が異なる場合を扱う。Recursive LSは、データを順次追加しながら回帰係数を推定し、係数の安定性を調べる方法である。
　線形回帰というときには、
説明変数$X$は、実験計画のように所与の値として扱う場合と、観察データのように確率変数として扱う場合がある。確率変数として扱う場合には、一般に外生性条件
$$E[\varepsilon_i\mid X]=0$$
を仮定する。

　一般化線形モデル(GLM)では、被説明変数の分布を指数型分布族から選び、その条件付き平均と説明変数の線形予測子をリンク関数によって結び付ける。一般化線形モデルに関連するモデルとして

• 一般化推定方程式
• 一般化混合モデル
• 一般化加法モデル

などがある。線形回帰ではOLSを用いるが、一般化線形モデルやその関連モデルでは、最尤法、準尤度法、推定方程式などによって回帰係数を推定する。

　また、時間の経過とともに記録したデータを時系列データという。時系列データを扱うときには従属変数$y$と独立変数$x$がともに確率変数であるときがある。 回帰モデルに説明変数の現在値と過去の値を含めるものを分布ラグモデルという。被説明変数の過去の値を説明変数とするものを自己回帰モデルという。被説明変数の過去の値と説明変数の現在値・過去値をともに含めるモデルには、自己回帰分布ラグモデル（ARDL）がある。

　複数の内生変数が相互に関係し、それらが複数の方程式によって同時に決定されるモデルを同時方程式モデルという。複数の時系列変数の相互関係を分析する方法としてVARがあり、変数間に共和分関係がある場合にはVECMが用いられる。

　時系列分析では、原系列が定常過程であるか、単位根をもつ非定常過程であるかを調べる必要がある。ランダムウォークは代表的な1次和分過程$I(1)$である。1次和分過程がすべて厳密なランダムウォークであるとは限らない。複数の$I(1)$過程について、ある非自明な線形結合が$I(0)$となるとき、それらの変数は共和分関係にあるという。共和分関係にある変数は共通する確率的トレンドをもち、共和分ベクトルによる線形結合ではそのトレンドが相殺される。VECMは、この長期的な共和分関係と短期的な変化を同時に表現するモデルだ。

108228091201082026

Python3ではじめるシステムトレード【第2版】環境構築と売買戦略