construct_cdnow_twohier_pnbd_models.qmd

---
title: "Construct Two Parameter Hierarchical P/NBD Model for CD-Now Transaction Data"
author: "Mick Cooney <mickcooney@gmail.com>"
date: "Last updated: `r format(Sys.time(), '%B %d, %Y')`"
editor: source
execute:
  message: false
  warning: false
  error: false
format:
  html:
    light: superhero
    dark: darkly
    anchor-sections: true
    embed-resources: true
    number-sections: true
    smooth-scroll: true
    toc: true
    toc-depth: 3
    toc-location: left
    code-fold: true
    code-summary: "Show code"
---


```{r import_libraries}
#| echo: FALSE
#| message: FALSE

library(conflicted)
library(tidyverse)
library(scales)
library(cowplot)
library(directlabels)
library(magrittr)
library(rlang)
library(fs)
library(purrr)
library(furrr)
library(glue)
library(cmdstanr)
library(brms)
library(posterior)
library(bayesplot)
library(tidybayes)


source("lib_utils.R")
source("lib_btyd.R")


conflict_lst <- resolve_conflicts(
  c("magrittr", "rlang", "dplyr", "readr", "purrr", "ggplot2", "MASS",
    "fitdistrplus")
  )

options(
  width = 80L,
  warn  = 1,
  mc.cores = parallel::detectCores()
  )


set.seed(42)
stanfit_seed <- 4203
n_sim        <- 2000

theme_set(theme_cowplot())
plan(multisession)
```

In this workbook we construct the non-hierarchical P/NBD models on the
synthetic data with the longer timeframe.




# Load and Construct Datasets


## Load CD-NOW Transaction Data

We now want to load the CD-NOW transaction data.


```{r load_cdnow_transaction_data}
#| echo: TRUE

customer_cohortdata_tbl <- read_rds("data/cdnow_customer_cohort_data_tbl.rds")
customer_cohortdata_tbl |> glimpse()

customer_transactions_tbl <- read_rds("data/cdnow_transaction_data_tbl.rds")
customer_transactions_tbl |> glimpse()

customer_subset_id <- read_rds("data/cdnow_customer_subset_ids.rds")
customer_subset_id |> glimpse()
```


We re-produce the visualisation of the transaction times we used in previous
workbooks.

```{r plot_customer_transaction_times}
#| echo: TRUE

plot_tbl <- customer_transactions_tbl |>
  group_nest(customer_id, .key = "cust_data") |>
  filter(map_int(cust_data, nrow) > 3) |>
  slice_sample(n = 30) |>
  unnest(cust_data)

ggplot(plot_tbl, aes(x = tnx_timestamp, y = customer_id)) +
  geom_line() +
  geom_point() +
  labs(
      x = "Date",
      y = "Customer ID",
      title = "Visualisation of Customer Transaction Times"
    ) +
  theme(axis.text.y = element_text(size = 10))
```





## Load Derived Data

```{r load_derived_data}
#| echo: TRUE

obs_fitdata_tbl   <- read_rds("data/cdnow_obs_fitdata_tbl.rds")
obs_validdata_tbl <- read_rds("data/cdnow_obs_validdata_tbl.rds")

customer_fit_stats_tbl <- obs_fitdata_tbl |>
  rename(x = tnx_count)
```


## Load Subset Data

We also want to construct our data subsets for the purposes of speeding up our
valuations.

```{r construct_customer_subset_data}
#| echo: TRUE

customer_fit_subset_tbl <- obs_fitdata_tbl |>
  filter(customer_id %in% customer_subset_id)

customer_fit_subset_tbl |> glimpse()


customer_valid_subset_tbl <- obs_validdata_tbl |>
  filter(customer_id %in% customer_subset_id)

customer_valid_subset_tbl |> glimpse()
```


We now use these datasets to set the start and end dates for our various
validation methods.


```{r set_start_end_dates}
dates_lst <- read_rds("data/cdnow_simulation_dates.rds")

use_fit_start_date <- dates_lst$use_fit_start_date
use_fit_end_date   <- dates_lst$use_fit_end_date

use_valid_start_date <- dates_lst$use_valid_start_date
use_valid_end_date   <- dates_lst$use_valid_end_date
```

We now split out the transaction data into fit and validation datasets.

