如何在R语言中为对数回归曲线添加置信区间？

Absolutely! The key issue here is that unlike linear models fit with lm(), the predict() function for nls() models doesn’t have a built-in interval argument to generate confidence intervals directly. But we can fix this in base R with two reliable approaches:

方法1：基于参数置信区间与方差传播 #

这种方法利用模型参数的协方差矩阵，通过数学推导计算拟合曲线的置信区间，计算速度快且基于经典统计假设。

步骤与代码 #

首先，重新加载你的数据并拟合nls模型：

# 加载数据
Treatment<-c(15,12,6,3,15,12,6,9,9,15,6,9,3,3,12) 
CSA<-c(70.32255036, 141.7157414, 185.6793193, 382.1145078, 51.68417543, 99.38527948, 164.1865075, 160.8565834, 43.85808722, 46.99992099, 75.46593836, 116.4039474, 250.3521394, 290.2665555, 122.2660148) 
data<-data.frame(Treatment, CSA)

# 拟合对数nls模型
model <- nls(CSA ~ a*log(Treatment)+b, start = list(a = -141,b = 437), data=data)

接下来，计算预测值的方差与置信区间：

# 获取参数估计值与协方差矩阵
params <- coef(model)
vc <- vcov(model)

# 生成预测用的x序列
newx <- seq(min(data$Treatment), max(data$Treatment), length.out=1000)

# 构造设计矩阵：对应模型公式的log(Treatment)和常数项1
X <- cbind(log(newx), rep(1, length(newx)))

# 计算每个预测值的方差（方差传播公式）
pred_var <- diag(X %*% vc %*% t(X))

# 计算拟合均值与95%置信区间（t分布临界值，自由度为样本量-参数数）
pred_mean <- predict(model, newdata = list(Treatment=newx))
df <- nrow(data) - length(params)
crit_val <- qt(0.975, df)
lower_ci <- pred_mean - crit_val * sqrt(pred_var)
upper_ci <- pred_mean + crit_val * sqrt(pred_var)

最后，将拟合曲线与置信区间添加到你的图中：

# 绘制原始图形（沿用你的代码）
par(mfrow=c(1,1)) 
par(mar=c(2.5,2.5,1,1)) 
plot(data$Treatment,data$CSA,ylim=c(0,400),xlim=c(3,15),pch=21, 
     xaxt="n",yaxt="n",cex=0.6,xlab=NA,ylab=NA,bty="l") 
axis(side=1,tck=-0.02,at=seq(3,15,3),cex.axis=0.6, mgp=c(0,0.3,0)) 
axis(side=2,tck=-0.02,at=seq(0,400,100),cex.axis=0.6, las=2,mgp=c(0,.5,0)) 
ylab<-expression("Total cross-sectional area (cm"^{2}~")") 
xlab<-c("Treatment") 
mtext(xlab,side=1,line=1.5,cex=0.7) 
mtext(ylab,side=2,line=1.5,cex=0.7)

# 添加拟合曲线与置信区间
lines(newx, pred_mean, col="grey23",lwd=1.5)
lines(newx, lower_ci, lty = 'dashed', col = "grey36",lwd=1)
lines(newx, upper_ci, lty = 'dashed', col = 'grey36',lwd=1)

方法2：自助法（Bootstrapping） #

如果担心正态假设不成立，自助法是更稳健的选择——它通过重复抽样数据并重新拟合模型来估计置信区间。

代码实现 #

set.seed(123) # 设置随机种子保证结果可重复
n_boot <- 1000 # 自助抽样次数
boot_preds <- matrix(NA, nrow = n_boot, ncol = length(newx))

# 循环进行自助抽样与模型拟合
for(i in 1:n_boot){
  # 有放回抽样生成自助数据集
  boot_data <- data[sample(nrow(data), replace = TRUE), ]
  # 尝试拟合模型（避免拟合失败中断循环）
  try({
    boot_model <- nls(CSA ~ a*log(Treatment)+b, start = coef(model), data=boot_data)
    boot_preds[i, ] <- predict(boot_model, newdata = list(Treatment=newx))
  })
}

# 移除拟合失败的样本
boot_preds <- boot_preds[complete.cases(boot_preds), ]

# 计算每个x的95%置信区间（基于分位数）
lower_boot <- apply(boot_preds, 2, quantile, 0.025)
upper_boot <- apply(boot_preds, 2, quantile, 0.975)

# 绘制自助法置信区间（可以和方法1的区间对比）
lines(newx, lower_boot, lty = 'dotted', col = "blue",lwd=1)
lines(newx, upper_boot, lty = 'dotted', col = 'blue',lwd=1)

两种方法的对比 #

• 方法1：计算快，依赖参数服从正态分布的假设，适合模型拟合良好、数据符合假设的场景。
• 方法2：不依赖分布假设，更稳健，但计算量稍大，适合数据偏离正态或模型假设存疑的情况。

内容的提问来源于stack exchange，提问作者A.Benson