GLM の基礎

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

モデルの定式化

GLM は正規線形モデルを指数型分布族に一般化した枠組みで、Nelder & Wedderburn (1972) が提唱しました。3つの要素で定義されます:

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

OLS は GLM の特殊ケース（Gaussian ファミリー + Identity リンク）です。この場合 IRLS は1回の反復で正規方程式の解に一致し、 $\phi$ をデータから推定する MIDAS の実装では Wald 統計量が $t_{n-p}$ に厳密に従うため、OLS の $t$ 検定と有限標本で等価になります。他のファミリー・リンクの組み合わせでは、Wald 検定は漸近的な近似にとどまります。

指数型分布族

GLM は分布族をこの形に限ることで、平均と分散の関係を $b(\theta)$ で統一的に記述し、ファミリーを問わない共通の推定アルゴリズム（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)$ は対数分配関数です。平均と分散は $b(\theta)$ から導かれます:

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

$b''(\theta)$ を $\theta$ ではなく $\mu$ の関数として書き直したものが分散関数 $V(\mu)$ です。つまり $\operatorname{Var}(Y) = V(\mu) \cdot a(\phi)$ が成り立ちます。たとえば Poisson では $b(\theta) = e^\theta$ なので $b'(\theta) = e^\theta = \mu$ 、 $b''(\theta) = e^\theta = \mu$ となり、 $V(\mu) = \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$	$\log(1+e^\theta)$	$\log\binom{n}{k}$
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$ は成功割合 $y = k/n$ , $0 \le y \le 1$ です。 $k$ は成功回数、 $n$ は試行回数、 $\mu$ は成功確率です。 $n=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$ の推定不確実性は含まれません（R の MASS::glm.nb と同じ扱い）

リンク関数

リンク関数は線形予測子 $\eta$ と応答変数の期待値 $\mu$ を結びつける単調関数 $\eta = g(\mu)$ です。 $g(\mu) = \theta$ （自然パラメータ）とするリンクを正準リンク（canonical link）と呼びます。

リンク関数	数式	正準リンクとなるファミリー
Identity	$\eta = \mu$	Gaussian
Logit	$\eta = \log\!\bigl(\mu / (1 - \mu)\bigr)$	Binomial
Log	$\eta = \log(\mu)$	Poisson
Inverse	$\eta = 1/\mu$	Gamma
Probit	$\eta = \Phi^{-1}(\mu)$	—

正準リンクには重要な性質があります。 $\eta = \theta$ となるため $X'y$ が $\beta$ の十分統計量になり、対数尤度が $\beta$ について凹になります。計画行列 $X$ がフルランクかつ最尤推定量が存在すれば、その解は一意になり、IRLS の収束も安定します。ただし完全分離（ある説明変数の線形結合で応答を完全に分離できるケース）では対数尤度が凹かつフルランクでも有限の最大値を持たず、最尤推定量は存在しません。二値応答のロジスティック回帰が典型例で、多項ロジットなど他の離散応答モデルでも類似のケースが起こります（GLM の操作方法の収束問題の項を参照）。

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

パラメータ推定（IRLS）

