ggplot2でRやrstanの分析結果を綺麗に魅せよう

明治大学Advent calender 2018，22日目です．
qiitaは前から利用していましたが投稿するのはこれが初めてです．稚拙な文章で恐縮ですが，頑張って書きますね．
ちなみに大学のアドベントカレンダーということで，学部生でこれから研究やっていくぜ！みたいな方に読んでいただくことを意識して書きたいと思います．

はじめに

皆さん，普段Rは使いますか？
使う方はご存知でしょうが，Rでは非常に簡単にグラフが描画できます．plot(x)の一言命令を実行すれば，折れ線グラフが描けるんですから，本当に便利ですよね．

しかし本日は「いつまでダサいグラフ描くつもり？」というお話．Rのグラフってすごく簡単に描けるんですが，見栄えがイマイチですよね．ダサいグラフではパワポやポスターまでダサくなります．折角の分析結果ですから，格好良く可視化したいですね．目指せインスタ映え！

というわけで本稿ではRで使えるggplot2というパッケージを紹介します．

そしてタイトルにrstanと載せましたが，Stanユーザーのベイズ屋さんに向けた可視化の話を後半に用意しました．応用編という形ですが，ggplot2をご存知の方は後半にジャンプしていただいて結構です．

実行環境

R(Version 3.4.4)の中でパッケージggplot2(Version 2.2.1)を呼び出して実行しました．開発にはRstudio (Version 1.1.456)を使いました．
後半の応用編ではrstan(Version 2.17.4)による計算結果を利用しました．

【基礎編】ggplot2を使おう

例として，irisという有名なデータセットを使います．setosa，versicolor，verginicaという３種類のアヤメについて，がく片長(Sepal Length)，がく片幅(Sepal Width)，花びら長(Petal Length)，花びら幅(Petal Width)を計測したデータセットです．具体的に見たい方はRのコンソールでirisと打って実行してください．

インストール

まずはお使いのPCのRでggplot2を使えるようにします．

install.packages("ggplot2")

以上を実行するだけです．一度インストールすれば，その後は使う前に以下を実行するだけでOKです．

library(ggplot2)

ヒストグラム

まずはggplot2を使わない場合と使う場合で比較してみましょう．以下を実行します．

# Rのデフォルト
hist(subset(iris,Species=="setosa")$Sepal.Length,breaks=seq(4,9,0.3),col="#0f717380",main="Histgram",xlab="Sepal Length")    # setosaのヒストグラムを描画
hist(subset(iris,Species=="versicolor")$Sepal.Length,breaks=seq(4,9,0.3),col="#48beff80",add=TRUE,main=NA,xlab=NA)    # versicolorを描画
hist(subset(iris,Species=="virginica")$Sepal.Length,breaks=seq(4,9,0.3),col="#e0bad780",add=TRUE,main=NA,xlab=NA)    # virginicaを描画
legend("topright",legend=c("setosa","versicolor","virginica"),col=c("#0f7173","#48beff","#e0bad7"),pch=15)    # レジェンド（凡例）を追加

# ggolot2を使った場合
ggplot(iris,aes(x=iris$Sepal.Length,fill=Species)) +    # 使うデータフレームと軸，分け方を宣言
  geom_histogram(binwidth = 0.2,position="identity",alpha=0.8) +    # ヒストグラムを描く
  labs(title="Histgram",x="Sepal Length") +    # タイトル，軸ラベルの指定
  theme(plot.title=element_text(hjust=0.5)) +    # タイトル位置を中央揃え
  scale_fill_manual(values = c("#0f7173","#48beff","#e0bad7"))    # 色の指定

実行結果は以下のようになります．左がRデフォルトで，右がggplot2による描画です．

今回は３種それぞれのSepal Lengthについてヒストグラムを重ねて描画しました．
どうでしょうか？同じ命令しかしていないのに，ggplot2の方がより綺麗な印象を受けませんか．

実は綺麗なのは見た目だけではありません．
Rデフォルトの文法では３行，ggplot2では５行で描きましたが，Rの文法は種ごとに毎回ヒストグラムを描いてオプションadd=TRUEで重ねて描画していくため，アヤメの種類が多ければ多いほど行が増えていきます．
それに対しggplot2の文法では，１行目にfill=Speciesと書けば「Speciesの値ごとに色分けしてね」という意味になるため，アヤメが何種類あってもこの５行よりも増えることはありません．
オプションも，Rデフォルトでは同じ様式なのに毎回指定する必要があり１行が長くなりがちですが，ggplot2なら１回ずつの宣言でよく，それらを足し込んでいく形なのでスッキリしますね．

散布図

次は散布図を描いてみたいと思います．

