Mixed Model ~混合モデル~

統計手法アラカルト
Mixed Model
~混合モデル~
行動計量学研究分野
学部三回 兼清道雄
内容
• GLMとは
• Mixed Modelとは
• SASでの使用例
– SPSSではver.11からサポート
2002/10/30
まぜまぜもでる
2
GLMとは
• General Linear Modelの略
• 一般線型モデル
• y を目的変数(従属変数)とし、p個の説明
変数(独立変数)x1,…,xpとの線型結合で
表される式(↓)のこと
2002/10/30
まぜまぜもでる
3
行
列
表
現
し
ま
す
と
・
・
・
2002/10/30
まぜまぜもでる
4
GLMによる解析
• 回帰分析(単回帰、重回帰)、分散分析(多
変量分散分析)、共分散分析など
• SASのGLMプロシージャでは上記の分析
が可能
2002/10/30
まぜまぜもでる
5
GLM風:分散分析
• 1要因(3水準(i=1,2,3))分散分析
• 各水準に2データ(j=1,2)、計6個のデータ
• データの構造式
全部書くと
2002/10/30
被験者1
まぜまぜもでる
被験者2
6
2002/10/30
まぜまぜもでる
7
GLMの欠点
• 誤差共分散の構造を柔軟に指定出来ない
– 誤差共分散:ε(誤差)同士の分散共分散行列
– 測定値間に何らかの関係があると考えて分析
できない
• ちゃんとした計画法を用いれば、まあ、それなりには
できるけれども・・・(分散分析)
– 反復測定や経時観察のデータでは致命的
2002/10/30
まぜまぜもでる
8
そこで・・・
• GLMでは満足に解析できない反復測定
データや経時観察データに対して・・・
• Mixed Modelの出番ってなわけです
• Mixed Modelって??
2002/10/30
まぜまぜもでる
9
Mixed Model
• General Linear Mixed Model
• 一般線型混合モデル
• 線型混合モデル、混合モデルともいう
• GLMの拡張
2002/10/30
まぜまぜもでる
10
どこが拡張??
• 変量効果の導入
– (対義語)固定効果
– GLMは固定効果しか扱えなかった
• 変量効果も固定効果として扱った
• 誤差共分散への様々な指定
– 後ほど(SASのところ)、詳しく説明
2002/10/30
まぜまぜもでる
11
固定効果?
• 固定効果(fixed-effect)
• ある要因を固定効果として解釈
– その要因に対して有限個の水準を想定
– 研究では、評価したい全ての水準を含んでい
ると考えることになる
– 解釈例
• 性別
• トレーニング有無
2002/10/30
まぜまぜもでる
12
変量効果?
• 変量効果(random-effect)=ランダム効果
• ある要因を変量効果として解釈
– その要因に対して無限個の水準(水準の母集
団)を想定
– 研究では、実際に検討する水準は母集団から
の標本であると考えることになる
– 解釈例
• 大学(実験や調査の場所として)や被験者
2002/10/30
まぜまぜもでる
13
固定?変量?具体例1
• 質問紙調査(複数の大学で)
• 因子分析後、尺度得点を使って分散分析
– コンピュータ不安の尺度得点とか
• 『大学』を1つの要因とする
– 他は性別とか
• つまり、大学や性別という要因でコン
ピュータに対する不安を説明しようとする
2002/10/30
まぜまぜもでる
14
固定?変量?具体例2
• 『大学』を・・・
• 固定効果と解釈
– 特定の大学の効果を考える(ex.阪大、神大、京大)
• 変量効果と解釈
– 大学間の全体的なばらつきを考える
• 次スライドを参照
2002/10/30
まぜまぜもでる
15
!図!
阪大は・・
神大は・・
固定効果だと
変量効果だと
35
コ
ン
尺ピ
度ュ
得ー
点タ
不
安
こ
の
ば
ら
つ
き
に
興
味
30
25
20
15
10
5
0
阪大
2002/10/30
神大
まぜまぜもでる
京大
16
本当はね(変量効果として解
釈)
こ
の
ば
ら
つ
き
に
興
味
35
コ
ン
尺ピ
度ュ
得ー
点タ
不
安
30
25
20
15
10
5
0
A大
2002/10/30
B大
C大
まぜまぜもでる
D大
E大
17
つまり
• 固定効果では個々の水準における効果に注目
ex.) ○阪大、神大、京大間に差はあるか?
○阪大と京大どちらがコンピュータ不安が高い
か?
• 変量効果では全体的なばらつきに注目
ex.) ○大学間で効果はどれだけばらつくか?
2002/10/30
まぜまぜもでる
18
変量効果?
• データの構造式ではよく「b」で表されます
– 固定効果の「β」に対して
• 実験毎に違うもの=変量効果
– 被験者(回答者)が毎回同じではない=変量効果
• つまり、確率変数です
– (誤差も確率変数です)
• 正規分布に従います(と仮定)
–あ
2002/10/30
まぜまぜもでる
19
モデルの構造式
• Mixed Model
• ちなみにGLMは
2002/10/30
まぜまぜもでる
20
モデル構造式(実際は)
• Mixed Model
• GLM
というわけで、パラメータの数は同じ(βG=β+b)
2002/10/30
まぜまぜもでる
21
具体例
対応のある1要因分散分析
• 1要因3水準(i=1,2,3)の分散分析
• 被験者を変量効果と解釈する(j=1,2)
• データの構造式は
全部書くと
2002/10/30
まぜまぜもでる
22
2002/10/30
まぜまぜもでる
23
平均、分散
• Mixed Model
• ちなみにGLMは
2002/10/30
まぜまぜもでる
24
以上より
• 測定値間の何らかの関係を、変量効果の
導入や誤差共分散の指定により、的確に
捉えることが出来る
• また、欠測値があっても解析可能
– ただし、MAR(Missing At Random)の場合
2002/10/30
まぜまぜもでる
25
Q.どうやってMixedを使う
の?
A.話はSASに飛びます
SAS
• PROC MIXEDで分析が可能
• 表1のデータをもとに、被験者効果を変量
効果とし、1要因分散分析を行う
• 要因Aは3水準
– 練習前、1週間練習、2週間練習
• データの構造式
2002/10/30
まぜまぜもでる
27
rとlとwの弁別実験(行動計量学講義資料より抜粋)
被験者
表
1
1
2
3
4
5
6
7
8
平均
練習前
0.40
0.42
0.42
0.42
0.47
0.38
0.44
0.40
0.41875
正答率
1週間練習 2週間練習
0.43
0.43
0.45
0.44
0.42
0.43
0.43
0.42
0.49
0.50
0.37
0.39
0.47
0.46
0.41
0.42
0.43375
0.43625
DATA rlw;
DO sub = 1 to 8;
DO a = 1 to 3;
INPUT y @@; OUTPUT;
SASプログラム1
データステップ編
END;END;
CARDS;
.40 .43 .43
.42 .45 .44
.42 .42 .43
.42 .43 .42
.47 .49 .50
.38 .37 .39
.44 .47 .46
.40 .41 .42
PRINTプロシージャ
で確認
;
PROC PRINT DATA=rlw;
RUN;
RUN;
SASプログラム2
PROC MIXED DATA=rlw;
CLASS a sub;
MODEL y = a;
RANDOM intercept/subject=sub;
RUN;
PROC MIXED DATA=rlw;
CLASS a sub;
MODEL y = a;
RANDOM sub;
どちらでも同じデータ
の構造になります
RUN;
2002/10/30
まぜまぜもでる
30
OUTPUT(付録OUT1参照)
• CLASSステートメントで指定した要因の水準
を表示
– Class Level Information
• 推定の為の反復計算の過程
– REML Estimation Iteration History
– 推定する分散成分が小さい場合、収束しないこと
もある
– 今回は収束(convergence criteria met.)
2002/10/30
まぜまぜもでる
31
OUTPUT(付録OUT1参照)
• 分散成分とその推定値
– Covariance Parameter Estimates(REML)
Cov Parm
Estimate
SUB
0.00098155
Residual
0.00007857
– 測定誤差によるばらつきは被験者によるばらつき
よりかなり小さい
• つまりばらつきのほとんどが被験者によるもの
– オプション‘COVTEST’を使えば、標準誤差や検定
等計量が表示される(後述)
2002/10/30
まぜまぜもでる
32
OUTPUT(付録OUT1参照)
• モデルの当てはまりに対する情報
– Model Fitting Information for Y
Description
Value
Observations
24.0000
Res Log Likelihood
53.5487
Akaike's Information Criterion
51.5487
Schwarz's Bayesian Criterion
50.5041
-2 Res Log Likelihood
-107.097
– おもにモデルを比較する時に使用
– 詳細は割愛
2002/10/30
まぜまぜもでる
33
OUTPUT(付録OUT1参照)
• 固定効果の検定
– Tests of Fixed Effects
Source
A
NDF
DDF
Type III F
Pr > F
2
14
9.12
0.0029
– 要因Aの効果あり
• 練習前と1週間練習した後と2週間練習した後では、正
答率が違う
• どことどこに有意な差?
• LSMEANSステートメントで多重比較可能(次スライド)
2002/10/30
まぜまぜもでる
34
PROC MIXED DATA=rlw;
PROC MIXED DATA=rlw;
CLASS a sub;
CLASS a sub;
MODEL y = a;
MODEL y = a;
RANDOM intercept/subject=sub;
RANDOM sub;
LSMEANS a/adjust=tukey;
LSMEANS a/adjust=tukey;
RUN;
RUN;
Least Squares Means
Effect A
LSMEAN
Std Error
DF
t Pr > |t|
A
1
0.41875000
0.01151151 7.74
36.38
0.0001
A
2
0.43375000
0.01151151 7.74
37.68
0.0001
A
3
0.43625000
0.01151151 7.74
37.90
0.0001
Differences of Least Squares Means
Effect A _A
Difference
Std Error
DF
t Pr > |t|
Adjustment
Adj P
A
1 2
-0.01500000
0.00443203
14
-3.38
0.0044
Tukey-Kramer 0.0170
A
1 3
-0.01750000
0.00443203
14
-3.95
0.0015 Tukey-Kramer 0.0068
A
2 3
-0.00250000
0.00443203
14
-0.56
0.5816
Tukey-Kramer 0.8417
• 練習することにより正答率は上がるが、練習の期
間によって正答率が変化するとはいえない
ステートメントやオプションの
紹介
PROC MIXED
• DATA=~
– 分析するデータ名を指定
• COVTEST
– 分散成分の標準誤差および検定統計量出力
– ただし、検定に関しては微妙
• METHOD=~
– 推定方法を指定 ex.METHOD=ML (最尤法)
– デフォルトはREML(制限付最尤法)
2002/10/30
まぜまぜもでる
37
PROC MIXED
• NOCLPRINT
– Class Level Informationを非表示
• NOITPRINT
– Estimation Iteration Historyを非表示
PROC MIXED DATA=rlw METHOD=ML COVTEST NOCLPRINT NOITPRINT;
… … …
RUN;
2002/10/30
まぜまぜもでる
38
CLASSステートメント
• 因子として考えるべき変数を指定
• 文字変数でも数値変数でもよい
2002/10/30
まぜまぜもでる
39
MODELステートメント
• 反応変数(1変数のみ)と固定効果を指定
– y = a とか y = a b a*b とか
• 切片は自動的に含まれる
– nointオプションで切片なしを指定出来る
• s または solution
– 固定効果の推定を表示
– 点推定、標準誤差、t統計量、p値
• ddfm=~
– 自由度の求め方を指定:デフォルトは以下
– ddfm=betwithin(REPEATEDステートメントのみ)
– ddfm=contain(RANDOMステートメント含む)
2002/10/30
まぜまぜもでる
40
MODELステートメント
Solution for Fixed Effects
Effect
A
INTERCEPT
Estimate
Std Error
DF
t
Pr > |t|
0.43625000
0.01151151
7
37.90
0.0001
A
1
-0.01750000
0.00443203
14
-3.95
0.0015
A
2
-0.00250000
0.00443203
14
-0.56
0.5816
A
3
0.00000000
.
.
.
.
• 前例:MODEL y=a/sの場合のアウトプット
• A3水準=切片となっている
2002/10/30
まぜまぜもでる
41
RANDOMステートメント
• 変量効果を指定
• 切片を入れるためにはinterceptを変数に
• subject=
– データセットにおける対象者を識別
• type=
– 変量効果の共分散行列の構造を指定
– 実用的にはUNかVC(UN、VCについては後述)
• g
– 変量効果の共分散行列の推定値を出力
• gcorr
– 変量効果の相関行列の推定値を出力
2002/10/30
まぜまぜもでる
42
RANDOMステートメント
• データの構造式
PROC MIXED DATA=rlw;
PROC MIXED DATA=rlw;
・・・・・・
・・・・・・
RANDOM intercept/SUBJECT=sub;
RANDOM sub;
RUN;
RUN;
43
REPEATEDステートメント
• 誤差共分散行列(共分散の構造)を指定
• type=
– 詳しくは次スライド
• subject=
– データセットにおける対象者を識別
• 反復効果を示す変数は名義変数
– 入れなくてもよい
2002/10/30
まぜまぜもでる
44
type=
• type=AR(1)
– 一次自己回帰
– 近いものには強い
関係
– 遠いものには弱い
関係
– 経時データ向き
2002/10/30
まぜまぜもでる
45
type=
• type=CS
– 複合対称性(Compound Symmetry)
– 測定値のばらつきが一定(定数分散)
– 測定値間の関係も一定(定数共分散)
2002/10/30
まぜまぜもでる
46
type=
• type=SIMPLE
• type=VC
– Variance
Components
– 測定値のばらつき
が一定(定数分
散)
– 測定値間の関係
なし(無相関)
2002/10/30
まぜまぜもでる
47
• 変量効果が複数個
• RANDOMステートメント
に‘type=SIMPLE(VC)’
で指定
– つまり変量効果同士の共
分散構造をtypeで指定し
た場合
• 変量効果ごとのばらつき
が一定ではない
(不定数分散)
• 変量効果間の関係なし
(無相関)
まぜまぜもでる
48
type=
• type=UN
– 無構造
– Unstructured
2002/10/30
まぜまぜもでる
49
その他にも
• たくさんの構造が指定できます
• 誤差についても
– 系列相関成分
– 測定誤差成分
– と分けることも出来ます
• 今回は割愛します
2002/10/30
まぜまぜもでる
50
LSMEANSステートメント
• 要因の各水準の平均値の推定
• 水準間の平均値の差の検定
– LSMEANS a / adjust=tukey;
– 要因Aの水準間の多重比較(tukey法を用いて)
• 単純主効果の検定
– LSMEANS a / slice = b;
– bの水準ごとの要因Aの効果を検定
2002/10/30
まぜまぜもでる
51
CONTRASTステートメント
• 固定効果パラメータの線型結合による仮説の
検定を行う
–あ
• CONTRASTステートメントでLを指定
– 次スライドを参照
• 様々な仮説に対して柔軟に検定
2002/10/30
まぜまぜもでる
52
具体例(前述:rlwの弁別実験、要因A、3水準)
intで指定
aで指定
→ int 1 a 1 0 0
→ a 1 -1 0
2002/10/30
まぜまぜもでる
53
CONTRASTステートメント
• 練習前と1週間練習後の平均に差がある
か検定(
) ‘ラベル名(20字以
内)’
• CONTRAST ‘before_1week’ a 1 –1
0;
CONTRAST Statement Results
Source
before_1week
NDF
DDF
F
Pr > F
1
14
11.45
0.0044
• 5%有意、差があると言える
– LSMEANSとの違いは多重性によるもの
2002/10/30
まぜまぜもでる
54
CONTRASTステートメント
• 練習前と2週間練習後の母平均の平均は、1週
間練習後の母平均と差があるかの検定(も可能)
• あ
• CONTRAST ‘(be+2w)/2=1w’ a 0.5 –1 0.5;
CONTRAST Statement Results
Source
(be+2w)/2=1w
NDF
DDF
F
Pr > F
1
14
2.65
0.1257
• 有意差無し
– 練習効果は非直線的ではないことがわかる
2002/10/30
まぜまぜもでる
55
CONTRASTステートメント
• あ
• CONTRAST ‘allmeans=0’ int 1 a 1 0 0,
int 1 a 0 1 0,
int 1 a 0 0 1;
CONTRAST Statement Results
Source
allmeans=0
2002/10/30
NDF
DDF
3
14
まぜまぜもでる
F Pr > F
494.41
0.0001
56
ESTIMATEステートメント
• 固定効果パラメータの線型結合による値
の推定を行う
–あ
• Lの指定の仕方はCONTRASTと同じ
– ラベルをつけるところも同じ
• 区間推定も可能
– ESTIMATE ‘~’ ~ / cl alpha=0.05;
2002/10/30
まぜまぜもでる
57
ESTIMATEステートメント
• 練習前の平均の推定(μ1)
• ESTIMATE ‘before’ int 1 a 1 0 0 / cl
alpha=0.05;
ESTIMATE Statement Results
Parameter
before
Estimate
Std Error
DF
t
0.41875000
0.01151151
14
36.38
2002/10/30
まぜまぜもでる
Pr > |t| Alpha
0.0001
0.05
Lower
Upper
0.3941
0.4434
58
レポート課題
1. 変量効果と考えられる要因を一つ挙げ、そ
の理由を述べてください(資料中のものも可)
2. 付録の表2のデータをMIXEDプロシージャを
用いて分析、解釈してください
•
プログラム、及び主要なアウトプットを添付のこと
3. 今回の講義でわからなかったところを書きた
いだけ書いてください
•
詳しく書いていただけるとなお嬉しい
〆切:11月13日(水) 提出先:狩野助教授室(北館304)
2002/10/30
まぜまぜもでる
59
参考文献
• 医学統計のための線型混合モデル-SAS
によるアプローチ-
– 松山 裕・山口拓洋(編訳) 2001
– サイエンティスト社
• 資料(なれよう、過去の発表)
2002/10/30
まぜまぜもでる
60
☆告知☆
• Mixed Modelになれよう
–
–
–
–
水曜4限 北館301でやってます
現在のお題は欠測データとMixed Modelです
以下のページに今までの資料などがあります
http://koko15.hus.osaka-u.ac.jp/~kanekiyo/mixed/
2002/10/30
まぜまぜもでる
61
なぜMixed Modelというのか?
• 固定効果(fixed-effect)
+変量効果(random-effect)
まぜまぜまぜまぜまぜまぜ
• 混合効果(mixed-effect)
• よってMixed Modelとなったわけです
2002/10/30
まぜまぜもでる
62