## ----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"))