```{r create_customer_transaction_splits}
#| echo: true

customer_fit_transactions_tbl <- customer_transactions_tbl |>
  filter(
    customer_id %in% customer_subset_id,
    tnx_timestamp >= use_fit_start_date,
    tnx_timestamp <= use_fit_end_date
    )
  
customer_fit_transactions_tbl |> glimpse()


customer_valid_transactions_tbl <- customer_transactions_tbl |>
  filter(
    customer_id %in% customer_subset_id,
    tnx_timestamp >= use_valid_start_date,
    tnx_timestamp <= use_valid_end_date
    )
  
customer_valid_transactions_tbl |> glimpse()
```

Finally, we want to extract the first transaction for each customer, so we
can add this data to assess our models.

```{r extract_customer_first_transaction}
#| echo: true

customer_initial_tnx_tbl <- customer_fit_transactions_tbl |>
  slice_min(n = 1, order_by = tnx_timestamp, by = customer_id)

customer_initial_tnx_tbl |> glimpse()
```

We now expand out these initial transactions so that we can append them to
our simulations.

```{r expand_initial_simulation_transactions}
#| echo: true

sim_init_tbl <- customer_initial_tnx_tbl |>
  transmute(
    customer_id,
    draw_id       = list(1:n_sim),
    tnx_timestamp,
    tnx_amount    = total_spend
    ) |>
  unnest(draw_id)

sim_init_tbl |> glimpse()
```


Before we start on that, we set a few parameters for the workbook to organise
our Stan code.

```{r setup_workbook_parameters}
#| echo: TRUE

stan_modeldir <- "stan_models"
stan_codedir  <-   "stan_code"
```


# Fit Hierarchical Mean-Lambda, Mean-Mu Model

Our first  model puts a hierarchical prior around the mean of both the
population $\lambda$, `lambda_mn`, and population $\mu$, `mu_mn`.

Once again we use a Gamma prior for both these parameters


## Compile and Fit Stan Model

We now compile this model using `CmdStanR`.

```{r compile_pnbd_twohiermean_stanmodel}
#| echo: TRUE
#| results: "hide"

pnbd_twohiermean_stanmodel <- cmdstan_model(
  "stan_code/pnbd_twohier_mean.stan",
  include_paths =   stan_codedir,
  pedantic      =           TRUE,
  dir           =  stan_modeldir
  )
```


We then use this compiled model with our data to produce a fit of the data.



```{r fit_pnbd_cdnow_twohiermean_stanmodel}
#| echo: TRUE

stan_modelname <- "pnbd_cdnow_twohiermean"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_lambda_mn_p1 = 0.25,
    hier_lambda_mn_p2 = 1.00,
    
    hier_mu_mn_p1     = 0.10,
    hier_mu_mn_p2     = 1.00,

    lambda_cv         = 1.00,
    mu_cv             = 1.00
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_cdnow_twohiermean_stanfit <- pnbd_twohiermean_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_cdnow_twohiermean_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  message(glue("Found file {stanfit_object_file}. Loading..."))
  
  pnbd_cdnow_twohiermean_stanfit <- read_rds(stanfit_object_file)
}

pnbd_cdnow_twohiermean_stanfit$print()
```


We have some basic HMC-based validity statistics we can check.

```{r calculate_pnbd_cdnow_twohiermean_hmc_diagnostics}
#| echo: TRUE

pnbd_cdnow_twohiermean_stanfit$cmdstan_diagnose()
```



## Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a
few different visualisations of the diagnostics.


```{r plot_pnbd_cdnow_twohiermean_traceplots_nowarmup}
#| echo: TRUE

parameter_subset <- c(
  "lambda_mn", "mu_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_cdnow_twohiermean_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))
```


We also check $N_{eff}$ as a quick diagnostic of the fit.


```{r plot_pnbd_cdnow_twohiermean_parameter_neffratio}
#| echo: TRUE

pnbd_cdnow_twohiermean_stanfit |>
  neff_ratio(pars = c("lambda_mn", "mu_mn", "lambda", "mu")) |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")
```


Finally, we want to check out the energy diagnostic, which is often indicative
of problems with the posterior mixing.

```{r plot_pnbd_cdnow_twohiermean_energy}
#| echo: true

pnbd_cdnow_twohiermean_stanfit |>
  nuts_params() |>
  mcmc_nuts_energy(binwidth = 50)
```



