## ---- echo = FALSE------------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(BTYDplus)
data("groceryElog")

## ---- fig.show="hold", fig.width=7, fig.height=2.8, fig.cap="Timing patterns for sampled grocery customers"----
library(BTYDplus)
data("groceryElog")
set.seed(123)
# plot timing patterns of 30 sampled customers
plotTimingPatterns(groceryElog, n = 30, T.cal = "2007-05-15",
                   headers = c("Past", "Future"), title = "")

## ---- echo=FALSE, results="asis"----------------------------------------------
cdnowElog <- read.csv(system.file("data/cdnowElog.csv", package = "BTYD"),
                      stringsAsFactors = FALSE,
                      col.names = c("cust", "sampleid", "date", "cds", "sales"))
cdnowElog$date <- as.Date(as.character(cdnowElog$date), format = "%Y%m%d")
knitr::kable(head(cdnowElog[, c("cust", "date", "sales")], 6), caption = "Transaction Log Example")

## ---- fig.show="hold", fig.width=7, fig.height=2.5, fig.cap="Weekly trends for the grocery dataset"----
data("groceryElog")
op <- par(mfrow = c(1, 2), mar = c(2.5, 2.5, 2.5, 2.5))
# incremental
weekly_inc_total  <- elog2inc(groceryElog, by = 7, first = TRUE)
weekly_inc_repeat <- elog2inc(groceryElog, by = 7, first = FALSE)
plot(weekly_inc_total, typ = "l", frame = FALSE, main = "Incremental")
lines(weekly_inc_repeat, col = "red")
# cumulative
weekly_cum_total  <- elog2cum(groceryElog, by = 7, first = TRUE)
weekly_cum_repeat <- elog2cum(groceryElog, by = 7, first = FALSE)
plot(weekly_cum_total, typ = "l", frame = FALSE, main = "Cumulative")
lines(weekly_cum_repeat, col = "red")
par(op)

## -----------------------------------------------------------------------------
data("groceryElog")
head(elog2cbs(groceryElog), 5)

## -----------------------------------------------------------------------------
data("groceryElog")
range(groceryElog$date)
groceryCBS <- elog2cbs(groceryElog, T.cal = "2006-12-31")
head(groceryCBS, 5)

## ---- fig.show="hold", fig.width=7, fig.height=2.2, fig.cap="Diagnostic plots for estimating regularity"----
data("groceryElog")
op <- par(mfrow = c(1, 2))
(k.wheat <- estimateRegularity(groceryElog, method = "wheat", 
                              plot = TRUE, title = "Wheat & Morrison"))
(k.mle <- estimateRegularity(groceryElog, method = "mle", 
                             plot = TRUE, title = "Maximum Likelihood"))
par(op)

## -----------------------------------------------------------------------------
# load grocery dataset, if it hasn't been done before
if (!exists("groceryCBS")) {
  data("groceryElog")
  groceryCBS <- elog2cbs(groceryElog, T.cal = "2006-12-31")
}
# estimate NBD parameters
round(params.nbd <- nbd.EstimateParameters(groceryCBS), 3)
# report log-likelihood
nbd.cbs.LL(params.nbd, groceryCBS)

## -----------------------------------------------------------------------------
# calculate expected future transactions for customers who've 
# had 1 to 5 transactions in first 52 weeks
est5.nbd <- nbd.ConditionalExpectedTransactions(params.nbd, 
              T.star = 52, x = 1:5, T.cal = 52)
for (i in 1:5) {
  cat("x =", i, ":", sprintf("%5.3f", est5.nbd[i]), "\n")
}

## -----------------------------------------------------------------------------
# predict whole customer cohort
groceryCBS$xstar.nbd <- nbd.ConditionalExpectedTransactions(
                    params = params.nbd, T.star = 52, 
                    x = groceryCBS$x, T.cal = groceryCBS$T.cal)
# compare predictions with actuals at aggregated level
rbind(`Actuals` = c(`Holdout` = sum(groceryCBS$x.star)), 
      `NBD`     = c(`Holdout` = round(sum(groceryCBS$xstar.nbd))))

