統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

Ad-Stock効果を推定しつつ回帰を回したい

MMM マーケティング 最適化 R マーケティングミックスモデル

最近ずっとMarketing Mix Modeling（MMM）をやっている。その中で広告効果（いわゆるROI）を推定しているのだけれど、広告の効果というものは出稿した時点だけでなく将来に渡って影響を及ぼすため、過去の広告の累積による影響（いわゆる残存効果・Ad-Stock効果）を取り上げる必要がある。ところが、このAd-Stock効果をモデルにより推定しようとすると少し難しい。

Ad-Stock効果として以下のような分布ラグ型：

{ \displaystyle y_t = \mu + \beta_{t} X_{t} + \beta_{t-1} X_{t-1} + \cdots + \beta_{t-m} X_{t-m} + \epsilon_t }

を仮定するのが最も単純なのだけれど、考慮するラグの次数mによってはモデルに投入する変数の数が多くなり、また多重共線性も強くなりがち。時系列データを用いているためデータポイントが少なく、このような方法は避けたい。

そこで「m期前の広告効果はλのm乗に従って減衰する」という仮定をおき（いわゆるKoyck型）、推定すべきパラメータの数を減らしている。この場合変数の数は減らせるが、非線形な効果であるλをどのように推定するかが鍵となる。

考えられる方法として、

  1. データ加工の段階で例えば相関係数などから推定する
  2. 過去の研究や業界での知見を導入する
  3. 非線形な効果を推定する方法を考える

といったものが思いついた。しかし1は季節性や他の広告の効果などを同時に考慮できないため✖️、2は使えそうな数値がなかった。そこで今回は回帰の中でλを同時に推定できないか、以下のように調査してみた。

set.seed(123)
# λに従って広告の効果が減衰するデータを生成する関数
simulate_y_01 <- function(pars) {
   X1 <- rnorm(100, 100, 2) * rep(c(0, rep(1, 7), 0, 1), 10)
   X2 <- filter(X1, pars[1], "recursive")
   Z1 <- rnorm(100, 5, 2) * rep(c(rep(0, 5), rep(1, 5)), 10)
   dat <- data.frame(
      "Y"  = 50 + pars[2] * X2 + pars[3] * Z1 + rnorm(100, 0, pars[4]),
      "X1" = X1,
      "Z1" = Z1)
   return(dat)
}

# λを与えるとAd-Stock変数を作成して回帰を実施し、AICを返す関数
return_AIC_01 <- function(a) {
   dat$Tmp <- filter(dat$X1, a, "recursive")
   AIC(lm(Y ~ Tmp + Z1, dat))
}

Koyck型の変数の作成にはstats::filterが便利だが、使う際には{dplyr}を読み込んでいないか注意しよう。これでreturn_AIC_01を最小化した時に得られるλが、simulate_y_01で用いた数値を一致するか確認する。実際にやってみると以下の結果が得られる。

## λを小さな値とした場合
pars      <- c(0.1, 0.5, 4, 5.0)
n         <- 500
res_01_01 <- matrix(0, n, 4)
for (i in 1:n) {
   dat <- simulate_y_01(pars)
   a   <- optimise(return_AIC_01, c(0, 1))$minimum
   X2  <- filter(dat$X1, a, "recursive")
   res_01_01[i, ] <- c(a, coef(lm(Y ~ X2 + Z1, data = dat)))
}

## λを大きな値とした場合
pars      <- c(0.9, 0.5, 4, 5.0)
n         <- 500
res_01_02 <- matrix(0, n, 4)
for (i in 1:n) {
   dat <- simulate_y_01(pars)
   a   <- optimise(return_AIC_01, c(0, 1))$minimum
   X2  <- filter(dat$X1, a, "recursive")
   res_01_02[i, ] <- c(a, coef(lm(Y ~ X2 + Z1, data = dat)))
}

> summary(res_01_01)
       V1                V2              V3               V4       
 Min.   :0.01545   Min.   :43.81   Min.   :0.4633   Min.   :3.486  
 1st Qu.:0.08780   1st Qu.:48.73   1st Qu.:0.4918   1st Qu.:3.877  
 Median :0.10450   Median :49.90   Median :0.4998   Median :3.999  
 Mean   :0.10324   Mean   :49.88   Mean   :0.4996   Mean   :3.997  
 3rd Qu.:0.11992   3rd Qu.:51.02   3rd Qu.:0.5081   3rd Qu.:4.109  
 Max.   :0.17782   Max.   :55.08   Max.   :0.5408   Max.   :4.474  

> summary(res_01_02)
       V1               V2              V3               V4       
 Min.   :0.8972   Min.   :40.15   Min.   :0.4768   Min.   :3.414  
 1st Qu.:0.8992   1st Qu.:47.93   1st Qu.:0.4946   1st Qu.:3.869  
 Median :0.9000   Median :50.04   Median :0.4995   Median :3.994  
 Mean   :0.9000   Mean   :49.94   Mean   :0.5002   Mean   :3.999  
 3rd Qu.:0.9007   3rd Qu.:52.10   3rd Qu.:0.5060   3rd Qu.:4.133  
 Max.   :0.9031   Max.   :58.68   Max.   :0.5236   Max.   :4.526

λが小さい時にはある程度バラついているものの、この誤差なら十分に実用的に思える。λが大きい時は全く問題がなさそうだ。
少し問題を難しくするため、2つの変数を共にAd-Stockに変更すると：

simulate_y_02 <- function(pars) {
   X1 <- rnorm(100, 100, 2) * rep(c(0, rep(1, 7), 0, 1), 10)
   X2 <- filter(X1, pars[1], "recursive")
   Z1 <- rnorm(100, 5, 2) * rep(c(rep(0, 5), rep(1, 5)), 10)
   Z2 <- filter(Z1, pars[2], "recursive")
   dat <- data.frame(
      "Y"  = 50 + pars[3] * X2 + pars[4] * Z2 + rnorm(100, 0, pars[5]),
      "X1" = X1,
      "Z1" = Z1)
   return(dat)
}

return_AIC_02 <- function(param) {
   a <- param[1]
   b <- param[2]
   dat$TmpX <- filter(dat$X1, a, "recursive")
   dat$TmpZ <- filter(dat$Z1, b, "recursive")
   AIC(lm(Y ~ TmpX + TmpZ, dat))
}

## λを小さな値とした場合
pars      <- c(0.1, 0.1, 0.5, 4, 5.0)
n         <- 100
res_02_01 <- matrix(0, n, 5)
for (i in 1:n) {
   dat <- simulate_y(pars)
   res <- optim(par = optim(par = init_par, fn = return_AIC)$par, 
                fn = return_AIC)$par
   X2  <- filter(dat$X1, res[1], "recursive")
   Z2  <- filter(dat$Z1, res[2], "recursive")
   res_02_01[i, ] <- c(res, coef(lm(Y ~ X2 + Z2, data = dat)))
}

## λを大きな値とした場合
pars      <- c(0.8, 0.7, 0.5, 4, 5.0)
n         <- 100
res_02_02 <- matrix(0, n, 5)
for (i in 1:n) {
   dat <- simulate_y(pars)
   res <- optim(par = optim(par = init_par, fn = return_AIC)$par, 
                fn = return_AIC)$par
   X2  <- filter(dat$X1, res[1], "recursive")
   Z2  <- filter(dat$Z1, res[2], "recursive")
   res_02_02[i, ] <- c(res, coef(lm(Y ~ X2 + Z2, data = dat)))
}

計算に少し時間がかかったので、反復は100回とした。結果は以下の通り。

> summary(res_02_01)
       V1                V2                 V3              V4               V5       
 Min.   :0.03261   Min.   :-0.02245   Min.   :44.54   Min.   :0.4685   Min.   :3.597  
 1st Qu.:0.08483   1st Qu.: 0.06107   1st Qu.:48.62   1st Qu.:0.4899   1st Qu.:3.902  
 Median :0.09598   Median : 0.09694   Median :50.23   Median :0.4997   Median :3.997  
 Mean   :0.09810   Mean   : 0.09459   Mean   :50.06   Mean   :0.4997   Mean   :4.017  
 3rd Qu.:0.11280   3rd Qu.: 0.11921   3rd Qu.:51.30   3rd Qu.:0.5081   3rd Qu.:4.122  
 Max.   :0.14426   Max.   : 0.24930   Max.   :55.33   Max.   :0.5419   Max.   :4.606  

> summary(res_02_02)
       V1               V2               V3              V4               V5       
 Min.   :0.7894   Min.   :0.6499   Min.   :41.16   Min.   :0.4630   Min.   :3.701  
 1st Qu.:0.7961   1st Qu.:0.6833   1st Qu.:48.19   1st Qu.:0.4905   1st Qu.:3.906  
 Median :0.7990   Median :0.7003   Median :50.51   Median :0.5016   Median :3.984  
 Mean   :0.7997   Mean   :0.7000   Mean   :49.99   Mean   :0.5006   Mean   :3.994  
 3rd Qu.:0.8030   3rd Qu.:0.7140   3rd Qu.:52.09   3rd Qu.:0.5093   3rd Qu.:4.065  
 Max.   :0.8137   Max.   :0.7496   Max.   :57.14   Max.   :0.5257   Max.   :4.287

先ほどと変わらず、λが小さいと推定の精度が落ちるようだ。また一部のデータセットではλが負となってしまった。

ついでなので、{rBayesianOptimization}パッケージによる最適化を通常のoptimと比較してみた。

library(rBayesianOptimization)
# BayesianOptimization関数では最大化を行うので、返り値となるAICに負の符号を加えてある
return_AIC_BO <- function(a, b) {
   dat$TmpX <- filter(dat$X1, a, "recursive")
   dat$TmpZ <- filter(dat$Z1, b, "recursive")
   out <- list(Score = -AIC(lm(Y ~ TmpX + TmpZ, dat)),
               Pred  = -AIC(lm(Y ~ TmpX + TmpZ, dat)))
   return(out)
}

pars   <- c(0.1, 0.1, 0.5, 4, 5.0)
dat    <- simulate_y(pars)
res_BO <- BayesianOptimization(return_AIC_BO, bounds = list(a = c(0, 1), b = c(0, 1)),
                               init_points = 20, n_iter = 1, acq = "ucb",
                               kappa = 2.576, eps = 0, verbose = TRUE)

init_par <- array(c(0.5, 0.5))
res_Opt  <- optim(par = optim(par = init_par, fn = return_AIC)$par, 
                  fn = return_AIC)$par

> res_BO$Best_Par
        a         b 
0.8274328 0.4757574 
> res_Opt
[1] 0.7959288 0.7930622

BayesianOptimizationの結果はあまり良くなく、通常のoptimの方が精度よく推定できた。以前の記事（optimってあんまり信用できないなぁ、って話 - 統計コンサルの議事メモ）ではoptimによるパラメータの推定はうまくいかなかったのだけど、lmと組み合わせているのが良いのかもしれない。

この方法でならAd-Stock効果の推定できそうだ。optimでは"L-BFGS-B"を指定すれば上限・下限を指定できたと思うので、もう少し検証してみよう。

dplyr::filterの注意点

R

最近すっかりHadley信者になってしまいデータ加工にもdplyrをよく使っているのだけれど、filterで少し躓いてしまったのでメモ。

まずは{dplyr}と{dtplyr}を読み込む：

library(dplyr)
library(dtplyr)

やりたかった処理とは以下のようなもので、irisを例とするとSpeciesごとにfilter内での条件を変えたかった。

> iris %>% 
+    filter(Species == "setosa", Sepal.Length >= 3 | 
+              Species == "versicolor", Sepal.Length >= 1 | 
+              Species == "virginica", Sepal.Length >= 2) %>% tbl_dt()

Source: local data table [50 x 5]

# tbl_dt [50 × 5]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
1           5.1         3.5          1.4         0.2  setosa
2           4.9         3.0          1.4         0.2  setosa
3           4.7         3.2          1.3         0.2  setosa
4           4.6         3.1          1.5         0.2  setosa
5           5.0         3.6          1.4         0.2  setosa
6           5.4         3.9          1.7         0.4  setosa
7           4.6         3.4          1.4         0.3  setosa
8           5.0         3.4          1.5         0.2  setosa
9           4.4         2.9          1.4         0.2  setosa
10          4.9         3.1          1.5         0.1  setosa
# ... with 40 more rows

しかしこのOutputにはsetosaしか残っておらず、どうやらOR条件("|")以下が評価されていないようなのである。。

{dplyr}のfilterでは"&"もAND条件として使用可能なので、そちらで試してみるとうまく行く。

> iris %>% 
+    filter(Species == "setosa" & Sepal.Length >= 3 | 
+              Species == "versicolor" & Sepal.Length >= 1 | 
+              Species == "virginica" & Sepal.Length >= 2) %>% 
tbl_dt()
Source: local data table [150 x 5]

# tbl_dt [150 × 5]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
1           5.1         3.5          1.4         0.2  setosa
2           4.9         3.0          1.4         0.2  setosa
3           4.7         3.2          1.3         0.2  setosa
4           4.6         3.1          1.5         0.2  setosa
5           5.0         3.6          1.4         0.2  setosa
6           5.4         3.9          1.7         0.4  setosa
7           4.6         3.4          1.4         0.3  setosa
8           5.0         3.4          1.5         0.2  setosa
9           4.4         2.9          1.4         0.2  setosa
10          4.9         3.1          1.5         0.1  setosa
# ... with 140 more rows

数値などは一切変えていないが、150行抽出されている。
しかし以下の2つは同じ結果を返すので、","と"&"だけの問題ではなさそう。

# ","でつなぐ
> iris %>% 
+    filter(Species == "setosa" , Sepal.Length >= 5) %>% tbl_dt()
Source: local data table [30 x 5]

# tbl_dt [30 × 5]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
1           5.1         3.5          1.4         0.2  setosa
2           5.0         3.6          1.4         0.2  setosa
3           5.4         3.9          1.7         0.4  setosa
4           5.0         3.4          1.5         0.2  setosa
5           5.4         3.7          1.5         0.2  setosa
6           5.8         4.0          1.2         0.2  setosa
7           5.7         4.4          1.5         0.4  setosa
8           5.4         3.9          1.3         0.4  setosa
9           5.1         3.5          1.4         0.3  setosa
10          5.7         3.8          1.7         0.3  setosa
# ... with 20 more rows

# "&"でつなぐ
> iris %>% 
+    filter(Species == "setosa" & Sepal.Length >= 5) %>% tbl_dt()
Source: local data table [30 x 5]

# tbl_dt [30 × 5]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
1           5.1         3.5          1.4         0.2  setosa
2           5.0         3.6          1.4         0.2  setosa
3           5.4         3.9          1.7         0.4  setosa
4           5.0         3.4          1.5         0.2  setosa
5           5.4         3.7          1.5         0.2  setosa
6           5.8         4.0          1.2         0.2  setosa
7           5.7         4.4          1.5         0.4  setosa
8           5.4         3.9          1.3         0.4  setosa
9           5.1         3.5          1.4         0.3  setosa
10          5.7         3.8          1.7         0.3  setosa
# ... with 20 more rows

"|"も組み合わせて使うと異なる結果が返るのだろうか。

> iris %>% 
+    filter(Species == "setosa" , Sepal.Length >= 5 | Species == "virginica") %>%
+    tbl_dt()

Source: local data table [30 x 5]

# tbl_dt [30 × 5]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
1           5.1         3.5          1.4         0.2  setosa
2           5.0         3.6          1.4         0.2  setosa
3           5.4         3.9          1.7         0.4  setosa
4           5.0         3.4          1.5         0.2  setosa
5           5.4         3.7          1.5         0.2  setosa
6           5.8         4.0          1.2         0.2  setosa
7           5.7         4.4          1.5         0.4  setosa
8           5.4         3.9          1.3         0.4  setosa
9           5.1         3.5          1.4         0.3  setosa
10          5.7         3.8          1.7         0.3  setosa
# ... with 20 more rows

> iris %>% 
+    filter(Species == "setosa" & Sepal.Length >= 5 | Species == "virginica") %>%
+    tbl_dt()
Source: local data table [80 x 5]

# tbl_dt [80 × 5]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
1           5.1         3.5          1.4         0.2  setosa
2           5.0         3.6          1.4         0.2  setosa
3           5.4         3.9          1.7         0.4  setosa
4           5.0         3.4          1.5         0.2  setosa
5           5.4         3.7          1.5         0.2  setosa
6           5.8         4.0          1.2         0.2  setosa
7           5.7         4.4          1.5         0.4  setosa
8           5.4         3.9          1.3         0.4  setosa
9           5.1         3.5          1.4         0.3  setosa
10          5.7         3.8          1.7         0.3  setosa
# ... with 70 more rows

ちなみにfilter内ではAND/OR条件の優先順位を付けるために括弧を使うことができるが、以下のように","を括弧で括る書き方はエラーとなる：

