階層ベイズモデルの直帰率分析への適用 with rstan

松浦さんの『StanとRでベイズ統計モデリング』の8章の階層ベイズがすごくわかりやすいなぁと思いつつも、自分の持っているデータで試していなかったので、これを機に実践してみようと思います。
やや変数を追加しているくらいで大した変更点はありませんが、題材としては当ブログのアクセスログにおける直帰率に関するデータで、どのような要素が直帰率に影響を与えるのかを分析します。

目次
・モデル概要
・前処理
・推定
・結果（非階層モデルとの比較）
・参考文献

モデル概要

モデルは8章のロジスティック回帰の階層モデルに一部変数を追加していますが、ほぼそのままです。記事ごとのパラメータやリファラーごとのパラメータを想定しています。
Nは記事数でnはそのインデックス、Cはリファラーの数でcはそのインデックス、Iはログとして残っているセッションの数でiはそのインデックスとなっています。hatebuは記事のはてぶ数、stringlineは記事の行数、holidayは休日or祝日ダミー変数、daytimeは12:00~18:00なら1をとるダミー変数、revisitedは再訪問ユーザーなら1を取るダミー変数となっています。記事ごと・リファラーごとに直帰のしやすさが違う（パラメータが従う正規分布のパラメータがそれぞれ異なる）という仮定のもとに立ったモデルとなります。

$$x[i] = b_{1} + x_{記事}[記事ID[i]] \\ + x_{リファラー}[リファラーID[i]] + x_{セッション}[i] $$

$$q[i] = inverselogit(x[i]) $$

$$Y[i] \sim Bernoulli(q[i]) $$

$$x_{記事}[n] = b_{2}hatebu + b_{3}stringline[n] \\ + b_{記事間の差}[n] $$

$$b_{記事間の差}[n] \sim Normal(0, \sigma_{記事番号}) $$

$$x_{リファラー}[c] = b_{リファラー間の差}[c] $$

$$b_{リファラー間の差}[c] \sim Normal(0, \sigma_{リファラー番号}) $$

$$x_{セッション}[i] = b_{4}holiday[i] 　+ b_{5}divice[i] \\ + b_{6}daytime[i] + b_{7}revisited[i] $$

前処理

GAのAPIからデータを取得して1セッション1記事になるようにデータを作成しています。数ヶ月で25000件ほどデータがあったのですが、計算に時間がかかるので、データ数を2400件くらいにサンプリングしています。

library(RGA)
library(tidyverse)
library(Nippon)

authorize()
prof <-list_profiles()

start_date <- "2017-04-01"
end_date   <- "2017-10-20"

accesslogdata <- get_ga(profileId = prof$id[2],
                        start.date = start_date,
                        end.date = end_date,
                        dimensions = "ga:pagePath,
                        ga:dateHourMinute,
                        ga:deviceCategory,
                        ga:userType,
                        ga:referralPath,
                        ga:fullReferrer",
                        sort = "-ga:sessions",
                        metrics = "ga:sessions,ga:bounces",
                        fetch.by = "day")

#reshape_url
accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))
accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))
accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

accesslogdata_filtered <- accesslogdata %>% filter(sessions == 1)
accesslogdata_filtered <- accesslogdata_filtered %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'t.co/'),'twitter', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'b.hatena.ne.jp/'),'b.hatena.ne.jp', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'d.hatena.ne.jp/'),'d.hatena.ne.jp', fullReferrer)) %>% 
                          mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'facebook.com/'),'facebook.com', fullReferrer))

#refferer list
refferer_cat <- c("google", "(direct)", "twitter", 
                  "yahoo", "b.hatena.ne.jp", "facebook.com",
                  "bing")

daytime_cat <- c("12","13","14","15","16","17","18")

