統計コンサルの議事メモ

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

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

先日、Marketing Mix Modelingにおける広告の累積効果(Ad-Stock効果)の推定について以下のような記事を書いた。

ushi-goroshi.hatenablog.com

その後も推定方法について調べていたところ、以下のような記事を発見した:

www.mm-lab.jp

要するに一期前の目的変数そのもの(y_{t-1})を説明変数として含めることで、過去の累積効果、すなわちAd-Stock効果を説明することができるらしい。これを発表者の名前からパルダ・モデルと呼ぶとのこと。
数式を用いた説明は以下で発見できた:

http://www-cc.gakushuin.ac.jp/~20130021/ecmr/chapter6.pdf

この記事において、コイック型の分布ラグを仮定したモデルはARX(1)と等価であることが解説されている。

それでは、このパルダ・モデルと先日のoptimによる推定ではどちらがより精度よく推定できるのだろうか?

シミュレーションの実行:パルダ・モデル

パルダ・モデルの推定精度を検証するためにシミュレーションを実行した。

まずは先日の記事同様にシミュレーション用のサンプルデータを発生させるための関数を作成する。真のモデルはstats::filterによって生成された説明変数にbetaを乗じたもので、このbetaを元データから推定したい。

simulate_y <- function(pars) {
   n         <- pars[1] # num of observation
   mu        <- pars[2] # intercept
   beta_01   <- pars[3] # regression coefficient of X1 to be esitmated
   beta_02   <- pars[4] # regression coefficient of X2 to be esitmated
   lambda_01 <- pars[5] # decay rate of Ad-Stock effect of X1
   lambda_02 <- pars[6] # decay rate of Ad-Stock effect of X2
   var_e     <- pars[7] # residual variance
   
   X_01_raw <- rnorm(n, 100, 2) * ifelse(runif(n) > 0.7, 0, 1)
   X_01_fil <- stats::filter(X_01_raw, lambda_01, "recursive")

   X_02_raw <- rnorm(n, 100, 2) * ifelse(runif(n) > 0.2, 0, 1)
   X_02_fil <- stats::filter(X_02_raw, lambda_02, "recursive")
   
   error <- rnorm(n, 0, sqrt(var_e))

   y     <- mu + beta_01 * X_01_fil + beta_02 * X_02_fil + error
   dat <- data.frame(
      "Y" = y,
      "X_01" = X_01_raw,
      "X_02" = X_02_raw,
      "Y_lag" = dplyr::lag(y, 1))
   return(na.omit(dat))
}

seedを固定し、500回の反復シミュレーションを実行する。lmの中にY_lagが含められていることに注意。betaの真の値は0.5および0.3と設定した。また(今回は推定の対象ではないが)減衰率はそれぞれ0.9および0.2である。

## simulation
set.seed(123)
init_par      <- array(c(100, 5, 0.5, 0.3, 0.9, 0.2, 5))
k             <- 500
results_palda <- matrix(0, k, 4)
for (i in 1:k) {
   dat <- simulate_y(init_par)
   res <- lm(Y ~ X_01 + X_02 + Y_lag, data = dat)
   results_palda[i, ] <- coef(res)
}
> summary(results_palda)
       V1                 V2               V3               V4        
 Min.   :-15.4731   Min.   :0.4322   Min.   :0.2117   Min.   :0.8300  
 1st Qu.: -1.8896   1st Qu.:0.4851   1st Qu.:0.2900   1st Qu.:0.8741  
 Median :  0.8618   Median :0.5000   Median :0.3049   Median :0.8820  
 Mean   :  1.0207   Mean   :0.5002   Mean   :0.3037   Mean   :0.8820  
 3rd Qu.:  4.0467   3rd Qu.:0.5146   3rd Qu.:0.3175   3rd Qu.:0.8907  
 Max.   : 16.4760   Max.   :0.5597   Max.   :0.3624   Max.   :0.9254  

V2およびV3に注目すると、それぞれ0.43〜0.56、0.21〜0.36の範囲で分布しており、Medianおよび平均は両方とも一致している。どうやらパルダ・モデルは過去の累積による影響をうまく取り除き、広告の当期に対する影響を精度よく推定できるようだ。
なおV4は一期前の目的変数に対するyの回帰係数であり、先に示したpdfによるとこの値は広告の減衰率λに相当するはずである。しかし変数が2つあるためか一致していない。

シミュレーションの実行:optimによる推定

では先日と同様にoptimを使いながら減衰率を同時に推定する方法はどうか。まずは先日と同じく、optimの目的関数を作成する。

return_AIC <- function(param) {
   a <- param[1]
   b <- param[2]
   dat$X_01_fil <- stats::filter(dat$X_01, a, "recursive")
   dat$X_02_fil <- stats::filter(dat$X_02, b, "recursive")
   AIC(lm(Y ~ X_01_fil + X_02_fil, dat))
}

