廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

ディリクレ・多項分布回帰による margarine データの分析

分析対象

R の bayesm パッケージに入っている margarine データを分析します。

このデータは家計ごとのマーガリンの購買が記録されている。

たとえば購買されたマーガリンのブランド(1〜10)ごとに家庭の収入の分布をみるとこんな感じ。

f:id:abrahamcow:20180804054450p:plain

マーガリンのブランド選択に影響を与える変数はなにか、ということを調べたいと思います。

説明変数として使えそうなものには以下がある。

  • Income: 収入
  • Fam_Size: 家族の構成人数
  • college: 教育(世帯主が大学を出たか否か?)
  • whtcollar: 職業(世帯主がホワイトカラー労働者か否か?)
  • retired: 職業(世帯主が退職してるか否か?)

よくわからない変数もあるけど、まあいいことにする。ソースは
Allenby, Greg and Peter Rossi (1991), "Quality Perceptions and Asymmetric Switching Between Brands," Marketing Science 10, 185–205.
らしいので気になる人は読んでください。

詳しい分析例は以下にあります:

ここではもっと単純なモデルを試します。

モデル

家庭 m がマーガリン k を購買した回数を w_{m,k} とし、w_{m,k} を要素とする行列を W=(w_{m,k}) とします。

説明変数の行列を X (M行D列)で表し、W の m 行めのベクトルを w_m、X の m 行めのベクトルを x_m として、以下のモデルを考えます。

 \phi_m \sim \mathrm{Dirichlet}(\exp(x_m\beta)')
w_m \sim \mathrm{Multinomial}(\phi_m)

\beta は回帰係数の行列でD行K列。
マーガリンの選択される確率は  \phi_m で家庭ごとに異なります。

モデルのパラメータを推定するための Stan のコードは以下のようになりました。

data {
   int<lower=1> M;
   int<lower=1> K;
   int<lower=1> D;
   matrix[M,D] X;
   int<lower=0> W[M,K];
}
parameters {
  simplex[K] phi[M];
  matrix<lower=0>[D,K] beta;
}
model {
   for (m in 1:M){
    phi[m] ~ dirichlet(exp(X[m,]*beta)');
    W[m,] ~ multinomial(phi[m]);
   }
}

簡単です。

結果

MCMCの収束はトレースプロットの目視とRhat指標で確かめた。

観測された w_m/\sum_k w_{mk}\hat \phi(事後期待値)はこんな感じ。

f:id:abrahamcow:20180804061024p:plain

まあまあ当てはまってる。

次の図は  \hat \beta のヒートマップです。

各属性を持つ層がどのブランドを選択しやすいかわかります。

f:id:abrahamcow:20180804061403p:plain

R のコード

data(margarine,package = "bayesm")

library(tidyverse)
library(rstan)
rstan_options(auto_write = TRUE)
require(bayesm)

choice_hh <-margarine$choicePrice %>% 
  select(hhid,choice) %>% 
  left_join(margarine$demos) %>% 
  mutate(choice=factor(choice))

ggplot(choice_hh,aes(x=Income))+
  geom_histogram(bins=20,colour="black",fill="white")+
  facet_wrap(~choice,scales = "free_y")

choice_s <-group_by(choice_hh, hhid, choice, Income, Fam_Size, college, whtcollar, retired) %>%
  tally() %>% 
  spread(choice,n,fill=0) %>% 
  ungroup() %>% 
  mutate(Income=Income/10)

W <-as.matrix(choice_s[,-c(1:6)])
X <-cbind(Intercept=1,as.matrix(choice_s[,2:6]))

DM <-stan_model("DM.stan")
dat4stan <- list(W=W,X=X,D=ncol(X),K=ncol(W),M=nrow(W))
outDM <- sampling(DM,dat4stan,seed=123,cores=4,chain=4,iter=1000)
#1215.21 seconds (Total)

traceplot(outDM)
all(summary(outDM)$summary[,"Rhat"]<1.1)

phi <- matrix(get_posterior_mean(outDM,"phi"),dat4stan$M,dat4stan$K,byrow = TRUE)
Wdf <-as.data.frame(W/rowSums(W)) %>% 
  mutate(row=row_number()) %>% 
  gather(col,value,-row) %>% 
  mutate(type="observed") %>% 
  mutate(col=factor(as.integer(col)))

phidf <-as.data.frame(phi) %>% 
  set_names(1:10) %>% 
  mutate(row=row_number()) %>% 
  gather(col,value,-row) %>% 
  mutate(type="fitted") %>% 
  mutate(col=factor(as.integer(col)))

jointdf <-bind_rows(Wdf,phidf)

ggplot(jointdf,aes(x=col,y=row,fill=value))+
  geom_tile()+
  scale_fill_continuous(low="white",high="orange")+
  facet_wrap(~type)

all(summary(outDM)$summay[,"Rhat"]<1.1)

betadf <-matrix(get_posterior_mean(outDM,pars="beta")[,"mean-all chains"],dat4stan$D,dat4stan$K,byrow = TRUE)[-1,] %>% 
  as.data.frame() %>%
  set_names(1:10) %>% 
  mutate(name=colnames(X[,-1])) %>% 
  gather(category,value,-name) %>% 
  mutate(category=factor(as.integer(category)))

ggplot(betadf,aes(x=category,y=name,fill=value))+
  geom_tile(colour="black")+
  scale_fill_continuous(low="white",high="orange")+
  theme_bw(base_size=16)

summary(outDM,pars="beta")$summary