#making weekday data and daytime data
accesslogdata_filtered <- accesslogdata_filtered %>% 
                            filter(fullReferrer %in% refferer_cat) %>% 
                            mutate(date=as.Date(format(substr(accesslogdata_filtered$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d")) %>% 
                            mutate(hourminutes = substr(accesslogdata_filtered$dateHourMinute,start = 9,stop = 10 )) %>% 
                            mutate(holiday = ifelse(is.jholiday(date), 1, 0)) %>% 
                            mutate(daytime = ifelse(hourminutes %in% daytime_cat, 1, 0)) %>% 
                            mutate(device = ifelse(deviceCategory == "desktop", 1, 0)) %>% 
                            mutate(revisited = ifelse(userType == "Returning Visitor", 1, 0))

selected_dataset <- accesslogdata_filtered %>% 
                      select(pagePath,fullReferrer,
                             device,revisited,holiday,daytime,bounces) %>% 
                      filter(!(pagePath %in% c("article_10", "article_53") )) %>% 
                      bind_cols(rand =runif(nrow(selected_dataset), min = 0, max = 1)) %>% 
                      filter(rand <= 2400/nrow(selected_dataset))

referer_cat <- selected_dataset %>% 
                  select(fullReferrer) %>% 
                  distinct() %>% 
                  mutate(referer_id = 1:n())

article_cat <- selected_dataset %>% 
                  select(pagePath) %>% 
                  distinct() %>% 
                  mutate(article_id = 1:n())

selected_dataset <- selected_dataset %>% 
                      left_join(article_cat, by="pagePath") %>% 
                      left_join(referer_cat, by="fullReferrer") %>% 
                      select(-pagePath, -fullReferrer)

article_data <- read_csv(file = "kamonohashiperry_text.csv")
article_data$url <- gsub(article_data$url,
                         pattern = "http://kamonohashiperry.com/archives/",replacement = "article_")

article_data <- article_data %>% left_join(article_cat, by=c("url"="pagePath"))
article_data <- article_data[!(is.na(article_data$article_id)),]

library(RGA)

library(tidyverse)

library(Nippon)

authorize()

prof <-list_profiles()

start_date <- "2017-04-01"

end_date <- "2017-10-20"

accesslogdata <- get_ga(profileId = prof$id[2],

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:referralPath,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:bounces",

fetch.by = "day")

#reshape_url

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))

accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

accesslogdata_filtered <- accesslogdata %>% filter(sessions == 1)

accesslogdata_filtered <- accesslogdata_filtered %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'t.co/'),'twitter', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'b.hatena.ne.jp/'),'b.hatena.ne.jp', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'d.hatena.ne.jp/'),'d.hatena.ne.jp', fullReferrer)) %>%

mutate(fullReferrer = ifelse(grepl(x = fullReferrer,'facebook.com/'),'facebook.com', fullReferrer))

#refferer list

refferer_cat <- c("google", "(direct)", "twitter",

"yahoo", "b.hatena.ne.jp", "facebook.com",

"bing")

daytime_cat <- c("12","13","14","15","16","17","18")

#making weekday data and daytime data

accesslogdata_filtered <- accesslogdata_filtered %>%

filter(fullReferrer %in% refferer_cat) %>%

