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")
# clear R's brain 一度参照したリストもしくは変数を次の作業で消す-----
rm(list = ls())
# Load libraries 必要なパッケージの読み込み-----
library(ggplot2)
library(ggfortify) #needed for autoplot() library(effectsize)
# Data input 今後以下のパスからデータを自動で読み込みする場合-----
library(readr)
CanRsample <- read.csv(file.choose()) #Run and choose CanRsample2021.csv # # ----------prep done 準備だん------------------------------------------------------------
aveTempY.lm <-lm(aveTemp_year ~ place, data = CanRsample)
oneway.test(CanRsample$aveTemp_year ~ CanRsample$place) #(データシート$従属変数 ~ データシート$独立変数) #post-hoc multiple comparison (Welch用にプーリングしたSDを用いる pool.sd=T を指定) pairwise.t.test(CanRsample$aveTemp_year, CanRsample$place, p.adj="holm", pool.sd=T)
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 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)
ggplot(CanRsample, aes(x = place, y = aveTemp_year)) +
scale_x_discrete(limit=c('Sendai', 'Tokyo', 'Kagoshima')) +
geom_jitter(shape = 21, size = 4, aes(fill = place), width = 0.2, alpha = 0.7) +
scale_fill_manual(values = c (Sendai = "deepskyblue3", Tokyo = "gray1", Kagoshima = "orange" )) +
stat_summary(fun.data = mean_se, geom = "errorbar", colour = "black", width = 0.2) + #widthは横線の長さ ggtitle("【4】comparison of means between 平均値の個体間比較") +
theme_classic(base_size = 22, base_family = "HiraKakuProN-W3")
【統計】
これはもう、「スクリプトを見て下さい」というダケで大丈夫でしょうか。変数指定の方法も簡単です。
その後、Welchの分散分析→後検定という流れです。
その後には、一応、2群間のt検定をするためのスクリプトを載せています。
実は、普通のt検定をする方がデータの指定方法がめんどくさくて、
t.test(群1のデータだけのオブジェクト$従属変数, 群2のデータだけのオブジェクト$従属変数)
という指定の仕方をしなければなりません。
なので、この検定の前に、各群のデータだけを抽出する filter を使っている部分がありますね。
さらに、t検定の効果量(Cohen's D)を計算する構文を最後につけています。
これは、2つの群(だけ)を含んでいるデータのオブジェクトを使わねばなりません。。。一緒に使うものはデータ指定形式を揃えて欲しいと思ってしまいますが、きっとパッケージや関数の作者さんが違うので、しょうがないですね。スクリプトの書き方で工夫するほかありません。
*最近、効果量を書くことが推奨される機会が増えていますが、検定したら自動でセットで出してくれるソフトや関数は少ない印象です。個人的には、SEMのエラーバーが群間・条件間で被っていないとか、箱ひげ図の箱が被っていないとか、dotplotを見れば個体分布がかなり離れているとか、そういうのを見る方が実は重要じゃないかと思っています。
※ 後検定について補足
Tukey-Kramer HSD 検定は通常のone-way ANOVAのあとの多重比較の検定として有名ですが、実は群間の等分散性を前提にしていて、不等分散の場合に使うのは不適切です。不等分散のときにも使える方法としてはGames-Howell法というものがあるそうなんですが、Games-Howell法も含む色んな多重比較をやってくれる PMCMRplusパッケージが、僕は全然うまくインストールできません。
このサイトでは、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.data = mean_se, geom = "errorbar", colour = "black", width = 0.2) + #widthは横線の長さ * バイオリンプロット
棒グラフにするくらいだったら、バイオリンプロットにするのが良いのではないかと思います。
以下では、冒頭に記したスクリプトに対し、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')) +
geom_violin(alpha = 0.2) +
geom_jitter(shape = 21, size = 4, aes(fill = place), width = 0.2, alpha = 0.7) +
scale_fill_manual(values = c (Sendai = "deepskyblue3", Tokyo = "gray1", Kagoshima = "orange" )) +
stat_summary(fun.data = mean_se, geom = "errorbar", colour = "black", width = 0.2) + #widthは横線の長さ ggtitle("【4】comparison of means between 平均値の個体間比較") +
theme_classic(base_size = 22, base_family = "HiraKakuProN-W3")