2012年12月18日火曜日

分散分析から回帰分析へ:回帰式による表現を理解する

 心理言語学を志す日本の学生の多くは、言語学科か心理学科に属しているかと思いますが、両者の違いの一つとして、次のことばに対する認識への違いを指摘できます。

「分散分析も所詮は回帰分析のひとつだよね。」

 統計学を必修科目としてお付き合いしてきた心理学出身の人がこれを聞くとたいてい
「そうだね。」「詳しくは知らないけど、授業で聞いたことはある。」という反応ですが、言語学出身の人のなかには
「そうなの?」「そんなの聞いたことないよ。」
という反応のひとがいるかもしれません。

今日は分散分析と回帰分析が同等であることを確認してみます。

1.条件の違いを示す棒グラフを線グラフに変更する。
以下は条件Aと条件Bについて何らかの成績の違いを示しているグラフです。一般的に条件間差はt検定で検定できます。













このグラフを線グラフに直してみます。また、「条件A」や「条件B」はラベルなので何と呼んでもいいはずです。ここで「条件A」を「0」、「条件B」を「1」と呼び変えます。


2.縦軸を移動する。
グラフ上の成績を示す縦軸は、横軸のどこについていてもいいはずです。ここで条件0の上に縦軸が引かれるように、縦軸を移動します。
















上の形を見てみると、中学校の数学で習った一次方程式を示す直線の形になっていることに気づきます。直線のy軸との交点は切片(intercept)、x軸が1変化するときのy軸の変化量を傾き(slope)と呼んでいました。ここで最初の二つの条件間差が傾きに一致していることに注意してください。

つまり、二つの条件はどちらも次のような方程式で書くことができます。

成績 = 切片 + 傾き x (条件の水準)

条件Aは条件の水準0なので
成績 = 切片 + 傾き x (0) = 切片

条件Bは条件の水準1なので
成績 = 切片 + 傾き x (1) = 切片 + 傾き

方程式上でも条件間差は傾きの大きさに一致しています。
こんな風に見方を少し変えるだけで、t検定で扱っていた条件間差が、方程式として表現できることが分かります。
 この考え方を理解できると「分散分析も回帰モデルの一種」という冒頭の言葉の意味がわかるのではないでしょうか。



2012年12月17日月曜日

R言語を用いた混合モデル解析の紹介その1:1要因デザインの場合(8)


既に出版されている論文を見てみると三つの対処法があるようです。

A:t値の絶対値が2より大きいときは十分に効果が大きいとみなす。
B:マルコフ連鎖モンテカルロシミュレーション(MCMCと呼ばれます)を実施して傾きの信頼区間を求めて有意性を調べる。
C:固定因子の傾きのあるモデルとないモデルによる推定を行い、二つのモデルを適合度で比較する。

 ではAの対処方法をとる場合ですが、今回の例ではt値が43.7と非常に大きいので、これを主張することができそうです。

Bの対処法ですが、languageRというパッケージに含まれるpvals.fnc()というコマンドを使うと
MCMCによる推定ができます。今回の例でMCMCを実施すると

MCMC(model.3)


$fixed
                  Estimate MCMCmean HPD95lower   HPD95upper  pMCMC Pr(>|t|)
(Intercept)       400.29   400.35              394.32           406.69     0.0001        0
CondDummy      48.15    48.14                45.14             51.25     0.0001        0

$random
    Groups        Name   Std.Dev. MCMCmedian MCMCmean HPD95lower HPD95upper
1  Subject (Intercept)    20.8991     8.4280               8.4860     6.8950    10.1637
2     Item (Intercept)       8.0714     5.9127               6.0536     3.8255     8.5196
3 Residual                     6.6118     9.3065               9.3425     7.9304    10.7475

上記のような出力結果が得られました。四行目がCondDummyの傾きの推定結果です。
95%信頼区間が45.14~51.25でした。pMCMCがいわゆるp値で0.0001となりました。

