廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

[RStan]差分方程式で呂布カルマのフリースタイルダンジョン出場を振り返る

モチベーション

[dlm]状態空間モデルでトレンドと広告の効果を分離して推定する - 廿TT で、「広告の効果測定において、残存効果、タイムラグをモデルに組み込みたい」というコメントを頂戴したので、それっぽいモデルを考えてみた。

提案モデルはちょっと変わった重回帰です。

はじめに

マスメディアの影響ってすごくて、「呂布カルマ」について Google トレンドをみると明らかにフリースタイルダンジョン出場を境に検索ボリュームが増えていることがわかる。

下図が「呂布カルマ」の Google トレンドです。

f:id:abrahamcow:20170502014124p:plain

R から Google トレンドのデータを読み込むには、gtrendsR というパッケージ(GitHub - PMassicotte/gtrendsR: R functions to perform and display Google Trends queries)が使える。

devtools::install_github("PMassicotte/gtrendsR")

CRAN からインストールしたバージョンだと「検索ボリュームが足りません」みたいなエラーが出て使えなかった。ギットハブ版は動いた。

上の折れ線グラフはこうやって書きました。単に plot(ryofu) とするだけでも、きれいな図がかけます。

library(gtrendsR)
library(cowplot)
ryofu <-gtrends("呂布カルマ",time = "2017-03-01 2017-04-29")

ggplot(ryofu$interest_over_time,aes(x=as.Date(date),y=hits))+
  geom_line()+
  geom_point()+
  theme_cowplot()+
  scale_x_date(date_labels = "%m/%d")+
  theme(legend.position="none")+
  xlab("")

モデル

折れ線グラフをよくみると、フリースタイルダンジョン放送日だけでなく、その前後も検索ボリュームが増えていることがわかる。

前後の日の検索ボリュームの増加も含めてモデル化したい。

直接効果

まず呂布カルマがフリースタイルダンジョンに出た回の放送当日の検索ボリュームの伸びをモデルに組み込もう。

放送当日の検索ボリュームの伸びを「直接効果」と呼ぶことにする。

FSD[t,j] という変数を t が j 回目のフリースタイルダンジョン放送日ならば 1、そうでなければ 0 の値をとるものとする。

FSD[t,j] を要素に持つ行列を FSD とする。

呂布カルマがバトルに出たのは 2017-04-05、2017-04-12、2017-04-19 だけれど、2017-03-29 も「次のチャレンジャーはこいつだ」みたいな感じでちらっとうつっていたとおもうので、2017-03-29、2017-04-05、2017-04-12、2017-04-19 を 1、他の日を 0 とする。

係数 beta = (beta[1], beta[2], beta[3], beta[4])' が直接効果を表すものとし、t 日目の直接効果は FSD[t,] beta である。

予習効果

次に放送前日の検索ボリュームの伸びをモデルに入れたい。

放送前日の検索ボリュームの伸びは「予習効果」と呼ぶことにしよう。

preFSD[t] という変数を、フリースタイルダンジョン放送前日ならば 1、そうでなければ 0 の値をとるものとする。

係数 gamma が予習効果を表すものとし、t 日目の予習効果は gamma preFSD[t] である。

予習効果は、ユーザーが「お、次回は呂布カルマが出るな」と知ることによって生まれると想像できる。

そのため preFSD[t] は次の日に呂布カルマが出るとわかっている日、2017-04-04、2017-04-11、2017-04-18 に 1、他の日は 0 である。

残存効果

さらに、放送後の検索ボリュームの伸びも評価したい。

変数 D[t] を差分方程式、

D[t] -D[t-1] = FSD[t,] beta - q D[t-1]

(0 < q < 1) で定義する。

ただし D[0] = 0 とする。

検索ボリュームの増加量 D[t] -D[t-1] は直接効果 FSD[t,] beta によって増え、一方で定期的に減衰するという関係を表している。

D[t] を移行すると、

D[t] = FSD[t,] beta + (1 - q) D[t-1]

であり、1 - q をあらためて p と置くと、

D[t] = FSD[t,] beta + p D[t-1]

と書ける。

D[t] から直接効果 FSD[t,] beta を引いたものを、「残存効果」と呼ぶことにしよう。

ベースライン

最後に「呂布カルマ」がフリースタイルダンジョンと関係なくもともと持っていた検索ボリュームを、切片項 base で表す。

これを「ベースライン」と呼ぼう。

誤差項

誤差項には正規分布を仮定しました。

Stan のコード

モデルのパラメータ推定には Stan を使いました。

わざわざ Stan 使う必要はあまりないかもしれないけれども、尤度の中で繰り返し計算をしているので、R だと遅いかなあというのと、目的関数が差分方程式のせいでやや複雑になるので、シンプル最尤法だと局所解にはまるのが少し不安というのが Stan 使った理由です。

