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); end

std(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=%.6f，t=%.4f\n', p, p_manual, t);

tcdf是 t 分布的累积分布函数，abs(t)取绝对值是因为双侧检验只看偏离零的幅度，最后乘 2 把两侧尾巴加起来。手算结果与corr输出应当一致到小数点后多位；如果对不上，八成是漏了abs或者自由度写成了 n-1。搞清这一层之后，p 值就不再是黑匣子里的魔法数字——它是 r 和 n 联合组装出来的一个概率，样本量一变，同样的 r 可以得出完全相反的结论。

3.2 为什么 n 只有 8 的时候 r=0.7 依然不显著

在小样本数学建模任务里，这种情况极其常见：跑出来 r=0.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 的上侧分位。n=8 时自由度是 6，临界 t 约为 2.447，对应的 r 临界值约 0.707。也就是说 r=0.70 恰好落在临界线内侧，p 值在 0.05 上方一点点，结论就是「没有足够证据认为相关」。同样的 r=0.70，如果 n=30，自由度 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.96，n-3 里的 3 同样来自自由度修正。手算区间与corr输出在小样本下会有细微差异，因为官方实现可能用了更精确的近似，但数量级和形态是一致的。置信区间的价值在于它直接暴露了估计的不确定性：r=0.75 配 n=12，区间可能宽到 [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 = i+1: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

内层循环遍历所有 i<j 的组合，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 三种系数的选择依据

选择路径我一般这么走：先画散点图，形状接近直线且没有明显的极端点，用皮尔逊；散点呈现单调曲线（指数、对数、饱和型），或者数据本身是等级、排名，用斯皮尔曼；样本量不大但并列值多、或者怀疑个别点会影响结论，用肯德尔并同时看 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 现象：两条互不相关的随机序列，算出 r=0.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 值，顺序颠倒一次就得返工一次。希望帮到你。

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