統計とR言語の勉強をしています。
「Rによるやさしい統計学」の写経で勉強。

#実行環境

Windows10 Pro 64bit
R：3.4.3
RStudio：1.1.383

RStudioのConsoleを使って実行しています

#一元配置分散分析(対応なし)を実行
3つ以上の平均値を比較するための統計的手法で分散分析があります。一元配置分散分析(対応なし)を例題から実施していきます。

項目	内容
前提	A, B, C, Dという4種類の指導法に分けた講座受講者に対して、同じテストを実施しました。テスト結果の得点に対して指導法の違いによる差異があるか
帰無仮説	4群の母平均は等しい(指導法の違いによる学習結果に差がない)
対立仮説	4群の母平均は等しくない(指導法の違いによる学習結果に差がある) ※4群のうち1つでも母平均が異なる場合も対立仮説が支持
検定統計量	$\frac{\frac{群間平方和}{郡間自由度}}{\frac{群内平方和}{郡内自由度}}$
有意水準	5%

Rで書いていくと、検定の手前まではこんな感じ。

> A <- c(15,9,18,14,18)
> B <- c(13,8,8,12,7)
> C <- c(10,6,11,7,12)
> D <- c(10,7,3,5,7)
> statistics_test2 <- c(A,B,C,D)
> statistics_test2
 [1] 15  9 18 14 18 13  8  8 12  7 10  6 11  7 12 10  7  3  5  7
> shidouhou <- c(rep("A",5),rep("B",5),rep("C",5),rep("D",5))
> shidouhou2 <- factor(shidouhou)
> shidouhou2
 [1] A A A A A B B B B B C C C C C D D D D D
Levels: A B C D

検定は3種類の関数があります。
まずは、oneway.test関数です。一元配置分散分析のみ実行可能な関数です。

> oneway.test(statistics_test2~shidouhou2, var.equal = TRUE)

	One-way analysis of means

data:  statistics_test2 and shidouhou2
F = 7.1111, num df = 3, denom df = 16, p-value = 0.002988

次にaov関数です。最も一般的らしいです。summary関数も付加しないと見にくいです。

> aov(statistics_test2~shidouhou2)
Call:
   aov(formula = statistics_test2 ~ shidouhou2)

Terms:
                shidouhou2 Residuals
Sum of Squares         184       138
Deg. of Freedom          3        16

Residual standard error: 2.936835
Estimated effects may be unbalanced
> summary(aov(statistics_test2~shidouhou2))
            Df Sum Sq Mean Sq F value  Pr(>F)   
shidouhou2   3    184   61.33   7.111 0.00299 **
Residuals   16    138    8.63                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

最後にanova関数です。複数のモデル比較など高度な分析にも対応できるようです。

> anova(lm(statistics_test2~shidouhou2))
Analysis of Variance Table

Response: statistics_test2
           Df Sum Sq Mean Sq F value   Pr(>F)   
shidouhou2  3    184  61.333  7.1111 0.002988 **
Residuals  16    138   8.625                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

どれもp値が5%以下なので、帰無仮説が棄却され、対立仮説が支持されます。

では、検定の中身についての解説です。1コマンドで検定しましたが、中では以下のような計算をしています。
まずは、全平均行列と群平均行列(列の平均行列)を作ります。

元データと全平均行列の偏差自乗を作り、合計を全体自由度で割ります。

群間行列(群平均と全平均行列の偏差自乗)を作り、合計を群間自由度で割ります。

郡内行列(元データと群平均の偏差自乗)を作り、合計を郡内自由度で割ります。

最後に群間平均平方を郡内平均平方で割ることで検定統計量を算出します。
今までの流れをRで実行すると下記のとおりです。

> all_data <- cbind(A,B,C,D)
> all_data
      A  B  C  D
[1,] 15 13 10 10
[2,]  9  8  6  7
[3,] 18  8 11  3
[4,] 14 12  7  5
[5,] 18  7 12  7

#全平均行列作成
> all_average <- mean(all_data)
> all_average
[1] 10
> all_average_matrix <- matrix(rep(all_average,20),nrow=5,ncol=4)
> all_average_matrix
     [,1] [,2] [,3] [,4]
[1,]   10   10   10   10
[2,]   10   10   10   10
[3,]   10   10   10   10
[4,]   10   10   10   10
[5,]   10   10   10   10

#群平均行列作成
> gun_average <- colMeans(all_data)
> gun_average
   A    B    C    D 
14.8  9.6  9.2  6.4 
> gun_average_matrix <- matrix(rep(gun_average,5),nrow=5,ncol=4,byrow=TRUE)
> gun_average_matrix
     [,1] [,2] [,3] [,4]
[1,] 14.8  9.6  9.2  6.4
[2,] 14.8  9.6  9.2  6.4
[3,] 14.8  9.6  9.2  6.4
[4,] 14.8  9.6  9.2  6.4
[5,] 14.8  9.6  9.2  6.4

