StanでPearsonの相関行列をベイズ推定

author: Unadon (見習い飯炊き兵) 動作環境:Mac OS Sierra 10.12.1; R version3.3.1; rstan 2.10.1

BayesCorr.png

 

はじめに

ベイズで相関係数を推定。

相関関係を考えたい変数がたくさんある場合を想定したコードをご紹介。毎度毎度お世話になっている こちらのスライドを参考にさせていただきました。 推定された相関係数の分布を可視化するところまでやっていきます。

本稿の目的

  1. 多変量の相関係数を推定する(ピアソンの積率相関係数)
  2. 相関係数の分布をマトリックス形式で可視化する
  3. MAP推定値を得る

パッケージの読み込み

irisのデータを使用します。コードには汎用性をつけていますので、ご自身のデータを読み込ませても構いません。

stanやggplot2以外にも、あとで使用するものがあるので、読み込んでおきます。

##################################################
#Stan : Correlations
##################################################
invisible({rm(list=ls());gc();gc()})

#library
library(rstan) #HMC sampling
library(ggplot2) #graphic
library(shinystan) #Stan results
library(reshape2) #data reshape
library(tidyverse)  #tidy data handling
library(plotly) #interactive graph
library(Rmisc) #summary correlation
library(GGally) #panel plot
library(ellipse) #なんやっけ
library(formattable)#HTML tableの出力

#Stan並列化
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

初期設定とデータの読み込み

stanのサンプリングに関する設定をし、パラメータを宣言しておきます。 パラメータ名はいじる必要がありません。

データと使用したい変数名の部分を皆さんのデータにあわせて変更してください。

#Initaial Settings-------------------------------------------
#サンプリング数
iterations =c(2000)
#Warm-up数
warmup=c(1000)
#チェーン数
chains=c(3)
#thinning
thin =c(1)
#推定に使用するパラメータ
parameters=c("rho","sigma","Sig") #log_likは必ず最後
#あとで結果を確認したいパラメータ
parameters2<-c("rho")

#データセットの読み込み(行は参加者,列は変数。Headerあり)
#dat<-read.csv("          .csv")
#read data サンプルデータ
dat<-iris

#変数間相関を求めたいデータの列名
varNames<-c("Sepal.Length","Sepal.Width","Petal.Length")

データ成形

上のブロックでデータと列名を指定したら、あとは自動的にデータ変換をしてStanに渡すデータセットにします。

Datが抽出した列のデータをmatrix形式にしたもの。Nvarで変数の数、Nsubで行数(参加者数)を定義しています。

spdという変数は、偏差平方和・積和行列と呼ばれる行列データです。

分散・共分散行列というワードを聞いたことが有る方もいらっしゃると思いますが、

cov()でそれを作成し、本来(N-1)で除算されていたものに(N-1)を乗算することで、

偏差平方和・積和行列を作っています。


#データ列の抽出
Dat<-data.matrix(dat[,varNames])

#Number of variables and Subjects
Nvar<-length(Dat[1,])
Nsub<-length(Dat[,1])

#sqrt(Cov)*SD matrix
spd<-cov(Dat)*(Nsub-1)

#data and parameters to be passed on stan
datastan = list(Nsub=Nsub,Nvar=Nvar,spd=spd)

stanの実行

stan model

こちらのスライドのStanレベル9あたりに詳細があります。 偏差平方和積和行列をデータとして読み込み、標準偏差sigmaと相関行列rhoをパラメータとして推定します。 相関係数ではなく、相関行列です。すごいですね。

transformed parameters{}のブロックで、推定された相関行列と標準偏差SDを使って共分散行列を作っています。 そして、その共分散行列とN-1をwishart分布の引数に入れて、偏差平方和積和行列spdをデータとしています。

wishart分布ってなんなのでしょうか。lkj_corrという分布は、もはや読み方からして謎です。

勉強しておきます。

modelCorr<-"
data {
    int Nsub;
    int Nvar;
    matrix[Nvar,Nvar] spd;//偏差平方和・積和行列
}

parameters{
    vector<lower=0>[Nvar]sigma;//標準偏差のベクトル
    corr_matrix[Nvar]rho;//相関行列
}

transformed parameters{
    cov_matrix[Nvar] Sig;//共分散行列
    Sig = quad_form_diag(rho,sigma);//相関行列とSDを共分散行列にする
}

model{
    rho ~ lkj_corr(1);//uninform prior
    spd ~ wishart(Nsub-1,Sig);
}

"

さて、上のモデルをつかってサンプリングを実行していきます。

