物質の毒性を調べる際などに，その指標としてLD50(50% Lethal Dose, 半数致死量)という値をしばしば求める*1．これはその濃度の物質を個体群に投与した場合に，半数の個体が死亡する濃度や量とされる．
素朴に考えると，濃度を説明変数，死亡率を従属変数として直線回帰を行い，死亡率=50%となる点の濃度をLD50とすれば良い様に思える．
ともあれ，まずはデータを用意する．今回データはRでGLMをやってみようから二硫化炭素ガスに5時間暴露されたカブトムシの死亡個体数データを用いる．
dose <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839) # ガス濃度[log_10 CS_2 mgl^-1] affected <- c(6, 13, 18, 28, 52, 53, 61, 60) # 死亡数 total <- c(59, 60, 62, 56, 63, 59, 62, 60) # 供試個体数
プロットするとこんな感じ．参考として死亡率10%, 50%, 90%にラインを引いた．
後のためこの図を描くための関数を作っておく．
baseplot <- function(){ plot(dose, affected/total, cex = total/max(total), ylim=c(0, 1), ylab="Mortality rate", xlab="Gas concentration" ) ## 50%ライン abline(0.5, 0, lty = 2, col="skyblue") ## 10%ライン abline(0.1, 0, lty = 2, col="pink") ## 90%ライン abline(0.9, 0, lty = 2, col="pink") }
一応以下ではLD10，LD50，LD90を求めることを目標とする．
直線回帰によるLD50
死亡率(死亡数/供試個体数)とガス濃度で単純にlm()してみる．また，predict()を使って観測値の95%信頼区間も求める．
baseplot() ### 直線回帰 fit.lm <- lm((affected/total)~dose) abline(fit.lm) coef(fit.lm) ## ED10, 50, 90 ed.lm <- function(r, p = c(ed10 = 0.1, ed50 = 0.5, ed90 = 0.9)){ (p - coef(r)[1]) / coef(r)[2] } ed.lm(fit.lm) ## 信頼区間 x <- data.frame(dose = seq(min(dose), max(dose), diff(range(dose))/100)) con <- predict(fit.lm, x, interval="confidence") lines(x$dose, con[,2], lty = 4) lines(x$dose, con[,3], lty = 4)
得られるLD10, 50, 90の値はこのようになる．
> ed.lm(fit.lm) ed10 ed50 ed90 1.699146 1.774264 1.849382
プロット結果を示す．
一見問題なさそうにも見えるのだが，予測値と信頼区間が0を下回ったり1を上回ったりしていて少し気持ち悪い．そもそもデータは死亡or生存の二値データなので誤差は二項分布しないとおかしい．
そこでこうしたデータを扱う場合は普通の線形モデルではなくて一般化線形モデル(Generalized Linear Model, GLM)というものを使ってパラメータの推定やらをする，というのが流行りらしい．
のだけれども，今ひとつ理解してないので詳しいことはググるか図書館で調べて下さい．
データをそのままロジスティック変換とかするのではなくて，データの期待値を変換して，観測値の誤差はまた別に考えるというものではないかと思ってはいるのですが，なにぶんちょっとググっただけのアレなのであまりアテにしないで下さい．下でやってることは合ってると思います．僕が作ったのではなくてもとからあるglm()を使ってるので大丈夫に決まってます．
logit法によるLD50
上のプロットから考えるに，死亡率とガス濃度の関係を説明するためには半分くらいの死亡率の時には傾きが急で，死亡率が1や0に近くなるほど傾きが緩やかになり漸近していくようなS字の曲線が良さそうに思える．
そのような性質を持つ関数として，ロジスティック関数と正規分布の累積関数が挙げられる．説明のための曲線にロジスティック関数を使うものはlogit法，正規分布累積関数を使うものはprobit法と呼ばれる．logitとはロジスティック関数の逆関数，probitとは正規分布累積関数の逆関数のことを指す．
Rではglm関数を使って次のようにするだけ．
### ロジット法 fit.logit <- glm(cbind(affected, total-affected) ~ dose, family=binomial(logit)) # 誤差は二項分布，リンク関数はロジット
glmの場合は最小二乗法でもってバシっとパラメータが決まる，ということはなくて，何らかの繰り返しアルゴリズムでもって尤度を最大化する必要があるのだけれどもその辺もやってくれるらしい．
lm()と目立って違うのはfamilyという引数の違いだろう．binomialというのは誤差分布が二項分布に従うということの指定で，(logit)というのは曲線がロジスティック関数に従うということの指定．(logit)はfamily=binomial指定時のデフォルトなので省略してもいい．
結果はlmオブジェクトと同様にsummary()で得られる．
> summary(fit.logit) Call: glm(formula = cbind(affected, total - affected) ~ dose, family = binomial(logit)) Deviance Residuals: Min 1Q Median 3Q Max -1.5941 -0.3944 0.8329 1.2592 1.5940 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -60.717 5.181 -11.72 <2e-16 *** dose 34.270 2.912 11.77 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 284.202 on 7 degrees of freedom Residual deviance: 11.232 on 6 degrees of freedom AIC: 41.43 Number of Fisher Scoring iterations: 4
さらにここから予測値，つまり期待値をつなげた曲線を得るにはpredict関数を使う．
baseplot() ## 予測値の計算 x <- data.frame(dose = seq(min(dose), max(dose), diff(range(dose))/100)) p.pred <- predict(fit.logit, x, type="response") lines(x$dose, p.pred)
また，誤差は二項分布に従うと考えているのだから，95%信頼区間はqbinom()を使って求めることができる．
## (doseに対するrateの)予測値信頼区間(95%) mean.N <- round(length(total)/sum(1/(total))) p.lw <- qbinom(0.025, mean.N, p.pred)/mean.N p.up <- qbinom(0.975, mean.N, p.pred)/mean.N lines(smooth.spline(x$dose, p.lw), lty = 2) lines(smooth.spline(x$dose, p.up), lty = 2)
さらに，LD50等を簡単に求めるためのdose.p()という関数がMASSパッケージに用意されている．ついでにSEまで求められる．至れり尽くせりとはこのこと．
> dose.p(fit.logit, p = c(0.1, 0.5, 0.9)) Dose SE p = 0.1: 1.707606 0.007126466 p = 0.5: 1.771721 0.003858052 p = 0.9: 1.835835 0.006192525
ただ，ここで求めたSEを1.96倍しても先程の95%信頼区間とは重ならないことに注意しておく必要がある．先程の信頼区間は死亡率の信頼区間で，例えば「濃度がこれこれの時に死亡率は95%の確率でここからここまでに入っていますよ」というようなことを言うときに使う値．一方，dose.p()により算出された標準誤差，またそれを1.96倍して得た95%信頼区間というのは「この死亡率を出すための濃度は95%の確率でここからここまでの間ですよ」というようなことを言うときに使う値(多分)．
実際の問題を考えると後者の方が有用なのだろう．他のLD50を求めるためのパッケージや関数(例えばdoByパッケージのdose.LD50()関数)で表示される信頼区間は後者の値となる．
そこで濃度の95%信頼区間曲線も求めてみよう．dose.p()関数を利用すれば比較的簡単に求められる(SEの値が妙な入れ方してあって少し戸惑ったが)．
## (rateに対するdoseの)予測値信頼区間(95%) d <- dose.p(fit.logit, p = (p <- seq(0.01, 0.99, 0.01))) dse <- attributes(d)$SE p2.lw <- d - dse*1.96 p2.up <- d + dse*1.96 lines(p2.lw, p, lty = 3) lines(p2.up, p, lty = 3)
ちなみにこちらの信頼区間の方が狭くなる．理由は知らない．プロット結果を示す．
さらに先程の直線回帰の結果と重ねて描いてみる．
probit法によるLD50
probit法だと正規分布を元にしているということもあって，logit法より自然な方法ととらえられるらしい．bioassayなんかの説明で「本当はprobit法を使うんだけど対応してる統計ソフトもあまりないしlogit法でも結果はほとんど一緒だからlogitでもいいよ」というようなことが言われる場合もあるらしい(僕が聞いたり読んだりしたわけじゃないですが)．
まあでもRなら簡単にできる．glm()のfamily=の指定でfamily=binomial(probit)と指定するだけ*2．あとはlogit法の場合と一緒．
なお，既に作成したモデルの一部だけ変えたいときはupdate()を使うことができる．
### プロビット法 fit.probit <- update(fit.logit, family = binomial(probit))
あとは一緒なのでざっと流して
baseplot() ## ED10, 50, 90 dose.p(fit.probit, p = c(0.1, 0.5, 0.9)) ## 予測値の計算 p.pred.pro <- predict(fit.probit, x, type="response") lines(x$dose, p.pred.pro) ## (rateに対するdoseの)予測値信頼区間(95%) mean.N <- round(length(total)/sum(1/(total))) p.lw.pro <- qbinom(0.025, mean.N, p.pred.pro)/mean.N p.up.pro <- qbinom(0.975, mean.N, p.pred.pro)/mean.N lines(smooth.spline(x$dose, p.lw.pro), lty = 2) lines(smooth.spline(x$dose, p.up.pro), lty = 2) ## (doseに対するrateの)予測値信頼区間(95%) d <- dose.p(fit.probit, p = p) dse.pro <- attributes(d)$SE p2.lw.pro <- d - dse.pro*1.96 p2.up.pro <- d + dse.pro*1.96 lines(p2.lw.pro, p, lty = 3) lines(p2.up.pro, p, lty = 3)
LD10, 50, 90は次の通り．
> MASS::dose.p(fit.probit, p = c(0.1, 0.5, 0.9)) Dose SE p = 0.1: 1.705891 0.006708782 p = 0.5: 1.770852 0.003803312 p = 0.9: 1.835814 0.005646886
logit法の場合と重ねて見ると分かるけど実際あまり変わらない．
なんか特徴のあるデータならともかく，今回の例のようなデータならprobitにするかlogitにするかは好みの問題だろうか．
一応線形回帰も合わせて重ねたやつを
あと，LD50を求める関数しかないけど，doByというパッケージを使うとdose.LD50()という関数でLD50と信頼区間を求められる．結局glmオブジェクトが必要なのだけど，helpを見る感じだと説明変数が複数の場合にも使えるらしい．
## LD50と信頼区間 library(doBy) dose.LD50(fit.logit, c(1, NA)) dose.LD50(fit.probit, c(1, NA))
出力は以下の通り．
> dose.LD50(fit.logit, c(1, NA)) ld50.numvec lower upper 1.771721 1.764159 1.779283 > dose.LD50(fit.probit, c(1, NA)) ld50.numvec lower upper 1.770852 1.763398 1.778307
