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

## ----echo=TRUE, message=FALSE, warning=FALSE----------------------------------
library(ccostr)
library(ggplot2)
library(knitr)
library(parallel)
library(msm)

## -----------------------------------------------------------------------------
set.seed(123)
sim <- simCostData(n = 1000, dist = "unif", censor = "light", cdist = "exp", L = 10)
est <- ccmean(sim$censoredCostHistory)
est

## ----fig.height=3, fig.width=7------------------------------------------------
plot(est) + geom_hline(yintercept = 40000, linetype = "dotted", size = 1)

## ----eval = FALSE-------------------------------------------------------------
# 
# nSim   <- 1000
# nYears <- 10
# indv   <- 1000 # increating individuals increases computing time exponential
# ## true mean for unif is 40000 and exp is 35956
# unif_light <- lapply(1:nSim, function(x) simCostData(n = indv, dist = "unif", censor = "light", cdist = "exp", L = nYears))
# unif_heavy <- lapply(1:nSim, function(x) simCostData(n = indv, dist = "unif", censor = "heavy", cdist = "exp", L = nYears))
# exp_light  <- lapply(1:nSim, function(x) simCostData(n = indv, dist = "exp",  censor = "light", cdist = "exp", L = nYears))
# exp_heavy  <- lapply(1:nSim, function(x) simCostData(n = indv, dist = "exp",  censor = "heavy", cdist = "exp", L = nYears))

## ----eval = FALSE-------------------------------------------------------------
# nCores <- parallel::detectCores() - 1
# cl <- makeCluster(nCores)
# clusterExport(cl = cl, c("unif_light", "unif_heavy", "exp_light", "exp_heavy"))
# invisible(clusterEvalQ(cl = cl, {library(dplyr)
#                                  library(ccostr)
#                                  library(data.table)
#                                  library(survival)}))
# est_unif_light <- parLapply(cl, unif_light, function(x) ccmean(x$censoredCostHistory, L = 10))
# est_unif_heavy <- parLapply(cl, unif_heavy, function(x) ccmean(x$censoredCostHistory, L = 10))
# est_exp_light  <- parLapply(cl, exp_light,  function(x) ccmean(x$censoredCostHistory, L = 10))
# est_exp_heavy  <- parLapply(cl, exp_heavy,  function(x) ccmean(x$censoredCostHistory, L = 10))
# stopCluster(cl)

## ----eval = FALSE-------------------------------------------------------------
# results_unif_light <- do.call(rbind, lapply(est_unif_light, function(x) x[[3]]))
# results_unif_heavy <- do.call(rbind, lapply(est_unif_heavy, function(x) x[[3]]))
# results_exp_light  <- do.call(rbind, lapply(est_exp_light,  function(x) x[[3]]))
# results_exp_heavy  <- do.call(rbind, lapply(est_exp_heavy,  function(x) x[[3]]))
# results_true <- data.frame("unif_light" = 40000,
#                            "unif_heavy" = 40000,
#                            "exp_light"  = 35956,
#                            "exp_heavy"  = 35956)
# results_bias <- data.frame("unif_light" = (colMeans(results_unif_light)),
#                            "unif_heavy" = (colMeans(results_unif_heavy)),
#                            "exp_light"  = (colMeans(results_exp_light)),
#                            "exp_heavy"  = (colMeans(results_exp_heavy)))
# results <- rbind(results_true, results_bias)
# row.names(results) <- c("true_mean", colnames(results_unif_light))
# 
# results_bias <- cbind(round(results[,c(1,2)] - 40000,2),
#                       round(results[,c(3,4)] - 35956,2))
# 
# cov_unif_light <- do.call(rbind, lapply(est_unif_light, function(x) ifelse(x[[4]][5,] >= 40000 & x[[4]][4,] <= 40000, 1, 0)))
# cov_unif_heavy <- do.call(rbind, lapply(est_unif_heavy, function(x) ifelse(x[[4]][5,] >= 40000 & x[[4]][4,] <= 40000, 1, 0)))
# cov_exp_light  <- do.call(rbind, lapply(est_exp_light,  function(x) ifelse(x[[4]][5,] >= 35956 & x[[4]][4,] <= 35956, 1, 0)))
# cov_exp_heavy  <- do.call(rbind, lapply(est_exp_heavy,  function(x) ifelse(x[[4]][5,] >= 35956 & x[[4]][4,] <= 35956, 1, 0)))
# results_coverage <- data.frame("unif_light" = (colMeans(cov_unif_light, na.rm = T)),
#                                "unif_heavy" = (colMeans(cov_unif_heavy, na.rm = T)),
#                                "exp_light"  = (colMeans(cov_exp_light,  na.rm = T)),
#                                "exp_heavy"  = (colMeans(cov_exp_heavy,  na.rm = T)))

## -----------------------------------------------------------------------------
load(system.file("extdata", "results.Rdata", package = "ccostr"))
kable(results)

## -----------------------------------------------------------------------------
kable(results_bias)

## -----------------------------------------------------------------------------
kable(results_coverage, digits = 3)