ランク統計量によるランダマイゼーション検定
平均値の比較をランクにより行う。Imbens & Rubin (2015)は、ランク統計量による検定は“attractive”であると書いている（pp. 74, 81）。ランダマイゼーション検定は、母集団の確率分布に関わらず行うことができるが、ランダマイゼーション検定の具体的方法は、Imbens & Rubin (2015; pp. 76-78)、岡本（2022; pp. 171-174）などを参照されたい。
２つの独立なデータが図１のように与えられたとする。図１のデータは、Imbens & Rubin (2015)のTable 5.3から抜き出したものである。
図１ Imbens & Rubin (2015), Table 5.3より。
第１列は、データのunit番号である。この列の値は、分析では用いられないので、任意の値でよい。第２列はデータ収集条件であり、実験群W=1、統制群W=0である。第３列がデータ値Yである。
データに対して、ランキングの値の平均が０になるように
が引かれた値を用いる。ここで、
である。
このとき、ランク統計量は次式で与えられる（Imbens & Rubin, 2015; Eq. (5.7)）。
式（5.7）の
によりランダマイゼーション検定を行うプログラムをリスト１のように作成した。プログラムのファイルなどはまとめて圧縮ファイルranktestfiles.zipとした。ダウンロード展開すれば自由に用いることができる。
リスト１のプログラムを実行すると、入力データファイル名を聞いてくる。
Input data
file (*.xlsx) =
図１の形式で用意したファイル名を設定する。図１のファイルを用いる場合は、以下のようになる。
Input data
file (*.xlsx) = HoneyStudyTable5_3.xlsx
入力データファイル名を設定すると、計算が始まる。
計算が終了すると、計算結果が以下のように出力される。
T_rank =
0.667, p-value = 0.800, s.e. = 0.00040
ランク統計量、ｐ値、ランダマイゼーションによる推定値のs.e.が表示されている。
直接計算されたｐ値は0.8であるが（Imbens & Rubin, 2015; p. 68）、これと比べてランダマイゼーションというシミュレーションで推定された値は満足できるものであると言えよう。
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
infl_nm = input('Input data file
(*.xlsx) = ')
#'HoneyStudyTable5_3.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]
# Treatments
Y = data[1:,2]
# Observed Outcome
print('W =\n',W)
print('Y =\n',Y)
#
#
Yの値により並べ換えてランク付け
#
s_idx = np.argsort(Y)
print('s_idx =\n', s_idx)
s_W = W[s_idx]
print('s_W =\n', s_W)
s_Y = Y[s_idx]
print('s_Y =\n', s_Y)
def ranking(x):
"""
ランク付け用関数
"""
r = np.empty(len(x))
s_p = 0
c_p = s_p
while True:
#print('s_p =', s_p)
while True:
if c_p < len(x) - 1:
if x[s_p] == x[c_p+1]:
c_p += 1
else:
break
else:
break
for i in
range(s_p, c_p+1):
r[i] = (s_p + c_p) / 2 + 1
s_p
= c_p + 1
c_p
= s_p
if
s_p >= len(x):
break
return r
Rk = ranking(s_Y)
print('Rk =\n', Rk)
Nt = np.sum(W)
# 実験群の数
Nc = len(W) - Nt
# コントロール群の数
print('Nt =', Nt)
print('Nc =', Nc)
Rk = Rk - (Nt + Nc + 1) / 2 # Normalized rank
print('sum_Rk =', sum(Rk))
sum_Rk_t = np.dot(Rk, s_W)
# 実験群のランク総和
sum_Rk_c = np.dot(Rk, 1 - s_W) # コントロール群のランク総和
T_rank = np.abs(sum_Rk_t/Nt -
sum_Rk_c/Nc) # データに対する統計量T
print('T_rank =', T_rank)
#
# ランダマイゼーション検定のためのシミュレーション
#
npg_shuffle =
np.random.default_rng().shuffle # Shuffling用関数
K = 1_000_000
cnt_p = 1
# データに対する統計量T以上の場合の数
for t in range(K-1):
npg_shuffle(s_W)
sum_Rk_t = np.dot(Rk, s_W)
sum_Rk_c = np.dot(Rk, 1 -
s_W)
vT_rank = np.abs(sum_Rk_t/Nt
- sum_Rk_c/Nc)
if vT_rank >= T_rank:
cnt_p += 1
if t % 10_000 == 0 and t
> 0:
print(f't
= {t}, p-value(temp.) =
{cnt_p/t:.5f}')
p_value = cnt_p / K
s_e_ = (p_value * (1 - p_value) / K)
** 0.5
print(f'\nT_rank = {T_rank:.3f}, p-value = {p_value:.3f}, s.e. = {s_e_:.5f}')