統計コンサルの議事メモ

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

Anscombeの例

二変数間の関連性の強さを評価する上で最も良く使われる指標は、間違いなくPearsonの積率相関係数、いわゆる相関だと思う。Excelでも関数により算出でき、関連の強さを-1〜1で評価できるというお手軽さが受け入れられやすいのだと思う。

さて、この相関は直線的な関係性しか評価できないため、必ず散布図とセットで見なくてはならない。表題の「Anscombeの例」とは、統計量が同じであってもデータセットが示す内容は大きく異なることがあるということを明確に教えてくれる良い教材である。実際に見てみよう。

## Load libraries
library(tidyverse)
library(broom)
library(forcats)
library(ggplot2)

まずはRに標準で格納されているデータセットanscombeを使う。

> anscombe
   x1 x2 x3 x4    y1   y2    y3    y4
1  10 10 10  8  8.04 9.14  7.46  6.58
2   8  8  8  8  6.95 8.14  6.77  5.76
3  13 13 13  8  7.58 8.74 12.74  7.71
4   9  9  9  8  8.81 8.77  7.11  8.84
5  11 11 11  8  8.33 9.26  7.81  8.47
6  14 14 14  8  9.96 8.10  8.84  7.04
7   6  6  6  8  7.24 6.13  6.08  5.25
8   4  4  4 19  4.26 3.10  5.39 12.50
9  12 12 12  8 10.84 9.13  8.15  5.56
10  7  7  7  8  4.82 7.26  6.42  7.91
11  5  5  5  8  5.68 4.74  5.73  6.89

anscombeにはこのように、4つのデータセットそれぞれがxとyを持つ。
このままだと少し使いにくいので整形する。

## Anscombe
dataAns <- select(anscombe, "x" = x1, "y" = y1) %>% 
   bind_rows(select(anscombe, "x" = x2, "y" = y2)) %>% 
   bind_rows(select(anscombe, "x" = x3, "y" = y3)) %>% 
   bind_rows(select(anscombe, "x" = x4, "y" = y4)) %>% 
   mutate("dataset" = c(rep("I", nrow(anscombe)), rep("II", nrow(anscombe)),
                        rep("III", nrow(anscombe)), rep("IV", nrow(anscombe))))

出来上がったデータに対して平均と標準偏差を求める。

dataAns %>% 
   group_by(dataset) %>% 
   summarise_each(funs(mean, sd), vars = c(x, y))

# A tibble: 4 × 5
  dataset vars1_mean vars2_mean vars1_sd vars2_sd
    <chr>      <dbl>      <dbl>    <dbl>    <dbl>
1       I          9   7.500909 3.316625 2.031568
2      II          9   7.500909 3.316625 2.031657
3     III          9   7.500000 3.316625 2.030424
4      IV          9   7.500909 3.316625 2.030579

平均と標準偏差が小数点以下2位まで揃っている。それぞれのデータセットについて相関を求めてみると:

dataAns %>% 
   group_by(dataset) %>% 
   summarise("COR" = cor(x, y))

# A tibble: 4 × 2
  dataset       COR
    <chr>     <dbl>
1       I 0.8164205
2      II 0.8162365
3     III 0.8162867
4      IV 0.8165214

これもほぼ同じ値となる。平均、標準偏差、相関が等しいデータセットなのでxとyも同じような値と取っているのだろうと思うのだが、プロットしてみると:

ggplot(dataAns, aes(x = x, y = y, label = dataset)) +
   facet_wrap(~ dataset, nrow = 2) +
   geom_point() +
   geom_smooth(method = "lm", formula = y ~ x, se = FALSE)


f:id:ushi-goroshi:20170708002403p:plain

このように曲線を描くものや外れ値に影響を受けているものがあることがわかる。
なお相関が同じなので当然と言えば当然だが、回帰係数も同じような値となる:

dataAns %>% 
   split(.$dataset) %>% 
   map(., ~ lm(y ~ x, data = .) %>% tidy()) %>% 
   bind_rows()

         term  estimate std.error statistic     p.value
1 (Intercept) 3.0000909 1.1247468  2.667348 0.025734051
2           x 0.5000909 0.1179055  4.241455 0.002169629
3 (Intercept) 3.0009091 1.1253024  2.666758 0.025758941
4           x 0.5000000 0.1179637  4.238590 0.002178816
5 (Intercept) 3.0024545 1.1244812  2.670080 0.025619109
6           x 0.4997273 0.1178777  4.239372 0.002176305
7 (Intercept) 3.0017273 1.1239211  2.670763 0.025590425
8           x 0.4999091 0.1178189  4.243028 0.002164602


