【R言語】ggplot2｜棒グラフに有意差を示すスターをつけてみる

こんにちは。プログラミング超初心者のえいこです。

前回、スターをつけるのに便利な関数パッケージであるgeom_signifをエラーメッセージと闘いながらインストールしました。

グラフを描こうと思ったらggpot2やdplyrが使えなくなったりと...紆余曲折があったのです。

が、今回はいよいよ作った棒グラフにスターをつけていこうと思います。

geom_signifなどの使い方を参考にしたのはこちらのサイト↓

ggplot2でお絵かき:棒グラフ - James' blog

振り返り：棒グラフの作り方

棒グラフを描く手順はこんな感じでした。

1. ggplotとdplyrを呼び出しておく
2. 作成したcsvファイルを読み込む
3. 読み込んだデータをもとに平均値・標準偏差などのリストを作っておく
4. geom_barを使って棒グラフを描く

今回は、スターをつけたりするので3.と4.の間に↓の3つの手順を追加しておきます。

a.  エラーバーをつけるための変数を事前に準備しておく
b.  二標本検定（t検定）を事前にしてpの値を文字列にしておく
c. yの最大値を定義しておく

a. エラーバーをつけるための変数を事前に準備しておく

私はerrorsという変数を準備して平均値±標準偏差のデータセットを入れておきます。

errors <- aes(ymax = mean + sd, ymin = mean - sd)

これを先に準備しておくことで、サクッとエラーバーが描けることがわかりました。

b. 二標本検定（t検定）を事前にしてpの値を文字列にしておく

pの値を文字列に変換する手順は、Rで棒グラフを描いていた時にまとめています。

今回も、スターを書いた隣にpの値を書いておこうと思っているのでpの値を事前に文字列にしました。

別にスターの数だけで良いならこんなことはしなくても良いかもしれません。

今回書いたコードはこちら↓

pvalue<-t.test(rawdata$value ~ rawdata$treatment)$p.value
c.pvalue<-as.character(round(pvalue,digits=4))
p.pvalue<-paste("p","=",c.pvalue,sep="")

c. yの最大値を定義しておく

これを棒グラフを描く前にしておくことで、データセットによって数値を変える場所を少なくできます。

#y軸の最大値を定義しておく
ylim=(max(rawdata$value)+ sd(rawdata$value))*1.2

では、棒グラフを描いていきます。

前回は一つ一つのステップごとに出力して、グラフの形を確認していました。

が、今回は前回書いたコードをすべて"+"でつなげて行っています。書いたコードはこちら↓

#棒グラフを作っていく
barplot <-
 ggplot(data_mean_sd, #使うデータセット
        aes(x = treatment,　#x軸
            y = mean,　#Y軸
            fill = treatment #塗りつぶしを処理によって変える
            )
        ) + 
 
 geom_bar(stat = "identity",　#統計処理しない数値を使用
          position = "dodge", #横に並べて表示
          colour = "black" #線の色は黒
          ) + 
 
 geom_errorbar(errors, #使うデータは"errors"
               width = 0.1 #幅は0.1
               ) +
 
 scale_fill_manual(values = c("gray", "white") #ぬりつぶしは灰色と白
                   ) +
 
 labs(x = "", y = "Expression level") +　#X軸とY軸のラベルを設定する
 
 theme_classic()+ #テーマを設定する
 
 scale_y_continuous(expand = c(0,0), limits = c(0,ylim))

前回と違うところは、"geom_errorbar"のところ。前回はこんなコードを書いていたのですが、

geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),width=.2)

今回は事前に変数で定義しているので、使うデータのところを"errors"のみで表記できるようになっています。

変数ってすごい！！！

やっと今回の本題geom_signifを使ってスターをつけていきたいと思います。

初めてgeom_signifを使ってみる

geom_signifってどんな関数なの？使うときにどんな引数が必要なの？

というのが最初に浮かぶ疑問ですよね。では、helpを見てみましょう。

Description
Create significance layer
Usage
stat_signif(mapping = NULL, data = NULL, position = "identity",
na.rm = FALSE, show.legend = NA, inherit.aes = TRUE,
comparisons = NULL, test = "wilcox.test", test.args = NULL,
annotations = NULL, map_signif_level = FALSE, y_position = NULL,
xmin = NULL, xmax = NULL, margin_top = 0.05, step_increase = 0,
tip_length = 0.03, size = 0.5, textsize = 3.88, family = "",
vjust = 0, parse = FALSE, manual = FALSE, ...)

geom_signif(mapping = NULL, data = NULL, stat = "signif",
position = "identity", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE, comparisons = NULL, test = "wilcox.test",
test.args = NULL, annotations = NULL, map_signif_level = FALSE,
y_position = NULL, xmin = NULL, xmax = NULL, margin_top = 0.05,
step_increase = 0, tip_length = 0.03, size = 0.5,
textsize = 3.88, family = "", vjust = 0, parse = FALSE,
manual = FALSE, ...)

んー...mappingってなんだかわからないし、何から入れた良いかさっぱりわかりません...

ということで、実際に使っている人のコードを参考にしてみました。

Xの座標、Yの座標、入れる文字列、線の装飾

だけでも良いみたい！

ということで、こんな感じにコードを書いてみました。

