廿TT

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

状態空間モデルで自然検索トラフィックの成長を予測する

場面設定

コンテンツを増やせばそれだけ自然検索にヒットするページが増え、ウェブサイトのトラフィックは増加します。

向こう一年間これだけ記事を書くぞ、というのが決まっていたとして、その計画から自然検索経由の訪問(セッション)数を予測できるでしょうか。

RGA パッケージでデータ抽出

RGA パッケージを使って Google アナリティクスから自然検索経由のセッション数を月次で引っ張ってきます。

library(RGA)
library(dplyr)
authorize()
prof <-list_profiles()

dat1 <-get_ga(profile.id = prof$id[1],
              start.date = "2013-06-01",
              end.date = "2016-01-31",
              dimensions = "ga:yearMonth,ga:landingPagePath",
              segment = "sessions::condition::ga:channelGrouping==Organic Search",
              metrics = "ga:sessions")

どれだけの記事を書いたか、ちゃんと把握してないので、(本当はあんまり良くないのですが)自然検索経由のランディングページの数を書いた記事の数として使うことにします。

monthly <-dplyr::group_by(dat1,year.month) %>>%
  dplyr::summarise(lp = length(landing.page.path),ss=sum(sessions)) %>>%
  dplyr::ungroup()

以上でデータの準備は完了です。

とりあえずプロット

Google アナリティクスの yearMonth を R の日付オブジェクトに変える関数、my_as.Date を定義しておきます。

my_as.Date <- function(year.month){
  tmp <-as.character(year.month)
  tmp <-paste(substr(tmp,1,4),substr(tmp,5,6),"01",sep="/")
  as.Date(tmp)
}

ランディングページの数(lp)とセッション数(ss)をプロットします。

monthly_gg <- gather(monthly,variable,value,-year.month)
monthly_gg$year.month <-my_as.Date(monthly_gg$year.month)

ggplot(monthly_gg)+
  geom_line(aes(x=year.month,y=value,colour=variable),size=1.2)+
  geom_point(aes(x=year.month,y=value,colour=variable),size=2)+
  facet_wrap(~variable,scales = "free",nrow = 2) +
  scale_x_date()

f:id:abrahamcow:20160206072133p:plain

ランディングページは増加を続けており、セッション数もそれにともなって増えているように見えます。

もう少しよくみると1月、7月あたりにぴょこっとセッションが増え、次の月減少、という傾向がありそうです。

このような季節性はあとで予測モデルを作るときに織り込みます。

次に、ランディングページとセッション数の関係を散布図で見てみます。

ggplot(monthly)+
  geom_point(aes(x=lp,y=ss),size=2)

f:id:abrahamcow:20160206072609p:plain

直線的に点がならびました。ランディングページのセッション数へ与える影響は線形になりそうです。

Stan でモデリング

最後の12ヶ月を予測期間、他を学習期間として使うことにします。

ランディングページの数は未来まで計画済みであり、向こう一年間わかっているものとします。

Next <- 12
N <-nrow(monthly) - Next
dat4stan <- list(y=monthly$ss[1:N],x=monthly$lp,N=N,Next=Next)

モデルは以下のようなものを考えました。

y=\alpha_t + \beta x +s_t+\varepsilon_3

y はセッション数、x はランディングページの数です。

\varepsilon_3 は誤差項で、正規分布に従うとします。

\alpha_t は切片ですが時間によって変化します。

その変化のしかたについて、「来月の切片は今月の切片に似ているだろう」という仮定をします。

\alpha_t = \alpha_{t-1} + \varepsilon_1 です。

\varepsilon_1 はノイズで正規分布に従うとします。

この仮定は今日体重が 60kg だったから、明日も 60kg 前後だろう、と予測するようなもので、ちっとも本質的でない、そんなモデルになんの意味があるのか、と思います(ぼくはそういうことを思う)。

しかし、こうすることで、時間的に変化していくもののダイナミクスを捉えることができるのだと考えます。

β は傾きで、これは時間的に変化はしません。

これも変化させたほうがおもしろいかもしれませんが、ランディングページとセッション数の関係は線形になりそうだったので固定しました。

s_t は季節調整項です。

セッション数は季節によって上がったり下がったりするけど、一年通してみればとんとんだ、と考えます。

 \sum_{k=1}^{12} s_{t-k} = 0 + \varepsilon_2 です。

\varepsilon_2 はやはり誤差で正規分布に従うとします。

式変形して s_t は、

 s_t =- \sum_{k=1}^{11} s_{t-k} + \varepsilon_2

と表せます。