GLM のパラメータ $\beta$ は最尤法で推定します。正則条件のもとで推定量は一致性・漸近正規性・漸近有効性を持ちます。一般には解析的に解けないため、IRLS（Iteratively Reweighted Least Squares）で数値的に求めます（Gaussian + Identity は例外で、 $V(\mu)=1$ ・ $d\eta/d\mu=1$ から下式の $W=I$ ・ $z=y$ となり、重みがデータに依存しないため任意の初期値から1回の反復で OLS 解 $\hat\beta = (X'X)^{-1}X'y$ に到達します）。

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

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

を解いて $\beta$ を更新します。 $W$ と $z$ は現在の $\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$ はリンク関数の導関数です。GLM に対する IRLS の定式化は Nelder & Wedderburn (1972) を参照してください。係数の変化量が収束閾値を下回ると終了です。

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

分散関数と過分散

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

ファミリー	$b''(\theta)$	$V(\mu)$	$a(\phi)$	$\operatorname{Var}(Y)$
Gaussian	$1$	$1$	$\phi$	$\phi$ （= $\sigma^2$ ）
Binomial	$\dfrac{e^\theta}{(1+e^\theta)^2}$	$\mu(1 - \mu)$	$1/n$	$\mu(1-\mu)/n$
Poisson	$e^\theta$	$\mu$	$1$	$\mu$
Gamma	$1/\theta^2$	$\mu^2$	$\phi$	$\mu^2 \phi$
Negative Binomial	$\dfrac{re^\theta}{(1-e^\theta)^2}$	$\mu + \mu^2/r$	$1$	$\mu + \mu^2/r$

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

Poisson で過分散が検出された場合、Negative Binomial に切り替えることで分散に $\mu^2/r$ の項が加わり、過分散を明示的にモデル化できます。 $r$ を推定する場合は過分散が分散関数で既にモデル化されているため $\phi = 1$ として扱い、 $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$ は推定分散パラメータを表します。Gaussian では残差 deviance を $n - p$ で割った値（Gaussian では Pearson $\chi^2$ と恒等的に等しい）、Gamma では Pearson $\chi^2$ を $n - p$ で割った値です（deviance ベースは Gamma の一致推定量にならないため Pearson $\chi^2/(n-p)$ を使用します）。Poisson/Binomial では $\hat\phi = 1$ です。 $h_i = x_\text{new}^T (X^T \hat W X)^{-1} x_\text{new}$ は予測点のレバレッジであり、予測点が学習データの説明変数空間の中心からどの程度離れているかを示します。

プラグイン法とは、推定したパラメータを真の値として扱い、その値に基づいて区間を計算する方法です。パラメータ推定の不確実性を含まない点が信頼区間との違いです。

予測区間の計算方法はファミリーに依存します:

Gaussian + identity link: 解析的公式 $\hat\mu \pm t_{n-p} \sqrt{\hat\phi(1 + h_i)}$ を使用します。新規観測の分散（ $\hat\phi$ ）と平均の推定不確実性（ $\hat\phi \cdot h_i$ ）の両方を含みます
Gaussian + 非 identity link: プラグイン法 $\hat\mu \pm t_{n-p} \sqrt{\hat\phi}$ を使用します。 $t_{n-p}$ は選択した信頼水準に対応する $t$ 分布の分位点です。非線形のリンク変換により、平均 $\mu$ スケールでの推定不確実性を閉じた形で組み込めません。代替として (a) リンクスケールで区間 $\hat\eta \pm t_{n-p}\sqrt{\hat\phi(1+h_i)}$ を作って $g^{-1}$ で応答スケールに逆変換する方法や、(b) delta 法による一次近似 $\operatorname{Var}(\hat\mu) \approx (d\mu/d\eta)^2 \cdot \hat\phi h_i$ を使う方法もありますが、いずれもリンクの非線形性や小標本で被覆確率の保証が弱くなるため MIDAS では採用していません。その結果、この組み合わせの予測区間は推定不確実性を反映しない簡易的なものとなり、データ中心の予測点と外挿点で同じ幅の区間になります
Poisson, Binomial, Gamma, Negative Binomial: プラグイン分位点法を使用します。推定した分布パラメータを真の値として扱い、分位点を直接計算します
- Poisson: 平均 $\hat\mu$ の Poisson $(\hat\mu)$ の分位点
- Binomial: 成功確率 $\hat\mu$ 、予測点で指定する試行回数 $n_\text{new}$ を使った Binomial $(n_\text{new}, \hat\mu)$ の分位点
- Gamma: 平均 $\hat\mu$ ・形状 $\alpha = 1/\hat\phi$ の Gamma 分布の分位点
- Negative Binomial: 平均 $\hat\mu$ と $r$ （自動推定モードでは $\hat r$ 、固定モードでは指定値）を使った Negative Binomial 分布の分位点

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

プラグイン法はパラメータ推定の不確実性を含まないため、実際の被覆確率（区間が真の値を含む確率）が設定した信頼水準を下回ることがあります。この傾向は小標本や観測データ範囲から離れた予測点で顕著になります。小標本で被覆確率が重要な場面では、標本サイズを増やすか、より計算コストの高いブートストラップ等の方法を別途検討してください（MIDAS では現在サポートしていません）。

参考文献

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