stanで線形モデル 視覚的事後予測チェック

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

 

今回は、これができるようになります: 線形モデルの切片と傾きの推定と事後予測チェック

はじめに

rstanで色々ベイズしてみたい!という方に向けて、基本的な解析とその周辺Tipsをご紹介しています。

前記事 をご参照ください。

今回は「線形モデルと事後予測チェック」です。

目的

  1. Rとrstanで線形モデルによる予測ができるようになる

  2. 事後予測チェックができるようになる

  3. 結果の図を作成する。

ライブラリの読み込み

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

#library and directory
library(rstan)
library(ggplot2)
library(shinystan)
library(reshape2)
library(loo)
library(tidyverse) 
library(plotly)
library(Rmisc)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#
# 目的: 線形モデルで傾きを推定する
#

#Options-----------------------------------------------------
# 1. WAIC&LOO: 情報量規準WAIC&LOOの算出
# 2. Parameters Visualization of beta
# 3. MAP estimate 
# 4. Posterior Predictive Check
#------------------------------------------------------------


推定の設定とデータの読み込み


#Initaial Settings-------------------------------------------
#サンプリング数
iterations =c(2000)
#Warm-up数
warmup=c(1000)
#チェーン数
chains=c(3)
#thinning
thin =c(1)

#推定に使用するパラメータ
parameters=c("beta","sigma","mu","log_lik") #取り出したいパラメータを最初に並べる
#あとで結果を確認したいパラメータ
parameters2<-c("beta")#推定に使用するパラメータの先頭から順になるように

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


#線形モデルの独立変数
varNames<-c("Sepal.Length","Sepal.Width","Petal.Length")

#線形モデルの従属変数
varY<-c("Petal.Width")


#Data shaping
Dat1<-data.matrix(dat[,varNames])#independent var
Y<-dat[,varY] #dependent var

#------------------------------------------------------------

データの加工

線形代数を使って高速演算を実現したいので、工夫を入れています。データの一列目に切片推定のためのデータを入れていきます。 何を入れるかというと、全部1の値。

データの加工とモデルの書き方についての詳細は、以下のスライドを参照させていただきました。

http://www.slideshare.net/simizu706/stan-62042940

#一列目の変数名をInterceptに
varNames<-c("Intercept",varNames)

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

#線形代数で解くので,切片用に1の列を作成して1列目に結合
InterceptDat<-c(rep(1,Nsub))
Dat<-cbind(InterceptDat,Dat1)

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

今回のデータは独立変数3列(+切片の疑似データ1列)が入ったDatと、従属変数Yが入ったYの2つがメインのデータです。

head(Y) #Petal.Widthのデータ
 [1] 0.2 0.2 0.2 0.2 0.2 0.4

head(Dat) #切片の擬似データと3つの独立変数 

     InterceptDat Sepal.Length Sepal.Width Petal.Length
[1,]            1          5.1         3.5          1.4
[2,]            1          4.9         3.0          1.4
[3,]            1          4.7         3.2          1.3
[4,]            1          4.6         3.1          1.5
[5,]            1          5.0         3.6          1.4
[6,]            1          5.4         3.9          1.7

Stanのモデルコード

以下はStanのモデルコードです。

LinearModel<-"
data {
    int <lower=0>Nvar;
    int <lower=0>Nsub;
    int <lower=0>NY;
    vector [Nsub]Y;
    matrix [Nsub,Nvar+NY] Dat;
}

parameters {
    vector [Nvar+NY]beta;
    real <lower=0>sigma;
}

transformed parameters {
    vector [Nsub]mu;
    mu = Dat*beta;
}

model {
    Y~ normal(mu,sigma);
}

generated quantities {
    matrix [Nsub,NY] log_lik;
    for (i in 1 : Nsub){
        for (j in 1 : NY){
            log_lik[i,j] = normal_lpdf(Y[i]|mu[i],sigma);
        }
    }
}

"

サンプリングの実行と推定結果の要約

サンプリングを実行します。結果の出力は、Summaryのところのコマンドです

