第4章は統計を扱います。
今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。
シロート統計学はEZRを使った統計分析をわかりやすく解説されています。
第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか？といった視点で進めていきます。
今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください！
今回はFriedman（フリードマン）検定を紹介します
まずFriedman（フリードマン）検定についてはハルさんのサイトをご参照ください。
また1.準備〜4.データの読み込みまでは【4-1】Rでt検定を行う方法と全く同じ流れになります。
もし1〜4まででわからない部分があれば確認してみてください。
1.準備
第4章は毎回ExcelデータをダウンロードしてRを使うのでプロジェクトで管理して行うことを勧めています。
ここではR練習というプロジェクトを作り、Excelファイルを入れるためのdataフォルダを作っています。
これを前提に次から進めていきます。
2.スクリプトファイルの作成
次にRのコードを書くためのスクリプトファイルを作ります。
3.データのダウンロード
今回もハルさんのサイトのデータを使わせていただきます。
デモデータ（Friedman検定）
この章ではRを使ってダウンロードしています。
download.file(url = “ファイルのURL”,destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。
実際のコードは以下になります。
前回のコードのURL（" "の中）とdestfileのdata/以降を変更するだけでOKです。
#データのダウンロード url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-friedman-rank-sum-test.xlsx" destfile = "data/demo-friedman-rank-sum-test.xlsx"
download.file(url, destfile)
dataフォルダにダウンロードできたことを確認します。
4.データの読み込み
データを読み込みます。
今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。
View Fileでデータを確認します。
データが入っているセルを確認します。
B2からC22までデータが入っています（B2:D32と表記）
次にImport Datasetでデータを取り込みます。
Import画面ではName, Sheet，Rangeを指定します。
Name：ハルさんのサイトではgripでしたがgrip_friedmanとします（大文字・小文字は別物とされます）
Sheet：このExcelは1つしかデータがないのでDefaultのままでOK
Range：先ほど確認したB2:D32（実は今回のように表以外に余計なものがなければ空欄でもOK）
Importボタンを押す前に右にあるコードをコピーしスクリプトファイルに貼り付けることも忘れずに行います。
library(readxl)
grip_friedman <- read_excel("data/demo-friedman-rank-sum-test.xlsx")
View(grip_friedman)データが正しく入っていることを確認します。
これでデータの取り込みは完成です。
Friedman検定は繰り返しのデータ（被験者内要因）なのでデータを横につなげています。
5.Friedman検定を行う（wideデータの場合）
ハルさんのサイトでは早速Friedman検定を行っていますのでこちらも早速行います。
今回はwideデータなのでまずwideデータの場合を紹介します。
Friedman検定はfriedman.test関数を使います。
friedman.test(cbind(繰り返しデータのある列))
friedman.test(cbind(grip_friedman$W0,grip_friedman$W1,grip_friedman$W2))
EZRと同じ結果が出ました。
6.Friedman検定を行う（longデータの場合）
次にlongデータでFriedman検定を行う方法も紹介します。
friedman.test(データ ~ グループ | 個人を識別するid, data = ○○)
Friedman検定は対応のあるデータになります。
今回は30人がW0,W1,W2の3回データを取っているので3*30=90データあります。
そのため作るlongデータは上記のイメージになります。
今回はid, time, gripとしてますが自分がわかれば好きな名前を付けてOKです。
id：30人の内誰のデータかわかるためのid（氏名でも可）wideデータをlongデータに変えるにはgather関数、もしくはpivot_longer関数を使います。
time：W0, W1, W2の3つ
grip：今回測定した握力のデータ
gather関数に関しては【2-5】Rでデータを集計するのに便利なlongデータとgather関数で紹介しています。
データの扱いはtidyverseパッケージに含まれていますので。まずlibrary関数でtidyverseパッケージを読み込みます。今回はgrip_friedman_longという名前にします。
library(tidyverse)
#gather関数の場合 grip_friedman_long <- grip_friedman %>% rowid_to_column("id") %>% gather(key = time, value = grip, W0:W2, factor_key = TRUE)
#pivot_longer関数の場合 grip_friedman_long <- grip_friedman %>% rowid_to_column("id") %>% pivot_longer(cols = c(W0,W1,W2), names_to = "time", values_to = "grip", names_ptypes = list(time = factor()))
ポイントはもともとのExcelデータにidにあたる列がないことです。
そこでrowid_to_column関数を使います。""の中身が列名になります。
これで準備が整いました。
longデータでFriedan検定を行う場合は改めてこうなります。
friedman.test(データ ~ グループ | 個人を識別するid, data = ○○)
friedman.test(grip ~ time | id, data = grip_friedman_long)
もちろん結果は同じです。
7.多重比較を行う
Friedman検定で多重比較を行いたいのですがwideデータだとEZR専用のパッケージでしかできません。
そのためRで行うにはlongデータで行う必要があります。
longデータで多重比較を行うにはpairwise.wilcox.test関数を使います。簡単に言うと全ての組み合わせでMann-Whitney U検定を行いその結果を元にp値の調整を行います。
どの方法を使うかはp.adjust.method = で指定します。
またタイのエラーを回避するためにexact = FALSEを入れています。
data =
pairwise.wilcox.test(目的変数, factor(グループ), p.adjust.method = "◯◯", paired = TRUE, exact = FALSE)
#Bonferroniの方法
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "bonferroni", paired = TRUE, exact = FALSE)
#Holmの方法
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "holm", paired = TRUE, exact = FALSE)
EZRとほぼ同じ結果になりました。
8.どうしてもwideデータのままなら手計算が必要
wideデータのまま多重比較を行う方法がEZRのパッケージしかありません。
そのためどうしてもwideデータのまま行うにはwilcox.test関数、もしくはcoinパッケージのwilcox_test関数を使い、1つずつ組み合わせを行い、BonferroniやHolm、Hochbergそれぞれの方法で自分でp値を直接計算することになります。
この方法は統計ERさんの記事に詳しく紹介されています。
更にBonferroni、Holm、Hochbergって何をしてるの？という疑問も解決するかもしれません。
仕組みがわかればp値に2とか3とかかけるだけなので意外と単純です。
8.まとめ
今回はFriedman検定を紹介しました。
もしRで行う場合はwideデータで分析するかlongで分析するかをあらかじめイメージしていた方がその後の解析もスムーズかもしれません。
次は相関係数について紹介します！
9.今回使ったコード
#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-friedman-rank-sum-test.xlsx"
destfile = "data/demo-friedman-rank-sum-test.xlsx"
download.file(url, destfile)
library(readxl)
grip_friedman <- read_excel("data/demo-friedman-rank-sum-test.xlsx")
View(grip_friedman)
friedman.test(cbind(grip_friedman$W0,grip_friedman$W1,grip_friedman$W2))
#pivot_longer関数の場合
grip_friedman_long <-
grip_friedman %>%
rowid_to_column("id") %>%
pivot_longer(cols = c(W0,W1,W2), names_to = "time", values_to = "grip", names_ptypes = list(time = factor()))
#gather関数の場合
grip_friedman_long <-
grip_friedman %>%
rowid_to_column("id") %>%
gather(key = time, value = grip, W0:W2, factor_key = TRUE)
friedman.test(grip ~ time | id, data = grip_friedman_long)
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "bonferroni", paired = TRUE, exact = FALSE)
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "holm", paired = TRUE, exact = FALSE)
記事一覧はこちらになります。
