StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

散布図行列を描くには (corrplot, pairs, GGally)

データが与えられた時にはまず可視化をします。そのデータがどのような仕組み(メカニズム)で作られてそうなったかを考えるために必須のプロセスです。しかしながら、どんな可視化がベストかははじめの段階では分からず、とにかくプロットしまくることになります。そのとっかかりに僕がよく使うのが散布図行列(scatter matrix,scatter plot matrix)です。

今回は3つほど紹介します。

{corrplot}パッケージの corrplot()関数

library(corrplot)

d <- iris
d <- d[,-ncol(d)]
M <- cor(d, method='spearman', use='pairwise.complete.obs')

png(file='corrplot.png', res=250, w=1500, h=1500)
corrplot.mixed(M, order = 'hclust')
dev.off()

5行目で相関係数行列を作ってそれを渡しておしまいです。相関係数行列の作り方は各自の自由です。上記ではSpearmanの順位相関係数を使っていますがMICとかでもいいと思います。

このcorrplotのデメリットとしましては散布図は表示できない点です。散布図行列と言っておきながらすみません。説明変数が100個以上あるときなどは散布図を描くのは苦しいのでこちらを使う場合が多いです。

{graphics}パッケージ(デフォルトで入っています)の pairs()関数

d <- iris
d <- d[,-ncol(d)]
hc <- hclust(as.dist(1-cor(d, method='spearman', use='pairwise.complete.obs')))
hc.order <- order.dendrogram(as.dendrogram(hc))
d <- d[ ,hc.order]
gr <- as.factor(iris$Species)

cols.key <- scales::muted(c('red', 'blue', 'green'))
cols.key <- adjustcolor(cols.key, alpha.f=1/4)
pchs.key <- rep(19, 3)

panel.hist <- function(x, ...) {
   usr <- par('usr'); on.exit(par(usr))
   par(usr=c(usr[1:2], 0, 1.5))
   h <- hist(x, plot=FALSE)
   breaks <- h$breaks
   nB <- length(breaks)
   y <- h$counts; y <- y/max(y)
   rect(breaks[-nB], 0, breaks[-1], y, col='gray', ...)
}

panel.cor <- function(x, y, ...){
   usr <- par('usr'); on.exit(par(usr))
   par(usr=c(0,1,0,1))
   r <- cor(x, y, method='spearman', use='pairwise.complete.obs')
   zcol <- lattice::level.colors(r, at=seq(-1, 1, length=81), col.regions=colorRampPalette(c(scales::muted('red'),'white',scales::muted('blue')), space='rgb')(81))
   ell <- ellipse::ellipse(r, level=0.95, type='l', npoints=50, scale=c(.2, .2), centre=c(.5, .5))
   polygon(ell, col=zcol, border=zcol, ...)
   text(x=.5, y=.5, lab=100*round(r, 2), cex=2, col='white')
   # pval <- cor.test(x, y, method='spearman', use='pairwise.complete.obs')$p.value
   # sig <- symnum(pval, corr=FALSE, na=FALSE, cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), symbols = c('***', '**', '*', '.', ' '))
   # text(.6, .8, sig, cex=2, col='gray20')
}

panel.scatter <- function(x, y){
   points(x, y, col=cols.key[gr], pch=pchs.key[gr], cex=1)
   lines(lowess(x, y))
}

png('pairs.png', res=250, width=1500, height=1500)
pairs(d,
   diag.panel=panel.hist,
   lower.panel=panel.scatter,
   upper.panel=panel.cor,
   gap=0.5,
   labels=gsub('\\.', '\n', colnames(d)),
   label.pos=0.7,
   cex.labels=1.4
)
dev.off()
  • 3~5行目: 変数の並びがhclustした結果の並びになるようにしています。panel.cor()のところは相関係数の下二桁を表示して、その値に応じて色と楕円の傾きを変えています。

オプションや使い方については各自ヘルプやstackoverflow見るなり、実際にいじって出力してみるなりして調べてください。コメントアウト部分につきましては「p-valueに応じたお星さまを頼むからつけてくれ。でないとお前を殺して俺も死ぬ。」と言われた時だけアンコメントして下さい。

