廿TT

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

負の二項分布を用いたリピート回数のモデル(Google アナリティクス)

負の二項分布は「試行が n 回成功するまでに失敗した回数 x の分布」として知られている( 負の二項分布 | r回の成功を得るのに必要な試行回数 )一方で, ポアソン分布のパラメータ λ がガンマ分布するような分布でもある( 可視化で理解する「負の二項分布」 - ほくそ笑む ).

ポアソン分布と捉えると無理がある計数データをモデル化するときに, 負の二項分布は便利そうだ.

ここに以下のグラフのようなデータがある.

f:id:abrahamcow:20161107041041p:plain

これは Google アナリティクスから持ってきたデータで, 今年の10月に当ブログにアクセスしたユーザーを, セッション(訪問)回数ごとに集計したものだ.

Google アナリティクスではユーザーのジェンダーも調べられるようなので, 信用してそのまま使うことにして性別で色分けしてある.

このようなデータを R から取得してグラフにするには以下のようにする.

#後で使うパッケージをまとめてロードしておく
library(RGA)
library(dplyr)
library(cowplot)
library(scales)
authorize()
prof <-list_profiles()

dat1 <-get_ga(profileId = prof$id[1],
              start.date = "2016-10-01",
              end.date = "2016-10-31",
              dimensions = "ga:sessionCount,ga:userGender",
              metrics = "ga:users")

#プロット
ggplot(dat1,aes(x=sessionCount,y=users,fill=userGender))+
  geom_bar(stat = "identity")+
  scale_y_continuous(labels = comma)

パラメタライゼーションと尤度関数

R の dnbinom は負の二項分布の確率関数で, 以下のようにパラメータを置いている.

\displaystyle p(x)=\frac{\Gamma(x+n)}{\Gamma(n) x!} p^n (1-p)^x

引数 mu を指定してやると, 以下のようにパラメータを置いた負の二項分布になり, 「ポアソン分布のパラメータがガンマ分布するような分布」として使うにはこちらのほうが便利である.

\displaystyle p(x)= \frac{\Gamma(x+\alpha)}{\Gamma(\alpha) x!} \left(\frac{\alpha}{\alpha+\mu}\right)^\alpha \left(\frac{\mu}{\alpha+\mu}\right)^x

α は整数でなくてよい.

この式から, 独立同分布に従うサイズ N のサンプル y_i が得られたときの大数尤度関数は,

l(\alpha,\mu) =\sum_{i=1}^{N}\log \Gamma(y_i+\alpha)-N\log\Gamma(\alpha)+\\
~~~~N \alpha \log\left(\frac{\alpha}{\alpha+\mu}\right)+\sum_{i=1}^{N}y_i \log \left(\frac{\mu}{\alpha+\mu}\right)

これは解析的には解けず, なんらかの方法で最大化してやってパラメータを求める.

モデル

今回, セッション回数をリピートの回数と再解釈することで負の二項分布を当てはめようと思う.

セッション回数 1 のユーザーはまだ一回しかサイトに訪れてないのでリピートの回数は 0, セッション回数 2 のユーザーはリピートの回数 1, セッション回数 3 のユーザーはリピートの回数 2,..., 以下同様である.

この「リピートの回数」に負の二項分布を仮定する.

R による実践

Google アナリティクスのデータは集計済み(足し算されたあと)のデータなので, 対数尤度関数を書くと, 以下のようになる.

ll_negbin <- function(par,dat){
  n <- sum(dat$users)
  size <- par[1]
  mu <- par[2]
  sum(dat$users*lgamma(dat$repeats+size))-n*lgamma(size)+
    size*n*log(size/(size+mu))+sum(dat$repeats*dat$users)*log(mu/(size+mu))
}

ここではまず, 男女の区別をせず負の二項分布を当てはめ, 男女別に負の二項分布を当てはめた結果と比較してどちらがいいかを選ぼうと思う.

R には optim というなんでも最大化(最小化)してくれる便利な関数があるのでそれを使う.

パラメータ推定からプロットまでは以下のコードで行う.

dat2 <-dat1 %>% mutate(sessionCount=as.integer(sessionCount)) %>% 
  mutate(repeats=sessionCount-1)

dat_p <-dat2 %>% group_by(repeats) %>%
  summarise(users=sum(users)) %>%
  ungroup()

dat_m <-dat2 %>% 
  filter(userGender=="male")

dat_f <-dat2 %>% 
  filter(userGender=="female")

fit_p <-optim(par=c(1,1),ll_negbin,dat=dat_p,control = list(fnscale=-1))
Theo_p <- sum(dat_p$users)*dnbinom(dat_p$repeats,size=fit_p$par[1],mu=fit_p$par[2])

fit_m <-optim(par=c(1,1),ll_negbin,dat=dat_m,control = list(fnscale=-1))
Theo_m <- sum(dat_m$users)*dnbinom(dat_m$repeats,size=fit_m$par[1],mu=fit_m$par[2])

fit_f <-optim(par=c(1,1),ll_negbin,dat=dat_f,control = list(fnscale=-1))
Theo_f <- sum(dat_f$users)*dnbinom(dat_f$repeats,size=fit_f$par[1],mu=fit_f$par[2])

dat_f$Theo <- Theo_f
dat_m$Theo <- Theo_m

dat_c <-rbind(dat_f,dat_m)

p1<-ggplot(dat_p,aes(x=repeats,y=users))+
  geom_bar(stat = "identity")+
  geom_point(aes(y=Theo_p),size=3,col="salmon")+
  geom_line(aes(y=Theo_p),size=1,col="salmon")
print(p1)
p2<-ggplot(dat_c,aes(x=repeats,y=users))+
  geom_bar(stat = "identity")+
  facet_wrap(~userGender)+
  geom_point(aes(y=Theo),size=2,col="salmon")+
  geom_line(aes(y=Theo),size=1,col="salmon")
print(p2)

実測値を棒グラフ, 負の二項分布による予測値を折れ線グラフで表した.

こちらが男女の区別をせずに当てはめた結果,

f:id:abrahamcow:20161107045525p:plain

こちらが男女別に当てはめた結果である.

f:id:abrahamcow:20161107045558p:plain


棒グラフより, 男女で明らかに傾向が違うので, これは別に扱ったほうがよさそうだぞ, と言ってもいいのであるが, 安心するために尤度比検定を行う.

l0 <-fit_p$value
l1 <-fit_f$value+fit_m$value
pchisq(-2*(l0-l1),2,lower.tail = FALSE) #p-value
#[1] 7.693881e-11

有意水準を 5% とすると, p 値は十分に小さく, 男女間には差があるといえる.

最後に, リピート回数が 0 より大きくなる確率を負の二項分布を使って求めてみる.

どのくらいの人がサイトのリピーターになってくれるのかは, 興味がある.

round(pnbinom(0,size=fit_f$par[1],mu=fit_f$par[2],lower.tail = FALSE),2)
#[1] 0.12
round(pnbinom(0,size=fit_m$par[1],mu=fit_m$par[2],lower.tail = FALSE),2)
#[1] 0.17

女性で 12%, 男性で 17% がリピーターになってくれそうなことがわかった.

モデルの発展としてはもっと他のデモグラフィック属性をいれるとかが考えられる.

abrahamcow.hatenablog.com