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

# はじめに

[前回の記事](https://qiita.com/toshi-ara/items/7a713276cb67daa5a10b)でRを使用して4係数ロジスティック曲線 (4-parameter logistic curve)をあてはめる方法を説明しました。その後、[この記事](https://www.pediatricsurgery.site/entry/2017/01/30/202953)から[`nplr`パッケージ](https://CRAN.R-project.org/package=nplr)の存在を知り、より簡単に目的の値を得ることができたので紹介します。

# 解析方法
## ELISAデータ

[前回の記事](https://qiita.com/toshi-ara/items/7a713276cb67daa5a10b)と同じデータを使用します。



```r
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`パッケージのインストール


```r
install.packages("nplr", repos = "https://cloud.r-project.org")
```
でインストールできます。依存パッケージはありません。

## 曲線あてはめ
### モデル式の作成

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


```r
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`関数を使用します。


```r
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
```

```r
res$x
```

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


```r
class(res)
```

```
## [1] "data.frame"
```
`getEstimates`の戻り値はデータフレームのようです。


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

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

```
## [1] 0.1512663 0.3958139 2.0115630
```

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

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

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


## あてはめ曲線の描画

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

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


```r
plot(model)
```

<img width=60% alt="unnamed-chunk-8-1.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373740/99745905-ff65-5689-1b45-7568e983aea9.png">


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

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

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


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

<img width=60% alt="unnamed-chunk-9-1.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373740/d9930dcd-dc9c-8354-6d22-5cc8def6b35a.png">


# 最終的なスクリプト


```r
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`パッケージを使用した方が関数を自分で定義しなくていいので楽かつ安心でした。