23. 正則化GLM：Lasso・Ridgeの一般化線形モデルへの拡張

なぜGLMに正則化が必要か：高次元・多重共線性の問題

一般化線形モデル（GLM）のパラメータ推定は通常、最尤推定によって行われます。最尤推定量は漸近的に有効であり、標本数 $n$ が十分大きく予測変数の数 $p$ が小さい古典的な設定では安定した推定が得られます。しかし現代のバイオインフォマティクスや医療データ解析では、遺伝子発現量やゲノム変異といった特徴量の数 $p$ が患者数 $n$ を大幅に上回る状況（$p \gg n$）が日常的に発生します。この設定において、最尤推定は理論的・数値的の両面で機能しなくなります。

GLMにおけるスコア方程式を解くためには、Fisher情報行列 $\mathcal{I}(\beta) = X^\top W X$（ここで $W$ はワーキング重みの対角行列）が逆行列を持つ必要があります。しかし $p \geq n$ のとき、$X^\top W X$ は最大で $n$ のランクしか持てないため、行列は特異または劣化となります。劣化した行列とは、固有値のなかに0に近い値が存在し、数値的な逆行列計算が不安定になる状態を指します。

$$
\mathcal{I}(\beta) = X^\top W X \in \mathbb{R}^{p \times p}, \quad \operatorname{rank}(\mathcal{I}) \leq \min(n, p)
$$

$p \geq n$ の場合 $\operatorname{rank}(\mathcal{I}) < p$ となり、最尤推定量 $\hat{\beta}$ は一意に定まりません。多重共線性が存在する場合も同様に情報行列は事実上特異となり、推定量の分散がCramér-Rao下界を通じて発散します。

このような状況では推定量の分散が爆発し、訓練データへの過適合が生じます。モデルは訓練標本には高い尤度を与えますが、新たなデータへの予測精度は著しく低下します。この現象はバイアス・分散トレードオフの観点から整理できます。最尤推定量はバイアスをゼロに保つ代わりに分散を無制限に許容しますが、予測精度を最大化するためには適度なバイアスを許容して分散を抑制することが有益な場合があります。

限界：正則化はバイアスを意図的に導入することで分散を制御する手法です。すなわち、推定精度の向上は情報の一部を犠牲にすることによってのみ達成されます。どの程度のバイアスが許容されるかはデータの構造と分析目的に依存します。

ペナルティ付き対数尤度の定式化

正則化GLMの基本的枠組みは、通常の最尤推定の目的関数である対数尤度 $\ell(\beta)$ にペナルティ項を加えることで得られます。対数尤度の一般形は

$$
\ell(\beta) = \sum_{i=1}^{n} \left[ y_i \theta_i – b(\theta_i) \right] / \phi + c(y_i, \phi)
$$

です（ここで $\theta_i$ は自然パラメータ、$b(\cdot)$ は対数分配関数、$\phi$ は分散パラメータ）。これにペナルティを加えたペナルティ付き負対数尤度の最小化問題として定式化すると次のようになります：

$$
\hat{\beta} = \arg\min_{\beta} \left\{ -\ell(\beta) + P_\lambda(\beta) \right\}
$$

ペナルティ項 $P_\lambda(\beta)$ の選択によってRidge、Lasso、Elastic Netの3手法が生まれます。

Ridge GLMではL2ペナルティを採用します：$P_\lambda(\beta) = \lambda \|\beta\|_2^2 = \lambda \sum_j \beta_j^2$。係数は連続的に縮小されますがゼロには達しません。Lasso GLMではL1ペナルティを採用します：$P_\lambda(\beta) = \lambda \|\beta\|_1 = \lambda \sum_j |\beta_j|$。L1ペナルティは解をスパースにします。スパースとは多くの係数が厳密にゼロとなる状態を指し、変数選択が自動的に行われることを意味します。Elastic Netは両者を混合します：

$$
P_{\lambda,\alpha}(\beta) = \lambda \left[ \alpha \|\beta\|_1 + \frac{1-\alpha}{2} \|\beta\|_2^2 \right], \quad \alpha \in [0,1]
$$

$\alpha=1$ がLassoに、$\alpha=0$ がRidgeに対応します。正則化パラメータ $\lambda \geq 0$ は縮小の強度を制御し、$\lambda \to 0$ では通常の最尤推定量に収束し、$\lambda \to \infty$ ではすべての係数がゼロに縮小されます。

