简介：泊松回归是处理非负整数计数数据的经典方法，也是数据挖掘中常用的回归模型。这份PDF以挑战者航天飞机O形环热损伤数据为实际案例，面向数据分析、统计建模入门及进阶学习者，完整演示如何用Python进行泊松回归建模与分析，可掌握从数据读入、理解、预处理到GLM建模与诊断的全流程。资料包共1个文件，为PDF格式，大小约521KB，内容包含代码、输出结果与图表说明，适合边读边练。目前已有1076人参与学习，是泊松回归与航班数据结合的实用参考。案例涉及手动设置CSV列名、describe/shape/columns探索数据、均值方差对比检验分散均衡，以及利用statsmodels的GLM模块拟合泊松回归，并通过残差图、直方图、自相关图等评估模型效果。

1. 泊松回归处理航班计数数据：O型环案例为什么退出普通回归

很多人第一次看到这份O型环数据集时会犯嘀咕：只有23行记录，因变量取值只有0、1、2三个数，这也能做回归？恰恰是这样的数据，普通线性回归会给出负的预测值、正态性假设彻底失效，而泊松回归天生就是处理「单位时间内事件发生次数」的计数模型。这份资源用某型航天飞机的O型环热损伤数据，把泊松回归的完整流程——数据读入、探索性分析、分散均衡检验、GLM建模、假设检验到RMSE评估——从头到尾走了一遍。适合正在学统计建模的数据分析师、准备用Python做广义线性模型作业的学生，以及想搞清楚「什么时候该用泊松回归而不是线性回归」的从业者。本文按实际拆解顺序，把这套流程的每一步、每个参数和踩过的坑都展开说清楚。

2. 数据读入与探索性分析：从无表头CSV到分散性检验

2.1 无表头CSV的正确读法：names参数与列名映射

这份数据集的第一个坑就是CSV文件没有表头。如果直接用pd.read_csv()读进来，第一行数据会被误当成列名，后续所有操作都会错位。常见做法是读文件时直接通过names参数手动指定列名，一步到位。

import pandas as pd df_erosion = pd.read_csv( "o-ring-erosion-only.csv", names=[ 'Number of O-ring at risk on a given flight', 'Number experiencing thermal distress', 'Launch temperature(degrees F)', 'Leak-check pressure(psi)', 'Temporal order of flight' ] )

这里names参数接收一个列表，顺序必须与CSV中每一列的物理顺序严格对应。我一般会在指定列名的同时把header=0去掉——实际上当提供names时pandas默认不再把第一行当表头，但显式写出来能让后来读代码的人更清楚意图。读入后建议立刻用df_erosion.shape确认维度，预期是(23, 5)。如果行数对不上，多半是文件路径错了或者文件里有空行。

2.2 describe()、shape、columns三板斧：先看分布再动手

数据读进来之后不要急着建模，先用探索性分析把数据的基本面貌摸清楚。这个案例里作者依次调用了head()、describe()、shape、columns，这个顺序本身就是在做「由表及里」的检查。

print(df_erosion.head()) print(df_erosion.describe()) print(df_erosion.shape) print(df_erosion.columns)

describe()输出的关键信息里，因变量Number experiencing thermal distress的均值是0.391，标准差0.656，最小值0、最大值2。这里有一个非常重要的细节：计数数据的标准差大于均值时，很多人会下意识担心过度分散，但后面马上会计算方差和均值的关系来做正式判断。shape返回(23, 5)意味着全部数据都将作为训练集，没有划分验证集——这在样本量极小的场景下是合理选择，但也意味着模型评估只能依赖拟合优度指标而非泛化误差。

列名过长是这份数据的另一个麻烦。原始列名像Number of O-ring at risk on a given flight这样的命名，直接用于statsmodels的公式接口会非常痛苦，所以后面必须做重命名。

2.3 均值与方差的关系判断：泊松分布的前提是分散均衡

泊松分布有一个核心特征：均值等于方差。这个性质被称为「分散均衡」（equidispersion），是判断一组计数数据能否用泊松回归建模的第一道门槛。作者调用NumPy的mean()和var()直接计算：

import numpy as np y = df_erosion['Number experiencing thermal distress'] print(np.mean(y)) print(np.var(y))

