6 Two-way ANOVA のグラフ

#basic の最後です。二元配置のグラフの書き方。これは、データの多さや水準の多さ、データの挙動によっても見やすさが変わってくるので、色んな描き方を掲載しています。

このページは、グラフが多いので、統計部分とグラフ部分のscriptを分けて載せておきます。

【統計】

今回も、5 対応ありパラメトリック群間差に引き続き、ANOVA君を使わせてもらっています。

パッケージのインストールから統計検定までのスクリプトは以下です。

code:r

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

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

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

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ウィンドウからインストール

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

rm(list = ls())

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

library(tidyverse)

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

library(ggplot2)

library(ggfortify) #needed for autoplot()

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

library(readr)

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

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

#横長データをtidyデータに変換したものも予め作っておく（多重比較や作図に便利）

CanRsample_tidy <- gather(CanRsample, key = "Season", value = "aveTemp", aveTemp_SpAut, aveTemp_summer, aveTemp_winter)

#keyは新しく作る列名（カテゴリカル変数名）、valueは値の新しい列名、その後に続くのが縦長に並べ替えたい列名

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

#【6】Line Graph　主に個体内比較用の折れ線グラフ-----------

source(file.choose())　#"anovakun_489.txt" を読み込み

#ワイド型データでやる場合（個体間要因の水準の命名順序を間違えそうなのでお勧めしない）

anovakun(dplyr::select(CanRsample, place, aveTemp_SpAut, aveTemp_summer, aveTemp_winter), "AsB", Place = c("Sendai", "Tokyo", "Kagoshima"), Season = c("SpAut", "Summer", "Winter"), holm = T, peta = T)

#後検定をウェルチ的なものでやる

#ロング型データでやる場合（お勧め）

anovakun(dplyr::select(CanRsample_tidy, ID_town, place, Season, aveTemp), "AsB", long = T, holm = T, peta = T)

#後検定をウェルチ的なものでやる

anovakun(dplyr::select(CanRsample_tidy, ID_town, place, Season, aveTemp), "AsB", long = T, holm = T, peta = T, welch = T)

aveTemp.mod <-lm(aveTemp ~ place*Season, data = CanRsample_tidy)

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

anova(aveTemp.mod) #交互作用などの結果を見るにはlm()の内容をanova()に渡す

summary(aveTemp.mod) #詳しいデータも知りたいときはこれもやる

ANOVA君の使い方は1つ前の 5 対応ありパラメトリック群間差のところでやったので、上記の通りにやって下さいという感じの説明になってしまいますが、ANOVAの形式の復習だけ書いておきます。ANOVA君では、使用するデータのオブジェクトの次に、

"個体間要因s個体内要因"

という書き方で要因のタイプと構造を指定します。前回は、繰り返しのある一元配置分散分析だったので "sA" という書き方で指定しました。

今回は、個体間要因1つ x 個体内要因1つの two-way なので、

"AsB"

と指定しています。

もう1つ、今回は、後検定のためのオプションについてだけ、補足します（figshare で共有してる script ではこのオプションを書いてません）。似たような anovakun() のコードを2つずつ、書いています。シャープで「後検定をウェルチ的なものでやる」と記したコメントのコードでは、"Welch = T" というものが追加されています。これは、後検定の多重比較を、Welch的な方法で行うオプションです。

※ figshare で共有してる script ではこのオプションを書いてません

ANOVA君このページの http://riseki.php.xdomain.jp/index.php?ANOVA君/ANOVA君の使い方、

多重比較の8に、以下のようにあります。

welch・・・「welch = T」とすると，多重比較の際にKeselman-Keselman-Shafferの統計量とWelch-Satterthwaiteの近似自由度を使用します。Welchの方法を複数の要因を含むデザインや反復測定要因を含むデザインに拡張した方法に基づいており，等分散性や多標本球面性が成立しない場合により適切な手法です（Keselman et al., 1991）。

Keselman, H. J., Keselman, J. C., & Shaffer, J. P. (1991). Multiple pairwise comparisons of repeated measures means under violation of multisample sphericity. Psychological Bulletin, 110(1), 162–170. https://doi.org/10.1037/0033-2909.110.1.162

このサイトでは基本的に、個体間要因の場合は、Welchの検定を激推ししてきました。Studentのt検定よりWelchのt検定、通常の one-way ANOVA よりも Welch's ANOVA をお勧めしています（詳しくははじめにこのサイトについて）。

であれば、two-way ANOVA の後検定でも、不等分散に強い方法で多重比較したくなるわけですが、ANOVA君はそれを Welch = T のオプションでやってくれるということで、僕も今後は使っていこうと思っています。

Welch = T が書いてあるもとと書いてないものの結果を試しに比較してみて下さい。

個体間要因（ここでは、Sendai・Tokyo・Kagoshima の違い、すなわち place の効果）の結果だけ異なっているはずです。

個体内要因の方は、普通にANOVAをする際も、球面性（等分散性みたいなもの）の仮定から逸脱していた場合に調整をする方法が市販ソフトとかでも実装されており、ANOVA君でも良い感じにやってくれるます。

Welch の検定がそもそも個体間要因の比較の方法なので、今回も個体間要因における結果だけが変わってきます。

球面性についても、ANOVA君のページの説明が分かりやすいんじゃないかと思います（井関先生すごい）。

http://riseki.php.xdomain.jp/index.php?ANOVA君/球面性検定の出力

【グラフの描画】

基本姿勢としては、反復測定の項目についてはライングラフにしたい、というものです。

https://scrapbox.io/files/6387a03678814900231fe9ad.png

code:r

#lineGraph <-

