IC50计算R代码正确性验证及纠错请求

问题概述 #

用户编写了一段用于计算半最大抑制浓度（IC50）的R代码，代码可运行但无法确认结果正确性，需要测试并修正潜在错误。

核心错误分析 #

• 变量对应关系完全颠倒：Hill方程中，抑制剂浓度C是自变量，观测到的响应值v1是因变量。原代码错误地将C代入V函数得到y_spa2，再以v1为x轴、y_spa2为y轴做插值，完全混淆了自变量和因变量的逻辑，导致IC50的计算对象错误。
• 重复添加噪声：代码中两次执行v1 = v1 + runif(length(v1), min = -0.1, max = 0.1)，过度叠加噪声会严重偏离原始观测数据。
• IC50定义误解：IC50是响应值降至0.5时对应的抑制剂浓度，原代码错误地将插值后的响应值序列当作浓度来取值，结果毫无意义。
• 插值范围不合理：插值序列x_dense2的范围（0-10）未覆盖原始浓度的最大值（19.5），且方向错误（应基于浓度范围插值）。

修正后的代码 #

# 抑制剂浓度（自变量）
C <- c(0, 0.3, 1.5, 3.5, 19.5)

# 观测响应值（因变量，代表活性保留比例，范围0-1）
v1 <- c(0.00, 0.00, 0.00, 0.25, 0.90)

# 仅添加一次随机噪声模拟实验误差
v1 <- v1 + runif(length(v1), min = -0.1, max = 0.1)
# 确保响应值在合理范围（0-1）内
v1 <- pmax(pmin(v1, 1), 0)

# Hill方程参数：H为IC50初始猜测，n为Hill系数（斜率）
H <- 1
n <- 1

# 定义Hill抑制模型：输入浓度C，返回理论响应值
V <- function(C, H, n) {
  1 / (1 + (C / H)^n)
}

# 生成高密度的浓度序列用于插值/绘图
x_dense <- seq(min(C), max(C), by = 0.01)

# 对观测数据进行线性插值，得到连续的响应曲线
# 这里以浓度C为x，观测响应v1为y，插值得到高密度浓度对应的响应值
y_dense <- approx(C, v1, xout = x_dense, method = "linear")$y

# 寻找响应值首次降到0.5以下对应的浓度（IC50）
# 先找到所有y_dense <= 0.5的索引，取第一个对应的浓度
if (any(y_dense <= 0.5)) {
  IC50 <- x_dense[which(y_dense <= 0.5)[1]]
} else {
  # 处理所有响应值都高于0.5的情况（无有效IC50）
  IC50 <- NA
  warning("所有观测响应值均高于0.5，无法计算IC50")
}

# 绘图验证
plot(x_dense, y_dense, type = "l", xlab = "抑制剂浓度", ylab = "响应值（活性保留比例）", main = "抑制曲线与IC50")
# 添加原始观测点
points(C, v1, pch = 16, col = "red")
# 添加y=0.5的水平线（半抑制点）
abline(h = 0.5, col = "blue", lty = 2)
# 添加IC50的垂直线
if (!is.na(IC50)) {
  abline(v = IC50, col = "green", lty = 2)
  text(IC50, 0.5, paste0("IC50 = ", round(IC50, 2)), pos = 4, col = "green")
}

验证方式 #

用已知参数的模拟数据测试：

• 当设置真实IC50为H=1，n=1时，生成理想响应值v1_true <- V(C, 1, 1)，添加少量噪声后运行代码，计算得到的IC50应接近1，以此验证代码逻辑正确性。

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