データフィッティングを用いた用量反応関係のモデリング

この記事は nAG 最適化コーナー シリーズの一部です。

原文：Modelling Dose-Response Relationships Using Data Fitting

---

nAG ライブラリの最適化モデリングスイートに含まれるデータフィッティングソルバーを用いた用量反応関係の数学的モデリングの調査を行いました。異なる出発点を調査した結果、出発点の違いがモデルの質に影響を与えることが分かりました。また、用量反応解析において、汎用性の高いデータフィッティングソルバー、つまり、損失関数と正則化関数を簡単に交換できるソルバーは、データに最適なモデルを特定するのに非常に役立つことが明らかになりました。また、この調査では、重回帰を実行するときに発生する共線性の問題についても検討し、正則化を用いて共線性を減らす方法も示します。

用量反応解析の背景

用量反応モデルは、投与された刺激に対する生物の反応を記述するために使用されます。化学物質を異なる濃度で生物に投与した後、その反応を測定します。人間は生涯にわたって様々な化学物質に常にさらされているため、用量反応モデルの一つの用途は、化学物質の毒性を低コストで決定することです（化学物質に対する生物の反応をモデル化することは、例えば、ネズミによる従来の生物学的試験を行うよりも格段に安価で迅速です [5] ）。用量反応解析のさらなる応用として、組織培養における腫瘍の増殖を減少させる薬剤の効果を調べることが挙げられます。また、用量反応関係の調査は、公共の政策において、薬剤や汚染物質の用量の、安全、危険、（また場合によって）有益、を決定する上で重要な役割を果たします。

用量反応解析の歴史は、紀元前８世紀にさかのぼります。その時代から多くの人々がこの関係を解釈しようとしてきました [9] 。しかしながら、毒物学の重要人物であるパラケルススが、化学物質が有毒かどうかを決定する上で最も重要な要因は、その投与量であることに気づき、彼が「The dose makes the poison（毒となるかは用量次第）」という格言を作ったのは 1500 年代のことでした [11] 。そして、この分野が純粋に統計的アプローチに移行し、用量反応モデルをデータに適合させること（この記事で扱うこと）が一般的になったのは 1970 年代となります。

nAG ライブラリによるデータフィッティング

用量反応データにモデルを適合させるプロセスは、非線形回帰によって行うことができます。この回帰により、未知の値の補間が可能となり、異なる実験の結果を比較することができます。しかしながら、このタイプの実験データには外れ値が含まれることがよくあります。これは非線形モデルにおけるパラメーター推定に深刻な影響を与える可能性があります。そのため、回帰がロバストであることが重要です。回帰のロバスト性は、異なる損失関数を使用することによって変更することができます。nAG ライブラリのソルバールーチン handle_solve_nldf は、様々な損失関数の使用が可能です。より具体的には、handle_solve_nldf は以下の形式のデータフィッティング問題を解くことができます。

$$\begin {aligned} &\min _{x \in \mathbb {R}^{n_{\mathrm {var}}}} \quad &&\sum ^{n_{\mathrm {res}}}_{i=1} \mathcal {X}(r_i(x)) + \rho \sum ^{n_{\text {var}}}_{i=1} \psi (x_i) \\ &\text {subject to} \quad &&l_g \leq g(x) \leq u_g, \\ & &&\tfrac {1}{2} x^T Q_i x + p_i^T x + s_i \leq 0, \quad 1 \leq i \leq m_Q, \\ & &&l_B \leq Bx \leq u_B, \\ & &&l_x \leq x \leq u_x. \end {aligned}$$

ここで、$\mathcal {X}$ は損失関数で、各データポイントの観測値と予測値の差である残差 $r_i(x)$, $i=1,\dots ,n_{\mathrm {res}}$ に作用します。さらに、$\psi$ は正則化関数で、正則化係数 $\rho$ を持ちます。使用できる正則化関数には、$l_1$-norm と $l_2$-norm の２つがあります。handle_solve_nldf では、様々な種類の制約を最小化問題に課すことができます。これには、非線形制約、二次制約、線形制約、ボックス制約があります（詳細については、handle_solve_nldf の ドキュメント をご参照ください）。先に述べたように、様々な損失関数を使用することができ、これらの異なる損失関数は回帰のロバスト性に影響します。例えば、Huber 損失関数は $l_2$-norm 損失関数（最小二乗）よりもロバストです。残差が大きくなると、Huber 損失関数や他のよりロバストな損失関数に比べて、 $l_2$-norm 損失はより高いペナルティを持ちます（図１を参照）。

