1
/
5

反復測定二元配置分散分析をR(carパッケージ)でやってみる(S. Kubota)

例えば,業務ソフトAとBを3カ月ずつ使用し,どちらが業務効率が良かったかを検証したいとする.

つまり,業務ソフト違いによる業務効率の影響と使い続けることでの時系列での業務効率の変化を検討する.このとき,業務ソフトの違い,時系列の変化をそれぞれ要因という.また,この検証を同一の社員で,AとBのソフトをそれぞれ3か月ずつ使い検証したこととする.

この場合は,同一の社員で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を組み合わせ,準備の2つ目が完了.

> 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
合同会社DILIGENCEでは一緒に働く仲間を募集しています
4 いいね!
4 いいね!

今週のランキング

窪田昂さんにいいねを伝えよう
窪田昂さんや会社があなたに興味を持つかも