construct_onlineretail_fixed_pnbd_models.qmd

---
title: "Construct Non-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(lobstr)
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 = parallelly::availableCores()
  )


set.seed(42)
stanfit_seed <- 5000

open_syslog("construct_onlineretail_fixed_pnbd_models")

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


In this workbook we construct our first 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 First P/NBD Model

We now construct our Stan model and prepare to fit it with our synthetic
dataset.

We also want to set a number of overall parameters for this workbook

To start the fit data, we want to use the 1,000 customers. We also need to
calculate the summary statistics for the validation period.


## Compile and Fit Stan Model

We now compile this model using `CmdStanR`.

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

pnbd_fixed_stanmodel <- cmdstan_model(
  "stan_code/pnbd_fixed.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_fixed1_stanmodel}
#| echo: TRUE

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

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(
    lambda_mn = 0.25,
    lambda_cv = 1.00,
    
    mu_mn     = 0.10,
    mu_cv     = 1.00,
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_onlineretail_fixed1_stanfit <- pnbd_fixed_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_fixed1_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_onlineretail_fixed1_stanfit <- read_rds(stanfit_object_file)
}

pnbd_onlineretail_fixed1_stanfit$summary()
```


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

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

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

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

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

pnbd_onlineretail_fixed1_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  as.numeric() |>
  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_fixed1_assessment}
#| echo: TRUE

pnbd_stanfit <- pnbd_onlineretail_fixed1_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_onlineretail_fixed1_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_onlineretail_fixed1",
  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         = 10
  )

pnbd_onlineretail_fixed1_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

simdata_tbl <- pnbd_onlineretail_fixed1_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_fixed1_valid_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_onlineretail_fixed1_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.

```{r delete_pnbd_onlineretail_fixed1_valid_assessment}
#| echo: false

rm(simdata_tbl)
rm(insample_plots_lst)
rm(outsample_plots_lst)

gc()
```



# Fit Alternate Prior Model.

We want to try an alternate prior model with a smaller co-efficient of variation
to see what impact it has on our procedures.


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

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

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(
    lambda_mn = 0.25,
    lambda_cv = 0.50,
    
    mu_mn     = 0.10,
    mu_cv     = 0.50,
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_onlineretail_fixed2_stanfit <- pnbd_fixed_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_fixed2_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_onlineretail_fixed2_stanfit <- read_rds(stanfit_object_file)
}

pnbd_onlineretail_fixed2_stanfit$summary()
```


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

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

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

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

pnbd_onlineretail_fixed2_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 want to check the $N_{eff}$ statistics also.


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

pnbd_onlineretail_fixed2_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  as.numeric() |>
  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_fixed2_assessment}
#| echo: TRUE

pnbd_stanfit <- pnbd_onlineretail_fixed2_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_onlineretail_fixed2_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_onlineretail_fixed2",
  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         = 20
  )

pnbd_onlineretail_fixed2_assess_data_lst |> glimpse()
```

### Check In-Sample Data Validation

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

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

simdata_tbl <- pnbd_onlineretail_fixed2_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_fixed2_valid_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_onlineretail_fixed2_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()
```




```{r delete_pnbd_onlineretail_fixed2_valid_assessment}
#| echo: false

rm(simdata_tbl)
rm(insample_plots_lst)
rm(outsample_plots_lst)

gc()
```


# Fit Tight-Lifetime Model

We now want to try a model where we use priors with a tighter coefficient of
variation for lifetime but keep the CoV for transaction frequency.


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

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

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(
    lambda_mn = 0.25,
    lambda_cv = 1.00,
    
    mu_mn     = 0.10,
    mu_cv     = 0.50,
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_onlineretail_fixed3_stanfit <- pnbd_fixed_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_fixed3_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_onlineretail_fixed3_stanfit <- read_rds(stanfit_object_file)
}

pnbd_onlineretail_fixed3_stanfit$summary()
```


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

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

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

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

pnbd_onlineretail_fixed3_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 want to check the $N_{eff}$ statistics also.


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

pnbd_onlineretail_fixed3_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  as.numeric() |>
  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_fixed3_assessment}
#| echo: TRUE

pnbd_stanfit <- pnbd_onlineretail_fixed3_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_onlineretail_fixed3_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_onlineretail_fixed3",
  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         = 30
  )

pnbd_onlineretail_fixed3_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

simdata_tbl <- pnbd_onlineretail_fixed3_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_fixed3_valid_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_onlineretail_fixed3_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()
```



```{r delete_pnbd_onlineretail_fixed3_valid_assessment}
#| echo: false

rm(simdata_tbl)
rm(insample_plots_lst)
rm(outsample_plots_lst)

gc()
```


# Fit Narrow-Short-Lifetime Model

We now want to try a model where we use priors with a tighter coefficient of
variation for lifetime but keep the CoV for transaction frequency.


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

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

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(
    lambda_mn = 0.25,
    lambda_cv = 1.00,
    
    mu_mn     = 0.20,
    mu_cv     = 0.30,
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_onlineretail_fixed4_stanfit <- pnbd_fixed_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_fixed4_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_onlineretail_fixed4_stanfit <- read_rds(stanfit_object_file)
}

pnbd_onlineretail_fixed4_stanfit$summary()
```


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

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

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

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

pnbd_onlineretail_fixed4_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 want to check the $N_{eff}$ statistics also.


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

pnbd_onlineretail_fixed4_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  as.numeric() |>
  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_fixed4_assessment}
#| echo: TRUE

pnbd_stanfit <- pnbd_onlineretail_fixed4_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_onlineretail_fixed4_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_onlineretail_fixed4",
  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         = 40
  )

pnbd_onlineretail_fixed4_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

simdata_tbl <- pnbd_onlineretail_fixed4_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_fixed4_valid_assessment}
#| echo: TRUE

simdata_tbl <- pnbd_onlineretail_fixed4_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()
```


```{r delete_pnbd_onlineretail_fixed4_valid_assessment}
#| echo: false

rm(simdata_tbl)
rm(insample_plots_lst)
rm(outsample_plots_lst)

gc()
```


# 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 <- obs_fitdata_tbl |>
  filter(tnx_count > 0) |>
  nrow()

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

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

obs_valid_total_count <- obs_validdata_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_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()
```

We have now constructed the simulation summary statistics and now reshape our
data to aid in our model assessment.


```{r calculate_model_assess_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()
```


We now use this data to construct model comparison plots for the different
models we have fit.


```{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 = 8)
    )
```


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





# R Environment


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

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