ほくそ笑む

R言語と統計解析について

実践 統計モデリング入門 【4. 三変数への拡張】

本シリーズの概要については、下記記事を最初に読んでください。

前回の記事では、2 種類のブランドから、消費者がどちらのブランドを選択するかという行動のモデル化を行いました。

今回は、ブランドを 3 種類に増やします。

三変数の場合のモデル

データとしては次のようなものを考えます。

y x1 x2 x3
2 58 53 51
3 48 34 61
2 56 50 46
1 51 51 60
3 54 52 69

y = 1 はブランド 1 を選択、y = 2 はブランド 2 を選択し、y = 3 はブランド 3 を選択したことを表します。
x1 はブランド 1 の商品の値引き額、x2 はブランド 2 の商品の値引き額、x3 はブランド 3 の商品の値引き額です。

このとき、次のようなモデルが考えられます。

  • 消費者は、より高い効用を得ることができるブランドを選択する
  • 効用は値引き額が大きいほど高くなる
  • 各ブランドは値引き額によらない固定された効用を持つ
  • 効用は消費者のその時の気分によってばらつきがある

これは、ブランドが 3 種類あるということ以外、前回のものと全く同じです。

これを数式で表すと、次のようになります。
任意の時点  t において、

  • ブランド 1 に対する効用 U_{1t}=x_{1t} + \varepsilon_{1t}
  • ブランド 2 に対する効用 U_{2t}=\alpha_2 + x_{2t} + \varepsilon_{2t}
  • ブランド 3 に対する効用 U_{3t}=\alpha_3 + x_{3t} + \varepsilon_{3t}
  •  \varepsilon_{1t}, \varepsilon_{2t},\varepsilon_{3t} \sim Norm(0, \sigma)

このとき、

  •  y=1 ( U_{1t} \geq U_{2t} かつ  U_{1t} \geq U_{3t} の場合)
  •  y=2 ( U_{2t} > U_{1t} かつ  U_{2t} \geq U_{3t} の場合)
  •  y=3 ( U_{3t} > U_{1t} かつ  U_{3t} > U_{2t} の場合)

このモデルにおける推定したいパラメータは、ブランド 2 および 3 の固定効用  \alpha_2, \alpha_3 と、気分によるばらつき  \sigma です。

ここまで表現できれば、あとは Stan によるモデル記述に移ります。

Stan によるモデルの記述

例のごとく、今回も increment_log_prob() を使用するため、上記モデルの対数尤度を求めていきます。

まず、 y_t=1 の場合、尤度は  U_{1t} \geq U_{2t} かつ  U_{1t} \geq U_{3t} となる確率と同じです。
ここで、各ブランド i に対して、 v_{it} = E[U_{it}] を定義します。ただし、 E[X] は確率変数  X の期待値です。
このとき、

  •  P(U_{1t} \geq U_{2t} \wedge U_{1t} \geq U_{3t})
  •  = P(U_{1t} \geq U_{2t})P(U_{1t} \geq U_{3t})
  •  = P(U_{1t} - U_{2t} \geq 0)P(U_{1t} - U_{3t} \geq 0)
  •  = (1-cdf_{12t}(0))(1-cdf_{13t}(0))

となります。ただし、 cdf_{ijt} は、 Norm(v_{it} - v_{jt}, \sqrt{2} \sigma) の累積分布関数です。
したがって、対数尤度は、

  •  LogLik_i = \log(1-cdf_{12t}(0)) + \log(1-cdf_{13t}(0))

となります。
同様に、 y_t=2 の場合、尤度は  U_{2t} > U_{1t} かつ  U_{2t} \geq U_{3t} となる確率ですので、

  •  LogLik_i = \log(1-cdf_{21t}(0)) + \log(1-cdf_{23t}(0))

 y_t=3 の場合、尤度は  U_{3t} > U_{1t} かつ  U_{3t} > U_{2t} となる確率ですので、

  •  LogLik_i = \log(1-cdf_{31t}(0)) + \log(1-cdf_{32t}(0))

となります。

Stan によるモデルの記述は次のようになります。

