偏相关（Partial Correlation）
偏相关是一种统计度量，用于衡量两个变量之间的相关性，同时控制一个或多个其他变量的影响。在多元变量分析中，偏相关可以帮助我们了解变量间的直接关系，而不受其他变量的干扰。
偏相关的数学定义
设有三个变量、和，偏相关系数 是控制了的条件下，和之间的相关程度。其计算公式如下：
其中：
是变量和之间的相关系数。 是变量和之间的相关系数。 是变量和之间的相关系数。
相关系数（Correlation Coefficient）的计算公式：
和 代表两个变量。 和 分别代表 和 的平均值。 代表 和 之间的相关系数。 代表样本数量。
偏相关（Partial Correlation）：
、 和 分别代表变量和、和、和之间的相关系数。 代表当作为控制变量时，X1和X2之间的偏相关系数。
APO与ENSO
2006年，赵平和张人禾用E0F分析研究东亚一太平洋区域海平面气压场时注意到在东亚与太平洋海平面气压之间存在个大尺度纬向遥相关，并称为东亚一太平洋偶极子 。也称亚洲太平洋涛动（Asian-Pacific Oscillation, APO），反映了亚洲大陆和太平洋地区大气压系统之间的相互作用和变化。APO对季风、降水和气候变化有着重要影响。
APO 的变化影响亚洲季风和太平洋地区的气候。正位相的 APO 通常与较强的亚洲夏季季风相关联，从而导致更多的降水，尤其是在东亚和南亚地区。 同时，APO 也影响太平洋地区的气候，夏季APO与太平洋年代际涛动(PDO) 和 厄尔尼诺-南方涛动(ENSO) 关系密切，例如它可以调节ENSO的强度。
Pingouin
这里我们计算偏相关使用Pingouin库实现，Pingouin是一个用 Python 3 编写的开源统计包，主要基于 Pandas 和 NumPy，目前已针对 Python 3.8-3.11 进行了测试。
可以使用 pip 轻松安装：
pip install pingouin
或者conda:
conda install -c conda-forge pingouin
Pingouin.partial_corr 参数说明
data: Pandas DataFrame 数据集。x,y: 字符串，指定需要计算偏相关的两个变量的列名。covar: 字符串或列表，指定一个或多个协变量的列名，这些协变量在计算偏相关时会被控制。x_covar: 字符串或列表，指定仅从 x 变量中去除影响的协变量。使用此参数时，即计算半偏相关。y_covar: 字符串或列表，指定仅从 y 变量中去除影响的协变量。使用此参数时，即计算半偏相关。alternative: 字符串，定义备择假设的尾部，可以是"two-sided"（默认）、"greater"或"less"。method: 字符串，指定相关性类型，可以是"pearson"（皮尔逊相关系数）或"spearman"（斯皮尔曼相关系数）。
返回值
返回一个包含以下内容的 Pandas DataFrame：
'n': 样本大小(去除缺测值后)。'r': 偏相关系数。'CI95': 95% 置信区间。'p-val': p 值。
注意事项
函数会自动去除含有缺失值的行。 偏相关系数的取值范围是 -1 到 1，表示完美的负相关到完美的正相关。 半偏相关与偏相关类似，但控制变量集仅从 x 或 y 中去除，而不是两者都去除。
计算NetCDF数据空间偏相关
这里我们使用上面介绍的夏季APO和ENSO指数数据，计算APO与海温(ERSSTv5) 的偏相关(去除掉ENSO信号)。
读取数据
import xarray as xr
import numpy as np
from matplotlib import pyplot as plt
import cartopy.crs as ccrs
from tqdm import tqdm
import pandas as pd
import pingouin as pg
import warnings
warnings.filterwarnings('ignore', category=RuntimeWarning)
def createFigure(figsize=(12, 8), dpi=300, subplotAdj=None, **kwargs):
figsize = figsize
figure = plt.figure(figsize=figsize, dpi=dpi, **kwargs)
if subplotAdj is not None:
plt.subplots_adjust(**subplotAdj)
return figure
apo = xr.open_dataset("apo.nc")['apo0'].values
month = [6, 7, 8]
nino = xr.open_dataset("NINO3.4_index.nc").where(
lambda m: m.time.dt.month.isin(month), drop=True).groupby(
"time.year").mean(
"time")['NINO3.4'].sel(year = slice(1980, 2014)).values
sst = xr.open_dataset("sst.mnmean.nc").where(
lambda m: m.time.dt.month.isin(month)).groupby(
"time.year").mean(
"time")["sst"].sel(year = slice(1980, 2014))
计算偏相关
lat, lon = sst.lat.values, sst.lon.values
sst = sst.values
r_matrix, p_matrix = [np.zeros_like(sst[0]) for _ in range(2)]
for i in tqdm(range(len(lat)), desc = 'Calc apo and sst partial corr...'):
for j in range(len(lon)):
if np.isnan(sst[:, i, j]).any():
continue
else:
df = pd.DataFrame({
'apo': apo,
'Nino3_4': nino,
'sst': sst[:, i, j]
})
pc = pg.partial_corr(data=df, x='sst', y='apo', covar='Nino3_4', method='pearson')
r_matrix[i, j] = pc['r'][0] # Partial correlation between IOB and precip
p_matrix[i, j] = pc['p-val'][0] # p_value for partial correlation
p_matrix[p_matrix == 0] = np.nan
绘图
fig = createFigure(figsize=(12, 8), dpi=300,
subplotAdj=dict(left=0.04, right=0.98,
top=0.9, bottom=0.05,
wspace=0.05, hspace=0.1))
ax = plt.subplot(111, projection=ccrs.Robinson(central_longitude=180))
ax.set_global()
ax.coastlines()
CF = plt.contourf(lon, lat,
r_matrix,
cmap= "RdBu_r",
levels=np.linspace(-0.5, 0.5, 21),
extend='both',
transform=ccrs.PlateCarree())
c1b = ax.contourf(lon, lat,
p_matrix,
levels=[0,0.05,0.5],
zorder=1,
hatches=['..', None],
colors="none", transform=ccrs.PlateCarree())
plt.title('Partial Correlation between APO and SST',
fontweight='bold',
fontsize=20)
position= fig.add_axes([0.31, 0.05, 0.4, 0.025])
fig.colorbar(CF,cax=position,orientation='horizontal',format='%.1f',)
plt.show()
Output
完整代码
import xarray as xr
import numpy as np
from matplotlib import pyplot as plt
import cartopy.crs as ccrs
from tqdm import tqdm
import pandas as pd
import pingouin as pg
import warnings
warnings.filterwarnings('ignore', category=RuntimeWarning)
def createFigure(figsize=(12, 8), dpi=300, subplotAdj=None, **kwargs):
figsize = figsize
figure = plt.figure(figsize=figsize, dpi=dpi, **kwargs)
if subplotAdj is not None:
plt.subplots_adjust(**subplotAdj)
return figure
apo = xr.open_dataset("apo.nc")['apo0'].values
month = [6, 7, 8]
nino = xr.open_dataset("NINO3.4_index.nc").where(
lambda m: m.time.dt.month.isin(month), drop=True).groupby(
"time.year").mean(
"time")['NINO3.4'].sel(year = slice(1980, 2014)).values
sst = xr.open_dataset("sst.mnmean.nc").where(
lambda m: m.time.dt.month.isin(month)).groupby(
"time.year").mean(
"time")["sst"].sel(year = slice(1980, 2014))
lat, lon = sst.lat.values, sst.lon.values
sst = sst.values
r_matrix, p_matrix = [np.zeros_like(sst[0]) for _ in range(2)]
for i in tqdm(range(len(lat)), desc = 'Calc apo and sst partial corr...'):
for j in range(len(lon)):
if np.isnan(sst[:, i, j]).any():
continue
else:
df = pd.DataFrame({
'apo': apo,
'Nino3_4': nino,
'sst': sst[:, i, j]
})
pc = pg.partial_corr(data=df, x='sst', y='apo', covar='Nino3_4', method='pearson')
r_matrix[i, j] = pc['r'][0] # Partial correlation between IOB and precip
p_matrix[i, j] = pc['p-val'][0] # p_value for partial correlation
p_matrix[p_matrix == 0] = np.nan
fig = createFigure(figsize=(12, 8), dpi=300,
subplotAdj=dict(left=0.04, right=0.98,
top=0.9, bottom=0.05,
wspace=0.05, hspace=0.1))
ax = plt.subplot(111, projection=ccrs.Robinson(central_longitude=180))
ax.set_global()
ax.coastlines()
CF = plt.contourf(lon, lat,
r_matrix,
cmap= "RdBu_r",
levels=np.linspace(-0.5, 0.5, 21),
extend='both',
transform=ccrs.PlateCarree())
c1b = ax.contourf(lon, lat,
p_matrix,
levels=[0,0.05,0.5],
zorder=1,
hatches=['..', None],
colors="none", transform=ccrs.PlateCarree())
plt.title('Partial Correlation between APO and SST',
fontweight='bold',
fontsize=20)
position= fig.add_axes([0.31, 0.05, 0.4, 0.025])
fig.colorbar(CF,cax=position,orientation='horizontal',format='%.1f',)
plt.show()
References
赵平,代玮 & 肖子牛.(2011).亚洲—太平洋涛动研究进展.气象科技进展(02),6-10.
