如何用贝叶斯统计拟合S型曲线并提取指定y对应x（brms改写）

我有多棵树木的多日测量数据集，测量值范围0到1。以Day为x轴、Value为y轴绘图时，每棵树的测量值呈S型曲线。目前用非线性最小二乘法（nls）拟合曲线，提取y=0.5对应的x值（日期）及该点的曲线梯度，后续将这些x值作为空间模型的响应变量。为考虑估计的不确定性，希望改用brms实现贝叶斯方法，作为贝叶斯统计新手找不到易懂教程，请求改写代码。

现有nls代码 #

results<-data.frame()
# 遍历每棵树
trees<-unique(all$TreeID)
for (tree in trees) {
  # 提取单棵树的数据
  subset_data<-filter(all, TreeID==tree)
  # 拟合S型逻辑曲线
  model<-nls(Value~SSlogis(Day, Asym, xmid, scal), data=subset_data)
  # 提取模型参数
  params<-coef(model)
  # 计算y=0.5对应的x值（萌芽日）
  bbday<-params["xmid"]+params["scal"]*log(0.5/(params["Asym"]-0.5))
  # 计算y=0.5处的曲线梯度
  bbgrad<-params["Asym"]*exp((params["xmid"]-bbday)/params["scal"])/(params["scal"]*(exp((params["xmid"]-bbday)/params["scal"])+1)^2) 
  # 存储结果
  results<-rbind(results, data.frame(TreeID=tree, BudBurstDay=bbday, BudBurstGrad=bbgrad))
}

示例数据 #

TreeID<-c(14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511, 14511)
Day<-c(84, 87, 89, 93, 96, 98, 103, 105, 108, 111, 114, 117, 120, 123, 127, 129, 132, 134, 138, 141, 144, 147, 150, 155, 158)
Value<-c(0.04166667, 0.06250000, 0.16666667, 0.26666667, 0.26666667, 0.30000000, 0.36666667, 0.42857143, 0.16666667, 0.43333333, 0.50000000, 0.41666667, 0.61111111, 0.77083333, 0.89583333, 1.00000000, 1.00000000, 1.00000000, NA, 1.00000000, 1.00000000, 1.00000000, 1.00000000, NA, 1.00000000)

brms改写方案 #

核心思路 #

1. 手动定义与SSlogis等价的非线性逻辑增长公式
2. 用分层建模语法为每棵树单独估计参数，替代循环拟合
3. 设置弱信息先验适配数据特性
4. 从模型后验样本中计算目标特征值，并量化不确定性

完整代码实现 #

library(brms)
library(dplyr)

# 处理示例数据（实际使用时替换为你的all数据集）
all_data <- data.frame(TreeID, Day, Value) %>% filter(!is.na(Value))

# 定义非线性逻辑增长模型公式
logistic_formula <- bf(
  Value ~ Asym / (1 + exp((xmid - Day)/scal)),
  Asym ~ 0 + (1|TreeID),  # 每棵树独立估计渐近值Asym
  xmid ~ 0 + (1|TreeID),  # 每棵树独立估计中点xmid
  scal ~ 0 + (1|TreeID),  # 每棵树独立估计缩放参数scal
  nl = TRUE  # 声明为非线性模型
)

# 设置弱信息先验
priors <- c(
  prior(normal(1, 0.5), nlpar = "Asym"),  # Asym接近1（Value上限为1）
  prior(normal(120, 30), nlpar = "xmid"), # xmid落在观测日期的中间范围
  prior(lognormal(0, 1), nlpar = "scal")  # scal为正数，用对数正态先验
)

# 拟合贝叶斯模型
brms_model <- brm(
  formula = logistic_formula,
  data = all_data,
  prior = priors,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  na.action = na.exclude
)

# 提取后验样本并计算目标特征值
posterior_samples <- posterior_samples(brms_model)
tree_ids <- unique(all_data$TreeID)
results_brms <- data.frame()

for (tree in tree_ids) {
  # 提取当前树的参数后验列
  asym_col <- paste0("Asym_", "TreeID[", tree, "]")
  xmid_col <- paste0("xmid_", "TreeID[", tree, "]")
  scal_col <- paste0("scal_", "TreeID[", tree, "]")
  
  asym_post <- posterior_samples[[asym_col]]
  xmid_post <- posterior_samples[[xmid_col]]
  scal_post <- posterior_samples[[scal_col]]
  
  # 计算每个后验样本对应的萌芽日和梯度
  bbday_post <- xmid_post + scal_post * log(0.5 / (asym_post - 0.5))
  bbgrad_post <- asym_post * exp((xmid_post - bbday_post)/scal_post) / (scal_post * (exp((xmid_post - bbday_post)/scal_post) + 1)^2)
  
  # 计算后验统计量（均值+95%可信区间）
  tree_result <- data.frame(
    TreeID = tree,
    BudBurstDay_mean = mean(bbday_post),
    BudBurstDay_lower = quantile(bbday_post, 0.025),
    BudBurstDay_upper = quantile(bbday_post, 0.975),
    BudBurstGrad_mean = mean(bbgrad_post),
    BudBurstGrad_lower = quantile(bbgrad_post, 0.025),
    BudBurstGrad_upper = quantile(bbgrad_post, 0.975)
  )
  
  results_brms <- rbind(results_brms, tree_result)
}

# 查看最终结果
print(results_brms)

关键说明 #

• 分层建模：通过(1|TreeID)实现单树参数独立估计，无需手动循环，效率更高
• 不确定性量化：结果包含均值和95%可信区间，直接反映估计的波动范围，满足后续空间模型的需求
• 模型诊断：拟合完成后可通过plot(brms_model)检查收敛性，pp_check(brms_model)验证拟合效果

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