ここでは、Kirk (1982, p.354) の Table 8.3.1 の完全無作為化２要因釣り合い型デ ザインデータ、及び竹内（監）・高橋ら (1990, p.297) の完全無作為化２要因非釣り 合い型デザインデータを例に取り、SAS による分析の方法を以下の４つのサブセクショ ンに分けて示す：
| 1.4.6.1 釣り合い型デザインプログラムとその 概要 |
| 1.4.6.2 釣り合い型デザインプログラムの出力結 果の概要 |
| 1.4.6.3 非釣り合い型デザインプログラムとそ の概要 |
| 1.4.6.4 非釣り合い型デザインプログラムの出力 結果の概要 |
また、この節には各プログラムの直後に以下のような SAS プログラムのダウンロー ドコーナーを用意してある：
| 1.完全無作為化２要因釣り合い型デザインのプログラム例 |
| 2.完全無作為化２要因非釣り合い型デザインのプログラム例 |
ここでは、２つのデータを例に取り、SAS による完全無作為化２要因デザインの 分析プログラムの作成方法について述べる。まず最初 は、Kirk (1982, p.354) の Table 8.3-1 のデータについて述べる。このデータは、 釣り合い型デザインデータである。もう１つは、竹内（監）・高橋ら (1990, p.297) の Table 15.7 の非釣り合い型でなおかつ因子相互が直交していないデザインで ある。２つのデータのうち、後者については高橋らも述べているように、1.3.7 節で 述べた SAS の４つの平方和がモデルにより微妙に異なるので、注意する必要がある。 彼らは、このような場合でかつサンプル数が極端にアンバランスな場合には、Type II 平方和を選択することを薦めている。また、後者データの場合、anova プロシジャ は用いない方がよい。これを用いると、計算結果は出力するが、ログ情報として SAS は glm を薦める。
1.3.6 節の場合と同様、それぞれのプログラムには、最初に長いコメントが付け てある。さらに、最後に分散分析の前提についての幾つかの検討のためのプログラムも 載せてある。その結果、長いプログラムになっているので、幾つかに分けて示す。 ユーザーは、これらを利用するに際しては、それらの中から必要なものを選んで 自分用に書き換えればよい。
それぞれのプログラムの利用に際しての最後の注意点は、つぎのようである。 すなわち、最初のデータの場合、全体的交互作用が有意になるので、anova もしくは glm プロシジャの means 文による対比検定は行わない方がよかろう。なぜならば、 そのような場合、各主効果の検定自身が 1.4.4 節で述べたように消極的な意味しか 持たせられない場合が多いであろうから。
一方、後者のデータの場合には、全体的交互作用は有意ではないので、部分的 交互作用の検定も行っていない点に、注意せよ。
プログラムは合計１０個の部分から成り立っており、以下に示すように順にそれぞれ について簡単な解説を施してある：
*-------------------------------------------------------------------------* | December 8, 1995 | | sas program -- crf33kk.sas | | | | Example of a balanced CRF-33 design analysis for Kirk (1982, p.354) | | Table 8.3-1 data. This data has 3 levels for both factor A and B. | | | | Cautions: | | | | (1) According to Winer (1971), Bartlett's test for homogeneity of vari- | | ance (Bartlett, 1937) is the most widely used test. However, the rou- | | tine use of it as a prelliminary test is not recommended. For, it is | | extremely sensitive to nonnormality (Scheffe, 1959, p.83). | | (2) Test for homogeneity of variance may be unnecessary for balanced | | designs like this data set. | | (3) Proc anova must be used only for balanced designs. | | (4) For pairwise comparisons, the tukey procedure is superior to the | | Scheffe procedure. | | (5) Use cldiff option in the mean statement for all pairwise compari- | | sons. | | (6) CLDIFF is the default for unbalanced designs unless DUNCAN, REGWF, | | REGWQ, SNK, or WALLer is specified. | | (7) Lsmeans option is not available with proc anova. | | (8) Contrast statement is not available with proc anova. | | (9) Lsmeans can be computed for any effect involving class variables | | as long as the effect is defined in the model. In other words, no | | effects should be specified in the Lsmeans statement unless they are | | defined in the model statement. | | (a) There is no limit to the number of CONTRAST statements, but thy must| | appear after the MODEL statement. | | (b) If overall interaction is statistically significant, it may neither | | be recommended to perform main effect tests nor to perform contrast | | tests of these effects. | | (c) Values of contrast statement, i.e., contrast 'label' effect values, | | for contrast-contrast interactions can be easily written by computing | | the Kronecker product of the contrast matrices corresponding to the | | two factors under study. | | | *-------------------------------------------------------------------------*; |
完全無作為化２要因デザインのプログラム crf33kk.sas についてのコメントから 成る。
options ps=60 ls=80;
/* (1) data input */
data work;
do a=1 to 3;
do b=1 to 3;
input y @;
output;
end;
end;
label a='treatment A'
b='treatment B'
y='dependent variable';
cards;
24 44 38 30 35 26 21 41 42
33 36 29 21 40 27 18 39 52
37 25 28 39 27 36 10 50 53
29 27 47 26 31 46 31 36 49
42 43 48 34 22 45 20 34 64
;
|
cards 文の下にある Kirk (1982, p.354) の Table 8.3-1 のデータは、２要因完全 無作為化要因デザインデータで、各要因とも３水準から成る。各水準は、５名の被験者 が無作為に割り付けられている。各行ともその９列は順に (A_1, B_1)、(A_1, B_2)、 (A_1, B_3)、(A_2, B_1)、(A_2, B_2)、(A_2, B_3)、(A_3, B_1)、(A_3, B_2)、 (A_3, B_3)、の各セルが対応する。
このように並べられたデータに対して、各セルの要因の組み合わせのコードを変数 a 及び変数 b の値として入力しつつ、それに対応する従属変数の値 y をも同時に読み 込むために、ここでは上のプログラムで do 文を用いた特殊な input 文によるデータ 入力形式を取っている。
/* (2-1) balanced CRF-33 ANOVA using proc anova */ title 'CRF-33 ANOVA for the Kirk (1982, p.354) data'; proc anova data=work; class a b; model y=a b a*b; /* model y=a|b; */ means a b/tukey alpha=0.01; /* means a b/tukey cldiff alpha=0.01; */ run; |
このデータは、釣り合い型デザイン、すなわち各セルの被験者数は等しいので、 SAS では上の例のように anova プロシジャを用いた分散分析が可能である。もちろん、 ２要因のそれぞれの主効果が有意の時、要因毎 tukey 法による対比較が可能なように、 means 文で指定がしてある。
/* (2-2) balanced CRF-33 ANOVA using proc glm */ options ps=60 ls=80; title 'CRF-33 ANOVA by proc glm'; proc glm data=work; class a b; model y=a|b; contrast 'c1(a)' a 1 -1 0; contrast 'c2(a)' a 0.5 0.5 -1; contrast 'c1(b)' b 1 -1 0; contrast 'c2(b)' b 1 0 -1; contrast 'c3(b)' b 0 1 -1; means a|b; means a b/tukey alpha=0.01; /* means a b/tukey cldiff alpha=0.01; */ lsmeans a|b/ out=glmout; run; |
ここでは、同一データを anova プロシジャでなく glm プロシジャで行うプログラム である。サンプル数が等しい場合には、どちらでも分析が可能であるし、結果も同一と なる。ただし、glm プロシジャでは anova ではできない非対比較を、上のプログラムの ように contrast 文を用いることにより検討できる。また、各セルのサンプル数が異 なる場合用の lsmeans 文によるセルの平均値の計算も可能となる。このデータのよう に各セルともサンプル数が等しい場合には、means 文による平均値と同一の結果と なる。
/* (3) means plot using proc glm/lsmeans output option */ options ps=30 ls=80; data work2; set glmout; if a=. or b=. then delete; run; title 'means plot utilizing proc glm/lsmeans'; proc plot data=work2; plot lsmean*a=b; run; |
このプログラムは、data ステップと plot プロシジャを用いて (2-2) の glm プロ シジャの lsmeans 文の out オプションにより書き出された一時ファイル glmout を呼 び出し、lsmeans 文により計算された各セルの平均値をプロットするためのものであ る。
/* (4) test for contrast-contrast interactions
using proc glm/contrast statement */
options ps=60 ls=80;
title 'test for contrast-contrast interactions';
proc glm data=work outstat=glmstat;
class a b;
model y=a|b;
contrast 'c1(a)-c1(b)' a*b 1 -1 0 -1 1 0 0 0 0;
contrast 'c1(a)-c2(b)' a*b 1 0 -1 -1 0 1 0 0 0;
contrast 'c1(a)-c3(b)' a*b 0 1 -1 0 -1 1 0 0 0;
contrast 'c2(a)-c1(b)' a*b 0.5 -0.5 0 0.5 -0.5 0 -1 1 0;
contrast 'c2(a)-c2(b)' a*b 0.5 0 -0.5 0.5 0 -0.5 -1 0 1;
contrast 'c2(a)-c3(b)' a*b 0 0.5 -0.5 0 0.5 -0.5 0 -1 1;
output out=temp residual=resid;
run;
|
このプログラムは、当該データの主効果と（全体的）交互作用の計算までは (2-2) の項の glm プロシジャによるプログラムと同一であるが、（全体的）交互作用が有意 の場合に備えて部分交互作用の１種である対比・対比交互作用の幾つかを計算するた めのものである。ここでの６つの対比・対比交互作用は、(2-2) の項での要因 a およ び要因 b のそれぞれ２通り及び３通りの対比の組み合わせによるもので、つぎの節 で述べるクロネッカー積を用いて上の contrast 文のような係数行列を計算する。
/* computation and print of the Scheffe adjusted F tests for
contrast-contrast interactions */
data work3;
set glmstat;
if _type_ ne 'CONTRAST' then delete;
ic=3; /* ic=number of levels of factor a */
jc=3; /* jc=number of levels of factor b */
n=45; /* n=total number of samples */
ndf=(ic-1)*(jc-1);
ddf=n-ic*jc;
fschef=f/ndf;
pschef=1-probf(fschef,ndf,ddf);
drop ic jc n ndf ddf;
run;
title 'sas original and Scheffe adjusted F tests';
proc print data=work3;
run;
|
(4) の glm での contrast 文により計算される対比・対比交互作用の検定では、 複数のそれらの危険率のコントロールが為されていない。このコントロールをシェ ッフェ方式でおこなうのが上のプログラムである。最後の print プロシジャは、単に これを出力するためのも野である。
/* (5-1) posterior analysis (1): histogram plot of error estimates */ options ps=30 ls=80; title 'histogram of the error estimates '; proc chart data=temp; vbar resid; run; |
上のプログラムは、(4) の項の glm プロシジャの最後に output 文を用いて一時 ファイルに出力したモデルの残差推定値を呼び出し、chart プロシジャによりその ヒストグラムを描かせるためのものである。
/* (5-2) posterior analysis (2): test for normality of error estimates */ options ps=60 ls=80; title 'test for global normality of the dependent variable'; proc univariate data=temp normal; var resid; run; |
これは、univariate プロシジャによるモデルの誤差推定値の分布の正規性の検定の ためのプログラムである。
/* (5-3) posterior analysis (3) : test for homogeneity of variances */
data work4;
set work;
ic=3; /* ic=number of levels of factor a */
jc=3; /* jc=number of levels of factor b */
do i=1 to ic;
do j=1 to jc;
if a=i and b=j then group=jc*(i-1)+j;
end;
end;
run;
title 'Bartlett test for homogeneity of variances';
proc discrim data=work4 pool=test slpool=0.05;
class group;
var y;
run;
|
ここでは、まずデータステップを使って、当該データの２要因の水準の組み合わせ による９個のセルに該当する被験者に group なる新たな変数とその値を定義しておき、 discrim プロシジャによるセルの分散の等質性の検定を行うためのプログラムを示し た。
| crf33kk.sas |
(2-1) のプログラムによる最初の主要な結果として、つぎのような分散分析表 が出力される：
CRF-33 ANOVA for the Kirk (1982, p.354) data
Analysis of Variance Procedure
Dependent Variable: Y dependent variable
Sum of Mean
Source DF Squares Square F Value Pr > F
Model 8 2970.0000000 371.2500000 5.94 0.0001
Error 36 2250.0000000 62.5000000
Corrected Total 44 5220.0000000
R-Square C.V. Root MSE Y Mean
0.568966 22.58770 7.9056942 35.000000
Source DF Anova SS Mean Square F Value Pr > F
A 2 190.0000000 95.0000000 1.52 0.2324
B 2 1543.3333333 771.6666667 12.35 0.0001
A*B 4 1236.6666667 309.1666667 4.95 0.0028
|
前節の１要因の場合とは対照的に、２種類の Source（変動因）の内容は異なる。前者 の変動因は Model として一括されているのに対して、後者のそれは２つの要因の主効果 と（全体的）交互作用の合計３つから成る。Pr > F の値から、要因 A は有意ではない が、要因 B は１パーセント以上の高い水準で統計的に有意であることがわかる。また、 交互作用も１パーセント以上の高い水準で統計的に有意である。
つぎに、(2-1) のプログラムの中の means 文による主効果の対比検定による結果と して、つぎのような tukey 法による要因 A 及び B についての対比検定結果が出力さ れる：
CRF-33 ANOVA for the Kirk (1982, p.354) data
Analysis of Variance Procedure
Tukey's Studentized Range (HSD) Test for variable: Y
NOTE: This test controls the type I experimentwise error rate, but
generally has a higher type II error rate than REGWQ.
Alpha= 0.01 df= 36 MSE= 62.5
Critical Value of Studentized Range= 4.396
Minimum Significant Difference= 8.9729
Means with the same letter are not significantly different.
Tukey Grouping Mean N A
A 37.333 15 3
A
A 35.333 15 1
A
A 32.333 15 2
Analysis of Variance Procedure
Tukey's Studentized Range (HSD) Test for variable: Y
NOTE: This test controls the type I experimentwise error rate, but
generally has a higher type II error rate than REGWQ.
Alpha= 0.01 df= 36 MSE= 62.5
Critical Value of Studentized Range= 4.396
Minimum Significant Difference= 8.9729
Means with the same letter are not significantly different.
Tukey Grouping Mean N B
A 42.000 15 3
A
B A 35.333 15 2
B
B 27.667 15 1
|
つぎに、(2-2) の項の glm を用いての同一データの分析結果が出力される。最初は、 (2-1) の項の anova プロシジャによるそれらと同様、要因の水準情報、サンプル数、 分散分析表が出力される。ここでは、分散分析表のみを示す：
CRF-33 ANOVA by proc glm
General Linear Models Procedure
Dependent Variable: Y dependent variable
Sum of Mean
Source DF Squares Square F Value Pr > F
Model 8 2970.0000000 371.2500000 5.94 0.0001
Error 36 2250.0000000 62.5000000
Corrected Total 44 5220.0000000
R-Square C.V. Root MSE Y Mean
0.568966 22.58770 7.9056942 35.000000
Source DF Type I SS Mean Square F Value Pr > F
A 2 190.0000000 95.0000000 1.52 0.2324
B 2 1543.3333333 771.6666667 12.35 0.0001
A*B 4 1236.6666667 309.1666667 4.95 0.0028
Source DF Type III SS Mean Square F Value Pr > F
A 2 190.0000000 95.0000000 1.52 0.2324
B 2 1543.3333333 771.6666667 12.35 0.0001
A*B 4 1236.6666667 309.1666667 4.95 0.0028
Contrast DF Contrast SS Mean Square F Value Pr > F
c1(a) 1 67.5000000 67.5000000 1.08 0.3056
c2(a) 1 122.5000000 122.5000000 1.96 0.1701
c1(b) 1 440.8333333 440.8333333 7.05 0.0117
c2(b) 1 1540.8333333 1540.8333333 24.65 0.0001
c3(b) 1 333.3333333 333.3333333 5.33 0.0268
|
うえの結果は、(2-1) の項の同一データに対する anova プロシジャによる結果と同一 であるが、anova プロシジャの場合と異なり、個々の要因の主効果や交互作用の検定結 果として、Type I 平方和と Type III 平方和、及び anova ではなかった contrast 文による対比検定結果も出力している。ただし、このデータは釣り合い型、すなわち 各水準のサンプル数が等しいので、２種類の平方和は全く同一である。
contrast 文による対比検定は、glm プロシジャに特有のもので、クライアントが自由 に指定する非対比較（対比較も含む）とその検定を行える。ただし、SAS の contrast 文の指定による対比検定では、いわゆる族あたりの危険率のコントロールはなされない ので注意が必要である。この場合、要因 A の主効果は有意でないので、最初の２つの 対比 c1(a) 及び c2(a) の検定結果が有意でないのは当然である。
(2-2) の項での means 文による平均値の出力と tukey 法による対比検定結果は、 anova プロシジャによるそれらと同一であるので、ここでは省略する。
CRF-33 ANOVA by proc glm
General Linear Models Procedure
Least Squares Means
A Y
LSMEAN
1 35.3333333
2 32.3333333
3 37.3333333
B Y
LSMEAN
1 27.6666667
2 35.3333333
3 42.0000000
A B Y
LSMEAN
1 1 33.0000000
1 2 35.0000000
1 3 38.0000000
2 1 30.0000000
2 2 31.0000000
2 3 36.0000000
3 1 20.0000000
3 2 40.0000000
3 3 52.0000000
|
(2-2) の項の lsmeans 文による、当該データの各水準での平均値は、水準間で サンプル数がすべて等しいので、(2-1) の項の means 文による各水準での平均値と 同一となる。
means plot utilizing proc glm/lsmeans
プロット : LSMEAN*A. 使用するプロット文字 : B の値．
LSMEAN |
60 +
|
|
| 3
50 +
|
|
|
40 + 2
|3
|2 3
|1
30 + 1
|
|
|
20 + 1
|
-+----------------------------------+----------------------------------+-
1 2 3
treatment A
NOTE: 1 オブザベーションを表示してません．
|
上の出力結果は lsmeans 文により得られた結果を (3) の plot プロシジャにより出 力したものである。
test for contrast-contrast interactions
General Linear Models Procedure
Class Level Information
Class Levels Values
A 3 1 2 3
B 3 1 2 3
Number of observations in data set = 45
..................................
..................................
..................................
Contrast DF Contrast SS Mean Square F Value Pr > F
c1(a)-c1(b) 1 1.2500000 1.2500000 0.02 0.8883
c1(a)-c2(b) 1 1.2500000 1.2500000 0.02 0.8883
c1(a)-c3(b) 1 5.0000000 5.0000000 0.08 0.7789
c2(a)-c1(b) 1 570.4166667 570.4166667 9.13 0.0046
c2(a)-c2(b) 1 1170.4166667 1170.4166667 18.73 0.0001
c2(a)-c3(b) 1 106.6666667 106.6666667 1.71 0.1997
|
上の結果は、上記プログラムの (4) の項の出力の一部である。(2-2) の項のプログラ ムと比較するとわかるように、両者は model 文のところまでは同一なので、出力結果も 当然そこまでは同一であり、それらはここでは省略した。(2-2) の項の glm プロシジャ と異なる点は、(4) の項での contrast 文による対比検定が、２要因の（全体的）交互 作用の部分交互作用の１種である対比・対比交互作用（６つ）の検定である点である。 ここでの対比・対比交互作用は、(2-2) の項における要因 A の２つの対比と同じく 要因 B の３つの対比から構成される６つの部分交互作用である点に注意されたい。
これらの対比・対比交互作用は、プログラムのところでも指摘したように、 (2-2) の項の要因 A の対比の係数行列（２行３列）と同要因 B の対比の係数行列（３ 行３列）の クロネッカー積を指定することにより、 (4) の項のような contrast 文にすることによって検討できる。クロネッカー 積の計算方法は、つぎの 1.4.7 節に定義と具体例を示してあるので、参照のこと。
sas original and Scheffe adjusted F tests
OBS _NAME_ _SOURCE_ _TYPE_ DF SS F PROB FSCHEF PSCHEF
1 Y c1(a)-c1(b) CONTRAST 1 1.25 0.0200 0.88833 0.00500 0.99995
2 Y c1(a)-c2(b) CONTRAST 1 1.25 0.0200 0.88833 0.00500 0.99995
3 Y c1(a)-c3(b) CONTRAST 1 5.00 0.0800 0.77892 0.02000 0.99918
4 Y c2(a)-c1(b) CONTRAST 1 570.42 9.1267 0.00462 2.28167 0.07947
5 Y c2(a)-c2(b) CONTRAST 1 1170.42 18.7267 0.00011 4.68167 0.00381
6 Y c2(a)-c3(b) CONTRAST 1 106.67 1.7067 0.19970 0.42667 0.78836
|
上の結果で、F から PROB までの項は、8. の SAS による contrast 文による 族あたりの危険率がコントロールされない場合の 結果である。その右側の FSCHEF と PSCHEF がシェッフェ方式で危険率をコントロール した場合の F 値と PROB の値である。例えば、第４番目の対比・対比交互作用は、 コントロール前では、１パーセント以上の高い水準で有意であるが、コントロール後 では、有意ではなくなっていることがわかる。
histogram of the error estimates
Frequency
| *****
| *****
| ***** *****
10 + ***** *****
| ***** *****
| ***** *****
| ***** ***** *****
| ***** ***** ***** *****
5 + ***** ***** ***** ***** *****
| ***** ***** ***** ***** *****
| ***** ***** ***** ***** ***** *****
| ***** ***** ***** ***** ***** *****
| ***** ***** ***** ***** ***** *****
-------------------------------------------------------------------------
-8 -4 0 4 8 12
RESID Midpoint
|
上の結果は、プログラム (4) の項の glm プロシジャの中の output 文で指定した 一時ファイルの出力した結果を読み込み、chart プロシジャで誤差推定値のヒストグラ ムを描いたものである。
test for global normality of the dependent variable
Univariate Procedure
Variable=RESID
Moments
N 45 Sum Wgts 45
Mean 0 Sum 0
Std Dev 7.150969 Variance 51.13636
Skewness 0.096003 Kurtosis -1.23423
USS 2250 CSS 2250
CV . Std Mean 1.066004
T:Mean=0 0 Pr>|T| 1.0000
Num ^= 0 39 Num > 0 19
M(Sign) -0.5 Pr>=|M| 1.0000
Sgn Rank 1 Pr>=|S| 0.9890
W:Normal 0.899259 Pr < W 0.0006
Quantiles(Def=5)
100% Max 12 99% 12
75% Q3 8 95% 10
50% Med 0 90% 10
25% Q1 -6 10% -10
0% Min -10 5% -10
1% -10
Range 22
Q3-Q1 14
Mode 0
Extremes
Lowest Obs Highest Obs
-10( 9) 10( 26)
-10( 25) 10( 33)
-10( 21) 10( 39)
-10( 20) 11( 34)
-10( 6) 12( 45)
|
上の出力結果の中程の、W:Normal の行の Pr < W の値を見ると 0.0006 となって おり、変数 resid のデータからは、データは正規母集団からの無作為標本である、 という帰無仮説が１パーセント以上の高い有意水準で棄却されることがわかる。
もっとも、1.2 節で指摘したように、通常の分散分析では正規性の仮定が多少崩れて いても、F-検定は頑健であることが知られているので、この結果にあまり深刻になる 必要はない。
Bartlett test for homogeneity of variances
Discriminant Analysis
45 Observations 44 DF Total
1 Variables 36 DF Within Classes
9 Classes 8 DF Between Classes
..............................
..............................
..............................
Discriminant Analysis Test of Homogeneity of Within Covariance Matrices
Notation: K = Number of Groups
P = Number of Variables
N = Total Number of Observations - Number of Groups
N(i) = Number of Observations in the i'th Group - 1
__ N(i)/2
|| |Within SS Matrix(i)|
V = -----------------------------------
N/2
|Pooled SS Matrix|
_ _ 2
| 1 1 | 2P + 3P - 1
RHO = 1.0 - | SUM ----- - --- | -------------
|_ N(i) N _| 6(P+1)(K-1)
DF = .5(K-1)P(P+1)
_ _
| PN/2 |
| N V |
Under null hypothesis: -2 RHO ln | ------------------ |
| __ PN(i)/2 |
|_ || N(i) _|
is distributed approximately as chi-square(DF)
Test Chi-Square Value = 1.360716
with 8 DF Prob > Chi-Sq = 0.9948
Since the chi-square value is not significant at the 0.05 level,
a pooled covariance matrix will be used in the discriminant function.
..............................
..............................
..............................
|
最後の結果は、データの誤差項の正規性と並ぶ 1.2 節で述べた分散分析のもう１つ の仮定である、セル間の分散の等質性の仮定の discrim プロシジャによる検定結果で ある。SAS では、分散分析にこのオプションを持っていないので、既に示したように discrim プロシジャを利用しなければならない。このプロシジャは、もともと判別 分析のためのものであるので、ここでのプログラムのような使い方をすると、当然 分散の等質性以外のいろいろな判別分析に特有な情報も出力するが、ここでは、それら は無視すればよい。クライアントは、うえのように途中の出力結果は飛ばして、カイ 二乗検定結果のみを見ればよい。SAS では、上のように、Bartlett の検定を用いて 分散の等質性の検定を行う。ここでは、Prob > Chi-Sq = 0.9948 となっているの で、分散の等質性の帰無仮説は採択される。
このプログラムは、以下のような７つの部分から成り立っている：
| 1.冒頭のコメント |
| 2.(1)データ入力部プログラム |
| 3. (2)glm プロシジャによる非釣り合い型 CRF-23 プログラム |
| 4. (3) 各セルの最小２乗平均値のプロットプログラム |
| 5. (4-1) モデルの誤差推定値のヒストグラムの作成 |
| 6. (4-2) モデルの誤差推定値の分布の正規性検定のためのプログラム |
| 7. (4-3) セルの分散の等質性の検定のためのプログラム |
*-------------------------------------------------------------------------* | December 9, 1995 | | sas program -- crf23toh.sas | | | | Example of an unbalanced and non-orthogonal CRF-pq ANOVA design | | analysis for Takeuchi (ed.)/Takahashi-Ohashi-Haga (1990) Table 15.7 | | data (p.297). This data has 2 levels for factor A, and 3 levels for | | factor B with 14 subjects, i.e., CRF-2 3 ANOVA design. | | | | Cautions: | | | | (1) According to Winer (1971), Bartlett's test for homogeneity of vari- | | ance (Bartlett, 1937) is the most widely used test. However, the rou- | | tine use of it as a prelliminary test is not recommended. For, it is | | extremely sensitive to nonnormality (Scheffe, 1959, p.83). | | (2) Test for homogeneity of variance may be unnecessary for balanced | | designs. | | (3) Proc anova must be used only for balanced designs. | | (4) For pairwise comparisons, the tukey procedure is superior to the | | Scheffe procedure. | | (5) Use cldiff option in the mean statement for all pairwise compari- | | sons. | | (6) CLDIFF is the default for unbalanced designs unless DUNCAN, REGWF, | | REGWQ, SNK, or WALLer is specified. | | (7) Lsmeans option is not available with proc anova. | | (8) Proc glm can be used for balanced or unbalanced designs. | | (9) Contrast statement is not available with proc anova. | | (a) Type II SS is recommended in case of extremely unbalanced designs. | | | *-------------------------------------------------------------------------*; |
options ps=60 ls=80;
/* (1) data input */
data work;
input (a b) (1.) +1 y 2. ;
label a='level of factor A'
b='level of factor B'
y='dependent variable';
cards;
11 10
11 13
21 15
21 14
12 14
12 12
12 15
12 11
22 16
22 18
13 22
13 19
23 21
23 18
;
|
/* (2) unbalanced & nonorthogonal CRF-23 ANOVA using proc glm */ options ps=60 ls=80; title 'CRF-23 ANOVA by proc glm'; proc glm data=work; class a b; model y=a|b/ss1 ss2 ss3 ss4; means a b/tukey alpha=0.01; /* means a b/tukey cldiff alpha=0.01; */ lsmeans a|b/ out=glmout; output out=temp residual=resid; run; |
/* (3) means plot using proc glm/lsmeans output option */ options ps=30 ls=80; data work2; set glmout; if a=. or b=. then delete; run; title 'means plot utilizing proc glm/lsmeans'; proc plot data=work2; plot lsmean*a=b; run; |
/* (4-1) posterior analysis (1): histogram plot of error estimates */ options ps=30 ls=80; title 'histogram of the error estimates '; proc chart data=temp; vbar resid; run; |
/* (4-2) posterior analysis (2): test for normality of error estimates */ options ps=60 ls=80; title 'test for global normality of the dependent variable'; proc univariate data=temp normal; var resid; run; |
/* (4-3) posterior analysis (3) : test for homogeneity of variances */
data work4;
set work;
ic=2; /* ic=number of levels of factor a */
jc=3; /* jc=number of levels of factor b */
do i=1 to ic;
do j=1 to jc;
if a=i and b=j then group=jc*(i-1)+j;
end;
end;
run;
title 'Bartlett test for homogeneity of variances';
proc discrim data=work4 pool=test slpool=0.05;
class group;
var y;
run;
|
| crf23toh.sas |
(2) のプログラムによる最初の主要な結果として、つぎのような分散分析表 が出力される：
CRF-23 ANOVA by proc glm
General Linear Models Procedure
Dependent Variable: Y dependent variable
Sum of Mean
Source DF Squares Square F Value Pr > F
Model 5 145.42857143 29.08571429 8.95 0.0039
Error 8 26.00000000 3.25000000
Corrected Total 13 171.42857143
R-Square C.V. Root MSE Y Mean
0.848333 11.57746 1.8027756 15.571429
Source DF Type I SS Mean Square F Value Pr > F
A 1 21.42857143 21.42857143 6.59 0.0332
B 2 108.80000000 54.40000000 16.74 0.0014
A*B 2 15.20000000 7.60000000 2.34 0.1586
Source DF Type II SS Mean Square F Value Pr > F
A 1 16.13333333 16.13333333 4.96 0.0565
B 2 108.80000000 54.40000000 16.74 0.0014
A*B 2 15.20000000 7.60000000 2.34 0.1586
Source DF Type III SS Mean Square F Value Pr > F
A 1 13.09090909 13.09090909 4.03 0.0797
B 2 105.20000000 52.60000000 16.18 0.0015
A*B 2 15.20000000 7.60000000 2.34 0.1586
Source DF Type IV SS Mean Square F Value Pr > F
A 1 13.09090909 13.09090909 4.03 0.0797
B 2 105.20000000 52.60000000 16.18 0.0015
A*B 2 15.20000000 7.60000000 2.34 0.1586
|
このデータは非釣り合い型デザインなので、上のプログラムでは model 文のオプ ションとして、ss1 ss2 ss3 ss4 のすべてを指定しているので、出力結果にもこれら ４つの平方和がすべて出ている。1.3.7 節で既に述べたように、非釣り合い型デザイン データで、交互作用を含むデザインの場合、４つの平方和は微妙に異なる。このような 場合、そこでも述べたように高橋ら (1990) は Type II 平方和を進めている。
上の結果では、Type II 平方和で見てみると、要因 A は傾向程度、要因 B は１パ ーセント以上の高い水準で有意となっているが、交互作用は有意ではない。
(2) の glm プロシジャの分析結果の後半は、means 文による対比検定結果である。 ここでは検定方法として tukey 法が指定されているので、tukey 法による出力がなさ れている：
CRF-23 ANOVA by proc glm
General Linear Models Procedure
Tukey's Studentized Range (HSD) Test for variable: Y
NOTE: This test controls the type I experimentwise error rate, but
generally has a higher type II error rate than REGWQ.
Alpha= 0.01 df= 8 MSE= 3.25
Critical Value of Studentized Range= 4.744
Minimum Significant Difference= 3.2659
WARNING: Cell sizes are not equal.
Harmonic Mean of cell sizes= 6.857143
Means with the same letter are not significantly different.
Tukey Grouping Mean N A
A 17.0000 6 2
A
A 14.5000 8 1
General Linear Models Procedure
Tukey's Studentized Range (HSD) Test for variable: Y
NOTE: This test controls the type I experimentwise error rate.
Alpha= 0.01 Confidence= 0.99 df= 8 MSE= 3.25
Critical Value of Studentized Range= 5.635
Comparisons significant at the 0.01 level are indicated by '***'.
Simultaneous Simultaneous
Lower Difference Upper
B Confidence Between Confidence
Comparison Limit Means Limit
3 - 2 1.030 5.667 10.304 ***
3 - 1 1.920 7.000 12.080 ***
2 - 3 -10.304 -5.667 -1.030 ***
2 - 1 -3.304 1.333 5.970
1 - 3 -12.080 -7.000 -1.920 ***
1 - 2 -5.970 -1.333 3.304
|
上の結果は、要因 B については、いわゆる総当たり方式による結果が出ている。 プログラムの冒頭のコメントの (6) で指摘したように、SAS では非釣り合い型デザイ ンデータの場合、２、３の特別なケースを除き、tukey 法は総当たり方式による出力 となる。要因 A については２水準であるので、もともと tukey 法などによる対比較 も不要で主効果を見れば足りる。
CRF-23 ANOVA by proc glm
General Linear Models Procedure
Least Squares Means
A Y
LSMEAN
1 15.0000000
2 17.0000000
B Y
LSMEAN
1 13.0000000
2 15.0000000
3 20.0000000
A B Y
LSMEAN
1 1 11.5000000
1 2 13.0000000
1 3 20.5000000
2 1 14.5000000
2 2 17.0000000
2 3 19.5000000
|
うえの結果から、例えば要因 A の２水準の平均値を上記 means 文による tukey 法の結果に付随して出ている値と lsmean 文により計算された平均値と比べるとわか るように、多少異なることがわかる。
means plot utilizing proc glm/lsmeans
プロット : LSMEAN*A. 使用するプロット文字 : B の値．
25 +
|
|
|
|
| 3
20 +
| 3
|
LSMEAN |
| 2
|
15 +
| 1
| 2
|
| 1
|
10 +
---+-----------------------------------------------+--
1 2
level of factor A
|
うえのプロット結果は、(3) の data ステップと plot プロシジャによるセル間の 最小２乗平均値のプロット結果である。
histogram of the error estimates
Frequency
5 + *****
| *****
| *****
| *****
4 + ***** *****
| ***** *****
| ***** *****
| ***** *****
3 + ***** ***** *****
| ***** ***** *****
| ***** ***** *****
| ***** ***** *****
2 + ***** ***** *****
| ***** ***** *****
| ***** ***** *****
| ***** ***** *****
1 + ***** ***** ***** ***** *****
| ***** ***** ***** ***** *****
| ***** ***** ***** ***** *****
| ***** ***** ***** ***** *****
--------------------------------------------------------------------
-2 -1 0 1 2
RESID 中間点
|
Univariate Procedure
Variable=RESID
Moments
N 14 Sum Wgts 14
Mean 0 Sum 0
Std Dev 1.414214 Variance 2
Skewness 0 Kurtosis -1.80638
USS 26 CSS 26
CV . Std Mean 0.377964
T:Mean=0 0 Pr>|T| 1.0000
Num ^= 0 14 Num > 0 7
M(Sign) 0 Pr>=|M| 1.0000
Sgn Rank -1.5 Pr>=|S| 0.9435
W:Normal 0.887012 Pr < W 0.0730
.............................
.............................
.............................
|
うえの結果から、誤差項についての正規性の帰無仮説は棄却出来ないことがわかる。
Bartlett test for homogeneity of variances
Discriminant Analysis
14 Observations 13 DF Total
1 Variables 9 DF Within Classes
5 Classes 4 DF Between Classes
...........................
...........................
...........................
Discriminant Analysis Test of Homogeneity of Within Covariance Matrices
Notation: K = Number of Groups
P = Number of Variables
N = Total Number of Observations - Number of Groups
N(i) = Number of Observations in the i'th Group - 1
__ N(i)/2
|| |Within SS Matrix(i)|
V = -----------------------------------
N/2
|Pooled SS Matrix|
_ _ 2
| 1 1 | 2P + 3P - 1
RHO = 1.0 - | SUM ----- - --- | -------------
|_ N(i) N _| 6(P+1)(K-1)
DF = .5(K-1)P(P+1)
_ _
| PN/2 |
| N V |
Under null hypothesis: -2 RHO ln | ------------------ |
| __ PN(i)/2 |
|_ || N(i) _|
is distributed approximately as chi-square(DF)
Test Chi-Square Value = 1.555951
with 4 DF Prob > Chi-Sq = 0.8167
Since the chi-square value is not significant at the 0.05 level,
a pooled covariance matrix will be used in the discriminant function.
...........................
...........................
...........................
|
セル間の分散の等質性の帰無仮説は、うえのカイ２乗検定の結果、棄却できないこと がわかる。