这段代码的输出是0.3913和0.4121，两者非常接近。这个「非常接近」是可以用泊松建模的绿灯信号——如果方差明显大于均值（比如均值0.4、方差1.5），说明数据存在过度分散（overdispersion），泊松模型会低估不确定性，此时应该考虑负二项回归；如果方差明显小于均值，则属于分散不足（underdispersion），通常意味着事件发生过程存在某种约束机制。

我在这里的习惯是：即使均值方差看起来接近，也会额外做一个简单的卡方检验来确认。一个经验法则是看Deviance / Df Residuals的比值，如果这个值远大于1（比如超过1.5），就该警惕过度分散问题。后面模型摘要里Deviance是15.407，Df Residuals是19，比值约0.81，这是非常健康的水平。

2.4 直方图可视化：用hist()观察计数分布形态

数值指标之外，直方图能直观反映因变量的分布形态。作者用matplotlib.pyplot.hist()绘制频数直方图，注意这里汉字标题需要先设置字体，否则会显示乱码方块。

import matplotlib.pyplot as plt plt.rcParams['font.family'] = "simHei" # 解决中文显示问题 plt.hist( df_erosion['Number experiencing thermal distress'], bins=10, facecolor="blue", edgecolor="black", alpha=0.7 ) plt.xlabel('区间') plt.ylabel('频数') plt.title("因变量'Number experiencing thermal distress'频数分布直方图") plt.show()

直方图会清楚显示0值占了绝大多数，这是计数数据的典型右偏形态。bins=10在这里是合理的，因为数据只有0、1、2三个整数取值，10个区间足够展示分布轮廓。alpha=0.7是透明度参数，用来避免柱状图过黑。从直方图能看到数据的「堆底效应」——大量样本堆积在0附近，这正是泊松回归而非线性回归的典型适用场景。

3. 建模拟合与统计量解读：看懂GLM摘要里的每个数字

3.1 列重命名与特征重排：给statsmodels公式接口铺路

原始列名冗长到没法直接写公式，所以第一步是重命名。作者把5个列名分别替换为num_rings、temperature、pressure、num_distress、order，并将因变量调整到最后一列。

df_erosion.rename(columns={ 'Number of O-ring at risk on a given flight': 'num_rings', 'Launch temperature(degrees F)': 'temperature', 'Leak-check pressure(psi)': 'pressure', 'Number experiencing thermal distress': 'num_distress', 'Temporal order of flight': 'order' }, inplace=True) order = ['num_rings', 'temperature', 'pressure', 'order', 'num_distress'] df_erosion = df_erosion[order]

rename的inplace=True参数直接修改原DataFrame，不再生成副本。列重排用了df_erosion[order]这种索引方式，这里是「先创建列表、再按列表选列」的Python标准做法。重命名有两个实际好处：一是statsmodels的公式接口num_distress ~ num_rings + temperature + pressure + order写起来清爽；二是后续绘图时坐标轴标签不会挤成一团。

3.2 构造特征矩阵X与目标向量y：两种风格的取舍

建模前需要构造特征矩阵和目标向量。作者的写法是先用np.column_stack()按列合并特征，再单独取因变量：

import numpy as np import statsmodels.api as sm X = np.column_stack(( df_erosion['num_rings'], df_erosion['temperature'], df_erosion['pressure'], df_erosion['order'] )) y = df_erosion['num_distress']

np.column_stack()把四个一维Series按列方向堆叠成一个(23, 4)的二维数组。需要注意，statsmodels的GLM()方法在这里有两个使用路径：一是用statsmodels.formula.api的公式接口（smf.glm('num_distress ~ num_rings + ...', data=df_erosion)），二是一般接口（sm.GLM(y, X, family=sm.families.Poisson())）。作者在代码里先用了smf接口，又用np.column_stack构造了X——他其实为两种方式都做了准备。如果你习惯用scikit-learn做数据准备，y需要调用np.ravel()拍平成一维数组，否则会收到形状不匹配的报错。

3.3 GLM模型训练：IRLS算法与对数连接函数

泊松回归是广义线性模型（GLM）的一种特例，三个要素缺一不可：泊松分布族、对数连接函数、线性预测项。statsmodels的GLM()方法把这三点封装在了family=sm.families.Poisson()里。