続いてシミュレーションを実行する。反復回数はパルダ・モデルと同じく500回で、betaの値も同じである。実行にやや時間がかかったので進捗を見られるように100回ごとにRoundを表示した。

set.seed(456)
init_par      <- array(c(100, 5, 0.5, 0.3, 0.9, 0.2, 5))
decay_par     <- array(c(0.5, 0.5))
k             <- 500
results_optim <- matrix(0, k, 5)
for (i in 1:k) {
   dat <- simulate_y(init_par)
   res <- optim(par = optim(par = decay_par, fn = return_AIC)$par, 
                fn = return_AIC)$par
   
   dat$X_01_fil <- stats::filter(dat$X_01, res[1], "recursive")
   dat$X_02_fil <- stats::filter(dat$X_02, res[2], "recursive")

   lm_out             <- lm(Y ~ X_01_fil + X_02_fil, data = dat)
   results_optim[i, ] <- c(res, coef(lm_out))
   if (i %% 50 == 0) cat("Round: ", i, "\n")
}
> summary(results_optim)
       V1               V2                V3               V4               V5        
 Min.   :0.8866   Min.   :0.06162   Min.   : 1.983   Min.   :0.4608   Min.   :0.2649  
 1st Qu.:0.8916   1st Qu.:0.17331   1st Qu.: 6.280   1st Qu.:0.4838   1st Qu.:0.2941  
 Median :0.8941   Median :0.19419   Median :35.109   Median :0.4923   Median :0.3001  
 Mean   :0.8949   Mean   :0.19289   Mean   :27.462   Mean   :0.4908   Mean   :0.3002  
 3rd Qu.:0.8996   3rd Qu.:0.21262   3rd Qu.:38.388   3rd Qu.:0.4989   3rd Qu.:0.3059  
 Max.   :0.9010   Max.   :0.29186   Max.   :46.506   Max.   :0.5200   Max.   :0.3286  

V4およびV5を見ると、それぞれの範囲は0.46〜0.52および0.26〜0.33と、パルダ・モデルと比較して精度よく推定できている。加えてoptimを用いる方法では各変数の減衰率λの値も推定することができ、その精度も十分高い。

まとめ

上記のシミュレーションの結果をまとめると:

  1. 目的変数の一期前の数値を説明変数として用いる(パルダ・モデル)ことで、Ad-Stock効果をモデルに取り入れることができる
  2. 簡単なシミュレーションではパルダ・モデルの有効性が示され、モデルを非線形にすることなく精度の良い回帰係数を推定することができた
  3. ただしoptimによる推定の方が回帰係数の推定精度が高く、またパルダ・モデルでは各変数の減衰率λを同時に得ることができない

となった。
パルダ・モデルではモデルを非線形にする必要がないためlmにより簡単に推定でき、またその推定精度も高い。さらに過去の累積効果をまとめて説明することができるため、Marketing Mix Modelの構築において最適なLagを気にする必要がなくなることは大きなメリットである。そのため一般的なMMMの構築においてはパルダ・モデルの方が望ましいのではないか。
一方で、パルダ・モデルは各変数のλを得ることができないため詳細な考察が必要な場合には向いておらず、また分布ラグの形もコイック型に限られてしまう。例えば広告間の効率を比較するためにλを示す必要があったり、使用期限が定められたクーポンのように幾何級数的な減衰が適当でない(期限ギリギリになると効果が上昇する)ことが想定されるような場合にはoptimによる推定の方が適当であると考えられる。

したがって分析者が必要とする数値および対象領域における仮定を踏まえて推定方法を選択する必要があるだろう。

Keras for R

RstudioがR上でKerasによるディープラーニングのモデルを構築するためのライブラリ{keras}を公開した。

R Interface to Keras • keras

以前から{tensorflow}を使えばtensorflow::import(module = "keras")でKerasを導入することができたようだが、{keras}を先にインストールすることでpythonさえ入っていればtensorflowのインストールもできるらしい。随分と楽にtensorflowを動かせることができそうなので早速試してみた。なお類似のプロジェクトがRstudio以外でも行われているようで、以下でもKerasをRから動かすためのライブラリを提供している。

github.com


実行環境はAWS上のWindowsマシン。Rのバージョンは3.3.2だけれど、Microsoft R Clientがインストールされている。

実行したスクリプトは上で紹介したページの内容なのだけれど、サンプルデータであるx_trainやx_testなどは自前で用意する必要があったのでirisを使って適当に作成した。なおこちらのページはKerasの公式ドキュメントをR用に書き換えただけの様子。

Keras Documentation

{Keras}のインストー

早速{keras}のインストールとtensorflowのインストールを実行してみる。なお私の環境では{keras}のインストールの際に{reticulate}が一緒に入らず、{keras}のインストールに失敗したため{reticulate}だけ先にインストールしている。