一点注意が必要なのは、languageRを利用したMCMCは一部のモデルではできないことです。
できないのは切片と傾きの両方を仮定し、かつそれらに共分散(相関)を仮定している場合です。
(上記のmodel.1とmodel.2がこれに該当します。)


Cの対処法をとる場合、CondDummyの固定因子を含まないモデルで推定してみます。

model.null = lmer(RT ~ (1 | Subject) + ( 1| Item))


次にこのモデルと先ほどのmodel.3を比較します。

anova(model.3, model.null)


Data: 
Models:
model.null: RT ~ (1 | Subject) + (1 | Item)
model.3: RT ~ CondDummy + (1 | Subject) + (1 | Item)
           Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)    
model.null  4 1390.7 1402.6 -691.35                             
model.3     5 1050.6 1065.5 -520.31 342.08      1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


出力結果を見ると、χ二乗値は342.08、自由度は1、p値は< 2.2e-16***と非常に小さくなることが示されています。

今回の例では、いずれの方法をとっても条件差が有意であると言えます。
解析は以上になります。


R言語を用いた混合モデル解析の紹介その1:1要因デザインの場合(7)


もう一度model.3の出力結果を見ながら、論文に載せるべき情報を整理します。


Linear mixed model fit by REML 
Formula: RT ~ CondDummy + (1 | Subject) + (1 | Item) 
  AIC  BIC logLik deviance REMLdev
 1043 1058 -516.5     1041    1033
Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 436.771  20.8991 
 Item     (Intercept)  65.147   8.0714 
 Residual              43.716   6.6118 
Number of obs: 144, groups: Subject, 12; Item, 12

Fixed effects:
            Estimate Std. Error t value
(Intercept)  400.292      6.514   61.45
CondDummy     48.153      1.102   43.70

Correlation of Fixed Effects:
          (Intr)
CondDummy -0.085

一行目は推定方法が述べられています。REMLは制限最尤法のことです。
二行目はモデルの回帰方程式です。RTという従属変数に対して、CondDummyという固定因子と切片の被験者ランダム因子、切片の項目ランダム因子がモデルに組み込まれていることを示しています。
十一行目からが固定因子の効果量を示しています。ここで重要なのはCondDummyの傾きで推定値が48.153、標準誤差が1.102、t値が43.70となっています。

さて、ここで一つ重要な問題が起こります。固定因子の効果量についてt値が報告されているのにp値が報告されていないんです。これでは有意がどうかを論じることができません。

(8)へ続く

R言語を用いた混合モデル解析の紹介その1:1要因デザインの場合(6)


ではさらにモデルを単純にできるかどうかを試してみます。以下のmodel.3では、model.2からさらに被験者のランダム因子のスロープを除いています。

model.3 = lmer(RT ~ CondDummy + ( 1 | Subject) + (1 + Item))

  出力は以下のようになります。
summary(model.3)

Linear mixed model fit by REML 
Formula: RT ~ CondDummy + (1 | Subject) + (1 | Item) 
  AIC  BIC logLik deviance REMLdev
 1043 1058 -516.5     1041    1033
Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 436.771  20.8991 
 Item     (Intercept)  65.147   8.0714 
 Residual              43.716   6.6118 
Number of obs: 144, groups: Subject, 12; Item, 12

Fixed effects:
            Estimate Std. Error t value
(Intercept)  400.292      6.514   61.45
CondDummy     48.153      1.102   43.70

Correlation of Fixed Effects:
          (Intr)
CondDummy -0.085

さらにmodel.2とmodel.3の適合度をanova()で検証してみます。
anova(model.2, model.3)

Data: 
Models:
model.3: RT ~ CondDummy + (1 | Subject) + (1 | Item)
model.2: RT ~ CondDummy + (1 + CondDummy | Subject) + (1 | Item)
        Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
model.3  5 1050.6 1065.5 -520.31                         
model.2  7 1054.4 1075.2 -520.22 0.1854      2     0.9115



