読者です 読者をやめる 読者になる 読者になる

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

R 統計

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

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

まず、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)

散布図は次のようになり、だいたい、本と近くなっています。 f:id:datascienceconsultant:20161227102812p:plain

モデルが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))

f:id:datascienceconsultant:20161227104203p:plain