## Assess the Model

As we intend to run the same logic to assess each of our models, we have
combined all this logic into a single function `run_model_assessment`, to 
run the simulations and combine the datasets.


```{r run_pnbd_cdnow_twohiermean_assessment}
#| echo: TRUE

pnbd_stanfit <- pnbd_cdnow_twohiermean_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_cdnow_twohiermean_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_cdnow_twohiermean",
  fit_end_dttm     = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm   = use_valid_end_date   |> as.POSIXct(),
  sim_seed         = 2210
  )

pnbd_cdnow_twohiermean_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

We first check the model against the in-sample data.

```{r run_pnbd_cdnow_twohiermean_fit_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_cdnow_twohiermean_assess_data_lst |>
  use_series(model_fit_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  bind_rows(sim_init_tbl) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)
```

This fit looks reasonable and appears to capture most of the aspects of the
data used to fit it. Given that this is a synthetic dataset, this is not
surprising, but at least we appreciate that our model is valid.


### Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

```{r run_pnbd_cdnow_twohiermean_valid_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_cdnow_twohiermean_assess_data_lst |>
  use_series(model_valid_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)
```

As for our short time frame data, overall our model is working well.



# Fit Hierarchical Lambda Model

Our next model puts a hierarchical prior around both $\lambda$ parameters,
`lambda_mn` and `lambda_cv`.

Once again we use a Gamma prior for both these parameters.


## Compile and Fit Stan Model

We now compile this model using `CmdStanR`.

```{r compile_pnbd_twohierlambda_stanmodel}
#| echo: TRUE
#| results: "hide"

pnbd_twohierlambda_stanmodel <- cmdstan_model(
  "stan_code/pnbd_twohier_lambda.stan",
  include_paths =   stan_codedir,
  pedantic      =           TRUE,
  dir           =  stan_modeldir
  )
```


We then use this compiled model with our data to produce a fit of the data.



```{r fit_pnbd_cdnow_twohierlambda_stanmodel}
#| echo: TRUE

stan_modelname <- "pnbd_cdnow_twohierlambda"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_lambda_mn_p1 = 0.25,
    hier_lambda_mn_p2 = 1.00,
    
    hier_lambda_cv_p1 = 1.00,
    hier_lambda_cv_p2 = 1.00,

    mu_mn             = 0.10,
    mu_cv             = 1.00
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_cdnow_twohierlambda_stanfit <- pnbd_twohierlambda_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_cdnow_twohierlambda_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  message(glue("Found file {stanfit_object_file}. Loading..."))
  
  pnbd_cdnow_twohierlambda_stanfit <- read_rds(stanfit_object_file)
}

pnbd_cdnow_twohierlambda_stanfit$print()
```


We have some basic HMC-based validity statistics we can check.

```{r calculate_pnbd_cdnow_twohierlambda_hmc_diagnostics}
#| echo: TRUE

pnbd_cdnow_twohierlambda_stanfit$cmdstan_diagnose()
```



## Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a
few different visualisations of the diagnostics.


```{r plot_pnbd_cdnow_twohierlambda_lambda_traceplots_nowarmup}
#| echo: TRUE

parameter_subset <- c(
  "lambda_mn", "lambda_cv", "r", "alpha",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_cdnow_twohierlambda_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))
```


We also check $N_{eff}$ as a quick diagnostic of the fit.


```{r plot_pnbd_cdnow_twohierlambda_parameter_neffratio}
#| echo: TRUE

pnbd_cdnow_twohierlambda_stanfit |>
  neff_ratio(
    pars = c("lambda_mn", "lambda_cv", "r", "alpha", "lambda", "mu")
    ) |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")
```


Finally, we want to check out the energy diagnostic, which is often indicative
of problems with the posterior mixing.

```{r plot_pnbd_cdnow_twohierlambda_energy}
#| echo: true

pnbd_cdnow_twohierlambda_stanfit |>
  nuts_params() |>
  mcmc_nuts_energy(binwidth = 50)
```



## Assess the Model

As we intend to run the same logic to assess each of our models, we have
combined all this logic into a single function `run_model_assessment`, to 
run the simulations and combine the datasets.


```{r run_pnbd_cdnow_twohierlambda_assessment}
#| echo: TRUE