4 パラメーターロジスティックモデル

米国国家毒性プログラム（NTP）の 1408 種類の化学物質ライブラリから、１つの化学物質の用量反応データを使用します [8] 。このライブラリは、従来の毒性試験（動物を用いたもの）から離れ、化学物質が細胞反応や毒性経路にどのような影響を与えるかを評価する、予測モデルの開発に重点を置く取り組みから生まれました [6] 。このライブラリには、溶剤、防腐剤、汚染物質など、多くの環境化合物の毒性プロファイルが含まれています [12] 。試験には 14 の用量があり、１用量につき３つの反復があります（図２参照）。用量反応曲線はシグモイド状で単調な傾向があるため、Hill [3] がヘモグロビンへの酸素の結合を記述するために用いた式から派生した、以下の 4 パラメーターロジスティックモデル [1] が一般的にデータに適合されます：

$$y = a + \frac {b-a}{1 + (\frac {c}{x})^d}.$$

この非線形回帰モデルでは、$x$ は化学物質の用量、$y$ は対応する反応を表し、４つのパラメーターを推定します。$a$ は最小漸近値、$b$ は最大漸近値、$c$ は曲線の変曲点であり、集団の半数が研究された反応を示すと予想されるときの薬剤の濃度（しばしば $\mathrm {EC}_{50}$ または $\mathrm {IC}_{50}$ と呼ばれる）を示します。そして、$d$ は曲線の傾き（Hill 勾配とも呼ばれる）です。

このモデルは、リガンド結合 [10]、定量的薬理学 [2]、生命科学の力学系 [7] などの分野で見られる、多くの非線形関係を表現できる Hill 方程式ファミリーの１つです。

各種回帰モデルによる結果

Python 用の nAG ライブラリを用いて、用量反応モデルを NTP データにフィットさせました。使用したソースコードは、弊社の Github Local Optimization ページでご覧いただけます。まず、このモデルを「良い」出発点（この場合、最適解に近い出発点を意味します）からフィットさせた場合を考えてみましょう。図３を見ると、すべての損失関数が、期待されるシグモイド形状に沿ったモデル曲線を示していることが分かります。しかし、最適化で使用する損失関数によって、これらの曲線にはばらつきがあります。例えば、Huber 損失による曲線は（Cauchy、quantile、$l_1$-norm、smooth $l_1$-norm と同様）、用量濃度が大きい場合、各用量の中央点を正確に通り、他の２つの反復を外れ値として扱います。しかし、$l_2$-norm 損失はロバスト性が低いため、このような外れ値の影響を受けやすく、曲線がやや上向きになります。さらに、$l_\infty$-norm 損失と arctan 損失による曲線は、最大用量に対して異なる値を外れ値として扱うため、大きな用量に対して異なるモデルフィットを持つことになります。

非線形回帰は局所的な手法であるため（局所的な最適値を見つけるだけで、必ずしも大域的な最適値を見つけるわけではない）、最終的な解は与えられた初期値の質に影響されることがあります。そこで次に、より悪い出発点（これを「素朴な」と呼ぶことにします）からフィッティングしたモデルを考えてみましょう。なお、ここでは話を簡単にするために、Huber と $l_2$-norm の２つの損失関数に限定して分析を行います。

正則化付き Huber 損失関数

図４では、Huber 損失関数を用いて、素朴な出発点からスタートし、正則化関数を変えて（None、$l_1$-norm、$l_2$-norm）モデルをフィッティングした結果が示されています。まず、正則化関数を使わず（None）、素朴な出発点でフィッティングしたモデルを考えてみましょう。図４で見られるように、この曲線は、良い出発点を用いてフィットした曲線よりも角ばっています。このモデルのパラメーターの標準誤差を分析したい場合、Hill 方程式のヤコビ行列（$J$）からなる行列 $J^TJ$ を計算し、この行列の逆行列を求める必要があります。ただし、この最適解では、この行列の逆行列を求めることができないため、モデルに共線性があることを示唆しています。これは他の分析を行ったときには起きなかったことなので、おそらくデータからではなく、モデル自体に起因する構造的な多重共線性です。

