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

前口上 | 目次 | 第1章 | 第2章 | 第3章 | 第4章 | 第5章 | 第6章 | 第7章 | 第8章 | 第9章 | 第10章
第11章 | 第12章 | 第13章 | 第14章 | 第15章 | 第16章 | 第17章 | 第18章 | 第19章 | 第20章 | 付録

1 | 2 | 3 | 4 | 5

## 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/kg | 0 | 20 | 20 | 0.0
0.1g/kg | 2 | 18 | 20 | 10.0
1g/kg | 16 | 14 | 30 | 53.3
10g/kg | 15 | 5 | 20 | 75.0
100g/kg | 19 | 1 | 20 | 95.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字状曲線)になり、それは累積正規分布で近似できることがわかります。

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

表13.4.2 プロビット分析の分散分析表
要因 | 平方和 | 自由度 | χ2値
--- | --- | --- | ---
回帰(直線性) | 38.6504 | 1 | 38.6504
ズレ(異質性) | 1.93636 | 3 | 1.93636
全体(用量) | 40.5868 | 4 |

○回帰
用量−プロビット直線：y=4.85495 + 0.975717x
用量反応曲線(プロビット曲線)：p=Φ(z)=Φ(y-5)=Φ(4.85795 + 0.975717x-5) … 図13.4.2の赤色の曲線
：標準正規分布
z=Φ-1(p)：Φ(z)の逆関数で、Φ(z)の関数値がpになる時の正規偏位
p：反応率　　y=z+5：pのプロビット　　x=log(D)：用量Dのメタメター
寄与率：r2=0.952291(95.2％)
χβ2=38.6504(p=5.06917×10-10)＞χ2(1,0.05)=3.841 … 有意水準5％で有意
○ズレ(LOF：Lack Of Fit)
χLOF2=1.93636(p=0.58572)＜χ2(3,0.05)=7.815 … 有意水準5％で有意ではない

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

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

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

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

表13.4.3 コクラン・アーミテージの
傾向分析による用量反応解析の分散分析表
要因 | 平方和 | 自由度 | χ2値
--- | --- | --- | ---
回帰(直線性) | 52.1752 | 1 | 52.1752
ズレ(異質性) | 1.79139 | 3 | 1.79139
全体(用量) | 53.9666 | 4 |

○回帰
用量反応直線：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値
--- | --- | --- | ---
回帰(直線性) | 26.5135 | 1 | 26.5135
ズレ(異質性) | 1.48434 | 3 | 1.48434
全体(用量) | 27.9979 | 4 |

○回帰
用量反応曲線(ロジスティック曲線)：図13.4.2の緑色の曲線

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

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

### (4) D50の推定

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

用量−プロビット直線：y=4.85495 + 0.975717x
Φ-1(p)=Φ-1(0.5)=0=zより、p=0.5の時z=0 → y=5

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

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

○コクラン・アーミテージの傾向検定による用量反応直線：p=0.472727 + 0.255x

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

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

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

表13.4.5 実用量の
プロビット分析の分散分析表
要因 | 平方和 | 自由度 | χ2値
--- | --- | --- | ---
回帰(直線性) | 15.0714 | 1 | 15.0714
ズレ(異質性) | 28.971 | 3 | 28.971
全体(用量) | 44.0424 | 4 |

○回帰
用量−プロビット直線：y=4.57215 + 0.0236603x
寄与率：r2=0.342202(34.2％)
χβ2=15.0714(p=0.0001)＞χ2(1,0.05)=3.841 … 有意水準5％で有意
○ズレ(LOF：Lack Of Fit)
χLOF2=28.971(p=2.71017×10-6)＞χ2(3,0.05)=7.815 … 有意水準5％で有意
D50=18.0829　　95％信頼区間　下限：D50L=7.13533　上限：D50U=36.2138

表13.4.6 実用量の傾向分析の分散分析表
要因 | 平方和 | 自由度 | χ2値
--- | --- | --- | ---
回帰(直線性) | 26.1624 | 1 | 26.1624
ズレ(異質性) | 27.8042 | 3 | 27.8042
全体(用量) | 53.9666 | 4 |

○回帰
用量−反応率直線：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

表13.4.7 実用量の
ロジット分析の分散分析表
要因 | 平方和 | 自由度 | χ2値
--- | --- | --- | ---
回帰(直線性) | 10.6076 | 1 | 10.6076
ズレ(異質性) | 17.3902 | 3 | 17.3902
全体(用量) | 27.9979 | 4 |

○回帰

表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の時に、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-1ｇk

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

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

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

期待プロビット初期値： 期待反応確率初期値：
実用プロビット初期値： 重み初期値：

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

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

これは重み付け最小2乗法に相当するため、対数用量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＞χ2(φLOF,α)の時、有意水準100α％で有意

なお用量を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.480818g/kg
t(∞,0.05)=1.96　　VR=1　　∑niwi=42.6067 　　mx=0.108156　　Sxx=40.5982