BayesMoFo-vignette.R
In BayesMoFo: Bayesian Mortality Forecasting

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
# install package
# Recommended installation
# install.packages("BayesMoFo")

# Development version (use only if needed)
# install.packages("devtools")
# devtools::install_github("jstw1g09/Rpackage-BayesMoFo")

#load package
library(BayesMoFo)

## -----------------------------------------------------------------------------
data(uk_mortalitydata)

## -----------------------------------------------------------------------------
head(uk_mortalitydata, n = 20)

## ----echo = FALSE, eval = FALSE-----------------------------------------------
#  death <- preparedata_fn(data_summarised[, c("Age", "Year", "Claim")],
#                                    ages = 35:65, years = 2016:2020)
#  expo <- preparedata_fn(data_summarised[, c("Age", "Year", "Exposure")],
#                                  ages = 35:65, years = 2016:2020)

## -----------------------------------------------------------------------------
death <- preparedata_fn(uk_mortalitydata[, c("Age", "Year", "Deaths")], 
                                  ages = 30:60, years = 2000:2020)
expo <- preparedata_fn(uk_mortalitydata[, c("Age", "Year", "Exposures")], 
                                ages = 30:60, years = 2000:2020)

## -----------------------------------------------------------------------------
data("dxt_array_product")
data("Ext_array_product")

# preview of death data the 1st insurance product called "ACI"
str(dxt_array_product["ACI",,,drop = FALSE])

## ----eval = FALSE-------------------------------------------------------------
#  
#  # inputting the data as a 3-way array
#  death <- preparedata_fn(dxt_array_product["ACI",,,drop = FALSE], ages = 35:65, years = 2016:2020)
#  expo <- preparedata_fn(Ext_array_product["ACI",,,drop = FALSE], ages = 35:65, years = 2016:2020)
#  
#  # specifying the name of the stratum to load using `strat_name`
#  death <- preparedata_fn(dxt_array_product,strat_name="ACI", ages = 35:65, years = 2016:2020)
#  expo <- preparedata_fn(Ext_array_product,strat_name="ACI", ages = 35:65, years = 2016:2020)

## -----------------------------------------------------------------------------
# preview of death data the 1st insurance product called "ACI"
str(dxt_array_product["ACI",,,drop = TRUE])

## ----eval = FALSE-------------------------------------------------------------
#  death<-preparedata_fn(dxt_array_product["ACI",,,drop = TRUE],data_matrix=TRUE,ages=35:65)
#  expo<-preparedata_fn(Ext_array_product["ACI",,,drop = TRUE],data_matrix=TRUE,ages=35:65)

## ----eval = FALSE-------------------------------------------------------------
#  fitmodel <- runBayesMoFo(death, expo,
#                           models = c("LC",
#                                      "CBD_M3",
#                                      "APCI")
#                           )

## ----eval = FALSE-------------------------------------------------------------
#  fitmodel <- fit_LC(death, expo)

## ----eval = FALSE-------------------------------------------------------------
#  fitmodel <- runBayesMoFo(death, expo, models = "LC")

## ----eval = FALSE-------------------------------------------------------------
#  fitmodel <- runBayesMoFo(death, expo,
#                           models = c("LC",
#                                      "CBD_M3",
#                                      "APCI"),
#                           family="poisson"
#                           )

## ----message = FALSE, warning = FALSE-----------------------------------------
fitmodel_forecast <- runBayesMoFo(death, expo,
                         models = c("LC",
                                    "CBD_M3",
                                    "APCI"),
                         forecast = TRUE,
                         h = 6, 
                         n.chain = 2
                         )

## -----------------------------------------------------------------------------
fitmodel_forecast$best_model
fitmodel_forecast$worst_model

## -----------------------------------------------------------------------------
fitmodel_forecast$DIC

## ----eval = FALSE-------------------------------------------------------------
#  fitmodel_forecast$result$best
#  fitmodel_forecast$result$worst

## ----fig.width = 10, fig.height = 10------------------------------------------
plot_param_fn(fitmodel_forecast)

## ----eval = FALSE-------------------------------------------------------------
#  plot_param_fn(fitmodel_forecast, pred_int = 0.80)

## ----fig.width = 10, fig.height = 10------------------------------------------
plot_param_fn(fitmodel_forecast, pred_int = 0.80, legends = FALSE)

## ----fig.width = 10, fig.height = 10------------------------------------------
plot_rates_fn(fitmodel_forecast)