mutate(date=as.Date(format(substr(accesslogdata_filtered$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d")) %>%

mutate(hourminutes = substr(accesslogdata_filtered$dateHourMinute,start = 9,stop = 10 )) %>%

mutate(holiday = ifelse(is.jholiday(date), 1, 0)) %>%

mutate(daytime = ifelse(hourminutes %in% daytime_cat, 1, 0)) %>%

mutate(device = ifelse(deviceCategory == "desktop", 1, 0)) %>%

mutate(revisited = ifelse(userType == "Returning Visitor", 1, 0))

selected_dataset <- accesslogdata_filtered %>%

select(pagePath,fullReferrer,

device,revisited,holiday,daytime,bounces) %>%

filter(!(pagePath %in% c("article_10", "article_53") )) %>%

bind_cols(rand =runif(nrow(selected_dataset), min = 0, max = 1)) %>%

filter(rand <= 2400/nrow(selected_dataset))

referer_cat <- selected_dataset %>%

select(fullReferrer) %>%

distinct() %>%

mutate(referer_id = 1:n())

article_cat <- selected_dataset %>%

select(pagePath) %>%

distinct() %>%

mutate(article_id = 1:n())

selected_dataset <- selected_dataset %>%

left_join(article_cat, by="pagePath") %>%

left_join(referer_cat, by="fullReferrer") %>%

select(-pagePath, -fullReferrer)

article_data <- read_csv(file = "kamonohashiperry_text.csv")

article_data$url <- gsub(article_data$url,

pattern = "http://kamonohashiperry.com/archives/",replacement = "article_")

article_data <- article_data %>% left_join(article_cat, by=c("url"="pagePath"))

article_data <- article_data[!(is.na(article_data$article_id)),]

推定

stanコードはこちらになります。

data {
  int N; #the number of article
  int C; #the number of referer
  int I; #the number of log
  int<lower=0> hatebu[N]; 
  real<lower=0> stringline[N];
  int<lower=1, upper=N> article_id[I];
  int<lower=1, upper=C> referer_id[I];
  real<lower=0, upper=1> holiday[I];
  real<lower=0, upper=1> device[I];
  real<lower=0, upper=1> daytime[I];
  real<lower=0, upper=1> revisited[I];
  int<lower=0, upper=1> Y[I];
}

parameters {
  real b[7];
  real b_P[N];
  real b_C[C];
  real<lower=0> s_P;
  real<lower=0> s_C;
}

transformed parameters {
  real x_P[N];
  real x_C[C];
  real x_J[I];
  real x[I];
  real q[I];
  for (n in 1:N)
    x_P[n] = b[2]*stringline[n] + b[3]*hatebu[n] + b_P[n];
  for (c in 1:C)
    x_C[c] = b_C[c];
  for (i in 1:I) {
    x_J[i] = b[4]*holiday[i] + b[5]*device[i] + b[6]*daytime[i] + b[7]*revisited[i];
    x[i] = b[1] + x_P[article_id[i]] + x_C[referer_id[i]] + x_J[i];
    q[i] = inv_logit(x[i]);
  }
}

model {
  for (n in 1:N)
    b_P[n] ~ normal(0, s_P);
  for (c in 1:C)
    b_C[c] ~ normal(0, s_C);
  for (i in 1:I)
    Y[i] ~ bernoulli(q[i]);
}

data {

int N; #the number of article

int C; #the number of referer

int I; #the number of log

int<lower=0> hatebu[N];

real<lower=0> stringline[N];

int<lower=1, upper=N> article_id[I];

int<lower=1, upper=C> referer_id[I];

real<lower=0, upper=1> holiday[I];

real<lower=0, upper=1> device[I];

real<lower=0, upper=1> daytime[I];

real<lower=0, upper=1> revisited[I];

int<lower=0, upper=1> Y[I];

}

parameters {

real b[7];

real b_P[N];

real b_C[C];

real<lower=0> s_P;

real<lower=0> s_C;

}

transformed parameters {

real x_P[N];

real x_C[C];

real x_J[I];

real x[I];

real q[I];

for (n in 1:N)

x_P[n] = b[2]*stringline[n] + b[3]*hatebu[n] + b_P[n];

for (c in 1:C)

x_C[c] = b_C[c];

for (i in 1:I) {

x_J[i] = b[4]*holiday[i] + b[5]*device[i] + b[6]*daytime[i] + b[7]*revisited[i];

x[i] = b[1] + x_P[article_id[i]] + x_C[referer_id[i]] + x_J[i];

q[i] = inv_logit(x[i]);

}

model {

for (n in 1:N)

b_P[n] ~ normal(0, s_P);

for (c in 1:C)

b_C[c] ~ normal(0, s_C);

for (i in 1:I)

Y[i] ~ bernoulli(q[i]);

}

rstanを用いたstan実行用のRコードです。ヴァイオリンプロットで主要な係数の分布を見る処理も書かれています。

library(rstan)

N <- nrow(article_cat)
C <- nrow(referer_cat)
I <- nrow(selected_dataset)

data <- list(N = N,
             C = C,
             hatebu = article_data$hatebu,
             stringline = article_data$str_lines,
             article_id = selected_dataset$article_id,
             referer_id = selected_dataset$referer_id,
             holiday = selected_dataset$holiday,
             device = selected_dataset$device,
             daytime = selected_dataset$daytime,
             revisited = selected_dataset$revisited,
             Y = selected_dataset$bounces)

fit <- stan(file = "model/model8-8_access_analysis.stan",
            data = data,
            pars = c("b", "b_P", "b_C", "s_P", "s_C", "q"),
            seed = 1234)

source('../common.R')

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', paste0('b', 1:7), 's_P', 's_C')
d_est <- data.frame(1:N_mcmc, ms$b, ms$s_P, ms$s_C)
colnames(d_est) <- param_names
d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])
d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')
d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()
p <- p + theme_bw(base_size=18)
p <- p + coord_flip()
p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')
p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)
p <- p + labs(x='parameter', y='value')
p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))
p