このanscombeと同様に、統計量が同程度で全く内容が異なるデータセットとしてdatasaurusというものがある。これを用いて同様に見てみよう。

## Datasaurus
datasaurus <- read_table2("../Data/DatasaurusDozen.tsv",
                          col_names = TRUE, col_types = "cnn") %>%
   mutate(dataset = as_factor(dataset))
> print(datasaurus)
# A tibble: 1,846 × 3
   dataset       x       y
    <fctr>   <dbl>   <dbl>
1     dino 55.3846 97.1795
2     dino 51.5385 96.0256
3     dino 46.1538 94.4872
4     dino 42.8205 91.4103
5     dino 40.7692 88.3333
6     dino 38.7179 84.8718
7     dino 35.6410 79.8718
8     dino 33.0769 77.5641
9     dino 28.9744 74.4872
10    dino 26.1538 71.4103
# ... with 1,836 more rows

このデータをプロットしてみるとどうか。

ggplot(datasaurus, aes(x = x, y = y)) +
   facet_wrap(~ dataset, nrow = 3) +
   geom_point()

f:id:ushi-goroshi:20170708002923p:plain

なんと、全く違う!
これらの全てが本当に同様の統計量を持っているのだろうか?

datasaurus %>% 
   group_by(dataset) %>% 
   summarise_each(funs(mean, sd), vars = c(x, y))
# A tibble: 13 × 5
      dataset vars1_mean vars2_mean vars1_sd vars2_sd
       <fctr>      <dbl>      <dbl>    <dbl>    <dbl>
1        dino   54.26327   47.83225 16.76514 26.93540
2        away   54.26610   47.83472 16.76982 26.93974
3     h_lines   54.26144   47.83025 16.76590 26.93988
4     v_lines   54.26993   47.83699 16.76996 26.93768
5     x_shape   54.26015   47.83972 16.76996 26.93000
6        star   54.26734   47.83955 16.76896 26.93027
7  high_lines   54.26881   47.83545 16.76670 26.94000
8        dots   54.26030   47.83983 16.76774 26.93019
9      circle   54.26732   47.83772 16.76001 26.93004
10   bullseye   54.26873   47.83082 16.76924 26.93573
11   slant_up   54.26588   47.83150 16.76885 26.93861
12 slant_down   54.26785   47.83590 16.76676 26.93610
13 wide_lines   54.26692   47.83160 16.77000 26.93790

anscombe同様に、平均と分散が小数点以下2位まで同じとなっている。

このように、データから得られる統計量が同一であっても内容が同じ傾向を示すとは限らない。そのため統計量を求める場合には必ずデータをプロットすることが重要である。平均や分散ならヒストグラムを、相関なら散布図を描くことでデータの分布について大まかに把握することができるだろう。

データの可視化というものは特に分析の後半に入ってくると忘れられがちであるが、データの特性を教えてくれる重要な手法である。分析を実施する時にはまず可視化を行うことを大切にしよう。

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

Marketing Mix Modelingにおける広告の累積効果の推定について、以下の記事を書いた。

ushi-goroshi.hatenablog.com
ushi-goroshi.hatenablog.com

本日はこの続きで、同じような条件で作成したデータに対してstanを用いたベイズ推定を実施したので、その内容を紹介する。
なお先日の記事でoptimによりAd-Stock効果の回帰係数および減衰率を十分に精度よく推定できそうなことはわかっているが、

  • 将来的なモデルの高度化
  • 多様な推定方法の確立(手持ちの武器を増やす)

ことを目的としている。

シミュレーションデータの作成 〜 stanによるフィッティングの実行

まずは事前準備として、必要なライブラリを読み込む。当然ながらstanを使うので、事前にインストールを済ませた上で{rstan}を読み込む。

## load libraries
library(dplyr)
library(tidyr)
library(tibble)

## load rstan
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

続いてシミュレーションデータ作成用の関数を用意する。これは先日の記事で紹介したものと大きな変更はないが、拡張性を考慮して変数の順序を入れ替え、また説明変数の作成に用いていたrnormの代わりにrgammaを使用している。

set.seed(123)
simulate_y <- function(pars) {
   n         <- pars[1] # num of observation
   mu        <- pars[2] # intercept
   var_e     <- pars[3] # residual variance
   beta_01   <- pars[4] # regression coefficient of X1 to be esitmated
   lambda_01 <- pars[5] # decay rate of Ad-Stock effect of X1
   beta_02   <- pars[6] # regression coefficient of X2 to be esitmated
   lambda_02 <- pars[7] # decay rate of Ad-Stock effect of X2
   
   X_01_raw <- rgamma(n, 3) * ifelse(runif(n) > 0.7, 0, 1)
   X_01_fil <- stats::filter(X_01_raw, lambda_01, "recursive")
   
   X_02_raw <- rgamma(n, 2) * ifelse(runif(n) > 0.8, 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,
      "X_01_Fil"   = X_01_fil,
      "X_02_Fil"   = X_02_fil,
      "Y_lag"      = dplyr::lag(y, 1),
      "True_Error" = error)
   return(dat)
}

init_par <- array(c(100, 5, 2, 0.5, 0.6, 0.8, 0.5))
dat_Ana  <- na.omit(simulate_y(init_par))

シミュレーションデータが作成できたので、stanを実行するための準備を実行する。stanに渡すためのパラメータを定義し、事前に準備したstanのスクリプトを用いてMCMCによるフィッティングを行う。
サンプルは10,000とし、そのうち最初の5,000サンプルをwarm-upとして破棄した。この辺りの設定は結果を見ながらの試行錯誤によるもので、明確な根拠に基づくものではない。traceプロットやヒストグラムを見る限り、もう少しサンプルを減らしても問題ないように思う。

dat_Ord <- list(N    = nrow(dat_Ana),
                Y    = dat_Ana$Y,
                X_01 = dat_Ana$X_01,
                X_02 = dat_Ana$X_02)

## fitting 
fit_01 <- stan(file = './AdStock_Simulation.stan', 
               data = dat_Ord, 
               iter = 10000,
               chains = 4,
               seed = 1234,
               warmup = 5000)

stanコードの記述

stan()で参照しているstanのコードは以下のように記述した。dataブロックおよびparameterブロックは大した工夫もないが、減衰率を示すlambda_X_01およびlambda_X_02はパラメータスペースを(0, 1)に制限している。

このコードで特筆すべきはtransformed parametersブロックで、ここでlambda_X_0xの値に基づきAd-Stock変数を作成している。Rであればstats::filterで簡単に作成できるのだが、今回の分析では減衰率も推定の対象であるためどうしてもstan内で記述する必要があった。ここで作成したX_01_filおよびX_02_filをmodelブロックで用いている。

data {
   int N;
   vector[N] Y;
   vector[N] X_01;
   vector[N] X_02;
}

parameters {
   real mu;
   real<lower=0> var_error;
   real beta_X_01;
   real<lower = 0, upper = 1> lambda_X_01;
   real beta_X_02;
   real<lower = 0, upper = 1> lambda_X_02;
}

transformed parameters {
   vector[N] X_01_fil;
   vector[N] X_02_fil;

   X_01_fil[1] = X_01[1];
   X_02_fil[1] = X_02[1];
   
   ## Create Ad-stock variable
   for(j in 2:N) {
      X_01_fil[j] = X_01[j] + X_01_fil[j-1] * lambda_X_01;
      X_02_fil[j] = X_02[j] + X_02_fil[j-1] * lambda_X_02;
   }
}

model {
   ## Sampling
   for(i in 1:(N)) {
      Y[i] ~ normal(mu + beta_X_01 * X_01_fil[i] + beta_X_02 * X_02_fil[i], 
      var_error);
   }
}

結果

stanを実行した結果は以下のように取得できる。

## Plot
pars <- c("mu", "var_error", "beta_X_01", "lambda_X_01", "beta_X_02",
          "lambda_X_02")
print(fit_01)
stan_trace(fit_01, pars = pars)
stan_hist(fit_01, pars = pars)
stan_dens(fit_01, separate_chains = T, pars = pars)
stan_ac(fit_01, separate_chains = T, pars = pars)

結果を見てみよう。まずは推定されたパラメータの事後分布から。

> print(fit_01)
Inference for Stan model: AdStock_Simulation.
4 chains, each with iter=10000; warmup=5000; thin=1; 
post-warmup draws per chain=5000, total post-warmup draws=20000.

               mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