## ----fig.width = 10, fig.height = 5-------------------------------------------
plot_rates_fn(fitmodel_forecast, plot_years = c(2016,2020,2024))

## ----fig.width = 10, fig.height = 5-------------------------------------------
plot_rates_fn(fitmodel_forecast, plot_type = "time", plot_ages = c(35,45,55))

## -----------------------------------------------------------------------------
summary_fitmodel<-summary_fn(fitmodel_forecast)

## ----eval = FALSE-------------------------------------------------------------
#  #posterior means
#  summary_fitmodel$rates_summary$mean
#  
#  #posterior standard deviations
#  summary_fitmodel$rates_summary$std

## ----eval = FALSE-------------------------------------------------------------
#  #posterior medians
#  summary_fitmodel$rates_pn$median
#  
#  #lower quantiles
#  summary_fitmodel$rates_pn$lower
#  
#  #upper quantiles
#  summary_fitmodel$rates_pn$upper

## ----eval = FALSE-------------------------------------------------------------
#  #posterior means
#  summary_fitmodel$param_summary$mean
#  
#  #posterior standard deviations
#  summary_fitmodel$param_summary$std
#  
#  #posterior medians
#  summary_fitmodel$param_pn$median
#  
#  #lower quantiles
#  summary_fitmodel$param_pn$lower
#  
#  #upper quantiles
#  summary_fitmodel$param_pn$upper

## ----fig.width = 10, fig.height = 10------------------------------------------
diagnostics_rates_result<-converge_diag_rates_fn(fitmodel_forecast)

## -----------------------------------------------------------------------------
diagnostics_rates_result$ESS

## ----fig.width = 10, fig.height = 10------------------------------------------
converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(35,45,55), plot_years = c(2016,2020,2024))

## ----eval = FALSE-------------------------------------------------------------
#  #for only trace plots
#  converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(35,45,55), plot_years = c(2016,2020,2024), trace = TRUE, density = FALSE)

## ----fig.width = 10, fig.height = 10------------------------------------------
#for only acf plots
converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(35,45,55), plot_years = c(2016,2020,2024), trace = FALSE, density = FALSE, acf_plot = TRUE)

## ----fig.width = 10, fig.height = 10------------------------------------------
diagnostics_param_result<-converge_diag_param_fn(fitmodel_forecast)

## -----------------------------------------------------------------------------
diagnostics_param_result$ESS

## ----fig.width = 10, fig.height = 10------------------------------------------
converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa","gamma","rho","phi","sigma2_kappa"))

## -----------------------------------------------------------------------------
fitmodel_forecast$result$best$param

## ----fig.width = 9, fig.height = 5--------------------------------------------
converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa[1,4]","gamma[1,2]"))

## -----------------------------------------------------------------------------
colnames(fitmodel_forecast$result$best$post_sample[[1]])[!startsWith(colnames(fitmodel_forecast$result$best$post_sample[[1]]),"q[")]

## ----eval = FALSE-------------------------------------------------------------
#  #for only trace plots
#  converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa[1,4]","gamma[1,2]"), trace = TRUE, density = FALSE)

## ----fig.width = 5, fig.height = 5--------------------------------------------
#for only acf plots
converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa[1,4]","gamma[1,2]"), trace = FALSE, density = FALSE, acf_plot = TRUE)

## ----fig.width = 9, fig.height = 5--------------------------------------------
converge_diag_result<-converge_diag_fn(fitmodel_forecast, plot_gelman = TRUE, plot_geweke = TRUE)

## -----------------------------------------------------------------------------
#Gelman's R
head(converge_diag_result$gelman_diag$psrf)

#Geweke's Z
head(converge_diag_result$geweke_diag$z)

#Heidel's Stationarity and Half-width tests
head(converge_diag_result$heidel_diag)

## -----------------------------------------------------------------------------
data(uk_deathscausedata)
head(uk_deathscausedata, n = 10)

## -----------------------------------------------------------------------------
death <- preparedata_fn(uk_deathscausedata[,c("Age","Year","Deaths","Cause")])
expo <- preparedata_fn(uk_deathscausedata[,c("Age","Year","Exposures","Cause")])

#or if require a subset of the data
death <- preparedata_fn(uk_deathscausedata[,c("Age","Year","Deaths","Cause")], 
                                 ages = seq(45,90,by=5), years = 2001:2020)
expo <- preparedata_fn(uk_deathscausedata[,c("Age","Year","Exposures","Cause")], 
                                ages = seq(45,90,by=5), years = 2001:2020)