fit<-stan(model_code=LinearModel,data=datastan,pars=parameters,
          iter = iterations,chains = chains,thin = thin)

#Summary
Summary<-summary(fit)$summary[c(1:c((Nvar+NY)*c(length(parameters2)))),c(1,3,4,8,9,10)]

beta[1]が切片。あとは独立変数に対応した傾きが順番に入っています。 beta[3]と[4]の95%信用区間がゼロをを超えており、[2]が下回っています。

 print(Summary)
 
              mean         sd       2.5%      97.5%     n_eff     Rhat
beta[1] -0.2396775 0.18797519 -0.6099599  0.1154222 1163.2953 1.001589
beta[2] -0.2062941 0.05119662 -0.3039755 -0.1051271  800.0843 1.002535
beta[3]  0.2216302 0.04840735  0.1302184  0.3134573 1016.1903 1.002512
beta[4]  0.5233562 0.02580370  0.4734307  0.5744832  804.6473 1.002729

MAP推定値の算出

切片と傾きについてのMAP推定値を計算しておきます。

#MAP estimate-----------------------------------------------------------------

#MAP estimate
MAPs  = c(rep(NA,c((Nvar+NY)*c(length(parameters2)))))
for (i in 1:c((Nvar+NY)*c(length(parameters2)))){
    MAPs[i]    <- density(df.ext[,i])$x[which(density(df.ext[,i])$y == max(density(df.ext[,i])$y))]
}
names(MAPs) <- c(colnames(df.ext))

こちらが結果になります。

MAPs
    beta.1     beta.2     beta.3     beta.4 
-0.2489881 -0.2102396  0.2255024  0.5223438 

収束診断と情報量基準の算出

shinystanで収束の確認をしましょう。DiagnoseからNeffやRhatのタブをクリックして、”None”が表示されていることを確認します。

launch_shinystan(fit)

LOOとWAICも一応計算しておきます。


#情報量規準WAICとLOOの算出
#WAIC and Leave one out cross validation--------------------------------------------
WAIC<-c()
for (i in 1:NY){
    res<-waic(extract_log_lik(fit)[,c(Nsub*(i-1)+1):c(i*Nsub)])
    WAIC[i]<-res$waic
}

#Leave-One-Out cross validation
LOO<-c()
for (i in 1:NY){
    res2<-loo(extract_log_lik(fit)[,c((Nsub)*(i-1)+1):c(i*Nsub)])
    LOO[i]<-res2$looic
}

結果です。

print(WAIC)
[1] -62.72323

 print(LOO)
[1] -62.67738

推定結果(事後分布)の可視化

推定された傾きと切片の事後分布を可視化します。コードはやや複雑ですが、ここはこういうものだ、というかたちで結構です。

やっていることは、MCMCサンプルをとりだして、データの形を整えてから、ggplot2でヒストグラムを書くという作業です。

#パラメータのggplot----------------------------------------

#Extract
d.ext<-rstan::extract(fit)    
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', 'varNames'))

#replace Num of parameter to varNames
for(i in 1: c(Nvar+NY)){
    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)

#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
)
)
g <- g + geom_histogram(position = "identity",alpha = 0.5)
g <- g + geom_density(stat = "density",position = "identity",alpha = 0.4)
#facet
g<-g+facet_wrap(~parameter,scales="free")

#Colour Variations
#Sequential(連続カラー): Blues BuGn BuPu GnBu Greens Greys Oranges OrRd PuBu PuBuGn PuRd Purples RdPu Reds YlGn YlGnBu YlOrBr YlOrRd
#Discrete(離散カラー)  : BrBG PiYG PRGn PuOr RdBu RdGy RdYlBu RdYlGn Spectral
g<-g + scale_color_brewer(palette="PuRd")
g<-g + scale_fill_brewer(palette="PuRd")
g<-g+ggtitle("Posterior Distributions")

#axis
g <- g + theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) 
g <- g + theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))   
g <- g + 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
g<-g+geom_vline(xintercept = 0)
g


