线性模型的置信区间和预测区间及R语言绘制

线性模型的置信区间和预测区间及R语言绘制 by 小白鱼的生统笔记

线性模型的置信区间和预测区间及R语言绘制

置信区间（Confidence Interval）和预测区间（Prediction Interval）是统计学中两个相关但不同的概念，它们用于描述参数估计的不确定性和未来观测值的变异性。

置信区间是对一个参数估计的范围的估计，通常以一定置信水平表示（例如95%置信水平）。置信区间提供了一个区间，我们可以合理地认为真实参数值落在这个区间内的概率较高。例如，95%置信水平的置信区间表示，我们有95%的置信度认为真实参数值在这个区间内。置信区间通常用于对参数估计的不确定性进行量化。

预测区间是用于对未来观测值的变异性进行估计的概念。预测区间通常比置信区间更宽，因为预测区间考虑了除了参数估计的不确定性之外，还考虑了未来观测值的不确定性。预测区间的含义是在给定的置信水平下（例如95%置信水平），未来观测值有多大的概率会落在这个范围内。

总的来说，置信区间用于估计参数的不确定性，用于统计推断；而预测区间用于考虑未来观测值的不确定性。

本篇以普通最小二乘（OLS）回归为例，展示如何在R语言中获取模型的置信区间和预测区间，以及绘制带置信区间和预测区间的回归散点图，如下所示的这样（如上所述，由于预测区间通常比置信区间更宽，因此外侧阴影是预测区间，内侧阴影是置信区间，中间是拟合线）。

示例数据

以R自带的iris数据集为例进行作图演示。该数据集记录了150朵鸢尾花的花朵性状测量值，这些鸢尾花来自三种物种，分别为setosa（n=50）、versicolor（n=50）和virginica（n=50）。包含四种性状，分别为萼片长度（sepal length，cm）、萼片宽度（sepal width，cm）、花瓣长度（petal length，cm）和花瓣宽度（petal width，cm）。

#数据集，详情 ?iris
data(iris)
head(iris)

线性拟合并计算置信区间和预测区间

作为演示，以上述数据集中的花瓣长度与花瓣宽度的部分数据为例，拟合一个花瓣长度与花瓣宽度关系的普通最小二乘（OLS）回归，并分别计算其置信区间和预测区间。

#构建一个花瓣长度与宽度的关系的普通最小二乘（OLS）回归
fit <- lm(formula = Petal.Width~Petal.Length, data = iris)
summary(fit)
 
#计算置信区间，默认 95% 置信区间
confidence <- cbind(iris, predict(fit, iris, interval = 'confidence', level = 0.95))
head(confidence)
 
#计算预测区间，默认 95% 预测区间
prediction <- cbind(iris, predict(fit, iris, interval = 'prediction', level = 0.95))
head(prediction)

如本文开头所述，置信区间反映了预测均值的不确定性，预测区间反映了单个数值的不确定性。例如，本示例展示的是花瓣长度与花瓣宽度的关系，当花瓣长度在1.4时，花瓣宽度的预测均值约为0.22，95%置信区间约为(0.16, 0.27)，95%预测区间约为(-0.19, 0.63)。这意味着当花瓣长度在1.4时，我们有95%的置信度认为花瓣宽度的平均范围在0.16到0.27之间，95%的花瓣宽度在-0.19（负值在这个数据中无生物意义，仅为数学上的统计值）到0.63之间。

使用ggplot2绘制散点图并添加拟合线、置信区间和预测区间

最后，我们绘制花瓣长度与花瓣宽度的散点图，并将上述得到的回归线、置信区间和预测区间添加在图中。

#使用 ggplot2 绘制所有物种的花瓣长度与宽度的散点图的简单示例
library(ggplot2)
library(ggpmisc)
 
ggplot(data = iris, aes(x = Petal.Length, y = Petal.Width)) +
geom_point(aes(color = Species)) +  #不同颜色表示不同物种
geom_ribbon(data = prediction, aes(x = Petal.Length, ymin = lwr, ymax = upr), fill = 'gray', alpha = 0.3) +  #添加上述获得的 95% 预测区间的阴影
geom_line(data = prediction, aes(y = lwr), color = 'red', linetype = 2) +  #添加上述获得的 95% 预测区间下边界虚线
geom_line(data = prediction, aes(y = upr), color = 'red', linetype = 2) +  #添加上述获得的 95% 预测区间上边界虚线
geom_ribbon(data = confidence, aes(ymin = lwr, ymax = upr), fill = 'red', alpha = 0.3) +  #添加上述获得的  95% 置信区间的阴影
geom_line(data = confidence, aes(y = fit), color = 'red') +  #添加回归线（上述获得的预测均值）
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., stat(p.value.label), sep = '~`,`~')), formula = y~poly(x, 1), 
    parse = TRUE, label.x.npc = 'left', label.y.npc = 'top', size = 2.7) +  #添加回归式、R2
theme(panel.grid = element_blank(), 
    panel.background = element_rect(color = 'black', fill = 'white'), 
    axis.ticks = element_line(color = 'black', size = 0.5), 
    axis.text = element_text(color = 'black', size = 9)) +
labs(x = 'Petal.Length', y = 'Petal.Width', color = 'Species')