## library install & load
devtools::install_github("rstudio/reticulate")
devtools::install_github("rstudio/keras")
library(keras)
library(dplyr) # for data handling

## install tensorflow
install_tensorflow()

Error: Installing TensorFlow requires a 64-bit version of Python 3.5

Please install 64-bit Python 3.5 to continue, supported versions include:
   
- Anaconda Python (Recommended): https://www.continuum.io/downloads#windows
- Python Software Foundation   : https://www.python.org/downloads/release/python-353/
   
Note that if you install from Python Software Foundation you must install exactly
Python 3.5 (as opposed to 3.6 or higher).

{keras}のインストールがうまくいけばinstall_tensorflow()という関数一つでtensorflowをインストールすることができる。とはいえpythonのインストールまではやってくれないようで私の環境ではエラーが出てしまった。指示に従い、Anaconda経由でpythonをインストールする。なお現時点でのpythonのバージョンは3.6だった。

Anacondaのインストールが済んだので再度tensorflowのインストールを試してみる。

## install tensorflow
install_tensorflow()

今度はうまくいった。

サンプルデータの作成

サンプルスクリプトでは分析用のデータが提供されていなかったため、irisを使ってサンプルデータを作成する。特別なことはしておらず、sampleによって乱数を振った上で学習用とテスト用にデータを分割した。なおデータはmatrixにする必要がある。

set.seed(123)
learn_id <- sample(1:nrow(iris), 100, replace = FALSE)

train_data <- iris %>%
   filter(row.names(.) %in% learn_id)
test_data  <- iris %>%
   filter(!row.names(.) %in% learn_id)

## training data
x_train <- train_data %>% 
   select(-Species) %>%
   as.matrix()
y_train <- as.matrix(model.matrix(~ train_data$Species - 1))

## test data
x_test <- test_data %>% 
   select(-Species) %>%
   as.matrix()
y_test <- as.matrix(model.matrix(~ test_data$Species - 1))

Kerasによるmodelオブジェクトの生成

サンプルデータが出来上がったので、{keras}によるフィッティングを行う。まずはmodelオブジェクトの生成から。

## generate model object
model <- keras_model_sequential() 

> Model
Model
_______________________________________________________________________________
Layer (type)                       Output Shape                    Param #     
===============================================================================
Total params: 0
Trainable params: 0
Non-trainable params: 0
_______________________________________________________________________________

入力層 → 隠れ層 → 出力層という簡単なニューラルネットワークを構築する。layer_denseで簡単にレイヤーを積み上げることができるようで、入力層のみinput_shapeを指定する必要がある。これはpython上でKerasを実行する際のmodel.addに相当するものと思われる。

model %>% 
   layer_dense(units = 10, input_shape = 4) %>% # Output & Input dimension
   layer_activation(activation = 'relu') %>% 
   layer_dense(units = 3) %>% 
   layer_activation(activation = 'softmax')

> model
Model
_______________________________________________________________________________
Layer (type)                       Output Shape                    Param #     
===============================================================================
dense_3 (Dense)                    (None, 10)                      50          
_______________________________________________________________________________
activation_3 (Activation)          (None, 10)                      0           
_______________________________________________________________________________
dense_4 (Dense)                    (None, 3)                       33          
_______________________________________________________________________________
activation_4 (Activation)          (None, 3)                       0           
===============================================================================
Total params: 83.0
Trainable params: 83.0
Non-trainable params: 0.0
_______________________________________________________________________________

パラメータの数は層ごとに(入力データの次元数+1) * 出力データの次元で決まる。+1は切片。

ここで出来上がったmodelオブジェクトの属性や構造などを見てみようとstrを実行してみたのだが、print(model)と結果は変わらなかった。やはりR的なオブジェクトではない様子。

ニューラルネットワークの構造が決まれば、最適化の際の条件を指定する。

## model parameter
model %>% 
   compile(
      loss = 'categorical_crossentropy',
      optimizer = optimizer_sgd(lr = 0.02),
      metrics = c('accuracy'))

フィッティング

以上で準備が整ったため、フィッティングを行う。

## fitting
model %>% 
   fit(x_train, y_train, epochs = 1000, batch_size = 32)

## 省略

