第4章は統計を扱います。今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。シロート統計学はEZRを使った統計分析をわかりやすく解説されています。第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!今回は一元配置分散分析を紹介しますまず一元配置分散分析(ANOVA)についてはハルさんのサイトをご参照ください。1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。1〜4まででわからない部分があれば確認してみてください。第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。これを前提に次から進めていきます。次にRのコードを書くためのスクリプトファイルを作ります。3.データのダウンロード今回もハルさんのサイトのデータを使わせていただきます。デモデータ(ANOVA)この章ではRを使ってダウンロードしています。リンクのアドレスをコピーdestfileは保存場所と保存のファイル名を指定します。実際のコードは以下になります。"の中)とdestfileのdata/以降を変更するだけでOKです。destfile)dataフォルダにダウンロードできたことを確認します。4.データの読み込みデータを読み込みます。今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。Fileでデータを確認します。データが入っているセルを確認します。Datasetでデータを取り込みます。Sheet,Rangeを指定します。Name:ハルさんのサイトではgripでしたがgrip_anovaとします(大文字・小文字は別物とされます)Sheet:このExcelは1つしかデータがないのでDefaultのままでOKRange:先ほど確認したB2:C62Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。"B2:C62")View(grip_anova)データが正しく入っていることを確認します。これでデータの取り込みは完成です。5.正規分布の確認ハルさんのサイトに沿って正規分布の確認を行います。正規分布の確認は【4-1】Rでt検定を行う方法で紹介しています。以下のコードがイメージできない場合はまずこちらの記事をご参照ください。ヒストグラムで分布を確認します。1つずつグラフを作るのは手間がかかるので、まとめてグラフを作ります。ヒストグラムはgeom_histgram関数を使います。グラフが重なっていると見にくいのでfacet.gridでカテゴリーごとにグラフを分けます。facet.gridに関しては【3-4】Rのggplot2でヒストグラムを作るgeom_histogram関数で紹介しています。次にハルさんのサイトでは正規性の検定を行っています。正規性の検討はshapriro.test関数を3つのカテゴリーそれぞれに行う必要があります。今回はsplit関数とmap関数を使いまとめて行ってみます。下の図は【4-1】Rでt検定を行う方法で使った方法ですが今回のコードと見比べてみてください。コードがほとんど同じです。コードを使うと他でも役に立つことがわかります。map(~shapiro.test(.$grip))まとめて結果が出ました。すべて正規分布であることが棄却されませんでした。次に等分散性の確認を行います。3群以上で等分散性を検定するにはBartlett検定を使います。使えば$は必要ありません。どちらでも大丈夫です。grip_anova)8.分散分析(ANOVA)を行うそれでは本題の一元配置分散分析と多重比較を行います。実は一元配置分散分析と重回帰分析(t検定も)は同じ分析方法です(詳しく知りたい方は一般線形モデルで検索してみてください)。出てきた結果を分散分析の形式で表示するか重回帰分析の結果で表示するかの違いだけです。そのため一元配置分散分析を行うには分散分析のaov関数か重回帰分析のlm関数で行います。その後に行う多重比較ではlm関数が使えないものもあるので基本はaov関数で大丈夫です。どちらの関数もまずモデルを作りサマリーを表示するというのが一般的な流れです。全体的な流れは以下の通りです。sumamry関数だけ結果の表示が変わるので注意が必要です。grip_anova)summary.aov(model_grip_lm)結果はどちらもEZRと同じになりました。9.多重比較を行うハルさんのサイトでは多重比較でTukeyを紹介していますが他も紹介します。Tukeyだけ関数が違います。Tukeyの場合Tukeyの場合は先程のaov関数の結果を使います。lm関数の結果は使えませんので注意が必要です。TukeyHSD(model_grip_anova)Bonferroniの場合pairwise.t.test関数はモデルの結果を使いません。"hochberg")10.まとめ今回は一元配置分散分析を紹介しました。行うときはlongデータである必要があります。Excelでデータ収集するときに注意するか、wideデータの場合はlongデータに直してから分析を行ってください。longデータとwideデータについてはで【2-5】Rでデータを集計するのに便利なlongデータとgather関数で紹介しています。次回はKruskal-Wallis検定を行います。(追記)グラフの作成EZRでは検定と同時にグラフを作成しますが、Rでは直接作成します。【4-1】t検定でグラフを使ったときのコードがそのまま使えます。グラフを作るために必要な要約のコードもそのまま使います。c(0,50))コードの違いはわかったでしょうか?t検定で使ったデータであるgripを今回のデータであるgrip_anovaに変更しただけです。それ以外のコードの中身は全く一緒です!11.今回のコード今回使ったコードのまとめになります。第4章は統計を扱います。今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。シロート統計学はEZRを使った統計分析をわかりやすく解説されています。第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!今回はMcNemar検定を紹介しますまずMcNemar検定についてはハルさんのサイトをご参照ください。1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。1〜4まででわからない部分があれば確認してみてください。第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。これを前提に次から進めていきます。次にRのコードを書くためのスクリプトファイルを作ります。3.データのダウンロード今回もハルさんのサイトのデータを使わせていただきます。デモデータ(McNemar検定)この章ではRを使ってダウンロードしています。リンクのアドレスをコピーdestfileは保存場所と保存のファイル名を指定します。実際のコードは以下になります。"の中)とdestfileのdata/以降を変更するだけでOKです。destfile)dataフォルダにダウンロードできたことを確認します。4.データの読み込みデータを読み込みます。今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。Fileでデータを確認します。データが入っているセルを確認します。Datasetでデータを取り込みます。Sheet,Rangeを指定します。Name:ハルさんのサイトではBWでしたがbwとします(大文字・小文字は別物とされます)Sheet:このExcelは1つしかデータがないのでDefaultのままでOKRange:先ほど確認したB2:C22Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。"B2:C22")View(bw)データが正しく入っていることを確認します。これでデータの取り込みは完成です。5.データテーブルの作成Fisherの正確検定と同様にまずはデータテーブルを作成してみます。データテーブルを作成するのはtable関数でした。bw$After)これでデータテーブルができました。後で使うのでtable_bwという名前を付けておきます。bw$After)6.McNemar検定を行う次にMcNemar検定を行います。McNemar検定はmcnemar.test関数をを使います。ncnemar.test関数は2つの記述方法があります。1集計前の表(Excelの形式)を使う今回のExcelの形でmcnemar.test関数を使う事ができます。bw$After)2集計後の表(テーブル形式)を使う先程のテーブル形式(table_bw)をそのまま使うこともできます。他にも既にExcelで集計してしまった場合です。mcnemar.test(tableデータ)mcnemar.test(table_bw)どちらもEZRと同じ結果になります。7.まとめ今回はRでMcNemar検定を行う方法を紹介しました。Fisherの正確検定、カイ二乗検定、McNemar検定は似た使い方ですので続けて記事を見られた方は行いやすかったかもしれません。このページの最後に連続補正について紹介しています。次回は分散分析に移ります。bw$After)mcnemar.test(table_bw)9.(追記)連続修正による正規近似McNemar検定はカイ二乗検定を元に近似値を求めています。mcnemar.test関数はデフォルトで連続修正を用いた正規近似を用いています。FALSEを付ける必要があります。FALSE)サイトによって連続補正を付けるか付けないか意見はあるようです。なし(甘めな検定、第1種の過誤が減る)総じて言うと上記の傾向のようですが、気になる方は下記をご参照ください。こちらでは連続補正を施したほうが近似値が良くなるとありました。この近似はかなり正確なことと、p値がだいたい0.01くらいまでは連続修正を施した方が近似が良くなることがわかると思います。そのためウィルコクソンの1標本検定と違って、通常はzoの値によらず全て連続修正を施します。FALSEと言っています。correct=FALSE)を追記して下さい。-charこちらの記事では連続補正をすると第1種の過誤が生まれやすいので無条件に使うには連続補正を行わないほうが良いといった記載がありありました。こちらの記事ではnデータ数によって連続補正あり・なしを指定しています。パソコンの性能上厳密な検定ができずカイ二乗検定をつかって近似値を求めていたようですが、二項検定やパッケージを使ってp値を厳密に計算することもできるようです。それに関しては奥村先生の記事をご参照ください。第4章では統計の中でも検定を扱っています。ここまでいろいろな検定を行ってきましたが、グラフを作る時はどうするでしょうか?【4-6】Rのt検定の結果からp値や信頼区間の数値を取り出す方法では検定した結果を取り出す方法を紹介しました。今回は取り出したデータで95%信頼区間のグラフを作る方法を紹介します。データとt検定のコードは以下になります。前回の記事でこのコードの解説をしていますのでわからない場合は先に確認をお願いします。dat_height)1.データを取り出すまずはp値、95%信頼区間のデータを取り出してみます。これも前回の記事で紹介しています。res_height$conf.int[2]pconf.lowconf.highこの数字をグラフに乗せると小数点が多すぎるのでround関数で丸めます。今回はp値は小数点第3位、95%信頼区間は小数点第2位までにしてみます。2)pconf.lowconf.highこれでグラフに貼り付けやすくなりました。2.グラフを作成する今回はうまくいかないグラフも載せることでグラフを作る流れもお見せできればと思います。目標は以下のようなグラフを作ることです。第3章で使ってきたtidyverseパッケージのggplot関数を使います。ggplotの使い方がわからない場合はこちらをご参照ください。最初にtidyverseパッケージを呼び出します。library(tidyverse)まずグラフを作ってみます。今回使うのはgeom_errorbar関数を使います。geom_errorbar関数はaes関数の中に3つの要素が必要です。x:x軸ymin:エラーバーの最小値ymax:エラーバーの最大値x軸に値するものは今回ないので空欄にしておきます。yminとymaxはconf.lowとconf.highになります。conf.high))なんだかすごいグラフができました。グラフを横向きにするのはcoord_flip関数です。()には何も入れません。coord_flip()これで横向きになりましたが線が幅長いので短くします。今回は0.1にしました。widthはaes関数の外に配置します。coord_flip()これでエラーバーっぽくなりました。今度は数値を入れます。文字を打ち込むのでgeom_text関数を使います。geom_text関数のaes()では3つ指定します。coord_frip()を使っているのでx軸が縦、y軸が横になっていることに注意します。x:文字を置くx軸の位置→今回は空欄y:文字を置くy軸の位置→conf.lowとconf.highの2つlabel:実際の文字→conf.lowとconf.highの2つconf.lowとconf.highの2つのデータを入れるのでc関数でつなげます。coord_flip()数字が出てきましたがグラフとぶつかってしまいます。上に上げるときはマイナスの値を入れます。今回は試してみて-1.5にしました。coord_flip()これでグラフができましたが、y軸(グラフを横にしたから下がy軸になっている)の名前を変えてみます。labs関数を使って「95%信頼区間」と入れてみます。xも消します(空欄にする)。coord_flip()Windowの場合はこれで完成かもしれませんが、Macだと日本語が□□□と豆腐になってしまいます。今回はヒラギノ角ゴproW3を指定します。coord_flip()今回のt検定は2群の差が0と仮定していました。0のラインで赤線を足してみます。y軸に垂直な線を引くにはgeom_hline関数を使います。付けなければ黒になります。aesは付けなくて大丈夫です。coord_flip()3.まとめ今回は検定の結果から値を取り出してグラフを作成しました。グラフの要素を1つずつ追加しました。ブログに書かれているコードは長くて読みにくいかもしれませんが、「今度はこの要素を付け加えたい」という順番にコードを書けばその通りにグラフができるのもRの特徴の1つです。今回の順番は1例ですので順番を変えて1つずつ作ってみていただければggplotの理解も深まりやすいと思います。第4章では統計の中でも検定を扱っています。ここまでいろいろな検定を行ってきましたが、グラフを作る時はどうするでしょうか?これは【4-1】Rでt検定を行う方法で行ったt検定の結果です。p値があり、下には95%信頼区間、その下には各グループの平均が記載されています。-1.745123A:28.12B:35.07この数値をグラフなどにするには手で打ち直すかコピペをすることが多いと思います。今回はこの数値がどこに入っているのか?どうやって取り出せばいいのかという話です。今回は架空のデータを使いますがもし第4章を実践された方はそのデータでも大丈夫です。1.使うデータこんなデータを作りました。)を指定すると疑似データを作ってくれます。heightの列に男女のデータをc関数でくっつけ100人のデータにします。sexの列は(男性,男性,男性,,,,女性,女性,女性,,,,,)としたいのでrep関数を使いました。tibble関数でheightとsexをつなげました。2.t検定をするこれをt検定してみます。t検定は【4-1】Rでt検定を行う方法で紹介しています。dat_height)ここまでできました。このままだとデータは取り出せません。結果を変数に入れます。dat_height)3.結果からデータを取り出すすると右上のEnvironmentタブに結果が格納されます。リストというのはExcelのシートが10個あって、それぞれにデータが入っているというイメージです。左がExcelでいうシート名です。p値はp.valueに入っていますのでこうなります。res_height$p.valueres_height[["p.value"]]もし小数点以下を四捨五入したい場合はround関数を使います。3)95%信頼区間はconf.intに入っています。そのまま出すと2つのデータが出てきます(95%信頼区間の最小と最大)。]]ではないので注意してください。]は【1-12】Rで特定の条件にあう要素を抜き出す方法で紹介しましたが、条件式ではなくただの数値をいれると◯番目のデータを取り出してくれます。[1]だと1つ目のデータ(つまり最小)[2]だと2つ目のデータ(つまり最大)が取り出せます。res_height$conf.int[1]res_height[["conf.int"]][2]そして点推定である各平均はestimateに入っています。res_height$estimateres_height$estimate[1]res_height[["estimate"]][2]4.データを取り出すメリット結果をコピペしても問題はないのですが、こうすることでのメリットがあります。それはもし後でデータが変わったとしてもコード自体は何も変えなくていいことです。考えたくはありませんが、後でデータが増えた、実は個々の数値が間違ってた等あると結果も変わります。データをコードで取り出しておくと、データが変わっても自動的に数値が更新されるのでミスする可能性が減ります。元データを修正するだけでいいというのはExcelやEZRでは簡単にはできません。5.まとめ今回は検定の結果からデータを取り出す方法を紹介しました。検定の種類によっては取り出し方が違うこともありますが、まずは結果が格納されたデータをみるとヒントがあるかもしれません。次回は取り出した値を使って95%信頼区間のグラフを作る方法を紹介します。第4章は統計を扱います。今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。シロート統計学はEZRを使った統計分析をわかりやすく解説されています。第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか?といった視点で進めていきます。今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください!今回はFisherの正確検定を紹介しますまずFisherの正確検定についてはハルさんのサイトをご参照ください。1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。1〜4まででわからない部分があれば確認してみてください。第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。これを前提に次から進めていきます。次にRのコードを書くためのスクリプトファイルを作ります。3.データのダウンロード今回もハルさんのサイトのデータを使わせていただきます。デモデータ(Fisherの正確検定)この章ではRを使ってダウンロードしています。リンクのアドレスをコピーdestfileは保存場所と保存のファイル名を指定します。実際のコードは以下になります。"の中)とdestfileのdata/以降を変更するだけでOKです。destfile)dataフォルダにダウンロードできたことを確認します。4.データの読み込みデータを読み込みます。今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。Fileでデータを確認します。データが入っているセルを確認します。Datasetでデータを取り込みます。Sheet,Rangeを指定します。Name:ハルさんのサイトと同じsexとします(大文字・小文字は別物とされます)Sheet:このExcelは1つしかデータがないのでDefaultのままでOKRange:先ほど確認したB2:C42Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。"B2:C42")View(sex)データが正しく入っていることを確認します。これでデータの取り込みは完成です。5.データテーブルの作成データの要約は【2-6】Rでgroup_by関数とsummarize関数を使ってグラフ作成に必要な統計量(平均や標準偏差など)を求めるで紹介しました。group_by関数とsummarize関数を使って要約するならこうなります。ExcelのSexの列のSが大文字なので注意が必要です。n())これでもいいのですが、EZRのようにtable形式にしたい場合はtable関数を使います。sex$category)6.Fisherの正確検定を行うFishsrの正確検定を行うにはfisher.test関数を使います。そのままでわかりやすいです。列名を2つ指定するだけです。sex$category)EZRで0.205なので今回の結果を四捨五入すると同じ結果です。95%信頼区間も出ているのでグラフを作ってみました。1、つまりオッズ比が1であるかどうかで判断してるので1で線を引きます。95%信頼区間が1を挟んでいますのでpは0.05以上と判断できます。今回の信頼区間はかなり広いことも読み取れます。95%信頼区間はデータ数が増えると幅が狭くなります。(追記)7.χ二乗検定を行うにはχ(カイ)二乗検定を行うにはchisq.test関数を使います。fisherをchisqに変えるだけで中身は同じです。sex$category)8.まとめ今回はFisherの正確検定を紹介しました。【4-1】から進めている方は少しずつ慣れてきたでしょうか。このサイトはそのため第1章から順に読むと徐々に知識が追加され、途中で復習できるよう構成しています。わからない箇所が多ければサイトマップを見ていただければ別の発見があるかもしれません。次回は検定の結果から(p値や信頼区間)のデータを取り出す方法を紹介します。9.今回使ったコード今回使ったコードをまとめて置いておきます。95%信頼区間のコードも置いています。