Note of Pediatric Surgery

腸内細菌、R、ときどき小児外科

ELISAの測定結果をRで計算する ( 2 ) 4パラメーターロジスティック曲線を描く

それでは、4パラメーターロジスティック曲線を描いて、回帰式からELISAの値を算出していきましょう。Rっで4パラメーターロジスティック曲線を描く方法はいくつかあるようですが、drcパッケージを使用します。こちらのサイトを参考にさせていただきました。ちなみにラボで使用している解析ソフトが、SoftMax Proですので、それに合わせて話を進めさせてもらいます。

まずdrcをダウンロードしてインストールしておきます。ついでにggplot2とreshape2も立ち上げておきましょう。

[code lang=r] install.packages("drc", dependencies = TRUE) library(drc) library(ggplot2) library(reshape2) [/code]

次にELISAの測定データの.txtファイルから96well分の吸光度データの部分をコピーした状態で下記のコマンドを入力します。

[code lang=r] X = data.frame(read.delim(pipe("pbpaste"))) head(X) [/code]

こんな感じのデータだとします。

1 2 3 4 5 6 7 8 9 10 11 12
0.9324 0.9336 0.2707 0.2687 0.4402 0.4419 0.4682 0.4682 0.6406 0.6379 0.5189 0.527
0.8577 0.8381 0.3158 0.316 0.455 0.4504 0.1769 0.1702 0.6805 0.6555 0.4735 0.5123
0.6447 0.6339 0.219 0.231 0.7963 0.7697 0.5783 0.5886 0.5299 0.5327 0.6142 0.669
0.3987 0.4057 0.1891 0.1932 0.5592 0.5629 0.4805 0.4732 0.5179 0.4979 0.2384 0.2657
0.2288 0.2384 0.2236 0.2154 0.5094 0.5118 0.6234 0.6319 0.4776 0.4641 0.2501 0.2783
0.134 0.1565 0.2721 0.2722 0.672 0.667 0.6897 0.7063 0.6585 0.66 0.1852 0.2077
0.0942 0.1035 0.3887 0.4028 0.7222 0.7013 0.5786 0.5832 0.6574 0.6511 0.2408 0.2732
0.0564 0.0581 0.2768 0.3425 0.7196 0.7106 0.4907 0.4986 0.6147 0.6229 0.2005 0.2149

僕は最初の2列 ( X1とX2 ) にduplicateでスタンダードを置くようにしています。4パラメーターロジスティック曲線をまず書くのに必要なのはスタンダードの値なので、それを抜き出してデータフレームにしておきます。

[code lang=r]

.最初の2列を抜き出す

Std1 <-X[,1:2]

.スタンダードのIgA濃度のベクトルを作成する

Concentration <- as.mumeric(c("1000","500","250","125","62.5","31.25","15.625","0"))

.スタンダードの吸光度のベクトルと結合する

Std2 <- cbind(Concentration,Std)

.Duplicateの吸光度を1列にする

Std3 <- melt(Std2,id.vars="Concentration",variable.name="variable",value.name="Absorbance",na.rm=FALSE, factorsAsStrings=TRUE)

.不要な2列目を削除する

Std <- Std3[,-2] Std [/code]

と入力すると

Concentration Absorbance 1 1000 0.9324 2 500 0.8577 3 250 0.6447 4 125 0.3987 5 62.5 0.2288 6 31.25 0.1340 7 15.625 0.0942 8 0 0.0564 9 1000 0.9336 10 500 0.8381 11 250 0.6339 12 125 0.4057 13 62.5 0.2384 14 31.25 0.1565 15 15.625 0.1035 16 0 0.0581

と出力されると思います。これで準備が整いました。それでは散布図を書いてみましょう。

[code lang=r]

.散布図を描画する

plot(Absorbance~Concentration,Std, subset=Conc!=0, #.Blankの値は除去するのが通常のようです log="x", #.X軸を対数表示に col="black",pch=19) [/code]

するとこんな感じになると思います。

<a href="http://pedsurgery.wp.xdomain.jp/wp-content/uploads/2015/11/ELISA_Std1.png"> class="aligncenter size-full wp-image-416" src="http://pedsurgery.wp.xdomain.jp/wp-content/uploads/2015/11/ELISA_Std1.png" alt="ELISA_Std1" width="900" height="900" /></a>

確かにシグモイド曲線になっていますね。それでは次に4パラメーターロジスティック回帰を行います。パッケージdrcの<strong>drm()</strong>を使用します。

[code lang=r] LL4 <- drm(Absorbance~Concentration, data=Std, fct=LL.4()) #.LL.4:4パラメータlogロジスティック [/code]

fctの部分を変更することで、他の回帰も使用できるようです。

余談ですが、実際はfctの設定は " L.4 " が4パラメーターロジスティックのoptionなのですが、ここでは " LL.4 "  の4パラメータlogロジスティックを使っています。" log " がつくのとつかないのでどういう違いがあるのかわかりませんでしたが、

<strong>・参考にしたサイトでもLL.4を使用している</strong> <strong>・Softmax Proで自動で計算してくれる、係数に近かったのはlogありの方</strong> <strong>・実際に吸光度から濃度計算した値もlogありの方が信用できる</strong>

というところから fct=LL.4のoptionを選択しました。あとは散布図を描けばOKです。

[code lang=r] plot(LL4, type="all",col="red",pch=19) [/code]

とするとこんな感じの曲線が得られます。

<a href="http://pedsurgery.wp.xdomain.jp/wp-content/uploads/2015/11/Rplot.png"> class="aligncenter size-full wp-image-417" src="http://pedsurgery.wp.xdomain.jp/wp-content/uploads/2015/11/Rplot.png" alt="Rplot" width="1139" height="1272" /></a>

実際に吸光度から濃度を計算する際に重要なのは、信頼性の高い予測値はシグモイド曲線が直線になっている範囲でしか得られない、ということです。この4パラメーターロジスティック回帰曲線は、直線部分の範囲が狭く、あまり良いものではなさそうです。

ちなみにこの回帰曲線が妥当なものかを調べる際には、<a href="">https://www.evernote.com/shard/s10/res/b4c58bf6-9d1d-4559-b39d-6ce1d9166c88/20131010_slide.pdf">こちらの参考</a>によると、Hosmer‐Lemeshow検定 ( HL検定 ) というものを用いるようです。これはロジスティック回帰モデルの適合度、つまり全体の当てはまりの良さを検定する時に使用します。この場合の帰無仮説は「当てはめたモデルが正しい」ことなので、帰無仮説を棄却しない、つまりHL検定はp値が大きい方がいいモデルということです。

<a href="">http://rstudio-pubs-static.s3.amazonaws.com/63117_47264132789f40f690a707ee5429f45a.html">こちらのサイト</a>や<a href="inside-R">http://www.inside-r.org/packages/cran/drc/docs/modelFit">inside-RによるとパッケージdrcのmodelFit()で良さそうなのですが、これでHosmer‐Lemeshow検定が行えているのかはわかりません。一応やってみると

[code lang=r] modelFit(LL4) [/code]

Lack-of-fit test

ModelDf RSS Df F value p value ANOVA 8 0.00061952 DRC model 12 0.00177176 4 3.7198 0.0538

と返って来ます。TMBの反応時間が長かったのだと思うのですが、p=0.05にかなり近くなり、いまいちな結果になってしまいました。いろいろと条件検討をして修正することが必要ですね…。それでは次回の記事で、回帰式を求めて、吸光度から測定濃度を計算していきたいと思います。