1_2 群ごとの相関と共分散分析
1_1では、dotに色分けはしてましたが、その色分けの基となる群も全部プールして、1つの集団としての相関をみるような図でした。今回は、複数の集団の集団内での相関を個別にみながら1つの図にまとめるような描画をやります。
https://scrapbox.io/files/63879263839d30001d26ad9e.png
最終的には、こんなのを描いてみましょう。
この図は、複数の図を組み合わせる関数を使っているので、まずはそれぞれの図を作るところから。
全体のscriptを先に載せておきます。
code:r
# ---------- preparation in advanca 事前準備--------------------------------------------------
# Install packages if not installed 必要なパッケージが無けれな自動でインストール-----
if(!require("ggplot2")) install.packages("ggplot2")
# clear R's brain 一度参照したリストもしくは変数を次の作業で消す-----
rm(list = ls())
# Load libraries 必要なパッケージの読み込み-----
library(ggplot2)
library(ggfortify) #needed for autoplot() # Data input 今後以下のパスからデータを自動で読み込みする場合-----
library(readr)
CanRsample <- read.csv(file.choose()) #Run and choose CanRsample2021.csv # # ----------prep done 準備だん------------------------------------------------------------
# 【1_2】scatter plots 散布図(相関など)----------------------------
temp.mod <-lm(aveTemp_summer ~ aveTemp_year*place, data = CanRsample)
Fig_smootheach <-
ggplot(CanRsample, aes(x = aveTemp_year, y = aveTemp_summer, color = place)) +
scale_color_manual(values = c (Sendai = "deepskyblue3", Tokyo = "gray1", Kagoshima = "orange" )) +
geom_smooth(method = 'lm', se = FALSE, mapping = aes(colour = place)) + #SEMを表示したいときは se = FALSE を消す ggtitle("【1_2】scatter plots 散布図(相関など)") +
theme_classic(base_size = 22, base_family = "HiraKakuProN-W3") #+ # white backgrond グラフの背景を白に ggplot(CanRsample, aes(x = aveTemp_year, colour = place, fill = place)) +
scale_color_manual(values = c (Sendai = "deepskyblue3", Tokyo = "gray1", Kagoshima = "orange" )) +
scale_fill_manual(values = c (Sendai = "deepskyblue3", Tokyo = "gray1", Kagoshima = "orange" )) +
ggtitle("【1_3】histogram ヒストグラム") +
theme_classic(base_size = 22, base_family = "HiraKakuProN-W3")
type = "histogram",
margins = "both",
size = 2, groupColour = TRUE, groupFill = TRUE,)
【複数集団の相関関係を1つのグラフに】
https://scrapbox.io/files/638792f5fc67870020a196ba.png
まず、図の説明から。Fig_smootheach というオブジェクトにした図は、こんな風に描かれるものです。
前回1_2では geom_smooth() をシャープでコメント化している方のscriptで実行しました。
前回との違いは上記スクリプト中で見比べて欲しいのですが、今回は、以下のように書いています。
code:r
geom_smooth(method = 'lm', se = FALSE, mapping = aes(colour = place))
mapping = aes(colour = place) を加えることで、placeを群として、群ごとに色を変えた回帰直線を引いています。
直線ごとにSEMのゾーンをグレーに塗ると、ちょっと目障りな感じだったので se = FALSE で消しています。
結構簡単に図が描けてしまいますね。
ここからは、図のスクリプトの前に書いてある統計の説明をちょっとします。難しいと感じる人は、別の本などで勉強してもらった方が良いと思います。前に紹介した、
でも解説されています。
各群ごとの相関をみるだけで良ければ、前回やったような普通の相関の分析をすればOKですが、各群の傾向の比較みたいなことをしたければ、共分散分析をするのが良いと思います。パラメトリックな検定なので、正規分布を前提としています。
正直、この解析をする場合は、「こんなscriptをで解析を指定してください」としか言いようがないように思うので、scriptのコメントを見て理解してください。ここでは、結果の理解のための要点を解説します。
今回入れているggfortify というパッケージで、autoplot() を実行出来ます。パラメトリック検定の前提となる正規性の確認ができます。正規性といっても、独立2群の比較の場合は各群が正規分布していることが前提なんですが、ANOVAで問題となってくる分布(正規性とか等分散性(球面性)、および球面性の補正)というのは、誤差の分布のことなので、計算するのが結構面倒ですね。その辺を、図解して確認させてくれるのが autoplot() です。まずは、実行結果の図を見ましょうか。
https://scrapbox.io/files/65e0ab1deae51900252ef902.png
僕も完全に理解しているわけではないんですが、、、dotがだいたい点線に対して線対象に均等に分布していればOK!みたいな理解です。外れ値と判定されたものに数字が振られています。誤差の正規性を表しているのは右上のQ-Qプロットで、斜めの直線の付近にdotが乗っかっていれば大丈夫で、端にいくほど逸脱しがちなデータが目立つようになります。今回のデータは、まぁ許容範囲だと思います。two-way ANOVAをするときとかも、一応見ておくと良いと思います。
で、共分散分析(ANCOVA)の結果ですが、anova(temp.mod) で出力されるものを見ます。
年間平均気温の効果が有意:全体的に年平均気温と夏の平均気温が相関する
placeの効果は有意ではない:3本の回帰直線をあてはめてますが、1_1でやったように1本引いて相関みるのと大差ない感じという理解
交互作用も有意ではない
交互作用が有意ではないというのは、どういうことでしょうか。「有意だった場合どうなるのか」を先にお伝えした方が早いと思います。交互作用が有意だった場合というのは、それぞれの線の傾きが異なるというときです。相関関係がある群とない群を比較した時や、どちらの群も相関があるが傾きが異なる場合(x軸の増加に対するy軸の値の応答性・反応性が異なる場合)に有意になると言えるでしょう(有意になるときに、そのように解釈できる)。
今回のサンプルデータではつまらない解析になってしまいますが、いろいろ使える解析方法だと多います。
いずれにせよ、複数の群の相関を1つの図で見ることができると、気持ちがいいですね。
【ヒストグラムを描いてみよう】
このページの最終目標は、散布図の傍にxとyの値のヒストグラムも描いてみよう!というわけで、このページの一番最初に掲げている図を作るのが目標です(統計の話は先ほどでもう終わり)。
で、その最終的な図を作るのには必要ではないのですが、ヒストグラムだけ書く時の方法もここに記しておきます。
---【1_3】histogram ヒストグラム--- 以下に書いてあるスクリプトは、年平均気温のヒストグラムを書くコードです。ggplot2特有のさまざまな設定は、前回1_1とだいたい同じですね。実行すると以下のようになります。
https://scrapbox.io/files/638793efc0d3f8001ff52e8c.png
いろんなところで使えるポイントは、
facet_wrap(~place)
というものです。この書き方だと、placeの水準ごとにそのスクリプトで指定した描画を個別にやってくれます。
試しに、facet_wrap(~place)+ の前にシャープをつけたり、この行を消したりしてから描画をすると、以下のようになります。
https://scrapbox.io/files/65e0b32594f65900242e1c91.png
目的に応じて使い分けてください。
【Marginal distribution を描いてみる】
いよいよ最終段階です。このページの冒頭の図は、Marginal distribution と言われているものだと思いますが、それを描くには、ggplot2 に加えて ggExtraというパッケージを使います。
冒頭に記したスクリプト通りにやってください、と言うだけなのですが、
R Graph Gallery に記載されているscriptに対して僕が工夫したのは、ヒストグラムも色分けして、「積み重ね」ではなく「重ね合わせ」にしたことです。コメント部分を読んで、自分でいろいろ設定変更して試すと、意味がわかると思います。
一応、注意点は、Fig_smootheach を先にオブジェクト化しておくことです。相関の図を先に描いて(書いて)オブジェクト化し、そのオブジェクトを ggExtra で読み込んで、傍にヒストグラムを追加で表示させています。
https://scrapbox.io/files/63879330a5e05f001d4e9966.png
スクリプトを実行すると、こんな図が出てくるはずです。凡例が邪魔ですね。その時は、Fig_smootheach の方で凡例を非表示にするように書き換えて(1_1 散布図と相関およびグラフ設定の基本 参照)、再度オブジェクト化しなおしてから、ggExtraのスクリプトを実行してください。 最後に、ちょっと遊んでみます。ggExtraでは、ヒストグラムの他に、バイオリンプロットのような密度かボックスプロットを描くことができます。
type = の部分を、
type = "boxplot"
type = "density"
などと、変更するだけで、以下のような図が描画出来ます。
https://scrapbox.io/files/65e0b84e9923e60024fcfc10.png
https://scrapbox.io/files/65e0b85bd22de50024ce1d06.png
こういう図が作れると、1つの図で相関も群間の差も表すことができるので、あとは、イラレなどで補足の文字などを入れれば、いろんな情報が一目で分かるを図を作ることができます。短報を投稿する際などにとても便利だと思います。
この他、
type = "violin"
type = "densigram" (densityとhistgramが重なったもの)
ちなみに、散布図に対して「ボックスプロットや密度のグラフがデカ過ぎだなぁ」と思ったら、
sizeの数値を上てください。上記の僕のスクリプトでは "size = 2" ですが、5にすると以下のように散布図の割合が増え(おそらく値はここを指定している)、相対的に傍のグラフが小さくなります。
https://scrapbox.io/files/65e0ba4801b0550023e14093.png
ちなみに、
分布もしくはboxplotで中央値などを表すだけじゃなく、平均値を表す方法もあれば良いなと思うのですが、簡単な方法はなさそうです。