読者です 読者をやめる 読者になる 読者になる

廿TT

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

ゼロ過剰負の二項分布によるセッションの間隔のモデル(Google アナリティクス)

Google アナリティクスではセッションの間隔(daysSinceLastSession)という指標を見ることができる.

これはサイトに訪れたユーザーの直前のセッションが何日前だったのかを示すものだ.

f:id:abrahamcow:20170106164331p:plain

ご覧の通り非常にゼロの多いデータで, しかもロングテールなので, これをどう読み解くのか, けっこう攻めあぐねていた.

そこでゼロ過剰負の二項分布(zero-inflated negative binomial distribution)に当てはめることで, 解釈しやすくすることを思いついた.

負の二項分布

離散時間の待ち時間の分布として有名なのが負の二項分布である.

確率関数は次の式で定義される.

\displaystyle p_{\rm NB}(x) = \frac{\Gamma(x+n)}{\Gamma(n) x!} \left(\frac{n}{n+\mu}\right)^n \left(\frac{\mu}{n+\mu}\right)^x

ここで \Gamma(x) はガンマ関数. x は非負の整数である.

R でグラフを書いて, 分布の形を見てみる.

plot(dnbinom(0:100,size=1,mu=20),pch=16,ylim = c(0,0.05),col="royalblue",xlab="",ylab="")
points(dnbinom(0:100,size=2,mu=20),pch=16,col="orange2")
points(dnbinom(0:100,size=0.5,mu=20),pch=16,col="forestgreen")
legend("topright",c("n=1","n=2","n=0.5"),col=c("royalblue","orange2","forestgreen"),
       pch=16)

f:id:abrahamcow:20170106170819p:plain

パラメータ n は分布の散らばり具合を表しており, 1 より大きくなればなるほど分布が平均の周りに集中する.

逆に 1 より小さくなればなるほど, 極端な値がでやすくなる.

ゼロ過剰負の二項分布

セッションの間隔には 0 が多い. このメカニズムを想像してみると, そこには 2 種類のゼロが混じっている可能性が考えられる.

  1. サイトを認識した上でその日のうちに再訪問した.
  2. ユーザー側は再訪問するつもりは特になかったが, セッション切れなどでその日のうちにセッションが 2 回カウントされた.

このような 2 種類のゼロの発生をモデル化したのがゼロ過剰モデル(zero inflated model)である.

ゼロ過剰負の二項分布の確率関数は次の形で与えられる.

\displaystyle p(x) = \begin{cases} \phi +(1-\phi)p_{\rm NB}(0) & x=0\\(1-\phi)p_{\rm NB}(x) & x=1,2,3,\ldots \end{cases}

この式では 2 種類のゼロの発生が混ぜ合わされている.

 \phi は 2 種類目のゼロが発生する確率である.
2 種類目のゼロの発生と 1 種類目のゼロの発生は背反なので足し算, 2 種類目のゼロの発生とゼロ以外の数の発生は独立としてかけ算してある.

R で書くとこう.

dZINB<-function(x,phi,n,mu){
  p=n/(n+mu)
  ifelse(x==0,phi+(1-phi)*dnbinom(x,size=n,mu=mu),(1-phi)*dnbinom(x,size=n,mu=mu))
}

実際に当てはめてみる

データ取得

まず RGA パッケージでデータを引っ張ってくる. 2016年一年間のデータを使うことにした.

Google アナリティクスでは, 新規ユーザーもセッションの間隔が 0 になってしまうので filters = "ga:userType==Returning Visitor" でリピーターのみを抽出している.

# 使うパッケージをまとめてロードする. 
library(RGA)
library(dplyr)
library(tidyr)
library(cowplot)

authorize()
prof <-list_profiles()

dat1 <-get_ga(profileId = prof$id[1],
              start.date = "2016-01-01",
              end.date = "2016-12-31",
              dimensions = "ga:userType,ga:month,ga:daysSinceLastSession",
              filters = "ga:userType==Returning Visitor",
              metrics = "ga:users")

#ディメンションは最初文字列型になっているので, 数値に直す. 
dat2 <- dat1 %>%  
  mutate(daysSinceLastSession=as.integer(daysSinceLastSession))

基礎集計

グラフをいくつか書いてみる.

月ごとのセッションの間隔ごとのユーザー数はこんな感じ.

ggplot(dat2,aes(x=daysSinceLastSession,y=users))+
  geom_bar(stat = "identity",fill="black")+
  facet_wrap(~month)

f:id:abrahamcow:20170106174803p:plain

