GLM の基礎

GLM タブで使われている統計理論の背景です。操作方法は GLM のページを参照してください。

モデルの定式化

GLM（一般化線形モデル）は、二値・カウント・正の連続値のように正規分布が合わない応答変数を回帰で扱うために、正規線形モデルを次節で定義する指数型分布族へ一般化した枠組みです。Nelder & Wedderburn (1972) が提唱しました。3つの要素で定義されます:

分布族: 応答変数 $Y$ の分布が指数型分布族に属する
線形予測子: $\eta = X\beta$ （説明変数の線形結合）
リンク関数: 単調関数 $g$ により $\eta = g(\mu)$ として線形予測子と平均 $\mu = E[Y]$ を結びつける

指数型分布族

応答変数の分布は、二値なら Binomial、カウントなら Poisson と、データの種類ごとに変わります。GLM が分布ごとに推定の理論を作り直さずに済むのは、分布を指数型分布族に限っているためです。この分布族では、分布ごとの違いが平均と分散の関係を表す 1 つの関数に集約され、それを差し替えるだけで同じ推定アルゴリズム（IRLS）が使えます。

確率密度（質量）関数が次の形で書ける分布族を指数型分布族と呼びます:

f(y \mid \theta, \phi) = \exp\!\left\{\frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi)\right\}

この式に現れる $\theta$ を自然パラメータ（canonical parameter）、 $\phi$ を分散パラメータ、 $b(\theta)$ を対数分配関数と呼びます。 $a(\phi)$ は各分布族ごとに定まる既知の関数です（下の表を参照）。この定義から、平均と分散が対数分配関数 $b(\theta)$ の微分として定まります:

$E[Y] = b'(\theta) = \mu$
$\operatorname{Var}(Y) = b''(\theta) \cdot a(\phi)$

$b''(\theta)$ を $\mu$ の関数として書き直したものが分散関数 $V(\mu)$ で、平均と分散の関係 $\operatorname{Var}(Y) = V(\mu) \cdot a(\phi)$ を 1 つの関数で表します。各分布族の $V(\mu)$ は分散関数と過分散の表に示します。

各分布族のパラメータを以下に示します。

分布族	$\theta$ （自然パラメータ）	$a(\phi)$	$b(\theta)$	$c(y, \phi)$
Gaussian	$\mu$	$\phi$	$\theta^2/2$	$-\dfrac{y^2}{2\phi} - \dfrac{\log(2\pi\phi)}{2}$
Binomial	$\log\!\bigl(\mu/(1-\mu)\bigr)$	$1/n_i$	$\log(1+e^\theta)$	$\log\binom{n_i}{k_i}$
Poisson	$\log\mu$	$1$	$e^\theta$	$-\log(y!)$
Gamma	$-1/\mu$	$\phi$	$-\log(-\theta)$	$(1/\phi - 1)\log y + (1/\phi)\log(1/\phi) - \log\Gamma(1/\phi)$
Negative Binomial	$\log\!\bigl(\mu/(\mu+r)\bigr)$	$1$	$-r\log(1-e^\theta)$	$\log\Gamma(y+r) - \log\Gamma(r) - \log(y!)$

Binomial の $y_i$ は成功割合 $y_i = k_i/n_i$ , $0 \le y_i \le 1$ です。 $k_i$ は観測 $i$ の成功回数、 $n_i$ は試行回数（サンプルサイズ $n$ ではなく各観測の試行回数）、 $\mu$ は成功確率です。 $n_i=1$ のとき Bernoulli 分布に帰着します
Negative Binomial の $r$ は MIDAS の操作画面では $\theta$ と表記されていますが、このページでは指数型分布族の自然パラメータとの混同を避けるため $r$ で統一しています。 $r$ が既知の場合のみ指数型分布族に属します。MIDAS の自動推定モードでは、 $\beta$ をプロファイルアウトした尤度 $L_p(r) = \max_\beta L(\beta, r)$ を外側ループで最大化して $r$ を推定します（GLM の操作方法を参照）。自動推定時に報告される $\hat\beta$ の標準誤差は $r = \hat r$ を既知として固定した情報行列から計算されるため、 $r$ の推定不確実性は含まれません

リンク関数の選択

リンク関数は線形予測子 $\eta$ と応答変数の期待値 $\mu$ を結びつける単調関数 $\eta = g(\mu)$ です。 $g(\mu) = \theta$ （自然パラメータ）とするリンクを正準リンク（canonical link）と呼びます。分布族が応答のデータ型でほぼ決まるのに対し、リンク関数には選択の余地があります。分布族の決め方は分布族の選び方を参照してください。