このpairs()は拡張性が高く、デフォルトのplotに詳しい人にはよいと思います。

{GGally}パッケージの ggpairs()関数

こちらは{ggplot2}に慣れている人向けです。

library(ggplot2)
library(GGally)

d <- iris
d <- d[,-ncol(d)]

N_col <- ncol(d)
ggp <- ggpairs(d, upper='blank', diag='blank', lower='blank')

for(i in 1:N_col) {
  x <- d[,i]
  p <- ggplot(data.frame(x, gr=iris$Species), aes(x))
  p <- p + theme(text=element_text(size=14), axis.text.x=element_text(angle=40, vjust=1, hjust=1))
  if (class(x) == 'factor') {
    p <- p + geom_bar(aes(fill=gr), color='grey20')
  } else {
    bw <- (max(x)-min(x))/10
    p <- p + geom_histogram(binwidth=bw, aes(fill=gr), color='grey20')
    p <- p + geom_line(eval(bquote(aes(y=..count..*.(bw)))), stat='density')
  }
  p <- p + geom_label(data=data.frame(x=-Inf, y=Inf, label=colnames(d)[i]), aes(x=x, y=y, label=label), hjust=0, vjust=1)
  ggp <- putPlot(ggp, p, i, i)
}

zcolat <- seq(-1, 1, length=81)
zcolre <- c(zcolat[1:40]+1, rev(zcolat[41:81]))

for(i in 1:(N_col-1)) {
  for(j in (i+1):N_col) {
    x <- as.numeric(d[,i])
    y <- as.numeric(d[,j])
    r <- cor(x, y, method='spearman', use='pairwise.complete.obs')
    zcol <- lattice::level.colors(r, at=zcolat,
      col.regions=colorRampPalette(c(scales::muted('red'), 'white', scales::muted('blue')), space='rgb')(81))
    textcol <- ifelse(abs(r) < 0.4, 'grey20', 'white')
    ell <- ellipse::ellipse(r, level=0.95, type='l', npoints=50, scale=c(.2, .2), centre=c(.5, .5))
    p <- ggplot(data.frame(ell), aes(x=x, y=y))
    p <- p + theme_bw() + theme(
      plot.background=element_blank(),
      panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
      panel.border=element_blank(), axis.ticks=element_blank()
    )
    p <- p + geom_polygon(fill=zcol, color=zcol)
    p <- p + geom_text(data=NULL, x=.5, y=.5, label=100*round(r, 2), size=6, col=textcol)
    ggp <- putPlot(ggp, p, i, j)
  }
}

for(j in 1:(N_col-1)) {
  for(i in (j+1):N_col) {
    x <- d[,j]
    y <- d[,i]
    p <- ggplot(data.frame(x, y, gr=iris$Species), aes(x=x, y=y, color=gr))
    p <- p + theme(text=element_text(size=14), axis.text.x=element_text(angle=40, vjust=1, hjust=1))
    if (class(x) == 'factor') {
      p <- p + geom_boxplot(aes(group=x), alpha=3/6, outlier.size=0, fill='white')
      p <- p + geom_point(position=position_jitter(w=0.4, h=0), size=1)
    } else {
      p <- p + geom_point(size=1)
    }
    ggp <- putPlot(ggp, p, i, j)
  }
}

png(file='ggally.png', res=250, w=1500, h=1500)
print(ggp, left=0.3, bottom=0.3)
dev.off()
  • 10~23行目: 対角線の部分を付け加えています。21行目で枠内にcolnameを表示するように工夫しています。
  • 28~47行目: 上三角の部分を付け加えています。さきほどのpairs()の時とほぼ同様です。
  • 65~67行目: ggpairs()で作られたオブジェクトは{ggplot2}パッケージのggsave()ができません。print()で出力します。

ディスプレイで見えにくい時は最後の手段として、A1サイズなどの大きな紙に印刷して閲覧するというのをオススメしています。