分类变量与连续变量相关性分析:方法选择、效应量与Python实战
1. 相关性分析真正在回答的问题
相关性分析这个词被用得太泛了。做数据的人嘴里说的"相关",有时候是统计学意义上的相关系数,有时候只是"这两列看起来有关系",还有时候是因果关系的委婉说法。我见过不少分析报告,拿一个 Pearson 系数就下结论说"A 越高 B 越大",但实际上 A 是个性别字段被编码成 0/1,B 是用户消费金额——这种"相关系数"算出来的数字在数学上没毛病,在业务解释上基本是一场灾难。
问题就出在变量类型上。相关性分析的本质,是衡量两个变量"共同变化"的程度和方向,而"共同变化"这四个字对不同类型的变量有不同的数学定义。连续变量之间的共同变化可以用协方差刻画,分类变量之间的共同变化只能靠频数分布去比对,而分类变量和连续变量之间的共同变化,既没有协方差可言、也没有频数表可列,必须借助"组间差异"这个中间桥梁来翻译。这就是为什么分类变量与连续变量的相关性分析是最容易翻车的一块:它不像连续对连续那样有唯一的标准答案,二分类、多分类、有序分类各自对应不同的方法,选错了数字还能跑出来,只是结论全错。
这篇东西就是想把这块讲透。适合谁看?如果你做过用户画像、AB 实验分析、问卷数据清洗、风控建模前的特征筛查,或者你只是单纯想知道"性别和收入到底有没有关系、关系多大"该怎么严谨地回答,那下面的内容应该能直接用上。我会按连续对连续、分类对分类、分类对连续三条线走,重点压在分类对连续上,因为那是方法最多、坑也最多的地方,最后给一套能跑的代码和一张排查表。
1.1 变量类型决定了你能问什么问题
先说一个容易被忽略的前提:变量类型不是数据自带的属性,而是分析者根据业务含义赋予的。同样一列"用户等级",取值 1 到 5,你可以把它当有序分类(等级高低有意义,但 1 到 2 的距离不等于 2 到 3 的距离),也可以当连续变量(假设间隔均匀),甚至当无序分类(只关心属于哪一档)。这三种身份会导向三种完全不同的分析方法,结论也可能不一样。
这个认识很重要。我见过太多人拿到数据第一件事就是df.corr(),然后对着一张全是数字的热力图发呆,完全没意识到里面有些列是标签编码后的省市区、有些是独热编码的残渣、有些是 0/1 的布尔值。corr()默认跑 Pearson,它对输入的唯一要求是"能转成数值",转完之后它不管语义,一律按等距变量处理。分类变量被编码成整数之后,Pearson 算出来的东西没有统计学意义——它衡量的是"编码值"的线性关系,而不是"类别"与"数值"的关系。省市区编码里"北京=1、上海=2",算出来的相关系数你跟谁都没法解释。
所以相关性分析的第一步永远是分类:把每一列贴上"连续""有序分类""无序分类"三个标签之一,然后按标签组合去查方法表。方法表长这样:
| 变量 A 类型 | 变量 B 类型 | 推荐方法 | 效应量指标 | 对称性 |
|---|---|---|---|---|
| 连续 | 连续 | Pearson / Spearman | r / ρ | 对称 |
| 无序分类 | 无序分类 | 卡方检验 / Fisher 精确 | Cramér's V、φ | 对称 |
| 二分类 | 连续 | 点二列相关 / 独立样本 t 检验 | r_pb、Cohen's d | 对称 |
| 多分类 | 连续 | 单因素方差分析 / Kruskal-Wallis | η²、ε² | 不对称 |
| 有序分类 | 连续 | Spearman / Jonckheere-Terpstra | ρ、趋势统计量 | 对称 |
| 有序分类 | 有序分类 | Spearman / Kendall's τ-b | ρ、τ | 对称 |
表里"对称性"那一列值得单独说一句。Pearson、Spearman、Cramér's V 都是对称的:A 和 B 换个位置,数字不变。但方差分析导出的 η² 不是——它算的是"分类变量能解释连续变量多少变异",反过来问"连续变量能解释分类变量多少"就没法用 η² 表达了。这种不对称在写报告时要注意措辞,别写成"两者互相影响",那已经滑向因果论了。
1.2 分类与连续之间的相关为什么总被算成"假的"
我总结过三类高频的假相关,几乎每个项目都会撞上一个。
第一类是编码当连续。把无序分类做 LabelEncoder 之后直接进corr()。比如城市编码、渠道 ID、商品类目 ID,这些都只是编号,编号大小没有任何顺序含义。算出来的 Pearson 系数纯属数字巧合,样本量一大还能"显著",特别能唬人。判断方法很简单:看这一列的取值是否等距、是否有序、0 值是否代表"没有"。任何一个不成立,就别用 Pearson。
第二类是把连续变量硬切成组再算相关。这属于信息损失型错误。有人为了"方便看",把年龄切成青年/中年/老年三组,然后算这三组和收入的关系。问题在于,切点是你自己定的,换个切点结论就可能变,而且分组会大幅损失变异信息,导致检验功效下降。原本连续对连续也许能看出 0.2 的相关,切完组之后组间差异不显著,你就误判成"没关系"。这个错误的方向是保守的,但保守也意味着漏掉真信号。
第三类是忽略分布形态。多分类对连续,标准做法是方差分析,而方差分析对正态性和方差齐性有要求,同时对极端值极其敏感。一个组里混进两个异常大的值,均值被拉高,F 值就上去了,你得到"组间差异显著",回头查数据发现是两个录错的 999999。所以做这类分析之前,箱线图必须画,Levene 检验必须跑,扭不过去的就换 Welch ANOVA 或者 Kruskal-Wallis。
这三类错误有个共同点:都能顺利跑出数字,都不报错,都看起来很像那么回事。这就是为什么我在团队里立了个规矩,凡是涉及分类变量的相关性分析,结果必须附上"方法选择说明"——用了什么方法、为什么用它、前提检查过没有。没有这段说明的结果一律打回。
2. 变量类型判定与数据前置体检
方法选对了一半,剩下一半靠数据本身干不干净。相关性分析对数据质量的要求比很多人想象的高:缺失值处理方式不同,相关系数能差出一大截;异常值不处理,Pearson 能被单个点带跑;类别频数太悬殊,卡方和方差分析都会失真。这一节讲的就是上桌前必须做的事。
2.1 连续、有序、名义:三种身份的边界
判定的核心问题是三个:取值能不能比大小、差值是否等距、零点是否有实在含义。三个都是"是",就是连续变量(收入、时长、温度、金额)。能比大小但差值不等距,是有序分类(满意度五级、学历层次、信用等级)。连大小都不能比,就是无序分类(城市、渠道、设备型号、是否点击)。
边界上有几个常被搞混的例子值得单独拎出来:
李克特量表。问卷里常见的"非常不满意到非常满意"五分制,绝大多数场合应该当有序分类处理。虽然实践中很多人直接当连续变量算均值,但那需要额外假设"相邻选项之间的心理距离相等",这个假设在题目设计良好时可以勉强成立,设计粗糙时就是瞎猜。稳妥起见,做相关分析时用 Spearman 而不是 Pearson,两种做法的结论一致时再放心用均值。
金额分段。业务报表里常见的"0-100 元、100-500 元、500-1000 元、1000 元以上"这种分箱,本质上是把连续变量切成有序分类。切完之后不该再当连续用,因为箱宽不一致,等距假设直接崩了。
布尔值。0/1 变量既可以看作二分类,也可以看作特殊的连续变量。这个双重身份很有用——正是因为二分类可以安全地编码成 0/1 而不损失信息,点二列相关才能在数学上等价于 Pearson 相关系数。但前提是它真的只有两个取值,而不是"金额为 0 或大于 0"这种被人工二分化的产物,后者属于人为切分,得用二列相关,后面会讲。
ID 类字段。订单号、用户 ID 这种,看着是数值,实际是名义分类,而且类别数等于样本量,任何基于频数的方法都会崩。这类字段在相关性分析里应该直接排除,不要试图分析。
判定完之后建议在数据字典里固化下来,写成枚举而不是靠记忆。项目一多人一多,同一个字段被两个人当两种类型用的情况太常见了,最后两份报告的结论对不上,排查成本极高。
2.2 上桌前先做的四件事
第一件:缺失值体检。相关性分析对缺失的处理有三种常见策略——按列删除、按对删除、插补。pandas.corr()默认是成对删除(pairwise),这意味着不同格子的相关系数是基于不同样本量算出来的,严格来说不可比。样本缺失有结构性的时候更麻烦:如果高收入人群普遍不填收入,那"收入与消费的相关"实际上只在中低收入群体里算出来的,解释时得打个大问号。我的习惯是先看缺失比例和缺失是否与其他变量相关,缺失超过 20% 的列直接标记为"谨慎使用"。
第二件:异常值体检。连续变量用箱线图或者 1.5 倍 IQR 规则扫一遍,重点看不合理的大值(比如年龄 200、时长负数)。异常值对 Pearson 的影响是压倒性的,对 Spearman 影响小很多,对基于秩的方法几乎无影响。所以如果你发现 Pearson 和 Spearman 结论不一致,八成是异常值在作祟,这时候要回到原始数据去看,而不是在两个系数之间挑一个好看的去汇报。
第三件:类别频数体检。分类变量看每一类的样本量。卡方检验的经验规则是期望频数不低于 5 的格子要占 80% 以上;方差分析虽然对组样本量的不平衡有一定容忍度,但一组只有 3 个样本、另一组有 300 个,均值的稳定性完全不在一个量级,F 检验的结果就不可靠了。遇到小类别,要么合并(注意合并要有业务依据,别把"高风险"和"中风险"合并了),要么改用精确检验或者非参数方法。
第四件:分布形态体检。连续变量看偏度和峰度,或者直接上直方图。严重右偏的收入、时长类变量,取对数之后往往更接近正态,这时候对数变换后的 Pearson 和原始的 Spearman 结论通常一致,可以互相印证。方差分析前要跑正态性检验(Shapiro 或 Anderson-Darling)和方差齐性检验(Levene),样本量大时正态性检验会过度敏感,这时候看 Q-Q 图比看 p 值靠谱。
这四件事做完,其实你已经对数据有了基本判断,该用哪个方法、结果能信几分,心里大概有数了。跳过这步直接算,后面返工的概率接近百分之百。
3. 连续对连续:拿来做对照的基线
分类对连续的方法,很多是从连续对连续的方法里"变形"出来的:点二列相关是 Pearson 的特例,方差分析本质上是组间方差与组内方差的比值,η² 是 r² 的多组推广。所以先把连续对连续这条基线捋清楚,后面的推导会顺畅很多。已经熟的人可以跳过这一节,但下面关于置信区间和样本量的部分,我建议还是扫一眼,因为那是分类对连续里同样存在的陷阱。
3.1 Pearson 与 Spearman 的取舍逻辑
Pearson 相关系数衡量的是线性关系,公式是协方差除以两个标准差的乘积,取值在 -1 到 1 之间。它对"线性"这两个字非常敏感:如果两个变量是完美的 U 型关系,Pearson 会算出接近 0 的结果,告诉你"没关系",这显然是错的。所以 Pearson 只能在肉眼确认散点图大致呈线性分布之后用。
Spearman 是 Pearson 在秩上的版本:把两个变量各自排序,用秩次代替原值算 Pearson。它的好处是只依赖顺序,不依赖具体数值,所以对单调非线性关系(比如对数关系)也能捕捉到,同时对异常值免疫。代价是它丢掉了"距离"信息,两个变量都是等距变量、关系确实线性时,Spearman 的功效略低于 Pearson。
实际项目里的取舍原则我一般这么定:
- 散点图线性、无明显异常值、双变量近似正态 → Pearson
- 单调但非线性、有异常值、变量是有序分类 → Spearman
- 样本量小(小于 30)→ Spearman 更稳,因为 Pearson 对正态假设的依赖在这个样本量下很明显
- 两边都要算,结论一致就放心报 Pearson,不一致就去看散点图找原因
"两个都算"这个做法看起来偷懒,其实是最实用的诊断手段。两者的差值本身就是信息:差了 0.1 以上,说明数据里要么有异常值、要么关系非线性,值得进一步查。
顺便说一句 Kendall's τ。它和 Spearman 都是基于秩的,但计算逻辑不同:Spearman 是"秩的 Pearson",Kendall 是"一致对减不一致对"的比例。Kendall 在小样本和大量并列值(ties)的情况下表现更稳,解释也更直观——τ = 0.3 大致意味着随机抽两对观测,同向的概率比反向高 30%。缺点是计算复杂度高,大样本时慢。样本量上万的数据,用 Spearman 就够了。
3.2 显著性不等于可靠:区间估计与样本量
这部分是我最想强调的。p < 0.05 只说明"在假设没有关系的前提下,观测到这个系数的概率很小",它不说明关系有多强。样本量一大,r = 0.05 也能显著;样本量一小,r = 0.6 也可能不显著。拿 p 值当"有没有关系"的判据,是相关性分析里最常见的误用。
正确的做法是同时看三样东西:系数大小、置信区间、样本量。
相关系数的置信区间可以用 Fisher z 变换来算。先把 r 转成 z = arctanh(r),z 的标准误是 1/√(n-3),据此算区间再反变换回 r。代码很简单:
import numpy as np from scipy import stats def corr_ci(r, n, alpha=0.05): z = np.arctanh(r) se = 1 / np.sqrt(n - 3) crit = stats.norm.ppf(1 - alpha / 2) lo, hi = np.tanh([z - crit * se, z + crit * se]) return lo, hi print(corr_ci(0.15, 2000)) # 窄区间,说明估计精细 print(corr_ci(0.15, 40)) # 宽区间,跨过 0,说明样本不足同样 r = 0.15,2000 个样本时区间大概是 [0.11, 0.19],40 个样本时区间可能从 -0.17 到 0.44,直接跨零。这两种情况下"r = 0.15"这句话的可信度天差地别,但很多报告只写一个数字。
效应量的经验参照(不同领域差异很大,仅作粗略参考):|r| 在 0.1 附近算小,0.3 中等,0.5 以上算大。社会科学里 0.2 已经值得讨论,物理实验里 0.99 才算合格。所以别照搬别人的阈值,要结合自己领域的常识。
还有一个常被忽略的点:相关系数会被测量误差稀释。如果两个变量都有随机测量误差,观测到的相关系数会系统性地低于真实相关(衰减偏误)。问卷数据尤其明显,一道题的测量信度只有 0.7,两个变量各有一半误差,观测相关就会被压低不少。所以拿到一个偏低的系数,别急着下"没关系"的结论,先看看变量本身测得准不准。
4. 分类对分类:卡方与它的效应量
分类对分类的情况在业务里非常常见:性别与是否流失、渠道与是否转化、会员等级与复购档位。这块的方法相对统一,但同样有一个高频错误——只报 p 值不报效应量。
4.1 卡方检验成立的前提
卡方检验的逻辑是:假设两个分类变量独立,按行合计和列合计推算每个格子的期望频数,然后比较观测频数与期望频数的偏离程度。偏离越大,越不像独立。统计量是 Σ(观测-期望)²/期望,自由度是 (行数-1)×(列数-1)。
用之前必须检查的三条:
期望频数。这是硬条件。任一格子的期望频数小于 1,或者超过 20% 的格子期望频数小于 5,卡方分布的近似就不可靠了,p 值会失真。处理方式:合并类别(合并要有业务逻辑,不能随便合),或者改用 Fisher 精确检验(2×2 表完全适用,更大的表可以用 Fisher-Freeman-Halton,计算量大但现代工具都能跑)。
样本独立性。每个观测只能落在一个格子里,不能重复计数。同一个用户既算"新客"又算"老客"的话,卡方就废了。这个在宽表拼接的时候特别容易出错,多表 join 之后行数膨胀,重复用户被数了两次。
类别定义一致。两个变量的分类标准要在时间上和口径上对齐。这个不是统计问题而是数据治理问题,但影响更大——去年按订单维度统计的"复购",今年按用户维度统计,两张表放一起做卡方,结果没有意义。
4.2 从 p 值走到 Cramér's V
卡方统计量本身不能当效应量用,因为它的数值随样本量线性增长:同样的分布偏离,样本量翻倍,卡方值翻倍。所以需要一个归一化的指标,这就是 Cramér's V:
V = sqrt(χ² / (n × (min(r, c) - 1)))取值 0 到 1,0 表示完全独立。2×2 表时 V 就等于 φ 系数,这时它的符号可以表示方向(正相关/负相关),但 Cramér's V 在更大的表里没有方向概念,只有强度。
这里有个细节要注意:Cramér's V 在小样本、大表的情况下会有向上的偏误,即真实无关的两个变量也可能算出 0.1 到 0.2 的 V。有个 Bergsma 校正公式可以修正,实践中如果表大于 3×3 且样本量不足 200,我建议用校正版:
def cramers_v_corrected(x, y): table = pd.crosstab(x, y) chi2 = stats.chi2_contingency(table, correction=False)[0] n = table.values.sum() r, k = table.shape phi2 = chi2 / n phi2corr = max(0, phi2 - (k - 1) * (r - 1) / (n - 1)) rcorr = r - (r - 1) ** 2 / (n - 1) kcorr = k - (k - 1) ** 2 / (n - 1) denom = min(kcorr - 1, rcorr - 1) return np.sqrt(phi2corr / denom) if denom > 0 else np.nan效应量大小的经验阈值,2×2 表时 φ 的 0.1/0.3/0.5 对应小/中/大;自由度更大时要按自由度调整,粗略的做法是把 V 除以 √(min(r,c)-1) 再比对。这块的争议一直存在,不同教材给的阈值不一样,我倾向的做法是:不要死守阈值,而是拿同类研究或历史数据做参照——"本次 V = 0.12,去年同口径是 0.11,基本稳定"这种表述比"V = 0.12 属于小效应"有用得多。
还有一点必须强调:卡方显著不等于有关系,更不等于有关系强度。n = 10000 的时候,V = 0.03 都能显著。这种"统计显著但业务无感"的结论,在流量数据里遍地都是。所以每次报卡方,我都会强制附上 V 和交叉表的百分比分布,让看图的人自己判断差异有没有业务意义。
5. 分类对连续:主线方法拆解
进入正题。分类变量与连续变量的相关性分析,核心思路是把"类别"翻译成"组",把"相关"翻译成"组间差异相对于组内波动的比例"。翻译得对不对,取决于分类变量有几个水平、这些水平有没有顺序。下面按二分类、多分类、有序分类三条支线拆,最后单独讲人为二分化的坑。
5.1 二分类对连续:点二列、t 检验与 Cohen's d
二分类变量(性别、是否会员、是否曝光)和连续变量(消费金额、停留时长)的关系,最标准的指标是点二列相关系数r_pb。它的定义和 Pearson 完全等价——只要把二分类编码成 0/1,两者数值一模一样。之所以单独立一个名字,是因为它的平方有明确的解释:r_pb² 等于"组间平方和占总平方和的比例",也就是这一组 0/1 标签能解释连续变量多少变异。
计算方式有两种,结果相同:
import numpy as np from scipy import stats # 方式一:直接调 scipy r_pb, p_pb = stats.pointbiserialr(df["is_member"], df["amount"]) # 方式二:手算,顺便验证公式 n1 = (df["is_member"] == 1).sum() n0 = (df["is_member"] == 0).sum() m1 = df.loc[df["is_member"] == 1, "amount"].mean() m0 = df.loc[df["is_member"] == 0, "amount"].mean() s = df["amount"].std(ddof=1) n = n1 + n0 r_manual = (m1 - m0) / s * np.sqrt(n0 * n1 / (n * (n - 1)))公式里的 (m1 - m0)/s 就是标准化均值差,后面的根号项是一个和组样本量有关的系数。这告诉我们一个重要事实:r_pb 的大小同时受"均值差有多大"和"两组样本量是否均衡"影响。同样差 0.5 个标准差,两组各 50 人的 r_pb 会明显大于 90 人 vs 10 人的情况。所以两组人数悬殊时,r_pb 会偏小,别误读成"关系弱"。
与之对应的、不受样本量比例影响的指标是Cohen's d:
d = (m1 - m0) / s_pooledd 只反映"差了几个标准差",0.2/0.5/0.8 对应小/中/大。在 AB 实验报告里,d 比 r_pb 更常用,因为它稳定、可比、和样本量解耦。我的习惯是两个都报:d 给业务看差距大小,r_pb² 给建模看解释力。
t 检验和点二列相关的关系也得说清楚。独立样本 t 检验检验的是"两组均值是否相等",它和 r_pb 的显著性检验在数学上等价:t² = (n-2) r_pb² / (1 - r_pb²)。换句话说,同一个问题,t 检验和点二列相关给出完全一致的 p 值,只是一个从"均值差"角度讲,一个从"共变"角度讲。别把两者当成两种独立证据去交叉验证,它们是同一个东西的两副面孔。
实操上要注意 t 检验的方差齐性。默认的 Student t 检验假设两组方差相等,不满足时会偏。直接用 Welch t 检验(equal_var=False),它在方差不齐时表现更好,方差齐时损失的功效很小,所以我现在基本默认用 Welch。样本量很小或者分布严重偏斜时,改用 Mann-Whitney U 检验,配套的效应量用秩二列相关:
u, p = stats.mannwhitneyu(g1, g0, alternative="two-sided") r_rb = 1 - 2 * u / (len(g1) * len(g0)) # 取值 -1 到 15.2 多分类对连续:方差分析、eta 方与 Welch 兜底
三个及以上类别(会员等级、城市线级、渠道来源)对连续变量,标准方案是单因素方差分析。它的核心思想是分解方差:总变异 = 组间变异 + 组内变异。如果组间变异相对于组内变异足够大,就认为各组均值不全相等。
F 统计量 = (组间均方) / (组内均方),自由度是 (k-1, n-k)。这里 k 是组数,n 是总样本量。F 显著只说明"至少有一组和其他组不同",不说明是哪一组,也不说明差异有多大——所以必须做事后多重比较,而且必须校正。
效应量用η²(eta squared):
η² = SS_between / SS_total解释是"分类变量解释了连续变量百分之多少的变异",取值 0 到 1。它和 r² 是同一个东西在多组情形下的推广,二分类时 η² 就等于 r_pb²。这个是很多人不知道的细节,但很有用:它意味着你可以把不同方法算出来的效应量放在同一把尺子上比较。
η² 有个已知问题:它是有偏的,样本量小时会高估。所以还有个ω²(omega squared)和偏 η²(partial eta squared)。偏 η² 是统计软件默认输出的那个(比如某些工具的 ANOVA 表),但它在单因素设计里和 η² 数值一样,多因素设计里才不同。报告时写清楚用的是哪个,别混。
手算 η² 很简单,比调包还快:
grand_mean = df["amount"].mean() ss_between = sum(len(g) * (g["amount"].mean() - grand_mean) ** 2 for _, g in df.groupby("level")) ss_total = ((df["amount"] - grand_mean) ** 2).sum() eta2 = ss_between / ss_total f_stat, p_val = stats.f_oneway(*[g["amount"] for _, g in df.groupby("level")])前提检查两条:正态性(每组残差近似正态)和方差齐性(Levene 检验 p > 0.05)。方差不齐时用Welch ANOVA,它不假设方差相等,自由度用 Welch-Satterthwaite 近似:
import pingouin as pg pg.welch_anova(data=df, dv="amount", between="level")严重偏态、异常值多、或者有序分类当无序处理时,换Kruskal-Wallis检验。它是方差分析的非参数版本,比较的是各组的秩均值。显著之后的事后比较用 Dunn 检验(配 Bonferroni 或 BH 校正),效应量可以用ε²(epsilon squared):
ε² = (H - k + 1) / (n - k)其中 H 是 Kruskal-Wallis 统计量,k 是组数,n 是总样本量。ε² 的解释和 η² 类似,但通常略小。
事后比较这一步千万别省,也千万别用 t 检验连打好几次还不校正。k 组两两比较有 k(k-1)/2 次,5 组就是 10 次,不做校正的话假阳性率会从 5% 飙到 40% 左右。校正方法用 Bonferroni(保守,比较次数少时够用)或者 Benjamini-Hochberg(控制 FDR,比较次数多时更合适)。
5.3 有序分类对连续:别把等级当哑变量
有序分类是三种类型里最尴尬的:它比无序分类多了一层顺序信息,比连续变量少了一层距离信息。处理方式取决于你想问什么问题。
如果只想问"有没有差异",把它当无序分类,走方差分析或 Kruskal-Wallis。这样做不损失正确性,但浪费了顺序信息,功效会低一些。
如果想问"有没有单调趋势",那就要用能利用顺序的方法。最直接的是把等级编码成 1、2、3、4、5,然后和连续变量算Spearman 秩相关。这个做法在等级数不少于 5 且分布大致均匀时表现不错;等级只有 3 档时会比较粗糙,但方向性判断仍然可用。
更严谨的做法是Jonckheere-Terpstra 趋势检验,它专门检验"随着等级升高,连续变量的分布是否有序移动",不要求等距假设,功效比 Kruskal-Wallis 高。这个检验在常规工具里不一定有,R 里有clinfun::jonckheere.test,Python 里可以用scikit-posthocs或者自己按公式实现。
绝对要避免的做法是把有序分类当无序分类拆成哑变量,然后逐个看系数。这样得到的是一堆"等级 2 比等级 1 高多少""等级 3 比等级 1 高多少"的比较,完全丢掉了"等级越高越好"这个整体趋势,解释起来支离破碎。如果确实要进回归模型,用正交多项式对比或者直接当连续变量处理(前提是等距假设能接受),都比拆哑变量强。
有一次做用户等级和客单价的分析,一开始按五个等级拆哑变量跑 ANOVA,F 检验显著但事后比较的结果很乱,看不出规律。后来换成 Spearman,ρ = 0.31,p 值极小,一句话就说清了"等级越高客单价越高"。同一个数据,换个方法,可解释性天差地别。
5.4 手工二分化连续变量的代价:二列相关
这个坑值得单独开一节,因为它太容易被当成"正常操作"。
场景是这样的:你有一个连续变量(比如用户活跃度得分),业务上想按某个阈值切成"高活跃/低活跃",然后分析它和另一个连续变量(消费金额)的关系。切完之后你算了个点二列相关,得到 r = 0.28。这个数字是低估的。
原因是二分化丢掉了组内的全部变异信息。原本活跃度从 10 分到 90 分是一个连续谱,切完之后只剩"高/低"两个标签,中间那么多信息全没了。统计上这叫人为二分化的信息损失,会导致相关系数系统性偏小(除非切点在均值处且分布正态)。
如果你手上就是二分后的数据,对应的正确指标是二列相关系数(biserial correlation),它估计的是"假如那个连续变量没被切,原本的相关会是多少"。公式是:
r_bis = r_pb × sqrt(p × q) / φ(t)其中 p 和 q 是两组的比例,φ(t) 是标准正态分布在切分点处的密度值(切分点 t 是使正态分布上侧面积为 p 的那个 z 值)。算出来 r_bis 通常会比 r_pb 大一些,因为它在补偿信息损失。
具体怎么算:
from scipy.stats import norm def biserial(r_pb, p): q = 1 - p t = norm.ppf(q) # 切分点 density = norm.pdf(t) # 该点处的正态密度 return r_pb * np.sqrt(p * q) / density p_high = (df["active_flag"] == 1).mean() print(biserial(0.28, p_high))注意二列相关的假设很强:它要求被二分的那个变量本身服从正态分布,且切分是人为的(不是天然的二分变量)。如果活跃度本身就是双峰的,或者"高/低"是天然存在的两类用户,那就该用点二列,不该用二列,硬套会得到错误结果。
我的实际建议是:能不切就不切。业务上需要"高/低"这种标签,可以在报告层做分组展示,但分析层始终保留连续形式。真要切,切点要有依据(比如业务阈值、分位数),而且要意识到所有基于二分数据的相关系数都是下界估计。
6. 完整实操:一套能跑的相关性矩阵
前面讲的是分散的方法,实际项目里真正有用的是"输入一张表,自动按变量类型分发到正确方法,输出统一的效应量和 p 值"。这一节给一套能直接抄的代码。
6.1 环境准备与模拟数据
依赖三个库:pandas 做数据处理,scipy 做核心检验,pingouin 补充 Welch ANOVA 和偏 η² 这类 scipy 里没有的实现。
pip install pandas numpy scipy pingouin scikit-posthocs造一份覆盖三种变量类型的模拟数据。为了让后面的结果可验证,我按已知的组间差异生成数据:
import numpy as np import pandas as pd from scipy import stats import pingouin as pg rng = np.random.default_rng(20240607) n = 480 # 无序分类:三个渠道 channel = rng.choice(["自然", "投放", "推荐"], size=n, p=[0.5, 0.3, 0.2]) # 有序分类:五个会员等级 level = rng.choice([1, 2, 3, 4, 5], size=n, p=[0.3, 0.25, 0.2, 0.15, 0.1]) # 二分类:是否会员 is_member = (level >= 3).astype(int) # 连续变量:消费金额,受渠道和等级共同影响,加噪声 base = np.where(channel == "自然", 120, np.where(channel == "投放", 145, 165)) amount = base + level * 12 + rng.normal(0, 45, n) amount = np.clip(amount, 0, None) df = pd.DataFrame({ "channel": channel, "level": level, "is_member": is_member, "amount": amount, "visits": rng.normal(8, 3, n) + level * 0.8 })这份数据里我埋了两条已知的关系:渠道影响金额、等级影响金额,后面就看方法能不能把它们捞出来。
6.2 逐类方法的代码与结果解读
先跑二分类对连续:
# 点二列相关 r_pb, p_pb = stats.pointbiserialr(df["is_member"], df["amount"]) print(f"点二列 r = {r_pb:.3f}, p = {p_pb:.4g}") # Cohen's d(Welch) g1 = df.loc[df["is_member"] == 1, "amount"] g0 = df.loc[df["is_member"] == 0, "amount"] sp = np.sqrt(((len(g1)-1)*g1.var(ddof=1) + (len(g0)-1)*g0.var(ddof=1)) / (len(g1) + len(g0) - 2)) d = (g1.mean() - g0.mean()) / sp print(f"Cohen's d = {d:.3f}") # Welch t 检验 t, p_t = stats.ttest_ind(g1, g0, equal_var=False) print(f"Welch t = {t:.3f}, p = {p_t:.4g}")跑出来应该是 d 在 0.9 到 1.2 之间,r_pb 在 0.4 上下。这里能验证前面说的等价性:t² 和 r_pb 换算出来的值应该严丝合缝对得上(verification:t² = (n-2)r²/(1-r²))。
再跑多分类对连续:
# 方差分析 groups = [g["amount"].values for _, g in df.groupby("channel")] f_stat, p_anova = stats.f_oneway(*groups) # eta 方 grand = df["amount"].mean() ss_between = sum(len(g) * (g["amount"].mean() - grand) ** 2 for _, g in df.groupby("channel")) ss_total = ((df["amount"] - grand) ** 2).sum() eta2 = ss_between / ss_total print(f"F = {f_stat:.3f}, p = {p_anova:.4g}, eta2 = {eta2:.3f}") # 方差齐性 lev = stats.levene(*groups) print(f"Levene p = {lev.pvalue:.4g}") # 事后比较(Dunn + BH 校正,非参数路线) import scikit_posthocs as sp_ print(sp_.posthoc_dunn(df, val_col="amount", group_col="channel", p_adjust="fdr_bh"))注意看 Levene 的结果。方差齐的话 ANOVA 没问题;不齐就换pg.welch_anova。这份模拟数据里三组的标准差都设成 45,理论上应该齐,但抽样波动可能让 p 值在 0.05 附近晃,这时候别纠结,直接上 Welch 更省心。
有序分类对连续:
rho, p_rho = stats.spearmanr(df["level"], df["amount"]) print(f"Spearman rho = {rho:.3f}, p = {p_rho:.4g}") # 对照:当连续变量处理算 Pearson r_pear, p_pear = stats.pearsonr(df["level"], df["amount"]) print(f"Pearson r = {r_pear:.3f}")两个数应该很接近(因为我在生成数据时让等级对金额的影响是线性的),这验证了"等距假设成立时两种方法一致"。如果真实数据里等级对金额是加速上升的(比如等级 5 突然贵很多),Spearman 会明显高于 Pearson,那时候更要信 Spearman。
6.3 组装成自动分发的相关矩阵
把上面的逻辑封装成函数,按变量类型自动选方法。核心是一个类型判定函数和一个分发函数:
def detect_type(s, max_cat=10): if s.dtype == object or s.dtype.name == "category": n_unique = s.nunique() return "ordinal" if is_ordered(s) else "nominal" n_unique = s.nunique() if n_unique <= 2: return "binary" if n_unique <= max_cat and (s == s.astype(int)).all(): return "ordinal" # 少数几个整数档,默认当有序 return "continuous" def pair_assoc(df, col_a, col_b): ta, tb = detect_type(df[col_a]), detect_type(df[col_b]) a, b = df[col_a], df[col_b] # 统一成"分类在前、连续在后"的顺序处理 if ta == "continuous" and tb != "continuous": a, b, ta, tb = b, a, tb, ta if ta == "continuous" and tb == "continuous": r, p = stats.pearsonr(a, b) return {"method": "pearson", "effect": r, "p": p} if ta in ("nominal", "ordinal") and tb == "nominal" and ta == "nominal": table = pd.crosstab(a, b) chi2, p, _, _ = stats.chi2_contingency(table) v = np.sqrt(chi2 / (len(a) * (min(table.shape) - 1))) return {"method": "cramers_v", "effect": v, "p": p} if tb == "continuous": if ta == "binary": r, p = stats.pointbiserialr(a, b) return {"method": "point_biserial", "effect": r, "p": p} if ta == "ordinal": r, p = stats.spearmanr(a, b) return {"method": "spearman", "effect": r, "p": p} f, p = stats.f_oneway(*[g[b.name].values for _, g in df.groupby(a.name)]) grand = b.mean() eta2 = sum(len(g) * (g[b.name].mean() - grand) ** 2 for _, g in df.groupby(a.name)) / ((b - grand) ** 2).sum() return {"method": "anova_eta2", "effect": eta2, "p": p} return {"method": "unsupported", "effect": np.nan, "p": np.nan}这个函数不完美,is_ordered需要你根据业务字典来实现,类型判定也有模糊地带,但它把方法论固化进了代码,比每次人工判断靠谱得多。真要用的话,补三件事:多重比较校正(循环里累积 p 值之后统一做 BH 校正)、样本量过滤(有效样本少于 30 的对直接跳过)、以及把结果写成带方法名和置信区间的表,方便回溯。
还有一点经验之谈:自动化分发的结果一定要抽查。我做过一次,自动跑完 200 个变量对,发现有几十对的 p 值全是 1.0,查下来是一列常量(某个渠道在某个时间段没有数据),方差为 0,检验直接失效。没有异常检测的自动化流程,出错的时候是静默的。
7. 常见问题速查与排查思路
方法学完了,剩下的就是踩坑。这一节把我在实际项目里反复遇到的问题整理成表,附上判断依据和处置方式。
7.1 问题现象与对应处置表
| 现象 | 可能原因 | 判断依据 | 处置方式 |
|---|---|---|---|
| Pearson 显著但 Spearman 不显著 | 异常值或非线性 | 散点图、箱线图 | 查异常值,或改报 Spearman |
| 相关系数极小但 p 极小 | 样本量过大 | 看 n 和置信区间 | 报效应量与区间,别只报 p |
| 卡方显著但看不出差异 | 期望频数或样本量问题 | 检查交叉表百分比 | 报 Cramér's V 和百分比分布 |
| ANOVA 显著但事后比较全不显著 | 未校正 / 组数太多 | 看比较次数 | 用 BH 校正,或减少比较 |
| 两组均值差很大但 t 检验不显著 | 组内方差极大 / 样本小 | 看标准差和 n | 增大样本或改用非参数 |
| 点二列相关系数偏低 | 两组样本量悬殊 | 看两组比例 | 改用 Cohen's d 判断 |
| 二分类编码后相关方向反了 | 编码顺序与业务方向不一致 | 看 0/1 对应什么 | 明确编码约定并写入文档 |
| 多个变量两两比较出现大量显著 | 多重比较膨胀 | 比较次数 | FDR 校正或降维后再看 |
表里第三条我要多嘴一句。卡方显著但事后看不出差异,八成是交叉表的百分比差异很小、但样本量大到能检出。这种时候别硬讲"有关系",直接写"存在统计关联但效应量很小(V = 0.04),业务上可视为无实质差异",比编一个故事强。
第七条也很常见。is_member如果 1 表示非会员,算出来的相关系数符号就是反的,解释的时候会得出完全相反的结论。这个错误极其低级但极其常见,尤其是多人协作、字段含义靠口口相传的项目。我的做法是在数据字典里强制写"取值含义"这一列,而且要求写具体("1 = 是会员",不是"1 = 是")。
7.2 我踩过的几个坑
第一个坑是把"相关"当"因果"用。这个错误不是统计层面的,是表达层面的。有一次做渠道与客单价的分析,算出来推荐渠道的客单价显著更高,报告里写了一句"推荐渠道带来更高客单价",被业务方理解成"多投推荐渠道就能提升客单价"。实际上很可能是高价值用户更倾向于走推荐路径,方向完全相反。后来我改了个规矩:所有相关性结论只写"关联""伴随""差异",动词一律不用"导致""带来""提升"。这个习惯救过我好几次。
第二个坑是忽略辛普森悖论。分组看和整体看结论相反的情况,在分类对连续的分析里特别容易发生。比如整体看会员的客单价更高,但按城市分层之后,每个城市里会员的客单价都更低——原因可能是会员在高线城市占比更高,而高线城市本身客单价就高。这个坑的排查方式很朴素:算完整体,一定要按关键维度分层再算一遍,两层结论不一致就说明有混杂变量,得用分层分析或者回归控制。
第三个坑是样本量不均衡导致的误判。多分类对连续,一组 400 人一组 20 人,ANOVA 的 F 值主要被大组驱动,小组的均值极不稳定却参与了组间平方和的计算。这种情况我一般会做两件事:一是看每组的置信区间,小组的区间宽得离谱就直接标注"样本不足,结论仅供参考";二是补一个 bootstrap 的组间差异估计,看看重抽样下结论稳不稳。
第四个坑是测量工具的刻度问题。问卷里的李克特量表,从"非常不满意"到"非常满意"标 1 到 5,看起来是等距的,但受访者对"满意"和"非常满意"的心理距离,跟"一般"和"满意"的距离往往不一样。这种不等距在算 Pearson 时会引入偏差,算 Spearman 就没问题。所以但凡数据来源是主观量表,我现在默认先跑 Spearman,Pearson 只作为对照。
第五个坑是把统计不显著当成"没有关系"。这个在样本量小的时候特别危险。n = 40 的时候,即使真实相关有 0.3,检出的概率也不到一半。所以不显著时正确的表述是"在当前样本量下未检出显著关联",而不是"两者无关"。如果业务上这个问题很重要,那就该去算需要多大样本量才能检出预期效应,而不是拿一个不显著的结果草草收场。
这套东西用熟了之后,回头看那些直接df.corr()一把梭的做法,会觉得挺粗糙的。但说实话,大多数业务场景也不需要多精细——知道该看哪个指标、知道 p 值不能单独用、知道分组之后要再检查一遍,就已经能避开 80% 的坑了。剩下的 20%,等你真的被一个辛普森悖论坑过一次,自然就记住了。