如何用R绘制Sigmoid曲线并从药物浓度数据拟合方程、获取R²值

针对你的药物浓度-响应数据，下面提供两种可靠的实现方法，包含拟合模型、获取方程及R²值、绘制美观曲线的完整步骤：

1. 推荐方案：使用drc包（专业剂量反应拟合） #

drc是专门用于剂量-反应曲线分析的R包，内置成熟的Sigmoid模型，拟合稳定性高，尤其适合小样本数据。

步骤1：准备数据 #

先将原始数据整理为数据框，方便后续处理：

# 构建数据框
drug_data <- data.frame(
  conc = c(0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30),  # 药物浓度（μM）
  response = c(97.01, 98.43, 98.29, 97.66, 96.51, 88.39, 38.54, 2.63)  # 响应值（%）
)

步骤2：拟合Sigmoid模型 #

安装并加载包后，用drm()函数拟合四参数逻辑斯蒂Sigmoid模型（最常用的药物抑制曲线模型）：

# 安装包（首次使用时运行）
install.packages("drc")
library(drc)

# 拟合四参数Sigmoid模型
sig_model <- drm(response ~ conc, data = drug_data, fct = LL.4())

该模型的方程形式为：
$$f(x) = d + \frac{c - d}{1 + e^{b(\log(x) - \log(e))}}$$
参数说明：

• c：响应渐近下限
• d：响应渐近上限
• b：曲线斜率（绝对值越大，曲线越陡）
• e：半数效应浓度（EC₅₀，即响应达到(c+d)/2时的药物浓度）

步骤3：获取模型方程及R²值 #

运行以下代码查看拟合参数，计算并输出R²：

# 查看模型拟合参数
summary(sig_model)

# 计算R平方值（drc包默认不输出，用残差偏差法计算）
r_squared <- 1 - deviance(sig_model)/deviance(lm(response ~ 1, data = drug_data))
cat("模型R平方值：", round(r_squared, 4), "\n")

步骤4：绘制美观的拟合曲线 #

用ggplot2绘制带对数刻度的清晰曲线（药物浓度通常用对数轴呈现更直观）：

# 安装并加载ggplot2（首次使用时运行）
install.packages("ggplot2")
library(ggplot2)

# 生成拟合曲线的预测数据
new_conc <- seq(min(drug_data$conc), max(drug_data$conc), length.out = 100)
predicted <- predict(sig_model, newdata = data.frame(conc = new_conc))
pred_df <- data.frame(conc = new_conc, response = predicted)

# 绘图
ggplot(drug_data, aes(x = conc, y = response)) +
  geom_point(size = 3, color = "#2c3e50", shape = 16) +  # 原始数据点
  geom_line(data = pred_df, aes(y = response), color = "#e74c3c", linewidth = 1.2) +  # 拟合曲线
  scale_x_log10(breaks = drug_data$conc, labels = drug_data$conc) +  # 对数刻度x轴
  labs(x = "药物浓度 (μM)", y = "响应值 (%)", title = "药物浓度-响应Sigmoid拟合曲线") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"))

2. 备选方案：用nls手动拟合基础Sigmoid模型 #

如果不想依赖专业包，可使用R内置的nls()函数拟合三参数Sigmoid模型，注意初始值设置要合理（否则可能拟合失败）：

# 拟合三参数Sigmoid模型：y = a / (1 + exp(-b(log(x) - c)))
sig_nls <- nls(response ~ a / (1 + exp(-b(log(conc) - c))), 
               data = drug_data,
               start = list(a = 100, b = -1, c = log(10)))  # 初始值需贴合数据趋势

# 查看拟合结果
summary(sig_nls)

# 计算R平方值
r_squared_nls <- 1 - deviance(sig_nls)/deviance(lm(response ~ 1, data = drug_data))
cat("手动拟合模型R平方值：", round(r_squared_nls, 4), "\n")

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