致死量LD-50を算出する

　LD-50とは,例えば10匹のゴキブリに殺虫剤をかけたときに50%のゴキブリ（つまり5匹のゴキブリ）が死に至る量のことをいいます。似たような言葉でED-50というのもありますが,これは10人の被験者にある睡眠薬を飲ませた場合に50%の人（つまり5人）が睡眠薬の効果を受けて眠りにつく量のことをいいます。これは特に心理学の分野というわけではありませんが,心理学でも例えば毒性行動学と呼ばれる分野で扱われることでもあります。

プロビット分析

　チェスター・ブリス(1899～1979)によって開発された致死量を推定するための統計分析方法のひとつ。現在ではLD-50（あるいはED-50）の値を推定するためのコンピュータ・プログラムを記述することによって容易に推定値を算出することができる。
　プロビット分析やそれに変わるモデル（例えばロジット分析など）を用いて,25%あるいは80%で死亡する分量を表すLD-25やLD-85といった致死量を推定することも可能である。しかし50%から離れた値を推定しようとするほど,よりよい推定値を得るためには大規模な実験が必要になる（例えば多くのマウスを必要とする）。
　なおプロビット分析は現在でも統計学でよく用いられている手法の1つである重回帰分析と似た手法であるが,従属変数が（例えば反応あり＝1,反応なし＝0などのように）2値をとる場合には重回帰分析では負の値や1よりも大きな値が出てきてしまうので,そのような場合に適用すべき手法である。

　下表のようなデータが得られた場合について実際にプロビット分析をしてみましょう。

投与量（mg/kg）	反応数（匹）	例数（匹）	p
0.30	1	10	0/10=0.0
1.50	2	10	1/10=0.1
4.00	3	10	3/10=0.3
10.50	6	10	6/10=0.6
28.00	9	10	10/10=1.0

　とはいっても,私も詳しい計算手順は知らないので群馬大学の青木先生が書いてくださったRプログラムを拝借しましょう。これはed50という関数名になっていますが,結果はLD-50であっても同じことです。

書式：ed50(投与量,例数,反応数,tr="ln")

> touyo<-c(0.30,1.50,4.00,10.50,28.00) #投与量
> hanou<-c(1,2,3,6,9) #反応数
> rei<-c(rep(10,5)) #例数
> rei
[1] 10 10 10 10 10
> ed50(touyo,rei,hanou,tr="ln")
P_hat = 4.00764 + 0.562399 * ln(x)
ED50 = 5.83873
95% confidence limit = [ 2.70472, 14.708 ]

   Dose       n       r          e        r/n        e/n
 -1.204      10       1    0.47512    0.10000    0.04751
0.40547      10       2     2.2234    0.20000    0.22234
 1.3863      10       3     4.1578    0.30000    0.41578
 2.3514      10       6     6.2932    0.60000    0.62932
 3.3322      10       9     8.1102    0.90000    0.81102

Chi square value = 1.74288, d.f. = 3, P value = 0.627

Observed log likelihood = -24.3444
Fitted   log likelihood = -25.1984

　この結果によれば,

LD-50（致死量）＝5.83873
95% Confidence Limit（95%信頼限界）＝{2.70472,14.708}

となります。つまり致死量は5.84mgで,この値が2.70から14.71の範囲に入る確率が95%であるということです。ということは今回使ったある毒薬を6mgほど投与すれば50%のマウスが死んでしまうわけですね。

　"Dose"とは薬物を「投与する」という意味で"ドウス"と発音します（Doの3人称単数現在形のDoesではないのでご注意を）。これは対数変換された値となっています。

> touyo
[1]  0.3  1.5  4.0 10.5 28.0
> log(touyo)
[1] -1.2039728  0.4054651  1.3862944  2.3513753  3.3322045

　横軸に投与量の対数,縦軸にプロビット値として投与量-反応曲線を描くと直線で近似させることができます(下図参照)。なお,プロビット値というのはp=(反応数)／(例数)と定義しているときに,pに対応する正規分布のパーセント点に5を足したものです。上の例ではp=1/10=0.1なので,正規分布においてp=0.1に対応するパーセント点は-1.281552となります。この値に5を足すと-1.281552+5=3.718448となり,下表のプロビット値になるわけです(この作業をプロビット変換といいます)。

> p<-hanou/rei
> p
[1] 0.1 0.2 0.3 0.6 0.9
> per<-sapply(p,qnorm)
> per
[1] -1.2815516 -0.8416212 -0.5244005  0.2533471  1.2815516
> per+5
[1] 3.718448 4.158379 4.475599 5.253347 6.281552

対数変換後の投与量	パーセント点	プロビット値
-1.2039728	-1.2815516	3.718448
0.4054651	-0.8416212	4.158379
1.3862944	-0.5244005	4.475599
2.3513753	0.2533471	5.253347
3.3322045	1.2815516	6.281552

　以上のデータを基にグラフを作成します（2回目のplotでx軸のメモリを細かく取れば直線になります）。

> probit<-c(per+5)
> probit
[1] 3.718448 4.158379 4.475599 5.253347 6.281552
> log.data<-log(touyo)
> log.data
[1] -1.2039728  0.4054651  1.3862944  2.3513753  3.3322045
> y=4.00764+0.562399*log(hanou)
> y
[1] 4.007640 4.397465 4.625498 5.015324 5.243357
> plot(log.data,probit,ylim=c(3.0,7.0))
> par(new=T)
> plot(log.data,y,ylim=c(3.0,7.0),ylab="",xlab="",type="l",col="red")