library(rstan)

N <- nrow(article_cat)

C <- nrow(referer_cat)

I <- nrow(selected_dataset)

data <- list(N = N,

C = C,

hatebu = article_data$hatebu,

stringline = article_data$str_lines,

article_id = selected_dataset$article_id,

referer_id = selected_dataset$referer_id,

holiday = selected_dataset$holiday,

device = selected_dataset$device,

daytime = selected_dataset$daytime,

revisited = selected_dataset$revisited,

Y = selected_dataset$bounces)

fit <- stan(file = "model/model8-8_access_analysis.stan",

data = data,

pars = c("b", "b_P", "b_C", "s_P", "s_C", "q"),

seed = 1234)

source('../common.R')

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

param_names <- c('mcmc', paste0('b', 1:7), 's_P', 's_C')

d_est <- data.frame(1:N_mcmc, ms$b, ms$s_P, ms$s_C)

colnames(d_est) <- param_names

d_qua <- data.frame.quantile.mcmc(x=param_names[-1], y_mcmc=d_est[,-1])

d_melt <- reshape2::melt(d_est, id=c('mcmc'), variable.name='X')

d_melt$X <- factor(d_melt$X, levels=rev(levels(d_melt$X)))

p <- ggplot()

p <- p + theme_bw(base_size=18)

p <- p + coord_flip()

p <- p + geom_violin(data=d_melt, aes(x=X, y=value), fill='white', color='grey80', size=2, alpha=0.3, scale='width')

p <- p + geom_pointrange(data=d_qua, aes(x=X, y=p50, ymin=p2.5, ymax=p97.5), size=1)

p <- p + labs(x='parameter', y='value')

p <- p + scale_y_continuous(breaks=seq(from=-2, to=6, by=2))

結果

係数を見る限りは、符号の向きが確かなのはb5(PCダミー)とb7(再訪問ユーザーダミー)なので、PCの方が直帰しにくく、再訪問ユーザーの方が直帰しにくいという傾向があると考えることができます。

教科書ではAUCを非階層モデルと比較していましたので、比較してみようと思います。
AUCの計算を行うためのコードもGithubに載っていましたのでそちらを使います。

#calculate auc
library(pROC)

ms <- rstan::extract(fit)
N_mcmc <- length(ms$lp__)
spec <- seq(from=0, to=1, len=201)
probs <- c(0.1, 0.5, 0.9)

auces <- numeric(N_mcmc)
m_roc <- matrix(nrow=length(spec), ncol=N_mcmc)
for (i in 1:N_mcmc) {
  roc_res <- roc(selected_dataset$bounces, ms$q[i,])
  auces[i] <- as.numeric(roc_res$auc)
  m_roc[,i] <- coords(roc_res, x=spec, input='specificity', ret='sensitivity')
}
quantile(auces, prob=probs)

#calculate auc

library(pROC)

ms <- rstan::extract(fit)

N_mcmc <- length(ms$lp__)

spec <- seq(from=0, to=1, len=201)

probs <- c(0.1, 0.5, 0.9)

auces <- numeric(N_mcmc)

m_roc <- matrix(nrow=length(spec), ncol=N_mcmc)