出力結果を見るとChisqが0.1854、自由度(Df)が2 Pr(>Chisq)が0.9115です。つまりmodel.2はmodel.3よりも有意に適合度が高いとは言えないことが示されました。
この場合は、model.3を採用します。そしてmodel.3が一番単純なモデルなので、これを最終モデルとします。
(7)へ続く

R言語を用いた混合モデル解析の紹介その1:1要因デザインの場合(5)


2-2. 簡略モデルによる推定と最大モデルとの比較

 混合モデル解析では、ランダム因子の構造について複数の可能性があります。そこでデータに対する適合度を指標として「適合度が最大」かつ「もっとも単純な」構造のモデルを最終モデルとして選択します。まずは上記の最大モデルから項目のランダム因子のスロープを取り除いたモデルで推定します。

model.2 = lmer(RT ~ CondDummy + ( 1 + CondDummy | Subject) + (1 | Item))

このモデルの推定結果は以下のようになります。
summary(model.2)


Linear mixed model fit by REML 
Formula: RT ~ CondDummy + (1 + CondDummy | Subject) + (1 | Item) 
  AIC  BIC logLik deviance REMLdev
 1047 1067 -516.3     1040    1033
Random effects:
 Groups   Name        Variance Std.Dev. Corr  
 Subject  (Intercept) 429.6527 20.7281        
          CondDummy     4.0436  2.0109  0.149 
 Item     (Intercept)    66.2690  8.1406        
 Residual                   42.6313  6.5293        
Number of obs: 144, groups: Subject, 12; Item, 12

Fixed effects:
                   Estimate Std. Error t value
(Intercept)      400.292      6.474   61.83
CondDummy     48.153      1.233   39.04

Correlation of Fixed Effects:
                    (Intr)
CondDummy -0.009


 次にmodel.maxとmodel.2を比較します。モデルの単純さという点では、model.2がランダム因子のスロープを除いた分(パラメータ数では2減っています)だけ優っています。あとは適合度を比べてmodel.2とmodel.maxが同程度の適合度なのか、それともmodel.maxの方が適当度が有意に高いのかを検証します。検証にはanova()というコマンドを用います。

anova(model.max, model.2)

Models:
model.2: RT ~ CondDummy + (1 + CondDummy | Subject) + (1 | Item)
model.max: RT ~ CondDummy + (1 + CondDummy | Subject) + (1 + CondDummy | 
model.max:     Item)
              Df    AIC    BIC  logLik      Chisq Chi Df Pr(>Chisq)  
model.2      7 1054.4 1075.2 -520.22                           
model.max  9 1053.7 1080.4 -517.84  4.7645      2    0.09234 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

 出力結果の下から3行目を見るとChisqが4.7645、自由度(Df)が2 Pr(>Chisq)が0.09234とあります。Chisqは各モデルの対数尤度の差の2倍の値です。自由度はパラメータ数の差になっています。P値はこの二つの値から求められており、model.maxはmodel.2よりも有意に適合度が高いとは
言えないことが示されました。ということで暫定的にmodel.2を採用します。

(6)へ続く

R言語を用いた混合モデル解析の紹介その1:1要因デザインの場合(4)


2-1. 最大モデルによる推定

それではモデルを構築して推定してみましょう。まずは混合モデル解析のためのパッケージlme4を起動します。

library(lme4)

  次にrawdataの中の変数を変数名で直接呼び出すためのコマンドを実行します。
attach(rawdata)

  準備ができたので、モデルを構築します。固定因子はCondDummyでSubjectとItemをランダム因子とします。ランダム因子の構造は切片と被験者内(項目内)要因のスロープすべてを仮定します。推定した結果をmodel.maxという変数に書きこむことにします。

model.max = lmer(RT ~ CondDummy + ( 1 + CondDummy | Subject) + (1 + CondDummy | Item))

  上の式を少しずつ説明します。
 混合モデル推定を実行するためのコマンドはlmer()です。()の中にモデルを書きこみます。
 モデルは 従属変数 ~ 固定因子 + (ランダム因子)のように書きます。 ~の左側に従属変数、
