一元配置分散分析

データの読み込み

サンプルデータは、「2-6b　実験の計画とデータ（2-6b-a.xls）」[1]を使わせていただきます。このデータは4水準の温度における何かの成分量がまとめられています。

R Console

> x <- read.table(pipe("pbpaste"), header=TRUE) # クリップボードからデータを読み込み、xに代入。ちなみにこちらはMac用のコードでWindowsは、read.table("clipboard", header=TRUE)。 
> str(x) # xの構造を確認。
'data.frame':	5 obs. of  4 variables:
 $ c50: num  77.4 78.2 78.1 77.8 77.9
 $ c55: num  78.3 78.2 78.4 77.3 79.1
 $ c60: num  79.2 79.3 79.1 78.2 79.3
 $ c65: num  78.9 78.8 78.1 78.1 78.9
> y <- stack(x) # アンスタックからスタックに変換し、yに代入。 
> str(y) # yの構造を確認。
'data.frame':	20 obs. of  2 variables:
 $ values: num  77.4 78.2 78.1 77.8 77.9 78.3 78.2 78.4 77.3 79.1 ...
 $ ind   : Factor w/ 4 levels "c50","c55","c60",..: 1 1 1 1 1 2 2 2 2 2 ...
> colnames(y) <- c("value", "temperature") # 変数名がピンとこないので、成分量をvalue、温度をtemperetureに変更。 
> str(y) # yの変数名が変わっていることを確認。
'data.frame':	20 obs. of  2 variables:
 $ value      : num  77.4 78.2 78.1 77.8 77.9 78.3 78.2 78.4 77.3 79.1 ...
 $ temperature: Factor w/ 4 levels "c50","c55","c60",..: 1 1 1 1 1 2 2 2 2 2 ...
      16

[1] 松原望の総合案内サイト – 2-6b　実験の計画とデータ

基本統計量の計算

R Console

> by(y$value, y$temperature, summary)
y$temperature: c50
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  77.40   77.80   77.90   77.88   78.10   78.20 
------------------------------------------------------------------------------------ 
y$temperature: c55
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  77.30   78.20   78.30   78.26   78.40   79.10 
------------------------------------------------------------------------------------ 
y$temperature: c60
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  78.20   79.10   79.20   79.02   79.30   79.30 
------------------------------------------------------------------------------------ 
y$temperature: c65
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  78.10   78.10   78.80   78.56   78.90   78.90

等分散性の検定

一元配置分散分析はそれぞれの水準が等分散を仮定して計算します。等分散でない場合はノンパラメトリック検定としてKruskal-Wallis検定を行います。等分散性の検定にはF検定の一種であるLevene検定を使うことにします。こちらはcarパッケージを使うと便利です。Levene検定は二つ以上の母集団が等分散であることを帰無仮説としますので、p>0.05なら帰無仮説は棄却されず、母集団は等分散であると考えます。p<0.05なら等分散でないことになります。

R Console

> library(car) # パッケージcarを読み込み。
> leveneTest(y$value, y$temperature, center = mean) # Levene検定の実行。
Levene's Test for Homogeneity of Variance (center = mean)
      Df F value Pr(>F)
group  3  0.3783 0.7699

pが0.05より大きいので等分散であり、一元配置分散分析を行うことができます。

対応のない一元配置分散分析の実行

このデータが対応があるのかないのか定かではありませんので、まずは対応のない方法で計算します。

MEMO

一元配置分散分析には以下の方法があり、どれも同じ結果が得られる。オプション機能は違うのかもしれないが、詳しくは分からない。

oneway.test (変数 ~ 要因に)
aov (変数 ~ 要因)
anova ( lm(変数 ~ 要因))

R Console

> anova(lm(y$value ~ y$temperature)) # 一元配置分散分析を実行します。
Analysis of Variance Table

Response: y$value
              Df Sum Sq Mean Sq F value  Pr(>F)  
y$temperature  3  3.482 1.16067    5.13 0.01124 *
Residuals     16  3.620 0.22625                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

p<0.05で4条件の温度の間に差が認められました。

多重比較検定の実行

以下の３つがあります。ここではHolmの方法の結果を貼り付けておきます。

MEMO

TukeyHSD(aov(変数 ~ 要因)) # テューキーの方法

pairwise.t.test(変数, 要因, p.adjust.method=”holm”) # Holm法

