ビッグデータ時代のマーケティングの事例をRでやってみた（１）

ビッグデータ時代のマーケティングは名著なのですが、事例に関するコードが公開されていないのが難点です。 逆に一つずつ自分でコードを書くと勉強になりそうなのでひとつやってみました。

２．４．「マーケティングにおけるモデリング例」です。

まず、P46の「データ生成の手順」に従ってデータフレームを作ります。

library(dplyr)
#データ生成
err <- rnorm(100, mean=0, sd=0.3)
Z1 <- runif(100, min=0.5, max=1)
u <- runif(100)
Z2 <- (u <= 0.3)*(Z1 <= 0.8)
y <- exp(0.5 - 2.5 * log(Z1) + 0.8 * Z2 + err)  %>% floor()
data <- data.frame (y, Z1, Z2)

散布図は次のようになり、だいたい、本と近くなっています。

モデルが10個あるのでリストを作って格納します。
#モデルリスト
mf <- list(y~1, y~Z1, y~Z1+Z2, y~log(Z1), y~log(Z1)+Z2,
           log(y)~1, log(y)~Z1, log(y)~Z1+Z2, log(y)~log(Z1), log(y)~log(Z1)+Z2)

AIC、ヤコビアンを計算します。

#AIC格納用リスト
AICb <- vector("double", 10)
AICa <- vector("double", 10)
#補正前AIC
for(i in 1:10){
  lm(mf[[i]], data=data) %>%
    AIC() -> AICb[i]
}
#ヤコビアン
jacob <- sum(log(numDeriv::grad(log, x=data$y)))
#補正後AIC
AICa <- AICb + c(rep(0, 5), rep(-2 * jacob, 5))

data.frame(AICb, AICa) %>%
  set_rownames(paste0("model", 1:10)) %>%
  knitr::kable(format="markdown")

	AICb	AICa
model1	621.08099	621.0810
model2	532.42325	532.4232
model3	502.95307	502.9531
model4	525.10030	525.1003
model5	493.74816	493.7482
model6	264.85676	553.5136
model7	136.97222	425.6290
model8	90.41715	379.0740
model9	139.38691	428.0437
model10	89.20195	377.8588

確かにmodel1～5では、1>2>4>3>5

model6～10では、6>9>7>8>10

model5と10を比べると5>10となり、model10が最適という結果が得られました。

最後に現況再現性の確認をします。

model[[10]] %>% predict() %>% plot(type="l", ylim=c(-1, 3.5))
par(new=T) 
log(y) %>% plot(type="l", lty="dashed", ylim=c(-1, 3.5))
par(new=T) 
(predict(model[[10]]) - log(y)) %>%
  plot(type="l", lwd=2, ylim=c(-1, 3.5))