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

# はじめに

競合的ELISAを行ったときに検量線を書く場合、4係数ロジスティック曲線(4-parameter logistic curve)をあてはめます。4係数ロジスティック曲線についてはリンク先（[参考URL](https://qiita.com/fogdiffusion/items/ac628df4cb97f2992114)）を参照して下さい。

この検量線を使用して実際のサンプルの濃度を求める方法を調べたので残しておきます。

**drc**パッケージを使用する方法と、使用しない方法の2種類を紹介します。

何かの参考になれば幸いです。

# 解析方法
## ELISAデータ

ここでは架空のデータを使用します。

### 検量線量のデータ

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

<img width=60% alt="unnamed-chunk-2-1.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373740/0e740a1d-e3c0-6f7c-a880-1585729b25c2.png">

### サンプルの吸光度

```r
Sample <- c(2.01, 1.53, 0.769)
```

## drcパッケージを使用した場合
### 曲線あてはめ

曲線をあてはめるために`drm`関数を使用します。`LL.4`関数を指定することで4係数ロジスティック曲線をあてはめることができます。


```r
## drcパッケージをインストールしていない場合
## install.packages("drc", repos = "https://cloud.r-project.org")
library(drc)
model1 <- drm(OD ~ conc, fct = LL.4(), data = dat)
```

参考までに係数です。

```r
coef(model1)
```

```
## b:(Intercept) c:(Intercept) d:(Intercept) e:(Intercept) 
##     1.0516584     0.4122003     2.5275349     0.4419621
```

### 予測値

`ED`関数で予測値を得ることができます。この場合には`type = "absolute"`を指定することが必要です。結果は`matrix`クラスで返ってきます。


```r
res1 <- ED(model1, Sample, type = "absolute", display = FALSE)
res1
```

```
##            Estimate  Std. Error
## e:1:2.01  0.1513043 0.002518151
## e:1:1.53  0.3966227 0.004727397
## e:1:0.769 2.0141149 0.055280669
```

```r
res1[,"Estimate"]  # 値だけ欲しいとき
```

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


## 基本パッケージのnls関数を使用した場合

`drc`パッケージをインストールするといくつかの依存パッケージも同時にインストールされます。ただし今回の目的のためにはそれらの依存パッケージは必要ないため、基本パッケージだけで解析できるかどうかを調べました。

### 曲線あてはめ

非線形回帰を行うための`nls`関数を使用すれば曲線のあてはめを行うことができました。

初期値として`b`は1（たぶんこれでうまく行きます）、`c`はOD値の最小値、`d`はOD値の最大値、`e`は濃度の中央値を使用しました。


```r
model2 <- nls(OD ~ d + (c - d) / (1 + (e / conc)^b),
              data = dat,
              start = list(b = 1,
                           c = min(dat$OD),
                           d = max(dat$OD),
                           e = median(dat$conc)))
```

係数です。

```r
coef(model2)
```

```
##         b         c         d         e 
## 1.0516508 0.4121948 2.5275357 0.4419650
```
`drc`関数を使用したときとほぼ同じ値が得られました。

### 予測値

4係数ロジスティック曲線の逆関数を求めて関数にしました。`parm`は先程のモデルの係数、`OD`はサンプルの吸光度です。


```r
calc_LL4 <- function(parm, OD) {
    exp(log((parm[3] - OD)/(OD - parm[2])) / parm[1] + log(parm[4]))
}
```


```r
calc_LL4(coef(model2), Sample)
```

```
## [1] 0.1513038 0.3966234 2.0141216
```

先程とほぼ同じ結果が得られました。実用上問題のない差です。

# 最終的なスクリプト


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

## drcパッケージを使用する場合
library(drc)
model1 <- drm(OD ~ conc, fct = LL.4(), data = dat)
coef(model1)

res1 <- ED(model1, Sample, type = "absolute", display = FALSE)
res1
res1[,"Estimate"]  # 値だけ欲しいとき


## nls関数を使用する場合
model2 <- nls(OD ~ d + (c - d) / (1 + (e / conc)^b),
              data = dat,
              start = list(b = 1,
                           c = min(dat$OD),
                           d = max(dat$OD),
                           e = median(dat$conc)))
coef(model2)

calc_LL4 <- function(parm, OD) {
    exp(log((parm[3] - OD)/(OD - parm[2])) / parm[1] + log(parm[4]))
}

calc_LL4(coef(model2), Sample)
```

# おわりに

競合的ELISAで測定を行ったときの解析方法について紹介しました。

drcパッケージを含むdockerコンテナを作成したときにサイズが大きくなってしまったため、drcパッケージを使用しない方法を試したところ、このようになりました。