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つです。

f:id:YamamotoLab:20180211082914p:plain

 

下の表が実験結果だとします。表のDeadは致死個体の数、Totalはビーカーに入れた総数、Concは曝露濃度をそれぞれ示しています。対照区Controlは便宜的に0と書いておきます。

表:96時間致死毒性試験結果

f:id:YamamotoLab:20180212155530p:plain

 

 

 

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 )  

f:id:YamamotoLab:20180211134315p:plain

Rでは目的変数y、説明変数xの関係をチルダ"~"を使って"y~x"のように表します。xは単一の変数でなく、"y~x1+x2"のようなモデル式でも構いません。

他に指定せず単純に"Dead/Total ~ Conc"と書くと、シンプルな線形回帰式 "Dead/Total = α × Conc + β" を想定していることになります。

 

その単純な式を図にすると下のようになります。濃度が20 μg/Lを超えると致死率が1 (= 100%) をオーバーするというモデルです。

f:id:YamamotoLab:20180211140237p:plain

さすがにこれは不自然なので、ある濃度以上だと致死率は1で一定となり (上限があり) 、ある濃度以下なら致死率はほぼ0である (下限がある) ことを表現できるモデル式が必要です。

 

 

 

4. log-logistic式 

ではどんなモデル式にあてはめましょう?毒性学分野で広く使われている式にlog-logistic式があります。log-logistic式は、

f:id:YamamotoLab:20180211144825p:plain

で表現できます。xはの濃度で、b,c,d,eは回帰係数です。

これらの係数を変化させてlog-logistic式を描いたものが下図です。log-logistic式はS字型(シグモイド)の曲線を表現できることが分かりますね。また図を見ると、bは傾き、cは下限値、dは上限値、eはf(x)が上限値の半分となる際のxの値(毒性学的にいうとED50, Effective Dose50)に対応していることが読み取れます。

f:id:YamamotoLab:20180212133347p:plain

 

ちなみに致死率を目的変数とする場合、致死率の下限は0です。なので、係数c = 0として下のように係数3つのlog-logistic式を書くこともできます。

f:id:YamamotoLab:20180211144852p:plain

同様に致死率の上限は1ですから、さらに係数d=1として係数2つの式にすることもできます。

f:id:YamamotoLab:20180212130819p:plain

 

都市工の学部生にもおそらく馴染み深い話をすると、係数3つのlog-logistic式はb=-1の時に酵素の反応速度論で用いられるミカエリス・メンテン式と等しくなります。式変形をしてみます。(bの値が間違っているとのご指摘を受けて修正しました。2018.02.28)

f:id:YamamotoLab:20180228183124p:plain

ちなみにミカエリスメンテン式は↓。Vmaxは基質濃度xが最大の時の反応速度で、Kmはv=Vmax/2の時の基質濃度ですね。

f:id:YamamotoLab:20180212112922p:plain

f:id:YamamotoLab:20180212113711p:plain

▲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)

f:id:YamamotoLab:20180212144601p:plain

これではダサいし見づらいので、色々付け加えていきます。

まず、線を太く、字と点を大きく、さらに全サンプルをプロットします。

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先生

f:id:YamamotoLab:20180212150111p:plain

 

対数軸に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" )

f:id:YamamotoLab:20180212150443p:plain

 

信頼区間を重ねて描きます。

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=TRUEconfidence.level=0.95, ylab="Mortality")# add=TRUEで重ね描きができる

f:id:YamamotoLab:20180212152713p:plain

 

縦軸をパーセント表記にしたかったのですが、drcパッケージのplot関数でそのような機能をサポートしてるか分からず、面倒なのでやめました。パーセント表記にしたい人は、普通のplot関数で致死データをプロットした後で、log-logistic式を重ねると良いかと思います。

 

 

 

最後にどのモデル式に当てはめるのが良いのか、という話を書こうかと思いましたが、さらに長くなりそうなのでここまでにします。気が向いたら続きを書くかも。

ある化学物質や汚染環境試料のLC50やEC50を知りたいだけならば、基本的にlog-logisticを使えば大きな問題はないと思います。たぶん。詳しくはdrcの作者Ritzさんの論文を読んでみましょう。

 

 

K. HIKI