1. 从一次翻车的拟合说起为什么Statsmodels的OLS值得单独拎出来讲很多人第一次接触线性回归都是从scikit-learn的LinearRegression开始的。几行代码.fit()一下.coef_和.intercept_就出来了干净利落。但只要你稍微往统计推断的方向走一步——比如想知道某个自变量的系数在统计上到底显不显著、置信区间是多少、残差是否满足正态性假设——scikit-learn就会显得力不从心。它是个优秀的机器学习工具但它不是为统计推断设计的。这就是Statsmodels存在的意义。它把计量经济学和统计学里那一整套推断框架搬到了Python里OLSOrdinary Least Squares普通最小二乘法是其中最基础也最常用的入口。标题里说的简单一元线性回归指的是只有一个自变量加一个截距项的模型形式上是 $y \beta_0 \beta_1 x \varepsilon$。别看它简单一元回归是所有回归分析的母体多元回归里遇到的共线性、异方差、自相关等问题在一元场景下都能找到最纯粹的影子。这篇内容适合两类人一类是刚学完统计学课程、想把公式落到代码上的学生或转行者另一类是用惯了scikit-learn、现在需要输出一份带p值和置信区间的分析报告的数据从业者。我会把Statsmodels做一元OLS的完整链路拆开讲——从数据准备、模型拟合、结果解读到残差诊断和常见报错处理尽量把每个参数背后的为什么说清楚而不是只丢一段能跑的代码。需要提前说明的是下面涉及的所有数据、案例、场景均为虚构构造仅用于演示方法本身。2. 环境准备与数据构造别急着fit先把地基打对2.1 Statsmodels的安装与版本选择Statsmodels的安装本身没什么门槛pip install statsmodels就能搞定。但这里有个容易被忽略的坑Statsmodels对numpy和pandas的版本比较敏感尤其是numpy2.0发布之后早期版本的Statsmodels会出现AttributeError或者矩阵运算结果异常。我的建议是直接装较新的稳定版并且让pip自己去解析依赖pip install -U statsmodels pandas numpy如果你用的是conda环境conda install -c conda-forge statsmodels通常更省心因为conda-forge渠道对科学计算栈的版本兼容性做过测试。装完之后用下面这段代码验证一下环境是否正常import statsmodels.api as sm import statsmodels.formula.api as smf import pandas as pd import numpy as np print(sm.__version__) print(pd.__version__) print(np.__version__)这里要区分两个导入路径statsmodels.api提供的是面向数组的接口sm.OLSstatsmodels.formula.api提供的是面向公式的接口smf.ols。两者底层是同一套引擎区别在于你怎么描述模型。一元回归用哪个都行但公式接口在变量名有意义的时候可读性更好后面我会两种都演示。2.2 构造一份有故事的演示数据为了让后面的诊断环节有东西可讲我特意构造一份带一点非线性趋势和异方差的数据而不是那种完美线性的教科书数据。这样你在跑残差图的时候才能看到真实的样子np.random.seed(42) n 120 x np.linspace(1, 20, n) # 真实关系带轻微二次项 随x增大的噪声异方差 y 3.0 2.1 * x - 0.03 * x**2 np.random.normal(0, 0.5 0.15 * x, n) df pd.DataFrame({x: x, y: y}) print(df.head()) print(df.describe())这段数据里y和x大体是正相关的但噪声的标准差随着x增大而增大这就是典型的异方差苗头。用一元线性回归去拟合它斜率估计仍然是无偏的但标准误会偏进而影响p值和置信区间的可靠性。我把它留着就是为了在诊断章节里让你亲眼看到残差图的喇叭口。提示真实项目里拿到数据后先做df.describe()和缺失值检查df.isnull().sum()比上来就fit要重要得多。缺失值不处理Statsmodels默认会整行丢弃样本量悄悄变少你都不知道。2.3 一元回归里那个容易被忽视的截距项Statsmodels的sm.OLS默认不包含截距你必须手动加一列常数项这是和scikit-learn最大的行为差异之一。scikit-learn默认带截距Statsmodels默认不带。很多人第一次用sm.OLS(y, x)跑出来截距是0还以为是数据问题其实是接口设计如此。正确做法是用sm.add_constant()X sm.add_constant(df[x]) print(X.head())加完之后X会多出一列名为const的全1列。为什么要有截距从几何上讲没有截距的回归线被强制过原点这只有在业务上明确知道x为0时y必然为0时才合理。绝大多数场景下强制过原点会引入系统性偏差。所以除非你有强理论依据否则永远加上const。3. 两种拟合姿势数组接口与公式接口的取舍3.1 数组接口sm.OLS控制力强适合程序化流程数组接口的写法是这样的X sm.add_constant(df[x]) model sm.OLS(df[y], X) results model.fit() print(results.summary())sm.OLS(endog, exog)里endog是因变量内生变量exog是自变量矩阵外生变量。注意exog必须是二维的如果你直接传一维的df[x]会报ValueError: endog and exog matrices are different sizes或者维度相关的错误。add_constant返回的是DataFrame天然是二维的所以没问题。数组接口的好处是变量可以是numpy数组、可以是动态拼接的矩阵适合写在循环里批量跑模型或者自变量是程序计算出来的中间结果。缺点是结果的可读性差一些summary里显示的是x1、x2这种列名变量一多就对不上号。3.2 公式接口smf.ols可读性优先贴近统计软件习惯公式接口长这样model smf.ols(y ~ x, datadf) results model.fit() print(results.summary())y ~ x这种写法直接借用了R语言的公式语法~左边是因变量右边是自变量。公式接口会自动帮你加截距这是它和数组接口的另一个重要区别。如果你确实想拟合过原点的模型要显式写成y ~ x - 1或者y ~ 0 x。两种接口怎么选我的经验是探索性分析、变量少、需要频繁看summary的时候用公式接口生产代码、变量动态生成、需要和其他矩阵运算混用的时候用数组接口。两者拟合出来的系数和统计量完全一致不存在哪个更准的问题。3.3 拟合结果对象里到底装了什么results model.fit()返回的是一个RegressionResults对象它是个信息富矿。除了summary()常用的属性包括属性/方法含义典型用途results.params系数含截距取斜率、截距做预测results.bse系数标准误手工算t值、置信区间results.tvaluest统计量判断显著性results.pvaluesp值显著性检验results.conf_int()置信区间报告区间估计results.rsquaredR方拟合优度results.resid残差残差诊断results.fittedvalues拟合值画拟合线results.mse_resid残差均方估计误差方差这些属性在写自动化报告的时候特别有用。比如你想把系数和p值导出成表格直接results.params和results.pvalues拼一下就行不用去summary的文本里正则匹配。4. 读懂summary那张表每个数字背后的统计含义4.1 系数、标准误与t检验的三角关系summary()里最核心的是中间那张系数表四列分别是coef、std err、t、P|t|。它们的关系是$$t \frac{\hat{\beta} - 0}{SE(\hat{\beta})}$$也就是说t值等于系数估计值除以它的标准误。这个t值检验的原假设是该系数等于0。p值就是在这个t分布下观测到比当前t值更极端情况的概率。p值小于0.05我们通常就说这个系数统计显著。但这里有个新手常犯的误解p值小不代表效应大。p值同时受效应大小和样本量影响。样本量足够大时一个微不足道的系数也能得到极小的p值。所以看结果要系数大小和p值一起看再结合业务判断这个效应有没有实际意义。4.2 R方、调整R方与F检验一元场景下的特殊之处一元回归里R方就是自变量和因变量的皮尔逊相关系数的平方。它表示因变量的变异中能被自变量解释的比例。R方0.8意味着80%的变异被解释了听起来不错但一元回归里R方高不代表模型就对——如果真实关系是非线性的线性模型也可能给出一个不低的R方但残差图会暴露问题。调整R方在多元回归里才有意义它惩罚了变量个数一元回归里调整R方和R方差别很小。F检验在一元回归里等价于对斜率的t检验因为只有一个自变量F统计量等于t统计量的平方。所以一元场景下你看到F检验显著和斜率t检验显著是一回事不用纠结。4.3 那几个常被忽略的诊断统计量summary底部还有几行Omnibus、Prob(Omnibus)、Durbin-Watson、Jarque-Bera、Skew、Kurtosis、Cond. No.。这些是残差诊断的快速指标Omnibus / Jarque-Bera检验残差是否正态。p值小于0.05说明残差偏离正态会影响小样本下的推断。Durbin-Watson检验残差自相关取值在0到4之间接近2说明无自相关。时间序列数据要特别关注这个。Skew / Kurtosis残差的偏度和峰度正态分布下分别接近0和3。Cond. No.条件数衡量自变量矩阵的病态程度。一元回归里一般不会有大问题多元回归里条件数过大意味着共线性严重。这些指标是体检报告不是判决书。它们提示你哪里可能有问题但要不要处理、怎么处理还得结合具体场景。5. 残差诊断实战从喇叭口到模型修正5.1 画残差图诊断的第一步永远是看图统计量是数字图是直觉。跑完回归后我习惯先画两张图残差对拟合值、残差对自变量。import matplotlib.pyplot as plt fitted results.fittedvalues resid results.resid fig, axes plt.subplots(1, 2, figsize(12, 4)) axes[0].scatter(fitted, resid, alpha0.6) axes[0].axhline(0, colorred, linestyle--) axes[0].set_xlabel(Fitted) axes[0].set_ylabel(Residual) axes[0].set_title(Residual vs Fitted) axes[1].scatter(df[x], resid, alpha0.6) axes[1].axhline(0, colorred, linestyle--) axes[1].set_xlabel(x) axes[1].set_ylabel(Residual) axes[1].set_title(Residual vs x) plt.tight_layout() plt.show()用我前面构造的数据跑你会看到残差随着x增大而张开形成喇叭口。这就是异方差的视觉信号。理想情况下残差应该像一条水平带均匀分布在0线两侧不呈现任何系统性形状。5.2 异方差、非线性、异常点三种典型病的识别残差图能暴露的问题主要有三类异方差残差的离散程度随拟合值或自变量变化。除了看图还可以用Breusch-Pagan检验或White检验做正式检验。Statsmodels里可以这样跑from statsmodels.stats.diagnostic import het_breuschpagan bp_test het_breuschpagan(results.resid, results.model.exog) labels [LM Statistic, LM-Test p-value, F-Statistic, F-Test p-value] print(dict(zip(labels, bp_test)))如果p值小于0.05就拒绝同方差假设说明存在异方差。处理办法包括对因变量做对数变换、使用稳健标准误results.get_robustcov_results(cov_typeHC3)、或者改用加权最小二乘。非线性残差图呈现U形或倒U形说明线性模型没抓住真实关系。解决办法是加二次项、做变量变换或者换非线性模型。我构造的数据里那个-0.03 * x**2项就是故意留下的非线性成分。异常点残差绝对值特别大的点。可以用标准化残差或学生化残差来识别绝对值大于3的通常值得关注。但要注意异常点不一定要删——它可能是真实的数据也可能是录入错误得回去查数据来源。5.3 用稳健标准误给推断上保险如果确认存在异方差但又不想改模型形式最省事的办法是用稳健标准误。Statsmodels支持多种协方差估计results_robust results.get_robustcov_results(cov_typeHC3) print(results_robust.summary())HC3是异方差稳健标准误里比较保守的一种小样本下表现较好。对比一下普通标准误和HC3标准误你会发现系数的点估计没变但标准误、t值、p值都变了。这就是为什么异方差下不能直接信普通summary的原因。注意稳健标准误只修正推断不改变系数估计。如果你的目的是预测异方差对点预测影响不大但如果目的是做假设检验、报告置信区间稳健标准误几乎是必须的。6. 预测与结果导出把模型用起来6.1 用get_prediction做区间预测拟合完模型预测新数据是常见需求。Statsmodels的get_prediction不仅能给点预测还能给预测区间和置信区间new_x pd.DataFrame({x: [5, 10, 15]}) new_X sm.add_constant(new_x) pred results.get_prediction(new_X) pred_summary pred.summary_frame(alpha0.05) print(pred_summary)summary_frame返回的列里mean是拟合值mean_ci_lower和mean_ci_upper是均值的置信区间obs_ci_lower和obs_ci_upper是单次观测的预测区间。这两个区间容易混淆置信区间描述的是平均响应的不确定性预测区间描述的是单个新观测的不确定性后者总是更宽因为它还包含了噪声的方差。6.2 把结果整理成可交付的表格实际工作里分析结果往往要交给别人看。把summary()的文本直接贴过去不够专业我一般会整理成结构化表格import pandas as pd coef_table pd.DataFrame({ coef: results.params, std_err: results.bse, t: results.tvalues, p_value: results.pvalues, ci_lower: results.conf_int()[0], ci_upper: results.conf_int()[1], }) print(coef_table.round(4))这样导出的表格可以直接进报告也可以存成CSV。如果要做多个模型的对比可以把每个模型的params和pvalues拼成一张宽表一眼看出系数随模型设定的变化。6.3 一元回归的边界什么时候该果断换模型一元线性回归再简单也有它的适用边界。以下几种情况我会果断放弃一元线性模型散点图明显呈现曲线关系且业务上无法用线性近似因变量是二分类或计数型该用逻辑回归或泊松回归残差诊断严重违反假设且变换后仍无改善存在明显的遗漏变量一元模型只是看起来显著。一元回归最大的价值在于它是理解回归思想的起点而不是终点。把一元场景下的每个环节吃透多元回归、广义线性模型的学习曲线会平缓很多。7. 那些年踩过的报错与坑7.1 维度不匹配与缺失值引发的报错最常见的报错是ValueError: endog and exog matrices are different sizes。原因通常是exog传了一维数组或者y和X的索引没对齐。用add_constant之后如果原DataFrame有缺失值add_constant不会自动对齐得先dropna()。另一个隐蔽的坑是df[x]和df[y]如果来自不同的DataFrame索引不一致sm.OLS会按位置对齐而不是按索引对齐结果就错了。稳妥做法是始终从同一个DataFrame里取列。7.2 公式接口里的变量名陷阱公式接口smf.ols(y ~ x, datadf)里x和y必须是df里真实存在的列名。如果列名带空格或特殊字符要用Q()包裹比如y ~ Q(my var)。另外公式接口默认会丢弃含有缺失值的行但不会告诉你丢了多少行样本量对不上时记得查一下。7.3 关于显著的执念最后一个坑是心态上的。很多人跑出一元回归看到p值小于0.05就松了口气觉得模型成立了。但统计显著只是最低门槛残差诊断、效应大小、业务合理性一个都不能少。我见过太多R方0.99但残差图惨不忍睹的模型也见过p值0.06但业务价值巨大的关系。工具是工具判断还得靠人。把Statsmodels的一元OLS跑通不难难的是养成拟合完先看残差、报告前先想假设的习惯。这个习惯一旦建立后面学多元回归、面板数据、时间序列都会顺很多。