相关系数代码篇

相关系数代码篇
- 1 R实现
- 2 Python实现

相关系数是对变量间相关程度的度量（我好像又在讲废话了😅），具体概念详见相关系数概念篇。我最近发现vscode和Jupyter Lab好像都可以实现对python和R的编译（vscode已经实现但是好像容易卡掉，具体的我还需要多了解了解【⭐有坑待填】）。

1 R实现#

1.1 选定相关系数#

怎么选相关系数，放在之前我的标准就是：看心情；或者会算哪个来哪个；又或者都算了，随便选一个😅。那现在，本着严谨认真的科学态度，选择标准就进化成了：首先看样本数据是否满足正态分布且没有离群值，这两项都满足ne，就选Pearson；要是不满足，看样本数量，要是小于30个ne，就选Kendall，否则就选Spearman。

1.1.1 正态分布检验#

Q-Q图

Copy
library(ggpubr)
ggqqplot(data$runoff.mode, ylab = "径流模数")
ggqqplot(data$NDVI, ylab = "NDVI")

Shapiro-Wilk检验

Copy
shapiro.test(data$NDVI)
shapiro.test(data$sa)

其他的检验还有Kolmogorov-Smirnov（KS）检验

1.1.2 离群值检测#

outliers包参考【译文】R语言中的离群值检测和处理 - 简书 (jianshu.com)

1.2 计算#

R里边相关系数矩阵的计算就非常方便la

Copy
cor_ma = cor(data, method = c("pearson"))
cor_k = cor(data, method = c("kendall"))
cor_s = cor(data, method = c("spearman"))

cor.mtest <- function(mat, ...) 
{mat <- as.matrix(mat)
n <- ncol(mat)
p.mat<- matrix(NA, n, n)
diag(p.mat) <- 0
for (i in 1:(n - 1)){
  for (j in (i + 1):n){
    tmp <- cor.test(mat[, i], mat[, j], ...)
    p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
  }
}
colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
p.mat
}
p.mat <- cor.mtest(data)

1.3 绘图#

Copy
# Cairo制图输出----
# AOE 特征向量角序
library(Cairo)
CairoPNG(width = 3800, height = 3800,
file = "cor_pearson.png",
dpi = 600,
bg = "white")
# CairoPDF(file = "cor_spearman.pdf",
#          bg = "white")
corrplot(cor_ma, order = "AOE", type="upper",
         tl.pos = "d", tl.col="black",  #tl 属性名称设置
         tl.cex = 0.75, cl.pos = "r",   #cl 图例设置
         cl.ratio = 0.2, cl.align.text = "l",cl.cex = 1,
         cl.offset = 0.1, addgrid.col = "black",
         insig = c("label_sig"), p.mat = p.mat, sig.level=0.05,
         mar = c(0, 0, 0, 0),
         col = COL2('PRGn', 10))
corrplot(cor_ma, add = TRUE, type = "lower", 
         method = "number", number.cex = 1,
         order="AOE", diag=FALSE,
         addgrid.col = "black", col = "black",
         tl.pos="n", cl.pos="n")
dev.off()

2 Python实现#

2.1 选定相关系数#

2.1.1 正态分布检验#

正态分布检验就是通过样本数据检验总体是否呈现正态分布（好像是个废话...）（统计检验见【⭐有坑待填】）。正态分布检验方法包括1️⃣图形法(q-q plot和p-p plot)；2️⃣统计量法。scipy包中有很多正态检验的方法（参考用 Python 检验数据正态分布的几种方法 - 知乎 (zhihu.com)）：shapiro、kstest、anderson 、normaltest ，本文就选取shapiro做一个示范。

Shapiro-Wilk检验H0为：样本所来自的总体分布服从正态分布。当p值小于0.05时，基本可以说明该样本所来自的总体不服从正态分布。

Copy
from scipy import stats
re = []
for i in range(0,len(df.columns)):
    x = df.iloc[:, i]
    shapiro_test = stats.shapiro(x)
    re.append(shapiro_test)
print(re)

2.1.2 离群值检测#