#if you need interactive plot
ggplotly(g)

よいしょ。

ゼロところに縦線を入れておきました。

こうすると、95%信用区間が0を超えているか否か。傾きに意味があるのかどうなのか、わかりやすいですね。

視覚的事後予測チェック

モデルによる予測が妥当かどうか、確かめる方法はいくつかあります。 その中の一つが視覚的事後予測チェック。推定された事後分布から乱数を生成し、実データとの重なりを見るものです。

方法は単純です。やりかたはいくつかあると思いますが、ここでは次の手順で行っていきます。

  1. 推定された切片や傾き、誤差の事後分布から乱数を生成し、予測分布のデータをつくる

  2. 事後分布からの乱数は、サンプリングの上から9つ(カラーパレットの色の上限が9だったから)をとりだし、それぞれを1回ずつ使用する

  3. 実データ(従属変数)のヒストグラムをかく

  4. 実データのヒストグラムに、予測分布の確率密度曲線(近似)を重ねてかく

#Posterior Predictive Checks---------------------------------------------------
#参加者数
Nsub<-Nsub

#何iterで事後予測チェックするか
repN<-9

#事後分布サンプルからの乱数生成##########################################
predY<-array(NA,c(repN,Nsub,NY))
for(i in 1:repN){
    for(j in 1: NY){
        #modelにあわせて変更--------------------------------------
        predY[i,,j] <-rnorm(Nsub,d.ext$mu[i,],d.ext$sigma[i])
        #---------------------------------------------------------
    }
}
#########################################################################


#ggplotようにデータ成形
#reshape "predictiveData"  to tidy data
m.pDat<-melt(predY,idvar=c())

#replace the parameters' number to varName 
for(i in 1: c(Nvar+1)){
    m.pDat$Var3<-replace(m.pDat$Var3,m.pDat$Var3==as.character(i), varY[i]) 
}
colnames(m.pDat)<-c("pred_y","Nsub","varNames","value")
m.pDat <- m.pDat[order(m.pDat$varNames),]

#reshape "rawData" to long data
m.Dat<-melt(Y,idvar=c())
colnames(m.Dat)<-c("value")


#ggplot---------------------
#data histogram
p<-ggplot(m.Dat,aes(x=value))+ 
    geom_histogram(aes(y=..density..
                      # ,fill=value#varNames
                      # ,colour="grey"#varNames
                       ),alpha=0.9)

#theme and graph title size
p<-p+theme_bw(base_size=18,base_family="HiraKakuProN-W3") # Graph Title Size and Font
p<-p+ggtitle("Posterior Predictive Checks")


#predictive data line added
m.pDat$pred_y<-as.factor(m.pDat$pred_y)
p<-p+geom_line(data=m.pDat,aes(x=value,y=..density..,
                               colour=pred_y
                               #,group=pred_y
                               ),stat="density")

#facet
p<-p+facet_wrap(~varNames,scales = "free")

#axis set 
p <- p + theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) 
p <- p + theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))   
p <- p + 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

#Colour Variations
#Sequential(連続カラー): Blues BuGn BuPu GnBu Greens Greys Oranges OrRd PuBu PuBuGn PuRd Purples RdPu Reds YlGn YlGnBu YlOrBr YlOrRd
#Discrete(離散カラー)  : BrBG PiYG PRGn PuOr RdBu RdGy RdYlBu RdYlGn Spectral
p<-p+scale_color_brewer(palette="PuRd") #Scale Fill

p
#if you need interactive plot
ggplotly(p)

ぃぃぃよいしょっと

まぁ、いいかんじですかね。二峰性の分布をちゃんと表現できています。

最後に

今回使用したコードを最後にまとめて記載しておきます。

それでは、EnjoyStan!

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

#library and directory
library(rstan)
library(ggplot2)
library(shinystan)
library(reshape2)
library(loo)
library(tidyverse) 
library(plotly)
library(Rmisc)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#
# 目的: 線形モデルで傾きを推定する
#

