［rstan］混合ロジスティック回帰を用いた検索キーワードのクラスタリング

Search Console のデータを使います。

まずは searchConsoleR でデータを読み込み、プロットします。

対象とするランディングページは
http://abrahamcow.hatenablog.com/entry/2015/01/17/064522
とします。

library(dplyr)
library(cowplot)
library(searchConsoleR)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
#######
#デフォルトでは、90日分のデータをとってくる
scr_auth()
sc_websites <- list_websites()

scdata <- search_analytics(sc_websites[1,1], 
                           dimensions = c("query"),
                           dimensionFilterExp = 
                             "page==http://abrahamcow.hatenablog.com/entry/2015/01/17/064522",
                           rowLimit = 20000)
#######
scdata <- scdata %>% 
  dplyr::filter(!is.na(query))

ggplot(scdata,aes(x=impressions,y=clicks,colour=position))+
  geom_point(size=3)

横軸にインプレッション、縦軸にクリック数をとった検索クエリの散布図は以下です。

当たり前ですがインプレッションが増えるとクリックが増えるという関係が伺えます。

ただし、下図のように傾きの違う線形の関係がふたつ混じっているように見えます。

傾きが大きいほうがクリック率が高くて効率のよいキーワード、傾きの小さいほうが効率の悪いキーワードだと考えられます。

効率のよいキーワードと悪いキーワードで層別して、傾向の違いを見たくなります。

そこで、インプレッション、クリック、掲載順位の数字を標準化して、k-means でクラスタリングしてみます。

標準化には、最小値を引いてレンジで割るという方法を用いました。

normalize <-function(x){(x-min(x))/diff(range(x))}
scdata_norm <- scdata %>% 
  mutate_each(funs(normalize(.)),-query)

cl <- kmeans(select(scdata_norm,clicks,impressions,position),2, iter.max = 200)
ggplot(scdata,aes(x=impressions,y=clicks,colour=factor(cl$cluster)))+
  geom_point(size=3,alpha=0.7)+
  labs(colour="class")

クラスタリングの結果は以下の図をご覧ください。

期待したような、効率のよいキーワードと悪いキーワードというわかれ方はしてくれませんでした。

別の方法を考えます。

モデル

考えたモデルは混合二項分布（rstan で混合二項分布のパラメータ推定 - 廿TT）とほぼ同じです。

クリック数が二項分布に従うと考えます。

試行回数パラメータをインプレッションとします。

成功確率パラメータには説明変数として、掲載順位（position）を入れることにします。

ふたつの二項分布で掲載順位にかかる係数 β は共通とし、切片項 α が異なると仮定します。

このふたつの二項分布が混合比率 φ で混じっているとします。

ここで、各クエリはだいたい半分ずつくらいにわかれてほしいので、φ の事前分布としてパラメータ3、3のベータ分布を仮定します。

パラメータ3、3のベータ分布は次のように中心 0.5 のあたりの密度が大きい分布です。

curve(dbeta(x,3,3),lwd=3)

Stanのコードは以下のようにしました。

data{
  int<lower=0> N;
  int<lower=0> CT[N];
  int<lower=0> IMP[N];
  real<lower=0> pos[N];
  real<lower=0> a;
  real<lower=0> b;
}
parameters{
  ordered[2] alpha;
  real beta;
  real<lower=0,upper=1> phi;
}
transformed parameters{
  matrix[N,2] lp;
  for(i in 1:N){
    lp[i,1] =log(phi)+binomial_logit_lpmf(CT[i]|IMP[i],alpha[1]+beta*pos[i]);
    lp[i,2] =log1m(phi)+binomial_logit_lpmf(CT[i]|IMP[i],alpha[2]+beta*pos[i]);
  }
}
model{
  phi ~ beta(a,b);
  for(i in 1:N){
    target += log_sum_exp(lp[i,]);
  }
}
generated quantities{
  matrix[N,2] label;
  for(i in 1:N){
    label[i,] = softmax(lp[i]')';
  }
}

2000回×4チェインで収束してくれたようです。

クラスタごとのクリック率（ctr）への当てはめプロットは以下のようになりました。

当てはまりはあまりよくなく、改善の余地がありそうですが、今回はこれでいいことにします。

クラスタで色分けして再度プロットしてみると、下図のように期待したわかれ方をしてくれました。

各クエリのクリック数をクラスタごとに棒グラフにしてみます。

効率の悪いクラスタ 1 は以下のようなキーワード（クリック数上位20）、

効率のよいクラスタ 2 は以下のようなキーワード（クリック数上位20）でした。

効率がいいクラスタ 2 は「花合わせ」関連語句を含むクエリが多いことがわかります。

分析対象としたページ
http://abrahamcow.hatenablog.com/entry/2015/01/17/064522
のタイトルは
花合わせの遊び方：よい子とよい大人のための花札入門 - 廿TT
であり、検索者のクエリとページタイトルがマッチしていたほうが、クリックされやすいと解釈できます。

自明といえば自明な結果になりました。

以下に R から Stan を操作して作図するまでのコードを貼ります。

mixture_binomil <-stan_model("mixture_binomial.stan")
datlist <- list(N=nrow(scdata),CT=scdata$clicks,IMP=scdata$impressions,pos=scdata$position,
                a=3,b=3)

smp1 <-sampling(mixture_binomil,datlist)
traceplot(smp1,par=c("alpha","beta","phi"))

s <-matrix(get_posterior_mean(smp1,pars=c("label"))[,"mean-all chains"],ncol=2,byrow=TRUE)
clabel_hat <-apply(s, 1, which.max)

table(clabel_hat)

alphahat <- get_posterior_mean(smp1,pars="alpha")[,"mean-all chains"]
betahat <- get_posterior_mean(smp1,pars="beta")[,"mean-all chains"]

ggplot(scdata[clabel_hat==1,],aes(x=position,y=ctr))+
  geom_point(alpha=0.5)+
  stat_function(fun=function(x){plogis(alphahat[1]+betahat*x)},size=1,colour="royalblue")

ggplot(scdata[clabel_hat==2,],aes(x=position,y=ctr))+
  geom_point(alpha=0.5)+
  stat_function(fun=function(x){plogis(alphahat[2]+betahat*x)},size=1,colour="royalblue")

theme_set(theme_cowplot(font_family = "HiraKakuProN-W3")) 

ggplot(scdata,aes(x=impressions,y=clicks,colour=factor(clabel_hat)))+
  geom_point(size=3,alpha=0.7)+
  labs(colour="class")

ggplot(head(scdata[clabel_hat==1,],20),aes(x=reorder(query,clicks),y=clicks))+
  geom_col()+
  xlab("")+
  coord_flip()

ggplot(head(scdata[clabel_hat==2,],20),aes(x=reorder(query,clicks),y=clicks))+
  geom_col()+
  xlab("")+
  coord_flip()

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