construct_cdnow_fixed_pnbd_models.qmd

---
title: "Construct Non-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(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 <- 5000

open_syslog("construct_cdnow_fixed_pnbd_models")

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

In this workbook we construct the non-hierarchical P/NBD models on the CD-NOW
transaction data.




# 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
```


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

stan_modelname <- "pnbd_cdnow_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_cdnow_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_cdnow_fixed1_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_cdnow_fixed1_stanfit <- read_rds(stanfit_object_file)
}

pnbd_cdnow_fixed1_stanfit$summary()
```


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

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

pnbd_cdnow_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_cdnow_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_cdnow_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_cdnow_fixed1_parameter_neffratio}
#| echo: TRUE

pnbd_cdnow_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.

We first run the assessment data.

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

pnbd_stanfit <- pnbd_cdnow_fixed1_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_cdnow_fixed1_assess_data_lst <- run_model_assessment(
  model_stanfit       = pnbd_stanfit,
  insample_tbl        = customer_fit_subset_tbl,
  fit_label           = "pnbd_cdnow_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(),
  precompute_rootdir  = "precompute",
  data_dir            = "data",
  summary_include_tnx = FALSE,
  sim_seed            = 10
  )

pnbd_cdnow_fixed1_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

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

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_subset_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_cdnow_fixed1_valid_assessment}
#| echo: TRUE

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

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_subset_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_cdnow_fixed1_assessment}
#| echo: false

rm(simdata_tbl)
```


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

stan_modelname <- "pnbd_cdnow_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_cdnow_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_cdnow_fixed2_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_cdnow_fixed2_stanfit <- read_rds(stanfit_object_file)
}

pnbd_cdnow_fixed2_stanfit$summary()
```


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

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

pnbd_cdnow_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_cdnow_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_cdnow_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_cdnow_fixed2_parameter_neffratio}
#| echo: TRUE

pnbd_cdnow_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.

We first run the assessment data.

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

pnbd_stanfit <- pnbd_cdnow_fixed2_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_cdnow_fixed2_assess_data_lst <- run_model_assessment(
  model_stanfit       = pnbd_stanfit,
  insample_tbl        = customer_fit_subset_tbl,
  fit_label           = "pnbd_cdnow_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(),
  precompute_rootdir  = "precompute",
  data_dir            = "data",
  summary_include_tnx = FALSE,
  sim_seed            = 20
  )

pnbd_cdnow_fixed2_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

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

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_subset_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_cdnow_fixed2_valid_assessment}
#| echo: TRUE

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

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_subset_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_cdnow_fixed2_assessment}
#| echo: false

rm(simdata_tbl)
```




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

stan_modelname <- "pnbd_cdnow_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_cdnow_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_cdnow_fixed3_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_cdnow_fixed3_stanfit <- read_rds(stanfit_object_file)
}
  
pnbd_cdnow_fixed3_stanfit$summary()
```


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

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

pnbd_cdnow_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_cdnow_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_cdnow_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_cdnow_fixed3_parameter_neffratio}
#| echo: TRUE

pnbd_cdnow_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.

We first run the assessment data.

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

pnbd_stanfit <- pnbd_cdnow_fixed3_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_cdnow_fixed3_assess_data_lst <- run_model_assessment(
  model_stanfit       = pnbd_stanfit,
  insample_tbl        = customer_fit_subset_tbl,
  fit_label           = "pnbd_cdnow_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(),
  precompute_rootdir  = "precompute",
  data_dir            = "data",
  summary_include_tnx = FALSE,
  sim_seed            = 30
  )

pnbd_cdnow_fixed3_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

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

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_subset_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_cdnow_fixed3_valid_assessment}
#| echo: TRUE

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

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_subset_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_cdnow_fixed3_assessment}
#| echo: false

rm(simdata_tbl)
```




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

stan_modelname <- "pnbd_cdnow_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_cdnow_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_cdnow_fixed4_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_cdnow_fixed4_stanfit <- read_rds(stanfit_object_file)
}

pnbd_cdnow_fixed4_stanfit$summary()
```


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

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

pnbd_cdnow_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_cdnow_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_cdnow_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_cdnow_fixed4_parameter_neffratio}
#| echo: TRUE

pnbd_cdnow_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.

We first run the assessment data.

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

pnbd_stanfit <- pnbd_cdnow_fixed4_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_cdnow_fixed4_assess_data_lst <- run_model_assessment(
  model_stanfit       = pnbd_stanfit,
  insample_tbl        = customer_fit_subset_tbl,
  fit_label           = "pnbd_cdnow_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(),
  precompute_rootdir  = "precompute",
  data_dir            = "data",
  summary_include_tnx = FALSE,
  sim_seed            = 40
  )

pnbd_cdnow_fixed4_assess_data_lst |> glimpse()
```


### Check In-Sample Data Validation

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

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

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

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_fit_subset_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_cdnow_fixed4_valid_assessment}
#| echo: TRUE

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

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = customer_valid_subset_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_cdnow_fixed4_assessment}
#| echo: false

rm(simdata_tbl)
```


# 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_cdnow_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_cdnow_(.*?)_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 = 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_cdnow_fixed_tbl.rds")
```




# R Environment

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

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