廿TT

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

EM アルゴリズムによるゼロ過剰二項分布のパラメータ推定と Google アナリティクスデータの分析

【GTM版】GoogleアナリティクスのClientIDを取得する方法 | SEM Technology を真似して Google アナリティクスでユニークユーザーごとの ID(より正確にはユニークブラウザごとの ID)を取れるようになった.

googleAnalyticsR でデータを取得すると、こんなふうだ.

> head(gadata)
             dimension1 deviceCategory channelGrouping sessions goal2Completions
1 1000012274.1497003540         mobile  Organic Search        1                0
2 1000415978.1471916875         mobile  Organic Search        1                0
3 1000440110.1442338475         mobile  Organic Search        2                0
4 1000533962.1497347896         mobile  Organic Search        1                0
5 1000594455.1498224061         mobile  Organic Search        1                0
6 1000851085.1432133759         mobile  Organic Search        1                0

セッション数が 2 以上のユーザーに関してコンバージョン率(CVR)の分布を見る.

f:id:abrahamcow:20170802012630p:plain

この棒グラフはちょっと特殊なことをやっていて, 一番左端がジャスト 0, それ以降は 0.1 刻みの幅で頻度を集計している.

ジャスト 0 が非常に多い.

ゼロ過剰二項分布(zero-inflated binomial distribution)で表せそう.

ゼロ過剰二項分布

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


\displaystyle p(y) = \begin{cases} \phi +(1-\phi)f(0;p,n) & y=0\\(1-\phi)f(y;p,n) & y=1,2,3,\ldots, n \end{cases}

f(y;p,n) は二項分布の確率関数. \phi は 0 から 1 の 範囲の実数.

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

1 種類目の 0 はふつうの二項分布で自然に発生する 0, 二種類目の 0 はふつうの二項分布とは別のメカニズムで発生する 0.

これを使うとそもそもコンバージョンする気が一切ないクライアントと, そうではないがたまたまコンバージョンしなかったクライアントが区別できそう.

EM アルゴリズム

各クライアントは確率 \phi で「そもそもコンバージョンする気が一切ないクライアント」, 1-\phi で「ふつうの二項分布に従ってコンバージョンするクライアント」になると考える.

潜在変数 z を考え, z=1 ならば確率 1 で 0 が出る, z=0 ならばふつうの二項分布に従う数が出る.

zy は独立とする.

コンバージョン数 y とゼロ過剰二項分布のパラメータ \phi, p が 与えられたときの, z の条件付き期待値は, 条件付き確率の定義から,


\displaystyle E[z|y]=f(z;p,n,\phi) = \frac{f(y,z;p,n, \phi)}{f(z;p,\phi)} \\ 
\displaystyle \quad= \frac{\phi I(y=0)}{\phi+(1-\phi)f(y;p,n)}

と表せる. ここで I は指示関数. (E-step)

独立性の仮定から, 潜在変数 z を含む, 完全な観測が得られたときの対数尤度関数は,


\displaystyle l^{C}(p,\phi) = \sum_{i=1}^{n}(1-z_i)\log p(y_i;p,n_i)+\\ \quad \sum_{i=1}^{n}\left\{z_i\log\phi +(1-z_i)\log(1-\phi) \right\}

となる.

p について微分して 0 と置いて解くと,  (\sum(1-z_i)y_i)/(\sum(1-z_i)n_i), \phi について微分して解くと, \sum z_i/N(N はサンプルサイズ)が得られ, これが EM アルゴリズムの更新式となる. (M-step)

R で書くとこんな感じだ.

doEM <- function(maxit,phi,p,n,Y){
  N = length(Y)
  # E-step
  probs <- function(m,phi,p,n,Y){
    if(m==1){(Y==0) * phi}else{dbinom(Y,n,p) * (1-phi)}
  }
  
  z.update <- function(phi,p,n,Y){
    unnorm <- sapply(1:2,function(i){probs(i,phi,p,n,Y)})
    norms <- apply(unnorm,1,sum)
    unnorm[,1] / norms
  }
  
  # M-step
  phi.update <- function(z,N){sum(z)/N}
  p.update <- function(z,n,Y){sum((1-z)*Y)/sum((1-z)*n)}
  
  for(i in 1:maxit){
    z <- z.update(phi,p,n,Y)
    phi2 <- phi.update(z,N)
    p2 <- p.update(z,n,Y)
    if(abs(phi2-phi)<1e-6 & abs(p2-p)<1e-6){break}
    phi <- phi2
    p <- p2
  }
  list("z"=z, "phi"=phi2, "p"=p2, "iter"=i)
}