選択の第一の基準は効果尺度、つまり係数を応答の差として読みたいか、比として読みたいかです。Identity リンクでは係数がそのまま応答の差を表します。Log リンクでは線形予測子上の足し算が平均の掛け算になるため、 $\exp(\beta)$ が比として読め、Poisson や Negative Binomial では率比（IRR）にあたります。Logit では $\exp(\beta)$ がオッズ比（OR）です。Logit と Log では、係数テーブルにこの解釈に対応する $\exp(\hat\beta)$ の列が表示されます。

リンク関数	数式	係数の解釈	正準リンクとなる分布族
Identity	$\eta = \mu$	$\beta$ は応答の差	Gaussian
Logit	$\eta = \log\!\bigl(\mu / (1 - \mu)\bigr)$	$\exp(\beta)$ はオッズ比	Binomial
Log	$\eta = \log(\mu)$	$\exp(\beta)$ は比（乗法的効果）	Poisson
Inverse	$\eta = 1/\mu$	$\beta$ は $1/\mu$ スケールの変化	Gamma
Probit	$\eta = \Phi^{-1}(\mu)$	$\beta$ は $\Phi^{-1}(\mu)$ （確率の z 値）の変化	—

どの効果尺度で読むのが適切かを最も強く裏づけるのは、データの生じ方を知っていることです。各要因の効果が発生率に倍率として掛かる過程なら、効果は log の尺度で足し算になります。観測されない連続量が閾値を超えると応答が 1 になる過程を想定できるなら、係数はその連続量の尺度上の効果として読め、連続量に正規分布を仮定した場合が Probit です（GLMM の基礎の閾値モデルを参照）。

正準リンクは、効果尺度で決まらないときの既定の選択肢にあたります。 $\eta = \theta$ となるため $X'y$ が $\beta$ の十分統計量になり、対数尤度が $\beta$ について凹になります。計画行列 $X$ がフルランクかつ最尤推定量が存在すれば、その解は一意になり、IRLS の収束も安定します¹。これらは推定の扱いやすさの性質であって、正準であることがモデルの正しさを保証するわけではありません。

非正準リンクではこれらの性質が保証されません。それでも効果尺度を優先して選ばれることがあります。たとえば Gamma 分布族の正準リンクは Inverse（ $\eta = 1/\mu$ ）ですが、係数が $1/\mu$ スケールになるため解釈が難しく、実務では $\exp(\beta)$ を乗法的効果として解釈できる Log リンクがよく使われます。

効果尺度で選んだリンクは、分布族の定義域と両立するとは限りません。たとえば率の差を読みたい Poisson + Identity では、線形予測子が負になる観測があると平均が定義域（正の値）を外れます。外れるかどうかは組み合わせではなくデータで決まり、外れた場合 MIDAS は平均を範囲内に押し込めず適合を中止してエラーを表示します（GLM の操作方法の注意事項を参照）。

パラメータ推定（IRLS）

GLM のパラメータ $\beta$ は最尤法で推定します。正則条件（対数尤度の微分可能性、真のパラメータがパラメータ空間の内点であることなど）のもとで推定量は一致性・漸近正規性・漸近有効性を持ちます。一般には解析的に解けないため、IRLS（Iteratively Reweighted Least Squares）で数値的に求めます²。

各反復で作業用重み $W$ と調整従属変数 $z$ を計算し、加重最小二乗:

