第4章は統計を扱います。今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。シロート統計学はEZRを使った統計分析をわかりやすく解説されています。第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!今回はKruskal-Wallis(クラスカル・ウォリス)検定を紹介しますまずKruskal-Wallis(クラスカル・ウォリス)検定についてはハルさんのサイトをご参照ください。1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。1〜4まででわからない部分があれば確認してみてください。第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。これを前提に次から進めていきます。次にRのコードを書くためのスクリプトファイルを作ります。3.データのダウンロード今回もハルさんのサイトのデータを使わせていただきます。デモデータ(Kraskal-Wallis検定)この章ではRを使ってダウンロードしています。リンクのアドレスをコピーdestfileは保存場所と保存のファイル名を指定します。実際のコードは以下になります。"の中)とdestfileのdata/以降を変更するだけでOKです。destfile)dataフォルダにダウンロードできたことを確認します。4.データの読み込みデータを読み込みます。今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。Fileでデータを確認します。データが入っているセルを確認します。Datasetでデータを取り込みます。Sheet,Rangeを指定します。Name:ハルさんのサイトではBWでしたがbw_kwとします(大文字・小文字は別物とされます)Sheet:このExcelは1つしかデータがないのでDefaultのままでOKRange:先ほど確認したB2:C92Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。"B2:C92")View(bw_kw)データが正しく入っていることを確認します。これでデータの取り込みは完成です。5.正規分布の確認ハルさんのサイトに沿って正規分布の確認を行います。正規分布の確認は【4-1】Rでt検定を行う方法で紹介しています。以下のコードがイメージできない場合はまずこちらの記事をご参照ください。ヒストグラムで分布を確認します。1つずつグラフを作るのは手間がかかるので、まとめてグラフを作ります。ヒストグラムはgeom_histgram関数を使います。グラフが重なっていると見にくいのでfacet.gridでカテゴリーごとにグラフを分けます。ちなみに【4-10】Rで分散分析(一元配置分散分析)を行う方法で作ったヒストグラムのコピペでできてしまいます。変数名を変えるだけです。次にハルさんのサイトでは正規性の検定を行っています。正規性の検討はshapriro.test関数を3つのカテゴリーそれぞれに行う必要があります。今回はsplit関数とmap関数を使いまとめて行ってみます。これも【4-10】Rで分散分析(一元配置分散分析)を行う方法で作ったコードと全く同じです。変数名を変えるだけです。コードを使うと他でも役に立つことがわかります。map(~shapiro.test(.$BW))まとめて結果が出ました。0.05ですので正規分布であるという仮説が棄却されました。7.グラフを作る今回はノンパラメトリックの検定を行うので箱ひげ図を作成します。箱ひげ図の作り方は【3-6】Rのggplot2で箱ひげ図を作るgeom_boxplot関数で紹介しています。BW))もしEZRと同じようなグラフであればboxplot関数も使えるので併せて消化しします。bw_kw)8.Kruskal-Wallis(クラスカル・ウォリス)検定を行う正規性が棄却されましたので今回はノンパラメトリックのKruskal-Wallis(クラスカル・ウォリス)検定を行います。Kruskal-Wallis検定はEZRではkruskal.test関数を使います。data=bw_kw)EZRと同じ関数を使っているので同じ結果になりました。検定のときに紹介したcoinパッケージでも行うことができます。coinパッケージを使ったことがなければインストールします。#coinパッケージのインストール#1度でもしていればしなくてOKinstall.packages("coin")coinパッケージのkruskal_test関数はkruskal.test関数と同じ使い方です。使う前にlibrary関数で呼び出します。グループにはfactor関数もしくはas.factor関数でfactor型に変更します。data=bw_kw)Asymptotic(漸近的:正規近似)とあるので正確検定ではないようです。そのためEZRと同じ結果になりました。実は3以上のグループではなく2変数だと正確検定を行うことができます。9.多重比較(Steel-Dwass法)を行うEZRで使っていたSteel-Dwass法はEZR専用のパッケージなので、Rで行うにはNSM3パッケージのpSDCFlig関数を使う必要があります。ちなみにEZRのSteel-Dwass法は正確検定ではなく正規近似を行います。NSM3パッケージではデータ数が少なければ正確近似ができますが、データ数が大きくなると計算量が膨大になりすぎて実施できません。モンテカルロ法もしくは正規近似を行うことになります。今回の90例では正確近似はできませんでした。検定が基準に詳しく紹介されています。factor(グループ))目的変数:bw_kw$BWグループ:factor(bw_kw$category)*グループはfactor型にする必要があるためfactor関数かas.factor関数を使います。factor(bw_kw$category))それぞれのp値は以下になりました。大きな差はありません。A-B:正規近似0.014、モンテカルロ0.0142A-C:正規近似0.0255、モンテカルロ0.0249B-C:正規近似0.9006、モンテカルロ0.90210.ひとまずまとめ今回はRでKruskal-Wallis(クラスカル・ウォリス)検定と多重比較であるSteel-Dwass法の紹介をしました。U検定でさんざん問題になったタイと正確検定に関してはどうなったんだ?という方は続きがありますので興味がある方は続きをご覧ください。次回は反復測定分散分析を紹介します。多重比較ではpairwise.t.test関数を使いましたが、Kruskal-Wallisではpairwise.wilcox.test関数を使います。U検定を行いその結果を元にp値の調整を行います。"hochberg")ただ、この方法だと警告が出ます。タイ(同順位)があるためです。U検定で紹介しました。FASEを追記します。すればエラーは出ません。U検定を行いその結果を元にp値の調整を行います。pairwise_wilcox_test関数が存在せず自分で手計算するしかありません。そのためどうしても正確検定を行いたい場合はcoinパッケージのwilcox_test関数を使い、1つずつ組み合わせを行い、BonferroniやHolm、Hochbergそれぞれの方法で自分でp値を直接計算することになります。この方法は統計ERさんの記事に詳しく紹介されています。Bonferroni、Holm、Hochbergって何をしてるの?という疑問も解決するかもしれません。仕組みがわかればp値に2とか3とかかけるだけなので意外と単純です。12.本当にまとめ今回はRでKruskal-Wallis(クラスカル・ウォリス)検定と多重比較であるSteel-Dwass法の紹介をしました。U検定のちゅんたろさんのツイートがここでも参考になると思います。10:40:35そしてそもそもp値がより小さければより差が大きいわけではありませんので、正確検定か正規近似か?よりも研究デザインが適切かどうかの方が大切になるかもと言うのが個人的なイメージです。次回は反復測定分散分析を紹介します!13.今回使ったコード今回使ったコードを紹介します。"hochberg")#正規近似を使えば警告は出ない。