#Options-----------------------------------------------------
# 1. WAIC&LOO: 情報量規準WAIC&LOOの算出
# 2. Parameters Visualization of beta
# 3. MAP estimate 
# 4. Posterior Predictive Check
#------------------------------------------------------------


#Initaial Settings-------------------------------------------
#サンプリング数
iterations =c(2000)
#Warm-up数
warmup=c(1000)
#チェーン数
chains=c(3)
#thinning
thin =c(1)

#推定に使用するパラメータ
parameters=c("beta","sigma","mu","log_lik") #取り出したいパラメータを最初に並べる
#あとで結果を確認したいパラメータ
parameters2<-c("beta")#推定に使用するパラメータの先頭から順になるように

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


#線形モデルの独立変数
varNames<-c("Sepal.Length","Sepal.Width","Petal.Length")

#線形モデルの従属変数
varY<-c("Petal.Width")


#Data shaping
Dat1<-data.matrix(dat[,varNames])#independent var
Y<-dat[,varY] #dependent var

#------------------------------------------------------------







#線形代数を使用するため、一列目の変数名を変更
varNames<-c("Intercept",varNames)

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

#線形代数で解くので,切片用に1の列を作成して1列目に結合
InterceptDat<-c(rep(1,Nsub))
Dat<-cbind(InterceptDat,Dat1)

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





LinearModel<-"
data {
    int <lower=0>Nvar;
    int <lower=0>Nsub;
    int <lower=0>NY;
    vector [Nsub]Y;
    matrix [Nsub,Nvar+NY] Dat;
}

parameters {
    vector [Nvar+NY]beta;
    real <lower=0>sigma;
}

transformed parameters {
    vector [Nsub]mu;
    mu = Dat*beta;
}

model {
    Y~ normal(mu,sigma);
}

generated quantities {
    matrix [Nsub,NY] log_lik;
    for (i in 1 : Nsub){
        for (j in 1 : NY){
            log_lik[i,j] = normal_lpdf(Y[i]|mu[i],sigma);
        }
    }
}

"

#sampling
fit<-stan(model_code=LinearModel,data=datastan,pars=parameters,
          iter = iterations,chains = chains,thin = thin)

#Summary
Summary<-summary(fit)$summary[c(1:c((Nvar+NY)*c(length(parameters2)))),c(1,3,4,8,9,10)]
print(Summary)




#情報量規準WAICとLOOの算出
#WAIC and Leave one out cross validation--------------------------------------------
WAIC<-c()
for (i in 1:NY){
    res<-waic(extract_log_lik(fit)[,c(Nsub*(i-1)+1):c(i*Nsub)])
    WAIC[i]<-res$waic
}


#Leave-One-Out cross validation
LOO<-c()
for (i in 1:NY){
    res2<-loo(extract_log_lik(fit)[,c((Nsub)*(i-1)+1):c(i*Nsub)])
    LOO[i]<-res2$looic
}




#パラメータのggplot----------------------------------------

#Extract
d.ext<-rstan::extract(fit)    
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', 'varNames'))

#replace Num of parameter to varNames
for(i in 1: c(Nvar+NY)){
    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)

#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
)
)
g <- g + geom_histogram(position = "identity",alpha = 0.5)
g <- g + geom_density(stat = "density",position = "identity",alpha = 0.4)
#facet
g<-g+facet_wrap(~parameter,scales="free")

#Colour Variations
#Sequential(連続カラー): Blues BuGn BuPu GnBu Greens Greys Oranges OrRd PuBu PuBuGn PuRd Purples RdPu Reds YlGn YlGnBu YlOrBr YlOrRd
#Discrete(離散カラー)  : BrBG PiYG PRGn PuOr RdBu RdGy RdYlBu RdYlGn Spectral
g<-g + scale_color_brewer(palette="PuRd")
g<-g + scale_fill_brewer(palette="PuRd")
g<-g+ggtitle("Posterior Distributions")

#axis
g <- g + theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) 
g <- g + theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))   
g <- g + 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

