ほくそ笑む

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

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

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

前回の記事では、ブランドを 1 種類に限定し、消費者が購入するかどうかの基準をモデル化しました。

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

二変数の場合のモデル

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

y x1 x2
1 54 47
2 56 51
1 52 44
1 50 38
2 42 41

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

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

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

ここで効用と言っているのは、商品を購入したときの消費者の満足度のことです。
各ブランドの固定効用というのは、値段に関わらないところで決まる効用のことです。例えば、明治おいしい牛乳とメグミルクを比較するときに、どちらのブランドがいいイメージを持っているかなどで決まる効用のことを指します。

これを数式で表すと、次のようになります。

任意の時点  t において、

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

このとき、

  •  y_t = 1 ( U_{t1} \geq U_{t2} の場合)
  •  y_t = 2 ( U_{t1} < U_{t2} の場合)

このモデルに対する推定したいパラメータはブランド 2 の固定効用  \alpha_2 と、気分によるばらつき  \sigma です。
固定効用について、なぜブランド 2 だけを考えるかというと、出力変数 y の決定について、効用 U_1U_2 のどちらが大きいかということだけから決まるため、各効用の相対的な値しか推定できないためです。ただし、ブランド 1 を基準としてブランド 2 のイメージがどれくらい高い(低い)かが分かるため、これで十分と言えます。

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

Stan によるモデルの記述

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

まず、 y_t=1 の場合、尤度は  U_{1t} \geq U_{2t} となる確率と同じです。
また、 U_{1t} \sim Norm(x_{1t}, \sigma),  U_{2t} \sim Norm(\alpha2 + x_{2t}, \sigma) です。
このとき、 P(U_{1t} \geq U_{2t}) = P(U_{1t}-U_{2t} \geq 0) より、正規分布の再生性を使って、

  •  P(U_{1t} \geq U_{2t}) = 1 - cdf_t(0)

となります。ただし、 cdf_t Norm(x_{1t} - (\alpha_2 + x_{2t}), \sqrt{2} \sigma) の累積分布関数です。
したがって、対数尤度は、

  •  LogLik_t = \log(1 - cdf_t(0))

となります。

同様に、 y_t=2 の場合、尤度は  U_{1t} < U_{2t} となる確率と同じですので、

  •  LogLik_t = \log(cdf_t(0))

となります。

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

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

以上により、モデルは完成しました。
次は、このモデルにシミュレーションデータを当てはめてパラメータをうまく推定できるかどうかを確認します。

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

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

generate_data <- function(N) {
  alpha2 <- 10 # parameter
  sigma  <-  5 # parameter
  x1 <- rpois(N, lambda = 50)
  x2 <- rpois(N, lambda = 45)
  u1 <-          x1 + rnorm(N, 0, sigma)
  u2 <- alpha2 + x2 + rnorm(N, 0, sigma)
  y  <- ifelse(u1 >= u2, 1, 2)
  data.frame(y=y, x1=x1, x2=x2)
}
data <- generate_data(N = 1000)
head(data, 5)
  y x1 x2
1 2 41 51
2 1 51 28
3 2 50 44
4 1 55 41
5 1 50 44

真のパラメータ  \alpha_2 = 10 \sigma = 5 のシミュレーションデータができました。

それでは、このデータに対して先程のモデルを当てはめてみましょう。

モデルの当てはめ

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

library(rstan)

standata <- list(N=nrow(data), x1=data$x1, x2=data$x2, y=data$y) # データの整形
model <- stan_model(file = "model2.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.38    0.01 0.39    9.65   10.11   10.37   10.64   11.16  1189    1
sigma     5.01    0.01 0.29    4.49    4.81    4.98    5.19    5.64  1478    1
lp__   -371.16    0.03 1.02 -373.83 -371.55 -370.85 -370.46 -370.21   910    1

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

まとめ

ブランドを 2 種類に増やした場合について、Stan によるモデルの記述を行い、シミュレーションデータに対してうまくパラメータ推定できることを確認しました。
次はブランドを 3 種類に増やした場合について考えます。