ほくそ笑む

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

実践 統計モデリング入門 【5. 異なる標準偏差への拡張】

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

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

これまでのモデルでは、全てのブランドに対して効用のばらつきの標準偏差  \sigma は同じという仮定を置いていました。
しかし、実際はこのような仮定は置けません。
そこで、今回は標準偏差がそれぞれ異なる場合を考えます。

異なる標準偏差に対するモデル

データは前回と同じです。

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 の商品の値引き額です。

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

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

これを数式で表すと、次のようになります。
任意の時点  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} \sim Norm(0, \sigma_1)
  •  \varepsilon_{2t} \sim Norm(0, \sigma_2)
  •  \varepsilon_{3t} \sim Norm(0, \sigma_3)

このとき、

  •  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_1, \sigma_2, \sigma_3 です。

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

Stan によるモデルの記述

今回も increment_log_prob() を使用するため、上記モデルの対数尤度を求めていきます。
といっても、前回との違いは標準偏差だけなので、あまり変わりません。

 y_t=1 の場合、

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

 y_t=2 の場合、

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

 y_t=3 の場合、

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

となります。ただし、 cdf_{ijt} は、 Norm(v_{it} - v_{jt}, \sqrt{\sigma_i^2 + \sigma_j^2}) の累積分布関数、 v_{it} U_{it} の期待値です。

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

// model4.stan
functions {
  real add_sigma(real sigma1, real sigma2) {
    return sqrt(pow(sigma1, 2) + pow(sigma2, 2));
  }
}
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> sigma1;
  real<lower=0> sigma2;
  real<lower=0> sigma3;
}
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, add_sigma(sigma1, sigma2)));
      increment_log_prob(normal_ccdf_log(0.0, v1 - v3, add_sigma(sigma1, sigma3)));
    } else if(y[i] == 2) {
      increment_log_prob(normal_ccdf_log(0.0, v2 - v1, add_sigma(sigma2, sigma1)));
      increment_log_prob(normal_ccdf_log(0.0, v2 - v3, add_sigma(sigma2, sigma3)));
    } else {
      increment_log_prob(normal_ccdf_log(0.0, v3 - v1, add_sigma(sigma3, sigma1)));
      increment_log_prob(normal_ccdf_log(0.0, v3 - v2, add_sigma(sigma3, sigma2)));
    }
  }
  alpha2 ~ normal(0, 1000);
  alpha3 ~ normal(0, 1000);
  sigma1 ~ cauchy(0, 5);
  sigma2 ~ cauchy(0, 5);
  sigma3 ~ cauchy(0, 5);
}

functions ブロックに再生性を考慮した標準偏差の加算関数 add_sigma() を作り、それを使うことで記述が簡略化されます。
以上により、モデルは完成しました。
次は、このモデルにシミュレーションデータを当てはめてパラメータをうまく推定できるかどうかを確認します。

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

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

