Rによる4係数ロジスティック曲線のあてはめ・予測値（その2）

はじめに

前回の記事でRを使用して4係数ロジスティック曲線 (4-parameter logistic curve)をあてはめる方法を説明しました。その後、この記事からnplrパッケージの存在を知り、より簡単に目的の値を得ることができたので紹介します。

解析方法

ELISAデータ

前回の記事と同じデータを使用します。

dat <- data.frame(
    conc = c(0.02, 0.15, 0.4, 1, 2, 5),
    OD = c(2.45, 2.01, 1.53, 1.04, 0.769, 0.567)
)

Sample <- c(2.01, 1.53, 0.769)

nplrパッケージのインストール

install.packages("nplr", repos = "https://cloud.r-project.org")

でインストールできます。依存パッケージはありません。

曲線あてはめ

モデル式の作成

nplr関数はformulaを使えません。nparsオプションでパラメータ数を指定することができます（2から5）。デフォルトでは最もあてはまりがよいものを選んでくれます。
getPar関数でパラメータを得ることができます。

library(nplr)
model <- nplr(dat$conc, dat$OD, silent = TRUE)
getPar(model)

## $npar
## [1] 5
## 
## $params
##      bottom      top      xmid       scal       s
## 1 0.4459926 2.540928 -0.212955 -0.9831879 1.30527

この場合、5-パラメータの式が最もあてはまりがよいという結果になりました。

パラメータのうちxmidがEC50です（常用対数として算出）。したがって10^xmidでEC50値を求めることができます。この場合にはEC50 = 0.612となりました。

予測値の算出

次に、OD値から濃度を求めます。この場合にはgetEstimates関数を使用します。

res <- getEstimates(model, Sample)
res

##       y     x.025         x     x.975
## 1 2.010 0.1492970 0.1501047 0.1509243
## 2 1.530 0.3976332 0.3992289 0.4007904
## 3 0.769 1.9782180 1.9918028 2.0059162

res$x

## [1] 0.1501047 0.3992289 1.9918028

予測値および95%信頼区間を算出してくれます。ブートストラップ法で信頼区間を求めているらしく、毎回信頼区間の値が異なりました（ごくわずかな違い）。ここではxカラムのデータだけを取り出しました。

class(res)

## [1] "data.frame"

getEstimatesの戻り値はデータフレームのようです。

次に4-パラメータの場合を試してみます

model4 <- nplr(dat$conc, dat$OD, npars = 4, silent = TRUE)
res <- getEstimates(model4, Sample)
res$x

## [1] 0.1512663 0.3958139 2.0115630

どの方法を使用してもほぼ同じ結果が得られました。

（参考）drcパッケージを使用したときの値

##  e:1:2.01  e:1:1.53 e:1:0.769
## 0.1513043 0.3966227 2.0141149

あてはめ曲線の描画

散布図にあてはめ曲線を追加した図です。

nplrパッケージで用意されている図

plot(model)

カスマタイズの方法が分かりません。

ggplot2パッケージを使用した図

nplr関数の戻り値に、あてはめ曲線用のデータが含まれているためそれを使用します。それぞれgetXcurve（ただし常用対数として）およびgetYcurve関数で値を得ることができます。

library(ggplot2)
ggplot(dat, aes(conc, OD)) +
    geom_line(data = data.frame(
                  conc = 10^getXcurve(model),
                  OD = getYcurve(model))) +
    geom_point() +
    scale_x_log10()

最終的なスクリプト

dat <- data.frame(
    conc = c(0.02, 0.15, 0.4, 1, 2, 5),
    OD = c(2.45, 2.01, 1.53, 1.04, 0.769, 0.567)
)

Sample <- c(2.01, 1.53, 0.769)

library(nplr)
model <- nplr(dat$conc, dat$OD, silent = TRUE)
## model <- nplr(dat$conc, dat$OD, npars = 4, silent = TRUE)
getPar(model)

res <- getEstimates(model, Sample)
res

library(ggplot2)
ggplot(dat, aes(conc, OD)) +
    geom_line(data = data.frame(
                  conc = 10^getXcurve(model),
                  OD = getYcurve(model))) +
    geom_point() +
    scale_x_log10()

おわりに

nplrパッケージを使用した方が関数を自分で定義しなくていいので楽かつ安心でした。