13.4 プロビット分析

(1) D₅₀

医学・薬学分野では、反応が計量尺度のデータではなく生／死や有効／無効といった名義尺度のデータでしか得られないことがよくあります。 この時、対象とする母集団の50％の個体が反応する用量を50％反応量といい、D₅₀とかD50などと書きます。 例えば反応が生／死の時は50％致死量または中央致死量(median lethal dose)といってLD₅₀と書き、反応が有効／無効の時は50％有効量といってED₅₀と書き、反応が毒性有／毒性無の時は50％毒性用量といってTD₅₀と書き、反応が阻害有／阻害無の時は50％阻害用量といってID₅₀と書きます。

薬剤の安全性の指標としてはLD₀(0％致死量)またはTD₀(0％毒性用量)が理想であり、有効性の指標としてはED₁₀₀(100％有効量)が理想です。 しかし現実にはLD₀やED₁₀₀を正確に推定するのは非常に困難です。 そこで普通はLD₅₀やED₅₀を推定し、これらの用量を薬剤の安全性や有効性の目安にしています。

(2) 閾用量

D₅₀については用量と反応率の関係を調べ、それに基づいて反応率が50％の時の用量を推定します。 そのためには用量反応直線における目的変数を、計量尺度の反応実測値の代わりに反応率にした手法を利用することが考えられます。 それは説明変数が計量尺度で目的変数が名義尺度の時の回帰分析である、コクラン・アーミテージの傾向分析に相当します。 (→5.3 計数値の相関分析と回帰分析)

しかし用量反応関係の場合、その関係を近似する関数を閾用量(いきようりょう、しきいようりょう、threshold dose)の分布に基づいて理論的に導くことができます。 閾用量とは、ある個体が反応を起こす最低の用量のことです。 例えば、ある人がお酒を呑んで酔いつぶれる最低の酒量は”酔いつぶれ反応”に関するその人の閾用量です。 この閾用量は、たった1合で酔いつぶれてしまう”酒に弱い人”もいれば、1升飲んでも大丈夫という”酒豪”もいるというように人によって様々です。

それと同様に特定の反応の閾用量も個体によって様々であり、閾用量はその反応特有の分布をしていると考えられます。 そして薬剤の閾用量分布は近似的に対数正規分布をすることが多く、用量を対数用量にすれば正規分布で近似できることが経験的にわかっています。 このことから用量反応率関数を理論的に導くことができます。

(3) 用量反応率関数とプロビット変換

具体的なデータに基づいて用量反応率関数を理論的に導いてみましょう。 反応が名義尺度の時の試験デザインは反応が計量尺度の時の試験デザインと同様です。 例えばマウスに、ある薬物について0.01g/kg、0.1g/kg、1g/kg、10g/kg、100g/kgの用量を無作為に割り付けて投与し、反応数を観測したところ表13.4.1のようになったとします。 なお1g/kg群だけ他の群よりも多く割りつけてあるのは、このあたりの用量がD₅₀と予想されるのでD₅₀推定値の精度を高くするためです。

表13.4.1において0.1g/kg群で反応した2匹は閾用量が0.1g/kg以下だったマウスであり、反応しなかった18匹は閾用量が0.1g/kgよりも多かったマウスのはずです。 つまり反応した2匹のマウスの閾用量は0g/kg〜0.1g/kgの範囲の色々な値のはずです。 同様に1g/kg群で反応した16匹のマウスの閾用量は0g/kg〜1g/kgの範囲の色々な値のはずです。 そうすると1g/kg群の反応数16というのは、1g/kgの薬剤を投与された30匹のうち閾用量が0g/kg〜1g/kgまでのマウスを合計した数ということになります。

このことから表13.4.1の5群のマウスが同じ特性を持ち、薬剤に対する反応性が同じなら、用量が多くなるにしたがって反応率は必ず増加または横ばいになり、減少することはないと予想されます。 そして閾用量が対数正規分布する時、対数用量と反応率の関係は正規分布を累積した関数つまり累積正規分布関数になるはずです。 このことから用量反応率関数はシグモイド曲線(S字状曲線)になり、それは累積正規分布で近似できることがわかります。

第10章の第2節で説明したように、累積正規分布曲線は反応率をプロビット変換することによって直線に変換できます。 この原理を利用して用量反応関係を累積正規分布曲線つまりプロビット曲線で近似し、D₅₀を逆推定する手法をプロビット分析(probit analysis)といいます。 表13.4.1のデータにその手法を適用すると次のようになります。 _(注1) (→10.2 各種のシグモイド曲線)