季節調整項はあとで推定値をプロットするので、それを見たほうが意味がわかりやすいと思います。

さて、これらを Stan でコーディングしてやります。とくにひねらず愚直にモデル通りに書いていきます。

data{
	int<lower=0> N;
	int<lower=0> Next;
	real x[N+Next];
	real y[N];
}
parameters{
	real alpha[N];
	real beta;
	real s[N];
	real<lower=0,upper=10000> sigma[3];
}
model{
	for(i in 2:N){
		alpha[i] ~ normal(alpha[i-1], sigma[1]);
	}
	for(i in 12:N){
		s[i] ~ normal(-s[i-1]-s[i-2]-s[i-3]-s[i-4]-s[i-5]-s[i-6]-s[i-7]-s[i-8]-s[i-9]-s[i-10]-s[i-11], sigma[2]);
	}
	for(i in 1:N){
		y[i] ~ normal(alpha[i]+beta*x[i]+s[i], sigma[3]);
	}
}
generated quantities{
	real alpha2[N+Next];
	real s2[N+Next];
	real predict[N+Next];
	for(i in 1:N){
		alpha2[i] <- alpha[i];
		s2[i] <- s[i];
	}
	for(i in (N+1):(N+Next)){
		alpha2[i] <- normal_rng(alpha2[i-1],sigma[1]);
		s2[i] <- normal_rng(-s2[i-1]-s2[i-2]-s2[i-3]-s2[i-4]-s2[i-5]-s2[i-6]-s2[i-7]-s2[i-8]-s2[i-9]-s2[i-10]-s2[i-11], sigma[2]);
	}
	for(i in 1:(N+Next)){
		predict[i] <- normal_rng(alpha2[i]+beta*x[i]+s2[i],sigma[3]);
	}
}

R から Stan を操作します。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
timevaryingmodel <- stan_model("~/Documents/dotR/timevaryingmodel.stan")
fit1 <- sampling(timevaryingmodel,data=dat4stan,iter=4000)

結果

観測値と推定値を重ねてプロットしました。

predict_ss <-get_posterior_mean(fit1,"predict")[,"mean-all chains"]

monthly_gg2 <-data.frame(year_month=my_as.Date(monthly$year.month),
                         predict_ss=predict_ss,
                         monthly_ss=monthly$ss)

p1 <-ggplot(monthly_gg2)+
  geom_point(aes(x=year_month,y=monthly_ss))+
  geom_line(aes(x=year_month,y=predict_ss))+
  geom_vline(xintercept = as.integer(as.Date("2015-01-01")),linetype=2)+
  scale_y_continuous(labels=comma)

p1

f:id:abrahamcow:20160206075456p:plain

点が観測値、線が推定値です。

縦の点線の右側が予測期間です。8月、9月のセッション減が予想以上でしたが、ほかは概ね予想に近い結果になっているのではないでしょうか。

95%ベイズ予測区間も合わせてみます。

ex1 <- extract(fit1)
CI <- apply(ex1$predict,2,function(x)quantile(x,c(0.025, 0.975)))
p1 +  geom_ribbon(aes(x=year_month,ymin=CI["2.5%",],ymax=CI["97.5%",]),alpha=0.2)

区間の幅はかなり広いです。しょうがないのかな?

f:id:abrahamcow:20160206075854p:plain

季節調整項の振る舞いを見てみます。

montheffect <- get_posterior_mean(fit1,"s2")[,"mean-all chains"]
year_month =my_as.Date(monthly$year.month)
ggplot()+
  geom_bar(aes(x=year_month,y=montheffect),stat="identity")

f:id:abrahamcow:20160206140239p:plain

多い月で、プラスマイナス 500 くらいの影響があるようです。

最後に、季節調整項を除いた \hat y= \hat \alpha_t + \hat \beta x_t をプロットします。

intercept <- get_posterior_mean(fit1,"alpha2")[,"mean-all chains"]
slope <- get_posterior_mean(fit1,"beta")[,"mean-all chains"]
adjust_ss <-intercept+slope*monthly$lp

ggplot()+
  geom_point(aes(x=year_month,y=adjust_ss))+
  geom_line(aes(x=year_month,y=adjust_ss))+
  scale_y_continuous(labels=comma)+
  geom_vline(xintercept = as.integer(as.Date("2015-01-01")),linetype=2)

f:id:abrahamcow:20160206080614p:plain

これが季節調整済みの正味の自然検索トラフィックの成長具合です。

最初のほうの月はマイナスになってしまいました。

1月、7月にちょっと季節性が残っているような気がしますが、しょうがないのかな?

参考文献

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1