相関の希薄化 - データサイエンス関連+αのメモ

相関の希薄化（attenuation of correlation） とは、データの測定誤差によって相関がゼロに近づく方向のバイアスをもつ現象。

信頼性は例えばテストやアンケートの設計に起因するものもあるし、連続変数を離散化したことによる情報の欠落によるものもある（→カテゴリカル変数の相関係数は希薄化を補正する考えから生まれた）

測定誤差による相関の希薄化¶

attenuation formula¶

相関の希薄化（attenuation）は Spearman (1904) によって発見された。またこの問題を修正するための式も提案された。

式を変形すると、「観測される相関は真の相関に信頼性の影響が掛かったものである」という考え方をする。

信頼性とは¶

再検査法（test-retest）的な考え方¶

例えばテストの得点においては「偶然高い得点（低い得点）になった」というようなランダム性の要素が少なくて、仮に同じ能力の人が同様のテストを繰り返し受験したら毎回同様の能力値の評価がなされるようなテストは信頼性が高い。逆にランダム性の影響が高いテストは信頼性が低いと考えられる。

なお、この考え方はテスト理論における再検査法の考え方。

再検査法的な考え方

r_{x x} = \operatorname{cor}(x_1, x_2)

(3)

$x_1$ ：1回目のテストの測定結果
$x_2$ ：2回目のテストの測定結果

なお Spearman (1904) も信頼性の部分は「同じ確率変数の独立に得られたサンプルの系列における相関係数の平均」という考え方をとっていたし、1940~50年代の文献でもこうした考え方が支持されていた。

Spearman (1904)の書き方

To ascertain the amount of this attenuation, and thereby discover the true correlation, it appears necessary to make two or more independent series of observations of both $p$ and $q$ . Then,

r_{pq} = \frac{r_{p'q'}}{\sqrt{r_{p'p'} \cdot r_{q'q'}}}

(4)

where

$r_{p'q'}$ = the mean of the correlations between each series of values obtained for $p$ with each series obtained for $q$ .

$r_{p'p'}$ = the average correlation between one and another of these several independently obtained series of values for $p$ .

$r_{q'q'}$ = the same as regards $q$ .

and $r_{pq}$ = the required real correlation between the true objective values of $p$ and $q$ .

信頼性係数¶

本質的に同様の考え方ではあるが、観測値が「真のスコア」と「測定誤差」から構成されて、真のスコアの割合が高いほど信頼性が高いという考えから次のように定式化することもできる

古典的テスト理論における信頼性係数

r_{x x}=\frac{\operatorname{Var}(T)}{\operatorname{Var}(X)}=1-\frac{\operatorname{Var}(\epsilon)}{\operatorname{Var}(X)}

(5)

$T$ ：真のスコア
$\epsilon$ ：測定誤差
$X=T+\epsilon$ ：観測値

離散化による相関の希薄化¶

連続変数を離散化することは測定の精度が下がることに相当する。例えば、ある連続変数を「0.5以上かどうか」で二値に離散化した場合、0.49も-999も「閾値である0.5より低い」という点では同じ値として観測されることになる。

こうした点では測定誤差の相関の希薄化と類似した問題が起こる。そのため離散化された質的変数を含む2変数間の相関関係をピアソンの積率相関係数で評価しようとすると希薄化が起こって過小評価バイアスが含まれることになる。

希薄化の例¶

Peters & Van Voorhis (1940) は相関係数が $\rho$ の 2 変量正規分布に従う確率変数 $X, Y$ があるとき、 $X$ を平均値や中央値で二値化した確率変数を $X_d$ とすると、 $X_d$ と $Y$ の間の相関係数は $0.798 \rho$ になる、つまり真の相関係数の約 0.8 倍へと過小評価する問題があることを報告している。0.798になるのは二値化のときであり、離散化のカテゴリ数や分割点を変えれば希薄化の強さは変わる。

モンテカルロシミュレーションしてみると、だいたいそのような結果になることが簡単に確認できる。例えば真の相関係数 $\rho=0.5$ の分布から疑似乱数でサンプルを生成して二値化して相関係数 $\hat{\rho}=\operatorname{cor}(x_d, y)$ を計算することを何回も繰り返してその標本分布を描くと、平均的には $\operatorname{avg}(\hat{\rho}) \approx 0.4$ となり、その比率 $\frac{\hat{\rho}}{\rho}$ の平均は

\operatorname{avg}\left( \frac{\hat{\rho}}{\rho} \right) \approx 0.798

(6)

となることがわかる。

# 1回だけテスト
import numpy as np
from scipy.stats import pearsonr

