はじめに

前回は、ゼロ過剰ポアソンを用いてCM視聴がゲームのプレイ時間に与える効果の測定を行いました。ゲームのプレイの有無に与える効果の部分はうまく推定できましたが、ゲームをプレイする時間に与える効果はポアソン分布ではうまく推定することができませんでした。

この部分もうまく推定できるよう、階層ベイズやゼロ過剰負の二項分布等を試してみましたが、そちらはうまくいかず、最終的にゲームのプレイ時間をスケールすると普通の負の二項分布に当てはまりが良くなったので、その結果を書いていきます。

負の二項分布について

負の二項分布については下記の資料がわかりやすかったです。
https://www.slideshare.net/simizu706/ss-50994149

負の二項分布は、下記のようなある成功率pで、r回正解するまでに必要な失敗回数を示す分布で、過分散対策としてよく使われるようです。

$\operatorname { Pr } ( y | p , r ) = \left( \begin{array} { c } { y + r - 1 } \\ { y } \end{array} \right) ( p ) ^ { r } ( 1 - p ) ^ { y }$ 　　（1）

これをガンマ分布を用いて連続変量に対応できるようにすると

$\operatorname { Pr } ( y | p , r ) = \frac { \Gamma ( y + r ) } { y ! \Gamma ( r ) } ( p ) ^ { r } ( 1 - p ) ^ { y }$ 　　（2）

また、ポアソン分布の $\lambda$ が確率的に変動すると考えると下記のような式になります。

$\operatorname { Pr } ( y | \phi , \mu ) = \frac { \Gamma ( y + \phi ) } { y ! \Gamma ( \phi ) } \left( \frac { \phi } { \mu + \phi } \right) ^ { \phi } \left( \frac { \mu } { \mu + \phi } \right) ^ { y }$ 　　（3）

これは成功率 $\phi/(\phi + \mu)$ で、 $\phi$ 回正解するまでの失敗回数yの分布を示します。
stanのモデルではこの $\phi, \mu$ を用いた式を使用します。

式(1)の負の二項分布についてパラメータを変えていった場合、グラフは以下のようになります。

f:id:rmizutaa:20190218203354p:plain

今回は前回のようにゲームのプレイ時間を0と１以上で分割することはせず、目的変数のゲームプレイ時間をスケールさせることで、分布を負の二項分布にあてはめます。

モデル

コードは以下のようになります。前回との主な変更点は以下の２点です。
・モデルを負の二項分布に変更
・分布に合うようにゲームのプレイ時間を1/10000にスケール。スケールによって0以上のプレイ時間が0になるのを防ぐために値は切り上げ

df=pd.read_csv("q_data_x.csv")

#スケーリング
df["TVwatch_day"]=df["TVwatch_day"]/10000
#スケール＋切り上げ
df["gamesecond"]=[math.ceil(i/10000) for i in df["gamesecond"]]

X=df[['cm_dummy','T','F1', 'F2', 'F3', 'M1', 'M2', 'M3', 'TVwatch_day',"child_dummy"]]

model="""
data {
  int N;
  int D;
  int<lower=0> Y[N];
  matrix[N,D] X;
}

parameters {
  vector[D] b1;
  real<lower=0> beta;
}

transformed parameters {
  vector[N] alpha_x;
  vector[N] alpha;
  
  alpha_x = X*b1;
  for (n in 1:N){
    alpha[n] = exp(alpha_x[n]);
    }
}

model {
  for (n in 1:N)
    Y[n] ~ neg_binomial_2(alpha[n], beta);
}
"""
data={'N': X.shape[0], 'D':X.shape[1] ,'Y':df["gamesecond"].values, 'X':np.array(X)}
fit = pystan.stan(model_code=model, data=data, iter=1000, chains=3)

結果

収束したので、各変数がゲームのプレイ時間に与える効果を確認します。

負の二項分布では $\phi,\mu$ の２つのパラメータを推定していますが、この分布の平均は $\phi$ 、分散は $\phi/(\phi + \mu)$ になります。なので各変数の効果を確認するために、ゲームのプレイ時間の平均時間に影響を与えるパラメータ $\phi$ について見ていきます。