## -----------------------------------------------------------------------------
# load grocery dataset, if it hasn't been done before
if (!exists("groceryCBS")) {
  data("groceryElog")
  groceryCBS <- elog2cbs(groceryElog, T.cal = "2006-12-31")
}
# estimate Pareto/NBD parameters
params.pnbd <- BTYD::pnbd.EstimateParameters(groceryCBS[, c("x", "t.x", "T.cal")])
names(params.pnbd) <- c("r", "alpha", "s", "beta")
round(params.pnbd, 3)
# report log-likelihood
BTYD::pnbd.cbs.LL(params.pnbd, groceryCBS[, c("x", "t.x", "T.cal")])

## -----------------------------------------------------------------------------
# calculate expected future transactions for customers who've 
# had 1 to 5 transactions in first 12 weeks, but then remained
# inactive for 40 weeks
est5.pnbd <- BTYD::pnbd.ConditionalExpectedTransactions(params.pnbd, 
               T.star = 52, x = 1:5, t.x = 12, T.cal = 52)
for (i in 1:5) {
  cat("x =", i, ":", sprintf("%5.3f", est5.pnbd[i]), "\n")
}

## -----------------------------------------------------------------------------
# predict whole customer cohort
groceryCBS$xstar.pnbd <- BTYD::pnbd.ConditionalExpectedTransactions(
                     params = params.pnbd, T.star = 52, 
                     x = groceryCBS$x, t.x = groceryCBS$t.x,
                     T.cal = groceryCBS$T.cal)
# compare predictions with actuals at aggregated level
rbind(`Actuals`    = c(`Holdout` = sum(groceryCBS$x.star)), 
      `Pareto/NBD` = c(`Holdout` = round(sum(groceryCBS$xstar.pnbd))))

## -----------------------------------------------------------------------------
# P(alive) for customers who've had 1 to 5 transactions in first 
# 12 weeks, but then remained inactive for 40 weeks
palive.pnbd <- BTYD::pnbd.PAlive(params.pnbd, 
                 x = 1:5, t.x = 12, T.cal = 52)
for (i in 1:5) {
  cat("x =", i, ":", sprintf("%5.2f %%", 100*palive.pnbd[i]), "\n")
}

## -----------------------------------------------------------------------------
# load grocery dataset, if it hasn't been done before
if (!exists("groceryCBS")) {
  data("groceryElog")
  groceryCBS <- elog2cbs(groceryElog, T.cal = "2006-12-31")
}
# estimate parameters for various models
params.bgnbd   <- BTYD::bgnbd.EstimateParameters(groceryCBS) # BG/NBD
params.bgcnbd  <- bgcnbd.EstimateParameters(groceryCBS)      # BG/CNBD-k
params.mbgnbd  <- mbgnbd.EstimateParameters(groceryCBS)      # MBG/NBD
params.mbgcnbd <- mbgcnbd.EstimateParameters(groceryCBS)     # MBG/CNBD-k
row <- function(params, LL) { 
  names(params) <- c("k", "r", "alpha", "a", "b")
  c(round(params, 3), LL = round(LL))
}
rbind(`BG/NBD`     = row(c(1, params.bgnbd), 
                         BTYD::bgnbd.cbs.LL(params.bgnbd, groceryCBS)),
      `BG/CNBD-k`  = row(params.bgcnbd, 
                         bgcnbd.cbs.LL(params.bgcnbd, groceryCBS)),
      `MBG/NBD`    = row(params.mbgnbd, 
                         mbgcnbd.cbs.LL(params.mbgnbd, groceryCBS)),
      `MBG/CNBD-k` = row(params.mbgcnbd, 
                         mbgcnbd.cbs.LL(params.mbgcnbd, groceryCBS)))

## -----------------------------------------------------------------------------
# calculate expected future transactions for customers who've 
# had 1 to 5 transactions in first 12 weeks, but then remained
# inactive for 40 weeks
est5.mbgcnbd <- mbgcnbd.ConditionalExpectedTransactions(params.mbgcnbd,
                  T.star = 52, x = 1:5, t.x = 12, T.cal = 52)
for (i in 1:5) {
  cat("x =", i, ":", sprintf("%5.3f", est5.mbgcnbd[i]), "\n")
}

## -----------------------------------------------------------------------------
# P(alive) for customers who've had 1 to 5 transactions in first 
# 12 weeks, but then remained inactive for 40 weeks
palive.mbgcnbd <- mbgcnbd.PAlive(params.mbgcnbd, 
                    x = 1:5, t.x = 12, T.cal = 52)