具体的に言えば、回帰モデルにおいて、変数間に相関がある場合、多重共線性が生じることがあります。これはいくつかの問題を引き起こす可能性があります：パラメーター推定値がモデルの小さな変化に敏感になり、パラメーター推定値の精度、および、回帰の検定力も低下し、さらに、独立変数の統計的有意性が損なわれる可能性があります。

回帰モデルにおいて、多重共線性を減らす一つの方法は、正則化を使用することです。正則化は、パラメーター推定値をゼロに向かって縮小し、オーバーフィッティングを避けるために、目的関数にペナルティを追加します。そこで、このモデルの多重共線性に対処するために、$l_2$-norm 正則化関数を用いてモデルを再フィットすることにしました [4] 。その結果は、見るからに印象の悪いモデルフィットになりました（図４）。しかし、$l_1$-norm 正則化関数を使用すると、結果の曲線はより滑らかになり、良い出発点から見つかる最適解ほどではありませんが、素朴な出発点から見つかる最良の解となります。さらに、この正則化により、行列 $J^TJ$ の逆行列を求めることができ、パラメーターの標準誤差を求めることができます。

正則化付き $l_2$-norm 損失関数

次に、素朴な出発点に対して、$l_2$-norm 損失関数を用い、様々な正則化関数を用いてフィッティングしたモデルを検討します。図５を見ると、$l_1$-norm 正則化と正則化なし（None）の場合の曲線は、良い出発点から求めた曲線と非常に似ているように見えます。しかし、$l_2$-norm 正則化を使用した場合の曲線は、著しくフィットが悪いように見えます。表１では、$l_2$-norm 損失関数を用い、２つの異なる出発点からモデルをフィットさせた場合の最適解を比較しています。良い出発点を使ってフィッティングしたモデルと、素朴な出発点と $l_2$-norm 正則化を使ってフィッティングしたモデルに焦点を絞ると、パラメーター $d$ の標準誤差は、良い出発点よりも、素朴な出発点と $l_2$-norm 正則化の方が、著しく大きいことが分かります。このパラメーターは Hill 勾配を表しており、その傾きは２つのモデルで顕著に異なっています。一般的に、素朴な出発点と $l_2$-norm 正則化を用いたモデルを除く、他のすべてのモデルのパラメーターの標準誤差は、より規則的であるように見えます。

表１：$l_2$-norm 損失関数、２つの異なる出発点（"1" は良い出発点、"2" は素朴な出発点）、異なる正則化関数を用いた、各回帰パラメーターの推定値（Estimate）、標準誤差（SE）、95% 信頼区間（CI）

Parameter	Start	Regularisation	Estimate	SE	95% CI
$a$	1	None	$-38.2$	$4.4$	$(-47.1, -29.3)$
	2	None	$-3.4$	$1.0$	$(-5.4, -1.3)$
	2	$l_1$-norm	$-3.3$	$1.0$	$(-5.4, -1.3)$
	2	$l_2$-norm	$-3.0$	$1.1$	$(-5.4, -0.7)$
$b$	1	None	$-3.4$	$1.0$	$(-5.4, -1.3)$
	2	None	$-38.2$	$4.4$	$(-47.1, -29.3)$
	2	$l_1$-norm	$-37.7$	$4.3$	$(-46.4, -29.0)$
	2	$l_2$-norm	$-28.0$	$3.4$	$(-34.8, -21.2)$
$c$	1	None	$32.5$	$5.6$	$(21.3, 43.8)$
	2	None	$32.5$	$5.6$	$(21.3, 43.8)$
	2	$l_1$-norm	$31.8$	$5.5$	$(20.7, 42.9)$
	2	$l_2$-norm	$21.3$	$3.6$	$(14.1, 28.5)$
$d$	1	None	$-3.4$	$1.3$	$(-5.9, -0.8)$
	2	None	$3.4$	$1.3$	$(0.8, 5.9)$
	2	$l_1$-norm	$3.5$	$1.3$	$(0.8, 6.1)$
	2	$l_2$-norm	$5.7$	$15.8$	$(-26.4, 37.7)$