for (i in 1:N_mcmc) {

roc_res <- roc(selected_dataset$bounces, ms$q[i,])

auces[i] <- as.numeric(roc_res$auc)

m_roc[,i] <- coords(roc_res, x=spec, input='specificity', ret='sensitivity')

}

quantile(auces, prob=probs)

＃ロジスティック回帰の階層ベイズ推定でのAUC
> quantile(auces, prob=probs)
      10%       50%       90% 
0.6683084 0.6770603 0.6846594


#ロジスティック回帰のベイズ推定でのAUC
> quantile(auces_non_hiral, prob=probs_non_hiral)
      10%       50%       90% 
0.5463287 0.5623697 0.5776898

＃ロジスティック回帰の階層ベイズ推定でのAUC

> quantile(auces, prob=probs)

10% 50% 90%

0.6683084 0.6770603 0.6846594

#ロジスティック回帰のベイズ推定でのAUC

> quantile(auces_non_hiral, prob=probs_non_hiral)

10% 50% 90%

0.5463287 0.5623697 0.5776898

80%が良いとされているAUCには程遠いですが、記事やリファラーごとの差を考慮しない非階層のものよりもAUCが高いと言えます。
ちなみに、教科書の例のAUCは80%ほどでした。
Webマーケのデータ分析においてロジットは汎用性が高いで、今回のコードを土台に色々と業務で試していこうと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル原著第2版

確率的プログラミングライブラリ「Edward」まとめ

Edwardで何ができるのか知らなかったので、忘備録として残しておきます。

目次
・Edwardとは
・Edwardでできること
・参考スライド
・参考文献

Edwardとは

・LDAで有名なコロンビア大のBlei先生の研究室で、2016年より開発されている確率的プログラミング((　プログラミング言語の変数をモデルの構成要素として使うプログラミング))のPythonライブラリ。
・積み木のように明快な形で確率的モデリングを行うことができる。（モデル→推論→評価　を一括でできる。）
・ベイズ統計と機械学習、深層学習、確率的プログラミングを融合させている。
・計算の際にTensorFlowを用いている。TensorBoardを可視化の際に用いることもできる。
・計算速度がStanやPyMC3よりも速い。GPUを用いた高速化も可能。((　pip install tensorflow-gpuでGPU版のTensolFlowを入れておく必要がある。))
・統計学者のGeorge Edward Pelham Boxから名前を取っている。

Edwardでできること

一般的なベイズ推定は当然ながら、深層学習向けのベイズ適用系の事例が豊富なようです。

・ベイズ線形回帰( Supervised Learning (Regression) )
・バッチトレーニング（巨大なデータセットにおける学習で用いる）( Batch Training )
・Tensorboardを用いた可視化（ Tensorboard ）
・Automated Transformations( Automated Transformations )
・線形混合効果モデル( Linear Mixed Effects Models )
・教師あり学習による分類( Supervised Learning (Classification) )
・教師なし学習( Unsupervised Learning )
・ニューラルネットワークの潜在空間モデル( Latent Space Models for Neural Data )
・混合密度ネットワーク( Mixture Density Networks )
・GAN( Generative Adversarial Networks )
・確率的デコーダー( Probabilistic Decoder )
・ネットワークの推論( Inference Networks )
・ベイジアンニューラルネットワーク( Bayesian Neural Network )
・確率的PCA(主成分分析)( Probabilistic PCA )

jupyterのコードたちはblei-lab/edwardのnotebookに載っています。

2層のニューラルネットワークへのベイズ推定の適用(
Bayesian Deep Learning with Edward (and a trick using Dropout) – Andrew Rowan – PyData London 2017)
Gounosyの方のブログによると、CTR予測などで扱うことができるようです。

参考スライド

確率的プログラミングライブラリEdward from Yuta Kashino

ベイズ統計の数理と深層学習 @db analytics showcase Sapporo 2017 from Takeshi Yoneda

（DL hacks輪読）Bayesian Neural Network from Masahiro Suzuki

参考文献

Edward: A library for probabilistic modeling, inference, and criticism

DEEP PROBABILISTIC PROGRAMMING