○回帰
用量−プロビット直線：y = 4.85489 + 0.975747x
用量反応曲線(プロビット曲線)：p = Φ(z) = Φ(y-5) = Φ(4.85489 + 0.975747x-5) … 図13.4.2の赤色の曲線
　　：標準正規分布
z=Φ^-1(p)：Φ(z)の逆関数で、Φ(z)の関数値がpになる時の正規偏位
p：反応率　　y=z+5：pのプロビット　　x=log(D)：用量Dのメタメター
寄与率：r² = 0.966806(96.7％)
χ_β² = 38.5139(p=5.4364×10^-10) ＞ χ²(1,0.05)=3.841 … 有意水準5％で有意
○ズレ(LOF：Lack Of Fit)
χ_LOF² = 1.94123(p=0.584694) ＜ χ²(3,0.05)=7.815 … 有意水準5％で有意ではない

表13.4.2の「全体(用量)」は5群の反応率の変動つまり用量の違いによる反応率の変動を表します。 「回帰(直線性)」は用量による変動の中で用量−プロビット直線つまり用量反応曲線で説明できる変動を表し、「ズレ(異質性)」は用量反応曲線では説明できない変動つまり用量反応曲線からのズレを表します。 そのため回帰の平方和とズレの平方和を合計したものが用量の平方和になり、回帰の自由度とズレの自由度を合計したものが用量の自由度になります。

また用量反応曲線の寄与率は回帰の平方和を用量の平方和で割った値であり、反応率の変動のうち用量反応曲線で説明できる変動の割合を表します。 データが出現率のため平方和がχ²値になることを除けば、これらは計量尺度のデータの用量反応解析と同様です。

このデータの場合、回帰の検定結果が有意でズレの検定結果は有意ではありません。 このことから反応率の変動には用量依存性があり、それはほぼ用量反応曲線で近似できる——曲線からのズレは5％程度——ということがわかります。 図13.4.2を見ると、確かに用量反応曲線は反応率のプロットとうまくフィットしていることがわかります。

ちなみに表13.4.1のデータにコクラン・アーミテージの傾向分析を適用すると次のようになります。 (→5.3 計数値の相関分析と回帰分析)

この場合も回帰の検定結果が有意でズレの検定結果は有意ではなく、用量反応直線の寄与率が約97％もあります。 したがってこのデータの場合、用量反応関係を直線で近似しても大差はないことがわかります。

またロジスティック回帰分析を用いた用量反応解析のことをロジット分析(logit analysis)といいます。 表13.4.1のデータにロジット分析を適用すると次のようになります。 (→10.3 ロジスティック回帰分析の計算方法)

この場合も検定結果はプロビット分析と同様であり、両者の用量反応曲線は非常によく似ています。 そして理論的にはプロビット分析の回帰係数とロジット分析の回帰係数の間には次のような関係があります。

しかしロジスティック曲線は質的に異なる2つの群があり、それぞれの群の説明変数が正規分布する時に説明変数と一方の群に属す確率の関係を表す曲線です。 したがって用量反応解析のように閾用量の分布に基づいて用量反応曲線を求める時はプロビット分析の方が正確です。

プロビット分析は計算が複雑なので、ロジット分析で近似的に代用する時があります。 しかしコンピュータが発達した現在ではどちらの手法も計算可能ですから、プロビット分析を用いた方が正確です。 計算が簡単という理由で代用するなら、むしろコクラン・アーミテージの傾向分析を用いた方が計算が簡単で結果の解釈も容易です。 なおこれら3種類の手法の検定は、たいてい有意性検定になります。 そのため検定結果よりも用量反応曲線(または直線)そのものや寄与率を科学的に検討する方が有意義です。

(4) D₅₀の推定

用量反応曲線を利用してD₅₀を推定する方法は、用量反応直線を利用して特定の反応の時の用量を逆推定する方法と原理的には同じです。 ただし用量反応曲線の式が複雑なので計算方法が少し複雑になります。 プロビット分析によってD₅₀とその95％信頼区間を求めると次のようになります。 _(注2)

D₅₀の用量メタメターx₀は閾用量分布の平均値になり、用量−プロビット直線の傾きの逆数は閾用量分布の標準偏差になります。 このことから用量−プロビット直線の傾きが大きければ閾用量分布の幅が狭くなり、用量の変化によって反応が敏感に変化することがわかります。 したがって閾用量分布の標準偏差つまり用量−プロビット直線の傾きの逆数は用量と反応の敏感性を表す指標として利用することができます。

ちなみにコクラン・アーミテージの傾向分析によって求めた用量反応直線と、ロジット分析によって求めた用量反応曲線を利用してD₅₀を推定すると次のようになります。

これらの計算結果から、どの方法でD₅₀を推定しても大きな違いはないことがわかると思います。

