基于R语言nls的鲁棒单调递增三相线性模型拟合及区间计算

Got it, let's walk through fitting this monotonic increasing three-phase linear model with breakpoints B1 and B2 using R's tools, then calculate 95% confidence and prediction intervals. We'll cover both standard nonlinear least squares and robust fitting to handle outliers.

1. 整理并探索数据 #

首先计算实际的x值，然后可视化数据来辅助选择初始参数：

# 输入原始数据
y <- c(4.5,4.3,2.57,4.40,4.52,1.39,4.15,3.55,2.49,4.27,4.42,4.10,2.21,2.90,1.42,1.50,1.45,1.7,4.6,3.8,1.9)
x <- 1500 - c(320,419,650,340,400,800,300,570,720,480,425,460,675,600,850,920,975,1022,450,520,780)

# 查看x的范围，辅助断点猜测
range(x)
# 输出: [1]  478 1200

# 绘制散点图观察趋势
plot(x, y, pch=16, col="steelblue", main="x vs y 原始数据", xlab="x", ylab="y")

从图中能看到数据有噪声，但我们需要强制模型单调递增——所以三段的斜率都必须非负，且满足B1 < B2。

2. 定义三相单调递增模型 #

我们的模型是分段线性函数，通过约束斜率非负保证递增：

• 当x < B1：y = a + b1*x（斜率b1 ≥ 0）
• 当B1 ≤ x < B2：y = a + b1*B1 + b2*(x - B1)（斜率b2 ≥ 0）
• 当x ≥ B2：y = a + b1*B1 + b2*(B2 - B1) + b3*(x - B2)（斜率b3 ≥ 0）

用ifelse()把模型写成nls()兼容的公式：

three_phase_formula <- y ~ a + 
  ifelse(x < B1, b1*x, 
         ifelse(x < B2, b1*B1 + b2*(x - B1), 
                b1*B1 + b2*(B2 - B1) + b3*(x - B2)))

3. 选择初始参数值 #

非线性拟合对初始值很敏感，结合散点图的趋势：

• a：截距，从最小y值附近取1.5
• b1：第一段斜率，初始设为0.005
• B1：第一个断点，猜测为600
• b2：第二段斜率，初始设为0.003
• B2：第二个断点，猜测为900
• b3：第三段斜率，初始设为0.004

start_params <- list(a=1.5, b1=0.005, B1=600, b2=0.003, B2=900, b3=0.004)

4. 拟合模型（标准+鲁棒） #

4.1 标准非线性最小二乘拟合（nlsLM） #

Base R的nls()收敛性较差，我们用minpack.lm包的nlsLM()——它更稳定，还支持参数边界约束来保证单调性：

# 安装/加载所需包
if (!require(minpack.lm)) install.packages("minpack.lm")
library(minpack.lm)

# 带边界约束的拟合，保证单调性和合理断点范围
fit_std <- nlsLM(three_phase_formula,
                 start = start_params,
                 lower = c(a=-Inf, b1=0, B1=min(x)+10, b2=0, B2=600, b3=0),
                 upper = c(a=Inf, b1=Inf, B1=900, b2=Inf, B2=max(x)-10, b3=Inf),
                 data = data.frame(x, y))

# 查看拟合结果
summary(fit_std)

4.2 鲁棒非线性最小二乘拟合 #

如果数据存在异常值（从散点图看可能有），用robustbase包的nlrob()做鲁棒M估计，降低异常值的影响：

# 安装/加载所需包
if (!require(robustbase)) install.packages("robustbase")
library(robustbase)

# 带边界约束的鲁棒拟合
fit_robust <- nlrob(three_phase_formula,
                    start = start_params,
                    lower = c(a=-Inf, b1=0, B1=min(x)+10, b2=0, B2=600, b3=0),
                    upper = c(a=Inf, b1=Inf, B1=900, b2=Inf, B2=max(x)-10, b3=Inf),
                    data = data.frame(x, y),
                    method = "M",
                    control = nlrob.control(maxit=1000))

# 查看鲁棒拟合结果
summary(fit_robust)

5. 计算95%置信区间 #

非线性模型没有线性模型那样简洁的置信区间公式，我们用两种常用方法：

5.1 剖面似然置信区间（适用于标准拟合） #

confint()可以为nlsLM拟合计算剖面置信区间：

# 95%剖面置信区间
ci_profile <- confint(fit_std, level=0.95)
print(ci_profile)

5.2 Bootstrap置信区间（通用，支持标准/鲁棒拟合） #

Bootstrap方法更灵活，尤其适用于鲁棒模型，用boot包实现：

# 安装/加载所需包
if (!require(boot)) install.packages("boot")
library(boot)

# 参数Bootstrap函数
boot_param_fun <- function(data, indices) {
  d <- data[indices,]
  fit <- nlsLM(three_phase_formula,
               start = start_params,
               lower = c(a=-Inf, b1=0, B1=min(x)+10, b2=0, B2=600, b3=0),
               upper = c(a=Inf, b1=Inf, B1=900, b2=Inf, B2=max(x)-10, b3=Inf),
               data = d)
  coef(fit)
}

# 运行Bootstrap（1000次保证稳定性）
set.seed(123) # 保证结果可重复
boot_results <- boot(data=data.frame(x,y), statistic=boot_param_fun, R=1000)

# 提取95% Bootstrap置信区间
ci_boot <- t(apply(boot_results$t, 2, function(z) quantile(z, c(0.025, 0.975))))
rownames(ci_boot) <- names(coef(fit_std))
colnames(ci_boot) <- c("2.5%", "97.5%")
print(ci_boot)

如果要对鲁棒拟合计算Bootstrap置信区间，只需把boot_param_fun里的nlsLM()替换为nlrob()即可。

6. 计算95%预测区间 #

预测区间同时考虑模型不确定性和新观测的随机噪声，我们同样用Bootstrap实现：

# 预测Bootstrap函数
boot_pred_fun <- function(data, indices, new_x_vals) {
  d <- data[indices,]
  fit <- nlsLM(three_phase_formula,
               start = start_params,
               lower = c(a=-Inf, b1=0, B1=min(x)+10, b2=0, B2=600, b3=0),
               upper = c(a=Inf, b1=Inf, B1=900, b2=Inf, B2=max(x)-10, b3=Inf),
               data = d)
  # 预测值 + 抽样残差模拟新观测噪声
  preds <- predict(fit, newdata=data.frame(x=new_x_vals))
  resids <- residuals(fit)
  preds + sample(resids, length(new_x_vals), replace=TRUE)
}

# 定义要预测的x值范围
new_x <- seq(min(x), max(x), length.out=50)

# 运行预测Bootstrap
set.seed(123)
pred_boot_results <- boot(data=data.frame(x,y), statistic=boot_pred_fun, R=1000, new_x_vals=new_x)

# 提取95%预测区间
pred_ci <- t(apply(pred_boot_results$t, 2, function(z) quantile(z, c(0.025, 0.975))))
colnames(pred_ci) <- c("95%预测区间下限", "95%预测区间上限")

# 可视化拟合结果+预测区间
plot(x, y, pch=16, col="steelblue", main="三相模型拟合结果 + 95%预测区间", xlab="x", ylab="y")
lines(new_x, predict(fit_std, newdata=data.frame(x=new_x)), col="darkred", lwd=2)
lines(new_x, pred_ci[,1], col="gray", lty=2)
lines(new_x, pred_ci[,2], col="gray", lty=2)
legend("bottomright", legend=c("拟合模型", "95%预测区间"), col=c("darkred", "gray"), lty=c(1,2), lwd=c(2,1))

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