str(death)
str(expo)

## ----eval = FALSE-------------------------------------------------------------
#  data("dxt_array_product");data("Ext_array_product")
#  str(dxt_array_product) # 3D data array
#  
#  death<-preparedata_fn(dxt_array_product,ages=35:65)
#  expo<-preparedata_fn(Ext_array_product,ages=35:65)

## ----message = FALSE, warning = FALSE-----------------------------------------
fitmodel_forecast <- runBayesMoFo(death, expo,
                         models = "MLiLee",
                         forecast = TRUE,
                         h = 5,
                         quiet = TRUE,
                         n.chain = 2
                         )

## ----eval = FALSE, echo = FALSE-----------------------------------------------
#  fitmodel_forecast$best_model
#  fitmodel_forecast$worst_model
#  fitmodel_forecast$DIC

## ----fig.width = 10, fig.height = 10------------------------------------------
plot_param_fn(fitmodel_forecast)

## ----fig.width = 10, fig.height = 5-------------------------------------------
plot_rates_fn(fitmodel_forecast, plot_years = c(2005,2020,2025))

## ----fig.width = 10, fig.height = 5-------------------------------------------
plot_rates_fn(fitmodel_forecast, plot_type = "time", plot_ages = c(45,65,85))

## ----fig.width = 10, fig.height = 10------------------------------------------
diagnostics_param_result<-converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(45,65,85), plot_years = c(2005,2020,2025))

## ----fig.width = 10, fig.height = 10------------------------------------------
#for only acf plots
converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(45,65,85), plot_years = c(2005,2020,2025), trace = FALSE, density = FALSE, acf_plot = TRUE)

## -----------------------------------------------------------------------------
diagnostics_param_result$ESS

## ----fig.width = 10, fig.height = 10------------------------------------------
converge_diag_param_fn(fitmodel_forecast, plot_params = c("beta","kappa","rho","phi","sigma2_kappa"))

## -----------------------------------------------------------------------------
fitmodel_forecast$result$best$param

## ----fig.width = 9, fig.height = 5--------------------------------------------
converge_diag_param_fn(fitmodel_forecast, plot_params = c("beta[1,3]","kappa[2,4]"))

## -----------------------------------------------------------------------------
colnames(fitmodel_forecast$result$best$post_sample[[1]])[!startsWith(colnames(fitmodel_forecast$result$best$post_sample[[1]]),"q[")]

## ----eval= FALSE--------------------------------------------------------------
#  #for only trace plots
#  converge_diag_param_fn(fitmodel_forecast, plot_params = c("beta[1,3]","kappa[2,4]"), trace = TRUE, density = FALSE)

## ----fig.width = 5, fig.height = 5--------------------------------------------
#for only acf plots
converge_diag_param_fn(fitmodel_forecast, plot_params = c("beta[1,3]","kappa[2,4]"), trace = FALSE, density = FALSE, acf_plot = TRUE)

## ----eval= FALSE--------------------------------------------------------------
#  converge_diag_result<-converge_diag_fn(fitmodel_forecast, plot_gelman = TRUE, plot_geweke = TRUE)

## ----eval= FALSE--------------------------------------------------------------
#  #Gelman's R
#  head(converge_diag_result$gelman_diag$psrf)
#  
#  #Geweke's Z
#  head(converge_diag_result$geweke_diag$z)
#  
#  #Heidel's Stationarity and Half-width tests
#  head(converge_diag_result$heidel_diag)

## ----eval = FALSE-------------------------------------------------------------
#  fitmodel_forecast <- runBayesMoFo(death, expo,
#                           models = c("APCI_sharegamma",
#                                      "RH_sharegamma"),
#                           forecast = TRUE,
#                           h = 5,
#                           quiet = TRUE
#                           )

## ----eval = FALSE-------------------------------------------------------------
#  fitmodel_forecast$DIC
#  
#  fitmodel_forecast$best_model
#  
#  plot_param_fn(fitmodel_forecast)

## ----eval = FALSE-------------------------------------------------------------
#  
#  plot_rates_fn(fitmodel_forecast, plot_years = c(2005,2020,2025))
#  
#  plot_rates_fn(fitmodel_forecast, plot_type = "time", plot_ages = c(45,65,85))
#  

## ----eval = FALSE-------------------------------------------------------------
#  converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(45,65,85), plot_years = c(2005,2020,2025))
#  
#  converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa","rho","phi","sigma2_kappa"))