Transition costs only;)

.Rbuildignore

@@ -23,3 +23,4 @@ vignettes/benchmarks*
 ^vignettes/precompile.R
 ^CRAN-SUBMISSION$
 ^Makefile$
+^CODEOWNERS$

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

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)
 }

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){

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

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

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
 }
 

src/RcppExports.cpp

@@ -261,6 +261,23 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// C_indiv_ctstm_trans
+std::vector<double> C_indiv_ctstm_trans(Rcpp::Environment R_CtstmTrans, Rcpp::DataFrame R_disease_prog, std::vector<int> strategy_idx, std::vector<int> 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<int> >::type strategy_idx(strategy_idxSEXP);
+    Rcpp::traits::input_parameter< std::vector<int> >::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<double> C_indiv_ctstm_los(Rcpp::DataFrame R_disease_prog, std::vector<int> strategy_idx, std::vector<int> 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},

src/indiv-ctstm.cpp

@@ -339,6 +339,68 @@ std::vector<double> 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<double> C_indiv_ctstm_trans(Rcpp::Environment R_CtstmTrans,
+					Rcpp::DataFrame R_disease_prog,
+					std::vector<int> strategy_idx,
+					std::vector<int> patient_idx,
+					Rcpp::Environment R_StateVal,
+					double dr, std::string type
+					){
+  std::unique_ptr<hesim::ctstm::transmod> transmod = hesim::ctstm::transmod::create(R_CtstmTrans);
+  hesim::ctstm::disease_prog disease_prog(R_disease_prog);
+  bool time_reset = Rcpp::as<bool>(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<double> out(N);
+  for (unsigned int i = 0; (int) i < N; ++i){
+    std::vector<int> ids = transmod->trans_mat_.trans_id(disease_prog.from_[i]);
+    std::vector<int> 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 

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)

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()