★手計算とRで学ぶ統計学: 分散分析

まず，対応のない一元配置分散分析のサンプルサイズ設計を行いましょう。以下のコードをR上で実行すると

# (初めてパッケージを使用するときのみ行う) パッケージをインストールする 
install.packages("pwr") 
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む 
library(pwr) 
# 一元配置分散分析 (対応なし)
# サンプルサイズ設計 
pwr.anova.test(k = 3,
            n = NULL,
            f = 0.5,
            sig.level = 0.05,
            power = 0.8)

それぞれのフレーバーを最低でも14人に食べてもらう必要があるとわかりました。

Balanced one-way analysis of variance power calculation
               k = 3
               n = 13.89521
               f = 0.5
       sig.level = 0.05
           power = 0.8
 NOTE: n is number in each group 

なお，サンプルサイズ設計で算出された最低サンプルサイズを使用することが不適切だと判断した場合，計算結果よりもやや大きいサンプルサイズを使うこともあります (※1)。
今回はそれぞれのフレーバーを20人ずつに食べてもらう設定で分析を進めます。

(※1) あまりにサンプルサイズが大きい場合 (例: このデザインで1000人からデータを取り分析する) には第一種の過誤が発生しやすくなりますが，今回のように14人→20人程度なら第一種の過誤が発生する可能性はそこまで高くありません。

それでは，実際にデータをRに読み込み，データの確認を行ってみましょう。

まずは作業ディレクトリを設定します。次にデータを読み込んで変数に代入し，列名でデータを指定できるようにします。
分散分析の場合，カテゴリとして数字以外の文字を含む列 (今回だとflavor列) がデータ中に含まれることが多いです。
そのため，文字を含む列の内容をカテゴリとして認識させる必要があります。

# 作業ディレクトリを設定 
setwd("csvファイルが入っているフォルダ名") 
# 使用するデータの設定 
data_anova_between <- read.csv("anova_oneway_between.csv") 
# 列名でデータを指定できるようにする
attach(data_anova_between)
# flavor列をfactor型に変換 (コンピュータがflavor列の中身を文字列ではなくカテゴリとして認識できるようにする)
flavor <- as.factor(flavor)

データを読み込んだら，必ずデータの中身を確認しましょう。解析を進めた後にデータの中身が違うと，それまでの作業が徒労に終わります。

#データ内容の確認
head(data_anova_between)

データの頭から6行を確認すると，今回行う解析に用いるデータを正しく読み込めたことがわかります。

> head(data_anova_between)
   ID  flavor value
 1  1 caramel     4
 2  2 caramel     3
 3  3 caramel     5
 4  4 caramel     4
 5  5 caramel     4
 6  6 caramel     5

データを読み込んだあとは，箱ひげ図を表示してカテゴリごとの分布の形状を確認します。

# 箱ひげ図の表示 
boxplot(value ~ flavor, data = data_anova_between)

箱ひげ図を見る限りでは，いずれのデータにも極端な外れ値などがないことがわかります。

対応のない一元配置分散分析を実施するコードを以下に示します。
このコードをRで実行すると

# 対応のない1要因分散分析の実施 
# (初めてパッケージを使用するときのみ行う) パッケージをインストールする
install.packages("rstatix")
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む
library(rstatix)

# 分散分析の実施 (効果量は偏イータ2乗) 
anova_result <-anova_test(value ~ flavor,data = data_anova,
                           between = c(flavor), effect.size = "pes")
# 分散分析の結果の表示 
anova_result_between

結果がたちまちに得られます。

> anova_result_between 
ANOVA Table (type II tests) 

  Effect DFn DFd      F        p p<.05   pes 
1 flavor   2  57 33.528 2.37e-10     * 0.541

# DFnは要因の自由度 (独立変数の数 - 1)，DFdは残差の自由度 (サンプルサイズ - 要因の自由度 -1)
# pesはpartial eta squareの略