for (i in 1:5) {
  cat("x =", i, ":", sprintf("%5.2f %%", 100*palive.mbgcnbd[i]), "\n")
}

## -----------------------------------------------------------------------------
# predict whole customer cohort
groceryCBS$xstar.mbgcnbd <- mbgcnbd.ConditionalExpectedTransactions(
                        params = params.mbgcnbd, T.star = 52,
                        x = groceryCBS$x, t.x = groceryCBS$t.x, 
                        T.cal = groceryCBS$T.cal)
# compare predictions with actuals at aggregated level
rbind(`Actuals`    = c(`Holdout` = sum(groceryCBS$x.star)),
      `MBG/CNBD-k` = c(`Holdout` = round(sum(groceryCBS$xstar.mbgcnbd))))

## ---- fig.show="hold", fig.width=4, fig.height=3.5, fig.cap="Weekly actuals vs. MBG/CNBD-k predictions"----
# runs for ~37secs on a 2015 MacBook Pro
nil <- mbgcnbd.PlotTrackingInc(params.mbgcnbd,
         T.cal = groceryCBS$T.cal,
         T.tot = max(groceryCBS$T.cal + groceryCBS$T.star),
         actual.inc.tracking = elog2inc(groceryElog))

## -----------------------------------------------------------------------------
# mean absolute error (MAE)
mae <- function(act, est) {
  stopifnot(length(act)==length(est))
  sum(abs(act-est)) / length(act)
}
mae.nbd <- mae(groceryCBS$x.star, groceryCBS$xstar.nbd)
mae.pnbd <- mae(groceryCBS$x.star, groceryCBS$xstar.pnbd)
mae.mbgcnbd <- mae(groceryCBS$x.star, groceryCBS$xstar.mbgcnbd)
rbind(`NBD`        = c(`MAE` = round(mae.nbd, 3)),
      `Pareto/NBD` = c(`MAE` = round(mae.pnbd, 3)),
      `MBG/CNBD-k` = c(`MAE` = round(mae.mbgcnbd, 3)))
lift <- 1 - mae.mbgcnbd / mae.pnbd
cat("Lift in MAE for MBG/CNBD-k vs. Pareto/NBD:", round(100*lift, 1), "%")

## -----------------------------------------------------------------------------
# load grocery dataset, if it hasn't been done before
if (!exists("groceryCBS")) {
  data("groceryElog")
  groceryCBS <- elog2cbs(groceryElog, T.cal = "2006-12-31")
}
# generate parameter draws (~13secs on 2015 MacBook Pro)
pnbd.draws <- pnbd.mcmc.DrawParameters(groceryCBS)
# generate draws for holdout period
pnbd.xstar.draws <- mcmc.DrawFutureTransactions(groceryCBS, pnbd.draws)
# conditional expectations
groceryCBS$xstar.pnbd.hb   <- apply(pnbd.xstar.draws, 2, mean)
# P(active)
groceryCBS$pactive.pnbd.hb <- mcmc.PActive(pnbd.xstar.draws)
# P(alive)
groceryCBS$palive.pnbd.hb  <- mcmc.PAlive(pnbd.draws)
# show estimates for first few customers
head(groceryCBS[, c("x", "t.x", "x.star",
                    "xstar.pnbd.hb", "pactive.pnbd.hb",
                    "palive.pnbd.hb")])

## -----------------------------------------------------------------------------
class(pnbd.draws$level_2)

# convert cohort-level draws from coda::mcmc.list to a matrix, with 
# each parameter becoming a column, and each draw a row
cohort.draws <- pnbd.draws$level_2
head(as.matrix(cohort.draws), 5)

# compute median across draws, and compare to ML estimates; as can be
# seen, the two parameter estimation approaches result in very similar
# estimates
round(
  rbind(`Pareto/NBD (HB)` = apply(as.matrix(cohort.draws), 2, median),
        `Pareto/NBD`      = BTYD::pnbd.EstimateParameters(groceryCBS[, c("x", "t.x", "T.cal")]))
, 2)

