リサンプリングを用いた多重比較

理論的な説明はとりあえずすっ飛ばして、
今回は、果たしてどういう計算をしているのか、という視点でまとめます。

今回使うのは ^[1] に書かれている有名な方法で、一部では min P とも呼ばれるらしいです。
リサンプリングによって全ての比較が帰無仮説の時 $H_{0}^{C}$ の帰無分布を生成し、帰無分布から各比較の $P$ 値の最小値の分布を求め、標本における各比較の $P$ 値がどの辺りに位置するのかを計算すると、調整 $P$ 値
\[
\tilde{p}_i=\Pr\left(\min_{1\leq j\leq k} P_j \leq p_i \mathrel{}\middle|\mathrel{} H_{0}^{C}\right)
\] が得られます。
$\tilde{p_i}$ が $i$ 番目の比較の調整 $P$ 値、$k$ はリサンプリング回数です。
FWER を押さえるとかそのあたりの話は、本を参照。

リサンプリングをして $k=20000$ の帰無分布を作ります。
※以下の例では、SAS で実装されている Dunnett 型の母平均の多重比較の結果との一致を確認するため、multtest procedure の結果を出しておきます。

data d;
call streaminit(642954121);
do g = 1 to 3;
  do r = 1 to 5;
  y = rand('Normal') + (g/2)**2;
  output;
end; end;
proc multtest data = d boot n = 20000
  outsamp = y seed = 74539234;
  class g;
  test mean(y); 
  contrast 'dunnett-1' 0.5 -0.5 0;
  contrast 'dunnett-2' 0.5 0 -0.5;
run;

                     p-Values
Variable    Contrast            Raw     Bootstrap
y           dunnett-1        0.0041        0.0111
y           dunnett-2        0.0009        0.0032

帰無分布から各比較の $P$ 値の最小値の分布を求め、標本の $P$ 値と比較します。

/* 標本のP値を得るためリサンプリング分布にくっつけます */
data y;
  set d(in=x) y;
  if x then do; _sample_ = -1;
    if g = 1 then _class_ = '1';
    if g = 2 then _class_ = '2';
    if g = 3 then _class_ = '3';
  end;
/* 統計量とP値の計算をして */
proc summary data = y;
  var y;
  by _sample_ _class_;
  output out = stat(keep=_class_ n m s _sample_) n = n mean = m css = s;
data stat;
  merge
    stat(where=(_class_="1") rename=(n=n1 m=m1 s=s1))
    stat(where=(_class_="2") rename=(n=n2 m=m2 s=s2))
    stat(where=(_class_="3") rename=(n=n3 m=m3 s=s3));
  retain p1obs p2obs minp 0;
  df = n1 + n2 + n3 - 3;
  s = (s1 + s2 + s3) / df;
  t1 = (m2 - m1) / sqrt(s * (1 / n1 + 1 / n2));
  t2 = (m3 - m1) / sqrt(s * (1 / n1 + 1 / n3));
  p1 = 2 * (1 - cdf('t', abs(t1), df));
  p2 = 2 * (1 - cdf('t', abs(t2), df));

  /* 標本のP値が、リサンプリング分布のmin Pを上回った回数をカウント */
  minp = min(p1, p2);
  if _sample_ = -1 then do; p1obs = p1; p2obs = p2; end;
  if p1obs >= minp then pt1 = 1; else pt1 = 0;
  if p2obs >= minp then pt2 = 1; else pt2 = 0;
  drop _class_;
proc means data = stat mean;
  where _sample_ ^= -1;
  var pt1 pt2;
run;

めでたし、めでたし。

   The MEANS Procedure
Variable            Mean
------------------------
pt1            0.0110500
pt2            0.0031500
------------------------

他の対比係数行列の場合も調べてみよう。
R だと $k$ 回 sample して統計量を求める関数に突っ込んで、pmin するだけ？重そう。

帰無分布を出すときに標準化したり、不等分散の場合だったり、不等標本数の場合だったりこの辺の細かいところをもう少し詰めたいな。