pairwise.t.test(変数, 要因, p.adjust.method=”bonferroni”) # ボンフェローニの方法

R Console

> pairwise.t.test(y$value, y$temperature, p.adjust.method = "holm")

	Pairwise comparisons using t tests with pooled SD 

data:  y$value and y$temperature 

    c50    c55    c60   
c55 0.4493 -      -     
c60 0.0096 0.1123 -     
c65 0.1524 0.4493 0.4373

P value adjustment method: holm

50℃と60℃の間に差が認められました。

対応のある場合と対応のない場合の一元配置分散分析

上述したのは対応のない一元配置分散分析でした。

anova(lm(変数 ~ 要因))

これに対し、対応のある一元配置分散分析では以下のようになります。

anova(lm(変数 ~ 要因 + サンプル))

ここでいうサンプルとは下の表の1〜5です。例えば、1から5までビーカーに番号が付けられており、そのビーカーを温めたときのそれぞれの温度での成分量を調べるといった実験であれば対応のあることになります。一方で、適当に集めてきた物をそれぞれの温度まで温めて成分量を計測し、捨てるというような実験であれば対応のないことになります。つまり、対応のある実験では、成分量を説明する要因として温度とビーカーがあることになります。なので、「+サンプル」を加えるそうです[1]。

R Console

   c50  c55  c60  c65 # 要因（温度）
1 77.4 78.3 79.2 78.9 # サンプル1
2 78.2 78.2 79.3 78.8 # サンプル2
3 78.1 78.4 79.1 78.1 # サンプル3
4 77.8 77.3 78.2 78.1 # サンプル4
5 77.9 79.1 79.3 78.9 # サンプル5

実際に対応のある一元配置分散分析を実行してみます。

R Console