これらの定式化にはベイズ統計学的な解釈があります。L2ペナルティによるRidgeは、$\beta_j \sim \mathcal{N}(0, 1/(2\lambda))$ という正規事前分布を仮定した最大事後確率推定に等しいと言えます。L1ペナルティによるLassoは、$\beta_j \sim \operatorname{Laplace}(0, 1/\lambda)$ というラプラス事前分布を仮定した最大事後確率推定に対応します。ラプラス分布は原点での確率密度が高く尾が厚いため、最大事後確率解が厳密にゼロをとりやすくなります。

前提条件：正則化パラメータ $\lambda$（およびElastic Netの $\alpha$）はデータから交差検証等によって選択するハイパーパラメータであり、モデル内で推定される固定パラメータではありません。

座標降下法によるペナルティ付きGLMの推定：glmnetアルゴリズム

L1ペナルティを含む目的関数は $\beta_j = 0$ において微分不可能であるため、Newton-Raphson法やFisher Scoringを直接適用できません。glmnetパッケージが採用するのは座標降下法であり、他の係数を固定した状態で1つの係数 $\beta_j$ について目的関数を最小化する操作をすべての $j$ にわたって繰り返します。

L1ペナルティ下における $\beta_j$ に関する最適化条件はsubgradient条件から導かれ、軟閾値作用素による更新式が得られます。軟閾値作用素とは、入力値の絶対値がしきい値以下のとき出力をゼロとし、それ以外では入力値からしきい値分だけ原点方向に縮小する関数です：

$$
\hat{\beta}_j \leftarrow S\!\left(\tilde{\beta}_j,\, \lambda\right) = \operatorname{sign}(\tilde{\beta}_j)\max\!\left(|\tilde{\beta}_j| – \lambda,\, 0\right)
$$

ここで $\tilde{\beta}_j$ は他の係数を固定したときの最小二乗的更新量です。$|\tilde{\beta}_j| \leq \lambda$ であれば $\hat{\beta}_j = 0$ となり、これがスパース解の数学的根拠となります。

GLMへの拡張では、各反復内で現在の推定値 $\hat{\beta}^{(t)}$ を中心に対数尤度を2次近似することでワーキング応答 $z_i$ とワーキング重み $w_i$ を構成します。ロジスティック回帰（リンク関数：logit）の場合は次のようになります：

$$
w_i = \hat{\mu}_i(1 – \hat{\mu}_i), \quad z_i = \hat{\eta}_i + \frac{y_i – \hat{\mu}_i}{w_i}
$$

ポアソン回帰（リンク関数：log）の場合は次のようになります：

$$
w_i = \hat{\mu}_i, \quad z_i = \hat{\eta}_i + \frac{y_i – \hat{\mu}_i}{\hat{\mu}_i}
$$

この操作はIRLS（反復重み付き最小二乗法）の枠組みと整合しており、各反復内で加重Lasso問題を座標降下法で解くという二重ループ構造をとります。正則化パス全体を効率的に計算するためにwarm startingが用いられます。これは $\lambda$ の値を大きい側から小さい側へ順に変化させ、前の $\lambda$ での解を次の $\lambda$ の初期値として使用する方法であり、収束の高速化に寄与します。

前提条件：座標降下法の収束保証は、目的関数が各座標について凸であることを必要とします。GLMの対数尤度は正則指数族のもとで凸であり、L1・L2ペナルティも凸であるため、ペナルティ付きGLMの目的関数全体の凸性は保証されます。

限界：カテゴリ変数のダミー変数群のようにグループ構造を持つ変数に座標降下法を単純適用すると、グループ内の任意の1変数のみが選択されるという恣意的な結果をもたらす可能性があります。このような場合にはグループLassoやSparse Group Lassoへの拡張が適切です。

正則化パスの解釈：ロジスティック回帰とポアソン回帰への適用

正則化パスとは、$\lambda$ を連続的に変化させたときの係数推定値の軌跡であり、横軸に $\log(\lambda)$、縦軸に各係数値をプロットすることで可視化されます（Fig.1）。$\lambda$ が大きい（右側）ではすべての係数がゼロに縮小されており、$\lambda$ を小さくするにつれて係数が順次ゼロから離脱します。最後まで非ゼロを保つ変数ほど応答変数との関連が強いと言えます。

   
   

（Fig1. ロジスティック回帰における正則化パス（Lasso）：log(λ)の減少とともに係数がゼロから非ゼロへ活性化される様子を示しています。上部の数字は非ゼロ係数数を示しています。）