#StanFit
fitCorr<-stan(model_code = modelCorr,data=datastan,pars=parameters,
          iter = iterations,chains = chains,thin = thin)

#Summary
Summary<-summary(fitCorr)$summary[c(1:c(Nvar*Nvar)),c(1,3,4,8,9,10)]

#結果のサマリーを確認
formattable(Summary)
Table. 出力結果の例 こちらの表の出力方法は、拙記事(5項目)を参照してください。

相関行列の可視化

今回の結果をggplot2で可視化していきます。

最近tidyverseを使うようになったので、以下のコードを改めてみると汚いなぁと思いますが、

変数列が増えても問題なく動きますので、そのまま使っていただいて構いません。

#Extract
d.ext<-rstan::extract(fitCorr)    
df.ext<-data.frame(d.ext[c(parameters2)])

#melting to long data
df.long<-melt(df.ext,idvar=c())

#col name separation
df.comp<-df.long %>% separate(variable, c('parameter', 'varNamesRow',"varNames"))

#replace Num of parameter to varNames
for(i in 1: Nvar){
    df.comp$varNamesRow<-replace(df.comp$varNamesRow,df.comp$varNamesRow==as.character(i), varNames[i]) 
    df.comp$varNames<-replace(df.comp$varNames,df.comp$varNames==as.character(i), varNames[i]) 
    }
df.comp <- df.comp[order(df.comp$varNames),]
df.comp$varNames<-as.factor(df.comp$varNames)
df.comp$varNamesRow<-as.factor(df.comp$varNamesRow)

#theme set
ggplot()+theme_set(theme_bw(base_size=18,base_family="HiraKakuProN-W3"))

#plot
g <- ggplot(df.comp,aes(x = value,y = ..density..,colour=varNames,fill=varNames))+
     geom_histogram(position = "identity",alpha = 0.5)+
     geom_density(stat = "density",position = "identity",alpha = 0.4)+
     facet_wrap(varNamesRow~varNames,scales="free")+
     scale_color_brewer(palette="Set2")+
     scale_fill_brewer(palette="Set2")+
     ggtitle("Posterior Distributions")+
     theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) +
     theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))+
     theme(legend.title = element_text(size=14), #legend size
               legend.text = element_text(size=10),  #legend text size 
               legend.key= element_rect(colour = "white")) # legend mark-box colour

plot(g)

よいしょ。

BayesCorr.png

少しですが左右対称の事後分布ではないですね。 これらの相関係数は、推定値が1または−1に近い場合、より非対称な分布になるため、最も確率密度の高いMAP推定値を代表値として採用する場合があります。

MAP推定値の算出

各相関係数について、MAP推定値を求めます。行列にして出力しましょう。出力はformattable()を使います。

#MAP estimate-------------------------------------------------------
MAPs  = c(rep(NA,c(Nvar*Nvar)))
for (i in 1:c(Nvar*Nvar)){
    MAPs[i]    <- density(df.ext[,i])$x[which(density(df.ext[,i])$y == max(density(df.ext[,i])$y))]
}
#置き換えるべき項目数についてのerrorが出ても構わない
#同じ変数の相関なので、MAP推定値”1”が2個以上選出されるため

names(MAPs) <- c(colnames(df.ext))
MAPcorMat<-matrix(MAPs,ncol = Nvar,nrow = Nvar)
rownames(MAPcorMat)<-varNames
colnames(MAPcorMat)<-varNames
formattable(data.frame(MAPcorMat))

相関行列(MAP推定値) CorrMat.png

ベイズ推定結果と通常の積率相関係数の対応性

Corr.png

念のため、通常の相関分析の結果と対応しているか見ておきました。大丈夫そうですね。


#RawData and PeasonCorrs(not Bayesian)--------------------------------------------------------------
d <- as.data.frame(Dat)
d.order<-order(colnames(Dat))
d <- d[ ,d.order]

N_col <- ncol(d)
ggp <- ggpairs(d, upper='blank', diag='blank', lower='blank',title = "RawData/Pearson-Corr(NotBayesian)")

for(i in 1:N_col) {
    x <- d[,i]
    p <- ggplot(data.frame(x), 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))+scale_colour_brewer(palette="Set2")
    } else {
        bw <- (max(x)-min(x))/10
        p <- p + geom_histogram(binwidth=bw, 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='pearson', 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))+scale_fill_brewer(palette="Set2")
        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), aes(x=x, y=y))
        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)
    }
}

ggp

ではでは、

Enjoy Stan!!

Written on February 17, 2017