SPMで複数の群の比較を行いたい時のDesign Matrix

SPM-MLに勉強になる話題が流れていたので、共有します。

出典はこちら。
https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=SPM;a65dd354.1307

コントロールAと疾患B、疾患Cという3群がある時、Design Matrixに全部を入れてしまった方がいいのか、
もしくは、コントロールと疾患A、コントロールと疾患Bというように、別々にDesign Matrixを作った方がいいのかという質問です。

これに対し、Cyril Pernetが非常にわかりやすいたとえを使って説明をしています。
ただ、少しだけtypoがあったのでそこを補足して説明します。

以下はMatlabを実際に動かしてみながらやってみるといいでしょう。

非常にわかりやすくするために、各群3人だとします。
そして、あるボクセルのデータが以下のようになっているとします。

A=[9; 10; 11];
B=[19; 20; 21];
C=[29; 30; 31];

3群を全部ひとつのモデルに入れることにしましょう。
この時、GLMで考えると、Yはデータ行列になりますので、以下のようになります。

Y=[A; B; C]
Y =
     9
    10
    11
    19
    20
    21
    29
    30
    31

この時、Design Matrixは次のようにあらわされます。(今はわかりやすく共変量はすべて1としています。)

X=[kron(eye(3), ones(3,1)) ones(9,1)]
X =
     1     0     0     1
     1     0     0     1
     1     0     0     1
     0     1     0     1
     0     1     0     1
     0     1     0     1
     0     0     1     1
     0     0     1     1
     0     0     1     1

GLMでは、Y=XB+Eであらわされます。
そして、Bを求めたい時は、Matlabで以下のコマンドを叩くことで簡単に求まります。
ただ、今、上でBという変数を使ってしまっているので、わかりやすくBETAとします。

BETA=pinv(X)*Y

そうすると、BETAは以下のように表示されます。

BETA =
   -5.0000
    5.0000
   15.0000
   15.0000

今、BETAは上からAの平均値、Bの平均値、Cの平均値、そして定数となります。

なので、AとBの差の絶対値は10となりますし、AとCの差の絶対値は20となります。

それでは、次にAとBの2群をモデルしてみましょう。

今の場合、Yは以下のようになります。

Y=[A;B]
Y =
     9
    10
    11
    19
    20
    21

このときのDesign Matrixは以下のようになります。

X = [kron(eye(2), ones(3,1)) ones(6,1)]
X =
     1     0     1
     1     0     1
     1     0     1
     0     1     1
     0     1     1
     0     1     1

同じようにBETAを求めましょう。

BETA=pinv(X)*Y
BETA =
   -0.0000
   10.0000
   10.0000

今の場合、BETAは上からAの平均値、Bの平均値、定数となります。

Aの平均値とBの平均値の差の絶対値は10となります。上のモデルと同じです。

では、AとCの場合はどうでしょうか。

Y=[A;C]
Y =
     9
    10
    11
    29
    30
    31

Xは同じものが使えますので、BETAを求めます。

BETA=pinv(X)*Y
BETA =
   -3.3333
   16.6667
   13.3333

今の場合、BETAは上からAの平均値、Cの平均値、定数となりますので、
Aの平均値とCの平均値の差の絶対値は(約ではありますが)20となります。

つまり、このような場合、Design Matrixを複数作るよりも、1つのDesign matrixを作ってしまった方がエレガントだということになります。

今までいつも迷っていたことだったので、すっきりしました。

コメントを残す

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください