4 対応のないパラメトリック群間差

一元配置分散分析です。

はじめにこのサイトについてと2 ノンパラ対応なし群間差を見ていただければと思いますが、基本的に、一元配置の分散分析は通常の one-way ANOVAではなく、Welch's one-way ANOVA をこのサイトでは使います。同様に、t検定もWelchを問答無用で使います。実は、Rでは、one-way ANOVAもt検定も、オプションの指定を特に書かなければ、デフォルトでWelchの検定を行いますので、このページのスクリプトでも等分散性を仮定するか否かの指定は書いてません。

https://scrapbox.io/files/63879c1a2b6e340020865a45.png

こんなグラフを描きます。まずは、スクリプト全体を以下に載せます。

code:r

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

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

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

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

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

if(!require("effectsize")) install.packages("effectsize", repos = "https://easystats.r-universe.dev")

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

rm(list = ls())

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

library(dplyr) #データ整形に使う

library(ggplot2)

library(ggfortify) #needed for autoplot()

library(effectsize)

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

library(readr)

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

#ここを実行（Run）するとfile選択画面になるので、どのfileを読み込んだか上に記録！

# # ----------prep done 準備だん------------------------------------------------------------

#【4】comparison of means between　平均値の個体間比較----------------------------

#データの分布など確認（診断プロット）

aveTempY.lm <-lm(aveTemp_year ~ place, data = CanRsample)

autoplot(aveTempY.lm, smooth.colour = NA) #ANOVAとか線形モデルを使う前提の診断プロット

#Welch's one-way ANOVA

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

#post-hoc multiple comparison （Welch用にプーリングしないSDを用いる pool.sd=F を指定）

pairwise.t.test(CanRsample$aveTemp_year, CanRsample$place, p.adj="holm", pool.sd=F)

#2群比較したい場合

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

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

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

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

t.test(dSendai$aveTemp_year, dKagoshima$aveTemp_year) #Welch's t-test

#特定の2群間でのみ比較したければ

dKago_Sen <-filter(CanRsample, place == "Kagoshima" | place == "Sendai") 　#2群をフィルターで選んでオブジェクト化

dKago_Tok <-filter(CanRsample, place == "Kagoshima" | place == "Tokyo")

dSen_Tok <-filter(CanRsample, place == "Sendai" | place == "Tokyo")

#エフェクトサイズの計算

cohens_d(aveTemp_year ~ place, data = dKago_Sen)

cohens_d(aveTemp_year ~ place, data = dKago_Tok)

cohens_d(aveTemp_year ~ place, data = dSen_Tok)

#Fig_mean <- #グリッドして並べる際はこのコメントを消すことでオブジェクト化

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

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

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

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

#ランダムなジッタープロット、position関数を使うとバグる、aesのfillにドットの色分け用グループ名を入れる

#geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8), aes(fill = 色分けグループ名), alpha = 0.7) +

#Prismのような揃ったプロット、aesのfillにドットの色分け用グループ名を入れる

#geom_point(size = 4, aes(colour = グループ変数), alpha = 0.7) + #alphaはdotの透明度、縦一直線のdot

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

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

stat_summary(fun.y = mean, geom = "crossbar", colour = "black", width = 0.4) + #棒グラフにしたい場合はcrossbarをbarに変える

stat_summary(fun.data = mean_se, geom = "errorbar", colour = "black", width = 0.2) + #widthは横線の長さ

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

ggtitle("【4】comparison of means between　平均値の個体間比較") +

#coord_flip() + #XとYの入れ替えに使う

theme_classic(base_size = 22, base_family = "HiraKakuProN-W3")

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

【統計】

これはもう、「スクリプトを見て下さい」というダケで大丈夫でしょうか。変数指定の方法も簡単です。

最初に、一応分布の確認をするプロットを作るスクリプト書いています（1_2 群ごとの相関と共分散分析参照）。

その後、Welchの分散分析→後検定という流れです。

その後には、一応、2群間のt検定をするためのスクリプトを載せています。

実は、普通のt検定をする方がデータの指定方法がめんどくさくて、

t.test(群1のデータだけのオブジェクト$従属変数, 群2のデータだけのオブジェクト$従属変数)

という指定の仕方をしなければなりません。

なので、この検定の前に、各群のデータだけを抽出する filter を使っている部分がありますね。

※ ただ、よくよく考えれば、2群に対して one-way ANOVA をすれば、結果は t検定と同一なので、

oneway.test(データセットオブジェクト名 ~ group, data = データセットオブジェクト名, var.equal = FALSE)

のように、Welchのt検定をかけてしまえば楽ですね。

さらに、t検定の効果量（Cohen's D）を計算する構文を最後につけています。

これは、2つの群（だけ）を含んでいるデータのオブジェクトを使わねばなりません。。。一緒に使うものはデータ指定形式を揃えて欲しいと思ってしまいますが、きっとパッケージや関数の作者さんが違うので、しょうがないですね。スクリプトの書き方で工夫するほかありません。

＊最近、効果量を書くことが推奨される機会が増えていますが、検定したら自動でセットで出してくれるソフトや関数は少ない印象です。個人的には、SEMのエラーバーが群間・条件間で被っていないとか、箱ひげ図の箱が被っていないとか、dotplotを見れば個体分布がかなり離れているとか、そういうのを見る方が実は重要じゃないかと思っています。

※ 後検定について補足

Tukey-Kramer HSD 検定は通常のone-way ANOVAのあとの多重比較の検定として有名ですが、実は群間の等分散性を前提にしていて、不等分散の場合に使うのは不適切です。不等分散のときにも使える方法としてはGames-Howell法というものがあるそうなんですが、Games-Howell法も含む色んな多重比較をやってくれる PMCMRplusパッケージが、僕は全然うまくインストールできません。

PMCMRplusパッケージ https://rdrr.io/cran/PMCMRplus/man/gamesHowellTest.html

このサイトでは、Welch's t-test をペアワイズに行って、Holm法でP値を補正する方法を採用しています。HADで Welch's ANOVAをして後検定までセットで行った場合と値が一致することを個人的に確認はしています。

【グラフの描画】

これまでに前の章で扱ってきたやり方とほとんど同じです。特に2 ノンパラ対応なし群間差に近いです。平均値の比較なので、中央値を表示させる箱ひげ図の代わりに、平均値を表示させる感じです。

geom_boxplotの代わりに、geom_barを使えば、いわゆる棒グラフが描けます。ただ、棒そのものの必要性をあまり感じないので、僕は普段平均値±SEMだけで描いています。

そのために使っているのが、stat_summary()という関数で、以下のようにこの関数を2回使って平均値±SEMを表示させています。

code:r

stat_summary(fun.y = mean, geom = "crossbar", colour = "black", width = 0.4) + #棒グラフにしたい場合はcrossbarをbarに変える

stat_summary(fun.data = mean_se, geom = "errorbar", colour = "black", width = 0.2) + #widthは横線の長さ

その他のポイントなどは、2 ノンパラ対応なし群間差を見てください。

＊バイオリンプロット

棒グラフにするくらいだったら、バイオリンプロットにするのが良いのではないかと思います。

以下では、冒頭に記したスクリプトに対し、ggplot2() の中の aes() に fill = place を追記して、その後に geom_violin() を足しています。データが膨大すぎる場合は、geom_jitter() を消してdotを表示させないという手もあるでしょう（もしくはdotはめっちゃ小さくする）。

https://scrapbox.io/files/65e61541c7f0450024d1cf8e.png

code:r

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

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

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

geom_violin(alpha = 0.2) +

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

#ランダムなジッタープロット、position関数を使うとバグる、aesのfillにドットの色分け用グループ名を入れる

#geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8), aes(fill = 色分けグループ名), alpha = 0.7) +

#Prismのような揃ったプロット、aesのfillにドットの色分け用グループ名を入れる

#geom_point(size = 4, aes(colour = グループ変数), alpha = 0.7) + #alphaはdotの透明度、縦一直線のdot

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

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

stat_summary(fun.y = mean, geom = "crossbar", colour = "black", width = 0.4) + #棒グラフにしたい場合はcrossbarをbarに変える

stat_summary(fun.data = mean_se, geom = "errorbar", colour = "black", width = 0.2) + #widthは横線の長さ

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

ggtitle("【4】comparison of means between　平均値の個体間比較") +

#coord_flip() + #XとYの入れ替えに使う

theme_classic(base_size = 22, base_family = "HiraKakuProN-W3")

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