ggplot(CanRsample_tidy, aes(x = Season, y = aveTemp, group = place, color = place)) +

geom_line(stat = "summary", fun = "mean", size = 1.0, alpha = 0.5)+

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

geom_point(size = 4, aes(colour = place), alpha = 0.7) + #alphaはdotの透明度、縦一直線のdot

#dotを消したい場合は上の関数を#で無効に

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

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

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

#色を自由に指定したい場合、上記先頭の#を取って色を指定

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

#色を自由に指定した際に凡例の色も揃える、上記先頭の#を取って色を指定

xlab("季節") + #X軸のラベル

ylab("季節ごとの平均気温") + #Y軸のラベル

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

ggtitle("【6_1】Line Graph　") +

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

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

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

僕自身は、あんまりこの描き方は使わない方針なのですが、同一個体には線を引かず、平均値を線で結ぶバージョンです。個体数や群の数が多すぎると、同一個体に引いた線が見た目に邪魔になってしまうことがあり、その場合、この描き方が良いかなと思います。

で、実は、この描き方の時に使う独特なコードがあり、

geom_line() の中に、

stat = "summary", fun = "mean"

というものを書いています。これにより、平均値を線で結んでいます。

https://scrapbox.io/files/6387a0459072ae00223826d4.png

code:r

#【6_2】Line Graph　主に個体内比較用の折れ線グラフ（個体ごとに線を引くver1）-----------

#lineGraph <-

ggplot(CanRsample_tidy, aes(x = Season, y = aveTemp, group = place, color = place)) +

geom_line(aes(group = ID_town, colour = place), alpha = 0.5)+

geom_point(size = 4, aes(colour = place), alpha = 0.7) + #alphaはdotの透明度、縦一直線のdot

#dotを消したい場合は上の関数を#で無効に

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

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

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

#色を自由に指定したい場合、上記先頭の#を取って色を指定

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

#色を自由に指定した際に凡例の色も揃える、上記先頭の#を取って色を指定

xlab("季節") + #X軸のラベル

ylab("季節ごとの平均気温") + #Y軸のラベル

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

ggtitle("【6_2】Line Graph　") +

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

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

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

あんまり解説はいらないですね。3 ノンパラ対応ありのときと比べると

geom_boxplot() が不要で、代わりに stat_summary() を2種つかっているというだけです。

この描き方は、交互作用があったときには、他の群と違う挙動を示す群の様子が分かりやすくて好きです。

https://scrapbox.io/files/6387a0512aaa410020a0a5c4.png

code:r

#【6_3】Line Graph　主に個体内比較用の折れ線グラフ（個体ごとに線を引くver2）-----------

#lineGraph <-

ggplot(CanRsample_tidy, aes(x = Season, y = aveTemp, group = place, color = place)) +

geom_line(aes(group = ID_town, colour = place), alpha = 0.5)+

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

geom_point(size = 4, aes(colour = place), alpha = 0.7) + #alphaはdotの透明度、縦一直線のdot

#dotを消したい場合は上の関数を#で無効に

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

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

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

#色を自由に指定したい場合、上記先頭の#を取って色を指定

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

#色を自由に指定した際に凡例の色も揃える、上記先頭の#を取って色を指定

xlab("季節") + #X軸のラベル

ylab("季節ごとの平均気温") + #Y軸のラベル

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

ggtitle("【6_3】Line Graph　") +

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

facet_wrap(~place) + #同じ操作を群ごとや個体ごとにする場合、場合わけ変数を~の後に指定

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

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

ポイントは、下の方に書いてある、

facet_wrap(~place)

です。~の後に書いてある変数について、水準ごとに分けて同じ作業を繰り返す指定です（便利）。

これを使って、反復あり一元配置を群の数だけ描くというやり方で、二元配置の図にしてます。

facet_wrap() については、1_2 群ごとの相関と共分散分析の【1_3】でヒストグラムを書く時にやりました。

https://scrapbox.io/files/6387a05d0f561e001dc536be.png

code:r

#【6_4】Line Graph　主に個体内比較用の折れ線グラフ（個体ごとに線を引くver2）-----------

#lineGraph <-

ggplot(CanRsample_tidy, aes(x = place, y = aveTemp, color = place)) + #9_3からxをplaceに変えたよ

#geom_line(aes(group = ID_town, colour = place), alpha = 0.5)+ #このケースでは機能しない

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

geom_point(size = 4, aes(colour = place), alpha = 0.5) + #alphaはdotの透明度、縦一直線のdot

#dotを消したい場合は上の関数を#で無効に

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

stat_summary(fun.data = mean_se, geom = "errorbar", aes(fill = place), width = 0.4) + #widthは横線の長さ

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

#色を自由に指定したい場合、上記先頭の#を取って色を指定

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

#色を自由に指定した際に凡例の色も揃える、上記先頭の#を取って色を指定

xlab("季節") + #X軸のラベル

ylab("季節ごとの平均気温") + #Y軸のラベル

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

ggtitle("【6_4】Line Graph（じゃないけど...）　") +

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

facet_wrap(~Season) + #同じ操作を群ごとや個体ごとにする場合、場合わけ変数を~の後に指定

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

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

最後は、個体間 x 個体内の描き方ではなく、個体間 x 個体間の two-way の図ですが、他のページでやってないけど（今後もページを追加するつもりもない）必要なときもあるので、やりました。

個体間 x 個体間の two-way は、不等分散に弱そうな気がして、そもそもこれを使わないといけなくなる実験計画をしないように僕はしているのですが、ANOVA君で Welch = T のオプションを使っておけば、だいぶ安全なのかなと思います。