剂量反应 - 使用 R 的全局曲线拟合答案

【问题标题】：Dose Response - Global curve fitting using R剂量反应 - 使用 R 的全局曲线拟合
【发布时间】：2017-02-26 07:35:36
【问题描述】：

我有以下剂量反应数据，并希望绘制剂量反应模型和全局拟合曲线。 [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))

请告诉我如何解决这个问题

【问题讨论】：

标签： r curve-fitting

【解决方案1】：

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

现在，进行此更改，我们得到了 drc package 的 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

以及以下图形拟合（单击图像将其展开）：

【讨论】：

e 是 log(EC50)。如果您使用LL.4（而不是LL2.4），那么e 将是EC50。请阅读提供的链接。
先生，感谢您让我开心。一个快速的可能愚蠢的问题。我将如何根据 EC50 值解释拟合数据？因为我需要所有 5 次拟合的 EC50 值并用它来绘制另一个图。我应该根据拟合值计算 EC50 吗？
当我在调整浓度范围 [ xdata = c(-5.00, -5.30, -5.60, -5.90, -6.20, -6.51, -6.81, - 7.11，-7.41，-7.71，-8.01，-8.31，-8.61，-8.91，-9.21，-9.52，-9.82，-10.12，-10.42，-10.72，-11.02，-11.32）] 我收到以下错误：lm.fit（x，y，offset=offset，singular.ok=singular.ok，...）中的错误：0（非NA）案例
答案中显示的代码使用问题中显示的df，而不是它的某些转换版本。请仔细阅读答案中提供的链接以充分了解模型。xdata 的日志转换由 LL2.4 或 LL.4 模型完成，而不是您。
非常感谢您的解释。

【解决方案2】：

只是上面 G.Grothendieck 答案的附录，用于叠加图，以防万一。

library(drc)
ys <- names(df)[-1]
for (i in 1:ys)  
 {fo <- as.formula(paste(ys[i], "~ xdata"))
  fit <- do.call("drm", list(fo, data = quote(df), fct = quote(LL2.4())))
    plot(fit, pch = 19+ x, ylim = c( min(df[,-1]),max(df[,-1])))
   par(new=TRUE)
  fit}

【讨论】：