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

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

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

install.packages("drc", dependencies = TRUE)
library(drc)
library(ggplot2)
library(reshape2)

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

X = data.frame(read.delim(pipe("pbpaste")))
head(X)

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

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パラメーターロジスティック曲線をまず書くのに必要なのはスタンダードの値なので、それを抜き出してデータフレームにしておきます。

#.最初の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

と入力すると

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

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

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

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

<a href="http://pedsurgery.wp.xdomain.jp/wp-content/uploads/2015/11/ELISA_Std1.png"><img 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>を使用します。

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

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

plot(LL4, type="all",col="red",pch=19)

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

<a href="http://pedsurgery.wp.xdomain.jp/wp-content/uploads/2015/11/Rplot.png"><img 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="http://www.inside-r.org/packages/cran/drc/docs/modelFit">inside-R</a>によるとパッケージdrcのmodelFit()で良さそうなのですが、これでHosmer‐Lemeshow検定が行えているのかはわかりません。一応やってみると

modelFit(LL4)

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にかなり近くなり、いまいちな結果になってしまいました。いろいろと条件検討をして修正することが必要ですね…。それでは次回の記事で、回帰式を求めて、吸光度から測定濃度を計算していきたいと思います。

Share on Facebook0Share on Google+0Share on Tumblr0Tweet about this on TwitterEmail this to someone

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

  1. ピンバック: ELISAの測定結果をRで計算する ( 3 ) 回帰式を使って濃度を求める | Note of Pediatric Surgery

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です