# ","を括弧でくくるとエラーになる
> iris %>% 
+    filter((Species == "setosa" , Sepal.Length >= 5)) %>%
Error: unexpected ',' in:
"iris %>% 
   filter((Species == "setosa" ,"

# "&"や"|"なら大丈夫
> iris %>% 
+    filter((Species == "setosa" & Sepal.Length >= 5)) %>%
+    tbl_dt()
Source: local data table [30 x 5]

# tbl_dt [30 × 5]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
1           5.1         3.5          1.4         0.2  setosa
2           5.0         3.6          1.4         0.2  setosa
3           5.4         3.9          1.7         0.4  setosa
4           5.0         3.4          1.5         0.2  setosa
5           5.4         3.7          1.5         0.2  setosa
6           5.8         4.0          1.2         0.2  setosa
7           5.7         4.4          1.5         0.4  setosa
8           5.4         3.9          1.3         0.4  setosa
9           5.1         3.5          1.4         0.3  setosa
10          5.7         3.8          1.7         0.3  setosa
# ... with 20 more rows
> 
> iris %>% 
+    filter((Species == "setosa" | Sepal.Length >= 5)) %>%
+    tbl_dt()
Source: local data table [148 x 5]

# tbl_dt [148 × 5]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
          <dbl>       <dbl>        <dbl>       <dbl>  <fctr>
1           5.1         3.5          1.4         0.2  setosa
2           4.9         3.0          1.4         0.2  setosa
3           4.7         3.2          1.3         0.2  setosa
4           4.6         3.1          1.5         0.2  setosa
5           5.0         3.6          1.4         0.2  setosa
6           5.4         3.9          1.7         0.4  setosa
7           4.6         3.4          1.4         0.3  setosa
8           5.0         3.4          1.5         0.2  setosa
9           4.4         2.9          1.4         0.2  setosa
10          4.9         3.1          1.5         0.1  setosa
# ... with 138 more rows

原因が明確になったわけではないが、ひとまずANDの指定には"&"を使ったほうが良さそうだ。

summaryの罠

R

年月を６桁の数値（YYYYMM）で表すために以下のように書いて何気なくsummaryを実行したところ、思わぬ挙動となった。

> Year  <- rep(2012:2016, each=12)
> Month <- rep(1:12, 5)
> YM    <- Year * 100 + Month
> summary(YM)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 201200  201300  201400  201400  201500  201600

上記の通り、月は1〜12の範囲なのにいずれも下2桁が00となっている。
念のため他の関数を使ってみると、ちゃんとした値が返る。

> fivenum(YM)
[1] 201201.0 201303.5 201406.5 201509.5 201612.0
> min(YM)
[1] 201201
> max(YM)
[1] 201612

summaryはgeneric functionなので、今回の例のようなnumeric型に対するsummaryの挙動を確認するためにmethodを調べてみると：

> methods("summary")
 [1] summary.aov                    summary.aovlist*               summary.aspell*               
 [4] summary.check_packages_in_dir* summary.connection             summary.data.frame            
 [7] summary.Date                   summary.default                summary.ecdf*                 
[10] summary.factor                 summary.glm                    summary.infl*                 
[13] summary.lm                     summary.loess*                 summary.manova                
[16] summary.matrix                 summary.mlm*                   summary.nls*                  
[19] summary.packageStatus*         summary.PDF_Dictionary*        summary.PDF_Stream*           
[22] summary.POSIXct                summary.POSIXlt                summary.ppr*                  
[25] summary.prcomp*                summary.princomp*              summary.proc_time             
[28] summary.srcfile                summary.srcref                 summary.stepfun               
[31] summary.stl*                   summary.table                  summary.tukeysmooth*

色々と出てきたが、どうやらsummary.defaultがそれっぽい。
summary.defaultの中身を調べてみる。

> summary.default
…
中略
…
      }) else if (is.numeric(object)) {
         nas <- is.na(object)
         object <- object[!nas]
         qq <- stats::quantile(object)
         qq <- signif(c(qq[1L:3L], mean(object), qq[4L:5L]), digits)
         names(qq) <- c("Min.", "1st Qu.", "Median", "Mean", "3rd Qu.", 
                        "Max.")
         if (any(nas)) 
            c(qq, `NA's` = sum(nas))
         else qq
      } 
…
中略
…

どうやらsummary.defaultの中身は{stats}のquantileで、間にmeanを挟んでいるだけの様子。そこでquantileを直接実行してみると：

> quantile(YM)
      0%      25%      50%      75%     100% 
201201.0 201303.8 201406.5 201509.2 201612.0

それっぽい値が返ってきている。どうやら次の処理のsignifの引数のdigitsが怪しい。
そこで今度はsummaryのdigitsを指定して実行してみる。

> summary(YM, digits = 8)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
201201.0 201303.8 201406.5 201406.5 201509.2 201612.0

狙った通りの挙動となった。どうやら6桁の整数では大きすぎて小数点以下が丸められてしまっていたようだ。
めでたしめでたし・・・・と思ったら、よく見るとsummary（つまりquantile）の返り値とfivenumの返り値が微妙に違っている。

> summary(YM, digits = 8)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
201201.0 201303.8 201406.5 201406.5 201509.2 201612.0 
> fivenum(YM)
[1] 201201.0 201303.5 201406.5 201509.5 201612.0

25パーセンタイルと75パーセンタイルの値が違う。調べてみると、fivenumは四分位数ではなくヒンジを返すのだそう。

http://cse.naro.affrc.go.jp/takezawa/r-tips/r/59.html

なのでデータが偶数個の場合に先ほどのような挙動となる様子。
こんな感じ。

> hoge <- 1:10
> summary(hoge)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    3.25    5.50    5.50    7.75   10.00 
> quantile(hoge)
   0%   25%   50%   75%  100% 
 1.00  3.25  5.50  7.75 10.00 
> fivenum(hoge)
[1]  1.0  3.0  5.5  8.0 10.0　# 合わない！

> foo <- c(0, hoge)
> summary(foo)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     2.5     5.0     5.0     7.5    10.0 
> quantile(foo)
  0%  25%  50%  75% 100% 
 0.0  2.5  5.0  7.5 10.0 
> fivenum(foo)
[1]  0.0  2.5  5.0  7.5 10.0　# 一致した！

SparklyrによるApache Sparkのインストールとロジスティック回帰の実行

R Spark