\hat\beta^{(t+1)} = (X'W^{(t)}X)^{-1}X'W^{(t)}z^{(t)}

を解いて $\beta$ を更新します。 $W$ は対角行列であり、その第 $i$ 対角要素 $W_{ii}$ と $z$ の第 $i$ 成分 $z_i$ は現在の $\hat\mu^{(t)}$ とリンク関数から次のように計算されます:

W_{ii} = \frac{1}{V(\mu_i)\,(d\eta/d\mu)_i^2}, \qquad z_i = \eta_i + (y_i - \mu_i)\,\Bigl(\frac{d\eta}{d\mu}\Bigr)_i

$V(\mu)$ は分散関数、 $d\eta/d\mu$ はリンク関数の導関数です。Grouped Binomial では、指数型分布族の表の $a(\phi) = 1/n_i$ のとおり観測 $i$ の分散が $V(\mu_i)/n_i$ になるため、 $W_{ii}$ にさらに試行回数 $n_i$ を掛けます。GLM に対する IRLS の定式化は Nelder & Wedderburn (1972) を参照してください。係数の変化量が収束閾値を下回ると終了です。

正準リンクを使う場合、対数尤度の凹性からこの反復は安定して収束します。非正準リンクでは収束が不安定になることがあるため、反復回数の増加や収束失敗に注意してください。

係数の標準誤差と信頼区間

GLM タブが報告する係数の信頼区間は Wald 型 $\hat\beta_j \pm c \cdot \operatorname{SE}(\hat\beta_j)$ です。標準誤差 $\operatorname{SE}(\hat\beta_j)$ は、推定共分散行列 $\hat\phi\,(X'\hat{W}X)^{-1}$ の第 $j$ 対角成分の平方根です。 $\hat W$ は IRLS の重み行列 $W$ を収束時の $\hat\mu$ で評価したものです。

推定分散パラメータ $\hat\phi$ は分布族ごとに決まります。Poisson・Binomial・ $r$ を推定する Negative Binomial では $\hat\phi = 1$ を使います。Gaussian では残差 deviance を $n-p$ で割った値を使います。Gaussian の deviance は $\sum_i (y_i - \hat\mu_i)^2$ なので、この値は $\text{Pearson }\chi^2/(n-p)$ と一致します。Gamma と $r$ を固定した Negative Binomial では $\text{Pearson }\chi^2/(n-p)$ を使います。Gamma で deviance ベースの推定量を使わないのは、一致推定量にならないためです（McCullagh & Nelder, 1989）。 $r$ 固定時に $\hat\phi$ を推定する意味は分散関数と過分散で述べます。

分位点 $c$ は、 $\hat\phi$ をデータから推定する分布族では自由度 $n-p$ の $t$ 分布から、 $\hat\phi = 1$ の分布族では標準正規分布から取ります。 $t$ 分布を使う分布族では、 $\hat\phi$ の推定に伴う不確実性が $c$ に反映されます。

Gaussian + Identity では、この区間の被覆確率（区間が真の値を含む確率）が有限標本で名目水準に一致します。それ以外の分布族・リンクの組み合わせでは、区間は最尤推定量の漸近正規性に基づく近似で、被覆確率が名目水準に一致するのは大標本に限られます。

分散関数と過分散

指数型分布族で述べたとおり、分散関数 $V(\mu) = b''(\theta)$ は対数分配関数の二階微分を $\mu$ で書き直したものです。 $\operatorname{Var}(Y) = V(\mu) \cdot a(\phi)$ の関係を通じて、各分布族の平均と分散の関係を規定します。

分布族	$V(\mu)$	$a(\phi)$	$\operatorname{Var}(Y)$
Gaussian	$1$	$\phi$	$\phi$ （= $\sigma^2$ ）
Binomial	$\mu(1 - \mu)$	$1/n_i$	$\mu(1-\mu)/n_i$
Poisson	$\mu$	$1$	$\mu$
Gamma	$\mu^2$	$\phi$	$\mu^2 \phi$
Negative Binomial	$\mu + \mu^2/r$	$1$	$\mu + \mu^2/r$

Poisson と Binomial では分散パラメータ $\phi = 1$ と仮定します。実データの分散がこの仮定より大きい場合を過分散（overdispersion）と呼びます。過分散があると標準誤差が過小推定され、信頼区間が狭くなりすぎます。過分散の診断には、View Diagnostics で開く GLM Diagnostics タブの Deviance Goodness-of-Fit に表示される Deviance/df 比を用います。仮定が正しければ 1 前後になるはずなので、1 から大きく離れる場合は過分散を疑います。 $n-p$ が小さいほど偶然のばらつきで 1 から離れやすい点には注意してください。

Poisson で過分散が検出された場合、Negative Binomial に切り替えることで分散に $\mu^2/r$ の項が加わり、過分散を明示的にモデル化できます。 $r$ を推定する場合、過分散は $r$ に吸収されるため $\phi = 1$ として扱います。 $r$ を固定する場合は、固定した $r$ では捉えきれない余剰分散を $\hat\phi = \text{Pearson }\chi^2/(n-p)$ で推定し、標準誤差と信頼区間に反映します。

ただし試行回数 $n_i = 1$ の二値データ（ロジスティック回帰）では、各観測が Bernoulli $(\mu_i)$ に従い、平均 $\mu_i$ が決まれば周辺分散 $\mu_i(1-\mu_i)$ も一意に決まります。観測レベルで分散に自由度がないため、「データの分散が理論分散より大きい」という対比を測る対象がそもそもありません。個体レベルのデータで Pearson $\chi^2$ や逸脱度から過分散を検出できないのはこのためです。これは「過分散がない」ことを意味するのではなく「同じデータからは検出できない」という意味で、クラスタや繰り返し測定（例: 同一施設に属する複数の患者を独立と扱う場合）に由来する余剰分散は別枠で存在することがあり、その場合は別途対応が必要です（用語集を参照）。過分散の検出と古典的な対処が意味を持つのは $n_i > 1$ の Grouped Binomial です。

Grouped Binomial で過分散が検出された場合、MIDAS には現在 quasi-binomial や Beta-Binomial などの対処法がありません。クラスタ構造由来の余剰分散であれば GLMM でランダム効果を導入する選択肢があります。過分散が疑われる場合は、分散パラメータの推定値を確認し、標準誤差や信頼区間が過小推定されている可能性を考慮してください。

平均応答の信頼区間と予測区間

GLM の予測機能が出力する、平均応答の信頼区間と新規観測の予測区間の数理的背景です。

以下の公式に現れる $\hat\phi$ は、係数の標準誤差と信頼区間で定義した推定分散パラメータです。分位点 $c$ も同節と同じ規則で選びます。 $h_i = x_\text{new}' (X'\hat W X)^{-1} x_\text{new}$ は予測点のレバレッジであり、予測に使う説明変数の値が学習データの中心からどの程度離れているかを示します。

平均応答の信頼区間は、リンクスケールの Wald 区間 $\hat\eta \pm c\sqrt{\hat\phi\, h_i}$ を $g^{-1}$ で応答スケールへ変換した区間です。 $g$ は単調なので、この変換は区間の被覆確率を変えません。リンクが非線形の場合、変換後の区間は $\hat\mu$ を中心に対称になりません。区間の正確さは係数の信頼区間と同じです。

予測区間は新規観測 $Y_\text{new}$ のばらつきを含む区間です。新規観測のばらつきはリンクスケールの区間の変換では表せないため、計算方法は分布族に依存します。Gaussian + Identity を除き、プラグイン法を使います。プラグイン法とは、推定したパラメータを真の値として扱い、その値に基づいて区間を計算する方法です。パラメータ推定の不確実性を含まない点が信頼区間との違いです。

Gaussian + Identity: 解析的公式 $\hat\mu \pm t_{n-p} \sqrt{\hat\phi(1 + h_i)}$ を使います。新規観測の分散（ $\hat\phi$ ）と平均の推定不確実性（ $\hat\phi \cdot h_i$ ）の両方を含みます
Gaussian + Identity 以外のリンク: プラグイン法 $\hat\mu \pm t_{n-p} \sqrt{\hat\phi}$ を使います。非線形のリンク変換により、平均 $\mu$ スケールでの推定不確実性を閉形式（closed form）で組み込めません。その結果、この組み合わせの予測区間は推定不確実性を反映しない簡易的なものとなり、データ中心の予測点と外挿点で同じ幅の区間になります。リンクスケールで区間を作って逆変換する方法も delta 法による一次近似も、リンクの非線形性や小標本で被覆確率の保証が弱くなるため採用していません
Poisson, Binomial, Gamma, Negative Binomial: プラグイン分位点法を使います。当てはめた分布の分位点を直接計算します
- Poisson: 平均 $\hat\mu$ の Poisson $(\hat\mu)$ の分位点
- Binomial: 成功確率 $\hat\mu$ と試行回数 $n_\text{new}$ （Grouped Binomial では予測データセットにある学習時の Trials 列の各行の値、二値データでは $1$ ）を使った Binomial $(n_\text{new}, \hat\mu)$ の分位点を $n_\text{new}$ で割った成功割合スケールの値
- Gamma: 平均 $\hat\mu$ ・形状 $\alpha = 1/\hat\phi$ ・スケール $\hat\phi \cdot \hat\mu$ の Gamma 分布の分位点
- Negative Binomial: 平均 $\hat\mu$ と $r$ （自動推定モードでは $\hat r$ 、固定モードでは指定値）を使った Negative Binomial 分布の分位点

離散分布（Poisson, Binomial, Negative Binomial）では分位点を保守側（ $P(X \le k) \ge \alpha$ を満たす最小の整数 $k$ ）に丸めます。このため、推定した分布を真の分布とみなしたときの被覆確率は名目信頼水準以上になります。個体 Binomial（ $n_i=1$ ）では分位点候補が $\{0, 1\}$ のみで、区間としての情報量は限られます。

プラグイン法はパラメータ推定の不確実性を含まないため、実際の被覆確率が設定した信頼水準を下回ることがあります。この傾向は小標本や観測データ範囲から離れた予測点で顕著になります。

参考文献

Nelder, J. A., & Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society: Series A, 135(3), 370-384. https://www.jstor.org/stable/2344614
McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman and Hall.

最尤推定量が存在しない例外が完全分離、つまりある説明変数の線形結合で応答を完全に分離できる状況で、対数尤度が凹かつフルランクでも有限の最大値を持ちません。二値応答のロジスティック回帰が典型例で、多項ロジットなど他の離散応答モデルでも類似のケースが起こります。MIDAS は分離を検知して警告します（GLM の操作方法の収束問題の項を参照）。 ↩
例外は Gaussian + Identity です。 $V(\mu)=1$ ・ $d\eta/d\mu=1$ から本文の式の $W=I$ ・ $z=y$ となり、重みがデータに依存しないため、任意の初期値から 1 回の反復で OLS 解 $\hat\beta = (X'X)^{-1}X'y$ に到達します。 ↩

このページの Markdown 版もあります。