## ---- fig.show="hold", warning=FALSE, fig.width=7, fig.height=3, fig.cap="MCMC traces and parameter distributions of cohort-level parameters"----
# plot trace- and density-plots for heterogeneity parameters
op <- par(mfrow = c(2, 4), mar = c(2.5, 2.5, 2.5, 2.5))
coda::traceplot(pnbd.draws$level_2)
coda::densplot(pnbd.draws$level_2)
par(op)

## ---- fig.show="hold", warning=FALSE, fig.width=7, fig.height=3, fig.cap="MCMC traces and parameter distributions of individual-level parameters for a specific customer"----
class(pnbd.draws$level_1)
length(pnbd.draws$level_1)

customer4 <- "4"
customer4.draws <- pnbd.draws$level_1[[customer4]]
head(as.matrix(customer4.draws), 5)
round(apply(as.matrix(customer4.draws), 2, median), 3)

# plot trace- and density-plots for customer4 parameters
op <- par(mfrow = c(2, 4), mar = c(2.5, 2.5, 2.5, 2.5))
coda::traceplot(pnbd.draws$level_1[[customer4]])
coda::densplot(pnbd.draws$level_1[[customer4]])
par(op)

## ---- eval = FALSE------------------------------------------------------------
#  # runs for ~120secs on a MacBook Pro 2015
#  op <- par(mfrow = c(1, 2))
#  nil <- mcmc.PlotFrequencyInCalibration(pnbd.draws, groceryCBS)
#  nil <- mcmc.PlotTrackingInc(pnbd.draws,
#           T.cal = groceryCBS$T.cal,
#           T.tot = max(groceryCBS$T.cal + groceryCBS$T.star),
#           actual.inc.tracking.data = elog2inc(groceryElog))
#  par(op)

## ---- eval = FALSE------------------------------------------------------------
#  # load CDNow event log from BTYD package
#  cdnowElog <- read.csv(
#                 system.file("data/cdnowElog.csv", package = "BTYD"),
#                 stringsAsFactors = FALSE,
#                 col.names = c("cust", "sampleid", "date", "cds", "sales"))
#  cdnowElog$date <- as.Date(as.character(cdnowElog$date),
#                            format = "%Y%m%d")
#  # convert to CBS; split into 39 weeks calibration, and 39 weeks holdout
#  cdnowCbs <- elog2cbs(cdnowElog,
#                       T.cal = "1997-09-30", T.tot = "1998-06-30")
#  
#  # estimate Pareto/NBD (Abe) without covariates; model M1 in Abe (2009)
#  draws.m1 <- abe.mcmc.DrawParameters(cdnowCbs,
#                mcmc = 7500, burnin = 2500) # ~33secs on 2015 MacBook Pro
#  quant <- function(x) round(quantile(x, c(0.025, 0.5, 0.975)), 2)
#  t(apply(as.matrix(draws.m1$level_2), 2, quant))
#  #>                        2.5%   50% 97.5%
#  #> log_lambda            -3.70 -3.54 -3.32
#  #> log_mu                -3.96 -3.59 -3.26
#  #> var_log_lambda         1.10  1.34  1.65
#  #> cov_log_lambda_log_mu -0.20  0.13  0.74
#  #> var_log_mu             1.44  2.62  5.05
#  
#  #' append dollar amount of first purchase to use as covariate
#  first <- aggregate(sales ~ cust, cdnowElog, function(x) x[1] * 10^-3)
#  names(first) <- c("cust", "first.sales")
#  cdnowCbs <- merge(cdnowCbs, first, by = "cust")
#  
#  #' estimate with first purchase spend as covariate; model M2 in Abe (2009)
#  draws.m2 <- abe.mcmc.DrawParameters(cdnowCbs,
#                covariates = c("first.sales"),
#                mcmc = 7500, burnin = 2500) # ~33secs on 2015 MacBook Pro
#  t(apply(as.matrix(draws.m2$level_2), 2, quant))
#  #>                         2.5%   50% 97.5%
#  #> log_lambda_intercept   -4.02 -3.77 -3.19
#  #> log_mu_intercept       -4.37 -3.73 -2.69
#  #> log_lambda_first.sales  0.04  6.04  9.39
#  #> log_mu_first.sales     -9.02  1.73  7.90
#  #> var_log_lambda          0.01  1.35  1.79
#  #> cov_log_lambda_log_mu  -0.35  0.22  0.76
#  #> var_log_mu              0.55  2.59  4.97

