このページは、令和２年４月２５日に一部更新しました。
ここでは、単純な RB-p デザイン ANOVA データの SAS による分散分析を以下の ２つのセクションに分けて示す：
| 1.5.4.1 単純な RB-p デザイン ANOVA プログ ラムとその概要 |
| 1.5.4.2 単純な RB-p デザイン ANOVA プログラム の出力結果の概要 |
| 単純なRB-pデザインのプログラム例 |
この節では、まず 1.5.1 節で述べた単純な乱塊法デザイン、とりわけ混合 模型の場合の SAS による ANOVA 分析プログラム例を示す。したがって、主 要因はこの場合反復測度からは構成されていないものとする。ただし、その ような単純な乱塊法デザインデータで適切なものが見あたらないので、ここでは Kirk (1982, p.88) の CR-4 デザインデータを RB-4 デザインデータとみなして SAS プログラム例として掲載することにする。もちろん、一般的に言って、 CR-4 デザインデータを RB-4 デザインデータとみなすこと自身は決して正しい 見方ではないが、あくまでもここでは単純な乱塊法デザインデータを手に した時、どのようなプログラムを書けばよいかを示すために架空の例として あげてあることに注意されたい。
プログラムは、以下に示すように：
の７項から成る。単に、RB デザインの分散分析のみを行うには、これらの内の (1) 及び (2-1) だけを実行すればよいが、RB デザインがらみの各種の分析を行うには この一連のプログラムを実行するとよい。
*-------------------------------------------------------------------------* | December 19, 1995 | | sas program -- rb4test.sas | | | Example of a RB-p ANOVA design analysis for Kirk (1982) data, p,88. | | exercises 27. This is not a RB-p design data, but a CR-p ANOVA design.| | We assume it to be an RB-p design data in order to examine the RB-p | | design data. This data has 4 levels of legibility of a dial in a air- | | plane. Twenty pilots and navigators were randomly assigned to one of | | four groups. We also assume that 5 blocks exist for each levels of | | legibility of a dial. Thus, these blocks are artificial. | | | | Cautions: | | | | (1) In general, RB-I ANOVA design requires IJ subjects, where J indi- | | cates the number of levels of block factor B. | | | *-------------------------------------------------------------------------*; |
この部分は、このプログラムにおけるデータのソースと内容についての説明等を 示したものである。
options ps=60 ls=80;
/* (1) data input */
data work;
do a=1 to 4;
do block=1 to 5;
input y @;
output;
end;
end;
label a='legibility factor with 4 levels'
block='artificial'
y='number of reading error';
cards;
5 11 7 7 5
3 6 3 5 3
15 8 10 6 13
10 6 6 12 6
;
|
この部分は、主要因子が４水準でブロック因子が５水準からなるデータが、うえの ように並べらているデータに対する do 文を用いた特殊な入力プログラムである。 また、２つ目の do ループでは、各列に対応するブロックの番号を block なる変数 に指定しながらデータを読むことを指示している。
ここで、データは必ずしも上のような並びになっている必要はない。最も単純な データ入力の方式は、各被験者毎例えば、ブロック因子の水準、主要因子の水準、 測定値、の順に並べたものを被験者分行に並べるもので、その場合にはもちろん、 上のような do 文は不要である。これに対して、上のようにデータを配列すると、 測定値とデザインとの関係が一見してわかりやすい、というメリットがある。
/* (2-1) RB-4 ANOVA using proc anova */ title 'RB-4 ANOVA for the Kirk (1982) data'; proc anova data=work; class block a; model y=a block; means a/tukey alpha=0.01; /* means a/tukey cldiff alpha=0.01; */ run; |
ここでは、単純な RB-4 デザインデータの ANOVA を指示している。class 文の中に ブロック因子の変数名 block を指定していることに注意せよ。また、ここでの主要 因子は a という変数名であることにも注意せよ。つまり、1.5.1 節の表 1.13 に 示したように、一般に乱塊法 (RB デザイン）では、主要な１要因の効果の有無 のみならず、ブロック因子の効果の有無も検討できることを思い出してほしい。 さらに、主効果が有意な場合、ここでは Tukey 法による対比較により対比検定を行 うよう指示している。
/* (2-2) Means plot using proc means output */ title 'means of each level of the factor'; proc means data=work; var y; class a; output out=temporal mean=aver; run; options ps=30 ls=80; proc plot data=temporal; plot aver*a; run; |
ここでは、means プロシジャによる平均値の計算と plot プロシジャによる平均 値プロットを行わせている。
/* (3-1) RB-4 ANOVA using proc glm */ options ps=60 ls=80; title 'RB-4 ANOVA by proc glm'; proc glm data=work; class a block; model y=a block; means a/tukey alpha=0.01; /* means a/tukey cldiff alpha=0.01; */ lsmeans a/ out=glmout; run; /* (3-2) means plot using proc glm/lsmeans output option */ options ps=30 ls=80; title 'means plot utilizing proc glm/lsmeans'; proc plot data=glmout; plot lsmean*a; run; |
ここでは、主要因子の４水準の平均値を anova プロシジャの means 文や means プロシジャを用いて計算するのではなく、glm プロシジャの lsmeans 文を用いて 計算させ、plot プロシジャで出力するプログラムを示した。means 文と lsmeans 文 の違いについては既に 1.3.6 節で述べたので、ここでは省略する。
/* (4) non-pairwise contrasts using proc glm/contrast statement */ options ps=60 ls=80; title 'pairwise and non-pairwise contrasts'; proc glm data=work outstat=glmstat; class a block; model y=a block; means a/tukey alpha=0.01; contrast 'level 1 vs the other' a 3 -1 -1 -1; contrast 'level 1-2 vs level 3-4' a 1 1 -1 -1; contrast 'level 1-2-3 vs level 4' a 1 1 1 -3; output out=temp residual=resid; run; |
ここでは、glm プロシジャの contrast 文を用いて、主要因子の４水準間の非対・ 比較による対比検定を指示している。同時に、後にモデルの誤差項の正規性の検定の ために、output 文によりモデルの残差項の推定値を resid なる変数名で一時ファイル temp に出力させている。
data work2; set glmstat; if _type_ ne 'CONTRAST' then delete; ic=4; /* ic=number of levels of the major factor, A */ jc=5; /* jc=number of levels of block factor, BL */ ndf=ic-1; ddf=ndf*(jc-1); fschef=f/ndf; pschef=1-probf(fschef,ndf,ddf); drop ndf ddf ic jc; run; title 'sas original and scheffe adjusted F tests'; proc print data=work2; run; |
ここでは、Scheffe 方式により、平均値の非対・比較による対比検定における族あ たりの危険率のコントロールを行い、その結果を出力する。既に、1.4.6 節で述べた ように、SAS の contrast 文では対比検定の検定に際して族あたりの危険率のコント ロールは行わないので、これを行うためには上のような新たなプログラムが必要と なることに注意せよ。
/* (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; /* (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; |
ここでは、誤差推定値のヒストグラムの作成と正規性の検定を行うために、あらか じめ (4) の項で一時ファイル temp に掃き出した各人の誤差推定値を chart プロシ ジャで読み込み、ヒストグラムを描かせ、さらに univariate プロシジャで誤差項の 正規性の検定を行わせている。
| rb4test.sas |
うえのプログラムを実行すると、以下のような結果が出力される：
(2-1) のプログラムによる最初の主要な結果として、つぎのような分散分析表 が出力される：
RB-4 ANOVA for the Kirk (1982) data 2
Analysis of Variance Procedure
Dependent Variable: Y number of reading error
Sum of Mean
Source DF Squares Square F Value Pr > F
Model 7 113.65000000 16.23571429 1.79 0.1794
Error 12 108.90000000 9.07500000
Corrected Total 19 222.55000000
R-Square C.V. Root MSE Y Mean
0.510672 40.98604 3.0124741 7.3500000
Source DF Anova SS Mean Square F Value Pr > F
A 3 105.35000000 35.11666667 3.87 0.0379
BLOCK 4 8.30000000 2.07500000 0.23 0.9170
RB-4 ANOVA for the Kirk (1982) data 3
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= 12 MSE= 9.075
Critical Value of Studentized Range= 5.502
Minimum Significant Difference= 7.4125
Means with the same letter are not significantly different.
Tukey Grouping Mean N A
A 10.400 5 3
A
A 8.000 5 4
A
A 7.000 5 1
A
A 4.000 5 2
|
最初の Model の項の変動因は、このデータの場合、主要因とブロック因子を合わせた 変動分を意味する。したがって、個々の変動因の効果の検定結果は、その数行下に 出力される Source（変動因）の項を見る必要がある。F-検定の結果の p-値より、 a 要因の効果は５パーセント水準で有意と言える。一方、ブロック因子の効果は有意 ではないと言える。
最後の数行は、(2-1) の項のプログラムにおける means 文による主効果の Tukey 法による対比検定結果である。この場合、Tukey Grouping の項のアルファベットの ラベルがすべての水準で A と同じものになっているので、どの対比較も統計的には 有意でないことがわかる。
(2-2) のプログラムにより、主要因子の４水準の平均値等の値とそのプロット結果 がつぎのように出力される：
means of each level of the factor 4
Analysis Variable : Y number of reading error
A N Obs N Mean Std Dev Minimum Maximum
-------------------------------------------------------------------------------
1 5 5 7.0000000 2.4494897 5.0000000 11.0000000
2 5 5 4.0000000 1.4142136 3.0000000 6.0000000
3 5 5 10.4000000 3.6469165 6.0000000 15.0000000
4 5 5 8.0000000 2.8284271 6.0000000 12.0000000
-------------------------------------------------------------------------------
means of each level of the factor 5
プロット : AVER*A. 凡例 : A = 1 OBS, B = 2 OBS, ...
AVER |
12 +
|
|
| A
10 +
|
|
|
8 + A
|
| A
|
6 +
|
|
|
4 + A
|
---+------------------+------------------+------------------+--
1 2 3 4
legibility factor with 4 levels
NOTE: 1 オブザベーションが欠損値です．
|
(3-1) のプログラムにより、主要因子の４水準の平均値等の値とそのプロット結果 が出力される。結果は、この例のようなデータでは、(2-1) 及び (2-2) の項のプログ ラムによるものと同一の結果になるので、ここでは省略する。もっとも、２要因以上 になり非釣り合い型デザインの場合では、anova の means 文による平均値と glm の lsmeans 文による平均値とは一般には異なることに注意したい。
(4) のプログラムにより、主要因子の４水準間の非対・比較による対比検定結果 と、Scheffe 方式によるそれらの検定の族あたりの危険率のコントロール結果が順に つぎのように出力される。なお、ここで対比検定結果の出力の前に、このデータの分散 分析結果と Tukey 法による対比検定結果が出力されるが、(2-1) や (3-1) の項と 同一の結果なので、ここではそれらは省略する：
pairwise and non-pairwise contrasts 14
General Linear Models Procedure
Dependent Variable: Y number of reading error
Contrast DF Contrast SS Mean Square F Value Pr > F
level 1 vs the other 1 0.81666667 0.81666667 0.09 0.7693
level 1-2 vs level 3 1 68.45000000 68.45000000 7.54 0.0177
level 1-2-3 vs level 1 2.81666667 2.81666667 0.31 0.5877
sas original and scheffe adjusted F tests 15
_
S
_ O _ F P
N U T S S
A R Y P C C
O M C P R H H
B E E E D S O E E
S _ _ _ F S F B F F
1 Y level 1 vs the other CONTRAST 1 0.8167 0.08999 0.76932 0.03000 0.99263
2 Y level 1-2 vs level 3 CONTRAST 1 68.4500 7.54270 0.01772 2.51423 0.10787
3 Y level 1-2-3 vs level CONTRAST 1 2.8167 0.31038 0.58769 0.10346 0.95647
|
ここでの出力結果では、前半に glm プロシジャによる族あたりの危険率のコント ロールなしの対比検定結果を示しており、後半に Scheffe 法に基づくこれら対比検 定の族あたりの危険率をコントロールした結果を出力している。
ただし、後半の出力部の左側の F 及び PROB に対応する値は前半の数値で少数点 以下の桁数のみ５桁まで示したものであることに注意されたい。右端の FSCHEF 及び PSCHEF が、これらに対応する族あたりの危険率をコントロールした結果である。両者 を比較すると、２つ目の非対・比較は、危険率をコントロールしない場合には、５ パーセントで有意であるが、コントロールした結果ではもはや有意でなくなることが わかる。
ここでは、(5-1) のchart プロシジャによる、モデルの誤差項の推定値のヒストグ ラムと、(5-2) 節の univariate プロシジャによるその正規性検定結果の一部を示す。 この出力については、
なる正規性検定結果の後の出力部分は省略した：
histogram of the error estimates 16
Frequency
7 + *****
| *****
6 + ***** *****
| ***** *****
5 + ***** *****
| ***** *****
4 + ***** ***** *****
| ***** ***** *****
3 + ***** ***** *****
| ***** ***** *****
2 + ***** ***** ***** *****
| ***** ***** ***** *****
1 + ***** ***** ***** ***** *****
| ***** ***** ***** ***** *****
--------------------------------------------------------------------
-4 -2 0 2 4
RESID 中間点
test for global normality of the dependent variable 17
Univariate Procedure
Variable=RESID
Moments
N 20 Sum Wgts 20
Mean 0 Sum 0
Std Dev 2.394072 Variance 5.731579
Skewness 0.109059 Kurtosis -0.66842
USS 108.9 CSS 108.9
CV . Std Mean 0.535331
T:Mean=0 0 Pr>|T| 1.0000
Num ^= 0 20 Num > 0 9
M(Sign) -1 Pr>=|M| 0.8238
Sgn Rank -2 Pr>=|S| 0.9488
W:Normal 0.961019 Pr < W 0.5681
|
正規性の検定結果である、W:Normal の項の p-値 (0.5681) から、このデータでは モデルの誤差項は正規分布に従う、と結論出来る。
