公開日:

【R】t検定を行い、有意差を示すアスタリスク付きのグラフを作成する【改】

以前、学部生の頃に作成した記事である【R】2群間においてt検定を実行し,アスタリスク付きのグラフを表示する がこのサイトのアクセス数第4位を維持し続けている。(1位はカナヘビの飼い方、2位はなぜ動物の交尾はバックばかり?~魚の妊娠の起源から陸上四肢動物まで~、3位は【Python】scipy.statsを使って正規分布の累積分布関数をプロットする)

今回論文を作成するにあたってggplot2を用いて有意差を示すアスタリスクとバーを含むきれいなグラフを作りたい、ということから勉強して作成した。

出来上がった図はこちら。

前回はライブラリ(ggplot2)を使わずに作成したが、以下のようなものだった。

ダサかったね。現時点で3000アクセス以上あったけど、ごめんなさい。

RとR studioをインストール

RR studioをインストールしておこう。特にR studioはR notebookを使って「ノートを取るように」解析結果を残しておくことができる。コマンドライン上でちまちまやっていると、その実行したコマンドが記録に残らなかったりしてあとから図の修正を行う際に面倒なことになる。

スクリプトファイルにしておくというのも手だが、いちいち図を出力し直したりするのが面倒くさい。

今回はRのバージョン4以上を対象に行っていく。

必要なライブラリのインストール

> library(ggplot2)
> library(ggsignif)

でエラー無く実行できるなら既にインストールされているので問題ない。もしインストールされていないのであれば以下のコマンドを入力してインストールする必要がある

> install.packages("ggplot2")
> install.packages("ggsignif")

これで1つ上のコマンドを実行すれば、問題なく実行できるはずである(何も表示されない)

インストールは完了。

ggplot2 で箱ひげ図なり棒グラフなりを作ろう

まずは自分の作りたい図をイメージしよう。

先の図は箱ひげ図(Box plot)であったが、場合によっては棒グラフなどを使うこともあるだろう。

ある生物種ごとにケラチンに含まれるアミノ酸残基数をカウントし、それをプロットしてやりたかった。与えてやるデータセットは以下のような形式である。

Irisでもできるが、種数が少ない
raw_data <- read.table(file='sample.txt', sep='\t', header = T, fileEncoding="UTF-8")

num <- transform(raw_data, Species= factor(Species, levels = c('A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'R', 'J')))

library(ggplot2)
plt <- ggplot(num, aes(y=CysteineResidues,x=Species)) +
  geom_boxplot(width=0.4) + 
  ylab("NUMBER") + expand_limits(y = 0) 
plt

このようにすると、

こんな感じのものがプロットされる。これだけでもなんとなく見栄えのする図となる。ggplot2スゲー。こだわる人間はもっと色々調節していくがそれは後ほど。

箱ひげ図の読み方については今私と同年代の人間は高校の数学Aだったかでやったはずなので教科書を見直そう。

検定を行って有意差のアスタリスクと”バー”をつける

次にこれらの種間(A-J)の間で数に統計的に有意な差が存在するかを明らかにすべく、検定を行う。通常、Rではデフォルトでt検定などの検定関数が入っているため、検定だけならば追加のライブラリは必要ないものが多い。

今回はこのggplot2のグラフに重ねる形でアスタリスクやバーを追加したいので、ggsignifというライブラリを使う。

使い方は非常にかんたん。

library("ggsignif")
plt2 <- plt +   geom_signif(comparisons = list(c("A", "B")),
              test = "t.test",
              na.rm = FALSE,
              map_signif_level = TRUE,
              y_position = 18,
              col = "red") + theme_minimal()


plt2

これによって以下の図が得られる。

わかりやすいように赤で示しているが、見た目が悪いので通常は黒に近い色にしたほうがいいだろう。

ggsignifの詳細

詳しくはこちらに書かれている(英語)

comparisons

comparisons = list(c("A", "B"))の部分で、比べる群を指定している。リスト形式にして渡す必要がある。

test

test = "t.test"で行う検定の種類を選んでいる。今回はt-test(t検定)を指定した。ただし、ここでt検定を用いることは正しくないかもしれないことに注意すること。t検定は2群ともに分布が正規分布に従っているかを検定した上でとか色々前提条件があったりする。他にはwilcox.testなどを指定することも可能だ。こちらはノンパラメトリック用の検定となっている。

na.rm

na.rm = FALSEはデフォルトでFALSEである。FALSEならば欠落データがあった場合に警告してくれる。Trueにすると警告なしに実行される。

map_signif_level

map_signif_levelはTRUEにすると上記のようにアスタリスクつき(p<0.05で*、p<0.01 **, p<0.001 ***)の表示にしてくれる。FALSEにすると以下のようにp値をバーの上に表示する形式になる。

また、マニュアルでアスタリスクとp値の関係を変えることができる。例えばもっと低いp値の記述をアスタリスク4個でしたい場合は以下のようになる。

plt2 <- plt +   geom_signif(comparisons = list(c("A", "B")),
              test = "t.test",
              na.rm = FALSE,
              map_signif_level = c("****" = 0.0001, "***" = 0.001, "**" = 0.01, "*" = 0.05),
              y_position = 18,
              col = "red") + theme_minimal()

分野によってこのアスタリスクとp値の関係が違ったりするものだが、必ず論文にはアスタリスクの数とp値の関係を示す必要がある。自明だと思わないこと。

y_position

y_positionによって、そのバーの場所を変えることができる。

y軸の値に対応しているので、一旦出力してどこに描画するかを調整しつつやっていくのが良さそうである。

col

col = "red"の部分を修正することでバーやアスタリスク、p値の色を変えられる。Rに登録されている色の単語の他に、カラーコードで指定してやることもできる。

***

今回はさらにtheme_minimal()を追加し、ggplot2によって標準出力されるグラフの見た目を変更した。こっちのほうがスッキリしていて見やすい。これは好みだと思う。

“+”でつないでいき、複数これを繰り返していくことでいくらでも2群間の検定をすることができる。

きれいにかけてうれしいね。

おわり。