使用ggplot绘制四参数逻辑回归模型并添加置信区间

我想用ggplot绘制四参数逻辑回归的结果，因变量y是鸟类数据集的多样性指标q1，自变量x是森林覆盖指标FC600。已经用mle2函数拟合好模型，用R基础绘图能成功画出原始数据和拟合曲线，但ggplot绘图始终失败，同时希望在图中添加标准误（SE）或95%置信区间。

数据 #

data$q1
 [1] 13.873635 13.377120 11.099290 11.481513  7.989731
 [6]  9.750019 12.737163 12.971196 12.413030 11.576204
[11] 13.909131 14.801793 11.270471 11.126155 13.620503
[16] 11.092484 11.651279 11.304653 13.736830

data$FC600
 [1]  42.70  23.83  12.08  20.46  29.91  27.69  25.87  30.60
 [9]  12.65  36.73  58.05  42.75  34.00  34.33  86.49  20.58
[17]  78.97  10.90 100.00

模型拟合代码 #

logip=function(p,mean,x){
  a=p[1]
  b=p[2]
  c=p[3]
  d=p[4]
  Riq1 = d+(a/(1+exp((b-(data$FC600))/c)))
  -sum(dnorm(x,mean=Riq1, log=TRUE))
}
parnames(logip)=c("a","b","c","d")

modTR.log6=mle2(minuslog=logip, start= c(a=max(data$q1),
                                         b=mean(data$FC600), 
                                         c=3,
                                         d=min(data$q1)),
                data=list(x=data$q1))

coefs=(cc <- coef(modTR.log6))
coefs=as.data.frame(coefs)

基础绘图代码 #

plot(data$FC600,data$q1, xlab="Forest cover", ylab="Richness",xlim=c(0,100))#original data
curve (coefs[4, 1]+(coefs[1, 1]/(1+exp((coefs[2, 1]-x)/coefs[3, 1]))), add=T,lwd=1.5)

解决方案 #

1. 修正模型拟合的硬编码问题

原模型函数直接调用data$FC600，导致后续无法灵活传入新的x序列，先重构模型：

# 重定义负对数似然函数，避免硬编码数据集
logip <- function(p, x, y) {
  a <- p[1]
  b <- p[2]
  c <- p[3]
  d <- p[4]
  # 四参数逻辑回归公式：y = d + a/(1 + exp((b - x)/c))
  pred <- d + (a / (1 + exp((b - x)/c)))
  -sum(dnorm(y, mean = pred, log = TRUE))
}
parnames(logip) <- c("a","b","c","d")

# 重新拟合模型，传入x和y参数
modTR.log6 <- mle2(minuslog = logip, 
                   start = c(a = max(data$q1),
                             b = mean(data$FC600), 
                             c = 3,
                             d = min(data$q1)),
                   data = list(x = data$FC600, y = data$q1))

2. 生成带置信区间的预测数据

ggplot要求所有绘图数据在数据框中，需先生成平滑x序列，再计算拟合值和95%置信区间：

# 生成0-100的平滑x序列
new_x <- seq(0, 100, length.out = 200)
new_data <- data.frame(FC600 = new_x)

# 定义四参数逻辑回归预测函数
predict_4pl <- function(coefs, x) {
  coefs[4] + (coefs[1]/(1 + exp((coefs[2] - x)/coefs[3])))
}

# 计算拟合值
fit_vals <- predict_4pl(coef(modTR.log6), new_x)

# 用delta方法计算标准误和95%置信区间
library(MASS)
vcov_mat <- vcov(modTR.log6)

# 定义梯度函数，用于delta方法
grad_fun <- function(x, coefs) {
  a <- coefs[1]
  b <- coefs[2]
  c <- coefs[3]
  d <- coefs[4]
  exp_term <- exp((b - x)/c)
  denom <- (1 + exp_term)^2
  
  da <- 1/(1 + exp_term)
  db <- (a * exp_term)/(c * denom)
  dc <- (a * (b - x) * exp_term)/(c^2 * denom)
  dd <- 1
  
  c(da, db, dc, dd)
}

# 计算每个x对应的标准误
se_vals <- sapply(new_x, function(x) {
  grad <- grad_fun(x, coef(modTR.log6))
  sqrt(t(grad) %*% vcov_mat %*% grad)
})

# 计算95%置信区间
lower_ci <- fit_vals - 1.96 * se_vals
upper_ci <- fit_vals + 1.96 * se_vals

# 合并为ggplot可用的数据框
pred_data <- data.frame(
  FC600 = new_x,
  fit = fit_vals,
  lower = lower_ci,
  upper = upper_ci
)

3. 用ggplot绘制图形

library(ggplot2)

ggplot() +
  # 原始数据散点
  geom_point(data = data, aes(x = FC600, y = q1), size = 2, alpha = 0.7) +
  # 拟合曲线
  geom_line(data = pred_data, aes(x = FC600, y = fit), linewidth = 1.5, color = "darkblue") +
  # 95%置信区间填充
  geom_ribbon(data = pred_data, aes(x = FC600, ymin = lower, ymax = upper), fill = "darkblue", alpha = 0.2) +
  # 坐标轴标签
  labs(x = "森林覆盖", y = "物种丰富度") +
  # 主题优化
  theme_bw()

关键说明

• 原模型硬编码data$FC600是ggplot绘图失败的核心原因，重构后可灵活传入新x序列；
• delta方法用于计算mle2拟合模型的预测置信区间，因为该模型无内置predict方法；
• ggplot要求所有绘图数据都在数据框中，因此必须提前生成包含拟合值和置信区间的预测数据框。

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