22. GLMの影響点診断：Cook距離・DFBETASと外れ値検出

影響点診断の動機：なぜGLMで特別な扱いが必要か

OLS（最小二乗法）における影響点診断は、ハット行列 $H = X(X’X)^{-1}X’$ の対角要素と標準化残差を組み合わせた手法として長年用いられてきました。しかし一般化線形モデルでは、応答変数の分布が正規分布以外（二項分布・ポアソン分布など）であり、最尤推定にIRLS（反復重み付き最小二乗法）が用いられるため、OLSの診断量をそのまま適用することはできません。

IRLSの各反復ステップでは、観測値ごとに異なる重み $w_i$ が付与されます。この重みはリンク関数の微分と分散関数から定まり、観測値の位置によって値が変化します。したがってOLSの等重みを前提とした診断量では、一般化線形モデルにおける各観測の実際の影響を正確に反映できません。

影響点診断では、以下の三類型を区別することが基本となります。外れ値とは応答変数の値がモデルの予測から大きく乖離している観測を指します。高レバレッジ点とは説明変数空間において他の観測から離れた位置にある観測を指し、予測値を強く引き寄せる潜在的な力を持ちます。影響点とはその観測を取り除いたときに推定係数が大きく変化する観測を指します。

限界：三類型は必ずしも排他的ではありません。ある観測が高レバレッジ点であると同時に影響点でもあることは頻繁に起こり得ます。逆に高レバレッジ点であっても回帰直線上に乗っている場合は影響点にならないため、単一の診断量のみで分類を完結させることは適切ではありません。

GLMのハット行列とレバレッジ

IRLSの収束後、一般化線形モデルにおける加重ハット行列は次のように定義されます。

$$
H_w = W^{1/2} X (X’ W X)^{-1} X’ W^{1/2}
$$

ここで $W = \mathrm{diag}(w_1, \ldots, w_n)$ は観測ごとの作業重みを対角に並べた行列であり、$X$ は計画行列です。OLSのハット行列は $W = I$ の特殊ケースとして包含されます。

$i$ 番目の観測のレバレッジ $h_{ii}$ は $H_w$ の $i$ 番目の対角要素として定義されます。$h_{ii}$ は作業応答変数の予測値に対する $i$ 番目の観測値の影響度を測り、$0 \le h_{ii} \le 1$ を満たします。高レバレッジの目安として $h_{ii} > 2p/n$（$p$：パラメータ数、$n$：サンプルサイズ）が広く用いられます。

一般化線形モデルにおける作業重み $w_i$ は、リンク関数の微分 $\left(d\mu/d\eta\right)^2$ と分散関数 $V(\mu_i)$ から

$$
w_i = \frac{\left(d\mu_i / d\eta_i\right)^2}{V(\mu_i)}
$$

として導かれます。したがって予測値 $\hat{\mu}_i$ が分布の裾（例：二項分布では確率が0または1に近い領域）に位置する観測ほど $w_i$ が小さくなり、レバレッジの値も変化します。これはOLSにはない挙動であり、診断値の解釈にモデル固有の文脈が必要な理由です。

仮定：加重ハット行列に基づくレバレッジはモデルが正しく特定されていることを前提とします。リンク関数や分布族の選択が誤っている場合、作業重みが不適切な値を取り、レバレッジの解釈が歪みます。

限界：レバレッジは説明変数空間における孤立度を測るに過ぎず、その観測が実際に推定係数を大きく動かしているかどうかはレバレッジ単独では判断できません。影響の程度を評価するにはCook距離等の追加診断量が必要となります。

GLMのCook距離：推定量への影響の総合指標

Cook距離は第 $i$ 観測を除外したときの係数ベクトルの変化を、Fisher情報行列で基準化した総合指標です。一般化線形モデルにおけるCook距離 $D_i$ の定義式は次の形をとります。

$$
D_i = \frac{(\hat{\boldsymbol{\beta}} – \hat{\boldsymbol{\beta}}_{(i)})’ \hat{J} (\hat{\boldsymbol{\beta}} – \hat{\boldsymbol{\beta}}_{(i)})}{p}
$$

ここで $\hat{\boldsymbol{\beta}}_{(i)}$ は第 $i$ 観測を除いて再推定した係数ベクトル、$\hat{J}$ はFisher情報行列の推定量、$p$ はパラメータ数です。$D_i$ が大きい観測ほど、その観測の有無によって係数ベクトル全体が大きく変化することを示します。

完全な再推定はコストが高いため、実用上は1点削除近似がよく用いられます。Pearson残差 $r_{P,i}$ とレバレッジ $h_{ii}$ を用いた近似式は次のように書けます。

