ggplot2とdrcによる線量応答曲線のプロット

Question

生物学では、しばしば線量応答曲線をプロットする必要があります。 Rパッケージ 'drc'は本当に便利で、基本グラフィックは 'drmモデル'を簡単に扱うことができます。しかし、私はdrmカーブをggplot2に追加したいと思います。

マイセット：ベースのグラフィックスを使用して

library("drc") 
library("reshape2") 
library("ggplot2") 
demo=structure(list(X = c(0, 1e-08, 3e-08, 1e-07, 3e-07, 1e-06, 3e-06, 
1e-05, 3e-05, 1e-04, 3e-04), Y1 = c(0, 1, 12, 19, 28, 32, 35, 
39, NA, 39, NA), Y2 = c(0, 0, 10, 18, 30, 35, 41, 43, NA, 43, 
NA), Y3 = c(0, 4, 15, 22, 28, 35, 38, 44, NA, 44, NA)), .Names = c("X", 
"Y1", "Y2", "Y3"), class = "data.frame", row.names = c(NA, -11L 
)) 

：

plot(drm(data = reshape2::melt(demo,id.vars = "X"),value~X,fct=LL.4(),na.action = na.omit),type="bars") 

は素敵な4パラメータ用量応答プロットを生成します。

ggplot2に同じプロットをプロットすると、2つの問題が発生します。

1. drmモデルカーブを直接追加する方法はありません。私は関数として4-PLを書き直し、stat_functionの形でそれを追加する必要があります。

   ggplot(reshape2::melt(demo,id.vars = "X"),aes(X,value)) + 
       geom_point() + 
       stat_function(fun = function(x){ 
       drm_y=function(x, drm){ 
        coef(drm)[2]+((coef(drm)[3]-coef(drm)[2])/(1+exp((coef(drm)[1]*(log(x)-log(coef(drm)[4])))))) 
       } 
   + drm_y(x,drm = drm(data = reshape2::melt(demo,id.vars = "X"), value~X, fct=LL.4(), na.action = na.omit)) 
   }) 
   

2. これで十分でない場合は、scale_xが連続している場合にのみ機能します。 scale_x_log10()を追加したい場合は、 Warning message: In log(x): NaNs producedとなります。

これを実現するには、これを処理する方法があります。いずれか（plot.drcの場合）のように、x = 0の値はx軸上にプロットされ、基本的には最低x値の1/100となります。 （demo$X[which.min(demo$X)+1]/100）またはGraphPad Prismの場合と同様に、0は線量応答曲線から完全に除外されます。

私の質問は以下のとおりです。

1. 直接ggplot2でDRMモデルをプロットする方法はありますか？

2. データセットを対応する4-PL曲線とリンクして、同じ色でプロットする方法を教えてください。

Answer 1

私は自分自身の疑問に答えるつもりです。これは他の人が同じ問題に直面するのを助けるでしょう。

scale_x_log10()場合、リニアスケールまたはgeom_またはstat_smooth (method=drm, fct=L.4(),se=FALSE)にプロットする場合geom_またはstat_smooth (method=drm, fct=LL.4(),se=FALSE)のいずれかの単純な足し算でggplot2とDRCパッケージに用量反応曲線をプロットすることはもちろん可能ですが追加されます。

私は私のデータを変換しましたlog10のスケールを使用することができるようにするために：この場合

demo <- demo %>% 
     mutate(X = 
     ifelse(X == 0, 
       yes = (sort(demo$X[which.min(sort(demo$X)) + 1]/100)), 
       no = X 
      ) 
      )   #looks for the pre-lowest value in X and divides it by 100 

を、私は= 1/100分のXでX = 0値を交換してきました最後のX値（この場合は1e-10）。しかし、Prismのようにデータセットから完全に省略することで、対数プロットを駄目にする0値を簡単に削除することができます。 私が知ったように、ggplotは最初に軸をスケーリングしてからデータを追加するので、コードがlog10（0）にしようとするとブレークします。

もう一つの微妙な点は、stat_smooth関数がmethod = drmを使用してdrmモデルを完全に処理できることですが、 'SE'信​​頼区間にどのように適合するかはわかりません。 se = FALSEを選択すると、プロットが可能になり、私の謙虚な意見では、とにかくあまり目立たないプロットになります。エラーバーを追加するだけです。

最後に、fct = LL.4()をfct = L.4()に変更すると、スケールが最初に選択され、後にフィットが行われるため、log10スケールでプロットすることができます。したがって、軸の値が非対数であっても、ggplotは実際にlog10にデータセットを変換しているため、log-logit-4P（LL（））の代わりにフィッティング関数をlogit-4P（L.4 .4（））。

geom_smooth（）関数とstat_smooth（）関数は当然データセットと同じ色を採用するため、フィットした関数の色をデータポイントの色に合わせて調整する必要はありません。要約する

：

demo <- demo %>% 
     mutate(X = 
     ifelse(X == 0, 
       yes = (sort(demo$X[which.min(sort(demo$X)) + 1]/100)), 
       no = X 
      ) 
      ) 
demo.long <- reshape2::melt(demo,id.vars = "X") #reshapes the demo dataset to long format 
ggplot(data = demo.long, 
     aes(x = X, y = value, col = variable) 
    ) + 
    geom_point() + 
    geom_smooth(method = drm, fct = L.4(), se = FALSE) + 
    scale_x_log10() #plots out the dataset with the corresponding 4-parameter log-logit dose response curves 

Answer 2

recent paperdrcパッケージの作者からはggplot2が使用するパラメータを抽出するための命令を含んでいました。これらはggplot2では機能しませんが、モデルからデータを抽出します。これはあなたのデータに適用されるソリューションです。

demo1 <- reshape2::melt(demo,id.vars = "X") # get numbers ready for use. 
demo.LL.4 <- drm(data = demo1,value~X,fct=LL.4(),na.action = na.omit) # run model. 

predict

機能はdrmモデルからパラメータを抽出することができます。 curveidを使用してフィットした複数の曲線と互換性がありません。

# predictions and confidence intervals. 
demo.fits <- expand.grid(conc=exp(seq(log(1.00e-04), log(1.00e-09), length=100))) 
# new data with predictions 
pm <- predict(demo.LL.4, newdata=demo.fits, interval="confidence") 
    demo.fits$p <- pm[,1] 
    demo.fits$pmin <- pm[,2] 
    demo.fits$pmax <- pm[,3] 

coord_transの問題を避けるため、ゼロ濃度をシフトすることをお勧めします。

demo1$XX <- demo1$X 
demo1$XX[demo1$XX == 0] <- 1.00e-09 

はその後geom_ribbonが描画されてからのエラーを停止して省略して、曲線をプロットしています。

ggplot(demo1, aes(x = XX, y = value)) + 
    geom_point() + 
    geom_ribbon(data=demo.fits, aes(x=conc, y=p, ymin=pmin, ymax=pmax), alpha=0.2) + 
    geom_line(data=demo.fits, aes(x=conc, y=p)) + 
    coord_trans(x="log") 

プロセスを繰り返すことができる複数の曲線を一緒にグラフ化します。各セットにIDを追加します。

demo.fits_1 <- data.frame(label = "curve1", demo.fits) 

次に、rbindを使用して、抽出されたすべてのパラメータを結合します。そこからggplotは色を扱うことができます。