2 ノンパラ対応なし群間差

今回から群間の差の検定とそのグラフです。

一元配置のデータで、個体間要因とか級間要因とか betwee-factor とか言われる場合です。

まずはノンパラから。

２群の比較と3群以上の比較は統計のかけ方が異なるのですが、、、

グラフを描く上では2群もそれより多い場合も何も変わらないので、同じページで解説しています。

2群の比較の検定は、サンプルデータから2群だけ抽出する方法の解説とともにやります（結局、そのやり方をしないとできないので...）。

ノンパラは、平均値の比較ではななくて、数値データを順位に変換して、順位が高いデータがどの群に多いか？ということを比較していると思うので、平均値のグラフを使わない方が妥当だと思っています。なので、分布を可視化する表現方法の箱ひげ図（boxplot）を僕は使うことにしています。

箱ひげの、箱がその群の順位の1/4区切りを示していて、箱の内部に収まるデータがその群の半数で、真ん中の線が中央値です。上下の棒（これがヒゲ）の先端が最大値と最小値になることが多いですが、外れ値がある場合はヒゲの外にもdotsがあります。

https://scrapbox.io/files/638798a69072ae002237dd59.png

グラフの描画はそんなに難しくないのですが、実は、群間差の場合のノンパラ検定は、不等分散の場合に使えない（第1種の過誤がメチャ起る）ので、紹介はするものの、あまり使わないで頂きたい感じ。

詳しくははじめにこのサイトについてのした方を再度見てください。

code:r

# ---------- preparation in advanca 事前準備--------------------------------------------------

# Install packages if not installed 必要なパッケージが無けれな自動でインストール-----

if(!require("ggplot2")) install.packages("ggplot2")

if(!require("dplyr")) install.packages("dplyr") #データ整形に使う（filterとか）

if(!require("NSM3")) install.packages("NSM3") #先に直接.tar.gzの圧縮fileをCRANからDLしてRstudioのPackgesウィンドウからインストール

#if(!require("partitions")) install.packages("partitions") #先に直接.tar.gzの圧縮fileをCRANからDLしてRstudioのPackgesウィンドウからインストール

#if(!require("gmp")) install.packages("gmp") #先に直接.tar.gzの圧縮fileをCRANからDLしてRstudioのPackgesウィンドウからインストール

if(!require("exactRankTests")) install.packages("exactRankTests")　#needed fo Exact Rank Tests

if(!require("lawstat")) install.packages("lawstat") #needed for Brunner-Munzel test

# clear R's brain 一度参照したリストもしくは変数を次の作業で消す-----

rm(list = ls())

# Load libraries 必要なパッケージの読み込み-----

library(ggplot2)

library(dplyr) #データ整形に使う（filterとか）

library(NSM3) #needed for Steel-Dwass test

#library(partitions) 　#needed for NSM3, maybe

#library(gmp) 　　　　 #needed for NSM3, maybe

library(exactRankTests) #needed fo Exact Rank Tests

library(lawstat) #needed for Brunner-Munzel test

# Data input 今後以下のパスからデータを自動で読み込みする場合-----

library(readr)

CanRsample <- read.csv(file.choose()) #Run and choose CanRsample2021.csv

#ここを実行（Run）するとfile選択画面になるので

# 【2】box plots between　箱ひげ図個体間----------------------------

kruskal.test(CanRsample$aveTemp_year ~ CanRsample$place) #(データシート$従属変数 ~ データシート$独立変数)

#下位検定の多重比較

pairwise.wilcox.test(CanRsample$aveTemp_year, CanRsample$place, paired=F, correct=F, exact=F, p.adj="holm")

#post-hoc multiple comparison, correct=F で連続性補正をしないようにしている、Rのデフォ設定でn50以上になると正規近似されるが、

##タイ値があってP値が計算されない場合があるため、正規近似するために exact=F にしている

pSDCFlig(CanRsample$aveTemp_year, as.factor(CanRsample$place), method = "Asymptotic")

#Steel-Dwass検定によるpost-hoc、漸近法（Asymptotic）の正規近似でやってます。

##n10以下ならmethod = "Exact" の方が良いけど計算重くてスタックするかも（頑張ってn11以上とりましょう）。

#ノンパラ2群比較したい場合---------

#各群だけのデータセット作成（オブジェクト化） genaration of data sets divided into each group

dSendai <-filter(CanRsample, place == "Sendai")　#(元データ, 抽出したい群名がある列名 == "群名")

dTokyo <-filter(CanRsample, place == "Tokyo")

dKagoshima <-filter(CanRsample, place == "Kagoshima")

wilcox.test(dSendai$aveTemp_year, dKagoshima$aveTemp_year, paired=F, exact=F, correct=F)　

#タイが無いとexact=Tにしても下の関数と同じ結果

wilcox.exact(dSendai$aveTemp_year, dKagoshima$aveTemp_year, paired=F)

#正規性や等分散性を仮定しないノンパラ2群検定

brunner.munzel.test(dSendai$aveTemp_year, dKagoshima$aveTemp_year) #今回のデータだと計算が収束しないぽい

#---------Figづくり--------------