# 正規分布からのサンプリング
rho = 0.5 # 真の相関係数
np.random.seed(0)
data = np.random.multivariate_normal(mean=[1, 5], cov=[[1, rho], [rho, 1]], size=100_000)
x, y = data[:, 0], data[:, 1]

# 平均値を閾値に二値化
yd = 1 * (y >= y.mean())
r = pearsonr(x, yd)[0]  # 二値化後の相関係数
print(f"真の相関={rho:.3f}, 二値化後の相関係数={r:.3f}, 比率={r/rho:.3f}")

真の相関=0.500, 二値化後の相関係数=0.399, 比率=0.798

# 複数回テスト
import numpy as np
from scipy.stats import pearsonr

def simulate_dichotomization(rho=0.5, size=100000, seed=0):
    np.random.seed(seed)
    data = np.random.multivariate_normal(mean=[0, 0], cov=[[1, rho], [rho, 1]], size=size)
    x, y = data[:, 0], data[:, 1]
    xd = 1 * (x >= np.mean(x))  # 平均で2値化
    r = pearsonr(xd, y)[0]
    return r

n = 1000
n_trials = 1500
true_rho = 0.5
estimates = []
for seed in range(n_trials):
    est_rho = simulate_dichotomization(rho=true_rho, size=n, seed=seed)
    estimates.append(est_rho)
estimates = np.array(estimates)


import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib_fontja
fig, axes = plt.subplots(ncols=2, figsize=[8, 3])

# fig1
hist = sns.histplot(estimates, ax=axes[0], color="#7abcff", edgecolor="#7abcff", kde=True, stat="density")
axes[0].set(title="積率相関係数の分布")

max_height = max([p.get_height() for p in hist.patches])
axes[0].axvline(true_rho, color="darkorange", linestyle="--", alpha=0.7)
axes[0].text(true_rho, max_height * 0.95, f"真のρ={true_rho:.1f} ", color="darkorange", alpha=0.9, ha="right")

r = estimates.mean()
axes[0].axvline(r, color="#007fff", linestyle=":", alpha=0.7)
axes[0].text(r, max_height * 0.95, f"平均={r:.3f} ", color="#007fff", alpha=0.9, ha="right")

# fig2
ratios = estimates / true_rho
hist = sns.histplot(ratios, ax=axes[1], color="#80f000", edgecolor="#80f000", kde=True, stat="density")
max_height = max([p.get_height() for p in hist.patches])
axes[1].set(title=r"比率 $\hat{\rho} / \rho$ の分布")

r = ratios.mean()
axes[1].axvline(r, color="#006400", linestyle=":", alpha=0.7)
axes[1].text(r, max_height * 0.95, f" 平均={r:.3f}", color="#006400", alpha=0.9, ha="left")

fig.show()

真の相関係数が異なる場合でも同様に一貫して0.798である

# 複数回テスト
import numpy as np
from scipy.stats import pearsonr

def simulate_dichotomization(rho=0.5, size=100000, seed=0):
    np.random.seed(seed)
    data = np.random.multivariate_normal(mean=[0, 0], cov=[[1, rho], [rho, 1]], size=size)
    x, y = data[:, 0], data[:, 1]
    xd = 1 * (x >= np.mean(x))  # 平均で2値化
    r = pearsonr(xd, y)[0]
    return r

import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib_fontja

def plot_hists(estimates, axes, set_title=True):
    # fig1
    hist = sns.histplot(estimates, ax=axes[0], color="#7abcff", edgecolor="#7abcff", kde=True, stat="density")

    max_height = max([p.get_height() for p in hist.patches])
    axes[0].axvline(true_rho, color="darkorange", linestyle="--", alpha=0.7)
    axes[0].text(true_rho, max_height * 0.95, f"真のρ={true_rho:.2f} ", color="darkorange", alpha=0.9, ha="right")

    r = estimates.mean()
    axes[0].axvline(r, color="#007fff", linestyle=":", alpha=0.7)
    axes[0].text(r, max_height * 0.95, r" $r$の"+ f"平均={r:.3f} ", color="#007fff", alpha=0.9, ha="right")

    # fig2
    ratios = estimates / true_rho
    hist = sns.histplot(ratios, ax=axes[1], color="#80f000", edgecolor="#80f000", kde=True, stat="density")
    max_height = max([p.get_height() for p in hist.patches])

    r = ratios.mean()
    axes[1].axvline(r, color="#006400", linestyle=":", alpha=0.7)
    axes[1].text(r, max_height * 0.95, r" $r$の"+ f"平均={r:.3f}", color="#006400", alpha=0.9, ha="left")

    if set_title:
        axes[0].set(title="二値化後の積率相関係数の分布")
        axes[1].set(title=r"比率 $r / \rho$ の分布")