ロジスティック回帰でのLassoは、スパースな分類モデルの構築に有効です。$p \gg n$ の高次元設定では非ゼロ係数を持つ変数のみがモデルに残るため、解釈可能性と予測精度を両立しやすくなります。一方、ポアソン回帰でのRidgeは、カウントデータの共線性が強い場面での係数安定化に有用です。Ridgeは係数を縮小しますがゼロにはしないため、すべての変数をモデルに保ちつつ分散を抑制します。

最適な $\lambda$ の選択には交差検証が用いられます（Fig.2）。cv.glmnetでは各 $\lambda$ についてCV逸脱度の平均と標準誤差を計算し、2種類の候補を提示します。$\lambda_{\min}$ はCV誤差を最小化する値であり、$\lambda_{1\text{se}}$ はCV誤差が最小値の1標準誤差以内に収まる最大の $\lambda$ です。後者はより強い正則化（より少ない変数）を与えます。

   
   

（Fig2. 交差検証によるλ選択：cv.glmnetの出力イメージです。λ_minとλ_1seの位置を示す二本の破線と、各λでのCV逸脱度（平均±1SE）を示しています。）

CV評価指標として逸脱度ベースの誤差が標準的に使用されます。二項分布の場合：

$$
\text{CV-deviance} = -\frac{2}{n} \sum_{i \in \text{val}} \left[ y_i \log \hat{\mu}_i + (1-y_i)\log(1-\hat{\mu}_i) \right]
$$

限界：正則化パスは係数の取捨選択の参考情報を提供しますが、統計的仮説検定の代替にはなりません。どの変数がゼロに縮小されるかという順序は探索的指標です。また、最適化された係数には縮小バイアスが残存するため、選択後の変数についてバイアス補正なしに信頼区間を構成することは統計的に問題があります。変数選択とパラメータ推定を同一データで行うことで推定量の分布が変化し、ナイーブな信頼区間は過狭になります（選択後推論問題）。

高次元2値応答変数（p>>n）での変数選択

ゲノムワイド関連解析（GWAS）やがんバイオマーカー研究では、遺伝子発現量や一塩基多型などの特徴量数 $p$ が数百から数万に達するのに対し、患者数 $n$ は数十〜数百にとどまることが多いです。このような $p \gg n$ の設定でLasso logistic回帰は変数選択の実用的な手段となります。

Lassoの変数選択一致性——真のゼロ係数を確率1でゼロと推定し、真の非ゼロ係数を正確に検出する性質——が成立するためには、irrepresentable conditionと呼ばれる設計行列の条件が必要です。この条件は、真のゼロ変数が真の非ゼロ変数の線形結合で近似できないこと、すなわちゼロ変数と非ゼロ変数の間の相関が十分弱いことを要求します。

実務的な高次元解析では、単一の $\lambda$ での選択結果に過度に依存することは推奨されません。安定性選択は、ブートストラップ標本を繰り返し生成し、各標本でのLasso選択結果から各変数の選択確率を推定する手法です。選択確率が高い変数を安定した重要変数とみなすことで、単一実行の偶然性を排除できます。

さらに、Lassoによる変数選択の後にその選択された変数で通常のGLMを再フィットすることが係数の推定として推奨されます。Lasso係数は縮小バイアスを含むため、推定値の大きさ自体を解釈する際には注意が必要です。選択後推論の問題は、最尤推定の信頼区間の構成手法がLasso選択後にはそのまま適用できないことを示しています——変数選択とパラメータ推定を同一データで行うことで推定量の分布が変化し、ナイーブな信頼区間は過信になります。

前提条件：Lassoの変数選択一致性はスパースな真のモデル（真に重要な変数が $p$ に比べて少数）を仮定します。また、観測は互いに独立であることが要求されます。反復測定データやクラスター構造を持つデータ（例：複数施設からの患者データ）に対しては、ペナルティ付きGEEや正則化GLMMへの拡張が必要です。

限界：Lassoは強相関を持つ変数群からは1つを恣意的に選ぶ傾向があります——どの変数が選ばれるかは標本の微小な変動に依存し、再現性が低いです。また、選択された変数の係数は縮小推定量であるため、効果の大きさを推論するには選択後の再推定が望ましいです。

Ridge・Lasso・Elastic Net の比較と選択指針

3手法の性質を理論・実用両面から整理すると以下のようになります（Table 1）。

