合同会社DILIGENCEでは一緒に働く仲間を募集しています
反復測定二元配置分散分析をR（carパッケージ）でやってみる（S. Kubota）
on 2017/08/31
例えば，業務ソフトAとBを３カ月ずつ使用し，どちらが業務効率が良かったかを検証したいとする．
つまり，業務ソフト違いによる業務効率の影響と使い続けることでの時系列での業務効率の変化を検討する．このとき，業務ソフトの違い，時系列の変化をそれぞれ要因という．また，この検証を同一の社員で，AとBのソフトをそれぞれ３か月ずつ使い検証したこととする．
この場合は，同一の社員で2つの要因について検証を行ったことになるので，反復測定二元配置分散分析というのを行う．
今回はRとcarパッケージを使ってみる．
データの並びは以下の通り．エクセルにあるデータを読みこんでいく前提で記述する．データは行に被験者，列に要因を分けて（各要因の水準ごと）作成する．employeeは社員のことで行ごとに社員（a～wの計23名）の業務効率（％）が記述されている．また，A_1st_monthの意味は，業務ソフトAを1か月使用した時の業務効率を示している．なお，業務効率というのは，勝手に作った指標であるため，実際はそれぞれの業務評価によって数値は変わるし，要因についても変わってくるはず．
このデータをコピーしてRに読み込ませる．datという名前でデータを読み込み．
> dat<-read.delim("clipboard",header = TRUE)
> dat
employee A_1st_month A_2_nd_month A_3rd_month B_1st_month B_2_nd_month B_3rd_month
1 a 30 43 50 42 53 67
2 b 38 46 51 50 57 65
3 c 37 44 52 43 59 60
4 d 36 43 59 50 50 70
5 e 33 42 57 45 57 65
6 f 39 41 58 45 51 69
7 g 37 49 57 50 50 62
8 h 39 46 56 44 53 68
9 i 37 45 58 47 58 63
10 j 39 47 59 48 51 63
11 k 33 44 55 40 55 66
12 l 36 40 59 43 51 68
13 m 37 44 53 48 51 70
14 n 39 48 54 41 53 61
15 o 35 46 58 43 56 64
16 p 31 42 59 40 55 68
17 q 34 41 57 50 54 69
18 r 34 50 55 47 60 63
19 s 40 45 52 44 50 70
20 t 40 42 60 40 60 61
21 u 32 42 59 48 51 67
22 v 33 50 53 47 53 67
23 w 38 48 52 50 54 68
> さらにデータを解析用に整え直す．datBindという名前の解析用データに変身させる．
> datBind<- cbind(dat$A_1st_month, dat$A_2_nd_month, dat$A_3rd_month, dat$B_1st_month, dat$B_2_nd_month, dat$B_3rd_month)
> datBind
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 30 43 50 42 53 67
[2,] 38 46 51 50 57 65
[3,] 37 44 52 43 59 60
[4,] 36 43 59 50 50 70
[5,] 33 42 57 45 57 65
[6,] 39 41 58 45 51 69
[7,] 37 49 57 50 50 62
[8,] 39 46 56 44 53 68
[9,] 37 45 58 47 58 63
[10,] 39 47 59 48 51 63
[11,] 33 44 55 40 55 66
[12,] 36 40 59 43 51 68
[13,] 37 44 53 48 51 70
[14,] 39 48 54 41 53 61
[15,] 35 46 58 43 56 64
[16,] 31 42 59 40 55 68
[17,] 34 41 57 50 54 69
[18,] 34 50 55 47 60 63
[19,] 40 45 52 44 50 70
[20,] 40 42 60 40 60 61
[21,] 32 42 59 48 51 67
[22,] 33 50 53 47 53 67
[23,] 38 48 52 50 54 68次に要因の組み合わせをつくる．factor1を業務ソフトAとB，factor2を時系列の流れとした．ようするに分散分析をするためのおまじない．idata_factでfactor1とfactor2を組み合わせ，準備の２つ目が完了．
> fact1<-factor(c("A","A","A","B","B","B"))
> fact2<-factor(c("1st","2nd","3rd","1st","2nd","3rd"))
> idata_fact<-data.frame(factor1=fact1,factor2=fact2)
> idata_fact
factor1 factor2
1 A 1st
2 A 2nd
3 A 3rd
4 B 1st
5 B 2nd
6 B 3rd解析用のモデルを作る．Model1という名前でモデルを作成．
> Model<- lm(datBind ~ 1)carパッケージを読み込む．なお，carパッケージがインストールされていなければ読み込めないので注意．パッケージのインストールは他のサイトを参考にしてほしい．
> library(car)でもって解析．Model_analysisという名前で分散分析を実行．summary関数を使って詳細を表示．いろいろな数値が出てきたが，Multivariate Tests: factor1とMultivariate Tests: factor2のPillaiの Pr(>F)を見ればいい．0.05よりも小さければ（未満）影響があるということになる．今回は，factor1のソフトAとB，factor2の時系列の影響ともに有るということが分かった．つまり業務ソフトによって業務効率は変わるし（今回はBの方が良い），両方のソフトとも使うほど業務効率が上がっていくということがわかった．
なお，Multivariate Tests: factor1:factor2のpillaiの結果も0.05を下回る場合にはちょっと厄介．これはそれぞれの要因が絡み合って影響しあっているってこと．たとえば，Aのソフトは使い込んでも業務効率が変わらないけど，Bのソフトは使い込めば業務効率が上がるなんていう可能性があったりするってこと．これらの要因ごとの絡み合った影響を検証するためには，されに解析を進めないといけない．この話はまた今度．
> Model_analysis <- Anova(Model,idata=idata_fact,idesign =~factor1*factor2)
Note: model has only an intercept; equivalent type-III tests substituted.
> summary(Model_analysis)
Type III Repeated Measures MANOVA Tests:
------------------------------------------
Term: (Intercept)
Response transformation matrix:
(Intercept)
[1,] 1
[2,] 1
[3,] 1
[4,] 1
[5,] 1
[6,] 1
Sum of squares and products for the hypothesis:
(Intercept)
(Intercept) 2093466
Multivariate Tests: (Intercept)
Df test stat approx F num Df den Df Pr(>F)
Pillai 1 0.9996 56243.7 1 22 < 2.22e-16 ***
Wilks 1 0.0004 56243.7 1 22 < 2.22e-16 ***
Hotelling-Lawley 1 2556.5319 56243.7 1 22 < 2.22e-16 ***
Roy 1 2556.5319 56243.7 1 22 < 2.22e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
------------------------------------------
Term: factor1
Response transformation matrix:
factor11
[1,] 1
[2,] 1
[3,] 1
[4,] -1
[5,] -1
[6,] -1
Sum of squares and products for the hypothesis:
factor11
factor11 19111.7
Multivariate Tests: factor1
Df test stat approx F num Df den Df Pr(>F)
Pillai 1 0.946827 391.741 1 22 1.6569e-15 ***
Wilks 1 0.053173 391.741 1 22 1.6569e-15 ***
Hotelling-Lawley 1 17.806408 391.741 1 22 1.6569e-15 ***
Roy 1 17.806408 391.741 1 22 1.6569e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
------------------------------------------
Term: factor2
Response transformation matrix:
factor21 factor22
[1,] 1 0
[2,] 0 1
[3,] -1 -1
[4,] 1 0
[5,] 0 1
[6,] -1 -1
Sum of squares and products for the hypothesis:
factor21 factor22
factor21 37201.09 21194.57
factor22 21194.57 12075.17
Multivariate Tests: factor2
Df test stat approx F num Df den Df Pr(>F)
Pillai 1 0.97709 447.7547 2 21 < 2.22e-16 ***
Wilks 1 0.02291 447.7547 2 21 < 2.22e-16 ***
Hotelling-Lawley 1 42.64331 447.7547 2 21 < 2.22e-16 ***
Roy 1 42.64331 447.7547 2 21 < 2.22e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
------------------------------------------
Term: factor1:factor2
Response transformation matrix:
factor11:factor21 factor11:factor22
[1,] 1 0
[2,] 0 1
[3,] -1 -1
[4,] -1 0
[5,] 0 -1
[6,] 1 1
Sum of squares and products for the hypothesis:
factor11:factor21 factor11:factor22
factor11:factor21 7.347826 9.608696
factor11:factor22 9.608696 12.565217
Multivariate Tests: factor1:factor2
Df test stat approx F num Df den Df Pr(>F)
Pillai 1 0.0155221 0.1655518 2 21 0.84852
Wilks 1 0.9844779 0.1655518 2 21 0.84852
Hotelling-Lawley 1 0.0157668 0.1655518 2 21 0.84852
Roy 1 0.0157668 0.1655518 2 21 0.84852
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
SS num Df Error SS den Df F Pr(>F)
(Intercept) 348911 1 136.48 22 56243.701 < 2.2e-16 ***
factor1 3185 1 178.88 22 391.741 1.657e-15 ***
factor2 9361 2 555.43 44 370.759 < 2.2e-16 ***
factor1:factor2 3 2 449.90 44 0.168 0.8459
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
factor2 0.89629 0.31676
factor1:factor2 0.91145 0.37775
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
factor2 0.90604 <2e-16 ***
factor1:factor2 0.91865 0.8284
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
HF eps Pr(>F[HF])
factor2 0.9827099 8.154032e-28
factor1:factor2 0.9983310 8.455884e-01
反復測定二元配置分散分析をR（carパッケージ）でやってみる（S. Kubota）
/assets/images/11207950/original/142f981a-6ea1-4ef8-93bb-320bf7fffe85?1668676934)
/assets/images/11207950/original/142f981a-6ea1-4ef8-93bb-320bf7fffe85?1668676934)
