diff --git a/.Rbuildignore b/.Rbuildignore
index 3cf63732..6f184bf9 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -23,3 +23,4 @@ vignettes/benchmarks*
^vignettes/precompile.R
^CRAN-SUBMISSION$
^Makefile$
+^CODEOWNERS$
diff --git a/DESCRIPTION b/DESCRIPTION
index a12e6ea6..86dafafc 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Package: hesim
Type: Package
Title: Health Economic Simulation Modeling and Decision Analysis
-Version: 0.5.5
+Version: 0.5.6
Authors@R: c(
person("Devin", "Incerti", , "devin.incerti@gmail.com", role = c("aut", "cre")),
person("Jeroen P.", "Jansen", role = "aut"),
diff --git a/R/RcppExports.R b/R/RcppExports.R
index ac42a7a6..fee8dce4 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -69,6 +69,10 @@ C_indiv_ctstm_starting <- function(R_disease_prog, strategy_idx, patient_idx, R_
.Call('_hesim_C_indiv_ctstm_starting', PACKAGE = 'hesim', R_disease_prog, strategy_idx, patient_idx, R_StateVal, dr, type)
}
+C_indiv_ctstm_trans <- function(R_CtstmTrans, R_disease_prog, strategy_idx, patient_idx, R_StateVal, dr, type) {
+ .Call('_hesim_C_indiv_ctstm_trans', PACKAGE = 'hesim', R_CtstmTrans, R_disease_prog, strategy_idx, patient_idx, R_StateVal, dr, type)
+}
+
C_indiv_ctstm_los <- function(R_disease_prog, strategy_idx, patient_idx, dr) {
.Call('_hesim_C_indiv_ctstm_los', PACKAGE = 'hesim', R_disease_prog, strategy_idx, patient_idx, dr)
}
diff --git a/R/ctstm.R b/R/ctstm.R
index 44dc8ba6..4b709d28 100644
--- a/R/ctstm.R
+++ b/R/ctstm.R
@@ -104,7 +104,8 @@ indiv_ctstm_sim_disease <- function(trans_model, max_t = 100, max_age = 100,
setattr(disprog, "class",
c("disprog", "data.table", "data.frame"))
setattr(disprog, "size",
- c(get_size(trans_model), n_states = nrow(trans_model$trans_mat)))
+ c(get_size(trans_model), n_states = nrow(trans_model$trans_mat),
+ n_transitions = max(trans_model$trans_mat, na.rm = TRUE)))
setattr(disprog, "absorbing", absorbing(trans_model$trans_mat))
return(disprog[, ])
}
@@ -577,8 +578,15 @@ IndivCtstm <- R6::R6Class("IndivCtstm",
private$disprog_idx$patient_idx,
stateval_list[[i]], dr[j],
sim_type)
+ } else if (stateval_list[[i]]$method == "transition"){
+ C_ev <- C_indiv_ctstm_trans(self$trans_model,
+ self$disprog_,
+ private$disprog_idx$strategy_idx,
+ private$disprog_idx$patient_idx,
+ stateval_list[[i]], dr[j],
+ sim_type)
} else{
- stop("The 'StateVals' 'method' must either be 'wlos' or 'starting'.")
+ stop("The 'StateVals' 'method' must be 'wlos', 'starting' or 'transition'.")
}
self$disprog_[, ev := C_ev]
if (lys){
diff --git a/R/sim-general.R b/R/sim-general.R
index 942d0473..dd9fb2b1 100644
--- a/R/sim-general.R
+++ b/R/sim-general.R
@@ -23,8 +23,8 @@
#'
#' The object also contains `size` and `absorbing` attributes.
#' The `size` attribute is a numeric vector with the elements `n_samples`,
-#' `n_strategies`, `n_patients`, and `n_states` denoting the number of samples,
-#' treatment strategies, patients, and health states The `absorbing` attribute
+#' `n_strategies`, `n_patients`, `n_states` and `n_transitions` denoting the number of samples,
+#' treatment strategies, patients, health states and transitions. The `absorbing` attribute
#' is a numeric vector containing the absorbing health states; i.e., the
#' health states that cannot be transitioned from. Operationally, an
#' absorbing state is a row in a transition matrix (as in the `trans_mat` field
@@ -433,7 +433,8 @@ sim_ev.NULL <- function(object, ...) {
#' the cohort. The method used to simulate expected values depends on the
#' `$method` field in the [`StateVals`] object(s) stored in `model(s)`. If
#' `$method = "starting"`, then state values represent a one-time value that
-#' occurs at time 0.
+#' occurs at time 0. If `$method = "transition"`, then the state values represent
+#' a value for the transition at the transition time.
#'
#' The more common use case is `$method = "wlos"`, or a "weighted length of stay".
#' That is, expected values for each health state can be thought of as state values
diff --git a/R/statevals.R b/R/statevals.R
index e9e043fc..f81a54c3 100644
--- a/R/statevals.R
+++ b/R/statevals.R
@@ -32,7 +32,9 @@ StateVals <- R6::R6Class(
#' used in a cohort model, the state values only accrue at time 0;
#' in contrast, in an individual-level model, state values
#' accrue each time a patient enters a new state and are discounted based on
- #' time of entrance into that state.
+ #' time of entrance into that state. When `transition`, then state values
+ #' represent a value that occurs at the transition time; only used for
+ #' individual-level models.
method = NULL,
#' @field time_reset If `FALSE` then time intervals are based on time since
@@ -49,7 +51,7 @@ StateVals <- R6::R6Class(
#' @param time_reset The `time_reset` field.
#' @return A new `StateVals` object.
initialize = function(params, input_data = NULL,
- method = c("wlos", "starting"),
+ method = c("wlos", "starting", "transition"),
time_reset = FALSE) {
self$params <- params
self$input_data <- input_data
@@ -362,9 +364,12 @@ create_StateVals.lm <- function(object, input_data = NULL, n = 1000,
}
#' @rdname create_StateVals
+#' @param method String matching a list of methods used for the `StateVals` object.
#' @export
-create_StateVals.stateval_tbl <- function(object, hesim_data = NULL, n = 1000, ...){
-
+create_StateVals.stateval_tbl <- function(object, hesim_data = NULL, n = 1000,
+ method = c("wlos","starting","transition"),
+ ...){
+ method <- match.arg(method)
grp_var <- attr(object, "grp_var")
# For backwards compatibility, use hesim_data attribute of object
@@ -505,21 +510,24 @@ create_StateVals.stateval_tbl <- function(object, hesim_data = NULL, n = 1000, .
} else{
time_intervals <- NULL
}
-
+
+ transp <- (method == "transition")
tparams <- new_tparams_mean(value = mu,
n_samples = n,
strategy_id = tbl$strategy_id,
n_strategies = length(unique(tbl$strategy_id)),
patient_id = tbl$patient_id,
n_patients = length(unique(tbl$patient_id)),
- state_id = tbl$state_id,
- n_states = length(unique(tbl$state_id)),
+ state_id = if (transp) NULL else tbl$state_id,
+ transition_id = if (transp) tbl$state_id else NULL,
+ n_states = if (transp) NULL else length(unique(tbl$state_id)),
+ n_transitions = if (transp) length(unique(tbl$state_id)) else NULL,
time_id = tbl$time_id,
time_intervals = time_intervals,
n_times = nrow(time_intervals),
grp_id = tbl$grp_id,
patient_wt = tbl$patient_wt)
- return(StateVals$new(params = tparams, ...))
+ return(StateVals$new(params = tparams, method=method, ...))
}
#' @export
diff --git a/R/utils.R b/R/utils.R
index f6482d0b..cd742d32 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -36,7 +36,7 @@ check_dr <- function(dr){
}
check_StateVals <- function(models, object,
- object_name = c("stateprobs", "disprog")) {
+ object_name = c("stateprobs", "disprog", "trans_model")) {
if(!is.list(models)){
stop("'models' must be a list", call. = FALSE)
@@ -84,14 +84,22 @@ check_StateVals <- function(models, object,
check_size(get_id_object(models[[i]])$n_patients,
expected_size[["n_patients"]],
z = "patients")
- check_size(
- get_id_object(models[[i]])$n_states,
- expected_states,
- msg = paste0("The number of states in each 'StateVals' model ",
- "must equal the number of states in the ",
- "'", object_name, "' object less the number of ",
- "absorbing states, which is ", expected_states, ".")
- )
+ if (models[[i]]$method == "transition") {
+ if (object_name %in% c("disprog", "trans_model")) {
+ check_size(get_id_object(models[[i]])$n_transitions,
+ expected_size[["n_transitions"]],
+ z = "transitions")
+ } else stop(paste0("Transition-specific 'StateVals' not available for object_name == '",
+ object_name,"'."))
+ } else # method != "transition"
+ check_size(
+ get_id_object(models[[i]])$n_states,
+ expected_states,
+ msg = paste0("The number of states in each 'StateVals' model ",
+ "must equal the number of states in the ",
+ "'", object_name, "' object less the number of ",
+ "absorbing states, which is ", expected_states, ".")
+ )
} # End loop over state value models
}
diff --git a/man/StateVals.Rd b/man/StateVals.Rd
index 365e57af..6f8889ff 100644
--- a/man/StateVals.Rd
+++ b/man/StateVals.Rd
@@ -57,7 +57,9 @@ that occurs when a patient enters a health state. When \code{starting} is
used in a cohort model, the state values only accrue at time 0;
in contrast, in an individual-level model, state values
accrue each time a patient enters a new state and are discounted based on
-time of entrance into that state.}
+time of entrance into that state. When \code{transition}, then state values
+represent a value that occurs at the transition time; only used for
+individual-level models.}
\item{\code{time_reset}}{If \code{FALSE} then time intervals are based on time since
the start of the simulation. If \code{TRUE}, then time intervals reset each
@@ -84,7 +86,7 @@ Create a new \code{StateVals} object.
\if{html}{\out{}}\preformatted{StateVals$new(
params,
input_data = NULL,
- method = c("wlos", "starting"),
+ method = c("wlos", "starting", "transition"),
time_reset = FALSE
)}\if{html}{\out{
}}
}
diff --git a/man/create_StateVals.Rd b/man/create_StateVals.Rd
index 8dbc362d..dc921f63 100644
--- a/man/create_StateVals.Rd
+++ b/man/create_StateVals.Rd
@@ -16,7 +16,13 @@ create_StateVals(object, ...)
...
)
-\method{create_StateVals}{stateval_tbl}(object, hesim_data = NULL, n = 1000, ...)
+\method{create_StateVals}{stateval_tbl}(
+ object,
+ hesim_data = NULL,
+ n = 1000,
+ method = c("wlos", "starting", "transition"),
+ ...
+)
}
\arguments{
\item{object}{A model object of the appropriate class.}
@@ -34,6 +40,8 @@ documentation in \code{\link[=create_params]{create_params()}}.}
\item{hesim_data}{A \code{\link{hesim_data}} object. Only required when \code{object} is of class
\code{\link{stateval_tbl}}. See "details".}
+
+\item{method}{String matching a list of methods used for the \code{StateVals} object.}
}
\value{
A \code{\link{StateVals}} object.
diff --git a/man/disprog.Rd b/man/disprog.Rd
index 89f91bf9..e410e886 100644
--- a/man/disprog.Rd
+++ b/man/disprog.Rd
@@ -27,8 +27,8 @@ state during the simulation and 0 otherwise.}
The object also contains \code{size} and \code{absorbing} attributes.
The \code{size} attribute is a numeric vector with the elements \code{n_samples},
-\code{n_strategies}, \code{n_patients}, and \code{n_states} denoting the number of samples,
-treatment strategies, patients, and health states The \code{absorbing} attribute
+\code{n_strategies}, \code{n_patients}, \code{n_states} and \code{n_transitions} denoting the number of samples,
+treatment strategies, patients, health states and transitions. The \code{absorbing} attribute
is a numeric vector containing the absorbing health states; i.e., the
health states that cannot be transitioned from. Operationally, an
absorbing state is a row in a transition matrix (as in the \code{trans_mat} field
diff --git a/man/sim_ev.Rd b/man/sim_ev.Rd
index ef966661..4e546545 100644
--- a/man/sim_ev.Rd
+++ b/man/sim_ev.Rd
@@ -89,7 +89,8 @@ the \code{\link{CohortDtstm}} and \code{\link{Psm}} classes) are mean outcomes f
the cohort. The method used to simulate expected values depends on the
\verb{$method} field in the \code{\link{StateVals}} object(s) stored in \code{model(s)}. If
\verb{$method = "starting"}, then state values represent a one-time value that
-occurs at time 0.
+occurs at time 0. If \verb{$method = "transition"}, then the state values represent
+a value for the transition at the transition time.
The more common use case is \verb{$method = "wlos"}, or a "weighted length of stay".
That is, expected values for each health state can be thought of as state values
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 37e3fb61..b462e3e0 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -261,6 +261,23 @@ BEGIN_RCPP
return rcpp_result_gen;
END_RCPP
}
+// C_indiv_ctstm_trans
+std::vector C_indiv_ctstm_trans(Rcpp::Environment R_CtstmTrans, Rcpp::DataFrame R_disease_prog, std::vector strategy_idx, std::vector patient_idx, Rcpp::Environment R_StateVal, double dr, std::string type);
+RcppExport SEXP _hesim_C_indiv_ctstm_trans(SEXP R_CtstmTransSEXP, SEXP R_disease_progSEXP, SEXP strategy_idxSEXP, SEXP patient_idxSEXP, SEXP R_StateValSEXP, SEXP drSEXP, SEXP typeSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject rcpp_result_gen;
+ Rcpp::RNGScope rcpp_rngScope_gen;
+ Rcpp::traits::input_parameter< Rcpp::Environment >::type R_CtstmTrans(R_CtstmTransSEXP);
+ Rcpp::traits::input_parameter< Rcpp::DataFrame >::type R_disease_prog(R_disease_progSEXP);
+ Rcpp::traits::input_parameter< std::vector >::type strategy_idx(strategy_idxSEXP);
+ Rcpp::traits::input_parameter< std::vector >::type patient_idx(patient_idxSEXP);
+ Rcpp::traits::input_parameter< Rcpp::Environment >::type R_StateVal(R_StateValSEXP);
+ Rcpp::traits::input_parameter< double >::type dr(drSEXP);
+ Rcpp::traits::input_parameter< std::string >::type type(typeSEXP);
+ rcpp_result_gen = Rcpp::wrap(C_indiv_ctstm_trans(R_CtstmTrans, R_disease_prog, strategy_idx, patient_idx, R_StateVal, dr, type));
+ return rcpp_result_gen;
+END_RCPP
+}
// C_indiv_ctstm_los
std::vector C_indiv_ctstm_los(Rcpp::DataFrame R_disease_prog, std::vector strategy_idx, std::vector patient_idx, double dr);
RcppExport SEXP _hesim_C_indiv_ctstm_los(SEXP R_disease_progSEXP, SEXP strategy_idxSEXP, SEXP patient_idxSEXP, SEXP drSEXP) {
@@ -685,6 +702,7 @@ static const R_CallMethodDef CallEntries[] = {
{"_hesim_C_ctstm_indiv_stateprobs", (DL_FUNC) &_hesim_C_ctstm_indiv_stateprobs, 11},
{"_hesim_C_indiv_ctstm_wlos", (DL_FUNC) &_hesim_C_indiv_ctstm_wlos, 7},
{"_hesim_C_indiv_ctstm_starting", (DL_FUNC) &_hesim_C_indiv_ctstm_starting, 6},
+ {"_hesim_C_indiv_ctstm_trans", (DL_FUNC) &_hesim_C_indiv_ctstm_trans, 7},
{"_hesim_C_indiv_ctstm_los", (DL_FUNC) &_hesim_C_indiv_ctstm_los, 4},
{"_hesim_C_psm_curves_summary", (DL_FUNC) &_hesim_C_psm_curves_summary, 4},
{"_hesim_C_psm_sim_stateprobs", (DL_FUNC) &_hesim_C_psm_sim_stateprobs, 1},
diff --git a/src/indiv-ctstm.cpp b/src/indiv-ctstm.cpp
index aa3f7358..741b7535 100644
--- a/src/indiv-ctstm.cpp
+++ b/src/indiv-ctstm.cpp
@@ -339,6 +339,68 @@ std::vector C_indiv_ctstm_starting(Rcpp::DataFrame R_disease_prog,
return out;
}
+/***************************************************************************//**
+ * @ingroup ctstm
+ * Transition-specific expected values.
+ * Simulate expected values in an individual patient simulation based on transition-
+ * specific values that occur when a patient follows a transition. Values
+ * accrue each time a patient has a transition and are discounted based on the
+ * transition.
+ * @param R_disease_prog An R object of simulating disease progression generated
+ * using C_ctstm_sim_disease.
+ * @param strategy_idx The strategy index starting at 0.
+ * @param patient_idx The patient index starting at 0.
+ * @param R_StateVal An R object of class @c StateVal.
+ * @param dr The discount rate.
+ * @param type @c predict for mean values or @c random for random samples.
+ * @return A vector of expected values in each row in R_disease_prog. These
+ * values are then summed by @c patient_id using @c data.table at the @c R level
+ * in the private member function @c IndivCtstm$sim_trans.
+ ******************************************************************************/
+// [[Rcpp::export]]
+std::vector C_indiv_ctstm_trans(Rcpp::Environment R_CtstmTrans,
+ Rcpp::DataFrame R_disease_prog,
+ std::vector strategy_idx,
+ std::vector patient_idx,
+ Rcpp::Environment R_StateVal,
+ double dr, std::string type
+ ){
+ std::unique_ptr transmod = hesim::ctstm::transmod::create(R_CtstmTrans);
+ hesim::ctstm::disease_prog disease_prog(R_disease_prog);
+ bool time_reset = Rcpp::as(R_StateVal["time_reset"]);
+ hesim::statmods::obs_index obs_index(hesim::statmods::get_id_object(R_StateVal));
+ hesim::statevals stvals(R_StateVal);
+
+ int N = disease_prog.sample_.size();
+ std::vector out(N);
+ for (unsigned int i = 0; (int) i < N; ++i){
+ std::vector ids = transmod->trans_mat_.trans_id(disease_prog.from_[i]);
+ std::vector tos = transmod->trans_mat_.to(disease_prog.from_[i]);
+ // find the transition index
+ int transition_idx = ids.size(); // default: no match
+ for (int j=0; j<(int)ids.size(); ++j) {
+ if (tos[j] == disease_prog.to_[i]) {
+ transition_idx = ids[j];
+ break;
+ }
+ }
+ if (transition_idx != (int) ids.size()) {
+ double time = time_reset ? (disease_prog.to_[i]-disease_prog.from_[i]) :
+ disease_prog.to_[i];
+ int time_idx = hesim::find_interval(time, obs_index.time_start_);
+ int obs = obs_index(strategy_idx[i],
+ patient_idx[i],
+ transition_idx,
+ time_idx);
+ double yhat = stvals.sim(disease_prog.sample_[i], obs, type);
+ out[i] = exp(-dr * time) * yhat;
+ } else {
+ out[i] = 0.0;
+ }
+ }
+ return out;
+}
+
/***************************************************************************//**
* @ingroup ctstm
* Simulate length of stay given simulated disease progression
diff --git a/tests/testthat/test-cpp-utils.R b/tests/testthat/test-cpp-utils.R
index e2276c7e..2d152e05 100644
--- a/tests/testthat/test-cpp-utils.R
+++ b/tests/testthat/test-cpp-utils.R
@@ -40,7 +40,7 @@ test_that("max_lt", {
})
# Test c++ function hesim_bound -----------------------------------------------------
-test_that("hesim_bound", {
+test_that("find_interval", {
vec = c(0,3,5)
expect_equal(hesim:::C_test_find_interval(-1,vec),0)
expect_equal(hesim:::C_test_find_interval(0,vec),0)
diff --git a/tests/testthat/test-ctstm.R b/tests/testthat/test-ctstm.R
index a86ee53b..06aa8402 100644
--- a/tests/testthat/test-ctstm.R
+++ b/tests/testthat/test-ctstm.R
@@ -552,6 +552,25 @@ test_that("Simulate costs and QALYs", {
}
test_starting(ictstm2)
+
+ ## Using method = "transition"
+ drugcost_tbl2 <- stateval_tbl(tbl = data.frame(state_id = 1:3,
+ est = c(10, 20, 30)),
+ dist = "fixed")
+ drugcostsmod2 <- create_StateVals(drugcost_tbl2, n = n_samples,
+ time_reset = FALSE, method = "transition",
+ hesim_data = hesim_dat)
+ ictstm2 <- ictstm$clone()
+ ictstm2$cost_models <- list(drugs=drugcostsmod2)
+ ictstm2$sim_costs(dr = 0, by_patient = TRUE)
+ test_transition_costs <- function(model){
+ model$disprog_[, costs := ifelse(from==1 & to==2, 10,
+ ifelse(from==1 & to==3, 20, 30))]
+ expect_equal(model$disprog_[, sum(costs), by=from]$V1,
+ model$costs_[, sum(costs), by=state_id]$V1)
+ }
+ test_transition_costs(ictstm2)
+
# Summarize costs and QALYs
## By patient = TRUE
ce_summary <- ictstm$summarize()