廿TT

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

[rstan]ランダム効果入り G-O モデルによるツイートのインプレッションの分析

今日の川柳

Twitter Analytics からツイートのインプレッション(何回表示されたか)などの情報を見ることができます。

ぼくのツイッターアカウントの最近のデータを以下に置いておきます。

tweet_activity_metrics_abiko_ushi_20180116_20180213_en.csv · GitHub

いま、ツイートが何日くらいでどのくらいの人に見られるかを知りたいとします。
@付きのツイートは除外して、横軸にツイートしてからの経過日数、縦軸に累計のインプレッションをとりプロットしてみます。

library(tidyverse)
library(lubridate)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
###
dat <-read_csv("~/Downloads/tweet_activity_metrics_abiko_ushi_20180116_20180213_en.csv")
Tmax <-as.Date("2018-02-12")
dat2 <- dat %>% 
  mutate(date=as.Date(time)) %>% 
  mutate(elapsed=as.integer(Tmax-date)) %>% 
  dplyr::filter(substr(dat$`Tweet text`,0,1)!="@") %>% 
  mutate(id=row_number())

p1<-ggplot(dat2,aes(x=elapsed,y=impressions))+
  geom_point()+
  ylim(c(0,NA))
print(p1)

f:id:abrahamcow:20180212110439p:plain

外れ値があるので、除外してプロット。

p2<-ggplot(dat2,aes(x=elapsed,y=impressions))+
  geom_point()+
  ylim(c(0,500))
print(p2)

f:id:abrahamcow:20180212110514p:plain

日数が経過するほど多くのインプレッションが得られるものの10日くらいからはほぼ横ばいになっています。

こういう動き方の背後にどういうメカニズムがあるのか想像します。

まず、ツイートを見てくれる人には限りがあるので、インプレッションには上限があることが考えられます。

そして、インプレッションの増加量は、インプレッションが上限に達する前の「まだツイートをみていない人」の量に比例すると考えることにします。

まだツイートをみていない人がいっぱいいる内はインプレッションの増加量が大きく、まだツイートをみていない人が少なくなってくると、インプレッションの増加量が小さくなるという仮説です。

この想像を微分方程式で書くと以下のようになります。

 \displaystyle \frac{d}{dt}A(t)=b(m-A(t))

t が経過日数、A(t) が t 日目の累計インプレッション 、m が上限、b は比例定数です。

この方程式はかんたんに解が求まり、初期条件  A(0)=0 の下で解くと、

 \displaystyle A(t)=m(1-\exp(-bt))

となります。

ポアソン過程のレートパラメータが時間の関数  \displaystyle \frac{d}{dt}A(t) であるとしたモデルは、G-O モデルとして知られています。

G-O モデルを使ってツイートのインプレッションを分析したいと思います。

ただし、ツイッターの場合リツイートとかで波及していくことがあるので、ツイートを見てくれる潜在的なユーザー層はツイートごとに異なることが考えられます。

そこで、ツイートごとに異なるランダム効果 r_i を入れて、ツイート i のインプレッションの上限が m\exp(r_i) だとします。

 \displaystyle A(t)=m\exp(r_i)(1-\exp(-bt))

インプレッションがパラメータ A(t)ポアソン分布に従うとし, ランダム効果 r_i が平均 0、分散 \sigma^2正規分布に従うとします。

このモデルを Stan で書くと以下のようになりました。

//hGO_poisson.stan
data{
  int<lower=0> T;
  vector<lower=0>[T] elapsed;
  int<lower=0> PV[T];
  real<lower=0,upper=1> p;
  real<lower=0> sigma;
}
parameters{
  real m;
  real r[T];
  real<lower=0> b;
}
model{
  for(t in 1:T){
   PV[t] ~ poisson_log(m+r[t]+log1m(exp(-b*elapsed[t])));
  }
  r ~ normal(0,sigma);
}
generated quantities{
  int predPV[T];
  real maxIMP[T];
  real ptime;
  for(t in 1:T){
    maxIMP[t] = exp(m+r[t]);
    predPV[t] = poisson_rng(exp(m+r[t])*(1-exp(-b*elapsed[t])));
  }
  ptime = -log1m(p)/b;
}

R から Stan にデータを渡して推定します。

hGO <-stan_model("hGO_poisson.stan")
dat4stan <- list(T=nrow(dat2),elapsed=dat2$elapsed,PV=dat2$impressions,
                 sigma=0.25,p=0.95)
smp <- sampling(hGO,dat4stan,iter=10000,thin=1,seed=1)

トレースプロットと Rhat を見る限り、収束は大丈夫そうです。

tp <- traceplot(smp,pars=c("m","b"))
print(tp)

f:id:abrahamcow:20180212114557p:plain

> all(summary(smp)$summary[,"Rhat"]<1.1)
[1] TRUE

ただし \sigma をかなり小さめに設定しないと、うまく収束しませんでした。

インプレッションの観測値を横軸に、モデルから算出した予測値を縦軸にとりプロットしてみます。両者が一致していれば、散布図の点が一直線にならびます。

predIMP <-get_posterior_mean(smp,pars="predPV")[,"mean-all chains"]
lower=summary(smp,pars="predPV")$summary[,"2.5%"]
upper=summary(smp,pars="predPV")$summary[,"97.5%"]

p_fit <- ggplot(dat2,aes(x=impressions))+
  geom_pointrange(aes(y=predIMP,ymin=lower,ymax=upper))+
  geom_abline(slope=1,intercept = 0,linetype=2)+
  xlim(c(0,NA))+ylim(c(0,NA))

print(p_fit)

f:id:abrahamcow:20180212113101p:plain

縦棒は95%ベイズ予測区間です。

当てはまっていそうです。

各ツイートのインプレッションの上限 m\exp(r_i) と現在のインプレッションを並べてプロットしてみます。

maxIMP=data_frame(current=dat2$impressions,
                  max=get_posterior_mean(smp,pars="maxIMP")[,"mean-all chains"]) %>% 
  mutate(id=row_number()) %>% 
  gather(type,IMP,-id)

p_IMP <-ggplot(maxIMP,aes(x=id,y=IMP,fill=type))+
  geom_col(position = "dodge")+
  theme(legend.title = element_blank())

print(p_IMP)

f:id:abrahamcow:20180212113318p:plain

インプレッションの上限と現在のインプレッションの差がインプレッションの伸びしろの見積もりです。

また、t 日目でインプレッションの上限の何パーセントに達したかは、(A(t)/(m\exp(r_i))) \times 100 でわかります。

A(t)/(m\exp(r_i)) を p と置いて、t について解くことで、ツイートのインプレッションが上限の  p \times 100 パーセントに達するのが何日目かを見積もることができます。

上限の95%に達する日を求めると、だいたい1週間前後になりました。

> summary(smp,pars="ptime")$summary
          mean    se_mean       sd     2.5%     25%      50%      75%    97.5%    n_eff
ptime 7.221732 0.02292474 1.317368 4.895182 6.29811 7.125261 8.062403 10.03966 3302.214
          Rhat
ptime 1.000074

以上です。

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

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