> sample <- factor(rep(c(1:5),4)) # 元のデータに1〜5までのサンプルが含まれていなかったので、加える。rep関数を使って1〜5までの数字を4回繰り返して書いて、それを順序なし要因と指定する。
> sample # データを確認。
 [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
Levels: 1 2 3 4 5
> anova(lm(y$value ~ y$temperature+sample)) # 要因にsampleを追加してanovaを実行。
Analysis of Variance Table

Response: y$value
              Df Sum Sq Mean Sq F value   Pr(>F)   
y$temperature  3  3.482 1.16067  8.8544 0.002279 **
sample         4  2.047 0.51175  3.9040 0.029571 * 
Residuals     12  1.573 0.13108                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

対応のない場合はp値が0.01124であったので、対応のある場合の方がp値が小さいことが分かります。

対応のある一元配置分散分析の多重比較検定では、paired=TRUEとオプションを指定して実行します。

R Console

> pairwise.t.test(y$value , y$temperature, p.adjust.method="holm", paired=TRUE ) # paired=TRUEと指定。

	Pairwise comparisons using paired t tests 

data:  y$value and y$temperature 

    c50   c55   c60  
c55 0.517 -     -    
c60 0.047 0.047 -    
c65 0.184 0.517 0.151

P value adjustment method: holm

対応のある場合では、55℃と60℃の間にも差が認められました。

[1] Introduction to Programming Language R – Rで統計学を学ぶ(6)

効果量の計算

一元配置分散分析の効果量η²は

η² ＝グループ間の平方和／平方和の合計

で計算することができます。分散分析の出力を利用する場合には以下のようにして計算します。

R Console

> result <- anova(lm(y$value~y$temperature)) #resultに結果を代入 
> result$"Sum Sq"[1]/(result$"Sum Sq"[1]+result$"Sum Sq"[2])
[1] 0.4902844

多重比較検定の効果量は組み合わせごとにt検定の効果量を計算します。

ANOVA君の利用

対応のない場合（被験者間）や対応のある場合（被験者内）が混在している二元配置分散分析などにも広く応用できる関数ANOVA君[1]はとても便利です。使い方は参考ページをご参照ください。参照ページからANOVA君のファイルをダウンロードしてrに読み込みます。

R Console

> source('~/Downloads/anovakun_482.txt', encoding = 'CP932') # 私はMac環境なので文字エンコードをCP932と指定。

ANOVA君はアンスタック形式のデータを用いて分析を行います。関数の書き方は参照ページを確認してください。

anovakun(データ, “要因計画の型”, 各要因の水準数，…)

と設定します。

ここでは対応のあるデータとして処理しようと思うので、参照ページで言うところの被験者内1要因のデータですので、要員計画の型を”sA”と記述し、各要因の水準は温度が4パターン（4水準）なので「4」と記述します。

R Console

> anovakun(x, "sA", 4, holm=TRUE, eta=TRUE) # 多重比較検定はデフォルトではShafferの方法ですが、ここでは上の結果と合わせるためholm=TRUEと指定。また、eta=TRUEでANOVAの効果量を算出するよう指定。

[ sA-Type Design ]

This output was generated by anovakun 4.8.2 under R version 3.4.3.
It was executed on Sun Feb  4 09:45:05 2018.

 
<< DESCRIPTIVE STATISTICS >>

---------------------------
  A   n     Mean    S.D. 
---------------------------
 a1   5  77.8800  0.3114 
 a2   5  78.2600  0.6427 
 a3   5  79.0200  0.4658 
 a4   5  78.5600  0.4219 
---------------------------


<< SPHERICITY INDICES >>

== Mendoza's Multisample Sphericity Test and Epsilons ==

-------------------------------------------------------------------------
 Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
-------------------------------------------------------------------------
      A  0.1356      2.7195   5 0.7575 ns  0.3333 0.6834 1.4104 0.6094 
-------------------------------------------------------------------------
                              LB = lower.bound, GG = Greenhouse-Geisser
                             HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller


<< ANOVA TABLE >>

---------------------------------------------------------
 Source     SS  df     MS  F-ratio  p-value      eta^2 
---------------------------------------------------------
      s 2.0470   4 0.5118                              
---------------------------------------------------------
      A 3.4820   3 1.1607   8.8544   0.0023 **  0.4903 
  s x A 1.5730  12 0.1311                              
---------------------------------------------------------
  Total 7.1020  19 0.3738                              
             +p < .10, *p < .05, **p < .01, ***p < .001


<< POST ANALYSES >>

< MULTIPLE COMPARISON for "A" >

== Holm's Sequentially Rejective Bonferroni Procedure ==
== The factor < A > is analysed as dependent means. == 
== Alpha level is 0.05. == 
 
---------------------------
  A   n     Mean    S.D. 
---------------------------
 a1   5  77.8800  0.3114 
 a2   5  78.2600  0.6427 
 a3   5  79.0200  0.4658 
 a4   5  78.5600  0.4219 
---------------------------

----------------------------------------------------------
  Pair     Diff  t-value  df       p   adj.p            
----------------------------------------------------------
 a2-a3  -0.7600   4.9472   4  0.0078  0.0467  a2 < a3 * 
 a1-a3  -1.1400   4.9241   4  0.0079  0.0467  a1 < a3 * 
 a3-a4   0.4600   3.0599   4  0.0377  0.1507  a3 = a4   
 a1-a4  -0.6800   2.5812   4  0.0612  0.1837  a1 = a4   
 a2-a4  -0.3000   1.3156   4  0.2587  0.5173  a2 = a4   
 a1-a2  -0.3800   1.2434   4  0.2816  0.5173  a1 = a2   
----------------------------------------------------------


output is over --------------------///

変数のラベルが、a1、a2…と変えられて表示されますが、結果自体は上の方法と同じです。

MEMO

ちなみに、Anova君でロング形式というまとめ方をされたデータを指定すれば、変数のラベルが結果に反映されるようです[3]。この方法は二元配置分散分析の方法をまとめるときに記載したいと思います。

ANOVA君は球面性の検定[2]の方法や効果量の算出方法を色々指定できるようです。

[1] 井関龍太のページ – ANOVA君

[2] 統計WEB – Mauchlyの球面性検定 | 統計用語集

[3] 井関龍太のページ – ANOVA君/より高度な入力方式

一元配置分散分析

データの読み込み

基本統計量の計算

等分散性の検定

対応のない一元配置分散分析の実行

多重比較検定の実行

対応のある場合と対応のない場合の一元配置分散分析

効果量の計算

ANOVA君の利用

1 COMMENT

コメントを残すコメントをキャンセル

データの読み込み

基本統計量の計算

等分散性の検定

対応のない一元配置分散分析の実行

多重比較検定の実行

対応のある場合と対応のない場合の一元配置分散分析

効果量の計算

ANOVA君の利用

1 COMMENT

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル