バナナでもわかる話

計量経済学・統計学が専門の大学院生です。統計学・経済学・投資理論・マーケティング等々に関する勉強・解説ブログ。ときどき趣味も。極力数式は使わずイメージで説明出来るよう心掛けていますが、時々暴走します。

階層ベイズ生存解析を用いたwebサイトの訪問者分析に関するStanでの実装

今回も少し難しめの話。
前回生存解析についてまとめたので、生存解析を使っているマーケティングモデルでも勉強してみようかと探してみたところ、次の論文があったので、論文モデルを整理しながらStanで実装してみようと思います。

前回の生存解析の記事
www.bananarian.net

今回まとめる論文
山口景子 『頻度の時間変化を考慮した階層ベイズモデルによるウェブサイト訪問行動の分析』
https://www.jstage.jst.go.jp/article/marketingscience/22/1/22_220102/_pdf/-char/ja

Stanについてよく知らない人はこの本がオススメ


目次

スポンサーリンク



webマーケティング研究の現状

webマーケティング研究

これは論文の内容ではありませんが、マーケティングサイエンス分野におけるwebマーケティング研究の現状について概観しておこうと思います。

まず、私が調べた限りですが、webマーケティングに関するマーケティング研究の関心は*1最近非常に高まっています。

しかし、その一方でwebマーケティング分析において欠かせないユーザー(or 顧客 or 消費者)のアクセスログデータは基本的に、各企業が管理しており、マーケティング研究者の手元にはありません。その為研究者がこれらの研究顧客分析を行うにあたっては、企業と共同研究を行う必要があり、ハードルが高く、報告も少ないのが現状です。

アクセスログの経営学的意義

経営学においては、消費者の購買行動(購買前、購買後含む)や購買心理を集計的な観点、個別的な観点から分析し、その結果を企業経営に役立つ形でインプリケーションとして提示することが求められます。

しかし購買前行動については、実際に使用できるデータ(消費者が商品購買前に何をしていたかに関するデータ)がほとんどなく、現象として観察することは出来ても計量が難しく、実際に役立てるのが難しいと言われてきました。

ただ今はというと、webサービスやECサイトの登場と共に、膨大なユーザーアクセスログが蓄積され、購買前の行動を十分に観察出来るようになりました。そのため、今経営学分野ではアクセスログを使った研究をしたいというニーズが高まりつつあるわけです。

one to one マーケティング

また、細かなユーザー情報を収拾できるようになったために、ウェブマーケティング分野ではone to oneマーケティングと言って、「消費者一人一人のニーズや状況を汲み取ったマーケティングコミュニケーションを行おう」という状況に移行しつつあります。

ウェブサイトの訪問頻度のモデル化

PV数

ウェブサイトのPV数はそのサイトの広告価値に影響します。そこでPV数を高めたり、今後の潜在的なリスクを見つけるためにどういった要因がPV数に影響を及ぼしているか分析出来たら嬉しいということになります。

さてここで、山口さんはとりあえずPV数を「訪問頻度」と「訪問一回当たりのPV数」に分解し、今回はそのうちの「訪問頻度」について考えてみたというわけですね。

モデル

尤度

 T日間において、ユーザー i=1,\cdots,Nがそれぞれ J_i回サイトに訪れたという状況を考えます。ユーザーiがj回目にサイトに訪れた日数を t_{(i,j)}と置くと、各回の訪問間隔 y_{(i,j)}=t_{(i,j+1)}-t_{(i,j)}はパラメータ \lambda_{ij}の指数分布に従うと仮定します。

ここで、 j=J_iの場合の y_{(i,j)}が本来は観測されないことに注意してください。観察期間 Tとの差を考え、指数分布の切断を考えます。


まず、 j≠J_iの場合の密度関数
 f(y_{(i,j)};\lambda_{ij})=\lambda_{ij}exp(-\lambda_{ij} y_{(i,j)})

 j=J_iの場合は次のように切断を考えます。

 y_{(i,j)}=T-t_{(i,j)}
この時、観測されていない変数 y_{(i,j)}^*=t_{(i,j+1)}-t_{(i,j)}は指数分布に従うため、
 f(y_{(i,j)}^*;\lambda_{ij})=\lambda_{ij}exp(-\lambda_{ij} y_{(i,j)}^*)