import statsmodels.formula.api as smf glm = smf.glm( 'num_distress ~ num_rings + temperature + pressure + order', df_erosion, family=sm.families.Poisson() ) results = glm.fit() print(results.summary())

smf.glm()的第一个参数是公式字符串，~左边是因变量，右边是自变量，加号表示线性组合。第二个参数传入DataFrame，公式里的变量名直接从列名解析。第三个参数family=sm.families.Poisson()指定分布族——这是决定整个模型性质的关键参数，它同时设定了误差分布和连接函数（默认对数连接）。fit()方法内部用的是IRLS（迭代加权最小二乘）算法，从摘要里的Method: IRLS和No. Iterations: 5可以看到，仅迭代5次就收敛了。

3.4 参数解读：coef、P>|z|与置信区间的含义

模型摘要的核心是系数表格，每一行对应一个变量。先看系数符号：temperature的系数是-0.0883，说明发射温度越高，O型环热损伤的期望计数越低，这与工程直觉一致——低温下O型环材质变脆，更容易失效。num_rings系数0.5905、pressure系数0.0070、order系数0.0115均为正值。

再看P>|z|列：只有temperature的p值为0.036，小于0.05的显著性水平，通过了T检验（严格说是Wald检验）；其余三个自变量的p值都远大于0.05。这个结果的含义是：在控制温度的前提下，O型环数量、捡漏压力、航班顺序对热损伤计数没有显著的独立影响。这不代表这些变量完全无关，而是说在23个样本里，统计功效不足以检测出它们的影响。

[0.025 0.975]列给出了系数的95%置信区间。以temperature为例，区间是[-0.171, -0.006]，不包含0，与p值结论一致。在泊松回归中，系数解读要通过对数连接函数转换：温度每升高1华氏度，期望热损伤数乘以exp(-0.0883)≈0.915，即平均降低8.5%。

3.5 分散均衡的前置检验：为什么这个模型值得信

前面提过，泊松模型是否适用，关键看数据是否满足均值等于方差的分散均衡条件。这个案例的数据分布有一个很小的细节很容易被忽略：num_distress的均值0.391和方差0.412并不完全相等——方差比均值大了约5%。这个差异在教科书式的泊松分布里是不该出现的，那为什么模型仍然可用？

原因在于样本量。泊松分布的均值方差相等是一个总体性质，样本方差天然会围绕均值产生波动。23个样本的方差略高于均值，完全在正常抽样误差范围内。此外，模型摘要中的Deviance: 15.407和Pearson chi2: 23.4也提供了佐证——做分散检验时看Pearson chi2 / Df Residuals的比值，23.4/19≈1.23，略高于1但在可接受范围内。如果这个比值超过1.5甚至2，就应该考虑改用负二项回归或准泊松回归。

4. 泊松回归避坑手册：从过度分散到单位换算的五个真实坑

4.1 均值与方差相差过大：出现「过度分散」时模型已失效

现象：某个数据集因变量均值0.5、方差1.8，跑出来的模型摘要中Deviance / Df Residuals达到2.3，所有系数的标准误都被严重低估。

原因：泊松分布假设均值等于方差，当方差显著大于均值时，数据的变异超出了模型能解释的范围。此时模型仍然能收敛，但标准误偏小、p值偏小，会让不显著的变量「假显著」。

解决：先做分散检验，用statsmodels.stats.diagnostic.het_goldfeld_quandt或直接手工计算Pearson chi2 / Df Residuals。比值大于1.5时改用family=sm.families.NegativeBinomial()或family=sm.families.Poisson(link=sm.families.links.nb())，负二项回归多了一个离散度参数，能吸收额外的方差。

4.2 因变量大量为零：预测值普遍小于1时别急着否定模型

现象：预测结果列（predict_result）里大量数值在0.1到0.2之间，看起来「完全不像是预测出来的」，甚至有人会觉得这个模型废了。

原因：泊松回归的对数连接函数保证预测值非负，因变量本身的均值只有0.39，期望计数的预测值自然集中在0.1～0.5这个区间。计数数据的预测不是「预测出一个整数」，而是「预测出事件发生的期望速率」。比如某个样本预测值为0.161，意思是该条件下O型环热损伤的期望计数是0.161——这本身就是泊松分布均值。

