★手計算とRで学ぶ統計学: ノンパラメトリック検定

まず，適合度の検定のサンプルサイズ設計を行いましょう。以下のコードをR上で実行すると

# (初めてパッケージを使用するときのみ行う) パッケージをインストールする  
install.packages("pwr")  
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む  
library(pwr) 
# 適合度の検定
# サンプルサイズ設計
pwr.chisq.test(N = NULL, # サンプルサイズを調べたいので設定しない
           w = 0.2, # 小さな効果量 (Cohen's w) を検出したい
           sig.level = 0.05, # 有意水準
           power = 0.8, # 検定力
           df = (3-1) # 自由度。カテゴリ数 (今回はSSR, SR, Rの3つ) - 1
　　　　　　)

最低でも241個のデータが必要であることがわかりました。

      Chi squared power calculation
               w = 0.2
               N = 240.8672
              df = 2
       sig.level = 0.05
           power = 0.8
 NOTE: N is the number of observations # Nは観測数 (データ全体の数) であることに留意

しかし，今回は10連ガチャを回すという試行の最大数が20回までという制約があります。つまり得られるデータ数は最大で200個しかないのです。
このような場合は，論文やレポートのMethod欄に「サンプルサイズ設計の結果，必要データ数は○○個であると分かった。しかし，システムの制約上今回は200個のデータを使用して分析を行った」と，サンプルサイズ設計の結果と実際のデータ数が異なる理由を記載しましょう。

適合度の検定を実施するコードを以下に示します。
これまでの検定では，csvファイルを読み込ませ，分布を確認するところから始めていました。
しかし，適合度の検定ではcsvファイルを作るよりも実測値や出現比率を数列として変数に代入したほうがより簡単であるため，今回はR上で直接データを打ち込むことにします。

# 適合度の検定の実施 
# 実際のデータを変数に格納 (gofはgoodness of fitの略称) 
data_gof <- c(5, 51, 144) 
# 公表されている排出率を変数に格納 
prob <- c(0.03, 0.18, 0.79) 
# 適合度の検定 
chisq.test(x = data_gof, p = prob)

適合度の検定の結果がたちまちに得られます。

    Chi-squared test for given probabilities 
data:  data_gof # 使用したデータ
X-squared = 7.6572, # カイ二乗値
df = 2, # 自由度
p-value = 0.02174 # p値

\(\chi^2(2) = 7.65, p = .022 \)  という結果から，5%有意水準のもとで帰無仮説が棄却され，公表された排出率と手元のデータの比率に統計的に有意な差があることがわかりました。

適合度の検定では，Cohen's wが効果の大きさの指標として用いられます。以下のコードでCohen's wを計算してみましょう。 

# 効果量 (Cohen's w) の計算 
# カイ二乗値をデータ数で割ったものの平方根を取る
effsize_gof <- sqrt(result_gof$statistic / sum(data_gof))
# 効果量の表示
effsize_gof

