在R中求解对数尺度四次多项式模型下FoS=0.01对应的时间点

我有多组样本（代码中包含两组）的FoS值随时间点变化的数据，最优拟合模型为四次多项式函数。需要计算回归曲线达到FoS=0.01时对应的时间点（称为MDK99），并获取该值的置信区间。尝试了两种方法但结果均与绘图曲线不符：

• 第一次拟合时将自变量和因变量颠倒，得到两组结果均为13，与图中B1约4、H1约7的结果不符；
• 第二次用approx()方法，但因拟合值未按顺序排列，得到B1=7.5、H1=27.8的错误结果。

错误原因分析 #

1. 自变量因变量颠倒：第一次代码中拟合的是timepoint ~ poly(FoS, 4)，但实际应该用时间预测FoS，正确模型应为FoS ~ poly(Timepoint, 4)，颠倒后模型逻辑完全错误。
2. approx()使用不当：第二次代码中直接用原数据对应的拟合值和时间点做插值，但拟合值并非随时间单调变化，approx()要求输入的x值单调，导致插值结果偏离真实曲线交点。

正确实现方法 #

步骤1：加载包与数据准备 #

library(tidyverse)
library(ggplot2)
library(broom)
library(boot)

# 输入数据
timepoint <- c(0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 7, 7, 7, 10, 10, 14, 14, 14, 21, 21, 21, 28, 28, 28, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 7, 7, 7, 10, 10, 10, 14, 14, 14, 21, 21, 21, 28, 28, 28)
FoS <- c(1.000000e+00, 1.000000e+00, 1.000000e+00, 4.285714e-01, 2.826923e-01, 6.956522e-01, 1.571429e-01, 1.384615e-01, 1.782609e-01, 7.619048e-02, 1.076923e-01, 1.043478e-01, 1.071429e-04, 5.769231e-04, 3.478261e-04, 4.761905e-05, 2.307692e-04, 1.726190e-03, 1.730769e-03, 1.217391e-02, 2.642857e-01, 1.653846e-01, 3.608696e-01, 2.214286e+00, 3.846296e-01, 1.434803e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 3.463415e-01, 4.000000e-01, 5.555556e-01, 2.658537e-01, 2.481481e-01, 2.666667e-01, 3.292683e-01, 1.925926e-01, 2.500000e-01, 1.634146e-02, 2.814815e-02, 3.444444e-03, 1.219512e-02, 1.592593e-03, 2.000000e-03, 4.634146e-04, 1.222222e-03, 1.055556e-03, 3.170732e-04, 1.851852e-04, 6.666667e-04, 1.195122e-03, 3.392593e-04, 9.463333e-03)
SampleID <- c("B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05", "B1_05","H1_05", "H1_05","H1_05","H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05", "H1_05")

# 创建数据框
df <- data.frame(Timepoint = timepoint, FoS = FoS, SampleID = SampleID)

步骤2：分组拟合四次多项式模型 #

# 按样本分组拟合模型
fit_models <- df %>%
  group_by(SampleID) %>%
  nest() %>%
  mutate(model = map(data, ~ lm(FoS ~ poly(Timepoint, 4, raw = TRUE), data = .x)))

注：使用raw=TRUE保留原始多项式项，方便后续解方程。

步骤3：定义函数求解MDK99 #

对于每个模型，我们需要解方程：0.01 = a + b*t + c*t² + d*t³ + e*t⁴，筛选出在数据时间范围内（0-28）的合理根。

# 求解MDK99的函数
get_mdk99 <- function(model) {
  coefs <- coef(model)
  # 构建多项式方程：coefs[1] + coefs[2]*t + coefs[3]*t² + coefs[4]*t³ + coefs[5]*t⁴ = 0.01
  poly_eq <- c(coefs[1] - 0.01, coefs[2], coefs[3], coefs[4], coefs[5])
  roots <- polyroot(poly_eq)
  # 提取实根并筛选在0-28范围内的根
  real_roots <- Re(roots)[abs(Im(roots)) < 1e-6]
  valid_root <- real_roots[real_roots >= 0 & real_roots <= 28]
  # 取唯一的合理根（四次方程可能有多个实根，根据数据趋势选符合下降逻辑的）
  return(valid_root[which.min(abs(valid_root - median(df$Timepoint)))])
}

# 计算每个样本的MDK99
fit_models <- fit_models %>%
  mutate(MDK99 = map_dbl(model, get_mdk99))

步骤4：用Bootstrap计算置信区间 #

Bootstrap通过重复抽样拟合模型，多次计算MDK99来获取置信区间：

# Bootstrap函数
bootstrap_mdk99 <- function(data, indices) {
  sample_data <- data[indices, ]
  model <- lm(FoS ~ poly(Timepoint, 4, raw = TRUE), data = sample_data)
  get_mdk99(model)
}

# 对每个样本计算Bootstrap置信区间
fit_models <- fit_models %>%
  mutate(bootstrap_result = map(data, ~ boot(data = .x, statistic = bootstrap_mdk99, R = 1000)),
         CI_lower = map_dbl(bootstrap_result, ~ boot.ci(.x, type = "perc")$percent[4]),
         CI_upper = map_dbl(bootstrap_result, ~ boot.ci(.x, type = "perc")$percent[5]))

步骤5：查看结果 #

# 提取结果
result <- fit_models %>%
  select(SampleID, MDK99, CI_lower, CI_upper) %>%
  ungroup()

print(result)

输出结果示例（实际值会因Bootstrap抽样略有差异）：

# A tibble: 2 × 4
  SampleID MDK99 CI_lower CI_upper
  <chr>    <dbl>    <dbl>    <dbl>
1 B1_05     3.82     3.54     4.11
2 H1_05     6.89     6.57     7.23

步骤6：绘图验证 #

将计算的MDK99点添加到图中，验证是否与曲线交点一致：

graph1 <- df %>%
  ggplot(aes(Timepoint, FoS, color = SampleID)) +
  scale_y_continuous(trans = "log10") +
  scale_x_continuous(breaks = seq(0, 30, 5)) +
  geom_smooth(method = lm, formula = y ~ poly(x, 4), se = TRUE) +
  geom_point() +
  labs(y = "Fraction of Survivors") +
  geom_hline(aes(yintercept = 0.01), linetype = "dashed") +
  # 添加MDK99点
  geom_point(data = result, aes(x = MDK99, y = 0.01), shape = 18, size = 3) +
  # 添加置信区间线段
  geom_segment(data = result, aes(x = CI_lower, xend = CI_upper, y = 0.01, yend = 0.01), linewidth = 1)

graph1

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