~の右側に説明変数を書きます。()でくくったものがランダム因子、くくらないものが固定因子と
認識されます。
 ランダム因子を示す()の中は|でしきられています。|の左側に設定する項を書き、|の右側にランダム変数を書きます。|の左側の数字1は切片を示しています。つまり
(1+ CondDummy | Subject)は被験者というランダム因子の変動を切片とCondDummyという因子のスロープの二つに仮定することを意味します。

結果を見るには先ほど使ったsummary()コマンドなどを使います。
summary(model.max)


Linear mixed model fit by REML 
Formula: RT ~ CondDummy + (1 + CondDummy | Subject) + (1 + CondDummy |      Item) 
  AIC  BIC logLik deviance REMLdev
 1045 1072 -513.6     1036    1027
Random effects:
 Groups   Name        Variance Std.Dev. Corr  
 Subject  (Intercept)  431.5118 20.7729        
          CondDummy      5.7161  2.3908  0.115 
 Item     (Intercept)     60.1954  7.7586        
          CondDummy    17.7028  4.2075  0.076 
 Residual                   37.5728  6.1297        
Number of obs: 144, groups: Subject, 12; Item, 12

Fixed effects:
                    Estimate Std. Error t value
(Intercept)      400.292      6.442   62.14
CondDummy     48.153      1.731   27.82

Correlation of Fixed Effects:
                  (Intr)
CondDummy 0.015 

出力結果の読みとりやどの情報を論文に載せるべきなのかは後述します。
(5)へ続く

R言語を用いた混合モデル解析の紹介その1:1要因デザインの場合(3)


1-2. ダミーコーディングを使って条件を符号化

データを見ていただいてわかるように、今のデータでは条件(Condition)がaまたはbとなっています。次はこれを数値へと符合化します。符合化の方法は複数ありますが、今回はダミーコーディングを使ってa条件を0、b条件を1とします。この数値を書き込むための別の変数(CondDuumy)を用意します。

まずは0をデフォルトの値として変数を設定します。
rawdata$CondDummy = 0

次にConditionの変数がbのケースの場合にCondDummyの値が1になるようにします。
rawdata[rawdata$Condition == "b","CondDummy"] = 1

確認としてもう一度、rawdataの要約を見てみましょう
summary(rawdata)



 Subject           Item             Condition       RT                  listnumber  
 Min.   : 1.00      Min.   : 1.00     a:72             Min.   :355.2       Min.   :1.00  
 1st Qu.: 3.75     1st Qu.: 3.75   b:72             1st Qu.:400.6      1st Qu.:1.75  
 Median : 6.50    Median : 6.50                     Median :423.1     Median :2.50  
 Mean   : 6.50    Mean   : 6.50                      Mean   :424.4     Mean   :2.50  
 3rd Qu.: 9.25     3rd Qu.: 9.25                     3rd Qu.:445.6      3rd Qu.:3.25  
 Max.   :12.00     Max.   :12.00                      Max.   :495.0      Max.   :4.00  

CondDummy  
 Min.   :0.0  
 1st Qu.:0.0  
 Median :0.5  
 Mean   :0.5  
 3rd Qu.:1.0  
 Max.   :1.0 


CondDummyという変数が追加できました。
(4)へ続く

R言語を用いた混合モデル解析の紹介その1:1要因デザインの場合(2)


1-1. データの読み込み

では、(1)のデータをRの中に読みこんでみましょう。(Rを起動してください)
一番簡単なのは、クリップボードを利用する方法です。
上のデータの部分(変数名ごと)ドラッグ&コピーをします。
コピーした状態で

rawdata = read.table("clipboard", header = TRUE)

とコマンドを書きこみます。
Mac版を使っている場合は次のように書きます。

rawdata = read.table(pipe("pdpaste"), header = TRUE)

これでrawdataというオブジェクトの中にすべての変数が名前付きで読み込まれました。

