如何使用R拟合多组累计和曲线的通用最优拟合方程？

你的数据属于多组纵向测量的S型增长序列，所有曲线都满足从0渐近到1的约束，适合用非线性混合效应模型（NLME）完成拟合，既可以得到所有组通用的整体拟合方程，也可以保留每组的个体变异，具体实现步骤如下： #

步骤1：选择候选拟合模型 #

由于累积和取值区间为[0,1]，优先选择带上下界的生长类模型作为候选：

• 逻辑模型（Logistic）：适合对称S型曲线，公式为cumul = 1 / (1 + exp(a*(time_diff_to_complete - b)))，其中a为斜率参数，b为拐点位置
• 龚珀兹模型（Gompertz）：适合前半段增长偏慢的非对称S型曲线，公式为cumul = exp(-a*exp(-b*time_diff_to_complete))
• 修正指数模型：适合无明显拐点的单调增长曲线，公式为cumul = 1 - exp(a*time_diff_to_complete)

步骤2：代码实现 #

以逻辑模型为例，完整代码如下：

# 加载依赖包
library(nlme)
library(ggplot2)

# 导入样本数据
run_group <- c('A_group', 'A_group', 'A_group', 'A_group', 'A_group', 'A_group', 'A_group', 'A_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'B_group', 'C_group', 'C_group', 'C_group', 'C_group', 'C_group', 'C_group', 'C_group', 'D_group', 'D_group', 'D_group', 'D_group', 'D_group', 'D_group', 'D_group', 'D_group', 'D_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group', 'E_group')
cumul <- c(0.052631579, 0.263157895, 0.342105263, 0.710526316, 0.868421053, 0.894736842, 0.973684211, 1, 0.0078125, 0.015625, 0.0390625, 0.0546875, 0.0703125, 0.1015625, 0.1640625, 0.3203125, 0.4921875, 0.734375, 0.875, 0.96875, 0.9921875, 1, 0.073529412, 0.220588235, 0.323529412, 0.507352941, 0.727941176, 0.970588235, 1, 0.006134969, 0.055214724, 0.141104294, 0.190184049, 0.349693252, 0.595092025, 0.858895706, 0.969325153, 1, 0.005649718, 0.011299435, 0.016949153, 0.039548023, 0.073446328, 0.124293785, 0.299435028, 0.451977401, 0.559322034, 0.728813559, 0.81920904, 0.960451977, 1)
time_diff_to_complete <- c(-155, -140, -125, -110, -95, -80, -65, -50, -270, -210, -195, -180, -165, -150, -135, -120, -105, -90, -75, -60, -45, -30, -130, -115, -100, -85, -70, -55, -40, -175, -160, -130, -115, -100, -85, -70, -55, -40, -225, -210, -195, -180, -150, -135, -120, -105, -90, -75, -60, -45, -30)
sample_data <- data.frame(run_group, cumul, time_diff_to_complete, stringsAsFactors=FALSE)

# 先拟合单组模型获取初始参数，保证混合模型收敛
init_nls <- nls(cumul ~ 1 / (1 + exp(a * (time_diff_to_complete - b))), 
                data = sample_data,
                start = list(a = 0.01, b = -100))

# 拟合非线性混合效应模型，允许参数在不同运行组之间存在随机变异
nlme_fit <- nlme(cumul ~ 1 / (1 + exp(a * (time_diff_to_complete - b))),
                 data = sample_data,
                 fixed = a + b ~ 1, # 固定效应：所有组共享的总体参数
                 random = a + b ~ 1 | run_group, # 随机效应：每组独有的参数偏移
                 start = coef(init_nls))

# 查看模型结果，fixed部分即为整体拟合方程的参数
summary(nlme_fit)

# 输出拟合值对比
sample_data$fit <- predict(nlme_fit)
ggplot(sample_data, aes(x = time_diff_to_complete, y = cumul, group = run_group)) +
  geom_line(alpha = 0.5, color = "gray") + # 原始曲线
  geom_line(aes(y = fit), color = "red") + # 拟合曲线
  theme_minimal()

步骤3：最优模型选择 #

分别拟合所有候选模型后，通过AIC()函数计算各模型的AIC值，AIC最小的即为最优拟合模型。如果不需要参数化的数学公式，仅需要平均拟合曲线，也可以直接用loess非参数拟合：ggplot(sample_data, aes(x = time_diff_to_complete, y = cumul)) + geom_smooth(method = "loess")。

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