如何手动计算线性回归拟合值置信带并在R中绘图？代码无线条求助

我尝试在R语言内置的trees数据集的线性回归周围绘制误差带，但绘图中未显示误差带线条，附上代码及生成的图，请求协助解决！

原始代码 #

xdev <- log(Volume) - mean(log(Volume))
xdev
ydev <- Girth - mean(Girth)
ydev
b1 <- sum((xdev)*(ydev))/sum(xdev^2)
b0 <- mean(Girth) - mean(log(Volume))*b1
plot(log(Volume) ~ Girth)
abline(coef(lm(log(Volume) ~ Girth)), col=2)
(paste(b0, b1))
y.hat <- b0 + b1*log(Volume)
ss.explained <- sum((y.hat - mean(Girth))^2)
ss.unexplained <- sum((y.hat - Girth)^2)
stderr.b1 <- sqrt((ss.unexplained/29)/sum(xdev^2))
stderr.b1

std.yhat <- sqrt((ss.unexplained/29)*((1/31) + (xdev^2/sum(xdev^2))))
std.yhat
upp <-  y.hat + qt(0.95, 29)*std.yhat
low <-  y.hat - qt(0.95, 29)*std.yhat
upp
low
library(lattice)
plot(log(Volume) ~ Girth, data = trees)
abline(c(b0, b1), lty=1)
points(upp ~ Girth, data=trees, type="l", lty=2)
points(low ~ Girth, data=trees, type="l", lty=2)  

原始生成图 #

错误原因分析 #

1. 回归预测逻辑颠倒：你的模型是log(Volume) ~ Girth（以Girth为自变量预测log(Volume)），但计算y.hat时写成了y.hat = b0 + b1*log(Volume)，相当于用因变量预测自变量，完全不符合模型逻辑，导致误差带计算错误。
2. 线条绘制方式错误：
   • 用points()加type="l"画线条不恰当，应该用专门的lines()函数。
   • 原始数据的Girth未排序，直接绘制线条会导致点之间乱连，视觉上看起来线条“消失”。
3. 置信区间分位数错误：95%双侧置信区间应该用qt(0.975, 29)，而不是qt(0.95, 29)，后者是单侧分位数，会导致区间宽度错误。

修正后的代码 #

# 加载数据集
data(trees)

# 定义自变量和因变量
x <- trees$Girth
y <- log(trees$Volume)

# 计算回归系数（对应y ~ x，即log(Volume) ~ Girth）
xdev <- x - mean(x)
ydev <- y - mean(y)
b1 <- sum(xdev * ydev) / sum(xdev^2)
b0 <- mean(y) - b1 * mean(x)

# 计算预测值与残差平方和
y.hat <- b0 + b1 * x  # 正确：用Girth预测log(Volume)
ss.unexplained <- sum((y - y.hat)^2)
df <- nrow(trees) - 2  # 自由度为n-2，动态计算更灵活

# 计算预测值的标准误与置信区间
std_yhat <- sqrt((ss.unexplained/df) * (1/nrow(trees) + (xdev^2)/sum(xdev^2)))
upp <- y.hat + qt(0.975, df) * std_yhat  # 95%双侧置信区间
low <- y.hat - qt(0.975, df) * std_yhat

# 按自变量Girth排序，保证线条连续
sorted_idx <- order(x)
x_sorted <- x[sorted_idx]
y_sorted <- y[sorted_idx]
upp_sorted <- upp[sorted_idx]
low_sorted <- low[sorted_idx]
yhat_sorted <- y.hat[sorted_idx]

# 绘图
plot(y_sorted ~ x_sorted, 
     xlab = "Girth", 
     ylab = "log(Volume)", 
     main = "Linear Regression with 95% Confidence Bands")
abline(b0, b1, col = "red", lwd = 1)
lines(x_sorted, upp_sorted, lty = 2, col = "blue")
lines(x_sorted, low_sorted, lty = 2, col = "blue")

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