玄関雑学の部屋雑学コーナー統計学入門

13.4 プロビット分析

(1) D50

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

薬剤の安全性の指標としてはLD0(0%致死量)またはTD0(0%毒性用量)が理想であり、有効性の指標としてはED100(100%有効量)が理想です。 しかし現実にはLD0やED100を正確に推定するのは困難です。 そこで普通はLD50やED50を推定し、これらの用量を薬剤の安全性や有効性の目安にしています。

(2) 閾用量

D50については用量と反応率の関係を調べ、それに基づいて反応率が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群だけ他の群よりも多く割りつけてあるのは、このあたりの用量がD50と予想されるのでD50推定値の精度を高くするためです。

表13.4.1 名義尺度の時の用量反応試験のデータ
用量反応数非反応数 計 反応率(%)
0.01g/kg020200.0
0.1g/kg2182010.0
1g/kg16143053.3
10g/kg1552075.0
100g/kg1912095.0

表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字状曲線)になり、それは累積正規分布で近似できることがわかります。

図13.4.1 閾用量の分布と用量反応率関数

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

表13.4.2 プロビット分析の分散分析表
要因平方和自由度χ2
回帰(直線性)38.5139138.5139
ズレ(異質性)1.9412331.94123
全体(用量)40.55524 
○回帰
用量−プロビット直線: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のメタメター
寄与率:r2 = 0.966806(96.7%)
χβ2 = 38.5139(p = 5.4364×10-10) > χ2(1,0.05) = 3.841 … 有意水準5%で有意
○ズレ(LOF:Lack Of Fit)
χLOF2 = 1.94123(p = 0.584694) < χ2(3,0.05) = 7.815 … 有意水準5%で有意ではない
図13.4.2 用量反応曲線と用量反応直線

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

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

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

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

表13.4.3 コクラン・アーミテージの
傾向分析による用量反応解析の分散分析表
要因平方和自由度 χ2値 
回帰(直線性)52.1752152.1752
ズレ(異質性)1.7913931.79139
全体(用量)53.96664 
○回帰
用量反応直線:p = 0.472727 + 0.255x … 図13.4.2の青色の直線
寄与率:r2 = 0.966806(96.7%)
χβ2 = 52.1752(p = 5.07069×10-13) > χ2(1,0.05) = 3.841 … 有意水準5%で有意
○ズレ
χLOF2 = 1.79139(p = 0.61681) > χ2(3,0.05) = 7.815 … 有意水準5%で有意ではない

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

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

表13.4.4 ロジット分析の分散分析表
要因平方和自由度χ2
回帰(直線性)27.8304127.8304
ズレ(異質性)1.3766231.37662
全体(用量)29.2074 
○回帰
用量反応曲線(ロジスティック曲線):図13.4.2の緑色の曲線
p = 1 1 + exp(0.134839 - 1.51748x)
寄与率:r2 = 0.952867(95.3%)
χβ2 = 27.8304(p = 1.32431×10-7) > χ2(1,0.05) = 3.841 … 有意水準5%で有意
○ズレ
χLOF2 = 1.37662(p = 0.711025) > χ2(3,0.05) = 7.815 … 有意水準5%で有意ではない

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

b(ロジット分析) ≒ b(プロビット分析) × π √(3)
※例題の場合
b(ロジット分析) = 1.51748 ≒ 0.975747× π √(3) = 1.76981

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

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

(4) D50の推定

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

用量−プロビット直線:y = 4.85489 + 0.975747x
Φ-1(p) = Φ-1(0.5) = 0 = z より、p = 0.5 の時 z = 0 → y = 5
x 0 = y - 4.85489 0.975747 = 5 - 4.85495 0.975747 = 0.148721 ← 閾用量分布の平均値
閾用量分布の標準偏差: SD x = 1 0.975747 = 1.0249
D50 = 100.148721 = 1.40838g/kg
95%信頼区間 下限:D50L = 0.673536g/kg (xL = -0.171639)  上限:D50U = 3.00642g/kg (xU = 0.47805)

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

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

○コクラン・アーミテージの傾向検定による用量反応直線:p = 0.472727 + 0.255x
x = 0.5-0.472727 0.255 = 0.106952
D50 = 100.1069529 = 1.27924g/kg
95%信頼区間 下限:D50L = 0.542033g/kg (xL = -0.265975)  上限:D50U = 3.13963g/kg (xU = 0.496879)
○ロジット分析による用量反応曲線: p = 1 1 + exp(0.134839 - 1.57148x)
0.134839 - 1.51748x = ln(1 - 0.5) - ln(0.5) = 0 → x = 0.134839 1.51748 = 0.0888574
D50 = 100.0888574 = 1.22704g/kg
95%信頼区間 下限:D50L = 0.518285g/kg (xL = -0.285431)  上限:D50U = 2.76428g/kg (xU = 0.441582)

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

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

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

