construct_onlineretail_twohier_pnbd_models.qmd

---
title: "Construct Single-Hierarchical P/NBD Model for Online Retail 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(rsyslog)
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 <- 4200

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

We start by modelling the P/NBD model using our synthetic datasets before we
try to model real-life data.



## Load Online Retail Data

We now want to load the online retail transaction data.


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

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

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

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




## Load Derived Data

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

customer_summarystats_tbl <- read_rds("data/onlineretail_customer_summarystats_tbl.rds")

obs_fitdata_tbl   <- read_rds("data/onlineretail_obs_fitdata_tbl.rds")
obs_validdata_tbl <- read_rds("data/onlineretail_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/onlineretail_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
```


Finally, we need to set our directories where we save our Stan code and the
model outputs.

```{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_onlineretail_twohiermean_stanmodel}
#| echo: TRUE

stan_modelname <- "pnbd_onlineretail_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_onlineretail_twohiermean_stanfit <- pnbd_twohiermean_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    adapt_delta     =                         0.95,
    save_warmup     =                        FALSE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_onlineretail_twohiermean_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_onlineretail_twohiermean_stanfit <- read_rds(stanfit_object_file)
}

pnbd_onlineretail_twohiermean_stanfit$summary()
```


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

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

pnbd_onlineretail_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_onlineretail_twohiermean_lambda_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_onlineretail_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_onlineretail_twohiermean_parameter_neffratio}
#| echo: TRUE

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


## 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_onlineretail_twohiermean_assessment}
#| echo: TRUE

pnbd_stanfit <- pnbd_onlineretail_twohiermean_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_onlineretail_twohiermean_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_onlineretail_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         = 4210
  )

pnbd_onlineretail_twohiermean_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

simdata_tbl <- pnbd_onlineretail_twohiermean_assess_data_lst |>
  use_series(model_fit_simstats_filepath) |>
  read_rds()

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_fitdata_tbl,
  simdata_tbl = simdata_tbl
  )

insample_plots_lst$multi_plot |> print()
insample_plots_lst$total_plot |> print()
insample_plots_lst$quant_plot |> 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_onlineretail_twohiermean_valid_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_onlineretail_twohiermean_assess_data_lst |>
  use_series(model_valid_simstats_filepath) |>
  read_rds()

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_validdata_tbl,
  simdata_tbl = simdata_tbl
  )

outsample_plots_lst$multi_plot |> print()
outsample_plots_lst$total_plot |> print()
outsample_plots_lst$quant_plot |> 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_onlineretail_twohierlambda_stanmodel}
#| echo: TRUE

stan_modelname <- "pnbd_onlineretail_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_onlineretail_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_onlineretail_twohierlambda_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_onlineretail_twohierlambda_stanfit <- read_rds(stanfit_object_file)
}

pnbd_onlineretail_twohierlambda_stanfit$summary()
```


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

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

pnbd_onlineretail_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_onlineretail_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_onlineretail_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_onlineretail_twohierlambda_parameter_neffratio}
#| echo: TRUE

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


## 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_onlineretail_twohierlambda_assessment}
#| echo: TRUE

pnbd_stanfit <- pnbd_onlineretail_twohierlambda_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_onlineretail_twohierlambda_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_onlineretail_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_onlineretail_twohierlambda_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

simdata_tbl <- pnbd_onlineretail_twohierlambda_assess_data_lst |>
  use_series(model_fit_simstats_filepath) |>
  read_rds()

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_fitdata_tbl,
  simdata_tbl = simdata_tbl
  )

insample_plots_lst$multi_plot |> print()
insample_plots_lst$total_plot |> print()
insample_plots_lst$quant_plot |> 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_onlineretail_twohierlambda_valid_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_onlineretail_twohierlambda_assess_data_lst |>
  use_series(model_valid_simstats_filepath) |>
  read_rds()

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_validdata_tbl,
  simdata_tbl = simdata_tbl
  )

outsample_plots_lst$multi_plot |> print()
outsample_plots_lst$total_plot |> print()
outsample_plots_lst$quant_plot |> 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.

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

calculate_simulation_statistics <- function(file_rds) {
  simdata_tbl <- read_rds(file_rds)
  
  multicount_cust_tbl <- simdata_tbl |>
    filter(sim_tnx_count > 0) |>
    count(draw_id, name = "multicust_count")
  
  totaltnx_data_tbl <- simdata_tbl |>
    count(draw_id, wt = sim_tnx_count, name = "simtnx_count")
  
  simstats_tbl <- multicount_cust_tbl |>
    inner_join(totaltnx_data_tbl, by = "draw_id")
  
  return(simstats_tbl)
}
```



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

obs_fit_customer_count <- customer_fit_subset_tbl |>
  filter(tnx_count > 0) |>
  nrow()

obs_valid_customer_count <- customer_valid_subset_tbl |>
  filter(tnx_count > 0) |>
  nrow()

obs_fit_total_count <- customer_fit_subset_tbl |>
  pull(tnx_count) |>
  sum()

obs_valid_total_count <- customer_valid_subset_tbl |>
  pull(tnx_count) |>
  sum()


obs_stats_tbl <- tribble(
  ~assess_type, ~name,               ~obs_value,
  "fit",        "multicust_count",   obs_fit_customer_count,
  "fit",        "simtnx_count",      obs_fit_total_count,
  "valid",      "multicust_count",   obs_valid_customer_count,
  "valid",      "simtnx_count",      obs_valid_total_count
  )


model_assess_tbl <- dir_ls("data", regexp = "pnbd_onlineretail_(one|two|fixed).*_assess_.*simstats") |>
  enframe(name = NULL, value = "file_path") |>
  filter(str_detect(file_path, "_assess_model_", negate = TRUE)) |>
  mutate(
    model_label = str_replace(file_path, "data/pnbd_onlineretail_(.*?)_assess_.*", "\\1"),
    assess_type = if_else(str_detect(file_path, "_assess_fit_"), "fit", "valid"),
    
    sim_data = map(
      file_path, calculate_simulation_statistics,
      
      .progress = "calculate_simulation_statistics"
      )
    )

model_assess_tbl |> glimpse()
```



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

model_assess_summstat_tbl <- model_assess_tbl |>
  select(model_label, assess_type, sim_data) |>
  unnest(sim_data) |>
  pivot_longer(
    cols = !c(model_label, assess_type, draw_id)
    ) |>
  group_by(model_label, assess_type, name) |>
  summarise(
    .groups = "drop",
    
    mean_val = mean(value),
    p10 = quantile(value, 0.10),
    p25 = quantile(value, 0.25),
    p50 = quantile(value, 0.50),
    p75 = quantile(value, 0.75),
    p90 = quantile(value, 0.90)
    )

model_assess_summstat_tbl |> glimpse()
```


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

ggplot(model_assess_summstat_tbl) +
  geom_errorbar(
    aes(x = model_label, ymin = p10, ymax = p90), width = 0
    ) +
  geom_errorbar(
    aes(x = model_label, ymin = p25, ymax = p75), width = 0, linewidth = 3
    ) +
  geom_hline(
    aes(yintercept = obs_value),
    data = obs_stats_tbl, colour = "red"
    ) +
  scale_y_continuous(labels = label_comma()) +
  expand_limits(y = 0) +
  facet_wrap(
    vars(assess_type, name), scale = "free_y"
    ) +
  labs(
    x = "Model",
    y = "Count",
    title = "Comparison Plot for the Different Models"
    ) +
  theme(
    axis.text.x = element_text(angle = 20, vjust = 0.5, size = 6)
    )
```


## 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_onlineretail_twohier_tbl.rds")
```




# R Environment

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

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