データ分析

EM アルゴリズムでパラメータを推定した結果, \phi は 0.804, p は 0.181 だった.

「そもそもコンバージョンする気が一切ないクライアント」が 959, 「そうではないクライアント」が 76 だった.

(これ, クライアント ID ごとにリマーケティングとかできたらおもしろいと思うんだけどできるのかな.)

CVR の分布を出すのはちょっと難しいので, シミュレーションでやったところ当てはまりはこんな感じだった.

f:id:abrahamcow:20170802022702p:plain

棒グラフが実測値, 点がゼロ過剰二項分布による予測値である.

一応縦棒で 95% 区間を示しているけれど, ほとんど見えない.

「そもそもコンバージョンする気が一切ないクライアント」を FALSE, 「そうではないクライアント」を TRUE として, クラスタごとの訪問経路を図示したのが次の帯グラフ,

f:id:abrahamcow:20170802023353p:plain

使用しているデバイスを図示したのが次の帯グラフである.

f:id:abrahamcow:20170802023514p:plain

FALSE はオーガニック検索経由のクライアントが多く, TRUE はソーシャル経由のクライアントが比較的多いことがわかる.

バイスに関しては, FALSE はデスクトップ使用のクライアントが少なく, TRUE は多いことがわかる.

スマホユーザーにとってこのブログはコンバージョンしにくい作りになっているのかもしれない.

R のコード

library(googleAnalyticsR)
library(tidyverse)
#####
ga_auth()
account_list <- ga_account_list()
ga_id <- account_list$viewId[3]
gadata <-
  google_analytics_4(ga_id,
                     date_range = c("2017-07-01","2017-07-31"),
                     metrics = c("sessions","goal2Completions"),
                     dimensions = c("dimension1","deviceCategory","channelGrouping"),
                     max = 50000)
#####
mycut<-function(x){
  ifelse(x==0,"0",as.character(cut(x,breaks =0:10/10)))
}

mylevels<-
  c("0","(0,0.1]","(0.1,0.2]","(0.2,0.3]","(0.3,0.4]","(0.4,0.5]",
    "(0.5,0.6]","(0.6,0.7]","(0.7,0.8]","(0.8,0.9]","(0.9,1]")

gadata2 <-gadata %>% 
  group_by(dimension1) %>% 
  summarise(SS=sum(sessions), CV=sum(goal2Completions)) %>% 
  dplyr::filter(SS>1) %>% 
  mutate(CVR=CV/SS) %>% 
  arrange(desc(CVR)) %>% 
  mutate(cut_CVR=mycut(CVR)) %>% 
  mutate(cut_CVR=factor(cut_CVR,mylevels))

p1 <-ggplot(gadata2,aes(x=cut_CVR))+
  geom_bar(fill="white",colour="black")+
  theme(axis.text = element_text(size=14))+
  xlab("CVR")+ylab("frequency")
print(p1)

###### EM ALGORITHM FUNCTIONS
doEM <- function(maxit,phi,p,n,Y){
  N = length(Y)
  # E-step
  probs <- function(m,phi,p,n,Y){
    if(m==1){(Y==0) * phi}else{dbinom(Y,n,p) * (1-phi)}
  }
  
  z.update <- function(phi,p,n,Y){
    unnorm <- sapply(1:2,function(i){probs(i,phi,p,n,Y)})
    norms <- apply(unnorm,1,sum)
    unnorm[,1] / norms
  }
  
  # M-step
  phi.update <- function(z,N){sum(z)/N}
  p.update <- function(z,n,Y){sum((1-z)*Y)/sum((1-z)*n)}
  
  for(i in 1:maxit){
    z <- z.update(phi,p,n,Y)
    phi2 <- phi.update(z,N)
    p2 <- p.update(z,n,Y)
    if(abs(phi2-phi)<1e-6 & abs(p2-p)<1e-6){break}
    phi <- phi2
    p <- p2
  }
  list("z"=z, "phi"=phi2, "p"=p2, "iter"=i)
}

