MATLAB相关性分析：皮尔逊、秩相关与显著性避坑

简介这份资料包面向使用MATLAB开展相关性分析数学建模的学习者与参赛者围绕变量关联程度的度量与实现展开涵盖皮尔逊相关系数、斯皮尔曼秩相关系数等常见方法适合需要完成统计分析、建模复现或课程作业的初中级读者。包内共11个文件以3个xls与1个xlsx体测数据、3个m脚本、2个mat数据文件为主并配有1个docx说明和1个jpg预览图压缩包约194KB便于直接导入MATLAB运行与核对结果。已有2409人学习下载。读者可从中获得相关系数计算、数据预处理、热力图绘制及结果解释的完整示例也能借助八年级体测数据案例理解线性关系与单调关系的差异并参考临界值文档判断显著性。脚本与数据文件便于复现实验说明文档可辅助梳理方法要点适合作为相关性分析建模的入门与查漏补缺材料。1. 相关系数不是「相关性」两列数就能骗过你的眼睛MATLAB 相关性分析几乎是数学建模里最先被调用的一类工具很多人打开编辑器敲下的第一行就是corrcoef(X)然后盯着那个落在 -1 到 1 之间的数字决定接下来往哪个方向建模。这个动作方便得像呼吸也危险得像呼吸——四组数据的散点图形状天差地别算出来却可以是同一个相关系数反过来两列数由同一个确定性公式生成依赖关系牢不可破相关系数却可能正好是 0。相关系数衡量的从来不是「两个变量有没有关系」而是「这两个变量在多大程度上能被一条直线解释」前者是客观世界的属性后者是我们自己选定的模型假设。搞清楚这条分界线是后面所有计算的立足点。2. 皮尔逊系数的分子与分母协方差被标准化之后还剩下什么2.1 从协方差到相关系数量纲在除法里消失协方差是相关系数的原料。它的定义是两个变量中心化之后的乘积和再取平均表达的是「x 偏离自己均值的时候y 是不是也朝着同一个方向偏离」。问题在于协方差的数值带着单位身高用米算一次、换成厘米再算一次得到的协方差相差 100 倍。这种量纲依赖让它没法跨数据集比较也没法判断「大」到底算多大。皮尔逊系数做的事情非常朴素——把协方差除上两个标准差。分子分母里的单位相互约掉结果变成一个无量纲的纯数取值范围被压缩到 [-1, 1]。这个变换看起来只是形式上的整理但它把「协同变化有多强」从绝对量变成了相对量从此不同数据集之间才有了可比性。% 同一组数据只改量纲看协方差和相关系数各自的表现 x [1 2 3 4 5 6]; % 单位米 y [2.1 4.3 5.9 8.2 9.8 12.1]; cov_raw cov(x, y) % 非对角元约为 6.88 r_raw corrcoef(x, y); % 非对角元约为 0.9989 x_cm x * 100; % 换成厘米 cov_cm cov(x_cm, y) % 非对角元被放大 100 倍 r_cm corrcoef(x_cm, y); % 非对角元纹丝不动cov返回的是 p×p 矩阵非对角元就是协方差单位是「米·目标单位」corrcoef返回同形状矩阵非对角元是皮尔逊系数对角线上恒为 1。上面对比的关键在于cov_cm的数值随量纲变化r_cm却与r_raw完全一致。这一点在建模里非常重要——做量纲统一是为了优化求解稳定但对相关性分析本身量纲不改结论。反过来说如果有人拿协方差大小来比较哪两个变量「更相关」那基本等于没做标准化就直接比长度血泪经验告诉我们这种结论很容易翻车。2.2 柯西-施瓦茨不等式|r| ≤ 1 不是巧合把中心化后的两个列向量记作 x̃ 和 ỹ那么皮尔逊系数可以写成内积除以两个范数之积r x̃, ỹ / (‖x̃‖ · ‖ỹ‖)这正是两个向量夹角余弦的定义。于是「|r| 不超过 1」就不再是一条需要记住的经验规律而是柯西-施瓦茨不等式的直接推论任意两个向量的内积绝对值永远不超过它们长度之积。等号成立的条件也很明确——两个向量共线也就是 ỹ k·x̃说明 y 与 x 之间是严格的线性关系。这个几何视角顺手解决了好几个常见困惑。第一r 的符号就是夹角的方向正数表示同向偏离负数表示反向偏离绝对值表示夹角接近 0 还是接近 90 度。第二r 为 0 只意味着两个向量正交也就是「没有线性成分」完全不排除二次、周期、分段等其他形式的依赖抛物线两端的数据点可以完美对称地分布在一条水平线两侧r 正好算出来是 0。第三r 对线性变换敏感——你在建模时对某个变量做了对数变换等于把它的向量方向掰弯了r 自然会变这不是算错了而是问题被重新定义了。2.3 矩阵视角corrcoef 一次算完全部变量对数据矩阵 X 是 n×p每一列是一个变量每一行是一次观测。corrcoef(X)返回 p×p 的对称矩阵第 (i,j) 个元素就是第 i 列和第 j 列的皮尔逊系数。这个函数的内部动作可以用三步概括按列中心化、求各列标准差、用中心化矩阵做一次矩阵乘法并逐元素除以标准差的乘积。function R naive_corr(X) % X: n×p每列一个变量返回 p×p 相关矩阵 Xc X - mean(X, 1); % 按列中心化n×p sd std(X, 0, 1); % 样本标准差分母 n-1 R (Xc * Xc) ./ (sd * sd) / (size(X,1) - 1); endstd(X, 0, 1)的三个参数分别是权重标志取 0 表示按 n-1 归一化与cov保持一致维度取 1 表示按列计算。Xc * Xc得到的是未除自由度的离差阵再逐元素除以sd*sd就完成了标准化。拿一个随机矩阵验证max(abs(naive_corr(X) - corrcoef(X)), [], all)误差通常在 1e-15 量级说明这几十行就是官方函数的核心。理解这个结构能带来两个实际好处。其一是复杂度矩阵乘法的代价是 O(n·p²)所以变量个数 p 上千时corrcoef会明显变慢而真正的瓶颈是 p² 而不是样本量 n这时候要么先降维要么只计算关心的变量对。其二是corrcoef和corr的分工——前者只做皮尔逊、只返回系数矩阵、遇到 NaN 会污染整行整列后者支持Type、Rows、Tail参数还能同时吐出 p 值和置信区间是建模时的主力。很多人混着用结果在缺失值场景下拿到一片 NaN 却找不到原因这类坑后面单独说。3. p 值那一层一个 r 和样本量怎么滚出显著性3.1 t 统计量与自由度 n-2 的来源r 本身只是一个描述统计量它告诉你样本里线性关系有多强但没法告诉你「这个强度是不是偶然」。要回答后者得给 r 配一个零假设总体相关系数为 0。在这个假设下统计量t r · sqrt( (n-2) / (1 - r²) )服从自由度为 n-2 的 t 分布。自由度为什么是 n-2不是 n-1因为计算 r 的过程中估计掉了两个均值x 的均值和 y 的均值每个估计消耗一个自由度剩下的信息量才是 n-2。这个细节在做小样本建模时影响很大n 本来就只有十几个的时候少一个自由度足以把显著性判反。rng(7); n 12; x randn(n,1); y 0.6*x 0.8*randn(n,1); [r, p] corr(x, y, Type, Pearson); t r * sqrt((n-2) / (1 - r^2)); % 检验统计量 p_manual 2 * (1 - tcdf(abs(t), n-2)); % 双侧 p 值 fprintf(corr 给 p%.6f手算 p%.6ft%.4f\n, p, p_manual, t);tcdf是 t 分布的累积分布函数abs(t)取绝对值是因为双侧检验只看偏离零的幅度最后乘 2 把两侧尾巴加起来。手算结果与corr输出应当一致到小数点后多位如果对不上八成是漏了abs或者自由度写成了 n-1。搞清这一层之后p 值就不再是黑匣子里的魔法数字——它是 r 和 n 联合组装出来的一个概率样本量一变同样的 r 可以得出完全相反的结论。3.2 为什么 n 只有 8 的时候 r0.7 依然不显著在小样本数学建模任务里这种情况极其常见跑出来 r0.70看着不算低兴冲冲写进论文说「显著相关」评审一句话就能戳穿。原因在于临界值随自由度急剧变化。把上面的 t 公式反解就能得到 r 的临界值r_crit t_crit / sqrt( n - 2 t_crit² )n 8; alpha 0.05; t_crit tinv(1 - alpha/2, n-2); % 双侧临界 t 值 r_crit t_crit / sqrt(n - 2 t_crit^2) % 对应的 r 临界值tinv返回指定累积概率下的 t 分位点1-alpha/2对应双侧 0.05 的上侧分位。n8 时自由度是 6临界 t 约为 2.447对应的 r 临界值约 0.707。也就是说 r0.70 恰好落在临界线内侧p 值在 0.05 上方一点点结论就是「没有足够证据认为相关」。同样的 r0.70如果 n30自由度 28临界 r 降到 0.36 左右显著性轻轻松松。这个数字对比说明了一件反直觉的事相关系数的大小不能脱离样本量单独解读。判断「相关强不强」和判断「相关是否可信」是两套逻辑前者看 r 的绝对值后者看 r 与 n 的组合。建模报告里两个都要写只给 r 不给 n 的说理是不完整的。3.3 Tail 参数与多重比较p 值的通胀corr的Tail参数控制检验方向。默认both做双侧适合「有没有关系」这种问题right只检验正相关left只检验负相关单侧检验的临界值更宽松但前提是方向必须在看数据之前就已经确定。见过太多人先算出 r 是正的再回头把Tail改成right去凑显著性这已经属于把结论倒推回假设统计上站不住。更隐蔽的问题是多重比较。一个 p 列的变量集corr返回的 p 值矩阵看起来是 p×p实际上由于对称性独立的检验只有 p(p-1)/2 个对角线上那 p 个 1 毫无信息量。做 10 个变量就有 45 次检验即使所有变量之间毫无关系按 0.05 的显著性水平平均也会冒出 2 个多「显著」的结果。R corr(X); m size(R,1); n_test m*(m-1)/2; % 实际独立检验次数 p_adj min(1, P .* n_test); % Bonferroni 校正P 为 corr 返回的 p 值矩阵上面是最粗糙也最保守的 Bonferroni 思路把每个 p 值乘以检验次数。更温和的做法是控制错误发现率MATLAB 里有mafdr可用。校正与否取决于你在做什么探索性筛变量可以放宽写进结论的显著性判断必须校正否则那 45 次检验里挑出来的两个「显著」早晚会变成模型里的噪声。3.4 Fisher z 变换置信区间从哪来corr还能返回置信区间第三、第四个输出参数就是上下界。它的算法不是直接对 r 加减标准误因为 r 的抽样分布是偏的接近 ±1 时会明显不对称。标准做法是先做 Fisher z 变换把 r 压到近似正态的尺度上算完区间再反变换回去。[r, p, rlo, rhi] corr(x, y); z atanh(r); % Fisher 变换 se 1 / sqrt(n - 3); % z 尺度上的标准误 ci tanh(z [-1 1] * norminv(0.975) * se); fprintf(corr 区间 [%.4f, %.4f]手算 [%.4f, %.4f]\n, rlo, rhi, ci(1), ci(2));atanh和tanh是变换与反变换norminv(0.975)给出标准正态的 97.5% 分位点即 1.96n-3 里的 3 同样来自自由度修正。手算区间与corr输出在小样本下会有细微差异因为官方实现可能用了更精确的近似但数量级和形态是一致的。置信区间的价值在于它直接暴露了估计的不确定性r0.75 配 n12区间可能宽到 [0.30, 0.93]这种宽度足以说明「点估计不能当结论用」。4. 秩与非线性斯皮尔曼、肯德尔换掉的到底是什么4.1 秩变换把「值」降级成「序」之后剩下什么皮尔逊系数吃的是数值本身所以对异常值和单调非线性都很敏感。斯皮尔曼的做法很干脆先把两列数据各自换成秩在秩上再跑一遍皮尔逊。秩变换只保留大小顺序扔掉具体数值和间隔信息于是任何严格单调的变换都不再影响结果。x (1:20); y exp(0.3*x) 0.5*randn(20,1); % 指数关系明显非线性 r_pear corr(x, y, Type, Pearson) % 被曲率拉低 r_sper corr(x, y, Type, Spearman) r_sper_manual corr(tiedrank(x), tiedrank(y), Type, Pearson)tiedrank返回平均秩没有并列时就是 1 到 n 的排序有并列时把并列位置的秩取平均保证秩和不变。r_sper与r_sper_manual应当完全相等这说明Spearman参数并不是另起炉灶的另一套公式而是在内部先做了秩变换。指数关系的例子里皮尔逊系数被曲率压到 0.85 以下斯皮尔曼系数却接近 1因为「x 越大 y 越大」这个序关系始终没变。这也解释了斯皮尔曼为什么在评分类、等级类数据上特别顺手。代价同样明确秩变换之后「差一点点」和「差很多」被一视同仁。如果业务上关心的是幅度而非顺序比如剂量与反应强度的关系斯皮尔曼会丢掉关键信息。4.2 肯德尔的一致对与不一致对肯德尔走的是另一条路它不比较秩的大小而是清点「对」。对任意两个观测 i 和 j如果 x 的排序方向和 y 的排序方向相同这一对叫一致对方向相反叫不一致对。统计量就是一致对减去不一致对再除以总对数。n 8; x [1 2 3 4 5 6 7 8]; y [1 2 3 4 5 7 6 9]; % 有一处局部交换 C 0; D 0; for i 1:n-1 for j i1:n s sign(x(i)-x(j)) * sign(y(i)-y(j)); C C (s 0); % 一致对 D D (s 0); % 不一致对 end end tau_a (C - D) / (n*(n-1)/2) % 无并列时的 tau-a tau_b corr(x, y, Type, Kendall) % MATLAB 默认算 tau-b内层循环遍历所有 ij 的组合sign的乘积为正说明两列同向为负说明反向。总对数是 n(n-1)/2本例中 28 对里有 1 对不一致tau 约 0.9286。MATLAB 的Kendall默认算的是 tau-b它在分母上对并列做了修正公式为(C-D) / sqrt((n0-n1)(n0-n2))其中 n0 是总对数n1、n2 分别是 x 列和 y 列内部并列对的数量。数据里没有并列时 tau-a 与 tau-b 相等一旦出现大量并列比如五级量表就必须用 tau-b否则分母会虚高。肯德尔与斯皮尔曼的差异在实际数据上能看出来。斯皮尔曼对「秩的差值的平方」敏感一个离群点被推到秩的极端影响会被平方放大肯德尔只看两个观测的相对顺序单个点的移动最多影响它与其余 n-1 个点的配对关系稳健性通常更好代价是朴素实现是 O(n²)样本上万时明显吃力。4.3 三种系数的选择依据系数捕捉的关系类型对异常值计算代价MATLAB 参数皮尔逊严格线性敏感O(n)Pearson斯皮尔曼任意单调中等O(n log n)Spearman肯德尔任意单调基于序对较稳健O(n²)Kendall选择路径我一般这么走先画散点图形状接近直线且没有明显的极端点用皮尔逊散点呈现单调曲线指数、对数、饱和型或者数据本身是等级、排名用斯皮尔曼样本量不大但并列值多、或者怀疑个别点会影响结论用肯德尔并同时看 tau-b。三种都算一遍然后挑最大的报这是典型的踩坑做法——挑最大等于把选择过程本身变成了一个未记录的检验结论的可信度就打了折扣。值得注意的是三种系数都只能捕捉单调关系。U 形、环形、周期性这些形态三者都可能给出接近 0 的结果此时需要的是距离相关或者分箱后重新审视这部分留到后面展开。5. 相关性建模的避坑清单五种翻车现场与排查路径5.1 现象分组内全是负相关合并后却成了强正相关这组数据的每一组内部都在下降整体一算却是 0.9 以上的正相关。G1 [1 2 3; 4 3 2]; % 组内 r -1 G2 [5 6 7; 10 9 8]; % 组内 r -1 G3 [9 10 11; 16 15 14]; % 组内 r -1 fprintf(组内 r: %.2f %.2f %.2f\n, corr(G1(:,1),G1(:,2)), ... corr(G2(:,1),G2(:,2)), corr(G3(:,1),G3(:,2))); X [G1(:,1); G2(:,1); G3(:,1)]; Y [G1(:,2); G2(:,2); G3(:,2)]; fprintf(合并 r: %.2f\n, corr(X, Y)); % 约 0.92方向完全相反原因在于分组变量同时影响了两列数据的水平。组与组之间 x 在升高、y 也在升高这种组间差异贡献了巨大的协方差把组内的负向关系整个盖住了。解决路径有两条一是分层分组各算各的系数二是把分组标签当作协变量用partialcorr(X, Y, G)计算偏相关剔除分组变量的影响后再看剩余关系。这条坑之所以值得单独列出来是因为它不会报错也不会警告输出就是一个漂亮的 0.92直接拿来写结论几乎必错。5.2 现象删掉一个观测相关系数从 0.97 掉到 0.6数据基本沿着一条直线分布只有一个点的 y 值明显偏离。原因在于皮尔逊系数的分子分母都被离差平方驱动。一个远离均值的点在中心化后离差很大它自己贡献的乘积项和平方项都能主导整个求和相当于给这个点加了杠杆。表现上就是系数被单个样本绑架。解决办法分三层先看散点图这一步没有替代品再用斯皮尔曼看秩相关是否稳定如果秩相关仍在 0.9 以上说明是数值异常而非结构变化最后追查这个点的来源是录入错误就修正是真实极值就保留并在报告中说明。rng(1); x (1:30); y 0.8*x randn(30,1); r_clean corr(x, y) % 接近 0.99 y(30) 5; % 一个点被录错 r_dirty corr(x, y) % 明显下降 r_spear corr(x, y, Type, Spearman)rng(1)固定随机种子保证可复现。对比三行输出就能看出皮尔逊掉得最狠斯皮尔曼只受轻微影响因为那个点的秩只挪动了一位。这个对比本身就是一次快速的稳健性诊断。5.3 现象相关矩阵求逆报错或者回归系数大得离谱两个自变量几乎是同一个变量的两份测量相关矩阵里对应元素达到 0.999。原因在于相关矩阵接近奇异最小特征值趋近 0条件数爆炸。求逆运算会把这个微小的特征值放大成巨大的数值回归系数因此失去物理意义符号还可能反直觉。n 50; x1 randn(n,1); x2 x1 0.01*randn(n,1); % 几乎就是 x1 的复制品 x3 randn(n,1); R corrcoef([x1 x2 x3]); fprintf(cond(R) %.1f\n, cond(R)); % 条件数极大 vif diag(inv(R)) % VIF_j 1/(1-R_j^2)cond返回最大与最小奇异值之比经验上超过 30 就要留意上千基本可以断定存在强共线性。diag(inv(R))给出各变量的方差膨胀因子超过 10 通常视为共线性严重。解决思路包括删掉冗余变量、做 PCA 取主成分、加岭回归的正则项。别指望corr会给你提示它只负责算出那个 0.999剩下的判断得自己做。5.4 现象corr 返回一片 NaN索引全是 NaN数据里零星有几个缺失但整行整列的结果都成了 NaN。原因是corrcoef和corr的默认行为Rows,all会丢掉任何含缺失值的整行如果某列恰好每个缺失都对上了其他列就会出现整列无效。解决办法是显式指定Rows参数。X randn(20, 3); X(3,1) NaN; X(7,2) NaN; R_all corr(X, Rows, all); % 缺一个就丢整行有效样本从 20 掉到 18 R_pair corr(X, Rows, pairwise); % 每对用各自可用的样本 eig(R_pair) % 可能出现负特征值all保证所有变量对用的是同一批样本矩阵正定、可比pairwise尽量多用数据但不同变量对的样本集不同矩阵可能出现负特征值拿去求逆或者做因子分析会直接崩。我的习惯是默认用all只有在样本损失过半时才考虑pairwise并且用完之后检查最小特征值是否为正。5.5 现象两条互不相关的随机序列算出 r0.85两条独立生成的序列理论上毫无关系相关系数却高得惊人。原因在于两条序列都带趋势。累积求和得到的随机游走其值会长时间停留在高位或低位两列都存在这种漂移时样本内的「共同上升段」会被算成线性相关而有效样本量远小于名义样本量。这类结果叫伪相关时间序列上尤其常见。rng(11); w1 cumsum(randn(200,1)); w2 cumsum(randn(200,1)); r_level corr(w1, w2) % 可能高达 0.8 上下 r_diff corr(diff(w1), diff(w2)) % 差分后回到 0 附近cumsum把白噪声变成随机游走diff取一阶差分还原成增量序列。排查路径是先看去趋势后的结果再画自相关函数确认序列的持续性必要时用有效样本量代替名义 n 重算显著性。只盯着 r 的数值不看序列结构这个坑会反复出现。6. 从相关矩阵到建模决策偏相关、距离相关与验证取舍相关矩阵算出来只是起点真正决定建模质量的是怎么用它。偏相关解决的是「控制住其他变量之后这两个还相关吗」——partialcorr(x, y, z)把 z 对 x 和 y 分别回归再对残差求相关本质上是把 z 解释掉的部分剔除。多变量场景下原始相关矩阵里几乎每一对都显著只有偏相关才能告诉你哪条边是真实的。距离相关则补上了单调关系之外的缺口。它把样本间距离构造成双中心化矩阵再算两个矩阵的协方差型统计量对任何形式的依赖都能给出非零响应。上面那个y x²的例子皮尔逊约等于 0距离相关却能明显跳出 0。rng(3); x randn(40,1); y x.^2 0.5*randn(40,1); r_pear corr(x, y) % 接近 0但显然有关系 B 5000; r_perm zeros(B,1); for b 1:B r_perm(b) corr(x, y(randperm(numel(y)))); end p_perm mean(abs(r_perm) abs(r_pear)); % 置换检验的 p 值randperm打乱 y 的顺序切断任何真实关联重复 B 次得到零假设下的分布再用实际统计量的位置估 p 值。置换检验的好处是完全不依赖正态假设而且可以套在任何统计量上——换成距离相关、换成肯德尔代码结构一模一样。真正的取舍在这里皮尔逊高效可解释、能直接进线性模型但对非线性和异常值脆弱秩相关稳健却丢掉了幅度信息没法直接代入回归距离相关最灵活但没有方向性、没有现成的系数可解释还要配置换检验才能定显著性。我的判断标准是看下游要什么——如果相关分析只是为了筛选特征喂给回归优先皮尔逊加 VIF 检查如果是为了在报告里说明变量间的关联结构秩相关配合散点图更稳妥如果怀疑存在被线性系数掩盖的非线性依赖才把距离相关和置换检验搬出来成本高但能兜底。从那以后我每做一次相关性分析都会强制先走一遍四张图原始散点、秩散点、分组散点、时间序列的自相关。四张图过了再谈系数系数过了再谈 p 值顺序颠倒一次就得返工一次。希望帮到你。本文还有配套的精品资源点击获取