#3群以上を前提にしてるので2群比較の図にしたいときはデータ整形したりオブジェクトを作ってください

#Fig_boxplot <-

ggplot(CanRsample, aes(x = place, y = aveTemp_year, group = place)) +

#データセットと使用する変数の指定、箱の色も帰る場合はaesに "color = グループ変数" を追加

scale_x_discrete(limit=c('Sendai', 'Tokyo', 'Kagoshima')) +

#X軸に表示したい順番の指定、不要なら消す

geom_boxplot(size = 1.5, width = 0.4) + #sizeは線の太さ、widthはboxの太さ

geom_jitter(shape = 21, size = 4, aes(fill = place), width = 0.2, alpha = 0.7) +

#geom_point(size = 4, colour = 'lightgrey', alpha = 0.7) + #alphaはdotの透明度

scale_fill_manual(values = c (Sendai = "deepskyblue3", Tokyo = "gray1", Kagoshima = "orange" )) +

#色をRで自動指定したい場合上記関数を消す

xlab("場所") + #X軸のラベル

ylab("年平均気温") + #Y軸のラベル

#ylim(max, min) + #Y軸の（最大値, 最小値）を指定する場合は記入して先頭の#を消す

ggtitle("【2】box plots between") +

#coord_flip() +　#X軸とY軸を逆にしたいときに使う

theme_classic(base_size = 22, base_family = "HiraKakuProN-W3") # white backgrond グラフの背景を白に

#theme(legend.position = "none") #凡例を消したい時は#を消して上のcodeの最後に+を足して有効にする

【統計の手法の確認】

パラメトリックな検定をする際には、

one-way ANOVA → Tukey-Kramer HSD 検定

という流れが用いられます。

ノンパラで3群以上の差を検定する際によく使われるのは、

Kruskal-Wallis test → Steel-Dwass test

という手順です。

Kruskal-Wallis test は、原理的には、マン・ホイットニーのU検定と同様だそうです。ちなみに、U検定とウィルコクソンの順位和検定（Wilcoxon rank sum test）は、同じ結果をもたらすものであると知られています。で、はじめにこのサイトについての最後の方で書いた通り（と上でも行った通り）、不等分散に非常に弱い。

ただ、このページでは（かつての僕がせっかく調べたので）一通りやっていきます。

Steel-Dwassをしてくれる pSDCFlig という関数が入っている NSM3 というパッケージをこのページではインストールしています。僕が初めてNSM3を使った時は、さらに、partitions と gmp というパッケージがないと動かないとRに言われたのでインストールしたのですが、今回、これらが無くても動きました。一応、誰かのためにコメントとして残してあります。

さて、統計を進める手順は上記の通りなので、上から実行してもらえればOKです。

後検定の多重比較としては、Steel-Dwass検定以外にも、ウィルコクソンの順位和検定を総当たり（ペアワイズ）にやった上でHolm法でP値を調整をする方法も載せています。ウィルコクソンの用法（オプションの指定方法）も、スクリプト中のコメントを読んでもらえれば分かるかと思います。

【2群比較をする場合】

2群の比較の検定をする前に、dplyrパッケージのfilterを使って、群（Sendai、Tokyo、Kagoshima）ごとのデータシートのオブジェクトを作っています。なんでこんなことしてるのかと言うと、こういう指定の仕方をしないと、検定をしてくれない関数だからです。初めて使った時は、ちょっと信じられませんでした。2群だけが含まれるデータシートのオブジェクトだったとしても、ダメです。

検定ごと（関数ごと）に群・条件（独立変数）の指定の仕方が全然違うので、最初はとまどったいうか、調べるのに時間がかかったものです。

上記 script 中では、SendaiとKagoshimaを比較する場合の書き方だけ示しています。

Brunner-Munzel 検定についてだけ、最後に。

これは、U検定と違って、等分散性を仮定しなくて良い検定です。むっちゃ便利ですね。使える時は使うと良いと思います。ただ、たまに計算が収束しません。それと、3群以上の場合に使える似たような検定も、たぶんまだ知られてないと思います。そんなにメジャーじゃないので、生物系の雑誌で使った際に、Reviewerが知らない可能性があるので、ちょっと怖い。心理系では、だいぶ知られているのではないでしょうか。

【グラフ描画】

すでに 1_1 と 1_2 を読んでくださった方なら、簡単だと思うので、それら前のページには書いてない特徴だけ説明します。

今回は、boxplot() という関数で描きます（その名の通りですね）。

横軸の並び

散布図じゃなくて群間の差なので、横軸にはカテゴリカル変数（群の名前）が並びます。

その並びを自分の好みの順番に指定するための

scale_x_discrete()

というものを入れてあります。不要なら消してください。

dotの表示：散らすか一直線にするか

geom_jitter()：dotを散らしたいときに使います（反復測定じゃないならこれが良いと思います、特に、nがでかいほどに）

geom_point()：dotを散らさず一直線にするときに使います（反復測定で同一個体とかに線を引くならこっちが良い）

このどちらか1つを使います。

dotの色指定

グラフの場合、色は scale_fill_manual() で決めます。

scale_color_manual() ではない！という点に注意。