set.seed(1)
generate_data <- function(N) {
  alpha2 <-  10 # parameter
  alpha3 <- -10 # parameter
  sigma1 <-   1 # parameter
  sigma2 <-   3 # parameter
  sigma3 <-   5 # parameter
  x1 <- rpois(N, lambda = 50)
  x2 <- rpois(N, lambda = 45)
  x3 <- rpois(N, lambda = 55)
  u1 <-          x1 + rnorm(N, 0, sigma1)
  u2 <- alpha2 + x2 + rnorm(N, 0, sigma2)
  u3 <- alpha3 + x3 + rnorm(N, 0, sigma3)
  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 2 47 47 55
2 2 44 45 55
3 3 37 52 66
4 2 50 49 55
5 1 55 38 55

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

モデルの当てはめ

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

library(rstan)
library(pforeach)

standata <- list(N=nrow(data), x1=data$x1, x2=data$x2, x3=data$x3, y=data$y) # データの整形
model <- stan_model("model4.stan")                                   # Stan コードのコンパイル
pforeach(i=1:3, .cores=3, .final=sflist2stanfit)({
  sampling(model, data = standata, chains = 1, seed=i)                               # MCMC サンプリングの実行
}) -> result

traceplot(result, inc_warm=FALSE)
print(result, digits = 2)

パラメータが増えて MCMC サンプリングに時間がかかるようになっため、pforeach パッケージを使って chain ごとに並列計算しています。

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


          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
alpha2   10.27    0.01 0.26    9.77   10.10   10.26   10.44   10.78  1353    1
alpha3   -9.33    0.01 0.41  -10.12   -9.61   -9.31   -9.04   -8.55  1235    1
sigma1    0.77    0.02 0.54    0.02    0.33    0.69    1.12    1.97   826    1
sigma2    3.14    0.01 0.30    2.57    2.95    3.15    3.34    3.71   976    1
sigma3    5.26    0.01 0.42    4.48    4.97    5.26    5.54    6.11  1061    1
lp__   -434.32    0.07 1.72 -438.40 -435.15 -434.00 -433.09 -432.10   674    1

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

さらなる拡張へ

標準偏差が異なる場合について、Stan によるモデルの記述を行い、シミュレーションデータに対してうまくパラメータ推定できることを確認しました。

以上により、今回の一連の記事で目的となる最終的なモデルを得ることができました。

ただし、シンプルなモデルからステップバイステップで拡張してきたことにより、さらなるモデルの拡張のアイデアがいくつか出てきたはずです。
ここでは、いくつかのモデル拡張アイデアを見てみましょう。

  • チラシに載ったかどうかの変数を追加する。どの牛乳を買うかは、折り込みチラシ等に載ったかどうかに強く影響を受けるかもしれません。それぞれのブランドがチラシに掲載されたかどうかのデータを確認して、利用可能ならばそれを変数として追加します。このとき、値引き額による効果とチラシによる効果は同等ではないため、それぞれの係数をパラメータとして推定する必要があります。また、値引き額が大きいとチラシに掲載されやすくなるという関係も考えられるため、交互作用項を含める必要があるかもしれません。
  • 個人に対する情報を含める。小売店ではポイントカードなどにより個人を特定できる情報を持っている場合があります。個々人に対するブランド効用や値引き効用は異なると考えられるため、このような情報を使って個人に対する効用を推定することができます。これを行うために、個人効用のばらつきに仮定を加え、階層モデルにしてパラメータ推定を行うという方法が考えられます。
  • 潜在変数の追加。今回、値引き額の値は定価からの差額で表されました。しかし、いつも安売りしているブランドに対しては、消費者は定価よりも低い価格を期待します。したがって、定価ではなく平均的な価格との差の方が変数として適切かもしれません。お客さんが値引き額の基準として採用している価格ははっきりとは分かりませんが、直近一か月間の平均値などを採用して近似的に表現することは可能です。このように、直接的には観測されないが仮定的に表現される変数のことを潜在変数といいます。この場合、値段が基準価格よりも高くなることも考えられ、さらに高い場合のマイナス効果と安い場合のプラス効果が非対称になると考えられます。このような場合には、インディケータ変数を用意して高い場合と低い場合で項を別にして係数の推定を行います。

このような拡張をモデルに組み込むのは、最初は複雑すぎると感じるかもしれません。

しかし、この一連の記事で見てきたように、複雑な拡張のアイデアに対しても、ステップバイステップでモデリングを行い、シミュレーションで逐一確認することによって、同じようにモデルを構築することができます。

おわりに

この一連の記事では、統計モデリングに対して、

  • シンプルなモデルから始めて段階的にモデルを拡張していく
  • 構築したモデルに対してシミュレーションを行う

という 2 つのコツに基づいて、実際のモデル構築の手順を説明しました。

複雑なモデルを作りたいとき、いきなり最終的なモデルを作り始めるのではなく、段階的なモデリングとシミュレーションを繰り返すことが、モデル構築の近道だと思います。

みなさまのモデリングライフのご参考になれば幸いです。

それでは。