表13.4.5 実用量の
プロビット分析の分散分析表
要因平方和自由度χ2
回帰(直線性)15.1168115.1168
ズレ(異質性)28.9702328.9702
全体(用量)44.0874 
○回帰
用量−プロビット直線:y = 4.57215 + 0.0236626x
寄与率:r2 = 0.342885(34.3%)
χβ2 = 15.1168(p = 0.0001) > χ2(1,0.05) = 3.841 … 有意水準5%で有意
○ズレ(LOF:Lack Of Fit)
χLOF2 = 28.9702(p = 2.2719×10-6) > χ2(3,0.05) = 7.815 … 有意水準5%で有意
D50 = 18.0812  95%信頼区間 下限:D50L = 7.13539 上限:D50U = 36.1706
表13.4.6 実用量の傾向分析の分散分析表
要因平方和自由度χ2
回帰(直線性)26.1624126.1624
ズレ(異質性)27.8042327.8042
全体(用量)53.96664 
○回帰
用量−反応率直線:p = 0.341832 + 0.00645033x
寄与率:r2 = 0.484789(48.5%)
χβ2 = 26.1624(p = 3.13873×10-7) > χ2(1,0.05) = 3.841 … 有意水準5%で有意
○ズレ(LOF:Lack Of Fit)
χLOF2 = 27.8042(p = 3.99246×10-6) > χ2(3,0.05) = 7.815 … 有意水準5%で有意
D50 = 24.5208  95%信頼区間 下限:D50L = 9.47437 上限:D50U = 41.0226
表13.4.7 実用量の
ロジット分析の分散分析表
要因平方和自由度χ2
回帰(直線性)9.7771419.77714
ズレ(異質性)11.0948311.0948
全体(用量)20.8724 
○回帰
用量−反応率曲線(ロジスティック曲線):p = 1 1 + exp(0.0501408 - 0.0332455x)
寄与率:r2 = 0.468434(46.8%)
χβ2 = 9.77714(p = 0.00176695) > χ2(1,0.05) = 3.841 … 有意水準5%で有意
○ズレ(LOF:Lack Of Fit)
χLOF2 = 11.0948(p = 0.0112239) > χ2(3,0.05) = 7.815 … 有意水準5%で有意
D50 = 1.5082  95%信頼区間 下限:D50L = -20.707 上限:D50U = 16.1064

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


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

表13.4.8 名義尺度の用量反応試験の一般的データ
用量対数用量反応個体数非反応個体数総個体数反応率
D1x1r1n1-r1n1p1
::::::
Dixirini-rinipi
::::::
Daxarana-ranapa

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

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


zi:標準正規分布Φ(zi)の関数値がpiになる時の正規偏位
yi = zi + 5 = β0 + β1xi:反応率piのプロビット

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

全体の尤度:
対数尤度:

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

傾斜ベクトル:k = L(k)   :ナブラ(ハミルトン演算子)   k = [b0k b1k]'
(j = 0,1)
ヘスの行列:k = k
(j = 0,1 l = 0,1)
ここで と zi = b0 + b1x1 - 5 より
  

ただし、重み: とする。
  
  
k+1 = k-k-1k

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

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

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

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

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

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


  


← 期待プロビットkから求めた値のため
収束条件: < ε = 0.01〜0.0001

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

  
  
  
SLOF = Syy - Sβ = χLOF2   φy = a - 1  φβ = 1   φLOF = φyy - φβ = a - 2
寄与率:
表13.4.9 プロビット分析の分散分析表
要因平方和自由度χ2
回帰(直線性)Sβφβχβ2
ズレ(異質性)SLOFφLOFχLOF2
全体(用量)Syyφy 
○回帰(直線性)の検定
帰無仮説 H0:対数用量−プロビット直線の傾きが0である(β1 = 0)
χβ2 > χ2β,α)の時、有意水準100α%で有意
○ズレ(異質性)の検定
帰無仮説 H0:対数用量−プロビット直線と観測プロビットが一致している
χLOF2 > χ2LOF,α)の時、有意水準100α%で有意

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

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

補正反応率: … pi = p0 + pi*(1 - p0)と考える
補正重み:

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

○観測プロビット:Y1 = Φ-1(0) + 5:↓(無限小のため計算から除外)
Y2 = Φ-1(0.1) + 5 = 3.71845   Y3 = Φ-1(0.533333) + 5 = 5.08365
Y4 = Φ-1(0.75) + 5 = 5.67449   Y5 = Φ-1(0.95) + 5 = 6.64485
○対数用量と観測プロビットの回帰分析:n = 4  mx = 0.5   my = 5.28036
Sxx = 5  Sxy = 4.68503  b00 = 4.81186   b10 = 0.937005
○期待プロビット初期値:
  
  
○実用プロビット初期値:
  
  
○対数用量と実用プロビットの重み付け回帰分析:∑niwi = 43.7251   mx = 0.141585  my = 4.99299
Sxx = 43.7763  Syy = 43.432  Sxy = 42.641   b01 = 4.85508  b11 = 0.974064
○期待プロビット1回目:
  
  
○実用プロビット1回目:
  
  
○対数用量と実用プロビットの重み付け回帰分析:∑niwi = 42.6067   mx = 0.108156  my = 4.96048
Sxx = 40.5982  Syy = 40.5868   Sxy = 39.6123  b02 = 4.85495   b12 = 0.975717
○期待プロビット2回目:
  
  

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

  SLOF = 40.5868 - 38.6504 = 1.9364  φy = 4   φβ = 1  φLOF = 3

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

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

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

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

用量−プロビット直線:y = 4.85495 + 0.975717x
  D50 = 100.148657 = 1.480818 g/kg
t(∞,0.05) = 1.960  VR = 1  ∑niwi = 42.6067   mx = 0.108156  Sxx = 40.5982