Rで致死試験の用量応答曲線を描く
用量応答曲線(dose-response curve)を描く必要に迫られた、Rで描けって言われたけど、Rを触ったことほぼないし…という学部生・M1向けの記事です。
This post is for undergraduate students who need to analyze dose-response relationships.
「そもそもRって何?」という方は、私が作った資料で恐縮ですが↓のスライドを見てみてください。この記事で紹介している教科書もおススメです。
Rで用量応答関係の解析をしているウェブサイトは、日本語でもいくつか見つかります(こちらとかこちら)。ただ、致死毒性試験の結果を解析した例が見当たらなかったので、この記事を作成することにしました。
記事作成にはRitz (2010, Environ Toxicol Chem 29: 220-229) とRitz et al. (2015, Plos One)、および上記の日本語で書かれたdrcについてのサイトを参考にしました。
結構長いです。とにかく図が描ければ良い、という忙しい方は「5. drcパッケージでlog-logistic式に当てはめる」と「6. drcパッケージで図を描く」だけ読んでいただければOKです。
1. 使用するデータ
下の図のような実験をおこなった場合を考えましょう。10匹のヨコエビをビーカーに入れ、96時間後の致死数をカウントします。1つの濃度区につきビーカーを3つ用意し、濃度区はControl・1 µg/L・2 µg/L・5 µg/L・10 µg/L・20 µg/Lの6つです。
下の表が実験結果だとします。表のDeadは致死個体の数、Totalはビーカーに入れた総数、Concは曝露濃度をそれぞれ示しています。対照区Controlは便宜的に0と書いておきます。
2. Rにデータを読み込ませる
Rに上のデータを取り込む方法はいくつかあります。
例えば Excelに上の表を作成し、該当部分をクリップポードに貼り付けた状態で次の関数を打ち込めばOKです。ここではtoxdataという名前でデータを格納しておきます。
toxdata <- read.delim ( "clipboard", header=TRUE ) ## headerは表の1行目を列の名前として認識するかどうかを示す
あるいはread.table関数でも同様です。
toxdata <- read.table ("clipboard", header=TRUE )
またはcsvファイルで上の表を保存してから、そのファイルを読み込むこともできます。ただこのread.csv関数を使う場合は、Rの作業ディレクトリとファイルのディレクトリ(保存場所)に注意してください。Rのディレクトリにファイルがなければ、csvファイルのディレクトリを指定する必要があります。
toxdata <- read.csv ( "toxdata.csv", header=TRUE )
toxdata <- read.csv ( "C:/home/user/toxdata.csv" ) #csvファイルの場所を指定
あるいは、Rにデータを直接打ち込んでもOKです。
toxdata <- data.frame ( Dead=c(1,0,0,2,0,1,2,1,1,3,2,4,7,5,9,10,9,10), Total =rep(10,18), Conc=c(rep(0, 3),rep(1,3),rep(2,3),rep(5,3),rep(10,3), rep(20,3)) )
3. とりあえず図を描いてみる
Rに取り込んだデータを図にしてみます。
plot( Dead/Total~Conc, data=toxdata )
Rでは目的変数y、説明変数xの関係をチルダ"~"を使って"y~x"のように表します。xは単一の変数でなく、"y~x1+x2"のようなモデル式でも構いません。
他に指定せず単純に"Dead/Total ~ Conc"と書くと、シンプルな線形回帰式 "Dead/Total = α × Conc + β" を想定していることになります。
その単純な式を図にすると下のようになります。濃度が20 μg/Lを超えると致死率が1 (= 100%) をオーバーするというモデルです。
さすがにこれは不自然なので、ある濃度以上だと致死率は1で一定となり (上限があり) 、ある濃度以下なら致死率はほぼ0である (下限がある) ことを表現できるモデル式が必要です。
4. log-logistic式
ではどんなモデル式にあてはめましょう?毒性学分野で広く使われている式にlog-logistic式があります。log-logistic式は、
で表現できます。xはの濃度で、b,c,d,eは回帰係数です。
これらの係数を変化させてlog-logistic式を描いたものが下図です。log-logistic式はS字型(シグモイド)の曲線を表現できることが分かりますね。また図を見ると、bは傾き、cは下限値、dは上限値、eはf(x)が上限値の半分となる際のxの値(毒性学的にいうとED50, Effective Dose50)に対応していることが読み取れます。
ちなみに致死率を目的変数とする場合、致死率の下限は0です。なので、係数c = 0として下のように係数3つのlog-logistic式を書くこともできます。
同様に致死率の上限は1ですから、さらに係数d=1として係数2つの式にすることもできます。
都市工の学部生にもおそらく馴染み深い話をすると、係数3つのlog-logistic式はb=-1の時に酵素の反応速度論で用いられるミカエリス・メンテン式と等しくなります。式変形をしてみます。(bの値が間違っているとのご指摘を受けて修正しました。2018.02.28)
ちなみにミカエリスメンテン式は↓。Vmaxは基質濃度xが最大の時の反応速度で、Kmはv=Vmax/2の時の基質濃度ですね。
▲Vmax=1・Km=1.5の時のミカエリス・メンテン式プロット
このような速度論とモデル式の話に興味がある人は、ヒルの式とかでググってください(適当)。
5. drcパッケージでlog-logistic式に当てはめる
ようやく本題です。
Rにはlog-ogistic式の係数を推定してくれて、描画まで楽ちんにしてくれるdrcというパッケージがあります。それを使ってみます。
始めてdrcを使う際は、まずパッケージをインストールしてみましょう。CRANのミラーサイトを選べと言われたら、とりあえずJapanのどこかを選んでください。
install.packages ("drc")
次にパッケージを読み込みます。
library(drc)
そして、「2. Rにデータを読み込ませる」で読み込んだtoxdataをlog-logistic式に当てはめます。drcパッケージでは、drm関数でモデルへのフィッティングができます。
model <- drm (Dead/Total ~ Conc, data=toxdata, type="binomial", weights=Total, fct = LL.2() )
Dead/Total ~ Concは上で書いた通り、"致死率を濃度で説明したい"ことを示しています。
次のtypeはどのような分布のデータを対象にしているかを示しています。致死試験の場合は、「死ぬ/生きている」の0/1で結果が出る二項分布のデータなので、"binomial"を指定します。デフォルトは正規分布になっていますが、他にもポアソン分布などが選べます。致死試験の場合はここで二項分布を指定して欲しい、ということがこの記事で書きたかった最大のポイントです。このあたりの議論の詳細は、「データ解析のための統計モデリング入門」などの教科書を参照してください。
そしてweights =Totalはビーカー当たりのヨコエビの個体数です。二項分布的に言うと試行回数ですね。
最後のfct=LL.2()は係数2つのlog-logistic式に当てはめる、という指示です。係数を3つにしたければLL.3、係数4つはLL.4にします。今回は、致死率データという性質上、係数2つのlog-logisticでOKでしょう。log-logistic式以外にも、対数をとらないlogistic (例えばL.4) やWeibull (W1.4やw2.4) などが選べます。
では、log-logistic式に当てはめた結果を見てみましょう。
summary (model)
Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 and upper limit at 1 (2 parms)
Parameter estimates:
Estimate Std. Error t-value p-value
b:(Intercept) -1.74314 0.26685 -6.53216 0
e:(Intercept) 5.96089 0.74685 7.98136 0
係数bとeの推定値が出力されます。drm関数は最尤推定法で係数を求めていると思います。詳しくはRitz et al., 2015 Plos Oneへ。eはED50 (ここではLC50) でしたが、LC50が5.9だとすると、「3. とりあえず図を描いてみる」の図となんとなく一致してますね。
LC50の信頼区間も計算できます。下の結果だと、95%信頼区間は4.50~7.42ですね。
ED(model,c(50), interval="delta", level= 0.95 ) #50をxにすればLCxを求められる
Estimated effective doses
Estimate Std. Error Lower Upper
e:1:50 5.96089 0.74685 4.49709 7.42469
6. drcパッケージで図を描く
また描画します。
plot (model)
これではダサいし見づらいので、色々付け加えていきます。
まず、線を太く、字と点を大きく、さらに全サンプルをプロットします。
par (mar =c(6,6,3,3),mgp=c(4,1,0), lwd=2) # 図の余白などの調節
plot (model, cex=2, cex.axis=2, cex.lab=2, lwd=2, type="all" ) # 詳しくはGoogle先生
対数軸に0が書かれているのは気持ち悪い、ということでControlを離して描きます。
par (mar =c(6,6,3,3),mgp=c(4,1,0), lwd=2)
plot (model, cex=2, cex.axis=2, cex.lab=2, lwd=2, type="all", broken=TRUE, conName="Control" )
信頼区間を重ねて描きます。
par (mar =c(6,6,3,3),mgp=c(4,1,0), lwd=2)
plot (model, cex=2, cex.axis=2, cex.lab=2, lwd=2, type="all", broken=TRUE, conName="Control",ylab="Mortality")
plot (model, cex=2, cex.axis=2, cex.lab=2, lwd=2, type="confidence", broken=TRUE, conName="Control", add=TRUE, confidence.level=0.95, ylab="Mortality")# add=TRUEで重ね描きができる
縦軸をパーセント表記にしたかったのですが、drcパッケージのplot関数でそのような機能をサポートしてるか分からず、面倒なのでやめました。パーセント表記にしたい人は、普通のplot関数で致死データをプロットした後で、log-logistic式を重ねると良いかと思います。
最後にどのモデル式に当てはめるのが良いのか、という話を書こうかと思いましたが、さらに長くなりそうなのでここまでにします。気が向いたら続きを書くかも。
ある化学物質や汚染環境試料のLC50やEC50を知りたいだけならば、基本的にlog-logisticを使えば大きな問題はないと思います。たぶん。詳しくはdrcの作者Ritzさんの論文を読んでみましょう。
K. HIKI