{sparklyr}というパッケージを使うことでWindowsであってもApache Sparkのインストールが簡単にできる。また{sparklyr}にはSpark MLlibの機械学習用の関数がラップされており、それを使ってみた結果を記しておく。

基本的にはRstudioの{sparklyr}の紹介ページをなぞってみただけだが、せっかくなので前回の記事で紹介した{MicrosoftML}のrxGLMと比較してみた。

http://spark.rstudio.com/

なお当然ながら環境はWIndows。

まずは{sparklyr}のインストール。普通にCRANからインストールできる。

install.packages("sparklyr")
library(sparklyr)

その他、必要なライブラリも合わせて読み込んでおく。

library(dplyr)
library(readr)
library(MicrosoftML)
library(RevoScaleR)

spark_installという関数でSparkがあっさりとインストールできる。なお私の環境（AWSに立ち上げたばかりのWindows Server）ではC++ SP1とJAVAが入っていなかったため以下のようなエラーがでた。

> spark_install(version = "1.6.2")
Error: Running Spark on Windows requires the Microsoft Visual C++ 2010 SP1 Redistributable Package. Please download and install from: 

  https://www.microsoft.com/download/en/details.aspx?id=13523

Restart your rsession after installation completes

もしくは

Error in shell_connection(master = master, spark_home = spark_home, app_name = app_name,  : 
  Java is required to connect to Spark. Please download and install Java from https://www.java.com/en/

エラーメッセージにしたがってC++ SP1とJAVAをそれぞれインストールする。インストール後はRのリスタートだけでなくWindowsの再起動をかけた方が良さそう。

なお{sparklyr}はまだまだ開発途上なので、なるべく最新版を用いた方が変なエラーに悩まされなくてすむようだ。以下のように{sparklyr}を最新版に更新する。なぜか{shiny}とRtoolsが一緒にインストールされた。

devtools::install_github("rstudio/sparklyr")


以上で準備が整ったので、Sparkに接続する。なお私はSparkに触ったことがないため、以降のコマンドは正直よくわからない。

sc <- spark_connect(master = "local")

copy_toによりRのオブジェクト（主にdata.frameが対象か？）をSparkに入れることができる。

install.packages(c("nycflights13", "Lahman"))
iris_tbl    <- copy_to(sc, iris)
flights_tbl <- copy_to(sc, nycflights13::flights, "flights")
batting_tbl <- copy_to(sc, Lahman::Batting, "batting")

Sparkに入れたテーブルはRのdata.frameと同様に、以下のように{dplyr}の関数が適用できる。最後のcollectは{dplyr}の関数でSpark上のテーブルをRに読み込むのだが、この辺りは{dplyr}のLazinessに関係しており、以下を参照すると良い。

Manipulating Data with dplyr

delay <- flights_tbl %>% 
   group_by(tailnum) %>%
   summarise(count = n(), dist = mean(distance), delay = mean(arr_delay)) %>%
   filter(count > 20, dist < 2000, !is.na(delay)) %>%
   collect

ローカルのCSVなどのファイルを読み込む時には、以下のようにspark_read_csvを使う。前回同様、SUSYをサンプルデータとして使用したが、ロードにかかった時間は4分ほどだった。

system.time(
   datTmp <- spark_read_csv(sc, "SUSY", 
                            path = "../MicrosoftML/Data/SUSY.csv",
                            header = FALSE)
)


Sparkに入れるのと、Rに直接読み込むのでは比較しても仕方ないかもしれないが、参考用に{readr}のread_csvで同様にSUSYを読み込んでみた。その結果こちらは10分程度かかった。
なお本来なら上記のようにcollectすれば良いのだが、RStudioがクラッシュしてしまったので仕方なくread_csvしている。もしかしたらas.data.frameの方が早いかもしれない。

system.time(
   datTmpDF <- read_csv("../MicrosoftML/Data/SUSY.csv",
                        col_names = FALSE)
)
   ユーザ   システム       経過  
     56.14      32.89     699.00

読み込んだデータを使い、ロジスティック回帰を実行する。本来ならGLMで比較するべきなのだが、ml_generalized_linear_regressionがSpark2.0以降しか使えないため、やむなくロジスティック回帰を用いている。

system.time(
   spark_logit <- datTmp %>%
      ml_logistic_regression(response = "V1",
                             features = c("V2", "V3", "V4", "V5"))
)
   ユーザ   システム       経過  
      0.07       0.01     270.53

結果は4分半で終了した。500万件の分析であることを考えると十分な速度と思われる。

system.time(
   res_rxGLM  <- rxGlm(X1 ~ X2 + X3 + X4 + X5, 
                       family = binomial("logit"), 
                       data = datTmpDF)
)
Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.271 seconds 

Starting values (iteration 1) time: 1.456 secs.
Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.180 seconds 

Iteration 2 time: 1.229 secs.
Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.556 seconds 

Iteration 3 time: 1.598 secs.
Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.229 seconds 

Iteration 4 time: 1.272 secs.
Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.291 seconds 

Iteration 5 time: 1.331 secs.
Rows Read: 5000000, Total Rows Processed: 5000000, Total Chunk Time: 1.263 seconds 

Iteration 6 time: 1.304 secs.

Elapsed computation time: 8.572 secs.
   ユーザ   システム       経過  
      0.89       1.02      27.23

一方、rxGlmは27秒と、相変わらず凄まじい速度で分析が実行できる。
それぞれの結果を、回帰係数で比較してみるとほぼほぼ同じ結果となった。

> cbind(spark_logit$coefficients,
        res_rxGLM$coefficients)

                     [,1]          [,2]
(Intercept) -1.7100304090 -1.7100301776
V2           2.6339536980  2.6339542674
V3          -0.0004056923 -0.0004060596
V4          -0.0009712714 -0.0009717918
V5          -0.9571219696 -0.9571225996


比較は以上である。前回の記事で書いた通りMRCを使えばかなりの規模のデータに対してRからモデリングを実行できるため、現時点でSparkあるいは{sparklyr}にはそれほど興味は引かれない。しかしRに読み込めないような規模のデータを、Sparkに格納しつつml_generalized_linear_regressionで分析が可能なのであれば、データの規模に応じて分析環境を変えることで分析の幅が広がると感じた。よって次の検証内容としては{sparklyr}を用いてのデータサイズの限界を探るようなものが必要となりそうだ。

なおSparkはH2Oに接続することでH2Oの機械学習用のライブラリが使え、これをH2O Sparkling Waterと言うが、Rでも{sparkling}というパッケージで対応しているのでそこも調査する必要がありそうだ。

最後に、{sparklyr}のチートシートを発見したのでついでにリンクを貼っておく。

http://spark.rstudio.com/images/sparklyr-cheatsheet.pdf

Microsoft R Clientによる大規模データの分析

R

Revolution Analyticsを買収したMicrosoftが、Revolution R Openに代わりMicrosoft R OpenというRのラッパーのようなものを出している。それに更に大規模データ分析用の独自開発パッケージを追加したMicrosoft R Client（MRC）というツールがあり、MRCの独自パッケージ{MicrosoftML}に実装されているrxGLMを用いて簡単な動作テストを実行したのでその結果を記しておく。MRCのインストールについてはこちらを参照：

bigdata4people.blogspot.jp

まずは大規模データを読み込むために{readr}パッケージをインストール。なお今回はAWSに新しく立てたWindows Serverを使うので、{tidyverse}を使って必要そうなパッケージをまとめてインストールすることにする。

install.packages("tidyverse")

続いてライブラリの読み込みと作業ディレクトリの変更を行う。

library(readr)
library(MicrosoftML)
library(RevoScaleR)
library(lattice)

wk_dir <-"C:/Users/Administrator/Desktop/MicrosoftML"
setwd(wk_dir)


データを読み込む。{MicrosoftML}がハイパフォーマンスやスケーラブルを謳っているので、そこそこの大規模データ（500万行×18列、2GB程度）を使用してみた。サンプルデータはこちらから取得した。

UCI Machine Learning Repository

### 大規模データの読み込みにreadr::read_csvを使用
### ヘッダーがないのでcol_namesはFALSEにする
datTmp <- read_csv("Data/SUSY.csv",
                   col_names = FALSE)

### X1が目的変数、X2〜X5を説明変数として使用した
### 全カラムを使用した時の結果は後日

### まずは結果の比較用にglmを実行。まぁ、動かないだろうと予想
start_time <- Sys.time()
res_GLM    <- glm(X1 ~ X2 + X3 + X4 + X5, 
                  family = binomial("logit"), 
                  data = datTmp)
end_time   <- Sys.time()
print(end_time - start_time)
# 予想通りメモリの割り当てに失敗した旨のエラーメッセージを受け取る。。。

### MicrosoftMLではどうか。rxGLMが{MicrosoftML}におけるGLMに相当する関数。
start_time <- Sys.time()
res_rxGLM  <- rxGlm(X1 ~ X2 + X3 + X4 + X5, 
                    family = binomial("logit"), 
                    data = datTmp)
end_time   <- Sys.time()

> print(end_time - start_time)
Time difference of 9.806952 secs

なんと、500万行のデータのglmが10秒足らずで実行できた。なおWindowsインスタンスはAWSでのt2.medium（64bit、RAM4GB）を使用しており、特別スペックが良いという訳ではない。

このままではglmとの比較ができないので、もう少しデータを減らして再挑戦してみる。

### gcを使ってメモリを解放しておく
gc()
gc()

### 使用データを100万行にする
datTmp_1MM <- datTmp[1:1000000, ]

### 再実行
start_time <- Sys.time()
res_rxGLM_02  <- rxGlm(X1 ~ X2 + X3 + X4 + X5, 
                      family = binomial("logit"), 
                      data = datTmp_1MM)
end_time   <- Sys.time()
print(end_time - start_time)
# Time difference of 1.839989 secs

gc()
gc()
start_time <- Sys.time()
res_GLM_02    <- glm(X1 ~ X2 + X3 + X4 + X5, 
                    family = binomial("logit"), 
                    data = datTmp_1MM)
end_time   <- Sys.time()
>print(end_time - start_time)
# Time difference of 7.471947 secs

今度はglmでも10秒かからず実行できた。
そこで両者の結果が一致するかを確認する。

coef(res_GLM_02) 
  (Intercept)            X2            X3            X4            X5 
-1.7128408655  2.6384939818  0.0002822391 -0.0012571389 -0.9573449310 

coef(res_rxGLM_02)
  (Intercept)            X2            X3            X4            X5 
-1.7128408655  2.6384939818  0.0002822391 -0.0012571389 -0.9573449310

回帰係数は両者で一致しており、問題なさそうだ。
続いてsummaryを実行する。

### まずはrxGLM
> summary(res_rxGLM_02)
Call:
rxGlm(formula = X1 ~ X2 + X3 + X4 + X5, data = datTmp_1MM, family = binomial("logit"))

Generalized Linear Model Results for: X1 ~ X2 + X3 + X4 + X5
Data: datTmp_1MM
Dependent variable(s): X1
Total independent variables: 5 
Number of valid observations: 1e+06
Number of missing observations: 0 
Family-link: binomial-logit 
 
Residual deviance: 1160233.9234 (on 999995 degrees of freedom)
 
Coefficients:
              Estimate Std. Error  z value Pr(>|z|)    
(Intercept) -1.7128409  0.0052066 -328.976 2.22e-16 ***
X2           2.6384940  0.0079644  331.285 2.22e-16 ***
X3           0.0002822  0.0022039    0.128    0.898    
X4          -0.0012571  0.0022356   -0.562    0.574    
X5          -0.9573449  0.0062846 -152.332 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Condition number of final variance-covariance matrix: 7.0668 
Number of iterations: 6

### 続いてglm
> summary(res_GLM_02)
Call:
glm(formula = X1 ~ X2 + X3 + X4 + X5, family = binomial("logit"), 
    data = datTmp_1MM)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-8.4904  -0.8894  -0.7140   1.0654   3.0185  

Coefficients:
              Estimate Std. Error  z value Pr(>|z|)    
(Intercept) -1.7128409  0.0052066 -328.977   <2e-16 ***
X2           2.6384940  0.0079644  331.287   <2e-16 ***
X3           0.0002822  0.0022039    0.128    0.898    
X4          -0.0012571  0.0022356   -0.562    0.574    
X5          -0.9573449  0.0062846 -152.333   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1379205  on 999999  degrees of freedom
Residual deviance: 1160234  on 999995  degrees of freedom
AIC: 1160244

Number of Fisher Scoring iterations: 5

summaryの結果も両者で大きく異なる様子はない。rxGLMの方でもう少し情報を追加したぐらいか。

一旦動作テストは以上である。{MicrosoftML}は少なくともオリジナルのRのglmと比較して大規模データに強いというのは間違いなさそうで、Rの弱点をカバーできる作りになっているようだ。従来MicrosoftはSQL ServerとRを連携させることで大規模データへの対応を行ってきたイメージがあったが、R単体（と言えないかもしれないが）でこれぐらいの規模のデータを処理できるのであれば、ウェブログを触るような状況とならない限りR一本で大抵の処理はできるように思う。RユーザーとしてはMRCを含めて近年のMicrosoftのRへの貢献はかなり大きいと感じている。

また今回の記事では割愛したが、{MicrosoftML}ではrxGLM以外にも様々な関数が用意されており、特に二項分類のための”Fast Tree"と呼ばれるアルゴリズムを実装した関数はかなり強力な感触を受けた。こちらの検証結果も後日書いてみようと思う。

人工知能と機械学習とDeep Learning

人工知能 機械学習 ディープラーニング

「人工知能（AI）による業務改革！」「Deep Learningを搭載した○○！」みたいな記事は巷に溢れかえっているが、そういった記事を見るにつけ思うのが「人工知能と機械学習、Deep Learningなどのアルゴリズムの区別ってどうなっているんだろう？」ということである。自分で機械学習のアルゴリズムを触ったことがあればわかるかと思うのだが、現在世間で言われているところの「人工知能」とは、「自律的な思考を行う」だとか「知性を獲得した」だとかそういったものでは全くなく、Deep Learningもそのための特殊なアルゴリズムというわけではない。

この辺りの関係性をもう少しスッキリと説明できないかなーと考えていたところ、以下の記事を発見した：

blog.algorithmia.com


またDeep Learningの大家であるモントリオール大学のBengioも同様に整理している：

http://www.deeplearningbook.org/contents/intro.html


私の人工知能に対する認識も上記と同じで（権威に阿ることにする）、Deep Learningは機械学習の中の1アルゴリズム（いわゆる「弱いAI」）として理解している。そのためDeep Learningを指して「人工知能」と称することにはかなりの抵抗がある。もしそれが成り立つなら、同じく機械学習の1アルゴリズムである（例えば）Random Forestを指して「人工知能」と呼んだって良いわけで、それならば決定木だってアリということになるだろうし、さらには線形回帰だって機械学習の1アルゴリズムとして捉えれば「人工知能」となってしまうだろう（もしかしたら実際に線形回帰程度のものを人工知能と謳っているところがあるかもしれない）。

「人工知能による〜〜」という言葉に惑わされず、内容を吟味した上で評価することを心がけたい。

RでTensorFlow

R ディープラーニング 機械学習 人工知能

知らない間にRでTensorFlowが使えるようになっていたので触ってみました。それにしてもRStudioは相変わらずイイ仕事をしますね。

まずは以下の通り、RStudioのGitHubからTensorFlowのライブラリをインストールします。なおこのライブラリはあくまでPCに事前にインストールされたTensorFlowをRから呼び出すためのものなので、TensorFlow自体は別途インストールする必要があります。

## TensorFlowのライブラリをインストールする
devtools::install_github("rstudio/tensorflow")

## ライブラリの読み込み
## データ加工用にdplyrとcaret、分析結果の比較用にglmnetとrandom forestのライブラリも読み込む
library(tensorflow)
library(dplyr)
library(caret)
library(glmnet)
library(ranger)


無事にインストールが完了したら、一旦TensorFlowのバージョンを確認しておきましょう。ちなみに私の環境は0.11です。

## TensorFlowのバージョン確認
tf$VERSION
[1] "0.11.0rc0"


では早速、TensorFlowによる「Hello,World!」です。この辺りは普通にマニュアルにある通りです。

## hello, tensorflow
sess  <- tf$Session()
hello <- tf$constant("Hello, TensorFlow!")
sess$run(hello)
[1] "Hello, TensorFlow!"

単回帰の実行

続いてTensorFlowによる単回帰を実行してみます。この書き方を理解できれば重回帰への拡張は容易でしょう。

## seedを固定しておく
set.seed(123)

## データ作成。回帰のパラメータは切片が-12、傾きが1.5（適当です）
x_data <- runif(100, min=0, max=1)
y_data <- x_data * 1.5 - 12

## TensorFlowによる変数の宣言
W <- tf$Variable(tf$random_uniform(shape(1L), -1.0, 1.0))
b <- tf$Variable(tf$zeros(shape(1L)))
y <- W * x_data + b

## 損失関数の定義と最適化のパラメータ設定
loss  <- tf$reduce_mean((y-y_data)^2)
opt   <- tf$train$GradientDescentOptimizer(0.2)
train <- opt$minimize(loss)


上記で分析環境が整ったので以下の通り実行しますが、その際の注意点として、TensorFlowのバージョンによって初期化に用いる関数が異なることに気をつけましょう。私の環境（0.11）では初期化の関数は「initialize_all_variables()」なのですが、新しいバージョン（0.12）では「global_variables_initializer()」が用いられるようです。

sess <- tf$Session()

## 変数の初期化。TensorFlowのバージョンによって、初期化に用いる関数が異なる。
# sess$run(tf$global_variables_initializer())
sess$run(tf$initialize_all_variables())

## 学習の実行
for (step in 1:1000){
  sess$run(train)
  if (step %% 20 == 0) {
    cat(step, "-", sess$run(W), sess$run(b), "\n")
  }
}
20 - -2.501043 -9.869881 
40 - -0.8908119 -10.72716 
60 - 0.0713779 -11.23942 
80 - 0.646331 -11.54552 
100 - 0.9898927 -11.72843
：
920 - 1.499992 -12 
940 - 1.499992 -12 
960 - 1.499992 -12 
980 - 1.499992 -12 
1000 - 1.499992 -12

LearningRateは小さめに設定してありますが、ちゃんと収束しているようですね。

多項分類

続けてIrisのデータを用いてSpeciesの分類を行います。

## Irisデータの多項分類
W  <- tf$Variable(tf$zeros(shape(4L, 3L)))
x  <- tf$placeholder(tf$float32, shape(NULL, 4L))
b  <- tf$Variable(tf$zeros(3L))
y  <- tf$nn$softmax(tf$matmul(x, W)+b)
y_ <- tf$placeholder(tf$float32, shape(NULL, 3L))

## 損失関数の定義
loss  <- tf$reduce_mean(-tf$reduce_sum(y_*y, reduction_indices = 1L))
opt   <- tf$train$GradientDescentOptimizer(0.9)
train <- opt$minimize(loss)

## 分析用データ（Iris）
trn_x <- iris %>% 
  select(-Species) %>% 
  as.matrix()
trn_y <- iris %>% 
  dummyVars(formula = ~ Species, sep = NULL) %>% 
  predict(object = ., newdata = iris) %>% 
  as.matrix()

sess <- tf$Session()
sess$run(tf$initialize_all_variables())

## 学習の実行
N    <- 20000
each <- 1000
resAll <- matrix(NA, N, 15)
for (i in 1:N){
  sess$run(train, dict(x=trn_x, y_=trn_y))
  resAll[i, ] <- c(sess$run(W), sess$run(b))
  if (i %% each == 0) {
    cat(i, "-", sess$run(W), sess$run(b), "\n")
    print(table(
      predict = sess$run(fetches = y, feed_dict = dict(x = trn_x, y_ = trn_y)) %>% apply(margen=1, FUN=which.max),
      true = trn_y %>% apply(margen=1, FUN=which.max)))
  }
}
1000 - 1.820591 3.826579 -5.199587 -2.53192 1.656648 0.04738447 -0.9168726 -2.692287 -3.47724 -3.873963 6.116461 5.224204 0.8330443 1.588295 -2.421341 
       true
predict  1  2  3
      1 50  0  0
      2  0 46  0
      3  0  4 50
2000 - 2.170959 4.40071 -6.013396 -2.967016 1.705534 0.348204 -1.100592 -3.198786 -3.876493 -4.748915 7.113993 6.165799 0.9800377 2.52502 -3.505058 
       true
predict  1  2  3
      1 50  0  0
      2  0 47  0
      3  0  3 50
：
19000 - 4.02061 6.796506 -9.696184 -4.990173 2.587524 0.9785848 -2.744421 -5.431025 -6.608002 -7.775012 12.4406 10.42116 1.755378 8.18106 -9.936396 
       true
predict  1  2  3
      1 50  0  0
      2  0 47  0
      3  0  3 50
20000 - 4.081785 6.860996 -9.807202 -5.050861 2.61721 0.9981934 -2.800498 -5.49942 -6.698854 -7.859092 12.60772 10.55024 1.781228 8.357506 -10.13868 
       true
predict  1  2  3
      1 50  0  0
      2  0 47  0
      3  0  3 50


20,000回ほど回してみたのですが、意外と分類が難しいのかパラメータは収束していません。また分類精度も最初の方で現在のレベルまで到達して以降、全く改善はありませんでした。この辺りはLearningRateや初期値の設定などに依存するのかもしれません。

最後に、以上の結果を他の分類器と比較してみます。比較対象にはglmnetとrandom forestを用いました。

## glmnetによる結果と比較してみる
resGLM <- glmnet(trn_x, iris$Species, "multinomial")
table(predict = predict(resGLM, newx = trn_x, type = "class", s = 0.01)[, 1],　true = iris$Species)
            true
predict      setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         47         1
  virginica       0          3        49

## random forestによる結果と比較してみる
resRF <- ranger(Species ~., data=iris, mtry=2, num.trees = 500, write.forest=TRUE)
table(predict = predict(resRF, data=iris)$predictions, true=iris$Species)
            true
predict      setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         50         0
  virginica       0          0        50

random forestは流石の精度ですね。Irisは分類が容易なデータなので、TensorFlow、glmnetもそれなりの精度で分類できています。

雑感

以上、RによるTensorFlowの流れを追ってみたのですが、正直に言えば「コレジャナイ」感が強いなーというのが私の感想です。一番引っ掛かりを感じる部分が「事前にTensorFlowのインストールが必要なこと」で、Rをメインに使用しているユーザーからはPythonの環境整備に対するハードルが高そうだなと思いました。他のライブラリのようにTensorFlow自体もinstall.packages()でインストールできると良いな、と。

もう一つ、これはPython用のTensorFlowをRから呼び出しているためなのでしょうが、書き方がRっぽくないなーと感じました。特に「tf$Session()」の部分などいかにもPythonっぽくてちょっと取っつきにくいなぁと思ってしまいます。

もちろんこれは私がそのように感じるというだけですので、Pythonにも馴染みのある方なら全く問題とはならないでしょう。ただ、私としてはmxnetの方が良さそうだなと思いましたので、もう少し検討を続けようと思います。