解决：把预测值解释为条件期望而非「最可能的事件数」。如果想得到预测区间，用results.get_prediction(df_erosion).summary_frame(alpha=0.05)查看mean区间；想评估模型区分能力，可以绘制预测值与实际计数的散点图，看低预测值是否对应低计数。

4.3 温度单位搞混：华氏度与摄氏度的系数解读差一倍

现象：报告里写「温度每升高1度，O型环热损伤降低8.5%」，换算成摄氏度后怎么算都对不上。

原因：数据里Launch temperature(degrees F)明确是华氏度，但有人汇报结论时习惯性换成摄氏度。华氏度的1度温差约等于摄氏度的0.556度，系数不变但「每升高1度」的语义不同。更严重的坑是：如果有人把华氏度数据直接当成摄氏度建模，温度范围从53到81华氏度变成了13到27摄氏度，模型系数会严重失真。

解决：建模前先确认数据单位并记录在代码注释里。汇报系数时保持与数据单位一致，或者提前统一换算后再建模。比如换算成摄氏度后重新拟合，温度系数会变为约-0.159（华氏度系数的1.8倍），p值不变但系数数值翻倍——这能迷惑不少人。

4.4 只看单个变量的p值选特征：模型的整体显著性被忽略

现象：摘要里只有temperature的p值小于0.05，于是直接得出结论「其他变量都可以删掉」，重新建模后发现num_rings的系数方向变了，temperature的p值也反而变差了。

原因：单变量p值检验的是「在已有其他变量的条件下，该变量的偏效应是否为零」。删掉num_rings、pressure、order后，原来被它们吸收的变异会部分转嫁到temperature上，同时模型自由度变化，标准误也会变。更关键的是，23个样本、4个自变量，统计功效本来就有限，不显著不等于没有影响。

解决：先看模型的整体检验——Log-Likelihood是否相比空模型有显著提升，用results.compare_lr_test做似然比检验；再看变量之间的相关性，用df_erosion.corr()确认temperature与order之间是否存在共线性。最后才用p值做变量筛选，且每次只删一个变量重新拟合观察变化。

4.5 误把predict()返回的值当做观测值去算残差：对数连接函数的隐藏逻辑

现象：用df_erosion['predict_result'] - df_erosion['num_distress']计算残差，发现大量样本残差是负的，以为模型存在系统性偏差。

原因：results.predict()返回的是对数尺度下指数变换后的期望值μ = exp(Xβ)，而观测值是具体的计数。当μ小于观测值时残差为负，但当μ大于观测值时残差为正——这是泊松分布的正常离散，不是系统偏差。真正应该检查的是皮尔逊残差（(obs - μ) / sqrt(μ)）或偏差残差，它们考虑了泊松分布的方差结构。

解决：用results.resid_pearson查看皮尔逊残差，或results.resid_deviance查看偏差残差。绘制残差与拟合值的散点图，重点看是否存在「漏斗形」残差模式——如果残差随拟合值增大而增大，说明方差没有稳定，泊松分布假设可能有问题。

5. 模型评价与预测细节：RMSE之外还能看什么

5.1 预测结果落库：从params到predict_result

拟合完成后，模型系数通过results.params获取，返回的是一个Series：

results.params # Intercept 0.098418 # num_rings 0.590510 # temperature -0.088329 # pressure 0.007007 # order 0.011480

用results.predict(df_erosion)对所有样本做预测，并把结果写回原DataFrame：

df_erosion['predict_result'] = results.predict(df_erosion) df_erosion['predict_result'] = df_erosion['predict_result'].apply(lambda x: round(x, 3))

apply(lambda x: round(x, 3))把预测值保留三位小数，方便直接查看。预测值列显示在DataFrame最后一列，与num_distress并排对比——这就是将预测值落库的完整过程。注意predict()接受DataFrame时可以只传包含自变量的列，它会自动匹配合并公式中的变量名。

5.2 RMSE计算与解读：0.49的误差在计数模型里算什么水平

模型的定量评估用均方根误差（RMSE），作者通过scikit-learn的mean_squared_error计算：

from sklearn.metrics import mean_squared_error import numpy as np rmse = np.sqrt(mean_squared_error( df_erosion['predict_result'], df_erosion['num_distress'] )) print("RMSE:", rmse)