//delayed_effect.stan
data{
  int<lower=1> T;
  matrix<lower=0,upper=1>[T,4] FSD;
  vector<lower=0,upper=1>[T] preFSD;
  vector[T] hits;
}
parameters{
  real<lower=0> sigma_obs;
  real<lower=0> p;
  vector<lower=0>[4] beta;
  real<lower=0> gamma;
  real<lower=0> base;
}
transformed parameters{
  vector[T] delayed_effect;
  delayed_effect[1] = FSD[1,] * beta;
  for(t in 2:T){
    delayed_effect[t] = FSD[t,]*beta+p*delayed_effect[t-1];
  }
}
model{
  hits ~ normal(base+delayed_effect+gamma*preFSD,sigma_obs);
}
generated quantities{
  vector[T] pred;
  for(t in 1:T){
    pred[t] = normal_rng(base+delayed_effect[t]+gamma*preFSD[t],sigma_obs);
  }
}

無事収束してくれたようです。

f:id:abrahamcow:20170502094658p:plain

結果のプロット

モデルによる予測値(期待値)と95%ベイズ予測区間は下の図です。

f:id:abrahamcow:20170502092840p:plain

まあまあよく当てはまっているんじゃないでしょうか。

直接効果の推定値は下の図です。

f:id:abrahamcow:20170502093102p:plain

最後の R-指定戦のインパクトがもっとも大きいことがわかる。

予習効果は、

f:id:abrahamcow:20170502093223p:plain

残存効果は、

f:id:abrahamcow:20170502093316p:plain

ベースラインは、

f:id:abrahamcow:20170502093333p:plain

です。

検索ボリュームの要因を分解することができました。

R のコード

#devtools::install_github("PMassicotte/gtrendsR")
library(gtrendsR)
library(cowplot)
ryofu <-gtrends("呂布カルマ",time = "2017-03-01 2017-04-29")
tail(ryofu$interest_over_time)
ggplot(ryofu$interest_over_time,aes(x=as.Date(date),y=hits))+
  geom_line()+
  geom_point()+
  theme_cowplot()+
  scale_x_date(date_labels = "%m/%d")+
  theme(legend.position="none")+
  xlab("")
######
#下処理
FSD1 <-as.integer(as.character(ryofu$interest_over_time$date) == "2017-03-29")
FSD2 <-as.integer(as.character(ryofu$interest_over_time$date) == "2017-04-05")
FSD3 <-as.integer(as.character(ryofu$interest_over_time$date) == "2017-04-12")
FSD4 <-as.integer(as.character(ryofu$interest_over_time$date) == "2017-04-19")
FSD=cbind(FSD1,FSD2,FSD3,FSD4)

preFSD <- as.integer(as.character(ryofu$interest_over_time$date) %in% 
  c("2017-04-04","2017-04-11","2017-04-18"))

Tim = nrow(ryofu$interest_over_time)
dat <-list(T=Tim,
           FSD=FSD,
           preFSD=preFSD,
           hits=ryofu$interest_over_time$hits)

#######
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
delayed_effect_model <-stan_model("delayed_effect.stan")

smp1 <-sampling(delayed_effect_model,dat,iter=2000,thin=1,
                control=list(adapt_delta=0.9))

#収束の確認
traceplot(smp1,pars=c("sigma_obs","beta","gamma","p"))
summary(smp1,pars=c("sigma_obs","beta","gamma","p"))

#結果の取り出し
pred <-get_posterior_mean(smp1,"pred")[,"mean-all chains"]
beta <-get_posterior_mean(smp1,"beta")[,"mean-all chains"]
base <-get_posterior_mean(smp1,"base")[,"mean-all chains"]
delayed <-get_posterior_mean(smp1,"delayed_effect")[,"mean-all chains"]
gamma <-get_posterior_mean(smp1,"gamma")[,"mean-all chains"]

lower=summary(smp1,pars="pred")$summary[,"2.5%"]
upper=summary(smp1,pars="pred")$summary[,"97.5%"]

#以下プロット
ggplot(ryofu$interest_over_time,aes(x=date,y=hits))+
  geom_point()+geom_line()+
  geom_line(aes(y=pred,colour="予測値"),size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2)+
  theme(legend.title = element_blank())
  
ggplot(ryofu$interest_over_time,aes(x=date,y=hits))+
  geom_point()+geom_line()+
  geom_line(aes(y=base,colour="ベースライン"),size=1)+
  theme(legend.title = element_blank())

ggplot(ryofu$interest_over_time,aes(x=date,y=hits))+
  geom_point()+geom_line()+
  geom_line(aes(y=FSD %*% beta,colour="直接効果"),size=1)+
  theme(legend.title = element_blank())

ggplot(ryofu$interest_over_time,aes(x=date,y=hits))+
  geom_point()+geom_line()+
  geom_line(aes(y=delayed-FSD %*% beta,colour="残存効果"),size=1)+
  theme(legend.title = element_blank())

ggplot(ryofu$interest_over_time,aes(x=date,y=hits))+
  geom_point()+geom_line()+
  geom_line(aes(y=gamma *preFSD,colour="予習効果"),size=1)+
  theme(legend.title = element_blank())

アマゾンアフィリエイトのコーナー

The Cool Core

The Cool Core

四次元HIP-HOP

四次元HIP-HOP

STRONG

STRONG

関連エントリ

abrahamcow.hatenablog.com