2段階t検定をSASでシミュレーション(1)

TeX方面で著名な奥村さんが、自身のブログでこんな記事をアップしていました。

2段階t検定の是非
http://oku.edu.mie-u.ac.jp/~okumura/blog/node/2262/
2群の平均値の差の検定で，まずF検定で分散が等しいかどうか検定してから，通常のt検定かWelchの検定かに振り分けることを勧める本やサイトがまだ多い。

なるほど、これは昔からちょっと変、と思っていました。明らかに等分散性が成立すると考えられる時以外は、Welch検定を行うべきでしょう。成立すると言えそうなのは、物理や化学の実験データくらいしかなさそうですが。

ということで、早速SASでやってみましょう。同じように100万回のシミュレーションを行って、SASの速さを見せつける、、、といきたいところですが、諸般の事情によりとりあえず1万回です。

%let loop=10000;
%let n1=10;
%let mean1=0;
%let sd1=1.5;
%let n2=30;
%let mean2=0;
%let sd2=1;

data test;
  call streaminit(12345);
  array n{2} (&n1 &n2);
  array mean{2} (&mean1 &mean2);
  array sd{2} (&sd1 &sd2);
  do loop=1 to &loop;
    do group=1 to 2;
	  do i=1 to n[group];
	    y=rand('normal',mean[group],sd[group]);
		output;
	  end;
	end;
  end;
run;

ods listing close;
proc ttest data=test;
  by loop;
  class group;
  var y;
  ods output ttests=t equality=e;
run;
ods listing;

data ttestout;
  merge t e;
  by loop;
  if first.loop then do;
    output;
	if probf>0.05 then do;
      variances="2step05";
	  output;
	end;
	if probf>0.20 then do;
	  variances="2step20";
	  output;
	end;
  end;
  else do;
	if probf<0.05 then do;
      variances="2step05";
	  output;
	end;
	if probf<0.20 then do;
	  variances="2step20";
	  output;
	end;
	variances="Unequal";
    output;
  end;
run;
proc format;
  value pfmt low-0.05="1, under 0.05"
             0.05-high="2, over 0.05";
  value $varfmt "Equal"="1, equal variance"
               "2step05"="2, 2step, 0.05"
               "2step20"="3, 2step, 0.20"
               "Unequal"="4, unequal variance";
run;

proc freq data=ttestout order=formatted;
  tables variances*probt/nocol nopercent;
  format probt pfmt. variances $varfmt.;
run;
/*
Variances         Probt(Pr > |t|)

度数             |
行のパーセント   |1, under|2, over |   合計
                 | 0.05   |0.05    |
-----------------+--------+--------+
1, equal varianc |   1080 |   8920 |  10000
e                |  10.80 |  89.20 |
-----------------+--------+--------+
2, 2step, 0.05   |    814 |   9186 |  10000
                 |   8.14 |  91.86 |
-----------------+--------+--------+
3, 2step, 0.20   |    662 |   9338 |  10000
                 |   6.62 |  93.38 |
-----------------+--------+--------+
4, unequal varia |    555 |   9445 |  10000
nce              |   5.55 |  94.45 |
-----------------+--------+--------+
*/

当然ながら、ほぼ同様の結果が得られています。
このように、DATAステップでサンプルデータを乱数から作成し、次に分析プロシジャでBYステートメントを使用して統計処理を行い、最後に得られるデータセットからシミュレーションの結果を導く、ということはよく行います。
以降の回では、このプログラムを改善して、より効率的なものを作成していきます。