Python与R实战:皮尔逊相关性假设检验全流程指南
当你拿到一组数据,兴奋地计算完皮尔逊相关系数准备发表结论时,是否想过这个数字可能完全失真?我曾在一个客户项目中犯过这个错误——在没有验证假设的情况下直接报告了0.85的高相关性,结果在复检时发现仅仅因为一个异常值就导致结论完全错误。本文将带你用代码完整走通皮尔逊相关性的五大假设检验流程,避开这些统计陷阱。
1. 环境准备与数据加载
工欲善其事,必先利其器。我们先配置好双语言环境:
# Python环境 import numpy as np import pandas as pd import seaborn as sns import matplotlib.pyplot as plt from scipy import stats# R环境 library(ggplot2) library(dplyr) library(car) # 用于QQ图假设我们有一组电商数据,包含用户浏览时长(分钟)和购买金额(元):
# Python生成模拟数据 np.random.seed(42) browse_time = np.random.normal(30, 5, 100) purchase_amt = 50 + 2*browse_time + np.random.normal(0, 10, 100) df = pd.DataFrame({'browse_time': browse_time, 'purchase_amt': purchase_amt})# R生成相同数据 set.seed(42) browse_time <- rnorm(100, 30, 5) purchase_amt <- 50 + 2*browse_time + rnorm(100, 0, 10) df <- data.frame(browse_time, purchase_amt)2. 变量类型验证
皮尔逊相关要求连续型数值变量。我们先做快速检查:
# Python类型检查 print(df.dtypes) # 输出应当显示: # browse_time float64 # purchase_amt float64# R类型检查 str(df) # 输出应显示: # 'data.frame': 100 obs. of 2 variables: # $ browse_time : num 31.4 29.8 32.2... # $ purchase_amt: num 107 116 128...注意:如果变量是分类数据,应该改用斯皮尔曼或肯德尔相关系数。例如客户评级(1-5星)与购买金额的相关性就不适合用皮尔逊。
3. 线性关系检验
散点图是最直观的检验工具。我们对比健康数据和存在非线性关系的数据:
# Python散点图 plt.figure(figsize=(12,5)) plt.subplot(121) sns.scatterplot(data=df, x='browse_time', y='purchase_amt') plt.title('健康数据') # 生成非线性数据示例 x = np.linspace(1, 20, 100) y = 50 + 0.5*x**2 + np.random.normal(0, 10, 100) plt.subplot(122) sns.scatterplot(x=x, y=y) plt.title('非线性数据');# R散点图 library(patchwork) p1 <- ggplot(df, aes(browse_time, purchase_amt)) + geom_point() + ggtitle("健康数据") x <- seq(1, 20, length.out=100) y <- 50 + 0.5*x^2 + rnorm(100, 0, 10) p2 <- ggplot(data.frame(x,y), aes(x,y)) + geom_point() + ggtitle("非线性数据") p1 + p2 # 并排显示当发现非线性模式时,解决方案包括:
- 尝试变量转换(对数、平方根等)
- 使用距离相关系数(dcor)等非线性相关度量
- 考虑回归模型而非简单相关
4. 正态性检验实战
正态性检验有可视化和统计检验两种方法。我们先看QQ图:
# Python正态检验 fig, axes = plt.subplots(1, 2, figsize=(12,5)) stats.probplot(df['browse_time'], plot=axes[0]) axes[0].set_title('浏览时长QQ图') stats.probplot(df['purchase_amt'], plot=axes[1]) axes[1].set_title('购买金额QQ图');# R正态检验 par(mfrow=c(1,2)) qqPlot(df$browse_time, main="浏览时长QQ图") qqPlot(df$purchase_amt, main="购买金额QQ图")统计检验推荐Shapiro-Wilk检验(样本量<5000):
# Python Shapiro检验 print("浏览时长检验:", stats.shapiro(df['browse_time'])) print("购买金额检验:", stats.shapiro(df['purchase_amt'])) # 输出示例: # ShapiroResult(statistic=0.992, pvalue=0.892) # p>0.05则接受正态假设# R Shapiro检验 shapiro.test(df$browse_time) shapiro.test(df$purchase_amt)当数据非正态时的处理方案:
| 情况 | 解决方案 |
|---|---|
| 轻度偏态 | 可继续使用皮尔逊(鲁棒性较强) |
| 严重非正态 | 改用斯皮尔曼相关系数 |
| 有明确分布 | 考虑基于特定分布的相关度量 |
5. 异常值检测与处理
异常值会显著扭曲相关系数。我们用三种方法检测:
IQR方法(最常用):
# Python异常值检测 def detect_outliers(series): Q1 = series.quantile(0.25) Q3 = series.quantile(0.75) IQR = Q3 - Q1 return ~((series < (Q1 - 1.5*IQR)) | (series > (Q3 + 1.5*IQR))) clean_df = df[detect_outliers(df['browse_time']) & detect_outliers(df['purchase_amt'])]# R异常值检测 is_outlier <- function(x) { Q <- quantile(x, probs=c(0.25, 0.75)) iqr <- IQR(x) !(x < (Q[1] - 1.5*iqr) | x > (Q[2] + 1.5*iqr)) } clean_df <- df[is_outlier(df$browse_time) & is_outlier(df$purchase_amt), ]可视化确认:
# Python箱线图 plt.figure(figsize=(10,4)) plt.subplot(121) sns.boxplot(data=df[['browse_time']]) plt.subplot(122) sns.boxplot(data=df[['purchase_amt']]);马氏距离(适用于多元异常值):
# Python马氏距离 from scipy.spatial.distance import mahalanobis cov = np.cov(df.values.T) inv_cov = np.linalg.inv(cov) mean = np.mean(df, axis=0) df['mahalanobis'] = [mahalanobis(x, mean, inv_cov) for x in df.values]处理异常值的策略优先级:
- 检查是否为数据录入错误
- 考虑业务背景判断是否保留
- 使用稳健统计量(如中位数代替均值)
- 进行变量转换
- 最后才考虑删除
6. 完整检验流程与自动化
我们将所有检验整合为可复用的函数:
def pearson_assumptions_check(df, var1, var2, alpha=0.05): """ 皮尔逊假设全自动检验 """ results = {} # 1. 变量类型检查 results['dtypes'] = df[[var1, var2]].dtypes.to_dict() # 2. 线性检验 r, p = stats.pearsonr(df[var1], df[var2]) results['linearity'] = {'r': r, 'p': p} # 3. 正态检验 shapiro1 = stats.shapiro(df[var1]) shapiro2 = stats.shapiro(df[var2]) results['normality'] = { var1: {'statistic': shapiro1[0], 'pvalue': shapiro1[1]}, var2: {'statistic': shapiro2[0], 'pvalue': shapiro2[1]} } # 4. 异常值检测 clean_df = df[detect_outliers(df[var1]) & detect_outliers(df[var2])] results['outliers'] = { 'original_size': len(df), 'clean_size': len(clean_df), 'removed': len(df) - len(clean_df) } # 综合判断 assumptions_met = all([ p < alpha, # 线性显著 shapiro1[1] > alpha and shapiro2[1] > alpha, # 正态性 results['outliers']['removed']/len(df) < 0.05 # 异常值<5% ]) results['assumptions_met'] = assumptions_met return resultspearson_assumptions_check <- function(df, var1, var2, alpha=0.05) { results <- list() # 1. 变量类型检查 results$dtypes <- sapply(df[c(var1, var2)], class) # 2. 线性检验 cor_test <- cor.test(df[[var1]], df[[var2]]) results$linearity <- list( r = cor_test$estimate, p = cor_test$p.value ) # 3. 正态检验 shapiro1 <- shapiro.test(df[[var1]]) shapiro2 <- shapiro.test(df[[var2]]) results$normality <- list( var1 = list(statistic = shapiro1$statistic, p.value = shapiro1$p.value), var2 = list(statistic = shapiro2$statistic, p.value = shapiro2$p.value) ) # 4. 异常值检测 clean_df <- df[is_outlier(df[[var1]]) & is_outlier(df[[var2]]), ] results$outliers <- list( original_size = nrow(df), clean_size = nrow(clean_df), removed = nrow(df) - nrow(clean_df) ) # 综合判断 results$assumptions_met <- all( cor_test$p.value < alpha, shapiro1$p.value > alpha && shapiro2$p.value > alpha, results$outliers$removed/nrow(df) < 0.05 ) return(results) }使用示例:
report = pearson_assumptions_check(df, 'browse_time', 'purchase_amt') print(report)最终决策流程可总结为:
如果所有假设满足:
- 报告皮尔逊相关系数及p值
- 建议补充95%置信区间
如果部分假设不满足:
- 考虑变量转换后重新检验
- 改用非参数方法(斯皮尔曼)
- 报告时注明限制条件
如果主要假设严重违反:
- 放弃相关性分析
- 考虑其他分析方法(回归、分类等)