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ステートメントを使用して統計処理を行い、最後に得られるデータセットからシミュレーションの結果を導く、ということはよく行います。
以降の回では、このプログラムを改善して、より効率的なものを作成していきます。