用量レスポンス - グローバルカーブフィッティングR

Question

以下の用量応答データがあり、用量応答モデルとグローバルフィット曲線をプロットすることを希望します。 [xdata =薬物濃度; ydata（0-5）=薬物の異なる濃度における応答値]。私はStd Curveを問題なくプロットしました。

STDカーブデータのフィット：

df <- data.frame(xdata = c(1000.00,300.00,100.00,30.00,10.00,3.00,1.00,0.30, 
          0.10,0.03,0.01,0.00), 
       ydata = c(91.8,95.3,100,123,203,620,1210,1520,1510,1520,1590, 
          1620)) 

nls.fit <- nls(ydata ~ (ymax*xdata/(ec50 + xdata)) + Ns*xdata + ymin, data=df, 
       start=list(ymax=1624.75, ymin = 91.85, ec50 = 3, Ns = 0.2045514)) 

用量応答曲線のデータフィット：私はRスクリプトを使用してフィットパラメータは以下に取得しようと

df <- data.frame(
     xdata = c(10000,5000,2500,1250,625,312.5,156.25,78.125,39.063,19.531,9.766,4.883, 
       2.441,1.221,0.610,0.305,0.153,0.076,0.038,0.019,0.010,0.005), 
     ydata1 = c(97.147, 98.438, 96.471, 73.669, 60.942, 45.106, 1.260, 18.336, 9.951, 2.060, 
        0.192, 0.492, -0.310, 0.591, 0.789, 0.075, 0.474, 0.278, 0.399, 0.217, 1.021, -1.263), 
     ydata2 = c(116.127, 124.104, 110.091, 111.819, 118.274, 78.069, 52.807, 40.182, 26.862, 
        15.464, 6.865, 3.385, 10.621, 0.299, 0.883, 0.717, 1.283, 0.555, 0.454, 1.192, 0.155, 1.245), 
     ydata3 = c(108.410, 127.637, 96.471, 124.903, 136.536, 104.696, 74.890, 50.699, 47.494, 23.866, 
        20.057, 10.434, 2.831, 2.261, 1.085, 0.399, 1.284, 0.045, 0.376, -0.157, 1.158, 0.281), 
     ydata4 = c(107.281, 118.274, 99.051, 99.493, 104.019, 99.582, 87.462, 75.322, 47.393, 42.459, 
        8.311, 23.155, 3.268, 5.494, 2.097, 2.757, 1.438, 0.655, 0.782, 1.128, 1.323, 0.645), 
     ydata0 = c(109.455, 104.989, 101.665, 101.205, 108.410, 101.573, 119.375, 101.757, 65.660, 35.672, 
        31.613, 12.323, 25.515, 17.283, 7.170, 2.771, 2.655, 0.491, 0.290, 0.535, 0.298, 0.106)) 

、私が手次のエラー：

nls（yda TA1は〜BOTTOM +（TOP - BOTTOM）/（1 + 10 ^（（logEC50 - XDATA）*：
単数勾配

nls.fit1 <- nls(ydata1 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata1), BOTTOM = min(df$ydata1),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

nls.fit2 <- nls(ydata2 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata2), BOTTOM = min(df$ydata2),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

nls.fit3 <- nls(ydata3 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata3), BOTTOM = min(df$ydata3),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

nls.fit4 <- nls(ydata4 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata4), BOTTOM = min(df$ydata4),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

nls.fit5 <- nls(ydata0 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata0), BOTTOM = min(df$ydata0),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

Answer 1

まずノートこの問題を解決する方法で私をアドバイスしてくださいxdataの最大値と最小値の比が200万であるため、xdataの代わりにlog(xdata)を使用したいと考えています。

この変更を加えると、の4つのパラメータlog-logistic LL2.4モデルが得られますが、質問とはわずかに異なるパラメータ設定が行われます。あなたがこれらの変更で大丈夫であると仮定すると、次のように最初のモデルに適合することができます。パラメータ化の詳細については?LL2.4を参照し、?ryegrassの下部に関連する例を参照してください。ここでdfはdfです。LL2.4モデル自体がlog(xdata)になります。

library(drc) 

fm1 <- drm(ydata1 ~ xdata, data = df, fct = LL2.4()) 
fm1 
plot(fm1) 

ここではすべての5つのモデルに適合し、視覚的にはフィットがかなり良いという最終的なプロットからわかります。

与える

library(drc) 

fun <- function(yname) { 
    fo <- as.formula(paste(yname, "~ xdata")) 
    fit <- do.call("drm", list(fo, data = quote(df), fct = quote(LL2.4()))) 
    plot(fit) 
    fit 
} 

par(mfrow = c(3, 2)) 
L <- Map(fun, names(df)[-1]) 
par(mfrow = c(1, 1)) 

sapply(L, coef) 

：

     ydata1 ydata2 ydata3 ydata4 ydata0 
    b:(Intercept) -1.37395 -1.1411 -1.1337 -1.0633 -1.6525 
    c:(Intercept) 0.70388 1.9364 1.5800 1.3751 5.7010 
    d:(Intercept) 101.02741 122.0825 120.8042 108.2420 107.9106 
    e:(Intercept) 6.17225 5.0686 4.3215 3.7139 3.2813 

と、次のグラフィカルなフィット（それを拡大するには画像をクリックしてください）：