离群值监测也可以通过图形法（箱型图）+Tukey fences进行初步判断。其他的判断方法有Isolation forest、Local outlier factor、DBScan、Elliptic envolope、Ensemble（参考Multivariate outlier detection in Python | by Philip Wilkinson | Towards Data Science）。看到DBScan，就会想到聚类，看到聚类就会想到距离【⭐有坑待填】。

2.2 计算#

计算这边就显得比较简单了，直接采用scipy中的相关系数计算函数

Copy
from scipy.stats import pearsonr, spearmanr, kendalltau
import pandas as pd
def corr_sci(df):
    corr_r_df =pd.DataFrame(index=df.columns, columns=df.columns)
    corr_p_df =pd.DataFrame(index=df.columns, columns=df.columns)
    for i in range(len(corr_p_df.columns)):
        for c in range(len(corr_p_df.columns)):
            # 下面的pearsonr可以变为kendalltau、spearmanr
            corr_r_df.iloc[i,c] = pearsonr(df[corr_r_df.index[i]], df[corr_r_df.columns[c]])[0]
            corr_p_df.iloc[i,c] = pearsonr(df[corr_p_df.index[i]], df[corr_p_df.columns[c]])[1]
    return corr_r_df,corr_p_df
df_r, df_p = corr_sci(df)
df_r = df_r.apply(pd.to_numeric,errors='ignore')
df_p = df_p.apply(pd.to_numeric,errors='ignore')

2.3 绘图#

绘图的主要思路可以先看正菜：就是采用seaborn中的heatmap包展示相关系数的颜色；采用掩膜将上三角的颜色抹掉，再用ax.text()在相应的位置加上属性名与相关系数。

首先，对热力图的上三角进行掩膜处理

Copy
import numpy as np
# Generate a mask for the upper triangle
mask = np.zeros_like(df_r, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

对绘图参数进行设置，Jupyter Notebook中加入%matplotlib widget食用感受更佳。

Copy
import matplotlib.pyplot as plt
%matplotlib widget
config={
    "font.family":'serif',
    "mathtext.fontset":'stix',
    "font.serif":['Times New Roman'],
    'axes.unicode_minus':False
}
plt.rcParams.update(config)

最后，上正菜...其中，配色可以参考网站tidyfriday_RStata

Copy
import matplotlib.colors as mcolors
import seaborn as sns
colors = ["#d7191c", "#fdae61", "#ffffbf", "#abd9e9", "#2c7bb6"]
# cmap = mcolors.ListedColormap(colors)
cmap1 = mcolors.LinearSegmentedColormap.from_list('cmap1', colors)
fig, ax = plt.subplots(figsize=(20,20),dpi=100)
# sns热力图
sns.heatmap(data = df_r, mask=mask, cmap=cmap1, vmin=-1, vmax=1,
            square=True, linecolor="black", linewidths=0.5, ax=ax,cbar = False)
# 展示属性名与相关系数数值
for i in range(len(df_r)):
    ax.text(i+0.5,(i+0.5), df_r.columns[i], 
            ha="center", va="center", rotation=0,fontdict={'family':'Times New Roman', 'size':13})
    for j in range(i+1, len(df_r)):
        s = "{:.2f}".format(df_r.values[i,j])
        ax.text(j+0.5,(i+0.5),s, 
            ha="center", va="center",fontdict={'family':'Times New Roman', 'size':16})
# p值小于等于0.05用*标记
for n in range(len(df_p)):
    for k in range(n+1, len(df_p)):
        l = "{:.3f}".format(df_p.values[len(df_p)-k-1,len(df_p)-n-1])
        if df_p.values[len(df_p)-k-1,len(df_p)-n-1] <= 0.05:
            ax.scatter(len(df_p)-(k+0.5),len(df_p)-(n+0.5),color="black",marker='*', s = 50)
# marker='o', 'v', '^', '<', '>', '8', 's', 'p', '*', 'h', 'H', 'D', 'd', 'P', 'X'
ax.axis("off")
cax = fig.add_axes([0.17, 0.06, 0.68, 0.028])#色标位置参数
cax.tick_params(labelsize=20)
fig.colorbar(ax.get_children()[0],orientation='horizontal',ax=ax,cax=cax)
plt.show()