コルモゴロフ・スミルノフ統計量によるランダマイゼーション検定
The Kolmogorov-Smirnov Two-Sample Test by Randomization
２つのデータセットの分布における何らかの特性、平均とか分散などの違いを、ランダマイゼーションによるコルモゴロフ・スミルノフ２サンプル検定によって調べるPythonプログラムを作成してみた。ランダマイゼーション検定は、母集団の確率分布に関わらず行うことができる。ランダマイゼーション検定の具体的方法は、Imbens & Rubin (2015; pp. 76-78)、岡本（2022; pp. 171-174）などを参照されたい。
２つの独立なデータセットが図１のように与えられたとする。図１のデータは、Gibbons
& Chakraborti (2003) のExample
2.1 (p. 237) におけるものをExcelデータとしたものである。
図１ Gibbons
& Chakraborti (2003) Example 2.1 (p. 237)より。
第１列は、データのunit番号である。この列の値は、分析では用いられないので、任意の値でよい。第２列はデータ収集条件であり、実験群W=1、統制群W=0である。第３列がデータ値であり、実験群W=1のデータは、標準正規分布からのサンプルである。統制群W=0は、自由度18のカイ2乗分布からのサンプルを標準化したものである。
データに対して経験分布（empirical distribution function）を次式で与える（Imbens & Rubin, 2015, p. 69）。
このとき、コルモゴルフ・スミルノフ検定統計量は次式で与えられる（Imbens & Rubin, 2015, Eq. (5.8)）
式（5.8）によるランダマイゼーション検定を行うプログラムをリスト１のように用意した。プログラムファイルなどは圧縮ファイルKS2smplsTestFiles.zipにまとめたので、ダウンロード展開すれば自由に利用できる。
コルモゴルフ・スミルノフ統計量による検定は、Imbens & Rubin (2015)に掲載されている結果（Table 5.9）では、共変量covariateを用いない検定の中で最も検定力の高いものであったことが示されている。
リスト１のプログラムを実行すると、入力データファイル名を聞いてくる。入力データファイルは、図１の形式のExcelファイルを用意する。次の例では、図１のファイルが設定されている。
Input data
file (*xlsx) = Example2_1.xlsx
ファイル名を設定して、Enterキーを押すと、データが読み込まれて計算が始まる。まず、実験群（treated）、統制群（control）の経験分布のグラフが表示される（図２）。
図２
経験分布は、次のコードのクラス型CumDistにより計算している。
class CumDist:
"""
Cumulative distribution
"""
def
__init__(self, y):
self.N = len(y)
self.y = sorted(y)
def
cdf(self, v):
pos = 0
while True:
if v < self.y[pos]:
break
if pos < self.N:
pos += 1
if pos >= self.N:
break
return pos / self.N
式（5.8）の統計量は、次のコードで求めている。
cum_yt = CumDist(Yt) # Cumulative distribution for Yt
cum_yc = CumDist(Yc) # Cumulative distributin
for Yc
CumVs_t = np.array([cum_yt.cdf(v) for v in Y])
CumVs_c = np.array([cum_yc.cdf(v) for v in Y])
Tks = np.max(np.abs(CumVs_t - CumVs_c))
図２のフォームの右上角のX印をクリックしてフォームを閉じると、ランダマイゼーションによる計算が始まる。
ランダマイゼーションは、Wのシャフリングにより、次のコードで行っている。
npg_shuffle(W)
Yt
= Y[W==1]
Yc
= Y[W==0]
シャフリングに対応する式（5.8）の統計量はdiff_maxとして、次のコードで求めている。
cum_yt
= CumDist(Yt)
cum_yc
= CumDist(Yc)
CumVs_t
= np.array([cum_yt.cdf(v) for
v in Y])
CumVs_c
= np.array([cum_yc.cdf(v)
for v in Y])
diff_max
= np.max(np.abs(CumVs_t - CumVs_c))
算出されたdiff_maxが、統計量Tks以上である回数を数えるための変数cntを以下のように用意している。
if diff_max
>= Tks:
cnt += 1
ランダマイゼーション・シミュレーションの終了後、diff_maxがTks以上であった比率をｐ値の推定値としている。
p_value = cnt / K
プログラムの実行終了後、端末には以下の結果が出力されている。
Tks = 0.250 p-value = 0.9801 s.e. =
0.00014
上記のp-value = 0.9801という結果は、Gibbons & Chakraborti (2003), p. 246に記載されているExact p-value = 0.9801と同じ値である。ただし、ランダマイゼーションによる結果は、シミュレーションごとに変動するので、常に一致するわけではないことに注意する必要がある。変動の大きさは、s.e. = 0.00014により推測できる。
Gibbons, J. D. & Chakraborti, S. (2003). Nonparametric statistical inference, fourth edition, revised and expanded. CRC Press.
Imbens, G. W. and Rubin, D. B. (2015). Causal
inference for statistics, social, and biomedical sciences: An introduction.
Cambridge University Press.
岡本安晴（2019）いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan, PyMC）―― 丸善出版
岡本安晴（2022）データを使いこなすための統計入門 22世紀アート
import pandas as pd
import numpy
as np
import matplotlib.pyplot
as plt
#
# Exc el形式データファイルの読み込み
#
infl_nm = input('Input data file (*xlsx) = ')
with pd.ExcelFile(infl_nm) as f_xlsx:
data_xlsx
= pd.read_excel(f_xlsx,
header=None)
# データ値をndarray型の配列に変換
data = data_xlsx.values
print('\n', data)
W = data[1:,1]
Y = data[1:,2]
print('W =\n',W)
print('Y =\n',Y)
Yt = Y[W==1]
# Treated group
Yc = Y[W==0]
# Control group
print('Yt
=\n', Yt)
print('Yc
=\n', Yc)
class CumDist:
"""
Cumulative distribution
"""
def
__init__(self, y):
self.N = len(y)
self.y = sorted(y)
def
cdf(self, v):
pos = 0
while True:
if v < self.y[pos]:
break
if pos < self.N:
pos += 1
if pos >= self.N:
break
return pos / self.N
min_y = np.min(Y)
max_y = np.max(Y)
cum_yt = CumDist(Yt) # Cumulative distribution for Yt
x_coord = np.arange(min_y-1,
max_y+1, 0.01)
yt_coord = [cum_yt.cdf(v)
for v in x_coord]
plt.plot(x_coord, yt_coord, label='treated')
cum_yc = CumDist(Yc) # Cumulative distributin
for Yc
yc_coord = [cum_yc.cdf(v)
for v in x_coord]
plt.plot(x_coord, yc_coord, label='control')
plt.legend()
plt.show()
CumVs_t = np.array([cum_yt.cdf(v) for v in Y])
CumVs_c = np.array([cum_yc.cdf(v) for v in Y])
Tks = np.max(np.abs(CumVs_t - CumVs_c))
npg_shuffle = np.random.default_rng().shuffle # Shuffling用関数
K = 1_000_000
cnt = 1
for i in
range(K-1):
npg_shuffle(W)
Yt
= Y[W==1]
Yc
= Y[W==0]
cum_yt
= CumDist(Yt)
cum_yc
= CumDist(Yc)
CumVs_t
= np.array([cum_yt.cdf(v) for
v in Y])
CumVs_c
= np.array([cum_yc.cdf(v)
for v in Y])
diff_max
= np.max(np.abs(CumVs_t - CumVs_c))
if diff_max
>= Tks:
cnt += 1
if i
% 10000 == 0 and i > 0:
print(f'{(i+1):,}/{K:,} {cnt/(i+1):.5f}, {diff_max:5f}')
p_value = cnt / K
s_e = (p_value
* (1.0 - p_value) / K) ** 0.5
print(f'\nTks
= {Tks:.3f} p-value =
{p_value:.4f} s.e. = {s_e:.5f}')