次に、素朴な出発点からフィッティングした $l_2$-norm 損失モデルのパラメーター推定値（表１参照）を、各パラメーターの元の定義と比較してみましょう。フィッティングするモデルは減少型なので、Hill 勾配であるパラメーター $d$ は負になることが期待されます。しかし、素朴な出発点から得られる最適解ではそうなっていません。さらに、最小漸近値を表すパラメーター $a$ と最大漸近値を表すパラメーター $b$ は、素朴な出発点から得られる最適解では意味が入れ替わっています。従って、回帰は基本的に減少するデータに増加する曲線を当てはめようとしており、パラメーターの推定値が本来の定義と一致した、良い出発点から始めた場合と比較して、用量反応曲線がわずかに異なっている理由を説明することができます。これは、良質な初期推定値（良い出発点）が、最終的な回帰モデルに影響を与えることを示しています。

結論

全般的に、用量反応解析には複数の要素を慎重に考慮する必要があります。非線形回帰モデルをフィットする際には、最適解に大きな影響を与える出発点を考慮することが重要です。また、最適解は損失関数の選択によっても影響を受けます。使用するデータに最も適したフィットを見つけるため、いくつかの損失関数（異なるロバスト性を持つ）を試すことが重要です。

参考文献

[1]   Gadagkar, S. R., & Call, G. B. $($2015$)$. Computational tools for fitting the Hill equation to dose-response curves. Journal of Pharmacological and Toxicological Methods, 71, 68-76.

[2]   Gesztelyi, R., Zsuga, J., Kemeny-Beke, A., Varga, B., Juhasz, B., & Tosaki, A. $($2012$)$. The Hill equation and the origin of quantitative pharmacology. Archive for History of Exact Sciences, 66$($4$)$, 427-438.

[3]   Hill, A. V. $($1910$)$. The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. Journal of Physiology, 40, 4-7.

[4]   Kesavulu, P., Vasu K., Bhupathi Naidu, M., Abbaiah, R., & Balasiddamuni, P. $($2016$)$. The effect of multicollinearity in nonlinear regression models. International Journal of Applied Research, 2, 506-509.

[5]   Lim, C., Sen, P. K., & Peddada, S. D. $($2013$)$. Robust analysis of high throughput screening $($HTS$)$ assay data. Technometrics, 55$($2$)$, 150-160.

[6]   National Toxicology Program. $($2004$)$. A national toxicology program for the 21st century: A roadmap for the future. National Toxicology Program: Research Triangle Park, NC, USA.

[7]   Pérez-Correa, J. R., Lefranc, G., & Fernández-Fernández, M. $($2015$)$. A new application of the Hill repressor function: Automatic control of a conic tank level and local stability analysis. Mathematical Problems in Engineering, 2015, Article ID 271216.

[8]   Tice, R. R., Fostel, J., Smith, C. S., Witt, K., Freedman, J. H., Portier, C. J., ...& Bucher, J. R. $($2007$)$. The National Toxicology Program high throughput screening initiative: Current status and future directions. Toxicologist, 46$($246$)$, 151.

[9]   Tsatsakis, A. M., Vassilopoulou, L., Kovatsi, L., Tsitsimpikou, C., Karamanou, M., Leon, G., ...& Spandidos, D. A. $($2018$)$. The dose response principle from philosophy to modern toxicology: The impact of ancient philosophy and medicine in modern toxicology science. Toxicology Reports, 5, 1107-1113.

[10]   Weiss, J. N. $($1997$)$. The Hill equation revisited: Uses and misuses. FASEB Journal, 11$($11$)$, 835-841.

[11]   Waddell, W. J. $($2010$)$. History of dose response. Journal of Toxicological Sciences, 35$($1$)$, 1-8.

[12]   Xia, M., Huang, R., Witt, K. L., Southall, N., Fostel, J., Cho, M. H., ...& Austin, C. P. $($2008$)$. Compound cytotoxicity profiling using quantitative high-throughput screening. Environmental Health Perspectives, 116$($3$)$, 284-291.