Eric's color bar icon

このページは、令和2年4月25日に一部更新しました。
under construction icon

1.5.4節 SAS による単純な RB-p デザインの ANOVA 分析

Eric's eye-bar icon

 ここでは、単純な RB-p デザイン ANOVA データの SAS による分散分析を以下の 2つのセクションに分けて示す:

1.5.4.1 単純な RB-p デザイン ANOVA プログ ラムとその概要
1.5.4.2 単純な RB-p デザイン ANOVA プログラム の出力結果の概要

 また、この節にはつぎの SAS プログラムのダウンロードコーナーを用意してある:

単純なRB-pデザインのプログラム例

Eric's eye-bar icon

1.5.4.1 単純なRB-pデザイン ANOVA プログラム とその概要

  この節では、まず 1.5.1 節で述べた単純な乱塊法デザイン、とりわけ混合 模型の場合の SAS による ANOVA 分析プログラム例を示す。したがって、主 要因はこの場合反復測度からは構成されていないものとする。ただし、その ような単純な乱塊法デザインデータで適切なものが見あたらないので、ここでは Kirk (1982, p.88) の CR-4 デザインデータを RB-4 デザインデータとみなして SAS プログラム例として掲載することにする。もちろん、一般的に言って、 CR-4 デザインデータを RB-4 デザインデータとみなすこと自身は決して正しい 見方ではないが、あくまでもここでは単純な乱塊法デザインデータを手に した時、どのようなプログラムを書けばよいかを示すために架空の例として あげてあることに注意されたい。

 プログラムは、以下に示すように:

1.冒頭のコメント
2.(1) データインプット
3.(2-1) anova プロシジャによる RB-4 デザインデー タの分散分析
4.(2-2) means プロシジャによる平均値の計算と plot プロシジャによる平均値プロット
5.glm プロシジャの contrast 文による平 均値の非対・比較による対比検定
6.Scheffe 方式による平均値の非対・比較 による対比検定における族あたりの危険率のコントロール
7.誤差推定値のヒストグラムの作成と正規性の 検定

の7項から成る。単に、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.                          |
|                                                                         |
*-------------------------------------------------------------------------*;

 この部分は、このプログラムにおけるデータのソースと内容についての説明等を 示したものである。

back icon


(1) データインプット

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
;

 この部分は、主要因子が4水準でブロック因子が5水準からなるデータが、うえの ように並べらているデータに対する do 文を用いた特殊な入力プログラムである。 また、2つ目の do ループでは、各列に対応するブロックの番号を block なる変数 に指定しながらデータを読むことを指示している。

 ここで、データは必ずしも上のような並びになっている必要はない。最も単純な データ入力の方式は、各被験者毎例えば、ブロック因子の水準、主要因子の水準、 測定値、の順に並べたものを被験者分行に並べるもので、その場合にはもちろん、 上のような do 文は不要である。これに対して、上のようにデータを配列すると、 測定値とデザインとの関係が一見してわかりやすい、というメリットがある。

back icon


(2-1) anova プロシジャによる RB-4 デザインデータの分散分析

/* (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 デザイン)では、主要な1要因の効果の有無 のみならず、ブロック因子の効果の有無も検討できることを思い出してほしい。 さらに、主効果が有意な場合、ここでは Tukey 法による対比較により対比検定を行 うよう指示している。

back icon


(2-2) means プロシジャによる平均値の計算と plot プロシジ ャによる平均値プロット

/* (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 プロシジャによる平均 値プロットを行わせている。

back icon


(3-1) glm プロシジャの lsmeans 文による平均値の計算と 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;

 ここでは、主要因子の4水準の平均値を anova プロシジャの means 文や means プロシジャを用いて計算するのではなく、glm プロシジャの lsmeans 文を用いて 計算させ、plot プロシジャで出力するプログラムを示した。means 文と lsmeans 文 の違いについては既に 1.3.6 節で述べたので、ここでは省略する。

back icon


(4) glm プロシジャの contrast 文による平均値の非 対・比較による対比検定

/* (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 文を用いて、主要因子の4水準間の非対・ 比較による対比検定を指示している。同時に、後にモデルの誤差項の正規性の検定の ために、output 文によりモデルの残差項の推定値を resid なる変数名で一時ファイル temp に出力させている。

back icon


Scheffe 方式による平均値の非対・比較による対比検 定における族あたりの危険率のコントロール

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 文では対比検定の検定に際して族あたりの危険率のコント ロールは行わないので、これを行うためには上のような新たなプログラムが必要と なることに注意せよ。

back icon


誤差推定値のヒストグラムの作成と正規性の検定

/* (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 プロシジャで誤差項の 正規性の検定を行わせている。

back icon


プログラムのダウンロード・コーナー

Eric's abar10 icon

rb4test.sas

Eric's eye-bar icon

1.5.4.2 単純なRB-pデザイン ANOVA プログラム の出力結果の概要

 うえのプログラムを実行すると、以下のような結果が出力される:

1. 要因の水準情報、サンプル数

2. anova プロシジャによる分散分析結果

 (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 要因の効果は5パーセント水準で有意と言える。一方、ブロック因子の効果は有意 ではないと言える。

 最後の数行は、(2-1) の項のプログラムにおける means 文による主効果の Tukey 法による対比検定結果である。この場合、Tukey Grouping の項のアルファベットの ラベルがすべての水準で A と同じものになっているので、どの対比較も統計的には 有意でないことがわかる。

3. means プロシジャによる平均値の出力及びそのプ ロット結果

 (2-2) のプログラムにより、主要因子の4水準の平均値等の値とそのプロット結果 がつぎのように出力される:

                       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 オブザベーションが欠損値です.

4. glm プロシジャ等による分散分析結果と、lsmeans 文に よる平均値の出力及び (3-2) の plot プロシジャによるそのプロット結果

 (3-1) のプログラムにより、主要因子の4水準の平均値等の値とそのプロット結果 が出力される。結果は、この例のようなデータでは、(2-1) 及び (2-2) の項のプログ ラムによるものと同一の結果になるので、ここでは省略する。もっとも、2要因以上 になり非釣り合い型デザインの場合では、anova の means 文による平均値と glm の lsmeans 文による平均値とは一般には異なることに注意したい。

5. contrast 文による主要因子の4水準間の非対・比較によ る対比検定結果と Scheffe 方式によるそれらの検定の族あたりの危険率のコントロ ール結果

 (4) のプログラムにより、主要因子の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 に対応する値は前半の数値で少数点 以下の桁数のみ5桁まで示したものであることに注意されたい。右端の FSCHEF 及び PSCHEF が、これらに対応する族あたりの危険率をコントロールした結果である。両者 を比較すると、2つ目の非対・比較は、危険率をコントロールしない場合には、5 パーセントで有意であるが、コントロールした結果ではもはや有意でなくなることが わかる。

6. モデルの誤差項の事後分析

 ここでは、(5-1) のchart プロシジャによる、モデルの誤差項の推定値のヒストグ ラムと、(5-2) 節の univariate プロシジャによるその正規性検定結果の一部を示す。 この出力については、

W:Normal  0.961019 Pr < W  0.5681

なる正規性検定結果の後の出力部分は省略した:

                       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) から、このデータでは モデルの誤差項は正規分布に従う、と結論出来る。

Eric's back icon

Eric's color bar icon