CM視聴の効果を推定する（ゼロ過剰ポアソン編）

Stan python

はじめに

最近は統計モデリングについて学んでいますが、具体的にこういうケースでは統計モデルを使うべきだというケースが自分の中で定まっていません。そのあたりをつけるために実験を行なっていきたいと思います。

データセット

今回は こちらの岩波データサイエンスvol.3のデータを使用しました。

このデータセットは市場調査のデータで、ユーザごとの年齢や性別などのデモグラ情報、CM視聴の有無とCMを実施したゲームのプレイ時間などが入っています。

CM視聴がゲームのプレイ時間に与える効果を推定したいのですが、 CM視聴者にバイアスがかかっているために集計や通常の回帰ではおかしな結果がでます。 そのため上記の本では傾向スコア等の因果推論の枠組みを用いて効果の推定を行なっています。

このデータセットに対する取組はこちらでも詳しく書かれています。 https://tjo.hatenablog.com/entry/2016/08/29/190000

方針

このデータセットにおいて、因果推論を行う必要性がある原因は以下の２点です。
1.CM視聴をしている人としていない人でデモグラ情報の偏りがある
2.目的変数をゲームのプレイ時間としているが、ほとんど(約93%)の人が全くプレイしておらず値に0が入っている

これらは統計モデリングをうまく使えば表現可能ではないかと考えました。1のケースは異なるデモグラ情報ごとに階層ベイズで、2のケースはゼロ過剰ポアソンを使えば対応可能なように見えます。

両方考慮するならゼロ過剰ポアソン+階層ベイズをする必要がありますが、かなり複雑化してしまうことが予測されるため、本記事ではまずゼロ過剰ポアソンのみを用いたCM視聴効果の推定を行います。

ゼロ過剰ポアソンについて

stanとrで統計モデリングの11.3章を参考にさせていただきました。 また。以下のサイトも参考にさせていただきました。 https://mrunadon.github.io/UnderShoesCensoredZIP/

ゼロ過剰ポアソンは、例えば飲食店への来店回数のように、ほとんどの人は0(利用しない)だが、１回以上来店する人はそれなりに来店する、というような場合に使えるモデルになります。

具体的には、まずベルヌーイ分布を用いて0か1以上かを判定し、そこで1と判定されたものはポアソン分布に従うという２段階のモデル構造をしています。

モデル

今回stanのモデルはアヒル本の例題そのままです。 説明変数としては、以下のものを使いました。

cm_dummy:CMの視聴有無
T：10代
F1：20-34歳の女性
F2：35-49歳の女性
F3：50歳以上の女性
M1：20-34歳の男性
M2：35-49歳の男性
M3：50歳以上の男性
TV_watch_day:１日のTV視聴時間(s)
child_dummy：子供の有無

コードは以下のようになります。

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

#スケールを合わせる
df["TVwatch_day"]=df["TVwatch_day"]/10000
df["gamesecond"]=(df["gamesecond"]/100).astype(int)

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

model="""
functions {
  real ZIP_lpmf(int Y, real q, real lambda) {
    if (Y == 0) {
      return log_sum_exp(
        bernoulli_lpmf(0 | q),
        bernoulli_lpmf(1 | q) + poisson_log_lpmf(0 | lambda)
      );
    } else {
      return bernoulli_lpmf(1 | q) + poisson_log_lpmf(Y | lambda);
    }
  }
}

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

parameters {
  vector[D] b[2];
}

transformed parameters {
  vector[N] q_x;
  vector[N] q;
  vector[N] lambda;

  q_x = X*b[1];
  lambda = X*b[2];
  for (n in 1:N)
    q[n] = inv_logit(q_x[n]);
}

model {
  for (n in 1:N)
    Y[n] ~ ZIP(q[n], lambda[n]);
}
"""

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)

うまく収束したので、結果を見ていきます。

まずは１段階目のゲームのプレイ有無への影響をみてみます。

# サンプル列を抽出
la  = fit.extract(permuted=True)
# ゲームプレイ有無への寄与度
plt.figure(figsize=(12,8))
plt.violinplot([la["b"][:,0,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:20190209111607p:plain

性年代でみると１０代がゲームをプレイする確率が高く、女性は年齢が変わっても同じ、男性は若い人の方が確率が高いですね。 CM視聴効果はプラスに出ており、平均値が0.17なので、オッズ比で考えると
np.exp(0.17)/np.exp(0) = 1.185
とCMを視聴していない人に対して約1.2倍ゲームをプレイする確率が上がるそうです。ちなみに子供がいる人はいない人と比べて1.5倍プレイする確率が上がるという結果が出ました。ゲームの詳細情報がないのでわかりませんが、子供と一緒に遊べるようなゲームなんでしょうか？

次に、２段階目のゲームをプレイしている人のプレイ時間への影響を確認していきます。

#ゲームプレイ時間への寄与度
plt.figure(figsize=(12,8))
plt.violinplot([la["b"][:,1,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:20190209111631p:plain

あまりうまく推定できている気がしません。 性年代が違ってもほぼ値が同じですし、CM視聴や子供の有無については効果がほぼないという結果になってしまっています。

ここをうまく推定でないのは心当たりがあって、実はプレイ時間が0ではない人のプレイ時間のヒストグラムは以下のようになっています。

df["gamesecond"].hist(bins=20,range=(1,300))

f:id:rmizutaa:20190209111650p:plain

ポアソン分布ではないのは確かですね。 ガンマ分布？負の二項分布？半コーシー分布？とかだったらうまくいったりするんでしょうか？

まとめ

ゼロ過剰ポアソンを用いてCM視聴がゲームのプレイ時間に与える影響の確認を行いました。

ゲームのプレイの有無に与える効果の部分はうまく推定できた気がしますが、ゲームをプレイする時間に与える効果はポアソン分布ではうまく推定することができませんでした。

今後はポアソン分布の部分を他の分布に変更してやってみるのと、今回できなかった階層モデルについて、stanがInitialization failedを吐かなかったら書きたいと思います。

コードは以下
https://github.com/rmizuta3/CMeffect