$$
D_i \approx \frac{r_{P,i}^2}{p} \cdot \frac{h_{ii}}{(1 – h_{ii})^2}
$$

この近似式はCook距離が「外れ度」と「レバレッジ」の積として表現されることを明示しています。残差が大きくレバレッジも高い観測が最も大きな $D_i$ を持ちやすくなります。閾値としては $4/n$ が簡便な基準としてよく使われるほか、$F(p, n-p)$ 分布の上側パーセント点を参照する方法もあります。

仮定：1点削除近似式はIRLSが収束しており、除外した1観測がモデル全体の重み構造を大きく変えないという条件のもとで精度を保ちます。強い影響点が複数存在する場合にはこの条件が崩れ、近似精度が低下します。

限界：Cook距離はすべての係数の変化を1つの数値に集約するため、どの係数が主に影響を受けているかの情報が失われます。また二項モデルとポアソンモデルとでは分散関数が異なるため、Pearson残差のスケールが異なり、閾値をモデル間で単純に比較することは適切ではありません。

   
   

（Fig1. レバレッジ（$h_{ii}$）とPearson残差の散布図（バブルサイズはCook距離）。三象限（安全・高レバレッジ・外れ値・影響点）を色分けして可視化しています。）

   
   

（Fig2. 観測番号ごとのCook距離棒グラフ。閾値ライン（$4/n$）を重ねて表示し、影響点として特定された観測を強調しています。）

DFBETAS：係数ごとの影響診断

Cook距離が係数ベクトル全体への影響を集約するのに対し、DFBETASは第 $i$ 観測が第 $j$ 係数に与える影響を個別に定量化します。$\mathrm{DFBETAS}_{j,i}$ の定義は次の通りです。

$$
\mathrm{DFBETAS}_{j,i} = \frac{\hat{\beta}_j – \hat{\beta}_{j(i)}}{\widehat{\mathrm{SE}}(\hat{\beta}_j)}
$$

分子は第 $i$ 観測の有無による第 $j$ 係数の変化量、分母は全データで推定した係数の標準誤差です。符号は変化の方向（正：除外により係数が増加、負：減少）を示します。絶対値の閾値として $2/\sqrt{n}$ がよく参照されます。

実計算では毎回 $n$ 回の再フィットを行うのではなく、Sherman-Morrison-Woodbury近似（あるいはIRLSの加重最小二乗更新式の逆行列更新）を利用することで計算量を削減できます。具体的には、情報行列 $(X’WX)$ の1行1列を除いた逆行列を既存の逆行列から近似的に求めることにより、係数変化量を閉形式で近似します。

Cook距離とDFBETASは相補的な診断ツールです。Cook距離で影響点の候補を絞り込んだ後、DFBETASにより「どの係数がその観測によって引き寄せられているか」を特定するという2段階のワークフローが実務上有効です。

仮定：DFBETASの解釈はモデルが正しく特定されており、標準誤差 $\widehat{\mathrm{SE}}(\hat{\beta}_j)$ が安定して推定されていることを前提とします。多重共線性が強い場合には標準誤差が膨張し、DFBETASの値が過大または過小に評価されます。

限界：説明変数の数 $p$ が多い場合、全観測×全係数のDFBETAS行列は $n \times p$ の規模になります。大規模なモデルではこの行列を一覧することは実用的でなく、事前にCook距離等で候補観測を絞った上でDFBETASを適用することが現実的です。

   
   

（Fig3. 観測×係数のDFBETASヒートマップ。色の濃淡で各観測が各係数に与える影響の方向と大きさを示します。）

マスキング効果と複数観測の同時影響

1点削除に基づくCook距離やDFBETASは、除外する観測が1つである場合の影響量を測ります。しかしデータに複数の影響点が存在するとき、これらが互いに影響を隠し合うマスキング効果が発生することがあります。具体的には、2つの影響点がほぼ同じ方向に係数を引っ張っているとき、どちらか一方だけを除外しても残りの影響点が引き続き係数を動かすため、除外した観測の1点削除指標は小さく見積もられます。

逆に、本来は外れ値でない観測が、近傍の影響点の存在によって外れ値として診断される現象をスワンピング効果と呼びます。マスキング効果とスワンピング効果はいずれも、1点削除近似の精度が低下する条件――すなわち複数の強い影響点がクラスターを形成している状況――で顕在化します。

