このページは、令和２年４月２５日に一部更新しました。
| 1.5.6.1 SAS/IML による反復測度 RB-p デザインの GMANOVA 分析 プログラムの概要 |
| 1.5.6.2 SAS/IML による反復測度 RB-p デザインの GMANOVA 分析 についての出力結果の概要 |
また、この節にはつぎの２つの SAS プログラムのダウンロードコーナーを用意 してあります：
| 1.反復測度RB-pデザインのIMLプログラム例 |
この節では、1.5.3 節で述べたホテリングの
T2 統計量を用いた反復測度 RB-p デザインデータの SAS プログラム
例を示す。前節の後半で指摘したように、
反復測度 RB-p デザインデータの SAS プログラム例の中の glm プロシジャに
よる repeated 文の contrast 項等を用いるか、manova 文を h=intercept と組
み合わせて用いると、MANOVA 出力における Wilks' Lambda などの検定基準に対
する F-値が出力されるが、４つの基準値はすべて一致し、この値はホテリングの
T2 統計量を F-値に直したものに相当しているのである。つまり、既
に前節で反復測度 RB-p デザインの MANOVA 分析（正確には GMANOVA 分析）の一
部を行っているのである。
しかし、前節のプログラムのみでは、contrast 項や manova 文の m= の項で
指定した対比の族あたりの危険率のコントロールは、Scheff\'e 修正でしかなか
った（より正確には、この修正に加え、SAS では反復測定要因の水準間の
contrast 文による対比検定においては、誤差分散の不偏推定値 UE
に関しては、個々の対比の分散の不偏推定値で置き換えている）。
実は、ホテリングの T2 統計量を用いた反復測度 RB-p デザインに
おける同時検定方式を用いれば、1.5.3 節の最後で述べたようにより単純ですっき
りとした族あたりの危険率のコントロールができるはずである。しかし、これを
glm プロシジャのみで実行することはできない。これを行うには、つぎのように
SAS/IML を用いればよい。
データは、前節の結果と比較可能なように、前節で用いた鏡映描写実験データ
である。ただし、プログラムを短くするため、ここでは前節のデータを一旦永久
SAS ファイル化し、それ (permfile.repmirr1) を入力しているので、注意せよ。
プログラムは、以下に示すように：
| 1.冒頭のコメント |
| 2.データ入力 |
| 3. 標本共分散行列・相関行列の計算モジュール covcor |
| 4. T2統計量の計算と検定モジュール hotcon |
| 5.モジュールの実行と終了 |
の５項から成る。
*----------------------------------------------------------------------* | December 1, 1998 | | | | file name: $HOME/sasprog/multivar/imlmirr.sas | | | | SAS/IML program for executing simultaneous contrast tests using | | Hotelling's T-squared statistics for a repeated measures RB-p design | | data. | | | *----------------------------------------------------------------------*; |
この部分は、このプログラムの簡単な内容紹介の部分である。
options pagesize=60 ls=80;
libname permfile '$HOME/sasset/multivar';
proc iml;
use permfile.repmirr1;
read all var{exer1 exer2 exer3 exer4 exer5 exer6 exer7}
into x;
print "Input raw-data matrix, X", x;
nr=nrow(x);
nc=ncol(x);
nr1=nr-1;
nc1=nc-1;
if nr1 <= 0 then do;
print "number of samples is too small"; return;
end;
|
この部分では、iml を起動し、あらかじめ作成してある永久 SAS ファイル permfile.repmirr1 （前節の２５名の被験者に対する鏡映描写実験での７回の練習 試行のデータを、SAS 名、exer1 から exer7 として提起してあるものとする） を行列 x に読み込みつぎからの解析のための若干の変数を定義するためのもの である。ここで、例えば nrow(x) や ncol(x) は iml の一種の組み込み関数で、 それぞれ、カッコ内の行列の行数や列数を左辺の変数に返すためのものである。
start covcor;
sum=x[+,];
sscp=x`*x-sum`*sum/nr;
scov=sscp/nr1;
sd=vecdiag(scov);
do j=1 to nc;
if sd[j]<=0.0001 then do;
print "sample unbiased sd is too small", j, sd; return;
end;
end;
print "Sample covariance matrix", scov[rowname=name colname=name];
vd=vecdiag(sscp);
d=diag(1/sqrt(vd));
r=d*sscp*d;
print "Correlation matrix", r[rowname=name colname=name];
finish;
|
この部分は、iml により入力データ行列から標本共分散行列と相関行列を計算する ための１種のサブルーチンプログラムであり、SAS ではモジュールと呼んでいる。 最初の sum=x[+,]; は、行列 x の行和を要素とする行ベクトルを計算するための ものである。また、sscp=x`*x-sum`*sum/nr; における x や sum の右側のバック クオテーションマークは、それぞれの行列やベクトルの転置行列や転置ベクトルを 指示する。関数 vecdiag(scov) は、カッコ内の正方行列 scov の対角要素を要素と する列ベクトルを生成するためのものである。また、一般に vec なる名前の変数を ベクトルとすると、関数 diag (vec) はベクトル vec の要素を対角要素とする対角 行列を生成するためのものである。
start hotcon;
xbar=sum/nr;
print "means of repeated measures", xbar;
c={1 -1 0 0 0 0 0,
0 1 -1 0 0 0 0,
0 0 1 -1 0 0 0,
0 0 0 1 -1 0 0,
0 0 0 0 1 -1 0,
0 0 0 0 0 1 -1};
d=c*xbar`;
d2=d##2;
scovt=c*scov*c`;
s=inv(scovt);
t2=nr*d`*s*d; /* Hotelling T-squared statistic */
vecs=vecdiag(scovt);
g=nr*diag(1/vecs);
ddf=nr-nc+1;
const=ddf/nr1/nc1;
fmain=const*t2;
pmain=1-probf(fmain,nc1,ddf);
print "Hotelling T-squared test for the main effect",
fmain[colname=namef],
nc1[colname=namedf1], ddf[colname=namedf2],
pmain[colname=namep];
h=const*g*d2;
pcon=vecdiag(I(nc1))-probf(h,nc1,ddf);
print "F values of simultaneous contrast test",
h[rowname=namecon colname=namef];
print "p-values of simultaneous contrast test",
pcon[rowname=namecon colname=namep];
finish;
name={ex1 ex2 ex3 ex4 ex5 ex6 ex7};
namecon={contrast1 contrast2 contrast3 contrast4
contrast5 contrast6 contrast7};
namef={F_value};
namep={Prob_ge_F};
namedf1={df_1};
namedf2={df_2};
|
この部分は、対比行列やデータからあらかじめ計算された共分散行列等を用いて ホテリングの T2統計量やそれを用いた検定の p-値を計算するための モジュールである。ここで、最初の方の変数 sum はそれ以前のモジュール covcor で定義された行ベクトルであることに注意したい。また、途中の d2=d##2; における ## は行列 d の各要素を２乗することを指示する「べき乗」作成のための 演算子である。最後の方の、vecdiag のカッコ内の I(nc1) は、次数 nc1 の単位行列 （対角要素のみ１で、残りはゼロの正方行列）の生成を指示する関数である。
run covcor; run hotcon; quit; |
「run モジュール名」で、当該モジュール（サブルーチン）の実行を指示する。 最後の quit; で iml を終了する。
| imlmirr.sas |
うえのプログラムを実行すると、以下のような出力結果が得られる：
SAS システム 5
Input raw-data matrix, X
X
81 48 51 47 31 41 28
48 36 28 31 26 22 22
47 62 35 35 34 25 25
55 28 35 31 28 23 23
103 112 67 70 68 42 39
53 33 28 28 22 23 20
38 44 25 22 21 19 20
55 49 42 35 38 37 33
56 34 39 27 27 20 26
62 61 47 45 43 44 40
127 102 96 102 76 70 56
37 35 29 38 35 31 28
36 36 32 29 34 23 24
41 42 42 35 34 27 29
38 29 29 26 22 20 18
50 34 27 22 23 15 18
29 28 26 29 32 25 25
91 52 33 34 40 28 33
89 51 40 41 35 32 26
41 44 41 60 39 34 30
28 29 27 25 24 23 22
34 31 28 27 24 26 19
33 24 25 21 23 16 19
41 35 35 32 30 27 25
36 38 29 30 32 24 24
|
この結果は、データ入力部のプログラム直後の print 文による。
Sample covariance matrix
SCOV EX1 EX2 EX3 EX4 EX5 EX6 EX7
EX1 638.45667 439.65333 328.435 346.87 258.94333 220.90333 163.57833
EX2 439.65333 455.22667 288.02167 321.335 258.04667 194.56 149.75167
EX3 328.435 288.02167 243.67333 261.63833 186.20667 166.23 120.01333
EX4 346.87 321.335 261.63833 315.44333 216.12167 190.71 134.81833
EX5 258.94333 258.04667 186.20667 216.12167 173.24 133.08833 104.28833
EX6 220.90333 194.56 166.23 190.71 133.08833 133.89333 91.71
EX7 163.57833 149.75167 120.01333 134.81833 104.28833 91.71 71.943333
Correlation matrix
R EX1 EX2 EX3 EX4 EX5 EX6 EX7
EX1 1 0.8155129 0.832683 0.7729303 0.778601 0.7555393 0.7632464
EX2 0.8155129 1 0.8647833 0.8479759 0.9188832 0.7880617 0.8274896
EX3 0.832683 0.8647833 1 0.943706 0.9062899 0.9202926 0.9064215
EX4 0.7729303 0.8479759 0.943706 1 0.9245141 0.9279687 0.8949381
EX5 0.778601 0.9188832 0.9062899 0.9245141 1 0.8738491 0.9341494
EX6 0.7555393 0.7880617 0.9202926 0.9279687 0.8738491 1 0.9344196
EX7 0.7632464 0.8274896 0.9064215 0.8949381 0.9341494 0.9344196 1
|
これらの結果、とりわけ後者から、このデータにおける反復測定水準間には相互に かなり高い相関があることがわかる。
means of repeated measures
XBAR
53.96 44.68 37.44 36.88 33.64 28.68 26.88
Hotelling T-squared test for the main effect
FMAIN F_VALUE
7.7963094
NC1 DF_1
6
DDF DF_2
19
PMAIN PROB_GE_F
0.0002481
F values of simultaneous contrast test
H F_VALUE
CONTRAST1 1.3251028
CONTRAST2 1.407374
CONTRAST3 0.0288628
CONTRAST4 0.6135276
CONTRAST5 1.9813895
CONTRAST6 0.4767658
p-values of simultaneous contrast test
PCON PROB_GE_F
CONTRAST1 0.2941361
CONTRAST2 0.2628
CONTRAST3 0.999867
CONTRAST4 0.7168598
CONTRAST5 0.119097
CONTRAST6 0.8172111
|
ここでは、モジュール hotcon の print 文に従い、まず反復測定各水準の平均値、
つぎにホテリングの T2統計量による主効果の F-検定統計量と自由度
及び p-値、危険率をコントロールした各対比検定のための F-統計量と対応する
p-値、が出力される。
ホテリングの T2統計量による主効果の F-検定における F-統計量の
値や自由度、p-値はすべて前節の repeated 文を用いた glm プロシジャによる出力
結果のうちの (G)MANOVA による主効果の検定結果と同一である。
一方、対比検定結果を見ると、同時検定方式で族あたりの危険率のコントロール
を行うと、前節で示した単純な SAS の対比検定では幾つかが高い水準で有意であった
にもかかわらず、５パーセント以上の水準で有意な対比は１つも無くなっていること
がわかる。ここでは、かろうじて Contrast 5 のみが有意な傾向を示すのみである。
もっとも、この結果は、前節で示した Scheffe 方式（で、かつ UEに
個別対比の誤差分散の不偏推定値をあてる方式）による結果に近い。