非線形回帰の基礎概念
線形回帰では応答変数を説明変数の線形結合で表現しますが、非線形回帰では応答変数をパラメータ$\boldsymbol{\theta}$に対して非線形な関数$f$で記述します。NLS(非線形最小二乗法)の出発点となる一般モデル形式は次のように定義されます。
$$y = f(\mathbf{x}; \boldsymbol{\theta}) + \varepsilon$$
ここで$\mathbf{x}$は説明変数ベクトル、$\boldsymbol{\theta} = (\theta_1, \ldots, \theta_p)^\top$は推定対象のパラメータベクトル、$\varepsilon$は誤差項です。「非線形」とは、$y$への説明変数ではなく、パラメータ$\boldsymbol{\theta}$に対して関数$f$が線形でないことを意味します。たとえば$f(x; \theta_1, \theta_2) = \theta_1 e^{\theta_2 x}$はパラメータに対して非線形であり、このモデルは説明変数$x$の線形変換だけでは表現できません。
一部の非線形モデルは変数変換によって線形化できます。対数変換$\log y = \log\theta_1 + \theta_2 x$がその例であり、このような変換可能な非線形性は変換後にOLS(最小二乗法)が適用可能です。一方、$f(x; \theta_1, \theta_2) = \theta_1 / (1 + e^{-\theta_2 x})$のような真の非線形性は変数変換だけでは線形化できず、解析的な閉形式解が存在しないため数値最適化が必要となります。非線形モデルが必要となる代表的な場面として、成長曲線(生物集団の増殖)、用量反応曲線(薬理学的応答)、指数減衰(薬物動態)が挙げられます。
誤差$\varepsilon_i$には独立性・均一分散性・正規性の仮定が課されます。これらはOLSと共通の出発点ですが、非線形回帰固有の制約として初期値依存性と局所最小値への収束リスクが存在します。$S(\boldsymbol{\theta})$が一般に非凸であるため、解析解が得られず初期値の設定が推定結果に本質的な影響を与えます。
非線形最小二乗法の数理的枠組み
NLSは残差平方和
$$S(\boldsymbol{\theta}) = \sum_{i=1}^{n} [y_i – f(\mathbf{x}_i; \boldsymbol{\theta})]^2$$
を最小化するパラメータ$\hat{\boldsymbol{\theta}}$を求める手続きです。目的関数$S(\boldsymbol{\theta})$の勾配ベクトルは
$$\nabla_{\boldsymbol{\theta}} S(\boldsymbol{\theta}) = -2\mathbf{J}^\top \mathbf{r}$$
と表されます。ここで$\mathbf{r} = \mathbf{y} – \mathbf{f}(\boldsymbol{\theta})$は残差ベクトルです。最小値における一階条件は$\mathbf{J}^\top\mathbf{r} = \mathbf{0}$であり、この方程式系を反復的に解くことがNLSアルゴリズムの本質です。
現在の推定値$\boldsymbol{\theta}^{(k)}$の近傍で$f(\mathbf{x}_i; \boldsymbol{\theta})$を1次テイラー展開すると局所線形近似
$$f(\mathbf{x}_i; \boldsymbol{\theta}) \approx f(\mathbf{x}_i; \boldsymbol{\theta}^{(k)}) + \mathbf{J}_i^\top \Delta\boldsymbol{\theta}$$
が得られます。ヤコビアン行列$\mathbf{J} \in \mathbb{R}^{n \times p}$は$(i, j)$要素が$\partial f(\mathbf{x}_i; \boldsymbol{\theta}^{(k)}) / \partial\theta_j$で定義される偏微分行列であり、$\Delta\boldsymbol{\theta} = \boldsymbol{\theta} – \boldsymbol{\theta}^{(k)}$です。この局所線形化のもとで$S(\boldsymbol{\theta})$の最小化問題は線形最小二乗問題に帰着されます。ヘッシアン$\nabla^2_{\boldsymbol{\theta}} S$の$2\mathbf{J}^\top\mathbf{J}$による近似(ガウス近似)は二階微分の積項を省略したものであり、この近似がOLS正規方程式$(\mathbf{X}^\top\mathbf{X})\hat{\boldsymbol{\beta}} = \mathbf{X}^\top\mathbf{y}$と形式的に対応しています。
この枠組みが有効であるためには$f(\mathbf{x}; \boldsymbol{\theta})$が$\boldsymbol{\theta}$に関して偏微分可能であること、および真のパラメータ近傍での局所線形近似が十分な精度を持つことが前提となります。高度な非線形性がある場合には線形近似の精度が低下し、初期値が真値から遠い場合にはアルゴリズムが発散するリスクがあります。
ガウス=ニュートン法
局所線形化されたOLS問題の正規方程式を解くと、ガウス=ニュートン更新式
$$\Delta\boldsymbol{\theta} = (\mathbf{J}^\top\mathbf{J})^{-1}\mathbf{J}^\top\mathbf{r}$$
が得られます。パラメータは次の反復更新則に従って改善されます。
$$\boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} + \Delta\boldsymbol{\theta}^{(k)}$$
$\mathbf{J}$が列フルランクであれば$\mathbf{J}^\top\mathbf{J}$は正定値となり更新方向が意味を持ちます。この更新式はOLS正規方程式と形式的に対応しており、NLSを局所的なOLS問題の繰り返し解法として解釈できます。本アルゴリズムにはステップサイズのダンピング機構がなく、各反復で算出した更新量をそのまま適用します。初期値が真値に十分近く残差が小さい条件下では二次収束が期待できますが、そのような好条件が常に保証されるわけではありません。
$\mathbf{J}^\top\mathbf{J}$が正則であること(ランク落ちしていないこと)と真値近傍での残差が十分小さいこと(小残差仮定)がアルゴリズムの安定動作の条件です。$\mathbf{J}^\top\mathbf{J}$が特異または近特異な場合には深刻な数値不安定性が生じ、大残差モデルでは収束が保証されません。収束速度は初期値とパラメータ空間の景観に依存します。
Levenberg-Marquardt法
ガウス=ニュートン法の数値不安定性を解消するために、ダンピングパラメータ$\lambda \geq 0$を導入したLevenberg-Marquardt法(LM法)では更新式が次のように修正されます。
$$\Delta\boldsymbol{\theta} = (\mathbf{J}^\top\mathbf{J} + \lambda\mathbf{I})^{-1}\mathbf{J}^\top\mathbf{r}$$
$\lambda > 0$のとき$\mathbf{J}^\top\mathbf{J} + \lambda\mathbf{I}$は常に正定値となり逆行列が必ず存在します。$\lambda \to 0$の極限ではガウス=ニュートン更新に一致し、$\lambda \to \infty$の極限では更新量が$\mathbf{J}^\top\mathbf{r}$方向の最急降下ステップに収束するという補間的性質があります。このため残差が大きく初期値が遠い状況では最急降下的に大域探索を行い、真値に近づくにつれてガウス=ニュートン的な高速収束に移行します。
$\lambda$は残差平方和の改善状況に応じて適応的に更新されます。残差が改善した反復では$\lambda$を縮小してガウス=ニュートン的な収束を促し、残差が悪化した反復では$\lambda$を拡大して探索を安定化します。この戦略は信頼領域法と等価であり、現在の推定値からの変化量を局所的な信頼領域内に制約する解釈が成立します。
$\lambda$の初期値と更新ルールは経験的に設定される点がLM法の制約です。また$\lambda\mathbf{I}$による等方的なレギュラライズはパラメータのスケールに依存するため、$\mathbf{J}^\top\mathbf{J} + \lambda\,\mathrm{diag}(\mathbf{J}^\top\mathbf{J})$による対角スケーリング変形がよりスケール不変な代替として用いられることがあります。
(Fig2. ガウス=ニュートン法とLM法の残差平方和の収束過程比較)
| 特性 | ガウス=ニュートン法 | Levenberg-Marquardt法 |
|---|---|---|
| 更新式の形式 | $\Delta\boldsymbol{\theta} = (\mathbf{J}^\top\mathbf{J})^{-1}\mathbf{J}^\top\mathbf{r}$ | $\Delta\boldsymbol{\theta} = (\mathbf{J}^\top\mathbf{J} + \lambda\mathbf{I})^{-1}\mathbf{J}^\top\mathbf{r}$ |
| 収束速度 | 真値近傍で二次収束(小残差条件下) | 超線形収束($\lambda$適応後) |
| 数値安定性 | $\mathbf{J}^\top\mathbf{J}$が近特異な場合に不安定 | $\lambda > 0$により常に安定 |
| 初期値への感度 | 高い(真値近傍が前提) | 比較的低い(大域的探索性あり) |
| ダンピング機構 | なし | あり(適応的$\lambda$更新) |
| 推奨する適用場面 | 良条件・小残差のモデル | 一般的な非線形回帰 |
収束診断とパラメータの不確実性評価
NLS推定の停止判定には3種の収束基準が用いられます。第1は勾配ノルム$\|\mathbf{J}^\top\mathbf{r}\| < \epsilon_g$(一階条件への近さ)、第2はパラメータ変化量$\|\Delta\boldsymbol{\theta}^{(k)}\| < \epsilon_\theta$(更新量の収束)、第3は残差変化量$|S(\boldsymbol{\theta}^{(k+1)}) - S(\boldsymbol{\theta}^{(k)})| < \epsilon_S$(目的関数の安定化)です。これら3基準のいずれかまたは複数の閾値への到達で収束を判定します。最大反復回数を超えても収束しない場合は初期値の変更または$\lambda$初期値の調整が検討されます。
NLS推定量$\hat{\boldsymbol{\theta}}$の漸近的分散共分散行列は
$$\widehat{\mathrm{Var}}(\hat{\boldsymbol{\theta}}) \approx s^2 (\mathbf{J}^\top\mathbf{J})^{-1}$$
で近似されます。残差標準偏差$s$は
$$s = \sqrt{\frac{S(\hat{\boldsymbol{\theta}})}{n – p}}$$
として定義され、$p$はパラメータ数です。この近似から得られる漸近信頼区間は$n \gg p$かつ誤差の正規性が成立する条件下で有効です。
初期値依存性の診断のために、複数の異なる初期値から最適化を実行して推定値の一致性を確認することが推奨されます。$\widehat{\mathrm{Var}}(\hat{\boldsymbol{\theta}})$を対角標準化して得られるパラメータ相関行列において高相関(例:$|r| > 0.99$)が観測される場合は、パラメータの識別性に問題がある可能性を示します。小標本($n < 20$程度)では漸近信頼区間がパラメータの不確実性を過小評価しやすく、ブートストラップ信頼区間による補完が必要です。さらにプロファイル尤度信頼区間はパラメータ空間の非対称性を直接反映できる点で漸近近似より優れており、精度が要求される場面での使用が推奨されます。
生体統計学への応用:用量反応モデルと成長曲線
薬理学的用量反応解析では4パラメータロジスティックモデルが標準的に用いられます。4PLモデルは
$$f(x; a, b, e, d) = d + \frac{a – d}{1 + (x/e)^b}$$
と定義されます。$a$は下漸近線(最小応答)、$d$は上漸近線(最大応答)、$e$はEC50(半数効果濃度、応答が$(a+d)/2$に達する濃度)、$b$はHill係数(濃度反応の急峻さ)です。EC50は薬効の効力指標として直接解釈可能であり、Hill係数の絶対値が大きいほど濃度に対する応答が急峻であることを定量的に示します。
ロジスティック成長曲線は個体群動態や細菌増殖の記述に用いられ、
$$f(t; K, r, t_0) = \frac{K}{1 + \exp(-r(t – t_0))}$$
と表されます。$K$は上限値(環境収容力)、$r$は成長率、$t_0$は変曲点(成長速度が最大となる時点)です。測定誤差が応答スケールに比例する場合には対数変換$\log(y)$を適用することで誤差の均一分散性を改善できます。2PLモデル(上下漸近線を固定)と4PLモデルの選択はAIC(赤池情報量規準)による数理的比較または先験的な生物学的情報に基づいて判断されます。
測定誤差の均一分散性は対数変換の適用根拠となる仮定であり、モデル形式の正しい特定(先験的情報に基づく選択)もパラメータ推定の有効性の前提です。外挿(上限・下限漸近値域外への予測)は曲線形状の仮定に強く依存するため避けるべきです。生物学的複製数が少ない実験($n < 20$)ではパラメータ推定精度が低く信頼区間の信頼性が下がります。また採用したモデル形式が真の生物学的機構と異なる場合(モデル誤特定)には系統的バイアスが生じ、EC50等の推定量の解釈が困難となります。
(Fig3. 4PLロジスティック用量反応曲線とEC50推定値)
(Fig1. ロジスティック成長曲線のNLS当てはめと95%信頼区間)
OLS・GLM・GAMとの比較と発展的手法
OLSは$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$に対して解析解$\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$を持ちますが、NLSには一般に解析解が存在せず数値反復が必要です。多項式回帰は入力変数の多項式変換によりOLSの枠内で曲線形状を表現できますが、パラメータに対して線形であるためパラメータ非線形性のある機構モデルとは区別されます。
GLM(一般化線形モデル)は指数型分布族の応答変数とリンク関数によるOLSの拡張であり、ポアソン回帰やロジスティック回帰が含まれます。応答変数が指数族に属し、かつ系統成分がパラメータに対してリンク関数を通じて線形に保たれる場合はGLMが適用可能です。一方、用量反応モデルのように応答がパラメータに対して真に非線形な構造(EC50やHill係数など実質的意味を持つパラメータ)を持つ場合にはNLSが必要です。
GAM(一般化加法モデル)は各説明変数に対してスプライン等の滑らかな非線形関数を推定する手法であり、非線形性の形状を事前に特定しない場合に有効です。NLSは機構的なモデル形式を前提としてパラメータが生物学的・物理的意味を持つ反面、GAMはデータ駆動的な柔軟な当てはめが可能ですがパラメータの実質的解釈が困難です。さらなる発展的手法として、ガウス過程回帰(カーネルによるベイズ的非線形回帰)やニューラルネットワーク回帰(大規模データでの表現学習)への接続が可能ですが、いずれもNLSと比較してパラメータ解釈可能性は低下します。
誤差分布仮定の違いも重要な選択基準です。OLSとNLSは正規誤差を前提とし、GLMは指数族の分散構造を、GAMは指定した分布族の分散構造をそれぞれ仮定します。離散応答変数(二値・カウント)に対してNLSを適用することは誤差構造の不整合を招くため不適切であり、その場合にはGLMまたはそれに準じた手法が選択されます。パラメータ解釈可能性と予測柔軟性は一般にトレードオフの関係にあり、NLSは機構モデルの枠内での解釈可能性を優先する場面で位置づけられます。