EdwardでBayesian DNN＋Variational Inferenceをやってみた話

DEEP PROBABILISTIC PROGRAMMING —”深層学習＋ベイズ”のライブラリ— Edwardの紹介

【Edward】MCMCの数学的基礎からStochastic Gradient Langevin Dynamicsの実装まで

Hello, world! Stan, PyMC3, and Edward
stanの開発者の方がstanとPyMC3とEdwardを比較しています。

Pythonで体験するベイズ推論 PyMCによるMCMC入門

RのFactoRizationMachinesパッケージを用いたFMのページビューデータへの適用

Googleタグマネージャーで集めたアクセスログデータを用いて、前回と同様に記事のレコメンドにチャレンジしてみようと思います。FactoRizationMachinesパッケージという便利そうなパッケージの存在も知れたことから、今回は以前から気になっていたFactorization Machineを扱います。

【目次】
・Factorization Machine(FM)の概要
・パッケージ紹介とインストール
・サンプルデータの構造把握と前処理
・FMの実行
・結果
・参考文献

Factorization Machine(FM)の概要

組み合わせ特徴量を扱う教師あり学習モデル。行列分解とSVMを合体させた手法。
スパースになりやすいデータの予測問題で扱う。
1ユーザーのある商品に対しての評価を、1評価1行として表して、ユーザーとアイテムの交互作用の特徴ベクトルを扱う。
相互作用項に関して、時間や文脈などを自由に入れられる。
相互作用項を次元圧縮する際の要素数を事前に決める必要がある。
Matrix Factorizationよりも精度が良いとされている。特徴量エンジニアリングなどで使われているようです。（Click-Through Rate Prediction）

パッケージ紹介とインストール

FactoRizationMachinesパッケージは線形SVMと2次のFMと高次のFMを実行することができ、引数で正則化項も加えることができます。現段階においては回帰のみで分類問題への適用は今後の開発となるようです。CRANから普通に
install.packages(‘FactoRizationMachines’)
でインストールします。libFMexeパッケージの場合は、libFMをインストールしてパスを指定しておく必要がありますが、このパッケージに関しては不要となります。

サンプルデータの構造把握と前処理

FactoRizationMachinesパッケージのサンプルコードにおいては、MovieLensのデータがサンプルデータとして載せられていました。ユーザーのID(整数)、映画のID(整数)、評価(整数、5段階)、日時(整数)からなるデータに対して、sparseMatrixに変換していました。

今回は、前回の投稿非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみると同じデータを使って、アクセスログデータに適用しようと思います。FactoRizationMachinesの形式に合わせるために、このブログのアクセスログも、クッキーのIDを整数に、記事のIDを整数に、閲覧回数を5段階(5以上を5に変換)に、日時を整数に変更しています。

library(RGA)
library(tidyverse)
library(reshape2)

authorize()
prof <-list_profiles()

start_date <- "2017-07-04"
end_date   <- "2017-10-07"

#demention until 7
accesslogdata <- get_ga(profileId = prof$id[2],
                        start.date = start_date,
                        end.date = end_date,
                        dimensions = "ga:pagePath,
                        ga:dimension1,
                        ga:dateHourMinute,
                        ga:deviceCategory,
                        ga:userType,
                        ga:fullReferrer",
                        sort = "-ga:sessions",
                        metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",
                        fetch.by = "day")

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))
accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))
accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

#ユーザーIDの整数化
user <- unique(accesslogdata$dimension1)
user <- data.frame(dimension1=user) %>% mutate(user_number=1:n())

#記事IDの整数化
articles <- unique(accesslogdata$pagePath)
articles <- data.frame(pagePath=articles) %>% mutate(article_number=1:n())