geom_signif(
   y_position=ylim*0.8,　#高さの設定
   xmin=1.0, #線の始点の設定
   xmax=2.0, #終点の設定
   annotation= #if文を使ってスターの数を変える
     if(pvalue<0.001){
       "***"
     }else if(pvalue<0.01){
       "**"
     }else if(pvalue<0.05){
       "*"
     }else{
       "n.s."
     },
   tip_length=0.05, #棒の下向きの線の長さ
   textsize=7　#スターの大きさ
   )

スターの数は、if文を使って表示させるというのがプログラミングをする醍醐味だと思っているので、以前の記事で作ったif文をそのままannotationのところに入れてみました。

スターの大きさは、”textsize”で変えられます。描いたグラフを見ながら微調整すると良いかもしれません。
（私の好みは7の時だったので、今回は7にしています）

このコードを、前の棒グラフのコードに"+"を使ってつなげます。

出力されたグラフは...

ちゃんと、スターが2つで表示されてる！！！

あとはp値を書き込むだけ♪

annotateを使ってp値を書き込む

グラフに何かを書き込みたいときには、"annotate"を使うようです。

今回は文字列を入れるので、annotate("text"~)で書き込む位置と文字列を指定していきます。書いたコードはこちら

 annotate("text", #グラフの中にテキストを入れるよ
          x=1.95, #xの場所
          y=ylim*0.85, #高さを調節
          label=p.pvalue, #書き入れる文字列
          fontface="italic", #書体はイタリックに変更
          size=3) #文字サイズ

fontfaceでボールドにしたり、イタリックにしたりできるようです。今回は、ちょっとしたこだわりでイタリックにしています。

このコードを例のごとく"+"で前のコードとつないでいきます。すると...

こんな感じのグラフになりました。

おまけ：作ったグラフを保存する

ggplot2では"ggsave"で指定すれば、作ったグラフを保存できます。

ファイルの名前、保存するグラフの名前、画質、大きさなどを指定できます。

今回はこんな感じのコードを最後に付け加えました。

ggsave(file="bar-plot.png",plot=barplot, dpi=100, width=4.0, height=5.7)

保存されるディレクトリは、カレントディレクトリ。保存場所に移動するコードを直前に入れておいた方が良いかもしれません。

保存できたグラフはこちら。

Rで１か月かけて地道にやってきたことが１週間ちょっとでできるなんてggplot2恐るべし！！

でも、１か月かけて勉強してきたことは全く無駄じゃなく、if文やpの値を文字列にするなどggplot2でも使っている知識はたくさんあります。

何事も積み重ねが大事ですね。

次回は、いよいよドットプロットと棒グラフを重ね合わせてみたいと思います。

それでは、また！

今回書いたコードの全体像はこちら↓

library(ggplot2) #ggplot2を呼び出す
library(ggsignif) #ggsignifを呼び出す
library(dplyr)　#dplyrを呼び出す

setwd("C:/Users/～～")#ファイルを呼び出すディレクトリを指定する

rawdata<- read.csv("xtreatment.csv",header=T,sep = ",") 

#平均値と標準偏差のリストを作る
data_mean_sd <- 
 rawdata %>%
 group_by(treatment) %>% 
 summarize(mean = mean(value), sd = sd(value))　
data_mean_sd

#エラーバーをつけるための変数を準備する
errors <- aes(ymax = mean + sd, ymin = mean - sd) 

#ttestするのpの値を文字列にする
pvalue<-t.test(rawdata$value ~ rawdata$treatment)$p.value
pvalue
c.pvalue<-as.character(round(pvalue,digits=4))
c.pvalue
p.pvalue<-paste("p","=",c.pvalue,sep="")
p.pvalue

#y軸の最大値を定義しておく
ylim=(max(rawdata$value)+ sd(rawdata$value))*1.2
ylim

#棒グラフを作っていく
barplot <-
 ggplot(data_mean_sd, #使うデータセット
        aes(x = treatment,　#x軸
            y = mean,　#Y軸
            fill = treatment #塗りつぶしを処理によって変える
            )
        ) + 
 
 geom_bar(stat = "identity",　#統計処理しない数値を使用
          position = "dodge", #横に並べて表示
          colour = "black" #線の色は黒
          ) + 
 
 geom_errorbar(errors, #使うデータは"errors"
               width = 0.1 #幅は0.1
               ) +
 
 scale_fill_manual(values = c("gray", "white") #ぬりつぶしは灰色と白
                   ) +
 
 labs(x = "", y = "Expression level") +　#X軸とY軸のラベルを設定する
 
 theme_classic()+ #テーマを設定する
 
 scale_y_continuous(expand = c(0,0), limits = c(0,ylim))+　#Y軸を設定する
 
 #有意差を表すスターを描いていく
 geom_signif(
   y_position=ylim*0.8,　#高さの設定
   xmin=1.0, #線の始点の設定
   xmax=2.0, #終点の設定
   annotation= #if文を使ってスターの数を変える
     if(pvalue<0.001){
       "***"
     }else if(pvalue<0.01){
       "**"
     }else if(pvalue<0.05){
       "*"
     }else{
       "n.s."
     },
   tip_length=0.05, #棒の下向きの線の長さ
   textsize=7　#スターの大きさ
   )+
 annotate("text", #グラフの中にテキストを入れるよ
          x=1.95, #xの場所
          y=ylim*0.85, #高さを調節
          label=p.pvalue, #書き入れる文字列
          fontface="italic", #書体はイタリックに変更
          size=3) #文字サイズ

barplot

ggsave(file="bar-plot.png",plot=barplot, dpi=100, width=4.0, height=5.7)