(5) 実用量を用いた場合

反応が計量尺度の時と同様に、反応が名義尺度の時も実用量を用いて解析することができます。 表13.4.1のデータに実用量を用いたプロビット分析と、コクラン・アーミテージの傾向分析と、ロジット分析を適用すると次のようになります。

表13.4.1のデータは対数用量用のデモデータですから、実用量を用いるとズレが大きくなっています。 しかし薬剤の有効率は用量と比例することが多いので、臨床試験などでは実用量を用いた方が良い結果になることも多いと思います。 そのため薬剤の用量反応関係を検討するための臨床試験では対数用量を用いた解析結果と実用量を用いた解析結果を比較し、医学的により妥当な方を選べば良いと思います。

(注1)　表13.4.1を一般化すると次のようになります。

表13.4.8 名義尺度の用量反応試験の一般的データ

用量	対数用量	反応個体数	非反応個体数	総個体数	反応率
D1	x1	r1	n1-r1	n1	p1
:	:	:	:	:	:
Di	xi	ri	ni-ri	ni	pi
:	:	:	:	:	:
Da	xa	ra	na-ra	na	pa

理論的な反応確率がπ_iの時に、n_i例中r_i例が反応する確率は二項分布から次のようになります。

プロビット分析では個体の閾用量が対数正規分布し、対数用量反応率関係が累積正規分布になると仮定します。 そのため理論反応確率π_iを次のように表すことができます。 この理論反応確率π_iをプロビット変換してプロビットy_iにすると、y_iと対数用量x_iの関係が直線になります。 (→10.2 各種のシグモイド曲線)

z_i：標準正規分布Φ(z_i)の関数値がp_iになる時の正規偏位
y_i=z_i + 5=β₀ + β₁x_i：反応率p_iのプロビット

ロジスティック回帰分析と同様に、この対数用量−プロビット直線のパラメーターβ₀とβ₁の推定値b₀とb₁を最尤法によって求めます。 (→10.3 ロジスティック回帰分析の計算方法 (注2))

全体の尤度：
対数尤度：

普通の最尤法では、次のように対数尤度を最大にする時のb₀とb₁をニュートン・ラプソン(Newton-Raphson)法によって近似的に求めます。

傾斜ベクトル：ｇ_k=∇L(ｂ_k) 　　∇：ナブラ(ハミルトン演算子) 　　ｂ_k=[b_0k b_1k]'
(j=0,1)
ヘスの行列：Ｈ_k=∇ｇ_k
(j=0,1 l=0,1)
ここでとz_i=b₀ + b₁x₁ - 5より
　　
∴
ただし、重み：とする。
　　
　　
ｂ_k+1=ｂ_k-Ｈ_k^-1ｇ_k

しかし次のような少し変則的なニュートン・ラプソン法を用いると、プロビットが正規偏位であることを利用して回帰とズレの検定を効率的に行うことができます。

観測プロビット(p_iをプロビット変換した値)：
理論プロビット(π_iをプロビット変換した値)：
z_iにおける確率密度(1正規偏位あたりの標準正規分布の積分値)：
∴

ここでb₀とb₁に適当な初期値b₀₀とb₁₀——例えば観測プロビットY_iとx_iの単純な回帰分析から求めた値——を代入して計算したプロビットを期待プロビットと呼び、期待プロビットを逆変換して求めた確率値を期待反応確率と呼ぶことにします。 これらの値を理論プロビットと理論反応率の代わりに用いて、観測プロビットを推測すると次のようになります。 これを実用プロビット(working probit)といいます。

期待プロビット初期値： 　　期待反応確率初期値：

実用プロビット初期値： 　　重み初期値：

この実用プロビットを用いてニュートン・ラプソン法を行うと、それは対数用量と実用プロビットの回帰直線y=β₀+β₁xについて、重み付け最小2乗法によって係数β₀とβ₁の推測値を求めることに相当します。

　　


← 期待プロビットはｂ_kから求めた値のため
収束条件：＜ε=0.01〜0.0001

収束後、対数用量x_iと収束後の実用プロビットy_iを用いて回帰とズレの分散分析を行うことができます。

　　
　　
　　
S_LOF=S_yy - S_β=χ_LOF² 　　φ_y=a - 1　　φ_β=1 　　φ_LOF=φ_yy - φ_β=a - 2
寄与率：

表13.4.9 プロビット分析の分散分析表

要因	平方和	自由度	χ2値
回帰(直線性)	Sβ	φβ	χβ2
ズレ(異質性)	SLOF	φLOF	χLOF2
全体(用量)	Syy	φy