複数影響点の検出には、逐次削除法（影響点を1点ずつ繰り返し除外して再フィットする）や前進探索法が有効です。前進探索法では、まず影響点を含まない可能性が高い小サブセットからフィットを開始し、残差の小さい観測を順次加えながらモデルを更新します。この過程でレバレッジや残差の推移を追うことで、影響点が加わる時点を検出できます。

   
   

（Fig4. マスキング効果の概念図：2つの影響点が互いに存在を隠す状況を、1点削除診断と2点同時削除診断の結果比較で図示しています。）

限界：逐次削除法は各ステップでモデルを再フィットするため計算コストが高く、観測数が多い場合には実行時間が問題となります。前進探索法はRではforwardパッケージ等で利用できますが、標準の一般化線形モデル関数には含まれません。

品質管理ドメインへの応用：製造不良率データの事例

製造業の品質管理では、ロットごとの検査個数と不良品数が記録されます。このデータへの自然なモデルは二項分布を仮定したロジスティック回帰であり、説明変数として製造条件（温度・圧力・原料ロット番号等）が用いられます。

影響点診断のワークフローは次の順序で実施されます。まずモデルを全データで推定し、Cook距離プロットを作成します。$D_i > 4/n$ を超える観測を影響点候補として特定したのち、レバレッジプロットで候補観測が説明変数空間での孤立点かどうかを確認します。続いてDFBETASプロットにより、どの製造条件パラメータが該当観測によって大きく変化しているかを特定します。

診断後のアクションは機械的な観測削除ではなく、まず当該ロットの製造記録・測定値の再確認から始まります。記録誤りや測定機器の異常が確認された場合は観測値の修正・除外を検討し、物理的に異常な製造条件が原因であれば、そのロットを特殊条件として説明変数に追加するか、サブグループごとのモデルを検討します。最終的なモデルの妥当性は、影響点除外後の診断量の再確認と、実務的な解釈の整合性の両面から評価します。

仮定：この分析では各ロットが互いに独立であることを前提とします。同一製造ラインの連続ロット間に系列相関がある場合や、複数ラインが混在する場合には、クラスター構造を無視した二項モデルは前提を満たしません。その場合はGEE（一般化推定方程式）や混合効果モデルへの移行を検討する必要があります。

限界：ロット数が少ない小標本では $4/n$ 等の閾値の統計的根拠が弱くなります。経験則に基づく閾値を機械的に適用することは避け、Cook距離の相対的な大きさや製造記録との照合を優先すべきです。

頑健GLMへの接続と診断手法の比較

これまで述べたCook距離・DFBETASに基づく診断は、いずれもpost-hocな手続きです。すなわちまず通常の一般化線形モデルを推定し、診断量を計算し、影響点を特定した後に観測を除外して再推定するという繰り返しを前提とします。この手続きには「何回除外すれば十分か」「除外後のモデルが安定しているか」を逐次確認する負担があります。

これに対して頑健GLMは、推定段階から影響点への感度を低減させる設計を持ちます。M推定に基づく頑健GLMでは、IRLSの各ステップで用いる重みを残差の大きさに応じて下方修正します。具体的には影響関数 $\psi(r)$（Huber関数・Biweight関数など）を用いて、大きな残差を持つ観測の重みを自動的に縮小します。Huber関数は閾値 $c$ 以内の残差には線形の重みを保ち、それを超える残差には一定の重みを割り当てます。Biweight関数はさらに極端な外れ値の影響をほぼゼロにします。

診断ベースの削除と頑健推定の選択は、目的と文脈によって異なります。診断に基づく削除は各観測の扱いが透明であり、削除理由を製造記録等の実務情報と照合できる点で解釈性が高いと言えます。一方、頑健推定は影響点の特定と除外の繰り返しを省き、安定した係数推定を一括して得られる点で効率的です。ただし係数の解釈が加重最尤推定量として標準的な一般化線形モデルと厳密には異なり、推定量の漸近性質も通常の最尤推定量とは異なります。

限界：頑健GLMはRの標準glm関数には実装されておらず、robustbaseパッケージのglmrob等の追加パッケージが必要です。また頑健推定を用いる場合であっても、分布族とリンク関数の選択はデータの生成過程に即して行わなければならず、誤ったモデル設定は頑健推定によっても補正されません。

OLS（最小二乗法）の影響診断（古典的Cook距離・DFFITS）との比較では、一般化線形モデル版診断量はいずれも加重構造を持つ点が根本的な相違です。OLSのDFFITSは予測値の変化を基準化した指標ですが、一般化線形モデルでは予測値がリンク関数を経た非線形スケールに存在するため、直接の数値比較には注意が必要です。