剂量响应 - 使用 R 的全局曲线拟合

debugcn 发表于 Dev

拉农·卡恩

我有以下剂量反应数据，并希望绘制剂量反应模型和全局拟合曲线。[xdata = 药物浓度；ydata(0-5) = 不同药物浓度下的响应值]。我没有问题地绘制了标准曲线。

标准曲线数据拟合：

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

剂量反应曲线数据拟合：

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

当我尝试使用下面提供的 R 脚本获取拟合参数时，出现以下错误：

nls(ydata1 ~ 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))

请建议我如何解决这个问题

G.格洛腾迪克

首先请注意，的最大值xdata与最小值的比率是 200 万，因此我们可能要使用log(xdata)代替xdata。

现在，进行此更改，我们得到了drc 包的 4 参数 log-logisticLL2.4模型，但参数化与问题中的略有不同。假设您对这些更改没有意见，我们可以按如下方式拟合第一个模型。有关参数化的详细信息，请参阅?LL2.4并参阅?ryegrass底部的相关示例。这是问题中显示的 - LL2.4 模型本身进行了转换。dfdflog(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

和以下图形适合（单击图像以展开它）：