// model3.stan
data {
  int<lower=0> N;
  real x1[N];
  real x2[N];
  real x3[N];
  int<lower=1, upper=3> y[N];
}
parameters {
  real alpha2;
  real alpha3;
  real<lower=0> sigma;
}
model {
  real v1;
  real v2;
  real v3;
  for(i in 1:N) {
    v1 <- x1[i];
    v2 <- alpha2 + x2[i];
    v3 <- alpha3 + x3[i];
    if(y[i]==1) {
      increment_log_prob(normal_ccdf_log(0.0, v1 - v2, sqrt(2) * sigma));
      increment_log_prob(normal_ccdf_log(0.0, v1 - v3, sqrt(2) * sigma));
    } else if(y[i]==2) {
      increment_log_prob(normal_ccdf_log(0.0, v2 - v1, sqrt(2) * sigma));
      increment_log_prob(normal_ccdf_log(0.0, v2 - v3, sqrt(2) * sigma));
    } else {
      increment_log_prob(normal_ccdf_log(0.0, v3 - v1, sqrt(2) * sigma));
      increment_log_prob(normal_ccdf_log(0.0, v3 - v2, sqrt(2) * sigma));
    }
  }
  alpha2 ~ normal(0, 1000);
  alpha3 ~ normal(0, 1000);
  sigma ~ cauchy(0, 5);
}

model ブロックの中で中間変数 v1, v2, v3 を定義することで式を簡潔に書くことができます。
以上により、モデルは完成しました。
次は、このモデルにシミュレーションデータを当てはめてパラメータをうまく推定できるかどうかを確認します。

シミュレーションデータの作成

Stan へ入力するためのシミュレーションデータを作成します。

generate_data <- function(N) {
  alpha2 <-  10 # parameter
  alpha3 <- -10 # parameter
  sigma  <-   5 # parameter
  x1 <- rpois(N, lambda = 50)
  x2 <- rpois(N, lambda = 45)
  x3 <- rpois(N, lambda = 55)
  u1 <-          x1 + rnorm(N, 0, sigma)
  u2 <- alpha2 + x2 + rnorm(N, 0, sigma)
  u3 <- alpha3 + x3 + rnorm(N, 0, sigma)
  y <- ifelse(u1 >= u2 & u1 >= u3, 1, ifelse(u2 > u1 & u2 >= u3, 2, 3))
  data.frame(y=y, x1=x1, x2=x2, x3=x3)
}
data <- generate_data(N = 1000)
head(data, 5)
  y x1 x2 x3
1 1 55 30 56
2 2 50 44 59
3 3 50 38 65
4 1 62 47 61
5 2 51 46 53

真のパラメータ \alpha_2 = 10\alpha_3 = -10 \sigma = 5 のデータができました。
それでは、このデータに対して先程のモデルを当てはめてみましょう。

モデルの当てはめ

Stan コードにより記述されたモデルをデータに当てはめるための R コードは次のようになります。

library(rstan)

standata <- list(N=nrow(data), x1=data$x1, x2=data$x2, x3=data$x3, y=data$y) # データの整形
model <- stan_model(file = "model3.stan")                                    # Stan コードのコンパイル
result <- sampling(model, data = standata, chains = 3)                       # MCMC サンプリングの実行

# 結果の表示
traceplot(result, inc_warm=FALSE)
print(result, digits = 2)

これを実行し、パラメータの推定結果を見てみます。


          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
alpha2   10.43    0.01 0.37    9.71   10.18   10.44   10.68   11.20  1517    1
alpha3   -9.87    0.01 0.44  -10.73  -10.18   -9.88   -9.58   -8.97  1310    1
sigma     4.89    0.01 0.22    4.48    4.74    4.88    5.03    5.33  1404    1
lp__   -638.88    0.04 1.22 -642.11 -639.45 -638.60 -638.00 -637.48  1042    1

Rhat < 1.1 であり、パラメータはちゃんと定常分布に収束していることがわかります。
alpha2 = 10、alpha3 = -10 および sigma = 5 がうまく推定できていることが確認できました。

まとめ

ブランドを 3 種類に増やした場合について、Stan によるモデルの記述を行い、シミュレーションデータに対してうまくパラメータ推定できることを確認しました。
これまでのモデルでは、全てのブランドに対して効用のばらつきの標準偏差  \sigma は同じという仮定を置いていましたが、次は標準偏差が異なる場合を考えます。