# Rのデフォルト
plot(subset(iris,Species=="setosa")$Sepal.Length,subset(iris,Species=="setosa")$Petal.Length,xlim=c(4,8),ylim=c(1,7),pch=20,col="#0f7173",main="Scatter Plot",xlab="Sepal Length",ylab="Petal Length")    # setosaの散布図を描画
par(new=T)    # 次の図を重ねて描画する宣言
plot(subset(iris,Species=="versicolor")$Sepal.Length,subset(iris,Species=="versicolor")$Petal.Length,xlim=c(4,8),ylim=c(1,7),pch=20,col="#48beff",main=NA,xlab=NA,ylab=NA)    # versicolorの散布図を描画
par(new=T)    # 次の図を重ねて描画する宣言
plot(subset(iris,Species=="virginica")$Sepal.Length,subset(iris,Species=="virginica")$Petal.Length,xlim=c(4,8),ylim=c(1,7),pch=20,col="#e0bad7",main=NA,xlab=NA,ylab=NA)    # virsinicaの散布図を描画
legend("bottomright",legend=c("setosa","versicolor","virginica"),col=c("#0f7173","#48beff","#e0bad7"),pch=20)    # レジェンド（凡例）を右下に描画

# ggplot2を使う場合
ggplot(iris,aes(x=iris$Sepal.Length,y=iris$Petal.Length,color=Species)) +    # 使うデータフレームと軸，分け方を宣言
  geom_point() +    # 散布図を描画
  labs(title="Scatter Plot",x="Sepal Length",y="Petal Length") +    # タイトルや軸ラベルを指定
  theme(plot.title=element_text(hjust=0.5)) +    # タイトル位置を中央揃え
  scale_color_manual(values = c("#0f7173","#48beff","#e0bad7"))    # 色の指定

Sepal LengthとPetal Lengthの２変量について，種別に色分けした散布図を描画しました．
versicolorとvirginicaはそれぞれ正の相関があるようです．

ggplot2では先ほどのコードを少し書き換えるだけで散布図が描けました．Rデフォルトはやや煩雑になりましたね．plot関数ではadd=TRUEのオプションが使えないらしく，par(new=T)によって重ねて描画する方法をとりました．
どちらのデザインがお好みでしょうか？

デザイン変更

ggplot2のデザイン，いかがでしたでしょうか．グレーの背景や目盛り線が勝手につくところが，個人的には嬉しくお気に入りです．
しかし残念ながらグラフにもTPOというものがありまして，ggplot2のデフォルトデザインが不適切な場もあるのだそうです．

渋々従うとこんな感じ．

ggplot(iris,aes(x=iris$Sepal.Length,y=iris$Petal.Length,color=Species)) +
  geom_point() +
  labs(title="Scatter Plot",x="Sepal Length",y="Petal Length") +
  theme_classic() +
  theme(plot.title=element_text(hjust=0.5)) +
  scale_color_manual(values = c("#0f7173","#48beff","#e0bad7"))

先ほどの散布図のコードに，４行目theme_classic() +を書き足しました．
このようなデザイン変更もggplot2なら簡単です．
とはいえ，個人的にはどんな時でもおしゃれしたいんですけどね…（これならRデフォルトの方が綺麗かも笑）

デザイン変更やTPO考慮時のTipsは以下の記事で紹介してくださっています．
＜ggplot2についてちょっと勉強した(3) -themeを利用した外観の変更-＞
＜学会発表のためのggplot2の設定めも＞

【応用編】rstanの結果を可視化する

ここまでggplot2で基本的な図の描画例をお見せしました．後半は，僕の研究における実際のggplot2利用シーンをご紹介したいと思います．

複雑な統計モデルのパラメータ推定にはMCMCというベイズ推定の手法が有効なのですが，近年統計的プログラミング言語によってこのMCMCの計算が手軽に実装できるようになっています．その代表がStanで，それをRで使えるようにしたパッケージがrstanです．

MCMCは事後分布からの乱数をマルコフ連鎖を用いて擬似的に大量生成する方法ですが，計算後得られたサンプルについて収束や自己相関をチェックする必要があります．その最も簡単なチェック方法は可視化だと言えます．

rstanデフォルトの可視化関数

前述したチェック用の関数は，rstanに標準搭載されています．
今回は題材としてStanとRでベイズ統計モデリングの12章から状態空間モデルのパラメータ推定の結果を用います．データ，モデル，実行用スクリプトはサポートサイトに掲載されています．

trace plot

パッケージrstanデフォルトの収束チェック関数を実行すると，

stan_trace(fit,pars=c('mu_all[1]'))

このようになります．

うーん，ちょっとダサい．急ぎ収束チェックしたい時などは便利そうですが，色が変えられないので，スライドやポスターに載せる用には向きませんね．
どうしてよりによって紫とオレンジなんでしょうか？？ハロウィン？？

自己相関関数

rstanデフォルトの関数で自己相関関数を可視化すると，

stan_ac(fit,pars=c("mu_all[1]"))

赤か〜…先ほどのtrace plotほどではないですが，個人的には気に入らないデザイン…

このように，rstanには事後チェック用の関数がついていて大変便利なのですが，デザインがいまひとつ．これじゃスライドやポスターには載せられません．
ちなみにこのようなラグと自己相関のプロットのことをコレログラムと呼びます．

そこでggplot2の出番です．

ggplot2でサンプル列のチェック