Ridgeの固有値収縮の観点からは、設計行列の特異値分解 $X = UDV^\top$ を用いると、線形モデルでのRidge推定量は $\hat{\beta}_{\text{Ridge}} = V(D^2 + \lambda I)^{-1}DU^\top y$ と表せます。分母に $\lambda$ が加わることで小さな固有値に対応する方向の係数が強く縮小され、情報行列の特異性が正則化されます。GLMではこの解析解は厳密には成立しませんが、2次近似の枠組みで類似の解釈が成立します。

Lasso解の非唯一性は、2変数が完全相関（$x_j = x_k$）の場合に生じます。このとき $\beta_j – \beta_k$ の任意の組み合わせが同じ残差を与えるため、解が一点に定まりません。Elastic NetはL2項がこの退化を防ぐため、強相関変数が存在してもより安定した解を返します。Elastic Netの $\alpha$ パラメータは二次元グリッド探索または交差検証によって $\lambda$ と同時に最適化されます。

非正則化GLMとの比較では、$p < n$ かつ共線性が軽微な設定では最尤推定量が最良不偏推定量としての性質を持ちます。正則化は不偏性を犠牲にして分散を削減するため、予測を目的とする場面や変数選択が必要な場面に適しています。ステップワイズ変数選択と比べると、正則化は連続的な縮小を通じて選択を行うため、段階的選択に伴う多重検定の問題を回避できる利点があります。

限界：3手法いずれも、予測変数が標準化されていない場合は係数のスケールに応じてペナルティの影響が変動します。実装上は予測変数を平均0・標準偏差1に標準化してからglmnetを適用することが標準的です。また、切片項はペナルティの対象外とするのがglmnetを含む主要実装の慣行です。

   
   

（Fig3. 同一データに対するRidge・Lasso・Elastic Netの係数推定値の比較（λは各手法でCV最適化）。Lassoのスパース性とRidgeの均等縮小を視覚的に対比しています。）

biostatistics応用例：高次元がんバイオマーカーデータでのLasso logistic回帰

ここでは正則化GLMの理論を具体的な分析フローに結びつけるため、高次元がんバイオマーカーデータへの適用を概説します。典型的な設定として、RNA-seqによる遺伝子発現量プロファイル（$p = 500 \sim 5000$ 遺伝子）を予測変数とし、腫瘍の悪性度や治療奏効の有無（二値応答）を $n = 50 \sim 200$ 症例で予測する問題を考えます。

分析の第一段階は前処理です。各遺伝子発現量を患者間で平均0・標準偏差1に標準化します。この操作を怠ると、発現量のスケールが大きい遺伝子が小さいペナルティで済むという不公平な縮小が生じます。発現量の分布が対数正規に近い場合は $\log_2$ 変換後に標準化することが多いです。

次に、cv.glmnetによる $\lambda$ の選択を行います。一般的な実装では10分割交差検証を用い、各フォールドで二項逸脱度を評価します。$\lambda_{\min}$ と $\lambda_{1\text{se}}$ の両方について選択変数のリストを確認し、安定して選択される変数を特定します。安定性選択を追加で実施する場合は、ブートストラップ標本（例：100回）ごとにLassoを適用し、各遺伝子の選択確率を算出します。

モデルの性能評価には、ホールドアウト検証データでのAUC（ROC曲線下面積）と逸脱度を用います。AUCは二値分類器の識別能力を要約する指標であり、閾値に依存しない比較が可能です。訓練データでのAUCは過楽観的になるため、報告には必ず検証データでの値を使用します。

結果の報告にあたっては、選択された変数数、使用した $\lambda$ の値（$\lambda_{\min}$ か $\lambda_{1\text{se}}$ かの明記を含む）、CV評価指標、および外部検証の有無を明記します。選択された変数の係数はLassoによる縮小推定量であるため、選択変数のみで再フィットした通常のロジスティック回帰の係数と比較することで、縮小バイアスの程度を確認することが望ましいです。

限界：内部交差検証のみによるAUC評価は、独立したコホートへの汎化性能を過大評価する傾向があります。これは $\lambda$ の選択自体がCVデータを用いるため、情報の漏洩が完全には排除されないためです。バイオマーカーとしての臨床的有用性の確認には、独立した外部コホートでの再現が不可欠です。また、発見された関連は予測的関連であり、個々の遺伝子との因果的関係を直接的に主張することはできません。