何かこのシリーズめちゃくちゃ久しぶりなんですが（汗）、ちょっと最近問題意識を抱いている話題があるのでそれに関連した形でStanでやってみようと思います。
それは時系列の「トレンド」の扱い。ビジネスの現場では、時系列を意識しなくても良い*1クロスセクションデータでは普通に線形*2モデルを組んだりしますが、ことパネルなど時系列データとなると途端にモデルを組んだりせずに「手なり」で適当に近似曲線を引いてしまったり、みたいなことが少なくない印象があります。特に見た目にそれっぽい「トレンド」がある場合は尚更。
ところがどっこい、最近このブログでも取り上げている動的線形モデルのように、時系列データであってもある程度はモデリングできるわけで、ならばいつまでも「見た目」でトレンドを憶測し「手なり」に近似曲線を引く、なんてことはせずにモデリングでズバリ推定するべきだと僕は思うわけです。
とは言え、内容としてはぶっちゃけこのシリーズ記事の2回前のものとほとんど同じなんですが（笑）、ちょっとデータを変えてやってみようと思います。
まずデータを見てみる
新しくGitHubに今回のサンプルデータセットを上げておいたので、落としてきてdata1という名前でインポートして下さい。その上で簡単にプロットしてみましょう。目的変数であるcv（コンバージョン数のつもり）、説明変数であるad1, ad2, ad3（広告投下量のつもり）とがデータフレームとして格納されています。
> plot.ts(data1) > plot(data1$cv,type='l')
見た感じ、目的変数であるcvには何となーくトレンドがありますね。これを例えば「手なり」で近似曲線を描く*3と以下のような感じになるかなと。
これが間違ってるとまでは思いませんが、何よりもこのデータは説明変数ad1, ad2, ad3を伴っていることが分かってます。その影響も加味した上で、トレンドを推定するということをしなければ真のトレンドはおそらく分からないと言って良いかと思います。けれども、「手なり」で近似してもそこは分からないわけです。そこで動的線形モデルを使ってみようというのが今回の記事の主旨です。
Stanで二階差分トレンドモデルを当てはめる
基本的には動的線形トレンドモデルなのでRなら{dlm}で推定できるよなぁと思うところですが、そこは{dlm}のダメなところで何と単回帰はあっても重回帰は実装されていないという。。。そこでStanでぶん回してやりましょう。まず、以下のようなStanコードを用意します。
data {
int<lower=0> N; // 日数
real<lower=0> cv[N]; // CV数
real<lower=0> ad1[N]; // 広告1
real<lower=0> ad2[N]; // 広告2
real<lower=0> ad3[N]; // 広告3
}
parameters {
real<lower=0> a;
real<lower=0> b;
real<lower=0> c;
real d;
real trend[N];
real s_trend;
real s_cv;
}
model {
real q[N];
real cum_trend[N];
for (i in 3:N)
trend[i]~normal(2*trend[i-1]-trend[i-2],s_trend);
cum_trend[1]<-trend[1];
for (i in 2:N)
cum_trend[i]<-cum_trend[i-1]+trend[i];
for (i in 1:N)
q[i]<-cv[i]-cum_trend[i];
for (i in 1:N)
q[i]~normal(a*ad1[i]+b*ad2[i]+c*ad3[i]+d,s_cv);
}
これをtrend_issue.stanという名前で保存した上で、以下のようなRコードでStanをキックし、結果をプロットしてみます。
> library(rstan) # {rstan}のインストールについては過去記事を参照のこと > dat1<-list(N=100,cv=data1$cv,ad1=data1$ad1,ad2=data1$ad2,ad3=data1$ad3) > cv.fit<-stan(file='trend_issue.stan',data=dat1,iter=1000,chains=4) # 中略 > smp<-extract(cv.fit) > t_a<-density(smp$a) > t_b<-density(smp$b) > t_c<-density(smp$c) > t_d<-density(smp$d) > a_w<-t_a$x[t_a$y==max(t_a$y)] > b_w<-t_b$x[t_b$y==max(t_b$y)] > c_w<-t_c$x[t_c$y==max(t_c$y)] > d_w<-t_d$x[t_d$y==max(t_d$y)] > trend_w<-rep(0,100) > for (i in 1:100) { + tmp<-density(smp$trend[,i]) + trend_w[i]<-tmp$x[tmp$y==max(tmp$y)] + } > cv.prd<-a_w*data1$ad1+b_w*data1$ad2+c_w*data1$ad3+d_w+cumsum(trend_w) > matplot(cbind(data1$cv,cv.prd,cumsum(trend_w)+d_w),type='l',lty=1,lwd=c(1,3,3),cex = 1.5,xlab="日数",ylab="CV数") > legend(0,700,text.width = 20,cex = 1.5,y.intersp = 0.4,legend = c("CV実測値","CV推定値","トレンド推定値"),col=c(1,2,3),lty=1,lwd=c(1,3,3)) > cor(data1$cv,cv.prd) [1] 0.9672645 # 推定合致率はかなり良さそう
これでだいぶはっきりしました。緩やかなS字カーブを描く大局的トレンドの上に、広告投下量によって説明されるCV増加分が都度乗っていく。。。その結果としての、やや凸凹した右肩上がりのCV実測値が得られているのだということが分かります。
重要なのは、最初に「手なり」で描いたトレンドとは異なり、動的線形モデルでStanによって推定された二階差分トレンドは最終的には「サチり」を迎えて下降傾向に転じていること。これは目視による「右肩上がり」という推測とは食い違う結果になっています。
理由は色々挙げられますが、やはり大きいのは「説明変数の効果をきちんと加味しているかどうか」の違いでしょう。目視でそこまで正確に見積もるのはやはり難しいわけで*4、そういう意味でもきちんと統計モデリングを行ってトレンドを推定することには意味があるということですね。
おまけ：パラメータの本当の値
以前の記事同様このサンプルデータにも当然本当の値というか正解がありまして、
> c(a_w,b_w,c_w,d_w) [1] 1.7772343 1.7304724 0.4276865 152.6054050 > c(a,b,c,d) [1] 2.0 1.3 0.5 150.0
また正解のトレンドの値をGitHubに置いといたので、落としてきてtrendとかいう名前にしてインポートしてから以下のようにして比較してみましょう。
> cor(trend,trend_w) [1] 0.9570091 > matplot(cbind(cumsum(trend),cumsum(trend_w)),type='l',lty=1,lwd=3)
まぁ大体うまく行ってるかなと。ちなみに、データ自体は以下のようにして生成しました。
> a<-2 > b<-1.3 > c<-0.5 > d<-150 > trend<-rep(0,100) > trend[1]<-0.1 > for (i in 3:100){ + trend[i]<-2*trend[i-1]-trend[i-2]+rnorm(n=1,mean=0,sd=0.05) + } > cv<-a*ad1+b*ad2+c*ad3+d+rnorm(n=100,mean=0,sd=50)+cumsum(trend) # 以下おまけ > t_s_cv<-density(smp$s_cv) > t_s_trend<-density(smp$s_trend) > s_cv_w<-t_s_cv$x[t_s_cv$y==max(t_s_cv$y)] > s_trend<-t_s_trend$x[t_s_trend$y==max(t_s_trend$y)] > s_cv_w [1] 44.3825 # 真の値は50 > s_trend [1] 0.05556993 # 真の値は0.05
実際にはサンプルサイズN=100ということでパラメータ推定精度自体はちょっと良くないんですが、モデル全体の合致度から言えばうまく推定できているのではないかと思います。ついでにMCサンプルの収束具合もチェック。
> library(coda) > cv.fit.coda<-mcmc.list(lapply(1:ncol(cv.fit),function(x) mcmc(as.array(cv.fit)[,x,]))) > plot(cv.fit.coda)
偏回帰係数の方は綺麗に収束してますね（ピークはずれちゃってますが）。ばらつきの標準偏差の方はやや収束が悪い気もしますが。。。ということでお後がよろしいようで。
追記（2016年6月20日）
Stanコードをベクトル表記にした時のやり方を最近になって真面目にやるようになったので、ここに備忘録として追記しておきます。
data {
int<lower=0> N; // 日数、すなわち行数
int<lower=0> M; // 計画行列の列数
real<lower=0> cv[N]; // CV数、目的変数だけベクトルで与える
matrix[N, M] X; // 計画行列として広告1-3を列方向に結合させて与える
}
parameters {
real trend[N];
real s_trend;
real s_cv;
vector<lower=0>[M] beta; // 偏回帰係数をベクトルで取得する
real d;
}
model {
real q[N];
real cum_trend[N];
for (i in 3:N)
trend[i]~normal(2*trend[i-1]-trend[i-2],s_trend);
cum_trend[1]<-trend[1];
for (i in 2:N)
cum_trend[i]<-cum_trend[i-1]+trend[i];
for (i in 1:N)
q[i]<-cv[i]-cum_trend[i];
for (i in 1:N)
q[i]~normal(dot_product(X[i],beta)+d,s_cv);
}
> library(rstan) > dvar <- data1[,-1] > dat1<-list(N=nrow(dvar), M=ncol(dvar), cv=data1$cv, X=dvar) > cv.fit<-stan(file='trend_issue.stan',data=dat1,iter=1000,chains=4) # 中略 > slength <- 2000 # iter=1000, chains=4だと500*4でサンプリング長が2000になる。ここは可変 > fit.smp<-extract(fit) > t_d<-density(fit.smp$d) > d<-t_d$x[t_d$y==max(t_d$y)] > beta<-rep(0,ncol(dvar)) > for (i in 1:ncol(dvar)) { tmp<-density(fit.smp$beta[(slength*(i-1)+1):(slength*i)]) beta[i]<-tmp$x[tmp$y==max(tmp$y)] } > trend<-rep(0,nrow(dvar)) > for (i in 1:nrow(dvar)) { tmp<-density(fit.smp$trend[,i]) trend[i]<-tmp$x[tmp$y==max(tmp$y)] } > beta_prod<-rep(0,ncol(dvar)) > for (i in 1:ncol(dvar)){beta_prod<-beta_prod + dvar[,i]*beta[i]} > pred <- d + beta_prod + cumsum(trend)
さらにseasonality項もtrendと同じやり方で追加できます。