\( F (2, 57) = 33.53, p < .005\) という結果から，5%有意水準のもとで帰無仮説が棄却され，少なくとも1組のフレーバーの組み合わせで美味しさの評価に統計的に有意な差があることがわかりました。

rstatixパッケージを使用すると，分散分析の結果とともに効果量 (偏イータ二乗) も出力されます。
しかし，この効果量は偏イータ二乗なので，サンプルサイズ設計の際に使用した効果量 \(Cohen's/;f \)に変換する必要があります。

# 偏イータ2乗からf (サンプルサイズ設計で使用) に変換した値を変数anova_effsize_fに代入 
anova_effsize_between_f <- sqrt((anova_result_between$pes) / (1-anova_result_between$pes)) 
# Cohen's fをコンソール上に表示 
anova_effsize_between_f

上記のコードを実行すると，\(Cohen's/;f = 1.09\)という値が変数anova_effsize_between_fに格納されます。
サンプルサイズ設計の時には\(Cohen's/;f f = 0.5\) という中程度の効果が得られるサンプルサイズを求めたので，今回のケースの効果量が非常に大きかったことがわかります。

分散分析を行なって統計的に有意な差がある，つまりキャラメルとクッキークリーム，キャラメルとわさび，クッキークリームとわさびというフレーバーの組み合わせのうち少なくとも1つの美味しさの母平均に統計的に有意な差があることがわかりました。
ここからはどの組み合わせに差があるのかを，Bonferroni法による\(p\)値の補正をした多重比較法を用いて検討します。

以下のコードを実行すると

# パイプ演算子を使用するためのパッケージのインストール (※1)
#　(初めてパッケージを使用するときのみ行う) パッケージをインストールする
install.packages("tidyverse")
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む
library(tidyverse) 
# 多重比較 
pwc_between <- data_anova %>%   
pairwise_t_test(     
RT_ave ~ Figure, paired = FALSE, # 対応がないのでFALSE   
p.adjust.method = "bonferroni"   
)
# 多重比較の結果の表示
pwc_between

以下のような結果が得られます。p.adjは補正後の\(p\)値を示しています。

> pwc_between
# A tibble: 3 × 9 
  .y.   group1  group2    n1    n2             p p.signif         p.adj p.adj.signif 
* <chr> <chr>   <chr>  <int> <int>         <dbl> <chr>            <dbl> <chr>       
1 value caramel cookie    20    20 0.198         ns       0.594         ns           
2 value caramel wasabi    20    20 0.00000000026 ****     0.00000000078 ****         
3 value cookie  wasabi    20    20 0.0000000382  ****     0.000000115   ****      

この結果から，キャラメルとわさび，クッキークリームとわさびの美味しさの評定値の平均値に有意な差があることがわかりました。

(※1) パイプ演算子を使用するためにtidyverseというパッケージ群 (dplyrを内包する) を読み込んでいます。多重比較法自体は分散分析を行う際に使用したrstatixパッケージのみで実行可能です。

どのような仮説のもとで手元のデータに対して統計的仮説検定を行なったか，その結果はどうだったか，を論文やレポートに記述する際は以下のように記載します。

まず，対応のある一元配置分散分析のサンプルサイズ設計を行いましょう。
これまでの分析ではRのpwrパッケージを使用してサンプルサイズ設計を行いました。しかし，対応のある分散分析のサンプルサイズ設計をRで行うことは難しいです。そのため，G*Powerというフリーソフトウェアを使用してサンプルサイズ設計を行います。

G*Powerをダウンロードして起動すると，以下のような画面が表示されます。
G*Powerを用いて検定力分析を行うには，まずTest family (検定で使用する分布による検定手法の分類)，Statistical test (検定手法)，Type of power analysis (何が目的か) を設定します。

 G*Powerにおいて，検定族，検定手法，検定力分析の目的まで設定した状態の画面。

対応のある一元配置分散分析のサンプルサイズ設計では，図1のように
①Test family: "F tests" (F分布を利用する検定)
②: Statistical test: "ANOVA: Repeated measures, within factor" (反復測定，要因内)
③: Type of power analysis: "A priori: Compute required sample size - given \(\alpha \), power, and effect size" (有意水準 \(\alpha \)と検定力，検出したい効果の大きさを与えて必要なサンプルサイズを調べる)
と設定します。

次に，有意水準 \(\alpha \)と検定力，検出したい効果の大きさを設定します。

G*Powerにおいて，検定力分析のパラメータを設定した状態の画面。

今回の場合は，
①Effect size f (効果の大きさ): 0.25 (Cohenの基準だと中程度の効果量)
②\(\alpha \) err prob (第一種の過誤を犯す確率をどこまで許容するか): 0.05
③Power (\(\left(1 - \beta \right)\) err prob, 統計的に有意な差がある時にどのくらいの確率で正しく差を検出できればよいか): 0.8
④Number of groups (要因内の水準数): 3
⑤Number of measurements (反復回数): 10
⑥Corr among rep measures (反復測定したデータ間の相関係数): 0.5
⑦Nonsphericity correction \(\epsilon \) (水準間で差の分散が等しいかの仮定): 1
と設定しました。 

ここまで入力して"Calculate"をクリックすると，Output parametersのTotal sample sizeに先ほど設定した基準を満たすために必要なサンプルサイズをはじめとした計算結果が表示されます。

 要因内分散分析の結果。
太枠で囲んだ部分が必要サンプルサイズを示している。

今回は15名の参加者が必要であることがわかりました。

今回の分析では，平均値を使用します。しかし，Excelで一人ずつ反応時間の平均値を計算して…を繰り返すのはお勧めできません。参加者15人がそれぞれ3つの図形について試行を行っているので45回同じ作業を繰り返すのは非常に面倒ですし，確実にミスが発生します。
そこで，dplyrというパッケージを使用して，刺激ごとに各参加者の平均値を計算して，CSVファイルとして保存してみましょう (※1)。

まずはこれまでと同様に作業ディレクトリを設定してCSVファイルを読み込み，変数型の変換を行います。

# 作業ディレクトリの設定
setwd("csvファイルが入っているフォルダ名")
# 使用するデータの設定
data_anova_within_raw <- read.csv("anova_oneway_within_raw.csv")

# ID, 図形の形状 (Figure列) をfactor型に， 反応時間 (RT列) をnumeric型に変換 
data_anova_within_raw$RT <- as.numeric(data_anova_within_raw$RT) 
data_anova_within_raw$ID <-as.factor(data_anova_within_raw$ID) 
data_anova_within_raw$Figure<-as.factor(data_anova_within_raw$Figure)

# 列名でデータを指定できるようにする
attach(data_anova_within_raw)

次に，以下のコードによって平均値を計算します。
コード中に出現する%>% という記号はパイプ演算子といい，左辺にあるオブジェクトを右辺の関数の第一引数 (通常はデータフレーム) として渡す動作を行います。

# (初めてパッケージを使用するときのみ行う) パッケージをインストールする
install.packages("dplyr")
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む
library(dplyr)

# 参加者， 刺激ごとに反応時間の平均値を計算し，その結果を変数data_anovaに代入
data_anova_within_processed <- data_anova_within_raw %>% # group_byの第一引数として変数data_anova_rawを渡す
  group_by(ID, Figure) %>% # ID, Figure列でデータをグループ化し， summarise関数の引数に計算結果を渡す
   summarise(RT_ave = mean(RT, na.rm = TRUE)) # RT列の平均値を計算する

# データフレーム型に変換
data_anova_within_processed <- as.data.frame(data_anova_within_processed)

head() 関数で変数data_anova_within_processedの中身を確認すると，以下のように表示されます。RT_ave列に，図形ごとの平均反応時間が記録されています。

> head(data_anova_within_processed)
   ID    Figure RT_ave
 1  1 dodecagon  780.7
 2  1   octagon  730.8
 3  1      star  397.1
 4  2 dodecagon  742.8
 5  2   octagon  704.9
 6  2      star  522.4

最後に，data_anova_within_processedをcsvファイルとして保存します。これでdplyrパッケージを使わないバージョンの解析で使用するものと同じ形のが完成しました。

write.csv(data_anova_within_processed, 
"anova_oneway_within_processed.csv", 
row.names = FALSE # 列名をつけずに保存。ここをTRUEとすると， IDの左に1から始まる連番の列ができてしまう
)

(※1) もちろん，ここで処理したデータフレームをそのまま使って分散分析を行っても構いません。

それでは，実際にデータをRに読み込み，データの確認を行ってみましょう。

まずは作業ディレクトリを設定し，以下のコードでデータを読み込み，列名でデータを指定できるようにします。そして，カテゴリ変数 (ID, Figure) をfactor型に変換し，この2つの列のデータを文字の集まりではなく1つ1つの水準として扱えるようにします。

# 作業ディレクトリの設定
setwd("csvファイルが入っているフォルダ名")
# 使用するデータの設定
data_anova_within <- read.csv("anova_oneway_within.csv")
# 列名でデータを指定できるようにする
attach(data_anova_within)
# ID, 図形の形状 (Figure列) をfactor型に変換  
data_anova_within$ID <-as.factor(data_anova_within$ID) 
data_anova_within$Figure<-as.factor(data_anovawithin$Figure)

データを読み込んだら，必ずデータの中身を確認しましょう。解析を進めた後にデータの中身が違うと，それまでの作業が徒労に終わります。

#データ内容の確認
head(data_anova_within)

データの頭から6行を確認すると，今回行う解析に用いるデータを正しく読み込めたことがわかります。

> head(data_anova_within)
  ID    Figure RT_ave 
1  1 dodecagon  780.7 
2  1   octagon  730.8 
3  1      star  397.1 
4  2 dodecagon  742.8 
5  2   octagon  704.9 
6  2      star  522.4

データを読み込んだあとは，箱ひげ図を表示してカテゴリごとの分布の形状を確認します。

# 箱ひげ図の表示 
boxplot(RT_ave ~ Figure, data = data_anova_within)

箱ひげ図を見る限りでは，いずれのデータにも極端な外れ値 (明らかなミスによる値) などがないことがわかります。

対応のある一元配置分散分析を実施するコードを以下に示します。
このコードをRで実行すると

# 対応のある1要因分散分析の実施 
# (初めてパッケージを使用するときのみ行う) パッケージをインストールする
install.packages("rstatix")
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む 
library(rstatix) 

# 分散分析の実施 (効果量は偏イータ2乗) 
anova_result_within <- anova_test(data = data_anova_within, 
                          dv = RT_ave, 
                          wid = ID, 
                          within = Figure, 
                          effect.size = "pes")
# 分散分析の結果の表示 
anova_result_within

結果がたちまちに得られます。

> anova_result_within
ANOVA Table (type III tests) 
$ANOVA   
  Effect DFn DFd      F        p p<.05   pes 
1 Figure   2  28 94.158 3.71e-13     * 0.871 
# DFnは要因の自由度 (独立変数の数 - 1)，　DFdは残差の自由度 (サンプルサイズ - 要因の自由度 -1)
# pesはpartial eta squareの略

# 球面性の仮定 (差の分散が等しいかどうか) の検定。この場合は帰無仮説 (差の分散は等しくない) が棄却されている
$`Mauchly's Test for Sphericity` 
  Effect     W     p p<.05 
1 Figure 0.919 0.577       

# 球面性の補正をおこなった場合の値。今回の例では球面性を満たすためそこまで重要ではないかもしれない
$`Sphericity Corrections` 
  Effect   GGe     DF[GG]    p[GG] p[GG]<.05  HFe      DF[HF]    p[HF] p[HF]<.05 
1 Figure 0.925 1.85, 25.9 2.52e-12         * 1.06 2.12, 29.67 3.71e-13         *

\(F (2, 28) = 94.16, p < .001\) という結果から，5%有意水準のもとで帰無仮説が棄却され，図形の組み合わせのうち少なくとも1つに統計的に有意な差があることがわかりました。

rstatixパッケージを使用すると，分散分析の結果とともに効果量 (偏イータ二乗) も出力されます。
しかし，この効果量は偏イータ二乗なので，サンプルサイズ設計の際に使用した効果量 \(f\)に変換する必要があります。

# 偏イータ2乗をCohen's f (サンプルサイズ設計で使用) に変換した値を変数anova_effsize_within_fに代入 (※1)
anova_effsize_within_f <- sqrt((anova_result_within$ANOVA$pes) / (1-anova_result_within$ANOVA$pes)) 
# Cohen's fをコンソール上に表示 
anova_effsize_within_f

上記のコードを実行すると，\(Cohen's\;f \fallingdotseq 2.60\)という値が変数anova_effsize_fに格納されます。
サンプルサイズ設計の時には\(Cohen's\;ff = 0.5\) という中程度の効果が得られるサンプルサイズを求めたので，今回のケースの効果量が非常に大きかったことがわかります。

(※1) 対応のない一元分散分析の場合はanova_result_between$pes と入力することで効果量の値を取り出しました。
今回は結果を格納した変数anova_result_withinの中にに分散分析の結果，球面性の仮定についての検定，球面性の補正という3つのリストが格納されています。そして，それぞれのリスト中に計算結果がオブジェクトとして格納されているので，anova_result_within$pes という表記では偏イータ2乗を取り出すことができません。
そのため，anova_result_withinというオプジェクトに格納されたANOVAというリストのpesというオブジェクトを指定するために，anova_result_within$ANOVA$pesという表記で偏イータ2乗の値を取り出しています。

分散分析を行って統計的に有意な差がある，つまり星型と正八角形，星型と正十二角形，正八角形と正十二角形という刺激の形状の組み合わせのうち少なくとも1つの組み合わせで平均反応時間の母平均に統計的に有意な差があることがわかりました。
ここからはどの組み合わせに差があるのかを，Bonferroni法による\(p\)値の補正をした多重比較法を用いて検討します。

以下のコードを実行すると

# パイプ演算子を使用するためのパッケージのインストール (※1) 
#　(初めてパッケージを使用するときのみ行う) パッケージをインストールする 
install.packages("tidyverse") 
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む 
library(tidyverse) 

# 多重比較 
pwc_within <- data_anova_within %>%  
 pairwise_t_test( 
    RT_ave ~ Figure, paired = TRUE, # paired = TRUEとすることで対応がある分散分析ということを示す
    p.adjust.method = "bonferroni" 
  )

以下のような結果が得られます。p.adjは補正後の\(p\)値を示しています。

> pwc_within
# A tibble: 3 × 10
   .y.    group1    group2     n1    n2 statistic    df             p       p.adj p.adj.signif
 * <chr>  <chr>     <chr>   <int> <int>     <dbl> <dbl>         <dbl>       <dbl> <chr>       
1 RT_ave dodecagon octagon    15    15      1.12    14 0.284         0.852       ns           
2 RT_ave dodecagon star       15    15     12.5     14 0.00000000533 0.000000016 ****         
3 RT_ave octagon   star       15    15     10.1     14 0.0000000803  0.000000241 ****           

この結果から，星型と正八角形，星型と正十二角形の反応時間の平均値に有意な差があることがわかりました。

(※1) パイプ演算子を使用するためにtidyverseというパッケージを読み込んでいます。多重比較法自体は分散分析を行う際に使用したrstatixパッケージのみで実行可能です。

どのような仮説のもとで手元のデータに対して統計的仮説検定を行なったか，その結果はどうだったか，を論文やレポートに記述する際は以下のように記載します。