月ごとのセッションの間隔の平均値(mean), ユーザー数(UU)とゼロの占める割合(zero) を計算してみる.

UU <-dat2 %>% 
  group_by(month) %>% 
  summarise(UU=sum(users),zero=users[daysSinceLastSession==0]/sum(users),
            mean=mean(users*daysSinceLastSession))

UU_g  <- UU %>% 
  gather(key,value,-month)

ggplot(UU_g,aes(x=month,y=value,group=key))+
  geom_line()+
  geom_point()+
  facet_wrap(~key,scales = "free_y",nrow = 4)+
  ylim(c(0,NA))

f:id:abrahamcow:20170106174933p:plain

結果

各月ごとにゼロ過剰負の二項分布を当てはめ, パラメータをプロットしてみる.

#尤度関数
llZINB <- function(par,x,users){
  phi <-boot::inv.logit(par[1])
  mu <- exp(par[2])
  n <- exp(par[3])
  p = n/(n+mu)
  zeropart <- log(phi+exp(log(1-phi)+dnbinom(0,size=n,prob=p,log = TRUE)))*users[x==0]
  nonzeropart <- users[x!=0]*(log(1-phi)+dnbinom(x[x!=0],size=n,prob=p,log = TRUE))
  zeropart+sum(nonzeropart)
}

#尤度関数を optim で最大化
estfunc <- function(data){
  optim(c(phi=1,mu=1,n=1),llZINB,users=data$users,x=data$daysSinceLastSession,control = list(fnscale=-1))
}

#month ごとに推定
fit_all <-by(dat2,dat2$month,estfunc)

#パラメータの取り出し
params <-sapply(fit_all,function(x)c(boot::inv.logit(x$par[1]),exp(x$par[2:3]))) %>% 
  t %>% as.data.frame() %>% 
  mutate(month=c(paste0("0",1:9),10:12))

params_g <-  params %>% 
  gather(parameters,estimates,-month)

#プロット
ggplot(params_g,aes(x=month,y=estimates,group=parameters))+
  geom_line()+
  geom_point()+
  facet_wrap(~parameters,scales = "free_y",nrow = 5)+
  ylim(c(0,NA))

f:id:abrahamcow:20170106180124p:plain

一番注目したいのは \mu (mu) で, これは分布の中心を表す. 20 から 30 くらいで推移している.

ユーザーの訪問間隔が 20日 から 30日を中心に分布しているというのは, それっぽい値である.

この数が小さいと, ユーザーの訪問間隔は短いということで, ユーザーのアクセスが頻繁であることになる.

リターゲティング(リマーケティング)広告やメルマガの配信間隔を変えたときにこれらの数字がどう変化するかを見るのはよさそうである.

\phi (phi) ばかり増えて \mu が増えていなかったら, 日を置いてのアクセスは増えていない, すなわちサイトへの関心は高くなっていないことがわかる.

当てはまりの確認はこのようにするとよいだろう.

まず確率関数で確率を計算する.

つぎに各月のユーザー数の総数で, セッションの間隔ごとのユーザー数を割り, 相対度数を出してやる.

このふたつを比べる.

pred <-left_join(dat2,params,by="month") %>%
  mutate(prob=dZINB(daysSinceLastSession,phi,n,mu)) %>% 
  left_join(UU,by="month")

ggplot(pred)+
  geom_point(aes(x=users/UU,y=prob),stat = "identity")+
  geom_abline(lty=2)+
  facet_wrap(~month)

f:id:abrahamcow:20170106193851p:plain

両者が一致していれば点は直線上にぴたっと並ぶ.

実際, 点はほとんど直線上に位置しており, 当てはまりは良さそうだ.

条件付き確率を使って, 2 番目のゼロ(ユーザー側は再訪問するつもりは特になかったが, セッション切れなどでその日のうちにセッションが 2 回カウントされた)以外の再訪問を計算するには, 普通の負の二項分布に推定されたパラメータを代入してやればよい.

12月のパラメータから, 累積分布関数を計算してみる.

pprob =pnbinom(0:100,size=params[12,"n"],mu=params[12,"mu"])
qprob =qnbinom(0.8,size=params[12,"n"],mu=params[12,"mu"])

ggplot()+
  geom_point(aes(x=0:100,pprob))+
  geom_hline(yintercept = 0.8,linetype=2)+
  geom_vline(xintercept = qprob,linetype=2)+
  xlab("")+ylab("")

f:id:abrahamcow:20170106200018p:plain

8 割のリピーターが38日以内に再訪問している.

これは例えば広告などの間接効果をみるのにどのくらいの期間さかのぼったらいいかとか, そういうのの参考になるかもしれない.