#全体データ作成
> zentai <- all_data - all_average_matrix
> zentai
      A  B  C  D
[1,]  5  3  0  0
[2,] -1 -2 -4 -3
[3,]  8 -2  1 -7
[4,]  4  2 -3 -5
[5,]  8 -3  2 -3

#群間データ作成
> gunkan <- gun_average_matrix - all_average_matrix
> gunkan
     [,1] [,2] [,3] [,4]
[1,]  4.8 -0.4 -0.8 -3.6
[2,]  4.8 -0.4 -0.8 -3.6
[3,]  4.8 -0.4 -0.8 -3.6
[4,]  4.8 -0.4 -0.8 -3.6
[5,]  4.8 -0.4 -0.8 -3.6

#郡内データ作成
> gunnai <- all_data-gun_average_matrix
> gunnai
        A    B    C    D
[1,]  0.2  3.4  0.8  3.6
[2,] -5.8 -1.6 -3.2  0.6
[3,]  3.2 -1.6  1.8 -3.4
[4,] -0.8  2.4 -2.2 -1.4
[5,]  3.2 -2.6  2.8  0.6
> zentai^2
      A B  C  D
[1,] 25 9  0  0
[2,]  1 4 16  9
[3,] 64 4  1 49
[4,] 16 4  9 25
[5,] 64 9  4  9
> gunkan^2
      [,1] [,2] [,3]  [,4]
[1,] 23.04 0.16 0.64 12.96
[2,] 23.04 0.16 0.64 12.96
[3,] 23.04 0.16 0.64 12.96
[4,] 23.04 0.16 0.64 12.96
[5,] 23.04 0.16 0.64 12.96
> gunnai^2
         A     B     C     D
[1,]  0.04 11.56  0.64 12.96
[2,] 33.64  2.56 10.24  0.36
[3,] 10.24  2.56  3.24 11.56
[4,]  0.64  5.76  4.84  1.96
[5,] 10.24  6.76  7.84  0.36

# 全体平方和(群間平方和 + 郡内平方和と同じ)
> all_square_sum <- sum(zentai^2)
> all_square_sum
[1] 322

# 群間平方和(aov関数結果のInputに対するSum of Squares)
> gunkan_square_sum <- sum(gunkan^2)
> gunkan_square_sum
[1] 184

# 郡内平方和(aov関数結果のResiduals(残差)のSum of Squares)
> gunnai_square_sum <- sum(gunnai^2)
> gunnai_square_sum
[1] 138

#自由度
> all_jiyudo <- length(all_data) - 1
> gunkan_jiyudo <- ncol(all_data) -1
> gunnai_jiyudo <- (nrow(all_data)-1)*ncol(all_data)

#群間平均平方(aov関数結果のInputに対するMean Sq)
> gunkan_average_square <- gunkan_square_sum / gunkan_jiyudo
> gunkan_average_square
[1] 61.33333

#郡内平均平方(aov関数結果のResiduals(残差)のMean Sq)
> gunnai_average_square <- gunnai_square_sum / gunnai_jiyudo
> gunnai_average_square
[1] 8.625

#全体平均平方
> all_square_sum/all_jiyudo
[1] 16.94737

# 全体平均平方(別解)
> var(statistics_test2)
[1] 16.94737

# 検定統計量(F分布)
> fstat <- gunkan_average_square/gunnai_average_square
> fstat
[1] 7.111111

##Tukey
ここまでの比較だと4つの指導法の母平均が等しくないことまでは検定でわかりました。さらに個別の指導法間の差を調べるために、Tukey(テューキー)の方法と呼ばれる多重比較をします。
検定統計量は以下の式で求めます。

{\frac{|比較する群の平均値差|}{\sqrt{\frac{郡内の平均平方}{各群のデータ数}}}}

検定統計量6.395651が棄却域4.046093以上なので対立仮説「母平均が等しくない」は支持されます。

> q <- abs(mean(A)-mean(D))/sqrt(gunnai_average_square/nrow(all_data))
> q
[1] 6.395651
> qtukey(0.95,4,16)
[1] 4.046093

関数TukeyHSDを使うことで一度に多重比較が出来ます。diffは平均値差、lwrとuprが95%信頼区間の下限・上限、p adjがp値です。

> TukeyHSD(aov(statistics_test2~shidouhou2))
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = statistics_test2 ~ shidouhou2)

$shidouhou2
    diff        lwr        upr     p adj
B-A -5.2 -10.514108  0.1141085 0.0562227
C-A -5.6 -10.914108 -0.2858915 0.0371222
D-A -8.4 -13.714108 -3.0858915 0.0017736
C-B -0.4  -5.714108  4.9141085 0.9963241
D-B -3.2  -8.514108  2.1141085 0.3446966
D-C -2.8  -8.114108  2.5141085 0.4561325