## ---- eval=FALSE--------------------------------------------------------------
#  # load grocery dataset, if it hasn't been done before
#  if (!exists("groceryCBS")) {
#    data("groceryElog")
#    groceryCBS <- elog2cbs(groceryElog, T.cal = "2006-12-31")
#  }
#  # estimte Pareto/GGG
#  pggg.draws <- pggg.mcmc.DrawParameters(groceryCBS) # ~2mins on 2015 MacBook Pro
#  # generate draws for holdout period
#  pggg.xstar.draws <- mcmc.DrawFutureTransactions(groceryCBS, pggg.draws)
#  # conditional expectations
#  groceryCBS$xstar.pggg  <- apply(pggg.xstar.draws, 2, mean)
#  # P(active)
#  groceryCBS$pactive.pggg <- mcmc.PActive(pggg.xstar.draws)
#  # P(alive)
#  groceryCBS$palive.pggg  <- mcmc.PAlive(pggg.draws)
#  # show estimates for first few customers
#  head(groceryCBS[, c("x", "t.x", "x.star",
#                      "xstar.pggg", "pactive.pggg", "palive.pggg")])
#  #>    x      t.x x.star xstar.pggg pactive.pggg palive.pggg
#  #> 1  0  0.00000      0       0.02         0.02        0.03
#  #> 2  1 50.28571      0       1.01         0.59        1.00
#  #> 3 19 48.57143     14      14.76         0.87        0.87
#  #> 4  0  0.00000      0       0.04         0.03        0.13
#  #> 5  2 40.42857      3       2.02         0.84        0.91
#  #> 6  5 47.57143      6       4.46         0.92        0.95
#  
#  # report median cohort-level parameter estimates
#  round(apply(as.matrix(pggg.draws$level_2), 2, median), 3)
#  #>     t gamma     r alpha     s  beta
#  #> 1.695 0.373 0.948 5.243 0.432 4.348
#  # report mean over median individual-level parameter estimates
#  median.est <- sapply(pggg.draws$level_1, function(draw) {
#    apply(as.matrix(draw), 2, median)
#  })
#  round(apply(median.est, 1, mean), 3)
#  #>    k lambda     mu    tau      z
#  #> 3.892  0.160  0.065 69.546  0.316

## ---- eval = FALSE------------------------------------------------------------
#  # compare predictions with actuals at aggregated level
#  rbind(`Actuals`    = c(`Holdout` = sum(groceryCBS$x.star)),
#        `Pareto/GGG` = round(sum(groceryCBS$xstar.pggg)),
#        `MBG/CNBD-k` = round(sum(groceryCBS$xstar.mbgcnbd)),
#        `Pareto/NBD (HB)` = round(sum(groceryCBS$xstar.pnbd.hb)))
#  #>                 Holdout
#  #> Actuals            3389
#  #> Pareto/GGG         3815
#  #> MBG/CNBD-k         3970
#  #> Pareto/NBD (HB)    4018
#  
#  # error on customer level
#  mae <- function(act, est) {
#    stopifnot(length(act)==length(est))
#    sum(abs(act-est)) / length(act)
#  }
#  mae.pggg <- mae(groceryCBS$x.star, groceryCBS$xstar.pggg)
#  mae.mbgcnbd <- mae(groceryCBS$x.star, groceryCBS$xstar.mbgcnbd)
#  mae.pnbd.hb <- mae(groceryCBS$x.star, groceryCBS$xstar.pnbd.hb)
#  rbind(`Pareto/GGG`      = c(`MAE` = round(mae.pggg, 3)),
#        `MBG/CNBD-k`      = c(`MAE` = round(mae.mbgcnbd, 3)),
#        `Pareto/NBD (HB)` = c(`MAE` = round(mae.pnbd.hb, 3)))
#  #>                   MAE
#  #> Pareto/GGG      0.621
#  #> MBG/CNBD-k      0.644
#  #> Pareto/NBD (HB) 0.688
#  
#  lift <- 1 - mae.pggg / mae.pnbd.hb
#  cat("Lift in MAE:", round(100*lift, 1), "%")
#  #> Lift in MAE for Pareto/GGG vs. Pareto/NBD: 9.8%

