重回帰分析で住宅価格の予測モデルを構築し、決定係数 $R^2 = 0.85$ を得たとしましょう。一見すると良いモデルに見えますが、もし残差に系統的なパターンが潜んでいたらどうでしょうか。例えば、高価格帯の物件ほど予測の誤差が大きくなっていたり、予測値に対して残差がU字型のカーブを描いていたりしたら、$R^2$ がいくら高くても、そのモデルの係数推定値やp値を信頼することはできません。
回帰分析は、誤差項に関するいくつかの前提条件（線形性、等分散性、正規性、独立性）が満たされているときにのみ妥当な結果を出します。これらの前提が破られると、推定量の最良線形不偏推定量（BLUE）としての性質が失われ、検定のp値や信頼区間が不正確になります。
残差分析（Residual Analysis）は、残差（実測値と予測値の差）のパターンを体系的に調べることで、これらの前提条件が満たされているかを検証する強力な診断手法です。残差は誤差項の推定値であり、残差に系統的なパターンが見られれば、モデルの前提に問題があることを示唆します。
残差分析を理解すると、以下のような応用が開けます。
- モデル改善: 非線形性や異分散の検出を通じて、変数変換やモデル構造の改良につなげる
- 外れ値・影響点の同定: モデルの推定結果に不当に大きな影響を与えているデータ点を発見する
- 報告の信頼性: 学術論文や報告書で回帰分析の結果を報告する際に、前提条件の検証結果を添えて分析の妥当性を保証する
- 予測精度の向上: 等分散性の違反を検出し、加重最小二乗法や一般化最小二乗法への切り替えを判断する
本記事では、残差の各種類（生残差、標準化残差、スチューデント化残差、削除残差）の定義と性質を数学的に導出し、残差プロットの読み方を体系的に解説します。さらに、てこ比やクックの距離による影響度分析の理論を示し、Pythonで包括的な回帰診断を実装します。
本記事の内容
- 回帰分析の4つの前提条件とその違反の影響
- 残差の4つの種類（生残差・標準化残差・スチューデント化残差・削除残差）の定義と導出
- ハット行列とてこ比の理論
- 残差プロット（6種類）の作成と読み方
- クックの距離と影響度分析の理論
- Pythonでの包括的な回帰診断の実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
回帰分析の4つの前提条件
正規線形モデルの仮定
重回帰分析の正規線形モデルは次のように定式化されます。
$$ \begin{equation} \bm{y} = \bm{X}\bm{\beta} + \bm{\varepsilon}, \quad \bm{\varepsilon} \sim N(\bm{0}, \sigma^2\bm{I}) \end{equation} $$
この式には暗黙のうちに4つの前提条件が含まれています。
前提1: 線形性（Linearity）
$$ E[\bm{y}|\bm{X}] = \bm{X}\bm{\beta} $$
目的変数の条件付き期待値が説明変数の線形関数であることを仮定しています。真の関係が $E[y|x] = \beta_0 + \beta_1 x + \beta_2 x^2$ のように非線形であれば、線形モデルでは系統的なバイアスが生じます。
前提2: 等分散性（Homoscedasticity）
$$ \text{Var}(\varepsilon_i) = \sigma^2 \quad \text{for all } i $$
誤差の分散が全てのデータ点で一定であることを仮定しています。例えば収入を予測するモデルで、高収入者ほど予測の誤差が大きくなる場合（分散が収入に比例する場合）、等分散性は成り立ちません。この状態を異分散性（heteroscedasticity）と呼びます。
前提3: 正規性（Normality）
$$ \varepsilon_i \sim N(0, \sigma^2) $$
誤差が正規分布に従うことを仮定しています。この仮定は、回帰係数のt検定やF検定の正確性に必要です。ただし、大標本では中心極限定理により正規性の仮定の重要性は低下します。
前提4: 独立性（Independence）
$$ \text{Cov}(\varepsilon_i, \varepsilon_j) = 0 \quad \text{for } i \neq j $$
誤差間が無相関であることを仮定しています。時系列データでは、誤差に自己相関があることが多く、この仮定が破られやすいです。
前提の違反がもたらす影響
各前提が破られた場合の影響を整理します。
| 前提の違反 | 推定量の不偏性 | BLUE性 | t検定・F検定 | 信頼区間 |
|---|---|---|---|---|
| 線形性の違反 | 失われる | 失われる | 無効 | 無効 |
| 等分散性の違反 | 保たれる | 失われる | 不正確 | 不正確 |
| 正規性の違反 | 保たれる | 保たれる（GLTにより） | 小標本で不正確 | 小標本で不正確 |
| 独立性の違反 | 保たれる | 失われる | 不正確（しばしば過度に有意） | 過度に狭い |
線形性の違反は最も深刻で、推定値自体にバイアスが生じます。等分散性と独立性の違反は推定値は不偏のままですが、標準誤差の推定が不正確になります。正規性の違反は大標本では問題が緩和されますが、小標本では検定結果に影響します。
これらの前提を検証するために、残差を詳しく分析します。まず、残差の数学的な性質を理解しましょう。
ハット行列と残差の性質
ハット行列の定義
最小二乗推定量 $\hat{\bm{\beta}} = (\bm{X}^T\bm{X})^{-1}\bm{X}^T\bm{y}$ による予測値は:
$$ \begin{equation} \hat{\bm{y}} = \bm{X}\hat{\bm{\beta}} = \bm{X}(\bm{X}^T\bm{X})^{-1}\bm{X}^T\bm{y} = \bm{H}\bm{y} \end{equation} $$
ここで定義された行列:
$$ \begin{equation} \bm{H} = \bm{X}(\bm{X}^T\bm{X})^{-1}\bm{X}^T \end{equation} $$
をハット行列（hat matrix）と呼びます。$\bm{y}$ に $\bm{H}$ を掛けると $\hat{\bm{y}}$（”hat” の付いた $\bm{y}$）になることからこの名前が付いています。
ハット行列の性質
ハット行列は以下の重要な性質を持ちます。
対称性: $\bm{H}^T = \bm{H}$
冪等性: $\bm{H}^2 = \bm{H}$（射影行列の性質）
トレース: $\text{tr}(\bm{H}) = p + 1$（切片を含む場合。$p$ は説明変数の数）
対角成分の範囲: $0 \leq h_{ii} \leq 1$
冪等性は重要な性質です。$\bm{H}$ が冪等行列であることは、$\hat{\bm{y}}$ が $\bm{X}$ の列空間への直交射影であることを意味しています。
てこ比（Leverage）
ハット行列の対角成分 $h_{ii}$ をてこ比（leverage）と呼びます。
$$ \begin{equation} h_{ii} = \bm{x}_i^T(\bm{X}^T\bm{X})^{-1}\bm{x}_i \end{equation} $$
てこ比の解釈は以下の通りです。
$h_{ii}$ は、$i$ 番目のデータ点の説明変数の値 $\bm{x}_i$ が、説明変数空間の「中心」からどれだけ離れているかを測る指標です。$\hat{y}_i = \sum_{j} h_{ij}y_j$ という式を見ると、$h_{ii}$ は自分自身の $y_i$ が自分の予測値 $\hat{y}_i$ にどれだけ影響するかの重みになっています。$h_{ii}$ が大きいデータ点は、回帰直線を自分の方に「引き寄せる」力が強く、結果として生残差が小さくなりがちです。
てこ比の平均値は $\bar{h} = (p+1)/n$ です。$h_{ii} > 2(p+1)/n$ のデータ点は高てこ比点と呼ばれ、回帰結果に大きな影響を与える潜在的な影響点です。
残差の共分散構造
生残差ベクトル $\hat{\bm{e}} = \bm{y} – \hat{\bm{y}} = (\bm{I} – \bm{H})\bm{y}$ の共分散行列を導出しましょう。
$\bm{y} = \bm{X}\bm{\beta} + \bm{\varepsilon}$ より:
$$ \hat{\bm{e}} = (\bm{I} – \bm{H})\bm{y} = (\bm{I} – \bm{H})(\bm{X}\bm{\beta} + \bm{\varepsilon}) = (\bm{I} – \bm{H})\bm{\varepsilon} $$
ここで $(\bm{I} – \bm{H})\bm{X} = \bm{X} – \bm{H}\bm{X} = \bm{X} – \bm{X} = \bm{0}$ を使いました。したがって:
$$ \begin{equation} \text{Cov}(\hat{\bm{e}}) = (\bm{I} – \bm{H})\text{Cov}(\bm{\varepsilon})(\bm{I} – \bm{H})^T = \sigma^2(\bm{I} – \bm{H}) \end{equation} $$
最後の等号では $(\bm{I} – \bm{H})^T = \bm{I} – \bm{H}$（対称性）と $(\bm{I} – \bm{H})^2 = \bm{I} – \bm{H}$（冪等性）を使いました。
この結果から重要なことがわかります。
第1に、$i$ 番目の残差の分散は $\text{Var}(\hat{e}_i) = \sigma^2(1 – h_{ii})$ です。つまり、てこ比が大きいデータ点の残差は、真の誤差よりも分散が小さくなります。高てこ比の点はモデルに大きな影響を与え、残差を見かけ上小さくしてしまうのです。
第2に、異なるデータ点の残差間には相関があります。$\text{Cov}(\hat{e}_i, \hat{e}_j) = -\sigma^2 h_{ij}$ です。真の誤差 $\varepsilon_i$ が互いに独立であっても、残差 $\hat{e}_i$ は互いに相関しています。
これらの性質を踏まえると、生残差をそのまま比較するのは適切でないことがわかります。次のセクションで、適切に修正された残差の定義を見ていきましょう。
残差の種類
生残差（Raw Residuals）
最も基本的な残差です。
$$ \begin{equation} \hat{e}_i = y_i – \hat{y}_i \end{equation} $$
前節で見たように、$\text{Var}(\hat{e}_i) = \sigma^2(1 – h_{ii})$ であり、てこ比によって分散が異なるため、生残差をそのまま比較するのは不公平です。
標準化残差（Standardized Residuals）
生残差をその標準偏差で割って標準化したものです。
$$ \begin{equation} r_i = \frac{\hat{e}_i}{\hat{\sigma}\sqrt{1 – h_{ii}}} \end{equation} $$
ここで $\hat{\sigma} = \sqrt{\text{RSS}/(n – p – 1)}$ は $\sigma$ の推定量です。前提条件が正しければ、標準化残差はおよそ $N(0, 1)$ に従います。$|r_i| > 2$ のデータ点は全体の約5%にあたり、$|r_i| > 3$ のデータ点は外れ値の候補です。
ただし、$\hat{\sigma}$ の計算に全データ（$i$ 番目を含む）を使っているため、外れ値の存在が $\hat{\sigma}$ を膨らませ、その結果として外れ値の標準化残差が「見えにくく」なる問題があります（マスキング効果）。
スチューデント化残差（Studentized Residuals / 外部スチューデント化残差）
マスキング効果を避けるため、$i$ 番目のデータ点を除外して $\sigma$ を推定する方法が考えられます。$i$ 番目を除いた残差標準偏差を $\hat{\sigma}_{(i)}$ と書きます。
$$ \begin{equation} t_i = \frac{\hat{e}_i}{\hat{\sigma}_{(i)}\sqrt{1 – h_{ii}}} \end{equation} $$
これを外部スチューデント化残差（externally studentized residual）または削除スチューデント化残差と呼びます。前提条件のもとで $t_i$ は自由度 $n – p – 2$ の $t$ 分布に従い、外れ値検定に直接使えます。
$\hat{\sigma}_{(i)}$ を直接計算するには全データ点について回帰を繰り返す必要がありますが、次の公式を使えば1回の回帰から計算できます。
$$ \begin{equation} \hat{\sigma}_{(i)}^2 = \frac{(n-p-1)\hat{\sigma}^2 – \hat{e}_i^2/(1-h_{ii})}{n-p-2} \end{equation} $$
この式は、$i$ 番目のデータ点を除いたRSSが元のRSSから更新式で計算できることを利用しています。
外れ値が1点だけであれば、スチューデント化残差で効果的に検出できます。しかし、複数の外れ値が存在する場合は、互いにマスキングし合う可能性があり、より頑健な手法（例えばロバスト回帰）が必要になることがあります。
各種残差の性質を理解した上で、具体的な残差プロットの作成と読み方を見ていきましょう。
残差プロットによる前提の検証
残差 vs 予測値プロット — 線形性と等分散性の検証
最も基本的な診断プロットは、横軸に予測値 $\hat{y}_i$、縦軸に残差 $\hat{e}_i$（または標準化残差 $r_i$）をプロットしたものです。
理想的なパターン: 残差がゼロの水平線の周りにランダムに散らばっている。
非線形性のシグナル: 残差がU字型やS字型のパターンを示す。これは、真の関係が非線形であるのに線形モデルを当てはめていることを意味します。対処法は、説明変数の二次項や交互作用項を追加する、あるいは変数変換（対数変換など）を行うことです。
異分散性のシグナル: 残差の散らばりが予測値に対して系統的に変化する。典型的には「漏斗形」（ファンネル形状）で、予測値が大きいほど残差の散らばりも大きくなります。対処法は、加重最小二乗法（WLS）を用いるか、目的変数に対数変換を施すことです。
Scale-Locationプロット — 等分散性の検証
横軸に予測値、縦軸に $\sqrt{|r_i|}$（標準化残差の絶対値の平方根）をプロットします。等分散性が成り立っていれば、$\sqrt{|r_i|}$ は予測値によらず一定の水準に分布します。右上がりの傾向があれば、正の異分散性（予測値が大きいほど分散が大きい）を示唆します。
Q-Qプロット — 正規性の検証
Q-Qプロット（Quantile-Quantile plot）は、標準化残差の順序統計量を標準正規分布の理論的分位点に対してプロットしたものです。
残差を昇順に並べて $r_{(1)} \leq r_{(2)} \leq \dots \leq r_{(n)}$ とし、対応する理論的分位点を $q_i = \Phi^{-1}(i/(n+1))$ とします。$r_{(i)}$ と $q_i$ をプロットしたとき:
理想的なパターン: 点が45度の直線上に並ぶ。
裾が重い分布（尖度が大きい）: 両端が直線から上下に逸脱する（S字型）。
右に歪んだ分布: 右端で直線から上に逸脱。
左に歪んだ分布: 左端で直線から下に逸脱。
Q-Qプロットは視覚的な判断であるため、補助的にShapiro-Wilk検定やAnderson-Darling検定などの統計的検定を行うことも推奨されます。
残差 vs 説明変数プロット
残差を各説明変数に対してプロットすることで、特定の変数との非線形関係を検出できます。あるいは、偏残差プロット（partial residual plot, component-plus-residual plot）を用いると、他の変数の影響を除いた上での各変数と目的変数の関係を可視化できます。
偏残差は次のように定義されます。
$$ \begin{equation} e_j^{(\text{partial})} = \hat{e}_i + \hat{\beta}_j x_{ij} \end{equation} $$
横軸に $x_j$、縦軸に偏残差をプロットしたグラフで曲線的なパターンが見られれば、$x_j$ との関係に非線形成分があることを示唆します。
残差プロットでモデルの前提を検証した後は、個別のデータ点がモデルにどれだけ影響を与えているかを評価することも重要です。次に影響度分析の理論を見ていきましょう。
影響度分析
影響点と外れ値の区別
外れ値（outlier）: 残差が大きいデータ点。$|r_i| > 2$ または $|t_i| > 2$ が目安。
高てこ比点（high leverage point）: 説明変数空間で中心から離れた点。$h_{ii} > 2(p+1)/n$ が目安。
影響点（influential point）: 回帰結果に大きな影響を与える点。外れ値かつ高てこ比であることが多い。
重要なのは、外れ値であることと影響点であることは必ずしも一致しないということです。てこ比が低い外れ値は回帰直線をあまり動かしませんが、てこ比が高い外れ値は回帰直線を大きく歪める可能性があります。逆に、てこ比が高くても残差が小さければ、その点は「外れ値」ではなくモデルの推定を安定化させる「良い影響点」です。
クックの距離の導出
クックの距離（Cook’s distance）は、$i$ 番目のデータ点を除外したときの回帰係数の変化を総合的に測る指標です。
$i$ 番目の点を除いた最小二乗推定量を $\hat{\bm{\beta}}_{(i)}$ とすると、クックの距離は次のように定義されます。
$$ \begin{equation} D_i = \frac{(\hat{\bm{\beta}}_{(i)} – \hat{\bm{\beta}})^T(\bm{X}^T\bm{X})(\hat{\bm{\beta}}_{(i)} – \hat{\bm{\beta}})}{(p+1)\hat{\sigma}^2} \end{equation} $$
この定義は直感的です。$\hat{\bm{\beta}}_{(i)} – \hat{\bm{\beta}}$ は $i$ 番目を除くことによる係数の変化であり、$(\bm{X}^T\bm{X})$ で重み付けすることで、$\hat{\bm{y}}$ の変化量を測っています。分母は自由度と分散で標準化しています。
$i$ 番目のデータ点を実際に除いて回帰をやり直す必要はありません。更新公式を用いると、クックの距離は1回の回帰の結果から計算できます。
$$ \begin{equation} D_i = \frac{r_i^2}{p+1}\cdot\frac{h_{ii}}{1-h_{ii}} \end{equation} $$
この式の構造は非常に示唆的です。クックの距離は:
- 標準化残差の二乗 $r_i^2$: そのデータ点が「目的変数の方向でどれだけ外れているか」
- てこ比に関する量 $h_{ii}/(1-h_{ii})$: そのデータ点が「説明変数の方向でどれだけ特異か」
この2つの積（を定数で割ったもの）として表されます。つまり、クックの距離は外れ値の程度とてこ比の両方を組み合わせた総合的な影響度指標です。
クックの距離の解釈
クックの距離の閾値については、いくつかの基準が提案されています。
- $D_i > 0.5$: 注意が必要
- $D_i > 1$: 強い影響力を持つ点
- $D_i > 4/n$: 比較的保守的な基準
- $D_i > F_{0.5}(p+1, n-p-1)$: F分布の50パーセント点（近似的に0.7〜1.0程度になることが多い）
DFFITS
DFFITSは、$i$ 番目のデータ点を除いたときの予測値の変化を測る指標です。
$$ \begin{equation} \text{DFFITS}_i = \frac{\hat{y}_i – \hat{y}_{(i)}}{\hat{\sigma}_{(i)}\sqrt{h_{ii}}} = t_i\sqrt{\frac{h_{ii}}{1-h_{ii}}} \end{equation} $$
DFFITSの閾値は $2\sqrt{(p+1)/n}$ とされています。
クックの距離とDFFITSは関連していますが、クックの距離が全ての予測値に対する影響を見るのに対し、DFFITSは $i$ 番目の予測値に対する影響のみを見る点が異なります。
これらの理論的指標をPythonで実装し、実際のデータで診断を行いましょう。
Pythonによる実装と可視化
問題のあるデータの作成と回帰診断
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n = 150
# 説明変数
x1 = np.random.uniform(0, 10, n)
x2 = np.random.randn(n)
# 非線形性 + 異分散性を含むデータ
y = 2 + 3*x1 - 0.25*x1**2 + x2 + (0.5 + 0.2*x1)*np.random.randn(n)
# 外れ値を追加
y[0] = 35 # 高い外れ値
y[1] = -15 # 低い外れ値
# 高てこ比点を追加
x1[2] = 20 # x1の範囲外
y[2] = 2 + 3*20 - 0.25*400 + 0.5*np.random.randn(1) # 適度な外れ
X = np.column_stack([np.ones(n), x1, x2])
p_vars = 2 # 説明変数の数（切片を除く）
# --- 最小二乗回帰 ---
beta_hat = np.linalg.lstsq(X, y, rcond=None)[0]
y_hat = X @ beta_hat
e = y - y_hat # 生残差
# --- ハット行列と診断統計量 ---
H = X @ np.linalg.inv(X.T @ X) @ X.T
h = np.diag(H) # てこ比
RSS = np.sum(e**2)
df = n - p_vars - 1
sigma_hat = np.sqrt(RSS / df)
# 標準化残差
r_std = e / (sigma_hat * np.sqrt(1 - h))
# スチューデント化残差（外部）
sigma_i_sq = ((df * sigma_hat**2 - e**2 / (1 - h)) / (df - 1))
sigma_i_sq = np.maximum(sigma_i_sq, 1e-10) # 数値安定性
t_stud = e / (np.sqrt(sigma_i_sq) * np.sqrt(1 - h))
# クックの距離
D_cook = r_std**2 / (p_vars + 1) * h / (1 - h)
# DFFITS
dffits = t_stud * np.sqrt(h / (1 - h))
print("=== 回帰結果 ===")
print(f"beta_hat = {beta_hat}")
print(f"R^2 = {1 - RSS/np.sum((y - y.mean())**2):.4f}")
print(f"sigma_hat = {sigma_hat:.4f}")
print(f"\n=== 上位5点の影響度 ===")
top5 = np.argsort(D_cook)[-5:][::-1]
for idx in top5:
print(f" i={idx:3d}: r={r_std[idx]:7.3f}, "
f"h={h[idx]:.4f}, D={D_cook[idx]:.4f}, "
f"DFFITS={dffits[idx]:.4f}")
この出力から、データ点0, 1は高い標準化残差（外れ値）を持ち、データ点2は高いてこ比（説明変数空間での特異性）を持つことが確認できます。クックの距離が最も大きい点が、モデルに最も強い影響を与えている点です。
6種類の診断プロット
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n = 150
x1 = np.random.uniform(0, 10, n)
x2 = np.random.randn(n)
y = 2 + 3*x1 - 0.25*x1**2 + x2 + (0.5 + 0.2*x1)*np.random.randn(n)
y[0], y[1] = 35, -15
x1[2] = 20
y[2] = 2 + 60 - 100 + 0.5*np.random.randn(1)
X = np.column_stack([np.ones(n), x1, x2])
p_vars = 2
beta_hat = np.linalg.lstsq(X, y, rcond=None)[0]
y_hat = X @ beta_hat
e = y - y_hat
H = X @ np.linalg.inv(X.T @ X) @ X.T
h = np.diag(H)
sigma_hat = np.sqrt(np.sum(e**2) / (n - p_vars - 1))
r_std = e / (sigma_hat * np.sqrt(1 - h))
D_cook = r_std**2 / (p_vars + 1) * h / (1 - h)
fig, axes = plt.subplots(2, 3, figsize=(17, 10))
# (a) 残差 vs 予測値
ax = axes[0, 0]
ax.scatter(y_hat, r_std, alpha=0.5, s=20, c='steelblue')
ax.axhline(0, color='red', linewidth=1, linestyle='--')
ax.axhline(2, color='gray', linewidth=0.5, linestyle=':')
ax.axhline(-2, color='gray', linewidth=0.5, linestyle=':')
# 外れ値にラベル
for idx in np.where(np.abs(r_std) > 2.5)[0]:
ax.annotate(str(idx), (y_hat[idx], r_std[idx]), fontsize=8, color='red')
# LOESS近似
from numpy.polynomial import polynomial as P
sort_idx = np.argsort(y_hat)
coef = P.polyfit(y_hat[sort_idx], r_std[sort_idx], 3)
y_fit_line = P.polyval(y_hat[sort_idx], coef)
ax.plot(y_hat[sort_idx], y_fit_line, 'r-', linewidth=1.5, alpha=0.7)
ax.set_xlabel('Fitted values $\\hat{y}$', fontsize=11)
ax.set_ylabel('Standardized residuals', fontsize=11)
ax.set_title('(a) Residuals vs Fitted', fontsize=12)
ax.grid(True, alpha=0.3)
# (b) Q-Qプロット
ax = axes[0, 1]
res_sorted = np.sort(r_std)
theoretical = stats.norm.ppf(np.arange(1, n+1) / (n+1))
ax.scatter(theoretical, res_sorted, alpha=0.5, s=20, c='steelblue')
ax.plot([-4, 4], [-4, 4], 'r--', linewidth=1.5)
for idx in np.where(np.abs(r_std) > 2.5)[0]:
rank = np.searchsorted(np.sort(r_std), r_std[idx])
if rank < n:
ax.annotate(str(idx), (theoretical[rank], res_sorted[rank]),
fontsize=8, color='red')
ax.set_xlabel('Theoretical quantiles', fontsize=11)
ax.set_ylabel('Standardized residuals', fontsize=11)
ax.set_title('(b) Normal Q-Q Plot', fontsize=12)
ax.grid(True, alpha=0.3)
# (c) Scale-Locationプロット
ax = axes[0, 2]
ax.scatter(y_hat, np.sqrt(np.abs(r_std)), alpha=0.5, s=20, c='steelblue')
sort_idx = np.argsort(y_hat)
coef_sl = P.polyfit(y_hat[sort_idx], np.sqrt(np.abs(r_std[sort_idx])), 3)
y_sl_line = P.polyval(y_hat[sort_idx], coef_sl)
ax.plot(y_hat[sort_idx], y_sl_line, 'r-', linewidth=1.5, alpha=0.7)
ax.set_xlabel('Fitted values $\\hat{y}$', fontsize=11)
ax.set_ylabel('$\\sqrt{|\\mathrm{Standardized\\ residual}|}$', fontsize=11)
ax.set_title('(c) Scale-Location', fontsize=12)
ax.grid(True, alpha=0.3)
# (d) てこ比 vs 標準化残差
ax = axes[1, 0]
sc = ax.scatter(h, r_std, c=D_cook, cmap='Reds', alpha=0.6, s=25,
edgecolors='gray', linewidths=0.5)
plt.colorbar(sc, ax=ax, label="Cook's distance")
ax.axhline(2, color='gray', linewidth=0.5, linestyle=':')
ax.axhline(-2, color='gray', linewidth=0.5, linestyle=':')
ax.axvline(2*(p_vars+1)/n, color='orange', linewidth=1, linestyle='--',
label=f'$2(p+1)/n = {2*(p_vars+1)/n:.3f}$')
for idx in np.argsort(D_cook)[-3:][::-1]:
ax.annotate(str(idx), (h[idx], r_std[idx]), fontsize=9, color='red',
fontweight='bold')
ax.set_xlabel('Leverage $h_{ii}$', fontsize=11)
ax.set_ylabel('Standardized residuals', fontsize=11)
ax.set_title('(d) Residuals vs Leverage', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (e) クックの距離
ax = axes[1, 1]
markerline, stemlines, baseline = ax.stem(range(n), D_cook, linefmt='steelblue',
markerfmt='o', basefmt='k-')
markerline.set_markersize(3)
ax.axhline(0.5, color='orange', linewidth=1.5, linestyle='--', label='$D=0.5$')
ax.axhline(1.0, color='red', linewidth=1.5, linestyle='--', label='$D=1.0$')
ax.axhline(4/n, color='gray', linewidth=1, linestyle=':', label=f'$4/n={4/n:.3f}$')
top_idx = np.argsort(D_cook)[-3:]
for idx in top_idx:
ax.annotate(str(idx), (idx, D_cook[idx]+0.02), fontsize=9,
ha='center', color='red')
ax.set_xlabel('Observation index', fontsize=11)
ax.set_ylabel("Cook's distance", fontsize=11)
ax.set_title("(e) Cook's Distance", fontsize=12)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
# (f) 残差のヒストグラムとKDE
ax = axes[1, 2]
ax.hist(r_std, bins=25, density=True, alpha=0.5, color='steelblue',
edgecolor='black', label='Residuals')
x_range = np.linspace(-5, 5, 200)
ax.plot(x_range, stats.norm.pdf(x_range), 'r-', linewidth=2, label='N(0,1)')
# KDE
from scipy.stats import gaussian_kde
kde = gaussian_kde(r_std)
ax.plot(x_range, kde(x_range), 'g--', linewidth=1.5, label='KDE')
ax.set_xlabel('Standardized residual', fontsize=11)
ax.set_ylabel('Density', fontsize=11)
ax.set_title('(f) Residual Distribution', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# Shapiro-Wilk検定
stat, p_value = stats.shapiro(r_std)
ax.text(0.05, 0.95, f'Shapiro-Wilk: W={stat:.4f}\np={p_value:.4f}',
transform=ax.transAxes, fontsize=9, va='top',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
plt.tight_layout()
plt.show()
6つの診断プロットから、このモデルの問題点が読み取れます。
-
残差 vs 予測値（図a）: 赤い曲線（3次多項式フィット）が水平線から系統的に逸脱しており、U字型のパターンが見えます。これは真の関係が $x_1^2$ の項を含む非線形であるのに、線形モデルを当てはめていることに起因する線形性の違反を示しています。データ点0と1は標準化残差が $\pm 2$ の点線を大きく超えており、外れ値として検出されています
-
Q-Qプロット（図b）: 大部分の点は45度線上にありますが、両端で直線から逸脱する点があります。特にデータ点0（右上の逸脱）と1（左下の逸脱）は外れ値の影響です。全体的な裾の重さは、異分散性や外れ値の存在を反映しています
-
Scale-Locationプロット（図c）: 赤い曲線が右上がりの傾向を示しており、予測値が大きくなるほど残差の散らばりが大きくなる異分散性を示唆しています。これはデータ生成時に $\sigma(x_1) = 0.5 + 0.2x_1$ として異分散性を組み込んだことと一致します
-
てこ比 vs 残差（図d）: 右上と右下の領域にあるデータ点は、てこ比と残差の両方が大きく、モデルへの影響が強い点です。クックの距離が大きい点ほど赤く表示されており、特にデータ点2は高いてこ比を持つことが確認できます
-
クックの距離（図e）: データ点0のクックの距離が突出して大きく、この点がモデル全体に不当に大きな影響を与えていることがわかります。$4/n$ の基準線（灰色点線）を超える点が複数あり、これらの点の妥当性を個別に検討する必要があります
-
残差分布（図f）: ヒストグラムとKDE（カーネル密度推定）を標準正規分布と比較しています。おおむね正規分布に近い形状ですが、外れ値により裾がやや重くなっています。Shapiro-Wilk検定の結果（p値）が小さければ、正規性の仮定は統計的に棄却されます
前提条件を満たすモデルとの比較
問題の原因を修正した改善モデルで、残差パターンがどう変わるかを確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n = 150
x1 = np.random.uniform(0, 10, n)
x2 = np.random.randn(n)
y = 2 + 3*x1 - 0.25*x1**2 + x2 + (0.5 + 0.2*x1)*np.random.randn(n)
# 外れ値なし、x1^2項を含む改善モデル
X_improved = np.column_stack([np.ones(n), x1, x1**2, x2])
beta_imp = np.linalg.lstsq(X_improved, y, rcond=None)[0]
y_hat_imp = X_improved @ beta_imp
e_imp = y - y_hat_imp
H_imp = X_improved @ np.linalg.inv(X_improved.T @ X_improved) @ X_improved.T
h_imp = np.diag(H_imp)
sigma_imp = np.sqrt(np.sum(e_imp**2) / (n - 4))
r_imp = e_imp / (sigma_imp * np.sqrt(1 - h_imp))
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# 残差 vs 予測値
ax = axes[0]
ax.scatter(y_hat_imp, r_imp, alpha=0.5, s=20, c='steelblue')
ax.axhline(0, color='red', linewidth=1, linestyle='--')
ax.axhline(2, color='gray', linewidth=0.5, linestyle=':')
ax.axhline(-2, color='gray', linewidth=0.5, linestyle=':')
ax.set_xlabel('Fitted values', fontsize=11)
ax.set_ylabel('Standardized residuals', fontsize=11)
ax.set_title('Residuals vs Fitted (improved)', fontsize=12)
ax.grid(True, alpha=0.3)
# Q-Qプロット
ax = axes[1]
res_sorted = np.sort(r_imp)
theoretical = stats.norm.ppf(np.arange(1, n+1) / (n+1))
ax.scatter(theoretical, res_sorted, alpha=0.5, s=20, c='steelblue')
ax.plot([-3, 3], [-3, 3], 'r--', linewidth=1.5)
ax.set_xlabel('Theoretical quantiles', fontsize=11)
ax.set_ylabel('Standardized residuals', fontsize=11)
ax.set_title('Q-Q Plot (improved)', fontsize=12)
ax.grid(True, alpha=0.3)
# Scale-Location
ax = axes[2]
ax.scatter(y_hat_imp, np.sqrt(np.abs(r_imp)), alpha=0.5, s=20, c='steelblue')
ax.set_xlabel('Fitted values', fontsize=11)
ax.set_ylabel('$\\sqrt{|r_i|}$', fontsize=11)
ax.set_title('Scale-Location (improved)', fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
R2_imp = 1 - np.sum(e_imp**2) / np.sum((y - y.mean())**2)
print(f"改善モデルの R^2 = {R2_imp:.4f}")
改善モデル（$x_1^2$ 項を追加し、外れ値を除外）の診断プロットでは、前のモデルで見られた問題が大幅に解消されています。残差 vs 予測値プロットの系統的なU字パターンが消え、残差がランダムに散らばっています。Q-Qプロットの点も直線に近づき、正規性が改善されています。ただし、Scale-Locationプロットにはまだ若干の傾向が残っている可能性があり、異分散性の完全な解消には加重最小二乗法が必要かもしれません。
まとめ
本記事では、残差分析とモデル診断の理論と実践を体系的に解説しました。
- 前提条件: 回帰分析は線形性・等分散性・正規性・独立性の4つの前提に基づく。これらの違反は推定量の性質やp値の信頼性に影響する
- ハット行列: $\bm{H} = \bm{X}(\bm{X}^T\bm{X})^{-1}\bm{X}^T$ の対角成分（てこ比 $h_{ii}$）がデータ点の影響力を決める
- 残差の種類: 生残差、標準化残差 $r_i = \hat{e}_i/(\hat{\sigma}\sqrt{1-h_{ii}})$、スチューデント化残差（外部スチューデント化残差）の3種があり、用途に応じて使い分ける
- 診断プロット: 残差 vs 予測値（線形性・等分散性）、Q-Qプロット（正規性）、Scale-Location（等分散性）、てこ比 vs 残差（影響点）の4種を基本とする
- 影響度指標: クックの距離 $D_i = r_i^2/(p+1)\cdot h_{ii}/(1-h_{ii})$ は外れ値度とてこ比を統合した総合指標。$D_i > 1$ は強い影響点
残差分析は回帰分析を行うたびに必ず実施すべき手続きです。$R^2$ やp値だけを報告するのではなく、残差プロットを添えてモデルの妥当性を確認することが、信頼できる統計分析の基本です。
次のステップとして、以下の記事も参考にしてください。