R で用量反応曲線
1. はじめに
用量反応曲線 (dose-response curve) を描く必要があったので、drc パッケージ を調査しました。
Bioassay Analysis using R が非常に分かりやすかったんですが、この資料が書かれてから drc パッケージは大幅に変更されていて、データセット名やメソッド名も変更されているので注意が必要です。
本記事は基本的に上記資料の “4. Fitting a single dose response curve” をもとに作成されています。
2. データセット
今回使用するデータセットは ryegrass データセットです。
library(drc)
data <- ryegrass
head(data)## rootl conc
## 1 7.580000 0
## 2 8.000000 0
## 3 8.328571 0
## 4 7.250000 0
## 5 7.375000 0
## 6 7.962500 0このデータはホソムギ (perennial ryegrass) をフェルラ酸 (ferulic acid) に暴露した時の濃度 conc と根の長さ rootl のデータだそうです。
plot(rootl ~ conc, data)
プロットしてみると、フェルラ酸濃度が上がると根の長さが下がっているようです。
また、濃度は対数にしたほうがよさそうです。
plot(rootl ~ conc, data, subset=conc!=0, log="x")
このデータに対してシグモイド曲線を当てはめ、用量反応曲線を描きたいと思います。
3. モデルの当てはめ
drc パッケージにはシグモイド曲線を表すことのできる関数がいくつか付属しています。
一番有名なのは、ロジスティック関数 (logistic)
\[L.4(x) = c + \frac{d-c}{1+\exp(b(x - e))}\] で、これは 4 つのパラメータ (b, c, d, e) を持つことから 4 パラメータロジスティック関数と呼ばれます。
\(c=0\) とした場合が 3 パラメータロジスティック関数です。
その他にも、ゴンペルツ関数 (Gompertz)
\[G.4(x) = c + (d-c)(\exp(-\exp(b(x-e))))\] ワイブル関数 (Weibull)
\[W1.4(x) = c + (d-c) \exp(-\exp(b(\log(x)-\log(e)))) \\ W2.4(x) = c + (d-c) (1 - \exp(-\exp(b(\log(x)-\log(e)))))\] などがあります。
用量-反応関係がホルミシス (hormesis) の場合は Brain-Cousens 関数などを使います。
drc パッケージ内の関数とモデルの対応を一部下記に示します。
| 関数 | モデル |
|---|---|
G.4() |
4 パラメータゴンペルツ |
L.4() |
4 パラメータロジスティック |
LL.4() |
4 パラメータ log ロジスティック |
W1.4() |
4 パラメータワイブル 1 |
W2.4() |
4 パラメータワイブル 2 |
BC.4() |
4 パラメータ Brain-Cousens |
これらの関数に対して非線形回帰モデル \[y_{i} = f(x_{i}) + \varepsilon_{i}\] を考えパラメータを推定します。
パラメータの推定は、drc では次のように行います。
model.LL4 <- drm(rootl ~ conc, data=data, fct=LL.4())
summary(model.LL4)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 2.98222 0.46506 6.41251 0.0000
## c:(Intercept) 0.48141 0.21219 2.26876 0.0345
## d:(Intercept) 7.79296 0.18857 41.32722 0.0000
## e:(Intercept) 3.05795 0.18573 16.46440 0.0000
##
## Residual standard error:
##
## 0.5196256 (20 degrees of freedom)drm() は fct で指定された関数に対してパラメータを推定します。
ここでは、4 パラメータ log ロジスティック関数 (LL.4()) に対して、パラメータ b, c, d, e を推定しています。
結果を見るには summary() を使います。
b は反応中間点での曲線の傾き、c は最小値、d は最大値、e は反応中間点での用量 (ED50) を表します。
c に対する p-value が大きい場合は、3 パラメータを検討します。
推定されたパラメータだけを見たいときは coef() も使えます。
coef(model.LL4)## b:(Intercept) c:(Intercept) d:(Intercept) e:(Intercept)
## 2.9822191 0.4814132 7.7929583 3.0579550モデルの適合度を調べたいときは modelFit() を使います。
modelFit(model.LL4)## Lack-of-fit test
##
## ModelDf RSS Df F value p value
## ANOVA 17 5.1799
## DRC model 20 5.4002 3 0.2411 0.8665ANOVA は、この場合(single curveの場合)、一元配置分散分析のことで、データに対する最小残差モデルになります。
したがって、当てはめたモデルが ANOVA と比べてどのくらい悪いかを残差を使って調べることで、モデルの適合度を知ることができます。
p-value が大きい場合、ANOVA と同程度の当てはめを持つことになります。
4. 用量反応曲線の描画
用量反応曲線の描画は plot() によって行います。
plot(model.LL4)
ここでプロットされている点は、各用量での平均です。
全データをプロットしたい場合は type="all" を指定します。type="none" で曲線のみ、type="obs" でデータのみをプロットすることもできます。
また、type="bars" でエラーバーを描画することもできます。
plot(model.LL4, type="all")
plot(model.LL4, type="bars")
plot(model.LL4, type="none")
plot(model.LL4, type="obs")
用量 0 は通常、コントロールとして分けて考えます。
引数 conName を指定することでコントロール名を変更できます。
また、broken=TRUE でコントロールを別扱いにできます。
plot(model.LL4, conName="Control", broken=TRUE)
別の曲線を追加するには plot() の引数 add=TRUE とします。
model.W14 <- drm(rootl ~ conc, data=data, fct=W1.4())
model.W24 <- drm(rootl ~ conc, data=data, fct=W2.4())
plot(model.LL4, conName="Control", broken=TRUE)
plot(model.W14, broken=TRUE, add=TRUE, type="none", col="green")
plot(model.W24, broken=TRUE, add=TRUE, type="none", col="blue")
5. Box-Cox 変換
Box-Cox 変換は boxcox() で行うことができます。
model.LL4.BC <- boxcox(model.LL4)
summary(model.LL4.BC)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 2.61839 0.39151 6.68795 0.0000
## c:(Intercept) 0.39083 0.10429 3.74744 0.0013
## d:(Intercept) 7.86633 0.29558 26.61364 0.0000
## e:(Intercept) 3.01662 0.21005 14.36124 0.0000
##
## Residual standard error:
##
## 0.2962958 (20 degrees of freedom)
##
## Non-normality/heterogeneity adjustment through optimal Box-Cox transformation
##
## Estimated lambda: 0.5
## Confidence interval for lambda: [0.269,0.949]ANOVA ベースで行う場合は method="anova" を指定します。
model.LL4.BC2 <- boxcox(model.LL4, method="anova")
summary(model.LL4.BC2)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 2.529649 0.377566 6.699894 0.000
## c:(Intercept) 0.372673 0.096817 3.849239 0.001
## d:(Intercept) 7.891138 0.327323 24.108080 0.000
## e:(Intercept) 2.997122 0.224248 13.365232 0.000
##
## Residual standard error:
##
## 0.2806173 (20 degrees of freedom)
##
## Non-normality/heterogeneity adjustment through optimal Box-Cox transformation
##
## Estimated lambda: 0.424
## Confidence interval for lambda: [0.124,0.782]plot(model.LL4, conName="Control", broken=TRUE)
plot(model.LL4.BC2, broken=TRUE, add=TRUE, type="none", col="red")