信頼区間

混合分布みたいなものを考えると, デルタ法(optim による glm:最尤推定 + 信頼区間, Wald 検定 - 廿TT でやったような方法)による信頼区間は近似精度が悪いとよく言われる.

そこでパラメトリックブートストラップによる信頼区間と, デルタ法による信頼区間を比べてみる.

パラメトリックブートストラップによる信頼区間とは, 乱数を発生させて実際に推定を何度もやってみて, 推定値がどのくらいばらつくのかを調べるという方法である.

ブートストラップをやるときはだいたい2000回以上回すとよいとされることが多い.

なので2000回回した. 計算時間は12ヶ月分で 5 分くらいだった. 並列処理とかはしていない.

#デルタ法
my_confint <-function(opt){
  my.Std <-sqrt(-diag(solve(opt$hessian)))
  up <- qnorm(0.975,opt$par,my.Std)
  low <- qnorm(0.025,opt$par,my.Std)
  c(lower=low,upper=up)
}

confint_delta <- as.data.frame(t(sapply(fit_all,my_confint)))

params2 <-sapply(fit_all,function(x)x$par) %>% 
  t %>% as.data.frame() %>% 
  mutate(month=1:12)
#####
#パラメトリックブートストラップ
rZINB <- function(length,n,mu,phi){
  flag <-rbinom(length,size = 1,prob = phi)
  ifelse(flag==1,0,rnbinom(length,size=n,mu=mu))
}

run_boot <-function(month){
  parmat <-matrix(,2000,3)
  for(i in 1:2000){
x <-rZINB(UU$UU[month],params$n[month],params$mu[month],params$phi[month])
simdat<-as.data.frame(table(x),stringsAsFactors = FALSE)
simdat$x <- as.numeric(simdat$x)
opt1 <- optim(c(phi=1,mu=1,n=1),llZINB,users=simdat$Freq,x=simdat$x,control = list(fnscale=-1))
parmat[i,] <- opt1$par
  }
  parmat
}

system.time({res_boot <-lapply(1:12,run_boot)})
#ユーザ   システム       経過  
#313.964      7.586    323.569 

lower <-t(sapply(res_boot,function(x)apply(x,2,quantile,0.025)))
upper <-t(sapply(res_boot,function(x)apply(x,2,quantile,0.975)))

ggplot(params2,aes(x=month,y=boot::inv.logit(phi)))+
  geom_line()+
  geom_point()+
  geom_ribbon(aes(ymax=boot::inv.logit(confint_delta[,"upper.phi"]),
                  ymin=boot::inv.logit(confint_delta[,"lower.phi"])),
              fill="blue",alpha=0.3)+
  geom_ribbon(aes(ymax=boot::inv.logit(upper[,1]),
                  ymin=boot::inv.logit(lower[,1])),
              fill="red",alpha=0.3)

ggplot(params2,aes(x=month,y=exp(mu)))+
  geom_line()+
  geom_point()+
  geom_ribbon(aes(ymax=exp(confint_delta[,"upper.mu"]),
                  ymin=exp(confint_delta[,"lower.mu"])),
              fill="blue",alpha=0.3)+
  geom_ribbon(aes(ymax=exp(upper[,2]),
                  ymin=exp(lower[,2])),
              fill="red",alpha=0.3)

ggplot(params2,aes(x=month,y=exp(n)))+
  geom_line()+
  geom_point()+
  geom_ribbon(aes(ymax=exp(confint_delta[,"upper.n"]),
                  ymin=exp(confint_delta[,"lower.n"])),
              fill="blue",alpha=0.3)+
  geom_ribbon(aes(ymax=exp(upper[,3]),
                  ymin=exp(lower[,3])),
              fill="red",alpha=0.3)

下図は青がデルタ法, 赤がパラメトリックブートストラップによる \phi の95%信頼区間である.

f:id:abrahamcow:20170106200606p:plain

こんなにずれるとは思わなかった. 信頼区間を出すときはブートストラップしたほうがいいだろう.

パラメータ n, \mu に関してはそこまで両者に差はない.

f:id:abrahamcow:20170106201134p:plain

f:id:abrahamcow:20170106201202p:plain

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

おもしろい本です. 買ってください.

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

ゼロ過剰ポアソン分布も出てきます.

余談

グラフを「かく」の「かく」は「描く」がただしいらしいけど, 今回からは「書く」と表記することにしました.

「描く」だとなんか大げさな感じがして恥ずかしいのだ.

関連エントリ

abrahamcow.hatenablog.com