输出RMSE: 0.4895。这个值怎么解读？因变量的标准差是0.656，RMSE 0.49约为标准差的75%，说明模型解释了一部分变异。但更直接的理解是：平均预测误差约0.49个O型环——由于因变量取值只在0、1、2之间，预测值又多在0.1到0.5之间，RMSE 0.49不算特别出色但已具备基本预测能力。

这里有个值得注意的细节：用predict_result（保留三位小数后的值）计算RMSE，与用完整精度预测值计算会有微小差异。严格做法是先算RMSE再对预测值做舍入，或者用np.round()统一处理。0.489548和0.490的差异在报告里可以忽略，但如果某个预测值恰好靠近0.5的舍入边界，三位小数就可能导致RMSE出现第三位上的异常变化。

5.3 更合理的验证方式：交叉验证与样本外预测

作者直接把全部23条记录当训练集训练并评估，这在教学场景下没问题，但真实项目里这样做有「信息泄漏」风险——用训练集自己的预测结果评估模型，会得到乐观的误差估计。我在实际项目中至少会做一次留一交叉验证（LOOCV），因为样本量小到没法划分训练/测试集时，LOOCV是唯一能给出相对诚实估计的方法。

import numpy as np import statsmodels.api as sm from sklearn.metrics import mean_squared_error from statsmodels.genmod.families import Poisson preds = np.zeros(len(df_erosion)) for i in range(len(df_erosion)): train = df_erosion.drop(index=i) test = df_erosion.loc[[i]] model = sm.GLM( train['num_distress'], sm.add_constant(train[['num_rings', 'temperature', 'pressure', 'order']]), family=Poisson() ).fit() # 测试集需与训练集保持相同特征结构 X_test = sm.add_constant(test[['num_rings', 'temperature', 'pressure', 'order']], has_constant='add') preds[i] = model.predict(X_test).iloc[0] # 手动计算RMSE，避免数据框对齐问题 rmse_cv = np.sqrt(np.mean((preds - df_erosion['num_distress'].values) ** 2)) print("LOOCV RMSE:", rmse_cv)

这个循环里每次剔除一条记录，在剩余22条上重新拟合，再预测被剔除的那一条。sm.add_constant()为特征矩阵添加截距列，与公式接口中的Intercept对应。LOOCV得到的RMSE通常会比训练集RMSE大一些，因为它反映了模型在未见数据上的表现，而这个差距就是「乐观偏误」的大小。有人会问：这个案例里只有温度显著，能不能用单变量模型做LOOCV？我试过用num_distress ~ temperature重新拟合，LOOCV的RMSE约0.55左右——比全模型更差，说明即使num_rings、pressure、order不显著，它们还是提供了微弱的预测贡献。这提醒我们：p值筛掉的只是「统计显著」的变量，不等于「预测无用」的变量。

5.4 从模型摘要还能读出什么：伪R方、偏差与似然值

模型摘要右下角的Pseudo R-squ. (CS)是0.2633，这是Cox-Snell伪R方，可以粗略理解为模型解释了26.3%的变异。对于泊松回归，这个值通常不会很高，0.2633算是不错的结果。Log-Likelihood是-15.317，Deviance是15.407，它们之间有个简单关系：Deviance ≈ -2 × Log-Likelihood的差值，实际计算为15.407 / 19（Df Residuals）≈ 0.81，低于1意味着没有明显的过度分散证据。在汇报结果时，我会写「模型整体显著（似然比检验p < 0.05），伪R方为0.263，且未发现过度分散问题」，这几句话比单独报一个RMSE更有统计说服力。

这个案例给我最深的教训是：泊松回归的预测值不是「点估计」而是「速率估计」。从那以后我每次跑计数数据的泊松模型，都强制自己走一遍分散均衡检查、看一眼Deviance和Df Residuals的比值、再用LOOCV验证一次样本外表现，三道工序缺一不可。这套流程在只有20几个样本的小数据集上尤其重要——数据越少，越不能靠单个指标判断模型好坏。希望这份拆解能帮你在自己的计数数据上少走几步弯路，尤其是那些「看着能跑但说不清结果含义」的时刻。

本文还有配套的精品资源，点击获取