○回帰(直線性)の検定
帰無仮説 H₀：対数用量−プロビット直線の傾きが0である(β₁=0)
χ_β²＞χ²(φ_β,α)の時、有意水準100α％で有意
○ズレ(異質性)の検定
帰無仮説 H₀：対数用量−プロビット直線と観測プロビットが一致している
χ_LOF²＞χ²(φ_LOF,α)の時、有意水準100α％で有意

この実用プロビットを用いた繰り返し計算法は反応確率が0または1のデータにも適用できます。 そこでこの方法を重み付け最小2乗法によるロジット分析に応用すると、出現率が0または1のデータがあっても、それらを除外せずにパラメータを求めることができます。 つまり出現率が0または1のデータがある時は、それらのロジットを期待値を用いて補完し、繰り返し計算によってパラメータを求めるわけです。 これを繰り返し補完法といます。 詳しい説明は「10.3 ロジスティック回帰分析の計算方法　(注1)」をご覧ください。 (→10.3 ロジスティック回帰分析の計算方法　(注1))

なお用量を0にしたコントロール群で反応があり、その反応率p₀が0よりも大きい時は次のようなアボット(Abbott)の式によって補正します。 この補正はp₀≦0.2なら適切に補正できると言われています。

補正反応率： … p_i=p₀ + p_i^*(1 - p₀)と考える
補正重み：

実用量を用いたプロビット分析はxとして用量そのものを用いるだけで、計算式は変わりません。 表13.4.1のデータについて実際に計算してみましょう。 実用量を用いたプロビット分析に興味のある人は自分で計算してみてください。

○観測プロビット：Y₁=Φ^-1(0) + 5：↓(無限小のため計算から除外)
Y₂=Φ^-1(0.1) + 5=3.71845 　　Y₃=Φ^-1(0.533333) + 5=5.08365
Y₄=Φ^-1(0.75) + 5=5.67449 　　Y₅=Φ^-1(0.95) + 5=6.64485
○対数用量と観測プロビットの回帰分析：n=4　　m_x=0.5 　　m_y=5.28036
S_xx=5　　S_xy=4.68503　　b₀₀=4.81186 　　b₁₀=0.937005
○期待プロビット初期値：
　　
　　
○実用プロビット初期値：
　　
　　
○対数用量と実用プロビットの重み付け回帰分析：∑n_iw_i=43.7251 　　m_x=0.141585　　m_y=4.99299
S_xx=43.7763　　S_yy=43.432　　S_xy=42.641 　　b₀₁=4.85508　　b₁₁=0.974064
○期待プロビット1回目：
　　
　　
○実用プロビット1回目：
　　
　　
○対数用量と実用プロビットの重み付け回帰分析：∑n_iw_i=42.6067 　　m_x=0.108156　　m_y=4.96048
S_xx=40.5982　　S_yy=40.5868 　　S_xy=39.6123　　b₀₂=4.85495 　　b₁₂=0.975717
○期待プロビット2回目：
　　
　　

1回目と2回目の期待プロビットの差の最大値が0.01よりも小さいため、ここで計算を終了します。 そして次のような統計量を求めると、表13.4.2の分散分析表を作成することができます。

　　S_LOF=40.5868 - 38.6504=1.9364　　φ_y=4 　　φ_β=1　　φ_LOF=3

(注2)　プロビット分析によってD₅₀を推定する方法は、用量反応直線を利用して特定の反応の時の用量を逆推定する方法とほぼ同様です。 そしてD₅₀の100(1-α)％信頼区間も、同様にフィーラーの式を利用して求めます。 (→13.1 用量反応直線 (注2))

用量−プロビット直線：y=b₀ + b₁x
y=5の時のxをx₀とすると
※閾用量分布の平均値はx₀、標準偏差は1/b₁になる。
x=log(D)の時：D₀=10^x0
x₀の100(1-α)％信頼区間：
ただし
※g＜0.1の時はg=0とした次のような近似式も用いられる。

t=t(∞,α)：正規分布の100α％点　　V_R=1
※ズレが大きい時はt=t(φ_LOF,α)、とすることもある

表13.4.1のデータについてD₅₀の95％信頼区間を実際に計算してみましょう。 実用量を用いた時はx₀の推定値をそのままD50の推定値にします。

用量−プロビット直線：y=4.85495 + 0.975717x
　　D₅₀=10^0.148657=1.480818g/kg
t(∞,0.05)=1.96　　V_R=1　　∑n_iw_i=42.6067 　　m_x=0.108156　　S_xx=40.5982

最終更新日：2016年3月10日 　第3節へ 第5節へ

　webmaster@snap-tck.com 　Copyleft (C) 2000 SNAP(Sugimoto Norio Art Production)