# Epoch 990/1000
# 100/100 [==============================] - 0s - loss: 0.0432 - acc: 0.9800     
# Epoch 991/1000
# 100/100 [==============================] - 0s - loss: 0.0442 - acc: 0.9800     
# Epoch 992/1000
# 100/100 [==============================] - 0s - loss: 0.0455 - acc: 0.9800     
# Epoch 993/1000
# 100/100 [==============================] - 0s - loss: 0.0508 - acc: 0.9800     
# Epoch 994/1000
# 100/100 [==============================] - 0s - loss: 0.0602 - acc: 0.9800     
# Epoch 995/1000
# 100/100 [==============================] - 0s - loss: 0.0437 - acc: 0.9800     
# Epoch 996/1000
# 100/100 [==============================] - 0s - loss: 0.0433 - acc: 0.9800     
# Epoch 997/1000
# 100/100 [==============================] - 0s - loss: 0.0432 - acc: 0.9800     
# Epoch 998/1000
# 100/100 [==============================] - 0s - loss: 0.1129 - acc: 0.9600     
# Epoch 999/1000
# 100/100 [==============================] - 0s - loss: 0.0529 - acc: 0.9700     
# Epoch 1000/1000
# 100/100 [==============================] - 0s - loss: 0.0647 - acc: 0.9600 

サンプルデータがものの100行程度しかないため一瞬で終わる。Accuracyが0.96とまずまず。続いてこちらのモデルを用いて、テストデータに対する精度評価を実行する。

## test
loss_and_metrics <- model %>% 
   evaluate(x_test, y_test, batch_size = 128)

> loss_and_metrics
[[1]]
[1] 0.1853937

[[2]]
[1] 0.96

テストデータに対しても精度が0.96と高く、過学習を起こしている様子はない。
テストデータに対するラベルの予測は以下のように実行する。

## predict
classes <- model %>% 
   predict(x_test, batch_size = 128)

> classes
[,1]         [,2]         [,3]
[1,] 9.983064e-01 1.693639e-03 6.461022e-20
[2,] 9.971897e-01 2.810343e-03 7.951551e-19
[3,] 9.983315e-01 1.668578e-03 3.816257e-20
[4,] 9.984826e-01 1.517428e-03 4.521292e-20
[5,] 9.994913e-01 5.087039e-04 3.988174e-23
## 省略

サンプルスクリプトは以上である。
注意点として、パイプ演算子(%>%)が使えるからといってdplyrが読み込まれている訳ではないので、データ加工を行うのであれば別途読み込む必要がある。

Rでディープラーニングを簡単に実行できる環境が整ってきているので、これを機会にRユーザーが増えると良いな。

帰無仮説は採択できない

統計的仮説検定は通常、以下のような手順に従って行われる:

  1. 帰無仮説を設定する
  2. 対立仮説を設定する
  3. 検定統計量を設定する
  4. 有意水準を設定する
  5. 実験やデータ解析によって検定統計量を求める
  6. 帰無仮説を真と仮定した時の検定統計量の得られる確率を求める
  7. 設定した有意水準に従い、判断を下す

このうち7番で行われる判断において、帰無仮説が棄却された場合は「対立仮説を採択する」という表現がなされる。では対立仮説が棄却された場合は「帰無仮説を採択する」と言うのかといえばそうではない。帰無仮説は採択できないからだ。

そもそも統計的仮説検定は「帰無仮説が真であるとき、この統計量が偶然に得られる確率はどの程度であろうか?」を定量的に評価するための手続きである。大事なポイントとしては「帰無仮説が真である」という点で、この前提の上でどの程度の検定統計量が得られれば帰無仮説を棄却できるか、ということを見ている。帰無仮説と対立仮説は排反であるため、帰無仮説が棄却された場合は自動的に対立仮説が採択される。

では帰無仮説が棄却されなかった場合はと言うと、これがややこしいところで、「帰無仮説を採択する」ことはできないのである。「棄却できない」と「採択する」は異なるからだ。

先にも述べたように、統計的仮説検定はあくまで「帰無仮説を棄却できるか」を検証するための手続きである。この手続きから得られる結果は「帰無仮説を棄却できる」「できない」の2つだけであり、採択に関する考慮は一切ない。そのため、例えば回帰分析を行った際にβ = 1.5, P < 0.01という回帰係数が得られたとしても、ここで言えることは「回帰係数が0であること(=帰無仮説)が棄却されたので、回帰係数は0ではない」ということだけであり、1.5という数値そのものに対しては何も言えないのである。

ちなみに検定統計量(t値やχ2値など)の大小をもって「極めて有意」「わずかに有意」などといった定量的は判断が下されることがあるが、厳密にはこれはよろしいものではない。伝統的な手続きに則った場合、棄却できるか否か(1か0か)の判断だけが可能なのである。
(念のため付け加えると、もちろん「t値が大きく、めったに起こらないことが発生した」といった表現は問題なく、「棄却できるか否か」について定量的な表現がよくない、ということである。最初に設定した有意水準を超えるか否かだけが焦点であり、5%と設定したのであれば、P値が0.04であろうと0.000001であろうと等しく「棄却できる」のである。)

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

最近ずっと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の注意点

最近すっかり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の罠

年月を6桁の数値(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のインストールとロジスティック回帰の実行

{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