データが正常に読み込まれたことを確認してみましょう。
データの要約を見るにはsummaryというコマンドを使用します。

summary(rawdata)

結果は次のようになります。
  Subject    Item       Condition     RT       listnumber
Min. : 1.00    Min. : 1.00    a:72    Min. :355.2    Min. :1.00
1st Qu.: 3.75  1st Qu.: 3.75   b:72     1st Qu.:400.6  1st Qu.:1.75
Median : 6.50  Median : 6.50         Median :423.1  Median :2.50
 Mean : 6.50   Mean : 6.50          Mean :424.4   Mean :2.50
 3rd Qu.: 9.25  3rd Qu.: 9.25         3rd Qu.:445.6  3rd Qu.:3.25
 Max. :12.00   Max. :12.00          Max. :495.0    Max. :4.00


こういったコマンドを利用してデータがきちんと読み込まれていることを確認していくことは意外に重要です。

(3)へ続く

R言語を用いた混合モデル解析の紹介その1:1要因デザインの場合(1)


本日はR言語で仮想データを実際に解析しながら、混合モデル解析の手続きを紹介します。
第一弾の今回は、もっとも単純な一要因のデザインを考えます。

仮想データは、次のような実験デザインを想定しています。

調べたい要因(Condition) 条件a vs. 条件b
従属変数 反応時間(RT)
被験者(Subject) 12名
刺激セット数(Item) 12セット
要因配置(listnumber) ラテン方格デザイン(4リストを作成)

データは以下の通りです。また以下では解説を黒字、R内でのコマンドを赤字、出力結果を青字で示します。

SubjectItemConditionRTlistnumber
11a363.51
12a362.751
13a3601
14a355.251
15a357.51
16a358.751
17b412.251
18b420.51
19b421.751
110b4231
111b434.751
112b4241
21b410.52
22b400.252
23b4072
24b408.252
25b420.52
26b429.252
27a386.252
28a378.52
29a382.752
210a3892
211a376.252
212a391.52
31a374.53
32b412.253
33a3643
34b435.253
35a366.53
36b423.253
37a376.253
38b448.53
39a394.753
310b4413
311a392.253
312b4443
41b422.54
42a367.754
43b4194
44a385.254
45b446.54
46a378.754
47b445.254
48a383.54
49b443.754
410a4004
411b444.754
412a392.54
51a382.51
52a386.751
53a3811
54a379.251
55a391.51
56a380.751
57b436.251
58b443.51
59b435.751
510b4351
511b440.751
512b4451
61b443.52
62b428.252
63b4342
64b442.252
65b444.52
66b444.252
67a403.252
68a397.52
69a411.752
610a4052
611a413.252
612a410.52
71a396.53
72b449.253
73a4053
74b450.253
75a391.53
76b443.253
77a412.253
78b455.53
79a409.753
710b4703
711a414.253
712b4623
81b454.54
82a400.754
83b4364
84a398.254
85b448.54
86a409.754
87b448.254
88a420.54
89b466.754
810a4154
811b475.754
812a422.54
91a396.51
92a409.751
93a4011
94a415.251
95a401.51
96a414.751
97b455.251
98b470.51
99b458.751
910b4651
911b463.751
912b4811
101b460.52
102b464.252
103b4582
104b465.252
105b479.52
106b461.252
107a414.252
108a432.52
109a434.752
1010a4252
1011a428.252
1012a425.52
111a418.53
112b465.253
113a4213
114b469.253
115a416.53
116b478.253
117a436.253
118b474.53
119a436.753
1110b4773
1111a436.253
1112b4953
121b479.54
122a424.754
123b4594
124a421.254
125b481.54
126a417.754
127b482.254
128a435.54
129b480.754
1210a4394
1211b491.754
1212a441.54



記述が長くなりそうなので、いくつかに分割して掲載することにします。
(2)へ続く

ワークショップへご参加いただいた皆様

 12月15日、仙台国際センターにおいてワークショップが開催されました。