mu             4.70    0.01 0.63   3.41   4.29   4.72   5.13   5.91 13079    1
var_error      1.45    0.00 0.11   1.26   1.37   1.44   1.52   1.68 17554    1
beta_X_01      0.50    0.00 0.08   0.35   0.44   0.50   0.55   0.65 15821    1
lambda_X_01    0.63    0.00 0.08   0.44   0.58   0.63   0.69   0.77 12174    1
beta_X_02      0.71    0.00 0.10   0.51   0.64   0.70   0.77   0.91 13692    1
lambda_X_02    0.56    0.00 0.07   0.41   0.52   0.57   0.61   0.69 14090    1
X_01_fil[1]    4.74    0.00 0.00   4.74   4.74   4.74   4.74   4.74     2    1
X_01_fil[2]    3.51    0.00 0.40   2.63   3.27   3.55   3.79   4.18 12174    1
X_01_fil[3]    4.94    0.00 0.52   3.87   4.59   4.96   5.30   5.91 12453    1
X_01_fil[4]    9.09    0.01 0.73   7.65   8.60   9.09   9.58  10.49 12477    1
…
省略
…

X_01_filは本来推定すべきパラメータではないが、stanでは値が決まっていない変数(確率変数と見なされるもの)は全てパラメータとして扱われるため、特に分布に関心はないものの結果として含まれている。

では実際に関心のあるパラメータに着目すると、muが設定した値からやや低く推定されているものの、50%ベイズ信頼区間には真の値が含まれている。またX_01の回帰係数および減衰率も同様に設定した値に近いものを得られている。
一方で残差分散およびX_02の回帰係数・減衰率については、大きく離れてはいないものの50%信頼区間からは外れてしまっている。特に残差分散は何度か試行を繰り返したものの、そのいずれにおいても過少に推定されてしまった。この原因についてはよくわからない。
なおX_02については、回帰係数および減衰率を高めに設定することでより収束しやすい傾向があった。

プロットにより作成されるグラフは、それぞれ以下の通り。

> stan_trace(fit_01, pars = pars)

f:id:ushi-goroshi:20170627001030p:plain

> stan_hist(fit_01, pars = pars)

f:id:ushi-goroshi:20170627001044p:plain

> stan_dens(fit_01, separate_chains = T, pars = pars)

f:id:ushi-goroshi:20170627001116p:plain

> stan_ac(fit_01, separate_chains = T, pars = pars)

f:id:ushi-goroshi:20170627001127p:plain

traceプロットを見ると安定した形となっており、Rhatの数値と合わせて考えても収束していると判断して良さそうだ。分布の形をchainごとで比較しても一部で歪んでいる様子は見当たらない。
なおヒストグラムを確認すると基本的には正規分布に近い形となっているが、減衰率については左に裾を引く形の分布となっている。

まとめ

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

  1. stanを用いることでAd-Stock効果をそれなりに精度よく推定できる
  2. 減衰率のような非線形の効果を素直にモデルに取り込むことができ、また分布の形もKoyckに限定されることなく柔軟なモデリングが可能
  3. しかしパルダ・モデルやoptimによる推定と比較すると精度に劣り、また処理時間およびコーディングに要する時間も比較的多くなる

と言えそうだ。分析準備に要するコストは習熟度が上がれば削減されるだろうが、一回の分析に要する処理時間が相対的に大きいため試行錯誤する場合には向いていないと思われる。
一方でベイズモデリングの大きなメリットである「モデルの柔軟性」は非常に魅力的に感じ、Marketing Mix Modelingの高度化を目指すに当たって積極的に使っていきたい技術である。

上記を踏まえてこれまでの分析結果をまとめると、

  • Ad-Stockの分布としてKoyck型を仮定でき、残存の規模(減衰率)に関心はない場合 → パルダ・モデル
  • 残存の規模や分布に関心があり、試行錯誤が必要 → optimによる推定
  • より柔軟なモデリング(階層ベイズや状態空間モデルへの拡張)を実施したい → stanによるモデリング

といった使い分けが良さそうだ。なおstanによるモデリングは分析者のアイディアおよび対象領域に対する知見の深さにかなり依存すると思われるため、少なくともMMM構築の初期では全く適していないだろう。
逆に言えば、当該領域について十分な経験を持った分析者がいるのであれば、アイディアの限りを用いて色々なモデルを試すことが出来るため大きな価値を提供することが出来るようになるだろう。

最後に、この分析に限らずMCMCによるもの全般に言えることだが、やはりパラメータの分布という形で結果を見せることができるというのは大きな強みであると感じた。統計解析に馴染みのない相手であっても、「パラメータがバラつきを持ち、その範囲を分布として示す」というのは直感的にわかりやすいだろうと思う。頻度論的なアプローチであってもパラメータの(誤差)分布を標準誤差として示すことは出来るが、理解の容易さという点では比較にならないだろう。

唯一の懸念は、「パラメータが分布を持つ」という仮定に耐えうるか、ということだろうか。

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の指定には"&"を使ったほうが良さそうだ。