# サンプル列を抽出
la  = fit.extract(permuted=True)
#ゲームプレイ有無への寄与度
plt.figure(figsize=(12,8))
plt.violinplot([la["b1"][:,i] for i in range(10)])
plt.xticks(list(range(1,len(X.columns)+1)), X.columns)
plt.xticks(rotation=45)
plt.tick_params(labelsize=18)
plt.show()

f:id:rmizutaa:20190218203414p:plain

性年代による影響が大きく出ています。 10代や20-49歳の男性は比較的ゲームのプレイ時間が長く、女性と50歳以上の男性はゲームのプレイ時間が短くなるという結果になっています。また、TV視聴時間はほぼ影響がなく、子供の有無は少しプラスの影響が出ています。１番見たかったCM視聴の効果は、この推定方法だとゲームのプレイ時間には与える影響は小さいという結果になりました。

最後に、今回のモデルが本当に当てはまりの良いモデルになっていたかを確認します。負の二項分布には２通りのパラメータの置き方があるため、numpyの負の二項分布に対応させるために少し式に変換を加えています。

負の二項分布の式参考
- https://mc-stan.org/docs/2_18/functions-reference/nbalt.html
- https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.negative_binomial.html

#numpyのパラメータに対応させる
mu = la["alpha"].mean()
phi = la["beta"].mean()
p = phi / (mu + phi)

まず、全体のヒストグラムを確認します。

#推定したパラメータからサンプリングしたヒストグラム
plt.figure(figsize=(12,4))
plt.subplot(1, 2, 1)
plt.ylim(0, 10000)
plt.hist(np.random.negative_binomial(phi,p, size=10000),bins=20)
#元のヒストグラム
plt.subplot(1, 2, 2)
plt.ylim(0, 10000)
plt.hist(df["gamesecond"],bins=20)
plt.show()

f:id:rmizutaa:20190218203442p:plain

左が推定したパラメータからランダムサンプリングで作成した分布で、右が元データの分布になります。 0の数はほとんど同じになっています。

１以上の形がわかりにくいので、１以上に制限したヒストグラムも確認します。

#推定したパラメータからサンプリングしたヒストグラム(1以上)
plt.figure(figsize=(12,4))
plt.subplot(1, 2, 1)
plt.ylim(0, 550)
plt.hist(np.random.negative_binomial(phi,p, size=10000),bins=20,range=(1,50))
#元のヒストグラム
plt.subplot(1, 2, 2)
plt.ylim(0, 550)
plt.hist(df["gamesecond"],bins=20,range=(1,50))
plt.show()

f:id:rmizutaa:20190218203456p:plain

１以上で見てもほとんど同じ形となっています。想像以上に当てはまりの良いモデルとなりました。

まとめ

負の二項分布を用いることで、各変数がゲームをプレイする時間に与える効果を測定することができ、ゲームのプレイ時間で考えるとCM視聴の有無はあまり効果がないという結果になりました。

前回ゲームのプレイの有無の部分に対してはCM視聴の効果が出ていたのに今回は出なかった原因としては、ゲームのプレイ時間で考えるとヘビーユーザが大きく平均値に影響を与えてしまうことが考えられます。ユーザ属性への依存度が高いため、多少CMの影響で新規プレイヤーが増えても、あまり結果に影響を与えないのではないかと思います。

今回の推定結果をみるとゲームのプレイ時間という指標はKPIの設定方法としてはあまり適切ではなく、新規者の数や休眠ユーザの復帰率などを目的としておいた方がよかったのではないかと思いました。

また、変数を2,3個入れ替えたりすると結果のプラスマイナスが反転するようなこともよくあるので、実際にやるならドメイン知識のある人が見て結果に違和感がないか確認する作業も必要となりそうです。複数人で同じ課題に対してモデルを作っても違う結果になると思うので、やってみても面白いかなと思います。

使用コードは以下
github.com