Tweedie分布族の動機と位置づけ
自動車保険の年間損害額データでは、契約者の大多数がその期間中に保険請求を行わないため、応答変数の相当割合が正確にゼロとなります。請求が発生した場合の損害額は連続的な正値をとり、右裾が重い分布を示します。この「点質量を持つゼロ」と「連続正値」が混在する構造は、正規分布を仮定する線形モデルでも、ポアソン回帰でも、ガンマ回帰でも、単独では適切にモデル化できません。
Tweedie分布族は、べき分散関数によって定義される指数型分布族の部分族であり、この混在構造を単一モデルとして扱います。べき指数$p$を連続的に変化させることで、複数の分布族が統一的に包含されます。
$$
\mathrm{Var}(Y) = \phi \cdot \mu^p
$$
ここで$\phi > 0$は分散パラメータ、$\mu = E[Y]$は平均です。$p$の値と対応する分布族は次の通りです。$p=0$のとき正規分布、$p=1$のときポアソン分布、$1 < p < 2$のとき複合ポアソン-ガンマ分布、$p=2$のときガンマ分布、$p=3$のとき逆ガウス分布に対応します。
Tweedie分布族は準尤度の枠組みと密接に関連します。準尤度は平均と分散の関係$V(\mu) = \mu^p$のみを指定して推定方程式を構成しますが、$p \in \{0\} \cup [1, \infty)$の範囲ではこの準尤度に対応する真の確率密度が存在し、Tweedie分布族として同定されます。したがってTweedie GLMは準尤度アプローチの特殊ケースではなく、完全な確率モデルとして最尤推定が可能です。
（Fig1. べき指数 p の連続変化に伴う Tweedie 分布形状の変遷（p=1.0, 1.5, 2.0 のサンプル密度比較））
指数型分布族としてのTweedie分布
Tweedie分布は指数型分布族の正準形式に収まります。$1 < p < 2$の場合、閉形式の確率密度は存在しませんが、ラプラス変換の逆変換を級数展開することで密度の数値評価が可能です。確率密度はラプラス変換経由で次のように表現されます。
$$
f(y;\theta,\phi) = \exp\!\left\{\frac{y\theta – \kappa(\theta)}{\phi} + c(y,\phi)\right\}
$$
ここで$\theta$は自然パラメータ、$\kappa(\theta)$はキュムラント母関数、$c(y,\phi)$は正規化定数に相当する項です。$p \neq 1, p \neq 2$のとき、自然パラメータ$\theta$と平均$\mu$の関係は次式で与えられます。
$$
\theta = \frac{\mu^{1-p}}{1-p}
$$
累積量母関数は$\kappa(\theta) = \mu^{2-p}/(2-p)$であり、$\kappa”(\theta) = V(\mu) = \mu^p$が分散関数に対応します。正準リンク関数は$\kappa'(\theta) = \mu$を$\theta$で解くことで導出され、次式となります。
$$
g(\mu) = \frac{\mu^{1-p}}{1-p} \quad (p \neq 1,\; p \neq 2)
$$
分散関数$V(\mu) = \mu^p$は指数型分布族における分散関数の条件を満たし、GLMの理論的枠組みに適合します。実務ではlogリンク$g(\mu) = \log(\mu)$が最も広く用いられます。logリンクは正準リンクではありませんが、係数を乗法的に解釈できる利点があり、IRLSによる数値的収束の安定性にも問題はありません。
仮定：分散パラメータ$\phi > 0$が必要です。また$p$が$0 < p < 1$の範囲では対応する確率分布が存在しないため、Tweedie分布族は$p \in (-\infty, 0] \cup [1, \infty)$の範囲でのみ定義されます。
| p値 | 分布名 | 正準リンク | 典型応用 |
|---|---|---|---|
| $p = 0$ | 正規分布 | $\mu$ | 連続量・計測値 |
| $p = 1$ | ポアソン分布 | $\log\mu$ | カウントデータ・事故件数 |
| $1 < p < 2$ | 複合ポアソン-ガンマ分布 | $\mu^{1-p}/(1-p)$（実務はlog） | 保険純保険料・電力消費量 |
| $p = 2$ | ガンマ分布 | $1/\mu$（実務はlog） | 損害額・待ち時間 |
| $p = 3$ | 逆ガウス分布 | $\mu^{-2}$ | 到達時間・信頼性工学 |
| $p > 2$ | 安定分布の部分族 | $\mu^{1-p}/(1-p)$ | 重裾分布・極値統計 |
複合ポアソン-ガンマ表現（1 < p < 2）
$1 < p < 2$のTweedie分布は、複合ポアソン過程として構成的に解釈できます。事故件数$N$がポアソン分布$N \sim \mathrm{Poisson}(\lambda)$に従い、各事故の損害額$X_1, X_2, \ldots$が独立かつ同一のガンマ分布$X_i \sim \mathrm{Gamma}(\alpha, \beta)$に従うとき、総損害額は次式で表されます。
$$
Y = \sum_{i=1}^{N} X_i \quad (Y = 0 \text{ if } N = 0)
$$
複合ポアソン-ガンマのパラメータ$(\lambda, \alpha, \beta)$とTweedieパラメータ$(\mu, \phi, p)$の変換式は次の通りです。
$$
\lambda = \frac{\mu^{2-p}}{\phi(2-p)}, \quad \alpha = \frac{2-p}{p-1}, \quad \beta = \phi(p-1)\mu^{p-1}
$$
ゼロ確率は$P(Y=0) = P(N=0) = e^{-\lambda}$であり、$\mu, \phi, p$を用いると次のように表されます。
$$
P(Y = 0) = \exp\!\left(-\frac{\mu^{2-p}}{\phi(2-p)}\right)
$$
$\phi$が大きいほど（または$\mu$が小さいほど）ゼロの割合が高くなります。この構成は保険数理における頻度と重大度の分解と自然に対応します。保険における純保険料は「一定期間あたりの期待損害額」であり、$E[Y] = \lambda \cdot \alpha/\beta = \mu$として乗法モデルに整合します。Tweedie GLMではlogリンクを通じて$\log(\mu_i) = \mathbf{x}_i^T\boldsymbol{\beta}$と表現することで、説明変数の効果が乗法的に解釈されます。
（Fig2. 複合ポアソン-ガンマ構造の模式図：ゼロ質量と連続正値部分の混在（p=1.6 のシミュレーション例））
仮定：個別損害額$X_i$は互いに独立かつ同一分布に従うことが前提です。また事故件数$N$と各損害額$X_i$は独立である必要があります。
限界：$p \in (1,2)$の範囲外（$p > 2$など）ではTweedie分布は逆ガウス分布など別の解釈が必要であり、複合ポアソン-ガンマの構成的解釈は成立しません。また分散パラメータ$\phi$が非常に小さい場合、$\lambda$が非常に大きくなり、数値計算上の密度評価が不安定になる場合があります。
パラメータ推定：p・φ・β の最尤推定
Tweedie GLMのパラメータ推定は、回帰係数$\boldsymbol{\beta}$、分散パラメータ$\phi$、べき指数$p$の三者を対象とします。$1 < p < 2$では閉形式の密度が存在しないため、対数尤度はラプラス変換の逆変換から得られる無限級数近似を用いて評価されます。
$$
\ell(\boldsymbol{\beta}, \phi, p \mid \mathbf{y}) = \sum_{i=1}^{n} \left[\frac{y_i \theta_i – \kappa(\theta_i)}{\phi} + \log W(y_i, \phi, p)\right]
$$
ここで$\theta_i = \mu_i^{1-p}/(1-p)$、$\kappa(\theta) = \mu^{2-p}/(2-p)$であり、$W(y, \phi, p)$は正規化定数に相当する無限級数項です。この級数の数値評価にはDunn & Smythによるアルゴリズムが用いられており、Rのtweedicパッケージに実装されています。$p$を既知として固定した場合、$\boldsymbol{\beta}$の推定は通常のGLMと同様にIRLS（反復加重最小二乗法）で実施されます。logリンクの場合、作業応答$z_i$と作業重み$W_i$は次式で与えられます。
$$
z_i = \hat{\eta}_i + (y_i – \hat{\mu}_i)\frac{d\eta_i}{d\mu_i}, \quad W_i = \frac{\hat{\mu}_i^{2-p}}{\phi}
$$
収束後の$\phi$はPearson推定量で求められます。
$$
\hat{\phi} = \frac{1}{n – q} \sum_{i=1}^{n} \frac{(y_i – \hat{\mu}_i)^2}{\hat{\mu}_i^p}
$$
ここで$q$は推定パラメータ数です。$p$が未知の場合、プロファイル対数尤度$\ell_p(p) = \max_{\boldsymbol{\beta}, \phi} \ell(\boldsymbol{\beta}, \phi, p \mid \mathbf{y})$を$p$の格子点上（例：$p \in \{1.01, 1.05, \ldots, 1.99\}$）で評価するgrid searchが実用的です。最適化アルゴリズムは計算コストを下げますが、プロファイル尤度面の非凸性から局所解に収束するリスクがあります。
（Fig3. 保険データを用いたプロファイル対数尤度による p 推定：最適 p とその 95% 信頼区間）
仮定：$p$が既知として固定された場合と未知として推定対象とする場合では、AIC計算時のパラメータ数カウントが異なります。後者では$p$の1自由度を追加します。推定戦略の選択はモデルの目的に応じて明示的に決定する必要があります。
限界：$p$が境界値1または2に近い場合、無限級数項$W(y, \phi, p)$の収束が遅くなり、尤度評価の数値精度が低下します。またサンプルサイズが小さい場合（$n < 200$程度）、$p$の推定値は不安定となり、プロファイル尤度曲線が平坦化する傾向があります。
モデル診断と適合評価
Tweedie GLMの適合評価には逸脱度、Pearson残差、および量子化残差を用います。逸脱度は飽和モデルとの対数尤度比の2倍として定義され、$p \neq 1, p \neq 2$の一般的なTweedieの場合は次式で表されます。
$$
D(y;\hat{\mu}) = 2\sum_{i=1}^{n}\left[\frac{y_i(y_i^{1-p} – \hat{\mu}_i^{1-p})}{1-p} – \frac{y_i^{2-p} – \hat{\mu}_i^{2-p}}{2-p}\right]
$$
ゼロ観測値（$y_i = 0$）については$D_i = 2\hat{\mu}_i^{2-p}/[(2-p)\phi]$として別途加算します。$p = 1.5$の場合に展開すると次のようになります。
$$
D_i\big|_{p=1.5} = 2\left[\frac{y_i(y_i^{-0.5} – \hat{\mu}_i^{-0.5})}{-0.5} – \frac{y_i^{0.5} – \hat{\mu}_i^{0.5}}{0.5}\right] = 2\left[-2y_i(y_i^{-0.5} – \hat{\mu}_i^{-0.5}) – 2(y_i^{0.5} – \hat{\mu}_i^{0.5})\right]
$$
Pearson残差と逸脱度残差はそれぞれ次のように定義されます。
$$
r_i^P = \frac{y_i – \hat{\mu}_i}{\sqrt{\hat{\mu}_i^p}}, \quad d_i = \mathrm{sign}(y_i – \hat{\mu}_i)\sqrt{D_i}
$$
Tweedie分布の非対称性から、通常の残差プロットでは系統的なパターンが誤解を招く場合があります。無作為化量子化残差を用いることで、理論的に正しいモデルのもとでこの残差が標準正規分布に従うため、QQプロットによる視覚的診断が有効です。$\phi$の過推定はモデルの分散を実際より大きく見積もり区間推定が過度に広くなります。$\phi$の過少推定は逆に標準誤差を小さく見積もり、検定の第一種過誤を膨らませます。AICは$p$固定時には$\mathrm{AIC} = -2\hat{\ell} + 2(q+1)$として計算されますが、$p$を推定対象とする場合はさらに1自由度を加えた$\mathrm{AIC} = -2\hat{\ell} + 2(q+2)$を用います。
（Fig4. Tweedie GLM の量子化残差 QQ プロット：モデル適合診断の例）
限界：標準的な残差プロットはTweedieの非対称な分布形状に対して誤解を招く評価を与える可能性があり、量子化残差との併用が推奨されます。
保険統計への応用：純保険料モデルの実装
自動車保険の純保険料予測において、Tweedie GLMは標準的なアプローチとして広く採用されています。応答変数$Y_i$は被保険者$i$の年間損害額（ゼロ含む）、エクスポージャー$e_i$は契約年数です。logリンクとオフセットを用いたモデルは次式で表されます。
$$
\log(\mu_i) = \log(e_i) + \mathbf{x}_i^T\boldsymbol{\beta}
$$
ここで$\log(e_i)$はオフセット項であり、係数は持たず固定値として扱われます。これにより$\mu_i / e_i$（単位エクスポージャーあたりの期待損害額）をモデル化することと等価になります。係数の解釈は乗法的であり、$j$番目の説明変数が1単位増加したとき、期待損害額の乗数として$\exp(\hat{\beta}_j)$が対応します。
$$
\frac{\mu_i(x_j + 1)}{\mu_i(x_j)} = \exp(\hat{\beta}_j)
$$
実務では$p \approx 1.5$–$1.8$の範囲で推定されることが多く観察されています。これは事故頻度と損害額の両面での変動が反映された結果であり、$p$が1に近いほど頻度効果が支配的、$p$が2に近いほど損害額の変動が支配的な状況に対応します。等級制度との整合性については、ゼロ確率と正値期待値が同一パラメータ$\mu$によって制御されるTweedieの性質が、等級割引の乗法的構造と自然に対応します。二段階モデルと比較すると、Tweedie単一モデルは推定パラメータ数が少なく実装が簡潔ですが、頻度と重大度に異なる説明変数を割り当てる柔軟性が失われます。
限界：損害額分布に極端な大規模損害が含まれる場合、ガンマ分布の仮定が成立せず、Tweedieモデルは裾の重い損害を過小評価する可能性があります。またTweedieは損害額ゼロ以外の連続部分の分布形状を暗黙的にガンマ分布に制約するため、実際の損害額分布がガンマから大きく逸脱する場合は注意が必要です。
関連モデルとの比較・拡張
Tweedie GLMは複数の隣接モデルと交差する位置に存在し、適用場面の明確化が重要です。二段階モデルとの比較では、Tweedieは単一モデルとして推定の効率性と実装の簡潔さで優位ですが、ゼロ発生プロセスと正値プロセスに異なる説明変数や分布を指定する柔軟性では二段階モデルが優れます。Tweedieは事実上、ゼロの多寡が$\mu$と$\phi$と$p$によって制御されるという制約を課していることに留意が必要です。
ZIPやZINBとの本質的な違いは、ゼロの生成機構の解釈にあります。ZIP/ZINBはゼロを「構造ゼロ」と「偶発ゼロ」の混合として明示的にモデル化します。Tweedieの$1 < p < 2$領域ではゼロは確率的ゼロ（複合ポアソン過程で$N=0$となった場合）として解釈され、構造ゼロの概念は含まれません。負の二項回帰との接続は$p \to 1$の極限で確認できます。$p=1$はポアソン分布に対応しますが、$p$を1からわずかに大きくすると過分散の連続的な導入が可能であり、負の二項回帰の分散関数$V(\mu) = \mu + \mu^2/k$とは異なる分散構造をとります。発展的な方向として、Tweedie過程（時系列版）は再保険数理や損害累積モデルへの応用が進んでいます。
| モデル | ゼロの扱い | 分散構造 | 推定複雑度 | 推奨場面 |
|---|---|---|---|---|
| Tweedie GLM（$1<p<2$） | 確率的ゼロ（$N=0$） | $\phi\mu^p$（べき） | 中（$p$推定が加わる場合は高） | 保険純保険料・ゼロ混在連続データ |
| ポアソン回帰 | 自然（$\lambda$依存） | $\phi\mu$（等分散） | 低 | 整数カウントデータ・過分散なし |
| ガンマ回帰 | 対応不可 | $\phi\mu^2$（CV一定） | 低 | 正値連続データのみ（ゼロなし） |
| 正規GLM | 対応不可（原理上） | $\phi$（定数） | 最低 | 対称連続データ・裾が軽い場合 |
| 負の二項回帰 | 自然（$\lambda$依存） | $\mu + \mu^2/k$ | 中（$k$の推定） | 過分散カウントデータ |
| ZIP | 構造ゼロ＋確率的ゼロ | 混合分布由来 | 高（2成分推定） | ゼロの生成機構が2種類存在する場合 |
| 二段階モデル | 明示的モデル化 | 2モデル独立 | 高（2モデル独立推定） | 頻度と重大度に異なる変数が必要な場合 |
限界：$p > 2$の領域（逆ガウス分布など）では複合ポアソン-ガンマの解釈は成立せず、分布の性質が大きく異なるため別途検討が必要です。また$\mu$が非常に小さい予測領域では、予測区間の計算に用いる無限級数評価のコストが増大する場合があります。
