import pandas as pd
import numpy as np
from scipy.stats import t,scoreatpercentile,norm
import matplotlib.pyplot as plt
%matplotlib inline
5.2 平均への回帰¶
5.2.1 概要¶
- 第1回の授業で取り上げた平均への回帰という現象に関してもう少し詳しく見ていきます。
- 平均への回帰は、確率的にある程度の幅で変動する値を測定する際に観察される現象です。
- 血圧などが最も良い例です。
- 教育の文脈では、TOEICやTOEFLなど問題が違っても能力が同じならだいたい同じ点数が得られるようなテストを2回受験した場合などが考えられます。
- ですが、例えば、ある一定の期間に100個の漢字を授業で学んでテストAを受験、別の時期に同じ長さの期間で異なる漢字を100個授業で学んだ場合も平均への回帰は見ることができそうです。
- テストの場合は同じ問題を出題することができないので、「似たような」問題で同じ能力を測定していると考えます。
5.2.2 導出¶
- 回帰分析を思い出してください。
$$ y = ax + b$$
- 上式の傾き($a$)、切片($b$)は以下のように求めることができました。
$$ a = \frac{\sigma_{xy}}{\sigma^2_x}, b = \bar{y} - a \bar{x} $$
$a$は$x$と$y$の共分散を$x$の分散で割った値、$b$は$y$の平均値から$x$を$a$倍した値を引いた値でした。
$\bar{y}$を用いて$ y = ax + b$を表してみます。
$$ y = ax + \bar{y} - a \bar{x} $$
- 係数$a$でまとめると、いかになります。
$$ y = \bar{y} + a (x - \bar{x}) $$
- $x$と$y$の相関係数$ r_{xy}$は以下の式で求められます。
$$ r_{xy} = \frac{\sigma_{xy}}{\sqrt{\sigma^2_x} \times \sqrt{\sigma^2_y}} $$
- ここで、$\sigma^2_x = \sigma^2_y$と仮定すると、相関係数は以下になります。
$$ r_{xy} = \frac{\sigma_{xy}}{\sqrt{\sigma^2_x} \times \sqrt{\sigma^2_x}} $$
$$ = \frac{\sigma_{xy}}{{\sigma^2_x}} $$
- これは傾き$a$です。$ y = \bar{y} + a (x - \bar{x}) $は以下のように書くことができます。
$$ y = \bar{y} + r_{xy} (x - \bar{x}) $$
- $\bar{y}$を左辺に移行すると以下が得られます。
$$ y - \bar{y} = r_{xy} (x - \bar{x}) $$
- $r_{xy}$は１から-1の間の値を取り、1になるのは稀ですので、多くの場合、以下の関係が言えます。
$$ |y - \bar{y}| < |x - \bar{x}| $$
5.2.3 平均への回帰とカットオフスコア¶
- 以下の論文が提案するモデルです。
Davis, C. E. （1976）．The effect of regression to the mean in epidemiologic and clinical studies. American Journal of Epidemiology, 104（5）, 493–498.
- ここではある標準化されたテストを２回受験した場合を考えます。
- 1回目ののスコアを$y_1$とし、２回目のスコアを$y_2$とする。
- 1回目、2回目のテストともに母平均$\mu$、母分散$\sigma^2$の正規分布に従うこととする。
- $\rho$を2回のテストスコアの相関係数とする。
- ここで、$k$をカットオフスコアとし、$k$を$z = \frac{|k - \mu|}{\sigma}$と標準化すると、1回目のテストにおいて$k$未満のスコアの平均値は以下となる。
$$ E[y_1 | y_1 < k] = \mu - C\sigma $$
- 同様に、$k$を超えた受験者の1回目のスコアの平均値は以下となる。
$$ E[y_1 | y_1 > k] = \mu + C\sigma $$
ここで、
$$ C = \frac{\phi (z)}{1- \Phi (z)} $$
$$\phi (z) = \frac{1}{\sqrt{2\pi}}exp^{-\frac{z^2}{2}}$$
$$\Phi (z) = \int_{\infty}^{z} \phi(x)dx $$
である。
また、これらの受験者の2回目($y_2$)の平均値は、それぞれ
$$ E[y_2 | y_1 < k] = \mu - \rho C\sigma $$
$$ E[y_2 | y_1 > k] = \mu + \rho C\sigma $$
- 母平均、母分散、母相関係数が未知の場合、それぞれ、平均、不偏分散、相関係数を用いる。
5.2.4 平均への回帰とカットオフスコアの例¶
- 以下のデータはあるテストの結果で、1回目のテスト(X)で100点以下であった受験者は特別プログラムを受講した。
data = pd.read_csv("../DATA01/IEDA_reg.csv",index_col=0)
data.describe()
| X | Y | |
|---|---|---|
| count | 30.000000 | 30.000000 |
| mean | 106.333333 | 110.666667 |
| std | 23.958777 | 21.684785 |
| min | 60.000000 | 70.000000 |
| 25% | 90.000000 | 96.250000 |
| 50% | 105.000000 | 105.000000 |
| 75% | 128.750000 | 130.000000 |
| max | 145.000000 | 150.000000 |
- 以下はこのデータのヒストグラムである。
data.plot(kind="hist",alpha=0.5)
<AxesSubplot: ylabel='Frequency'>
- 以下の散布図で、x軸は１回目のテストのスコア、y軸は２回目のテストのスコアを示している。
- 赤の実線がこのデータの回帰直線、赤の点線が$y = x$のグラフである。
- ここで示されている実線と点線の違いが平均への回帰の大きさである。
x = data["X"].values
y = data["Y"].values
x_i = np.arange(60,160,1)
a,b = np.polyfit(x,y,1)
y_hat = a*x_i +b
plt.scatter(x,y)
plt.plot(x_i,y_hat,color="red")
y_i = x_i
plt.plot(x_i,y_i,linestyle="dotted",color="red")
[<matplotlib.lines.Line2D at 0x179190460>]
- ここで、１回目のスコアと２回目のスコアの差得点(diff)を算出し、x軸を１回目のスコア、y軸を差としてプロットする。
- 赤の実線は回帰直線である。
- １回目のスコアと差得点には負の相関がある。
data["diff"] = data["Y"] - data["X"]
d = data["diff"]
x_i = np.arange(60,160,1)
a,b = np.polyfit(x,d,1)
y_hat = a*x_i +b
plt.scatter(x,d)
plt.plot(x_i,y_hat,color="red")
[<matplotlib.lines.Line2D at 0x1792177c0>]
- 1回目のスコアが平均より低い人は２回目のスコアが高い傾向がある。
- 差得点と１回目のスコアには負の相関がある。
data.corr()
| X | Y | diff | |
|---|---|---|---|
| X | 1.000000 | 0.899224 | -0.425399 |
| Y | 0.899224 | 1.000000 | 0.013400 |
| diff | -0.425399 | 0.013400 | 1.000000 |
- １回目で100点以下の受験者の１回目の平均値
data[data["X"] < 101].mean()[0]
83.84615384615384
- １回目で100点以下の受験者の2回目の平均値
data[data["X"] < 101].mean()[1]
91.53846153846153
5.2.5 カットオフスコアを用いた場合の平均への回帰の算出¶
- 100点をカットオフスコアとし、以下の式を用いて平均への回帰を算出する。
- $E[y_1 | y_1 < k]$が1回目のテストが$k$点以下の受験者の1回目の平均値の期待値(1)
- $ E[y_2 | y_1 < k]$が1回目のテストが$k$点以下の受験者の２回目の平均値の期待値(2)
- (1)と(2)の差が差得点の期待値、すなわち、平均への回帰である。
- カットオフスコアを標準化する。
$$ z = \frac{x - \mu}{\sigma} $$
## １回目のスコアをx、2回目のスコアをyに保存する。
x = data["X"].values
y = data["Y"].values
mu = np.average(x)
sigma = np.std(x,ddof=1)
z = (100 - mu)/sigma
z
-0.26434293248823776
- $C$を求める。
$$ C = \frac{\phi (z)}{1- \Phi (z)} $$
ここで
$$\phi (z) = \frac{1}{\sqrt{2\pi}}exp^{-\frac{z^2}{2}}$$
$$\Phi (z) = \int_{\infty}^{z} \phi(x)dx $$
norm.cdf()は$1- \Phi(z)$の値。
phi = norm.pdf(z)
Phi = norm.cdf(z)
C = phi / Phi
C
0.9734348716713084
- 以下を求めて(2)の値から(1)の値を引く。
- 1回目の平均値の期待値$ E[y_1 | y_1 < k] = \mu - C_\sigma $ (1)
- 2回目の平均値の期待値$ E[y_2 | y_1 < k] = \mu + \rho C_\sigma $ (2)
rho = np.corrcoef(x,y)[0][1]
E_y2 = mu - rho * C * sigma
E_y1 = mu - C * sigma
E_y2 - E_y1
2.3503179132614207
- この引き算は以下と同じ。 $$C\sigma(1 - \rho)$$
C * sigma * (1 - rho)
2.350317913261422
- 実際の点数は約８点上昇しているがそのうち2.35点は平均への回帰と言える。
- ここで実際に$ E[y_1 | y_1 < k] = \mu - C_\sigma $となるか確認してみよう。
# 平均50、標準偏差10のデータを生成する。
x = np.random.normal(50, 10, 1000)
# カットオフスコアを30とし、30以下のデータの平均値を算出する
L = []
for i in x:
if i <= 30:
L.append(i)
np.average(L)
26.496855111254774
# 30を標準化する
z = (30 - np.average(x))/np.std(x,ddof=1)
z
-2.0371031574101486
- $ C = \frac{\phi (z)}{1- \Phi (z)} $を求める。
phi = norm.pdf(z)
Phi = norm.cdf(z)
C = phi /(Phi)
C
2.4061190974227036
- $ \mu - C_\sigma $を求める。
np.average(x) - C * np.std(x,ddof=1)
26.38062490557109
練習問題¶
"../DATA01/IEDA_reg2.csv"には500人の受験者が２回のテスト（1回目がX、2回目がY）を受験した結果が保存されている。このデータに関して、
- ヒストグラムを書きなさい。
- １回目のスコアをx軸、２回目のスコアをyとして散布図を書きなさい。また、これら2つのテストの回帰直線と１回目のテストと２回目のテストで同じスコアだった場合を示す直線($y = x$)を挿入しなさい。
- 1回目のテストのスコアと２回目のテストのスコアとの差を算出し、１回目のテストのスコアをx軸、差得点をy軸とした散布図を書き、これら２つの値の回帰直線を挿入しなさい。
- 3で求めた差得点と１回目のテストのスコアの相関係数を算出しなさい。
- 下位25%のカットスコアを算出し、この受験者群の１回目のテストのスコアの平均値の期待値(A)および２回目のテストのスコアの平均値の期待値(B)を求め、$B - A$（平均への回帰）を算出しなさい。