使用drm绘制量效曲线时出现Non finite value error的技术问询

错误原因 #

出现non-finite value supplied by optim和收敛失败的问题，通常由以下几种情况导致：

• 数据无明确量效趋势：你的AR231数据可能没有呈现出随浓度变化的典型响应规律（比如示例数据中浓度升高响应值明显下降），导致优化算法找不到合理的拟合方向。
• 初始值不合适：drm()默认的参数初始值可能与你的数据分布不匹配，使得优化过程陷入局部极值或无法计算有效梯度。
• 数据异常：存在离群值、缺失值或浓度范围不足（比如缺少低浓度/高浓度的平台期数据），干扰了模型拟合。
• 模型选择不当：LL.4()四参数逻辑模型可能不适合你的数据（比如如果响应的下限平台接近0，三参数模型LL.3()更合适）。

解决步骤 #

1. 先诊断数据 #

首先可视化AR231的原始数据，确认是否存在量效趋势：

library(ggplot2)
ggplot(dosage, aes(x = conc, y = AR231)) +
  geom_point() +
  scale_x_log10() + # 浓度通常取对数更直观
  labs(x = "浓度", y = "响应值")

同时检查数据的基本统计信息，排查异常值：

summary(dosage$AR231)
sum(is.na(dosage$AR231)) # 检查缺失值

2. 调整模型拟合策略 #

尝试简化模型

如果数据的下限平台（比如高浓度时的响应值）接近0，改用三参数对数逻辑模型：

AR231 <- drm(AR231 ~ conc, data = dosage, fct = LL.3())
summary(AR231)

提供自定义初始值

根据数据的实际范围手动指定初始参数，帮助优化算法收敛：

• b：斜率（示例中约1.4，可根据数据趋势调整）
• c：下限平台（数据最小值附近）
• d：上限平台（数据最大值附近）
• e：EC50（数据响应下降一半对应的浓度）

# 假设AR231的最大值是4，最小值是0，EC50大概在20左右
AR231 <- drm(AR231 ~ conc, data = dosage, fct = LL.4(),
             start = c(b = 1, c = 0, d = 4, e = 20))
summary(AR231)

更换优化方法

默认的优化方法可能不适用，尝试BFGS或Nelder-Mead：

AR231 <- drm(AR231 ~ conc, data = dosage, fct = LL.4(),
             method = "BFGS")
summary(AR231)

3. 批量获取所有受试对象的EC50值 #

如果需要为所有列（除了conc）拟合模型并提取EC50，可以用循环或lapply实现：

library(drc)

# 获取所有受试对象的列名
subject_cols <- colnames(dosage)[!colnames(dosage) %in% "conc"]

# 定义拟合函数
fit_drc <- function(col_name) {
  model <- tryCatch(
    drm(get(col_name) ~ conc, data = dosage, fct = LL.4()),
    error = function(e) NULL # 跳过拟合失败的对象
  )
  if (!is.null(model)) {
    ec50 <- ED(model, 50)[[1]]
    return(data.frame(Subject = col_name, EC50 = ec50))
  } else {
    return(data.frame(Subject = col_name, EC50 = NA))
  }
}

# 批量拟合并整理结果
ec50_results <- do.call(rbind, lapply(subject_cols, fit_drc))
print(ec50_results)

注意事项 #

• 如果某些受试对象的数据确实没有量效关系，拟合失败是正常的，可考虑排除这些样本或调整模型。
• 拟合前建议对浓度取对数（比如log10(conc+1)避免0值问题），有时能改善收敛效果。

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