Rを使って解析方法をお伝えする。今回はCochranQ検定と、あと解析でMcNemer検定の多重比較になります。これをlibrary：”rstatix”を使って解析を行います。
解析について知りたい場合には下記のサイトを確認してください。
対象となるデータは？
1.2項分布 ： yes/noなどの
2.関連のあるデータ
3.3項目以上になります。
yes/no、あり/なしなど、関連のあるデータの検定としてCockranQ検定・McNemer検定があります。今回はその検定をRで行う方法を解説したいと思います。
手順は？
1.rstatixを使えるようにする
2.データを縦型のデータに変える
3.cochranQ検定を行う
4.cochranQ検定で有意差があれば、McNemer検定（holm）を行う
の4つの工程を行っていきます
まずはrstatix・を使えるようにする
今回の解析では解析では”rstatix”を使っていきますが、それには縦型のデータにする必要があります。そこで”tidyverse”の”pivot_longer”を使って縦型のデータへ変更します。
そのため、まずはインストールしていなければインストールをしていきます。
install.packages(“rstatix”)
install.packages(“tidyverse”)
それから、２つともlibraryを使えるようにします。
library(rstatix)
library(tidyverse)
データを縦型のデータにする
次に解析するにあたって縦型のデータにしていきます。まずはデータを。このデータはアイスの趣向調査のつもりで作成しました。idは人の名前、つまりAさん・Bさん・・・
次の項目はアイスの味を記載しており、バニラアイス・チョコアイス・抹茶アイスが”１：好き”・”0:嫌い”という形で作成しました。
データ名は”icePreference”とします。
| id | vanilla | chocolate | Matcha |
| A | 1 | 0 | 0 |
| B | 1 | 1 | 1 |
| C | 0 | 0 | 1 |
| D | 1 | 0 | 0 |
| E | 0 | 1 | 0 |
| F | 1 | 0 | 0 |
これを縦型のデータにすると下記になります。
| id | taste | preference |
| A | vanilla | 1 |
| A | chocolate | 0 |
| A | Matcha | 0 |
| B | vanilla | 1 |
| B | chocolate | 1 |
| B | Matcha | 1 |
| C | vanilla | 0 |
| C | chocolate | 0 |
| C | Matcha | 1 |
| D | vanilla | 1 |
| D | chocolate | 0 |
| D | Matcha | 0 |
| E | vanilla | 0 |
| E | chocolate | 1 |
| E | Matcha | 0 |
| F | vanilla | 1 |
| F | chocolate | 0 |
| F | Matcha | 0 |
味はtaste・好みはpreferenceと変数名を作成しました。
最初から縦型のデータであればそのまま解析ですが、最初のような形であれば、変換の必要があります。
ここで使用するのが”pivot_longer”。書き方は下記。
l_icePreference <-icePreference %>% pivot_longer（cols = c（”vanilla”, “chocolate”,”Matcha”）,names_to = “taste” , value_to=”preference”）
上記を解説して書くと下記。アンダーラインが引いてあるところは、データに合わせて書き換えてください。
作成したデータの名前を指定 <- 元のデータを指定 %>% pivot_longer（cols = c（”まとめるデータ①”, “まとめるデータ②”,”まとめるデータ③”）,names_to = “まとめるデータの名前” , value_to=”結果のデータのまとめ”）
name_to・value_toはそれぞれの名前を指定しています。記載しなければ”name”と”value”となります。
これでできたデータを”l_icePreference”へ入れるように記載しました。本当は”icePreference”を上書きしてもいいのですが、これならミスした時に直せるので私はこうしています。
cochranQ検定を行う
縦型のデータに変更したら、次はcochranQ検定を行います。その方法は、cochran_qtestを使って行います。数式は下記。アンダーラインが引いてあるところは、データに合わせて書き換えてください。
結果 <- cochran_qtest(調べるデータ , 検索するグループ ~ 結果 |人や検体 )
実際に上記データで行ってみると下記になります
c_icePreference <- cochran_qtest(l_icePreference, taste ~ preference | id)
print(c_icePreference)
１列目でcochranQ検定を行い、結果をc_icePreferenceに入れる。print(c_icePreference)で表示します。
結果の各項目の見方は下記
.y.（Outcome変数）：
検証した結果件数の名前が表示されます
n（サンプル数）：
分析に含まれた独立した被験者（参加者）の総数です。
データに欠損値（NA）がある場合、その被験者は自動的に除外されてカウントが減ります。
statistic（検定統計量Q）
コクランのQ検定における計算統計量（Q値）です
この値が大きいほど各グループ感のh自立の偏りや差が大きいことを意味します
df（自由度）
「比較したグループ数-1」になります
p（p値）
p値になります。ここで有意差をみていきます
最後にMcNemer検定の多重検定をおこなう
cochranQ検定で有意差があれば、その中のどこかで差があることがわかりました。ですが、実際には差がどこにあるのかわかりません。
そこでMcNemer検定をすべての項目で行いどこに有意差があるか、確認していきます。この検定を事後検定（post-hoc test）といいます。
ただ、検定を何度も行うため有意差がないのに有意差があるとしてしまう可能性が高くなってしまいますので、それに補正を行っていきます。
数ある補正の中で、この条件で使用するのはBonferroni法とholm法があります。
今回の数式は下記。アンダーラインが引いてあるところは、データに合わせて書き換えてください。
結果 <- pairwise_mcnemer_test (調べるデータ , 検索するグループ ~ 結果 |人や検体 ,type=”exact” , p.adjust.method = 補正方法 )
上記のデータで記載してみると下記
p_icePreference <- pairwise_mcnemer_test(l_icePreference, taste ~ preference | id , p.adjust.method = “holm”)
多重解析を１行で記載できると言うのが”R”のすごいところ。しかも、cochranQ検定の書き方とほぼ一緒という素晴らしい計算式です。
type=”exact”の記載はMcNemer検定の正確検定（Exact test）を行う場合に記載します。通常のMcNemer検定を行うのであれば、削除して構いません。
holm法で行っていますが、holm→bonferroniに記載を変えることでbonferron法に変わります
結果の見方は下記
- group1 / group2（比較対象のグループ）
- 現在どのグループとどのグループを比較しているかを示します。
- 例えば、1行目は「グループA」と「グループB」の比率を直接比較しています。
- p（修正前のp値）
- 多重比較の補正を行う前の、単発のマクネマー検定のp値です。
- 基本的にはこの値は無視して構いません。
- p.adj（修正後のp値）
- 最も重要チェック項目です。
- 指定した補正方法（holmやbonferroni）によって、グループを何度も比較したことによるエラーを修正したp値です。
- p.adj.signif（有意性のマーク）
- p.adj の結果を視覚的に分かりやすく星（アスタリスク）などで表したものです。
- ns: 有意差なし (not significant, p.adj ≥ 0.05)
- *: 有意差あり (p.adj < 0.05)
- **: 強い有意差あり (p.adj < 0.01)
- ***: 非常に強い有意差あり (p.adj < 0.001)
- method（使用した検定手法）
- 実行された検定名が表示されます（通常は McNemar test、または引数で正確確率を指定した場合は Exact binomial test となります）。
これで一通り解析できました。
McNemar準正確検定（Mid-p検定）については多重解析の方法はわかりませんでした。行うのであれば１個づつ検定して補正（holmやbonferroni）をかけることになると思います。
コメント