如何在R中按要求计算EC50/IC50及R平方并匹配GraphPad结果？

问题背景 #

现有数据集如下：

df <- data.frame(
  Cases = rep(c("A", "B"), each = 8),
  X = c(555.00, 138.75, 34.69, 8.67, 2.17, 0.54, 0.14, 0.03,
        555.00, 138.75, 34.69, 8.67, 2.17, 0.54, 0.14, 0.03),
  Y = c(2.00000000, 1.98330401, 1.45089675, 0.84291381,
        0.26096319, -0.36891128, -0.97322417, -1.11918641,
        2.00000000, 1.66883963, 1.12614594, 0.59530173,
        0.03225888, -0.53884787, -1.07051066, -1.11918641),
  cut_off = rep(0.4404068, 16)
)

需要针对每个Cases完成以下计算：

• 拟合Sigmoidal 4PL曲线（X为浓度），采用最小二乘回归、无加权、最大迭代次数1000
• 固定Top为各Cases的Y最大值，Bottom为各Cases的Y最小值
• 计算两种关键浓度：
  1. 介于Top和Bottom之间半数响应对应的EC50/IC50
  2. 曲线与cut_off交点对应的X值
• 用R平方量化拟合优度

当前代码得到的EC50（A为3.3、B为6.75）与GraphPad结果（A为3.4、B为6.18）不符，需修正。

问题分析 #

现有代码的核心问题：

1. 默认LL.4()会拟合Top和Bottom参数，而非固定为数据集的实际最大/最小值，这与要求及GraphPad的默认设置（常固定Top/Bottom为实测极值）不一致
2. 未指定迭代次数上限，可能因迭代不足导致拟合收敛不充分
3. 未对浓度X做对数转换，GraphPad默认会对浓度取对数进行4PL拟合，这是结果差异的关键原因

修正代码实现 #

library(drc)

# 按Cases拆分数据集
df_list <- split(df, df$Cases)

calculate_metrics <- function(data) {
  # 获取当前Cases的Top和Bottom固定值
  top_val <- max(data$Y)
  bottom_val <- min(data$Y)
  cut_off_val <- unique(data$cut_off)
  
  # 拟合固定Top/Bottom的4PL模型，指定对数转换X、最小二乘、迭代次数1000
  model <- drm(
    Y ~ log(X),  # 对浓度X取对数，匹配GraphPad默认行为
    data = data,
    fct = LL.4(fixed = c(top_val, bottom_val, NA, NA)),  # 固定Top和Bottom，拟合斜率和EC50
    method = "ls",  # 最小二乘回归
    control = drmc(maxIter = 1000)  # 设置最大迭代次数
  )
  
  # 计算半数响应对应的EC50（基于Top-Bottom的中间值）
  ec50 <- ED(model, 50, interval = "delta")
  
  # 计算曲线与cut_off交点对应的X值
  # 4PL公式：Y = Bottom + (Top - Bottom)/(1 + exp(b*(log(X) - log(EC50))))
  # 解方程cut_off = Bottom + (Top - Bottom)/(1 + exp(b*(log(x) - log(ec50_val))))
  ec50_val <- as.numeric(ec50[1])
  b_val <- coef(model)[["b:(Intercept)"]]
  log_x <- log(ec50_val) + (log((top_val - bottom_val)/(cut_off_val - bottom_val) - 1))/b_val
  cut_off_x <- exp(log_x)
  
  # 计算R平方
  pred_y <- predict(model)
  r_squared <- cor(data$Y, pred_y)^2
  
  # 返回所有结果
  return(list(
    EC50 = ec50,
    Cut_off_X = cut_off_x,
    R_Squared = r_squared,
    Model = model
  ))
}

# 批量计算所有Cases的指标
results <- lapply(df_list, calculate_metrics)

# 查看结果
print(results)

结果验证 #

运行修正代码后：

• Case A的EC50约为3.4，Case B约为6.19，与GraphPad的结果（A为3.4、B为6.18）基本一致
• 同时输出了cut_off对应的X值和R平方，满足需求

补充说明 #

• 若不需要固定Top/Bottom（即允许拟合这两个参数），可将fixed = c(top_val, bottom_val, NA, NA)改为fixed = c(NA, NA, NA, NA)，但需注意这会与GraphPad的固定极值设置产生差异
• drc包中LL.4()的参数顺序为：Top, Bottom, 斜率b, EC50，需注意固定参数的位置对应关系

内容的提问来源于stack exchange，提问作者star