n = 10000
n_trials = 10000
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=[8, 9])
for i, true_rho in enumerate([-0.5, 0.25, 0.75]):
    # シミュレーション
    estimates = []
    for seed in range(n_trials):
        est_rho = simulate_dichotomization(rho=true_rho, size=n, seed=seed)
        estimates.append(est_rho)
    estimates = np.array(estimates)
    # plot
    plot_hists(estimates, axes[i, :], set_title = i == 0)
    fig.show()

希薄化の係数¶

Cohen (1983)によれば、Peters and Van Voorhis (1940)が報告しているような二値化の希薄化の係数は次のように一般化できる

e = \frac{h}{\sqrt{p(1 - p)}}

(7)

$h$ ：標準正規分布の密度関数の縦軸（the ordinate of the standard unit normal curve）
$p$ ：二値化した変数の比率

例えば標準正規分布に従う連続変数 $X$ を平均値で二値化した場合は分割点が0なので、標準正規分布の密度関数

\phi(z)=\frac{1}{\sqrt{2 \pi}} e^{-z^2 / 2}

(8)

より、

$h = \phi(0) \approx 0.3989$
$p = 0.5$

となり、

e = \frac{h}{\sqrt{p(1 - p)}} = \frac{0.3989}{\sqrt{ 0.5 \cdot 0.5}} = \frac{0.3989}{0.5} \approx 0.798

(9)

となる。

from scipy.stats import norm
p = 0.5
norm.pdf(0) / np.sqrt(p * (1 - p))

np.float64(0.7978845608028654)

biserial相関係数¶

なお、二値化変数と連続変数の間の相関係数の一種である Biserial相関係数 はCohen (1983)が述べたような希薄化の係数の逆数をピアソンの積率相関係数に掛けて補正する（Spearmanの希薄化の修正式を使う）タイプの相関係数である

r_{\text{bi}} = \frac{ r_{\text{Pearson}} }{ e } = r_{\text{Pearson}} \cdot \frac{\sqrt{p (1-p)}}{h}

(10)

# 二変量標準正規分布からのサンプルとする
import numpy as np
n = 10000
rho = 0.5
np.random.seed(0)
data = np.random.multivariate_normal(mean=[0, 0], cov=[[1, rho], [rho, 1]], size=n)
x, y = data[:, 0], data[:, 1]

# 平均値を閾値として二値化されたとする
yd = 1 * (y >= np.mean(y))  # 二値化

# Cohen (1983) にも記載の希薄化率 e
from scipy.stats import norm
p = np.mean(yd)  # 離散化後に y_d == 1 の比率
z = norm.ppf(p)  # 閾値の位置（= y_d == 1の比率、の標準正規分布上の値）
h = norm.pdf(z)  # 標準正規分布の閾値での確率密度
e = h / np.sqrt(p * (1 - p))

# biserial 相関係数
from scipy.stats import pearsonr
r = pearsonr(x, yd).statistic
r_bi = (1 / e) * r

print(f"真の相関ρ={rho:.3f}, 観測値r={r:.3f}, 誤差e={e:.3f}, 補正後(biserial)={r_bi:.3f}")

真の相関ρ=0.500, 観測値r=0.396, 誤差e=0.798, 補正後(biserial)=0.496

Bedrick (1995) は離散化の仕方や離散化前の潜在変数が従う分布などの条件によって希薄化の強さが異なることを報告している。

参考文献¶

Spearman, C. (1904). The proof and measurement of association between two things. Am J Psychol, 15, 72-101.
Bedrick, E. J. (1995). A note on the attenuation of correlation. British Journal of Mathematical and Statistical Psychology, 48(2), 271-280.
Cohen, J. (1983). The cost of dichotomization. Applied psychological measurement, 7(3), 249-253.
Drasgow, F. (1986). Polychoric and polyserial correlations. In S. Kotz & N. Johnson (Eds.), The Encyclopedia of Statistics.
Peters, C. C., & Van Voorhis, W. R. (1940). Statistical procedures and their mathematical bases. McGraw-Hill.
Trafimow, D. (2016). The attenuation of correlation coefficients: A statistical literacy issue. Teaching Statistics, 38(1), 25-28.
- テスト理論的な観点から説明
Muchinsky, P. M. (1996). The correction for attenuation. Educational and psychological measurement, 56(1), 63-75.
- 歴史的な経緯も含め説明