g<-g+geom_vline(xintercept = 0)
g
#if you need interactive plot
#ggplotly(g)





#MAP estimate-----------------------------------------------------------------

#MAP estimate
MAPs  = c(rep(NA,c((Nvar+NY)*c(length(parameters2)))))
for (i in 1:c((Nvar+NY)*c(length(parameters2)))){
    MAPs[i]    <- density(df.ext[,i])$x[which(density(df.ext[,i])$y == max(density(df.ext[,i])$y))]
}
names(MAPs) <- c(colnames(df.ext))
MAPs





#Posterior Predictive Checks---------------------------------------------------
#参加者数
Nsub<-Nsub

#何iterで事後予測チェックするか
repN<-9

#事後分布サンプルからの乱数生成##########################################
predY<-array(NA,c(repN,Nsub,NY))
for(i in 1:repN){
    for(j in 1: NY){
        #modelにあわせて変更--------------------------------------
        predY[i,,j] <-rnorm(Nsub,d.ext$mu[i,],d.ext$sigma[i])
        #---------------------------------------------------------
    }
}
#########################################################################

#reshape "predictiveData" array to long data
m.pDat<-melt(predY,idvar=c())

#replace the parameters' number to varName 
for(i in 1: c(Nvar+1)){
    m.pDat$Var3<-replace(m.pDat$Var3,m.pDat$Var3==as.character(i), varY[i]) 
}
colnames(m.pDat)<-c("pred_y","Nsub","varNames","value")
m.pDat <- m.pDat[order(m.pDat$varNames),]

#reshape "rawData" to long data
m.Dat<-melt(Y,idvar=c())
colnames(m.Dat)<-c("value")




#ggplot---------------------
#data histogram
p<-ggplot(m.Dat,aes(x=value))+ 
    geom_histogram(aes(y=..density..
                      # ,fill=value#varNames
                      # ,colour="grey"#varNames
                       ),alpha=0.9)

#theme and graph title size
p<-p+theme_bw(base_size=18,base_family="HiraKakuProN-W3") # Graph Title Size and Font
p<-p+ggtitle("Posterior Predictive Checks")


#predictive data line added
m.pDat$pred_y<-as.factor(m.pDat$pred_y)
p<-p+geom_line(data=m.pDat,aes(x=value,y=..density..,
                               colour=pred_y
                               #,group=pred_y
                               ),stat="density")

#facet
p<-p+facet_wrap(~varNames,scales = "free")

#axis set 
p <- p + theme(axis.title.x = element_text(size=14), #axis title letter size
               axis.title.y = element_text(size=14)) 
p <- p + theme(axis.text.x = element_text(size=8),  #axis text letter size 
               axis.text.y = element_text(size=8))   
p <- p + 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

#Colour Variations
#Sequential(連続カラー): Blues BuGn BuPu GnBu Greens Greys Oranges OrRd PuBu PuBuGn PuRd Purples RdPu Reds YlGn YlGnBu YlOrBr YlOrRd
#Discrete(離散カラー)  : BrBG PiYG PRGn PuOr RdBu RdGy RdYlBu RdYlGn Spectral
p<-p+scale_color_brewer(palette="PuRd") #Scale Fill

p
#if you need interactive plot
#ggplotly(p)



#Final----------------------------------------------------------
# Print and Writing Results
AnalysisTime <- Sys.time() 
Results<-list(AnalysisTime=AnalysisTime,
              iterations=iterations,chains=chains,warmup=warmup,thin=thin,
              varY=varY,varNames=varNames,Summary=Summary,MAPs=MAPs)
InformationCriteria<-list(WAIC=WAIC,LOO=LOO)

print(Results)
print(InformationCriteria)

#All Graphics
multiplot(g,p, cols=1)

#writing 
write.csv(list(Results=Results),"LinearModel.csv")
ggsave("LinearModel.png",g)
ggsave("LinearModelPPcheck.png",p)
ggsave("LinearModelFullImages.png",multiplot(g,p, cols=1))


Written on January 13, 2017