上記のコードを実行すると，効果量が\(Cohen's\,w = 0.1956677 ...\) であったことがわかります。
5%有意水準のもとでは公表されている排出率と統計的な有意差はあります。しかし，効果としてはそこまで大きくない，というイメージを持っておくとよいでしょう。

> effsize_gof
X-squared 
0.1956677 

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

まず，独立性の検定のサンプルサイズ設計を行いましょう。以下のコードをR上で実行すると

# (初めてパッケージを使用するときのみ行う) パッケージをインストールする  
install.packages("pwr")  
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む  
library(pwr) 
# 独立性の検定
# サンプルサイズ設計 
pwr.chisq.test(N = NULL, # サンプルサイズを調べたいので設定しない
                w = 0.5, # 比較的大きな効果量 (Cohen's w) を検出したい
                sig.level = 0.05, # 有意水準
                power = 0.8, #　検出力
                df = (2-1)*(2-1) # (行数 (好き, 嫌い) -1) * (列数 (普通ピーマン，品種改良ピーマン) -1)
                )

最低でも32個のデータが必要であることがわかりました。

     Chi squared power calculation
               w = 0.5
               N = 31.39544
              df = 1
       sig.level = 0.05
           power = 0.8
 NOTE: N is the number of observations　# Nは観測数 (データ全体の数) であることに留意

今回は食育イベントの一環として行うことを想定して，参加者全員 (50名) からデータを集めることにします。

独立性の検定を実施するコードを以下に示します。
適合度の検定と同様，R上で直接データを打ち込んでから検定を行います。

# 独立性の検定の実施 
# データを変数に格納 
data_ind <- matrix(c(32, 23, 18, 27), nrow = 2, byrow = T) 
# 独立性の検定 
result_ind <- chisq.test(x = data_ind)

上記のコードを実行すると，独立性の検定の結果がたちまちに得られます。
Rのchisq.test関数では，2\(\times \)2のクロス集計表の検定を行う際はデフォルトでイェーツの補正がなされます。

    Pearson's Chi-squared test with Yates' continuity correction 
data:  data_ind # 使用したデータ
X-squared = 2.5859, # カイ2乗値
df = 1, # 自由度
p-value = 0.1078 # p値

\(\chi^2(1) = 2.59, p = .108 \)  という結果から，5%有意水準のもとで帰無仮説に関する判断が保留され，ピーマンの品種と好き嫌いの区分に関係があるとはいえないことがわかりました。

適合度の検定と同様，独立性の検定でもCohen's wが効果の大きさの指標として用いられます。以下のコードでCohen's wを計算してみましょう。 

# 効果量 (Cohen's w) の計算 
# カイ二乗値をデータ数で割ったものの平方根を取る
effsize_ind <- sqrt(result_gof$statistic / sum(data_gof))
# 効果量の表示
effsize_ind

上記のコードを実行すると，効果量が(Cohen's\,w = 0.1608061 ...\) であったことがわかります。
5%有意水準のもとでも帰無仮説についての判断が保留されているので，効果の大きさについての議論も一旦保留しておきましょう。

> effsize_ind
X-squared 
0.1608061 

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

まず，作業ディレクトリを設定し，データを読み込みます。
今回のデータには文字を含む列 (major列) があります。そのため，major列をカテゴリとして扱うようRに認識させる必要があります。

# 作業ディレクトリの設定
setwd("csvファイルを入れたフォルダの名前")
# データの読み込み
data_Utest <- read.csv("Utest.csv")
# major列をfactor型に変更
data_Utest$major <- as.factor(data_Utest$major)

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

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

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

> head(data_Utest)
   ID major value
 1  1   sci     4
 2  2   sci     5
 3  3   sci     1
 4  4   lit     2
 5  5   sci     5
 6  6   lit     2

データを読み込んだあとは，棒グラフによって文系，理系ごとの分布の形状を確認します。

# データの図示
# 文系のデータ (major列がlit) を変数Utest_data_litに代入
Utest_data_lit <- subset(data_Utest,major=="lit") 
# 文系のデータのみを棒グラフとして図示
barplot(table(Utest_data_lit$value)) 

# 理系のデータ (major列がsci) を変数Utest_data_sciに代入
Utest_data_sci <- subset(data_Utest,major=="sci") 
# 理系のデータのみを棒グラフとして図示
barplot(table(Utest_data_sci$value))

文系と理系それぞれのデータを確認すると，分布が大きく異なることがわかります。

文系の授業満足度の分布

理系の授業満足度の分布

それでは，マン・ホイットニーのU検定を実施するコードを以下に示します。

# マン・ホイットニーのU検定 
wilcox.test(Utest_data_lit$value, Utest_data_sci$value)

上記のコードを実行すると，マン・ホイットニーのU検定の結果がたちまちに得られます。
実行結果に"Warning message"とありますが，データ中に同じ値が存在すると必ずこの警告メッセージが現れます。そのため，この警告についてはそこまで気にする必要はありません。

> wilcox.test(Utest_data_lit$value, Utest_data_sci$value)
     Wilcoxon rank sum test with continuity correction 

data:  Utest_data_lit$value and Utest_data_sci$value 
W = 147, p-value = 0.00178 
alternative hypothesis: true location shift is not equal to 0 

Warning message: In wilcox.test.default(Utest_data_lit$value, Utest_data_sci$value) : 
  cannot compute exact p-value with ties

\(W = 147, p = .002 \)  という結果から，5%有意水準のもとで帰無仮説が棄却され，学生の文系 / 理系によって授業の満足度に統計的に有意な差があることがわかりました。

マン・ホイットニーのU検定ではCliff's deltaが効果の大きさの指標として用いられます。以下のコードでCliff's deltaを計算してみましょう。 

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

# 効果量の計算
cliff.delta(Utest_data_lit$value, Utest_data_sci$value)

上記のコードを実行すると，\(Cliff's\,\Delta = -0.5288462 ...\) と大きな効果量であったことがわかりました。

Cliff's Delta 

delta estimate: -0.5288462 (large) 
95 percent confidence interval:
      lower      upper
 -0.7610288 -0.1765589 

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