実行用スクリプトを見ていただくと，rstanの結果はfitという変数に格納されることがわかります．
可視化するためにfitの中身を自由に扱いたいのですが，この変数はS4クラスという珍しい変数型で，扱いが少し難しいです．

変数fitからサンプル抽出

まずfitの正体を暴くためにstr(fit)を実行すると，死ぬほど長文が返されます．（長すぎるので載せませんが）こんなに色々抱えていたとは…

そこで，使う部分だけ抽出したいと思います．

ext_fit <- extract(fit)

すると今度は

> str(ext_fit)
List of 4
 $ mu_all: num [1:1600, 1:24] 11.2 11.3 10.9 11.2 11.3 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ iterations: NULL
  .. ..$           : NULL
 $ s_mu  : num [1:1600(1d)] 0.412 0.538 0.332 0.43 0.403 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ iterations: NULL
 $ s_Y   : num [1:1600(1d)] 0.22 0.0422 0.1824 0.1002 0.0574 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ iterations: NULL
 $ lp__  : num [1:1600(1d)] 30.4 59.3 38.2 46.6 58.5 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ iterations: NULL

かなり短くなりましたね．
今回は4チェーンをそれぞれ4000回繰り返しサンプルして，初めの2000回分は初期値依存期間として捨てています．また，サンプル系列間の自己相関をなくすために5倍に薄めて使っています．すなわち得られたサンプル数は$4 \times (4000-2000)/5 =1600$個．
ですから，ext_fit$mu_allは得られた1600個のサンプルを24個の変数ごとに格納した行列形式のデータということになりますね．このように，extract関数によっていとも簡単にサンプル抽出ができました．

trace plot

では，trace plotを書きましょう．

library(ggplot2)    # ggplot2を読み込み
trace <- function(x,Nchains)    # trace plot描画用の関数
{
  xdf <- data.frame(Index=rep(rep(1:(length(x)/Nchains)),Nchains),
                    chain=rep(1:Nchains,each=length(x)/Nchains),
                    data=x)    # ベクトルで入力されるサンプル列をデータフレーム化
  ggplot(xdf,aes(x=Index,y=data,colour=factor(chain))) +
    geom_line() +
    labs(x="iteration",y="state") +
    scale_colour_manual(values = c("#0f7173","#48beff","#e0bad7","#70f8ba"),name="chain")    # 好きな色を指定
}
trace(ext_fit$mu_all[,1],4)    # 描画

実行するとこんな感じ．

綺麗に描けました．
定常分布に収束していそうですね．

自己相関関数

自己相関関数も描画しましょう．

library(ggplot2)    # ggplot2を読み込み
aveac <- function(x,Nchains,maxlag=20)    # 自己相関関数のチェーン平均を描画する関数
{
  len <- length(ext_fit$mu_all[,1])/Nchains    # 各チェーンのサンプル長
  y <- array(NA,dim=c(Nchains,(maxlag+1)))
  for(i in 1:Nchains)    # チェーンごとに自己相関を求める（acf関数を利用）
  {
    y[i,] <- acf(ext_fit$mu_all[((i-1)*len+1):(i*len),1],lag.max=maxlag,plot=FALSE)$acf[,1,1]
  }
  ave_y <- apply(y,2,mean)    # チェーン平均を求める
  ggplot(data.frame(ac=ave_y,lag=c(0:maxlag)),aes(x=lag,y=ac,fill=factor(0))) +
    geom_bar(stat="identity") +
    labs(title="average autocorrelation",y="average autocorrelation") +
    scale_fill_manual(values = c("#70f8ba"), guide = "none")    # 好きな色を指定
}
aveac(ext_fit$mu_all[,1],4)

実行するとこんな感じ．

きっちり自己相関が消えていますね．自己相関が強いと独立なサンプルとみなせないのでもっと強い倍率で薄めなくてはならないのですが，これだけ弱ければ十分でしょう．
ちなみにStanはHMC法というMCMCアルゴリズムを積んでいて，他のMCMCサンプラーのWinBUGSやJAGSよりも自己相関が弱くなりやすいです．素晴らしい！

おわりに

前半と後半で内容の温度差が大きくて申し訳ありませんでした．
Rやrstanのデフォルト関数もお手軽で悪くはありませんが，今回比較したようにggplot2はスマートなコーディングで綺麗なグラフを描ける便利パッケージです．
これから研究やってくぜ！という学部生の方には特に触れていただきたいツールです．

後半ではStanの出力結果からサンプル列を抽出する方法，そしてそれらのサンプル列のチェックをggplot2で可視化する方法をご紹介しました．既存の関数も便利ですが，グラフを書くならデザインをカスタマイズできた方が良いですよね．
ただ実は，ggplot2でMCMCのグラフを作成するためのパッケージ，ggmcmcというツワモノがおりまして．自分はggmcmcは利用していないため詳しくありませんが，今回紹介した関数よりも実はそちらの方が便利かもしれません．

参考文献

ggplot2に関すること
＜ggplot2を使ったデータのビジュアライゼーション＞

stanに関すること
＜StanとRでベイズ統計モデリング　１２章＞