ここで、少なくとも y_{(i,j)}^*≧y_{(i,j)}であることはわかっているので、

 Prob(y_{(i,j)}^*≧y_{(i,j)})=1-Prob(y_{(i,j)}^*<y_{(i,j)})=1-F(y_{(i,j)}) = exp(-\lambda_{ij} y_{(i,j)})



更に、顧客の情報を織り込みたいので、次のようにパラメータを規定します。
1を含む顧客情報に関する(K+1)次元ベクトル x_{ij}を用いて、 \lambda_{ij}を次のように考えます。

 \lambda_{ij} = exp(\beta_i x_{ij})

ただし \beta_iは(K+1)次元パラメータベクトル。

事前分布と階層パラメータ

更に事前分布として次のように設定する。

 \beta_i ~ MultiNormal(\delta,\Sigma_{\beta})

 \delta~  MultiNormal(適宜設定)
 \Sigma_{\beta} ~  InvWishart(適宜設定)

共変量の選択

この論文では、目的上経時情報が知りたいので、そのような共変量を選択し、セミパラメトリックなcox比例ハザードモデルを行い、妥当な共変量を判断した後に、再度実行を行っている。

cox比例ハザードモデルについてはこちらをどうぞ
www.bananarian.net

Stanで実装してみた

R で擬似アクセスログを作ってみた - あらびき日記
access_logs.R · GitHub

アクセスログやユーザー情報のデータなんて当然私は持ち合わせていないので、上記サイトさんを参考にアクセスログデータもどきを作らせていただきました。

アクセスログ作成プロセスはこの記事では省略しますので、上記サイトさんを参考にしてください。

model="

functions{
	real TruncatedExp_log(real Y,real lam){
		real pr;
		pr = exp(-lam*Y);
		return(log(pr));
	}
}

data{
	//個人ID数
	int N;
	//特性数
	int CNUM;
	//各IDに対応する訪問回数
	int T[N];
	//T[N]の和
	int SUMT;


	//個人の特性
	matrix[SUMT,CNUM] CX;
	//Truncate情報
	int<lower=0,upper=1> Tru[SUMT];
	//DiffTime=Y
	vector<lower=0>[SUMT] Y;
}

parameters{
	//係数パラメータ
	vector[N] beta0;
	matrix[N,CNUM] beta;
	
	//事前分布のパラメータ
	vector[CNUM+1] BETA;
	cov_matrix[CNUM+1] SIGMA;
}

transformed parameters{
	vector<lower=0>[N] lambda0;
	vector<lower=0>[SUMT] lambda;
	vector[CNUM+1] Ibeta[N];

	for(i in 1:N){
		lambda0[i] = exp(beta0[i]);
		Ibeta[i][1]=beta0[i];
		for(c in 1:CNUM){
			Ibeta[i][c+1]=beta[i,c];
		}
		if(i==1){
			for(t in 1:T[i]){
				lambda[t]=lambda0[i]*exp(dot_product(beta[i],CX[t]));
			}
		}
		else{
			for(t in 1:T[i]){
				lambda[T[i-1]+t]=lambda0[i]*exp(dot_product(beta[i],CX[T[i-1]+t]));
			}
		}
	}
}

model{
	//尤度関数
	for(i in 1:SUMT){
		if(Tru[i]==0){
			Y[i] ~exponential(lambda[i]);
		}
		else{
			Y[i] ~TruncatedExp(lambda[i]);
		}
	}			
	//事前分布
		for(i in 1:N){
			Ibeta[i] ~ multi_normal(BETA,SIGMA);
		}

	//階層事前分布
	//適宜調整(仮に下のような感じにしてみた。下のよりももっと無情報に近づけるのが良い)
	BETA ~ multi_normal(rep_vector(0,CNUM+1),diag_matrix(rep_vector(10,CNUM+1)));
	inverse(SIGMA) ~ wishart(10,diag_matrix(rep_vector(10,CNUM+1)));
}

"

library(rstan)
fit=stan(model_code=model,data=data,seed=1234,par=c("Ibeta"))

雑感

かなり解釈しやすいモデルで、使える場面が多々ありそうだな~と思いました。
論文にもあるように、各人のパラメータを使ってセグメンテーションを行うと直感的にもわかりやすいかもしれないです。

*1:理由は後述