当日は、企画者側の予想をはるかに超える多数の方にご参加いただき、
ありがとうございました。 また、配布資料が足りなかったことについて、ご迷惑おかけしましたことをお詫びします。

 企画者としての感想ですが、国内研究者がこの問題へ高い興味・関心を抱いていることを実感しました。学会中、何人かの先生とお話ししましたが、統計手法の勉強は「今度時間ができたときに」と思いつつもなかなかその時間が取れないというご意見をいただきました。学会中はある意味、日常の雑務から離れている時間かと思いますので、そういう意味で認知科学会大会のイベントとして企画できたことはよかったのかもしれないと思います。その一方で、もっと時間があれば丁寧に説明して、活発な議論に導けたのではないかという反省もあります。

 当日、配布資料を希望された方には順次送付していきますので、よろしくお願いいたします。感想やコメントがありましたら、mixed.model.jp@gmail.com
までお願いいたします。

 ワークショップで取り上げた話題の一部について、今後当ブログでもご紹介していく予定です。

2012年12月13日木曜日

論文誌でのLME使用率(3)

文理解に関する研究でどの程度混合モデル解析が使用されているかを把握するため、論文誌における混合モデル解析の使用率を計算しています。

 第3弾ではCognitionを取り上げます。計算の対象となったのは2011年、2012年の電子版が発行されている論文(全部で310本もありました)です。

 このうち、線形混合モデルを使用している論文が35本ありました。使用率はわずか11%です。これまでに調べた3つの論文誌の中で最低となりました。

 つぎに研究内容(文理解研究かそれ以外か)で分類してみました。その結果、文理解研究の論文がが36本、それ以外の研究領域(さまざまなトピックが含まれています)の論文が274本ありました。

 この分類別の線形混合モデルの使用率を計算してみました。すると、文理解研究の論文のうち線形混合モデルを使用していは20本、使用率56%でした。また、他の研究領域では論文数が15本、使用率5.5%でした。線形混合モデルの使用が研究領域によってばらつくという傾向が顕著に見られました。
 

2012年12月10日月曜日

Rを用いたデータ解析ミニチュートリアル応募締め切りました

認知科学会大会二日目にR言語を使用した混合モデル解析のミニチュートリアルを行います。
参加希望者数が当初の予定に達しましたので、12月9日に募集を締め切りました。
本日以降、メールをいただきましてもお受けできません、どうぞ悪しからずご了承ください。

2012年12月3日月曜日

論文誌でのLME利用率(2)


文理解に関する研究でどの程度混合モデル解析が使用されているかを把握するため、論文誌における混合モデル解析の使用率を計算しています。

 第2弾ではLanguage and Cognitive Processesを取り上げます。計算の対象となったのは2011年、2012年の電子版が発行されている論文(全部で127本ありました)です。

 このうち、線形混合モデルを使用している論文が27本ありました。使用率は21%です。この値は前回調べたJournal of Memory and Language(29%)に比べて若干低いです。

 つぎに研究内容(文理解研究かそれ以外か)で分類してみました。その結果、文理解研究の論文がが49本、それ以外の研究領域(単語レベルの処理の研究や発話についての音響、音韻的な要因の研究が多く見られました)の論文が78本ありました。

 分類別の線形混合モデルの使用率を計算してみました。すると、文理解研究の論文のうち線形混合モデルを使用していは19本、使用率39%でした。また、他の研究領域では論文数が8本、使用率10%でした。線形混合モデルの使用が研究領域によってばらつくという傾向は、前回のJournal of Memory and Languageの場合と同様です。そして文理解研究は比較的使用率が高いという傾向も見られました。

 さらに、今回は、年度による違いを見ることができました。2011年の文理解研究では混合モデルの使用率が30%だったのですが、2012年ではこれが45%まで上昇しました。この上昇を考慮すると、Journal of Memory and Languageと比較してこの論文誌の使用率が低いとは単純に言えず、今後同程度にまで上昇するかもしれません。