『マーケティングの統計モデル』第8章の階層ベイズモデルをStanでモデリングしてみた

はじめに
階層ベイズモデルによるブランド選択行動のモデル化

はじめに

『マーケティングの統計モデル』の第8章で登場する階層ベイズモデル (階層ベイズモデルによるブランド選択行動のモデル化) を Stan でモデリングしてみたという記事になります。

同書ではサポートページ経由で R によるサンプルコードも提供されていますが、第8章のベイズモデリングに関しては bayesm や dlm パッケージを用いたモデリングになっており、Stan を用いたコードはありません。

そこで、本稿では Stan を用いて同書の分析を再現してみたいと思います。なお、分析は Google Colaboratory 上で Pystan を用いて行いました。

マーケティングの統計モデル (統計解析スタンダード)

作者:佐藤忠彦
朝倉書店

Amazon

階層ベイズモデルによるブランド選択行動のモデル化

分析内容

同書で取り上げられている階層ベイズモデルによる分析は、同書4章で取り上げられたロジットモデルによるブランド選択行動のモデル化をさらに発展させてたものになります。

まず初めに上記のロジットモデルによる分析について簡単に触れたいと思います。同書4章では3つの醤油の商品(ナショナルブランド2つ、プライベートブランド1つ) の ID 付 POS データを用いて、ブランド選択行動のモデル化を行っています。ID 付 POS データのデータ構造は以下の通りです。

パネル番号
購買日付
商品1選択の有無 (1選択0非選択)
商品2選択の有無 (1選択0非選択)
商品3選択の有無 (1選択0非選択)
選択商品番号 (1or2or3): Y
商品1価格掛率: X1
商品2価格掛率: X2
商品3価格掛率: X3
商品1山積み陳列実施の有無(1実施0非実施): X4
商品2山積み陳列実施の有無(1実施0非実施): X5
商品3山積み陳列実施の有無(1実施0非実施): X6
商品1チラシ掲載の有無(1掲載0非掲載): X7
商品2チラシ掲載の有無(1掲載0非掲載): X8
商品3チラシ掲載の有無(掲載実績なし．すべて0）

このモデルをさらに発展させて、同書第8章では階層ベイズモデルを用いることによってパネルごとの消費者異質性を考慮したモデル化を行っています。その際に以下のようなパネルごとの属性データを連結させ、ブランド選択に影響する各説明変数の係数が個々のパネルによって異なるというモデルを構築しています。

パネル番号
年齢: Z1
家族人数: Z2

つまり、醤油のブランド選択行動に対する各商品の価格掛率やプロモーション (山積み陳列 / チラシ掲載) の影響度合いを分析する、さらにその際その影響度合いが個々のパネルごとに異なるであろうということを踏まえてその消費者異質性を評価するという形になります。以下では、それを Stan を用いてモデリングしていきたいと思います。

モデル式の記述

上記をモデル式にすると以下のようになります。これを Stan のコードに落としていきます。

$\displaystyle{\eta_{i} \sim Normal(0, s_{b})}$

$\displaystyle{b_{1}[m] = \theta_{1,1}Z_{1}[m] + \theta_{1,2}Z_{2}[m] + \eta_{1}}$

$\displaystyle{b_{2}[m] = \theta_{2,1}Z_{1}[m] + \theta_{2,2}Z_{2}[m] + \eta_{2}}$

$\displaystyle{b_{3}[m] = \theta_{3,1}Z_{1}[m] + \theta_{3,2}Z_{2}[m] + \eta_{3}}$

$\displaystyle{b_{4}[m] = \theta_{4,1}Z_{1}[m] + \theta_{4,2}Z_{2}[m] + \eta_{4}}$

$\displaystyle{b_{5}[m] = \theta_{5,1}Z_{1}[m] + \theta_{5,2}Z_{2}[m] + \eta_{5}}$

$\displaystyle{U_{1}[n] = b_{1}[N2ID[n]]log(X_{1}[n]) + b_{2}[N2ID[n]]X_{4}[n] + b_{3}[N2ID[n]]X_{7}[n] + b_{4}[N2ID[n]]}$

$\displaystyle{U_{2}[n] = b_{1}[N2ID[n]]log(X_{2}[n]) + b_{2}[N2ID[n]]X_{5}[n] + b_{3}[N2ID[n]]X_{8}[n] + b_{5}[N2ID[n]]}$

$\displaystyle{U_{3}[n] = b_{1}[N2ID[n]]log(X_{3}[n]) + b_{2}[N2ID[n]]X_{6}[n]}$

$\displaystyle{Y[n] \sim Categorical(softmax(\vec{U[n]}))}$

Stan コード

上記の階層ベイズモデルを Stan コードで再現したのが以下になります。

data {
    int N;  // POS データのレコード件数
    int K;  // 離散選択肢の数 = ブランド数
    int D;  // 共変量の次元
    matrix [N, D] X;  // 共変量行列
    int<lower=1, upper=K> Y[N];  // 結果変数 = 選択したブランド

    int M;  // ユニークな ID 数
    int<lower=1, upper=M> N2ID[N];  // パネル n の属性テーブルにおけるインデックス
    int D2;  // パネルの属性に関する共変量の次元
    matrix[M, D2] Z;  // パネルの属性に関する共変量行列
}

parameters {
    vector<lower=0>[5] s_b;
    matrix[5, 2] theta;
    vector[5] eta;
}

transformed parameters {
    matrix[N, K] U;
    matrix[M, 5] b;

    for (m in 1:M) {
        b[m, 1] = theta[1, 1] * Z[m, 1] + theta[1, 2] * Z[m, 2] + eta[1];
        b[m, 2] = theta[2, 1] * Z[m, 1] + theta[2, 2] * Z[m, 2] + eta[2];
        b[m, 3] = theta[3, 1] * Z[m, 1] + theta[3, 2] * Z[m, 2] + eta[3];
        b[m, 4] = theta[4, 1] * Z[m, 1] + theta[4, 2] * Z[m, 2] + eta[4];
        b[m, 5] = theta[5, 1] * Z[m, 1] + theta[5, 2] * Z[m, 2] + eta[5];
    }

    for (n in 1:N) {
        U[n, 1] = b[N2ID[n], 1] * log(X[n, 1]) + b[N2ID[n], 2] * X[n, 4] + b[N2ID[n], 3] * X[n, 7] + b[N2ID[n], 4];
        U[n, 2] = b[N2ID[n], 1] * log(X[n, 2]) + b[N2ID[n], 2] * X[n, 5] + b[N2ID[n], 3] * X[n, 8] + b[N2ID[n], 5];
        U[n, 3] = b[N2ID[n], 1] * log(X[n, 3]) + b[N2ID[n], 2] * X[n, 6];            
    }
}

model {
    eta ~ normal(0, s_b);

    for (n in 1:N) {
        Y[n] ~ categorical_logit(U[n, ]');
    }
}