phi0 <- 0.8
p0 <- 0.5
EM.result <- doEM(1000,phi0,p0,gadata2$SS,gadata2$CV)
print(EM.result$iter)
print(round(EM.result$phi,3))
print(round(EM.result$p,3))
C.hat <- EM.result$z<0.5
print(table(C.hat))
rzinfbinom <- function(N,M,p,phi){
  g <- rbinom(N,1,phi)
  ifelse(g==1,0,rbinom(N,M,p))
}
N=nrow(gadata2)
pred <- matrix(NA,2000,11)
for(i in 1:2000){
  simy1 <-rzinfbinom(N,gadata2$SS,EM.result$p,EM.result$phi)
  cut_pred <- mycut(simy1/gadata2$SS)
  cut_pred <- factor(cut_pred,levels=mylevels)
  tab_pred <- table(cut_pred)
  pred[i,] <- as.numeric(tab_pred)
}

pred_df <- data_frame(cut_CVR=mylevels,
                  pred_lower=apply(pred,2,quantile,0.025),
                  pred_upper=apply(pred,2,quantile,0.975),
                  pred_mean=apply(pred,2,mean))

obs_df <- gadata2 %>%
  group_by(cut_CVR) %>%
  tally()

out_df <- left_join(pred_df, obs_df, by="cut_CVR") %>% 
  mutate(cut_CVR = factor(cut_CVR,levels=mylevels))
  
p2 <-ggplot(out_df,aes(x=cut_CVR,y=n))+
  geom_col(fill="white",colour="black")+
  geom_pointrange(aes(y=pred_mean,ymin=pred_lower,ymax=pred_upper))+
  theme(axis.text = element_text(size=14))+
  xlab("CVR")+ylab("frequency")
print(p2)

gadata3 <- gadata2 %>% 
  mutate(C.hat) %>% 
  select(dimension1,C.hat)
gadata4 <- gadata %>% 
  dplyr::filter(sessions>1) %>% 
  left_join(gadata3,by="dimension1")
ga_device <-gadata4 %>% 
  group_by(deviceCategory,C.hat) %>% 
  summarise(SS=sum(sessions)) %>% 
  ungroup()

p3 <-ggplot(ga_device,aes(x=C.hat,y=SS,fill=deviceCategory))+
  geom_col(position = "fill")+
  scale_y_continuous(label = scales::percent)+
  theme(axis.text = element_text(size=14))+
  xlab("cluster")+ylab("sessions")
print(p3)

ga_channel <-gadata4 %>% 
  group_by(channelGrouping,C.hat) %>% 
  summarise(SS=sum(sessions)) %>% 
  ungroup()

p4 <-ggplot(ga_channel,aes(x=C.hat,y=SS,fill=channelGrouping))+
  geom_col(position = "fill")+
  scale_y_continuous(label = scales::percent)+
  theme(axis.text = element_text(size=14))+
  xlab("cluster")+ylab("sessions")
print(p4)

付録:途中計算

p について

l^{C}(p,\phi)p について微分して 0 と置いて解く.

\displaystyle l^{C}(p,\phi) = \sum_{i=1}^{N} (1-z_i) \left( \log \binom{n_i}{y_i} + y_i \log(p) +(n_i-y_i)\log(1-p)\right) \\ \quad+\sum_{i=1}^{N}( z_i \log(\phi) + (1-z_i) \log(1-\phi)

\displaystyle \frac{\partial}{\partial p}l^{C}(p,\phi) =\sum_{i=1}^{N} (1-z_i) \left( \frac{y_i}{p}  \frac{n_i-y_i}{1-p}\right) =0

\displaystyle \sum_{i=1}^{N} (1-z_i) \left\{ (1-p)y_i  - p(n_i-y_i)\right\} =0

\displaystyle \sum_{i=1}^{N} (1-z_i) \left\{ y_i  - n_i p \right\} =0

\displaystyle \sum_{i=1}^{N} (1-z_i) y_i  - \sum_{i=1}^{N}(1-z_i) n_i p =0

 p=(\sum(1-z_i)y_i)/(\sum(1-z_i)n_i).

φ について

l^{C}(p,\phi)\phi について微分 0 と置いて解く.

\displaystyle \frac{\partial}{\partial \phi}l^{C}(p,\phi) = \sum_{i=1}{N}\left( \frac{z_i}{\phi}- \frac{1-z_i}{1-\phi}\right) =0

\displaystyle \sum_{i=1}^{N}\left( z_i(1-\phi)- (1-z_i) \phi \right) =0

\displaystyle \sum_{i=1}^{N}( z_i-\phi) =0

\phi=\sum z_i/N .