如何在R中实现BLI实验数据的全局贝叶斯拟合与局部拟合

先明确：你用的geom_smooth(method = "loess")不是目标拟合 #

Loess是局部加权回归，仅做曲线平滑，完全不涉及BLI的结合-解离动力学模型，和你需要的局部/全局动力学拟合不是一回事。我们需要基于BLI的核心动力学模型来实现。

第一步：整理数据为长格式 #

你的当前数据是宽格式（每个浓度一列响应值），这种格式不适合批量拟合，先转成长格式，方便后续分组处理：

# 模拟完整数据（替换成你真实的response2.5_y、response1.25_y数据）
mydata <- data.frame(
  time_x = rep(0:12, 3),
  conc_B = rep(c(5, 2.5, 1.25), each = 13),
  response_y = c(0, 0.1, 0.25, 0.6, 0.65, 0.65, 0.63, 0.62, 0.6, 0.59, 0.57, 0.55, 0.53,
                 0, 0.05, 0.12, 0.3, 0.32, 0.32, 0.31, 0.3, 0.29, 0.28, 0.27, 0.26, 0.25,
                 0, 0.02, 0.06, 0.15, 0.16, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.1, 0.09)
)

# 标记结合/解离阶段（x=3.9为解离开始时间）
mydata$phase <- ifelse(mydata$time_x <= 3.9, "association", "dissociation")

第二步：局部拟合（每个浓度单独拟合动力学参数） #

局部拟合是给每个浓度的结合、解离阶段分别拟合动力学模型，用基础R的nls（非线性最小二乘）即可实现：

1. BLI核心动力学模型 #

• 结合阶段：R(t) = R_max * (1 - exp(-k_on * conc * t))，其中k_on为结合速率常数，R_max为最大结合响应
• 解离阶段：R(t) = R0 * exp(-k_off * (t - t_dissoc))，其中k_off为解离速率常数，t_dissoc=3.9，R0为解离开始时的响应值

2. 批量执行局部拟合 #

library(dplyr)
library(purrr)

# 按浓度分组拟合
local_fits <- mydata %>%
  group_split(conc_B) %>%
  map(function(df) {
    # 分离结合/解离数据
    assoc_df <- df %>% filter(phase == "association")
    dissoc_df <- df %>% filter(phase == "dissociation")
    
    # 拟合结合阶段
    fit_assoc <- nls(response_y ~ R_max * (1 - exp(-k_on * conc_B * time_x)),
                     data = assoc_df,
                     start = list(R_max = max(assoc_df$response_y), k_on = 0.1))
    
    # 提取结合阶段的R_max作为解离阶段的初始R0
    R0 <- coef(fit_assoc)["R_max"]
    # 拟合解离阶段
    fit_dissoc <- nls(response_y ~ R0 * exp(-k_off * (time_x - 3.9)),
                      data = dissoc_df,
                      start = list(k_off = 0.01))
    
    # 生成拟合预测值
    pred_assoc <- predict(fit_assoc, newdata = assoc_df)
    pred_dissoc <- predict(fit_dissoc, newdata = dissoc_df)
    
    return(data.frame(
      time_x = c(assoc_df$time_x, dissoc_df$time_x),
      conc_B = df$conc_B[1],
      local_pred = c(pred_assoc, pred_dissoc)
    ))
  }) %>%
  bind_rows()

第三步：全局贝叶斯拟合（跨浓度共享核心参数） #

全局拟合要求所有浓度共享k_on和k_off参数，利用多浓度数据提升参数估计的可靠性，用brms包（基于Stan的贝叶斯建模工具）实现，适合新手：

1. 安装并加载包 #

install.packages("brms")
library(brms)

2. 定义全局贝叶斯模型 #

我们把结合和解离阶段合并成分段模型，共享k_on和k_off，每个浓度保留独立的R_max：

# 定义分段动力学模型公式
bli_model <- bf(
  response_y ~ ifelse(time_x <= 3.9, 
                      R_max[conc_B] * (1 - exp(-k_on * conc_B * time_x)),
                      R_max[conc_B] * exp(-k_off * (time_x - 3.9))),
  k_on ~ 1, k_off ~ 1, R_max ~ 0 + factor(conc_B),
  nl = TRUE
)

# 设置参数初始值（基于数据大致估计）
init_values <- list(
  list(k_on = 0.1, k_off = 0.01, R_max_factorconc_B5 = 0.65, 
       R_max_factorconc_B2.5 = 0.32, R_max_factorconc_B1.25 = 0.16)
)

# 运行贝叶斯拟合（迭代数设小以加快速度，正式分析可增大）
global_bayes_fit <- brm(
  formula = bli_model,
  data = mydata,
  init = init_values,
  chains = 2, iter = 2000,
  cores = 2  # 用2个核心加速
)

# 生成全局拟合的预测值（取后验均值）
global_pred <- predict(global_bayes_fit, newdata = mydata) %>%
  as.data.frame() %>%
  bind_cols(mydata %>% select(time_x, conc_B)) %>%
  rename(global_pred = Estimate)

第四步：整合数据与拟合曲线绘图 #

用ggplot2把原始数据、局部拟合、全局贝叶斯拟合整合到一张图中：

library(ggplot2)

ggplot(mydata, aes(x = time_x, y = response_y, color = factor(conc_B))) +
  # 原始数据曲线
  geom_line(size = 1) +
  # 局部拟合曲线（灰色虚线）
  geom_line(data = local_fits, aes(y = local_pred), color = "grey", linetype = "dashed", size = 1) +
  # 全局贝叶斯拟合曲线（黑色实线）
  geom_line(data = global_pred, aes(y = global_pred), color = "black", size = 1) +
  # 解离阶段分隔线
  geom_vline(xintercept = 3.9, linetype = "solid", color = "grey50", linewidth = 1.2) +
  # 样式与标签设置
  scale_color_manual(values = c("sienna1", "steelblue", "forestgreen"),
                     labels = c("1.25 nM", "2.5 nM", "5 nM"),
                     name = "B concentration") +
  labs(x = "Time [s]", y = "Wavelength shift [nm]") +
  theme_bw() +
  theme(
    axis.line.x.bottom = element_line(size = 1.5, color = 'black'),
    axis.line.y.left   = element_line(size = 1.5, color = 'black'),
    panel.grid.minor.x = element_blank(),
    panel.border       = element_blank(),
    axis.title.x = element_text(size = 20),
    axis.title.y.left = element_text(size = 20),
    axis.text = element_text(color = "black", size = 16),
    legend.title = element_text(size = 16),
    legend.text = element_text(size = 14)
  )

关键说明 #

• 局部拟合：每个浓度独立估计k_on、k_off、R_max，适合观察单个浓度的拟合效果，但参数不共享，精度较低。
• 全局贝叶斯拟合：所有浓度共享k_on和k_off，利用多浓度数据提升参数估计的可靠性，贝叶斯框架还能给出参数的不确定性区间。
• 如果真实数据中解离阶段的初始响应与结合阶段的R_max存在偏差，可在模型中加入偏移项调整。

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