当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.29の「不完全データの統計処理」に関して演習問題を中心に解説を行います。欠損値の取り扱いは実際によく出てくる課題であるので、演習などを通して一通りの手法を確認しておくと良いと思われました。
本章のまとめ
欠測メカニズム
欠測が生じる原理は主にMCAR(Missing Completely At Random)、MAR(Missing Completely Random)、MNAR(Missing Not At Random)の3通りに分類される。$(x,y)$が観測されるとき、$y$の欠測はそれぞれ下記のように考えることができる。
・MCAR
-> yの欠測は完全にランダムに生じる。平均、標準偏差、相関係数、回帰の際の係数は欠測が生じなかった場合に概ね一致する。
・MAR
-> yの欠測はxの値に基づいて生じる。xとyに相関がある場合、回帰係数は概ね同じ一方で$y$の平均や標準偏差には変化が生じる場合がある。
・MNAR
-> yの欠測はyの値に基づいて生じる。xとyに相関がある場合、回帰係数やyの平均や標準偏差に変化が生じる場合がある。削除法・代入法
削除法(deletion method)は「欠測のある個体に対応する行を削除する手法」である。削除法はCC(Complete Case)解析とAC(Available Case)解析に分けられ、CC解析では欠測がある行を全て取り除く一方で、AC解析では使える値は全て用いて計算を行う。
欠測数があまり多くない場合や分散共分散行列を計算する場合はCC解析が妥当である一方で、変量数と欠測が多い場合はAC解析の方が妥当である。
補間法・代入法(imputation method)は「欠測箇所に何らかの値を代入する手法」である。主に下記の手法が用いられる。
i) 平均値代入
-> 着目する変数の観測値の平均値を代入する
ⅱ) 回帰代入
-> 回帰式によって欠測部分を予測する
ⅲ) Hot Dech法
-> 欠測のある個体と類似の個体を探し、その観測値を欠測部分に代入する演習問題解説
例29.1
下記を実行することで、打ち切りとトランケートされた場合の二つの場合の最尤推定値を計算できる。
import numpy as np
from scipy import stats
x_all = np.array([35., 38., 42., 45., 49., 52., 53., 56., 63., 65.])
x_observed = np.array([35., 38., 42., 45., 49., 52., 53., 56.])
c = 60.
sigma = 10.
mean_x_all = np.mean(x_all)
mean_x_observed = np.mean(x_observed)
mean_mle_censoring = np.mean(x_observed)
mean_mle_truncate = np.mean(x_observed)
for i in range(10):
a1 = (c-mean_mle_censoring)/sigma
mean_mle_censoring = (np.sum(x_observed)+(x_all.shape[0]-x_observed.shape[0])*(mean_mle_censoring+(stats.norm.pdf(a1)*sigma)/(1-stats.norm.cdf(a1))))/x_all.shape[0]
a2 = (c-mean_mle_truncate)/sigma
mean_mle_truncate = mean_x_observed + (stats.norm.pdf(a2)/stats.norm.cdf(a2))*sigma
print("Result_Censoring: {}".format(mean_mle_censoring))
print("Result_Truncate: {}".format(mean_mle_truncate))・実行結果
> print("Result_Censoring: {}".format(mean_mle_censoring))
Result_Censoring: 50.0523586591
> print("Result_Truncate: {}".format(mean_mle_truncate))
Result_Truncate: 48.6544391675どちらの結果からも標本平均が実際の平均よりも小さいことを改善していることが確認できる。
例29.2
下記を実行することで、$(29.5)$式〜$(29.8)$式の値の計算を行うことができる。
・$(29.5)$式
import numpy as np
x = np.array([35., 37., 40., 45., 47., 54., 57., 60., 62., 63])
print("mu_x: {:.1f}".format(np.mean(x)))
print("sigma_x^2: {:.1f}".format(np.sum((x-np.mean(x))**2)/x.shape[0]))・実行結果
mu_x: 50.0
sigma_x^2: 100.6・$(29.6)$式
import numpy as np
x = np.array([35., 37., 40., 45., 47., 54., 57.])
y = np.array([40., 38., 48., 41., 60., 46., 57.])
print("Mean x: {:.1f}".format(np.mean(x)))
print("s_x^2: {:.2f}".format(np.sum((x-np.mean(x))**2)/x.shape[0]))
print("Mean y: {:.2f}".format(np.mean(y)))
print("s_y^2: {:.2f}".format(np.sum((y-np.mean(y))**2)/y.shape[0]))
print("r: {:.2f}".format(np.sum((x-np.mean(x))*(y-np.mean(y)))/np.sqrt(np.sum((x-np.mean(x))**2)*np.sum((y-np.mean(y))**2))))・実行結果
Mean x: 45.0
s_x^2: 59.71
Mean y: 47.14
s_y^2: 62.41
r: 0.64・$(29.7)$式
import numpy as np
x = np.array([35., 37., 40., 45., 47., 54., 57.])
y = np.array([40., 38., 48., 41., 60., 46., 57.])
s_xy = np.sum((x-np.mean(x))*(y-np.mean(y)))/x.shape[0]
s_x2 = np.sum((x-np.mean(x))**2)/x.shape[0]
s_y2 = np.sum((y-np.mean(y))**2)/y.shape[0]
beta = s_xy/s_x2
alpha = np.mean(y)-beta*np.mean(x)
tau2 = s_y2 - s_xy**2/s_x2
print("alpha: {:.2f}".format(alpha))
print("beta: {:.2f}".format(beta))
print("tau^2: {:.2f}".format(tau2))・実行結果
alpha: 17.65
beta: 0.66
tau^2: 36.75・$(29.8)$式
import numpy as np
x_all = np.array([35., 37., 40., 45., 47., 54., 57., 60., 62., 63])
x = np.array([35., 37., 40., 45., 47., 54., 57.])
y = np.array([40., 38., 48., 41., 60., 46., 57.])
sigma_x2 = np.sum((x_all-np.mean(x_all))**2)/x_all.shape[0]
s_xy = np.sum((x-np.mean(x))*(y-np.mean(y)))/x.shape[0]
s_x2 = np.sum((x-np.mean(x))**2)/x.shape[0]
s_y2 = np.sum((y-np.mean(y))**2)/y.shape[0]
beta = s_xy/s_x2
alpha = np.mean(y)-beta*np.mean(x)
tau2 = s_y2 - s_xy**2/s_x2
print("mu_y: {:.2f}".format(alpha+beta*np.mean(x_all)))
print("sigma_y: {:.2f}".format(tau2+beta**2*sigma_x2))
print("sigma_xy: {:.2f}".format(beta*sigma_x2))
print("rho: {:.2f}".format(beta*sigma_x2/np.sqrt(sigma_x2*(tau2+beta**2*sigma_x2))))・実行結果
mu_y: 50.42
sigma_y: 79.98
sigma_xy: 65.94
rho: 0.74$(29.8)$式の結果の計算にあたって有効数字の取り扱いを行わなかったことで、ワークブックの解答と少々異なる結果が得られた。
問29.1
$[1]$
⑤
$2$回の試験の間の相関係数が正であり、$1$回目の試験の合否で$2$回目の試験を受験するかが決まることから、$2$回目の試験結果が良いと思われる生徒が$2$回目の試験を受験していないと考えることができる。
よって、$\bar{y}_B$は$\mu_Y$を過小評価し、同様に$r_{B}$は$\rho$を過小評価すると考えられる。
$[2]$
②
標本平均には偏りがないが、回帰による予測値が回帰直線に一致することから相関係数を過大に評価すると考えることができる。
参考
・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1