#データの結合
accesslogdata <- accesslogdata %>% mutate(date=as.numeric(as.POSIXlt(as.Date(format(substr(accesslogdata$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d"))))
accesslogdata <- accesslogdata %>% left_join(articles,by="pagePath")
accesslogdata <- accesslogdata %>% left_join(user,by="dimension1")

#最終接触日の抽出
last_date <- accesslogdata %>% select(user_number,article_number,date) %>% 
                arrange(user_number,desc(date))
last_date <- last_date[!duplicated(last_date[c("user_number","article_number")]),]

#各ユーザーの各記事に対するページビュー数の集計
accesslogdata_sum <- accesslogdata %>% group_by(article_number,user_number) %>% summarise(pageview=n()) %>% 
                        select(user_number,article_number,pageview)
accesslogdata_sum <- accesslogdata_sum %>% mutate(pageview=ifelse(pageview < 5,pageview,5))

accesslogdata_sum <- accesslogdata_sum %>% left_join(last_date, by = c("user_number" = "user_number", "article_number" = "article_number"))

library(RGA)

library(tidyverse)

library(reshape2)

authorize()

prof <-list_profiles()

start_date <- "2017-07-04"

end_date <- "2017-10-07"

#demention until 7

accesslogdata <- get_ga(profileId = prof$id[2],

start.date = start_date,

end.date = end_date,

dimensions = "ga:pagePath,

ga:dimension1,

ga:dateHourMinute,

ga:deviceCategory,

ga:userType,

ga:fullReferrer",

sort = "-ga:sessions",

metrics = "ga:sessions,ga:goal1Completions,ga:pageviews",

fetch.by = "day")

accesslogdata$pagePath <- vapply(strsplit(accesslogdata$pagePath,"\\?"), `[`, 1, FUN.VALUE=character(1))

accesslogdata <- accesslogdata %>% filter(grepl(x = pagePath,"/archives/[0-9]+$"))

accesslogdata$pagePath <- gsub(accesslogdata$pagePath,pattern = "/archives/",replacement = "article_")

#ユーザーIDの整数化

user <- unique(accesslogdata$dimension1)

user <- data.frame(dimension1=user) %>% mutate(user_number=1:n())

#記事IDの整数化

articles <- unique(accesslogdata$pagePath)

articles <- data.frame(pagePath=articles) %>% mutate(article_number=1:n())

#データの結合

accesslogdata <- accesslogdata %>% mutate(date=as.numeric(as.POSIXlt(as.Date(format(substr(accesslogdata$dateHourMinute,start = 1,stop = 12), format="%Y%m%d%"),format = "%Y%m%d"))))

accesslogdata <- accesslogdata %>% left_join(articles,by="pagePath")

accesslogdata <- accesslogdata %>% left_join(user,by="dimension1")

#最終接触日の抽出

last_date <- accesslogdata %>% select(user_number,article_number,date) %>%

arrange(user_number,desc(date))

last_date <- last_date[!duplicated(last_date[c("user_number","article_number")]),]

#各ユーザーの各記事に対するページビュー数の集計

accesslogdata_sum <- accesslogdata %>% group_by(article_number,user_number) %>% summarise(pageview=n()) %>%

select(user_number,article_number,pageview)

accesslogdata_sum <- accesslogdata_sum %>% mutate(pageview=ifelse(pageview < 5,pageview,5))

accesslogdata_sum <- accesslogdata_sum %>% left_join(last_date, by = c("user_number" = "user_number", "article_number" = "article_number"))

FMの実行

デフォルトの設定c(1, 10)では線形のウェイトが有効で、2次の項の要素数が10で正則化項なしのFMを実行することになります。引数に関する詳しい情報はPackage ‘FactoRizationMachines’に書かれています。今回はサンプルを参考に正則化項ありでモデルを実行します。まず、アクセスログデータに対して、ユーザーのIDからなる整数ベクトル、記事のIDからなる整数ベクトル、セッションのあった日時のデータからなる整数ベクトルを作成し、sparseMatrix関数を用いて元データを変形し、80%のデータをトレーニングに、20%のデータをテストに割り当てます。さらに、テストデータに関して、予測値との平均二乗誤差を計算します。

library(FactoRizationMachines)
library(Matrix)

user = accesslogdata_sum$user_number
items = accesslogdata_sum$article_number + max(user)
wdays = (as.POSIXlt(accesslogdata_sum$date,origin="2017-07-04")$wday+1)+max(items)

# Transform access log to feature form
data = sparseMatrix(i=rep(1:nrow(accesslogdata_sum),3),j=c(user,items,wdays),giveCsparse=F)
target = accesslogdata_sum$pageview

# Subset data to training and test data
set.seed(123)
subset = sample.int(nrow(data),nrow(data)*.8)
data.train = data[subset,]
data.test = data[-subset,]
target.train = target[subset]
target.test = target[-subset]

# Predict ratings with Support Vector Machine with linear kernel
model = SVM.train(data.train,target.train)

# RMSE resulting from test data prediction
sqrt(mean((predict(model,data.test)-target.test)^2))

# Predict ratings with second-order Factorization Machine
# with second-order 10 factors (default) and regularization
model = FM.train(data.train,target.train,regular=0.1)
# RMSE resulting from test data prediction
sqrt(mean((predict(model,data.test)-target.test)^2))

library(FactoRizationMachines)

library(Matrix)

user = accesslogdata_sum$user_number

items = accesslogdata_sum$article_number + max(user)

wdays = (as.POSIXlt(accesslogdata_sum$date,origin="2017-07-04")$wday+1)+max(items)

# Transform access log to feature form

data = sparseMatrix(i=rep(1:nrow(accesslogdata_sum),3),j=c(user,items,wdays),giveCsparse=F)

target = accesslogdata_sum$pageview

# Subset data to training and test data

set.seed(123)

subset = sample.int(nrow(data),nrow(data)*.8)

data.train = data[subset,]

data.test = data[-subset,]

target.train = target[subset]

target.test = target[-subset]

# Predict ratings with Support Vector Machine with linear kernel

model = SVM.train(data.train,target.train)

# RMSE resulting from test data prediction

sqrt(mean((predict(model,data.test)-target.test)^2))

# Predict ratings with second-order Factorization Machine

# with second-order 10 factors (default) and regularization

model = FM.train(data.train,target.train,regular=0.1)

# RMSE resulting from test data prediction

sqrt(mean((predict(model,data.test)-target.test)^2))

結果

各モデルについての平均二乗誤差を計算しています。
線形モデルや高次元モデルよりも、2次の項を持つFMが精度が高いようです。

> model = SVM.train(data.train,target.train)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.514266
> model = FM.train(data.train,target.train,regular=0.1)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.369983
> model　=　HoFM.train(data.train,target.train,c(1,3,1),regular=0.1)
> sqrt(mean((predict(model,data.test)-target.test)^2))
[1] 1.388247

> model = SVM.train(data.train,target.train)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.514266

> model = FM.train(data.train,target.train,regular=0.1)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.369983

> model　=　HoFM.train(data.train,target.train,c(1,3,1),regular=0.1)

> sqrt(mean((predict(model,data.test)-target.test)^2))

[1] 1.388247

こちらは、この中で性能の良かった2次の項を持つFMの予測結果とテストデータの結果をプロットしたものです。4点を超える値をあまり予測できていないようです。今回はサンプルを回しただけなので、本来であれば次元の数kや正則化のセッティングをいろいろいじったり、相互作用項を新しく追加するなどして精度を高めることが必要です。

結果の比較だけでは仕事で使えないので、実際に予測した結果を取り出したいと思います。
実際に運用するとなると、ページIDを所与として、ページビュー数を0とおいて（型をそろえるため。NULLだとエラーになった）、任意のタイミング(date)を想定して、モデルにデータを適用し、評価の高いものをサジェストするスタイルになるのではないでしょうか。

この結果だと、ユーザー98に記事1を見せることに対して4.02点が与えられています。

参考文献

Factorization Machinesを今更読みました
 Factorization Machines
High-order factorization machines with R #tokyor 61
Factorization Machinesのおはなし。
libFMexeを動かすまで (R Wrapper for the libFM Executable参照記事)
一歩Matrix Factorization、二歩Factorization Machines、三歩Field-aware Factorization Machines…『分解、三段突き！！』
[論文] Factorization Machines (ICDM 2010) 読んだ 22:41
Factorization machines with r
Factorization Machinesについて調べてみた