pnbd_stanfit <- pnbd_cdnow_twohierlambda_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_cdnow_twohierlambda_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_cdnow_twohierlambda",
  fit_end_dttm     = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm   = use_valid_end_date   |> as.POSIXct(),
  sim_seed         = 2220
  )

pnbd_cdnow_twohierlambda_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

We first check the model against the in-sample data.

```{r run_pnbd_cdnow_twohierlambda_fit_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_cdnow_twohierlambda_assess_data_lst |>
  use_series(model_fit_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  bind_rows(sim_init_tbl) |>
  arrange(customer_id, draw_id, tnx_timestamp)

assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)
```

This fit looks reasonable and appears to capture most of the aspects of the
data used to fit it. Given that this is a synthetic dataset, this is not
surprising, but at least we appreciate that our model is valid.


### Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

```{r run_pnbd_cdnow_twohierlambda_valid_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_cdnow_twohierlambda_assess_data_lst |>
  use_series(model_valid_index_filepath) |>
  read_rds() |>
  use_series(sim_file) |>
  map_dfr(read_rds) |>
  select(customer_id, draw_id, sim_data) |>
  unnest(sim_data) |>
  arrange(customer_id, draw_id, tnx_timestamp)


assess_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_transactions_tbl,
  simdata_tbl = simdata_tbl
  )

assess_plots_lst |> map(print)
```

As for our short time frame data, overall our model is working well.


# Compare Model Outputs

We have looked at each of the models individually, but it is also worth looking
at each of the models as a group.

We now want to combine both the `fit` and `valid` transaction sets to
calculate the summary statistics for both.

```{r calculate_observed_statistics}
#| echo: true

obs_summstats_tbl <- list(
    fit   = customer_fit_transactions_tbl,
    valid = customer_valid_transactions_tbl
    ) |>
  bind_rows(.id = "assess_type") |>
  group_by(assess_type) |>
  calculate_transaction_summary_statistics() |>
  pivot_longer(
    cols      = !assess_type,
    names_to  = "label",
    values_to = "obs_value"
    )

obs_summstats_tbl |> glimpse()
```



```{r load_model_assessment_data}
#| echo: TRUE

model_assess_transactions_tbl <- dir_ls("data", regexp = "pnbd_cdnow_(fixed|one|two).*_assess_.*index") |>
  enframe(name = NULL, value = "file_path") |>
  mutate(
    model_label = str_replace(file_path, "data/pnbd_cdnow_(.*?)_assess_.*", "\\1"),
    assess_type = if_else(str_detect(file_path, "_assess_fit_"), "fit", "valid"),
    
    assess_data = map(
      file_path, construct_model_assessment_data,
      
      .progress = "construct_assess_data"
      )
    ) |>
  select(model_label, assess_type, assess_data) |>
  unnest(assess_data)

model_assess_transactions_tbl |> glimpse()
```

We now want to calculate the transaction statistics on this full dataset, for
each separate draw.

```{r calculate_model_assessment_transaction_statistics}
model_assess_tbl <- model_assess_transactions_tbl |>
  group_by(model_label, assess_type, draw_id) |>
  calculate_transaction_summary_statistics()

model_assess_tbl |> glimpse()
```

We now combine all this data to create a number of different comparison plots
for the various summary statistics.


```{r construct_model_comparison_plot}
#! echo: TRUE

create_multiple_model_assessment_plot(
  obs_summstats_tbl, model_assess_tbl,
  "total_count", "Total Transactions"
  )

create_multiple_model_assessment_plot(
  obs_summstats_tbl, model_assess_tbl,
  "mean_count", "Average Transactions per Customer"
  )

create_multiple_model_assessment_plot(
  obs_summstats_tbl, model_assess_tbl,
  "p99", "99th Percentile Count"
  )
```


## Write Assessment Data to Disk

We now want to save the assessment data to disk.

```{r write_model_assessment_data}
#| echo: TRUE

model_assess_tbl |> write_rds("data/assess_data_pnbd_cdnow_twohier_tbl.rds")
```



# R Environment {.unnumbered}

```{r show_session_info}
#| echo: TRUE
#| message: TRUE

options(width = 120L)
sessioninfo::session_info()
options(width = 80L)
```