rstanで個人のパラメーターを推定した話(JapanR2015のLT補足)

昨年の12/5に開催されたJapanRでLTをしました。

rstanで個人のパラメーターを推定した話 from Yuya Matsumura

www.slideshare.net

5分しかなかったのであんまりちゃんと説明できなかったのですが

要約すると以下のような感じです。

rstanで階層ベイズモデルを書いて推定した。
「個人の」というのは「一人一人の」という意味
パネルデータを用いた。
目的変数は「シャンプーのブランドスイッチ」、説明変数は「価格」、「CM放送数」、「そのブランドがSNSでつぶやかれた数」
階層構造の説明変数は「性別」「年齢」「未既婚ダミー」「世帯年収」

モデル式などはスライドを見てください。

問題なのはスライドではすべてのコードを載せられなかったことです。

これじゃいかんということですべてのRコードとStanコードを載せておきたいと思います。

まずはRコード

setwd("E:/TokyoR/JapanR2015")
library(rstan)

#stanファイル名
scr <- "standraft.stan"
d <- read.csv("sample.csv")

#ユニークな個人の情報
d1 <- d[!duplicated( data.frame(d$id) ), ]
#ファイル内の変数をリストにまとめる
data = list(N=nrow(d1),
            NT=nrow(d),
            y = d$離反,
            x = cbind(rep(1,length=nrow(d)),matrix(scale(d$価格),nrow=nrow(d), ncol=1)),
            w = cbind(
                      matrix(scale(d$ucm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$scm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$qcm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$kcm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$pcm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$usns),nrow=nrow(d), ncol=1),
                      matrix(scale(d$ssns),nrow=nrow(d), ncol=1),
                      matrix(scale(d$qsns),nrow=nrow(d), ncol=1),
                      matrix(scale(d$ksns),nrow=nrow(d), ncol=1),
                      matrix(scale(d$psns),nrow=nrow(d), ncol=1)
                    ),
            id = d$id,
            z = cbind(
　　　　　　　　　　　matrix(scale(d1$年齢),nrow=nrow(d1), ncol=1),
　　　　　　　　　　　matrix(scale(d1$性別),nrow=nrow(d1), ncol=1),
                      matrix(scale(d1$未既婚),nrow=nrow(d1), ncol=1),
                      matrix(scale(d1$世帯年収),nrow=nrow(d1), ncol=1),
                     )
            )


par <- c("alpha","beta") #サンプリング結果を出力する変数名
war <- 10 #バーンイン期間
ite <- 100 #繰り返し回数
see <- 1234 #乱数の種類
dig <- 3 #有効数字
cha <- 1 #連鎖構成数

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

#モデル
fit <- stan(file="mitaronstandraft.stan",
            data=data,
            pars=c("alpha","beta"),
            verbose=F,
            seed=71,
            chains=1,
            warmup=1000,
            iter=10000)


#結果の出力

fit.summary <- data.frame(summary(fit)$summary)
write.csv(fit.summary, "fit_summary1.csv", quote=F)


pdf("fit_plot1.pdf", width=600/72, height=600/72)
plot(fit)
dev.off()

pdf("fit_traceplot1.pdf", width=600/72, height=600/72)
traceplot(fit)
dev.off()

library("coda")
fit.coda<-mcmc.list(lapply(1:ncol(fit),function(x) mcmc(as.array(fit)[,x,])))
pdf("fit_density1.pdf", width=600/72, height=600/72)
plot(fit.coda) #事後分布をプロット
dev.off()