Skip to content

Instantly share code, notes, and snippets.

@mclements

mclements/code.R

Created July 21, 2024 10:17
Show Gist options
  • Save mclements/ba5ce9964f838c25e07da13d487db0dd to your computer and use it in GitHub Desktop.
Save mclements/ba5ce9964f838c25e07da13d487db0dd to your computer and use it in GitHub Desktop.
Download ZIP
Extension to hesim for cCTSTMs
Raw
code.cpp
/**
Requires the transition_costs branch
Mark Clements
2024-07-08
*/
#include <boost/numeric/odeint.hpp>
#include <RcppArmadillo.h>
#include <hesim.h>
#include <hesim/ctstm/ctstm.h>
#include <hesim/statevals.h>
#include <hesim/dtstm.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(hesim)]]
// [[Rcpp::depends(BH)]]
// Boilerplate code to ensure that arma plays nicely with boost::numeric::odeint
namespace boost {
namespace numeric {
namespace odeint {
template <>
struct is_resizeable<arma::vec>
{
typedef boost::true_type type;
const static bool value = type::value;
};
template <>
struct same_size_impl<arma::vec, arma::vec>
{
static bool same_size(const arma::vec& x, const arma::vec& y)
{
return x.size() == y.size();
}
};
template<>
struct resize_impl<arma::vec, arma::vec>
{
static void resize(arma::vec& v1, const arma::vec& v2)
{
v1.resize(v2.size());
}
};
}
}
} // namespace boost::numeric::odeint
#define CTSTM_OUT(OBJECT,COUNTER) \
OBJECT.sample_[COUNTER] = s; \
OBJECT.strategy_id_[COUNTER] = k; \
OBJECT.patient_id_[COUNTER] = i; \
OBJECT.grp_id_[COUNTER] = transmod->obs_index_.get_grp_id(); \
OBJECT.patient_wt_[COUNTER] = transmod->obs_index_.get_patient_wt(); \
OBJECT.state_id_[COUNTER] = h;
/***************************************************************************//**
* @ingroup ctstm
* Simulate disease progression (i.e., a path through a multi-state model).
* This function is exported to @c R and used in @c CohortCtstmTrans$sim_stateprobs() and
* in @c CohortCtstm$sim_disase().
* @param R_CtstmTrans An R object of class @c CohortCtstmTrans.
* @param start_state The starting health state for each patient and random sample
* of the parameter sets. (TODO: allow for probabilities in the initial states.)
* @param start_age The starting age of each patient in the simulation.
* @param start_time The starting time of the simulation.
* @param max_t The maximum time to simulate the model for.
* @param max_age The maximum age that a patient can live to.
* @return An R data frame of the same format as stateprobs_out.
******************************************************************************/
// [[Rcpp::export]]
Rcpp::List C_cohort_ctstm_sim(Rcpp::Environment R_CtstmTrans,
Rcpp::List R_CostsStateVals,
Rcpp::Environment R_QALYsStateVal,
std::vector<int> start_state, // by patient
std::vector<double> start_age, // by patient
std::vector<double> times, // common
std::string clock,
std::vector<int> transition_types,
int progress = 0,
double dr_qalys = 0.0,
double dr_costs = 0.0,
std::string type="predict",
double eps = 1e-100) {
using namespace boost::numeric::odeint;
typedef arma::vec state_type;
enum TransitionType {tt_time, tt_age};
// Initialize
std::unique_ptr<hesim::ctstm::transmod> transmod = hesim::ctstm::transmod::create(R_CtstmTrans);
int n_costs = R_CostsStateVals.size();
std::vector<hesim::statevals> costs_lookup;
std::vector<hesim::statmods::obs_index> obs_index_costs;
for (int cost_=0; cost_<n_costs; ++cost_) {
costs_lookup.push_back(hesim::statevals(Rcpp::as<Rcpp::Environment>(R_CostsStateVals[cost_])));
obs_index_costs.push_back(hesim::statmods::obs_index(hesim::statmods::get_id_object(Rcpp::as<Rcpp::Environment>(R_CostsStateVals[cost_]))));
}
// use unique_ptr because statevals and obs_index can be null and do not have null constructors
std::unique_ptr<hesim::statevals> qalys_lookup;
std::unique_ptr<hesim::statmods::obs_index> obs_index_qalys;
bool valid_R_QALYsStateVal = R_QALYsStateVal.exists("method");
if (valid_R_QALYsStateVal) {
qalys_lookup = std::unique_ptr<hesim::statevals>(new hesim::statevals(R_QALYsStateVal));
auto id_object = hesim::statmods::get_id_object(R_QALYsStateVal);
obs_index_qalys =
std::unique_ptr<hesim::statmods::obs_index>(new hesim::statmods::obs_index(id_object));
}
int n_samples = transmod->get_n_samples();
int n_strategies = transmod->get_n_strategies();
int n_patients = transmod->get_n_patients(); // actually, the number of cohorts
int n_states = transmod->trans_mat_.n_states_;
int n_times = times.size();
int N = n_samples *
n_strategies *
n_patients *
n_states*
n_times;
hesim::stateprobs_out out(N);
int N2 = n_samples *
n_strategies *
n_patients *
n_states *
(1+(valid_R_QALYsStateVal ? 1 : 0)+n_costs); // LYs, QALYs, cost categories
hesim::ev_out out2(N2);
int counter = 0, counter2 = 0;
// Loop
for (int s = 0; s < n_samples; ++s){
if (progress > 0){
if ((s + 1) % progress == 0){ // R-based indexing
Rcpp::Rcout << "sample = " << s + 1 << std::endl;
}
}
for (int k = 0; k < n_strategies; ++k){
transmod->obs_index_.set_strategy_index(k);
for (int i = 0; i < n_patients; ++i){
transmod->obs_index_.set_patient_index(i);
arma::vec p0 = arma::zeros(n_states*(3+R_CostsStateVals.size()));
p0[start_state[i]] = 1.0;
int n_trans = transmod->trans_mat_.n_trans_;
arma::uvec from(n_trans);
arma::uvec to(n_trans);
for (int state_ = 0; state_<n_states; ++state_) {
std::vector<int> trans_ids = transmod->trans_mat_.trans_id(state_);
std::vector<int> tos = transmod->trans_mat_.to(state_);
int n_trans_state = trans_ids.size();
for (int trans_ = 0; trans_ < n_trans_state; ++trans_){ // NB: within a state
int trans_id = trans_ids[trans_];
from(trans_id) = state_;
to(trans_id) = tos[trans_];
}
}
arma::mat report(n_times,p0.size());
report.row(0) = p0.t();
auto stepper = make_dense_output(1.0e-10, 1.0e-10, runge_kutta_dopri5<state_type>());
for (size_t j=1; j<n_times; j++) {
size_t n = integrate_adaptive(stepper,
[&](const state_type &Y , state_type &dYdt, const double t)
{
arma::vec rates(n_trans);
arma::vec utilities(n_states);
arma::mat costs_by_state(n_states,n_costs);
arma::mat costs_by_transition(n_trans,n_costs);
double age = start_age[i]+t;
double tstar = std::max(eps,t);
// rate calculations
for (int state_ = 0; state_<n_states; ++state_) {
for (int trans_ = 0; trans_ < n_trans; ++trans_){
if (clock == "forward"){
rates[trans_] = transmod->summary(trans_, s, {tstar}, "hazard")[0];
} else { // clock == "mixt"
if (transition_types[trans_] == tt_age){
rates[trans_] = transmod->summary(trans_, s, {age}, "hazard")[0];
} else { // transition_types[trans_] == tt_time (tt_reset is *not* currently available)
rates[trans_] = transmod->summary(trans_, s, {tstar}, "hazard")[0];
}
}
} // end loop over transitions
} // end loop over states
// utility input calculations
if (valid_R_QALYsStateVal) {
int t_index_qalys = hesim::hesim_bound(t, obs_index_qalys->time_start_);
for (int state_ = 0; state_<n_states; ++state_) {
int obs_qalys = (*obs_index_qalys)(k, // strategy
i, // patient
state_,
t_index_qalys);
utilities(state_) = qalys_lookup->sim(s, obs_qalys, type);
} // end loop over states
}
// cost input calculations
for (int cost_=0; cost_<n_costs; ++cost_) {
int t_index_costs = hesim::hesim_bound(t, obs_index_costs[cost_].time_start_);
if (costs_lookup[cost_].method_ == "wlos") {
for (int state_ = 0; state_<n_states; ++state_) {
int obs_costs = obs_index_costs[cost_](k, // strategy
i, // patient
state_,
t_index_costs);
costs_by_state(state_, cost_) += costs_lookup[cost_].sim(s, obs_costs, type);
} // end loop over states
} else if (costs_lookup[cost_].method_ == "starting") {
for (int trans_=0; trans_<n_trans; ++trans_) {
int obs_costs = obs_index_costs[cost_](k, // strategy
i, // patient
to(trans_),
t_index_costs);
costs_by_state(to(trans_), cost_) = costs_lookup[cost_].sim(s, obs_costs, type);
} // end loop over transitions
} else { // costs_lookup[cost_].method_ == "transition"
for (int trans_=0; trans_<n_trans; ++trans_) {
int obs_costs = obs_index_costs[cost_](k, // strategy
i, // patient
trans_, // assumes transition
t_index_costs);
costs_by_transition(trans_, cost_) = costs_lookup[cost_].sim(s, obs_costs, type); // by trans_
} // end loop over transitions
} // end case: "transition"
} // end loop over costs
dYdt = dYdt*0.0;
arma::vec delta = Y(from) % rates;
dYdt(to) += delta;
dYdt(from) -= delta;
double drr_utilities = std::exp(-dr_qalys*t);
double drr_costs = std::exp(-dr_costs*t);
// undiscounted life-years
dYdt(arma::span(n_states,2*n_states-1)) = Y(arma::span(0,n_states-1));
// discounted utilities
dYdt(arma::span(2*n_states,3*n_states-1)) = utilities % Y(arma::span(0,n_states-1))*drr_utilities;
// discounted costs
for (int cost_=0; cost_<n_costs; ++cost_) {
if (costs_lookup[cost_].method_ == "wlos")
dYdt(arma::span((3+cost_)*n_states,(4+cost_)*n_states-1)) += costs_by_state.col(cost_) % Y(arma::span(0,n_states-1))*drr_costs;
else if (costs_lookup[cost_].method_ == "starting") {
for (int trans_=0; trans_<n_trans; ++trans_)
dYdt(to[trans_] + (3+cost_)*n_states) += costs_by_state(to[trans_],cost_) * Y(from[trans_]) * rates(trans_) *drr_costs;
} else { // costs_lookup[cost_].method_ == "transition"
for (int trans_=0; trans_<n_trans; ++trans_)
dYdt(from[trans_] + (3+cost_)*n_states) += // costs arbitrarily assigned to the state of origin (from state)
costs_by_transition(trans_,cost_) * Y(from(trans_)) * rates(trans_) * drr_costs;
}
}
},
p0,
times[j-1],
times[j],
times[j]-times[j-1]);
report.row(j) = p0.t();
} // end j times_ loop
// Output
for (int h = 0; h < n_states; ++h){
for (int ti = 0; ti < n_times; ++ti){
CTSTM_OUT(out,counter);
out.t_[counter] = times[ti];
out.prob_[counter] = report(ti, h);
++counter;
} // end cycles loop
// report for life-years
CTSTM_OUT(out2,counter2);
out2.dr_[counter2] = 0.0; // assume no discounting for life-years
out2.outcome_[counter2] = "ly";
out2.value_[counter2] = report(n_times-1,n_states+h);
++counter2;
// report for qalys
if (valid_R_QALYsStateVal) {
CTSTM_OUT(out2,counter2);
out2.dr_[counter2] = dr_qalys;
out2.outcome_[counter2] = "qaly";
out2.value_[counter2] = report(n_times-1,2*n_states+h);
++counter2;
}
// report for costs
for (int cost_=0; cost_<n_costs; ++cost_) {
CTSTM_OUT(out2,counter2);
out2.dr_[counter2] = dr_costs;
out2.outcome_[counter2] = "Category " + std::to_string(cost_+1);
out2.value_[counter2] = report(n_times-1,(3+cost_)*n_states+h);
++counter2;
} // end cost category loop
} // end state loop
} // end patient loop
} // end strategy loop
} // end parameter sampling loop
// Return
using namespace Rcpp;
return(List::create(_["stateprobs"]=out.create_R_data_frame(),
_["ev"]=out2.create_R_data_frame()));
}
#undef CTSTM_OUT
// test example that does not use hesim
class TestODE {
public:
typedef arma::vec state_type;
size_t n_states;
double discount_rate;
TestODE(double discount_rate = 0.0)
: n_states(4), discount_rate(discount_rate) { }
inline double hweibull(const double t, const double shape, const double scale) {
return R::dweibull(std::max(1e-100,t),shape,scale,0)/R::pweibull(t,shape,scale,0,0);
}
void operator() (const state_type &Y , state_type &dYdt, const double t)
{
run_step(Y, dYdt, t);
}
void run_step(const state_type &Y , state_type &dYdt, const double t)
{
using namespace arma;
// conv_to<uvec>::from(an_ivec)
// vec rates = {0.1, 0.2, 0.3};
vec rates = {hweibull(t,0.8,1.0), hweibull(t,0.8,1.0),
hweibull(t,0.8,2.0), hweibull(t,0.8,3.0)};
vec utilities = {0.9, 0.8, 0.7, 0.0};
vec costs_per_unit_time = {10000.0, 20000.0, 30000.0, 0.0}; // by state
vec costs_per_transition = {10.0, 15.0, 20.0, 30.0}; // by transition
vec starting_costs = {0.0, 2000.0, 3000.0, 0.0}; // by state
uvec from = {0, 1, 1, 2};
uvec to = {1, 0, 2, 3};
dYdt = dYdt*0.0;
// state occupation/transition probabilities
vec delta = Y(from) % rates;
dYdt(to) += delta;
dYdt(from) -= delta;
// double dr = 1.0/pow(1.0+discount_rate,t);
double dr = std::exp(-discount_rate*t);
// undiscounted life-years
dYdt(span(n_states,2*n_states-1)) = Y(span(0,n_states-1));
// discounted utilities
dYdt(span(2*n_states,3*n_states-1)) = utilities % Y(span(0,n_states-1))*dr;
// discounted costs
dYdt(span(3*n_states,4*n_states-1)) = costs_per_unit_time % Y(span(0,n_states-1))*dr;
// discounted transition costs
for (int j=0; j<4; ++j)
dYdt(from[j] + 4*n_states) += costs_per_transition[j] * Y(from[j]) * rates(j) * dr;
// discounted starting costs
for (int j=0; j<4; ++j)
dYdt(to[j] + 5*n_states) += starting_costs[to[j]] * Y(from[j]) * rates(j) * dr;
}
arma::mat run(arma::vec p0, arma::vec times) {
using namespace boost::numeric::odeint;
size_t n_times = times.size();
// combine the results
arma::mat combined(n_times,p0.size());
combined.row(0) = p0.t();
auto stepper = make_dense_output(1.0e-10, 1.0e-10, runge_kutta_dopri5<state_type>());
for (size_t i=1; i<n_times; i++) {
size_t n = integrate_adaptive(stepper,
[this](const state_type &Y , state_type &dYdt, const double t) {
this->run_step(Y,dYdt,t);
},
p0,
times[i-1],
times[i],
times[i]-times[i-1]);
combined.row(i) = p0.t();
}
return combined;
}
};
// [[Rcpp::export]]
Rcpp::List runTestODE(arma::vec p0, arma::vec times, double discount_rate = 0.0) {
using namespace Rcpp;
return List::create(_("times")=times,
_("Y")=TestODE(discount_rate).run(p0,times));
}
Raw
code.R
library("RcppArmadillo")
library("Rcpp")
library("hesim")
library("flexsurv")
sourceCpp("code/**
Requires the transition_costs branch
Mark Clements
2024-07-08
*/
#include <boost/numeric/odeint.hpp>
#include <RcppArmadillo.h>
#include <hesim.h>
#include <hesim/ctstm/ctstm.h>
#include <hesim/statevals.h>
#include <hesim/dtstm.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(hesim)]]
// [[Rcpp::depends(BH)]]
// Boilerplate code to ensure that arma plays nicely with boost::numeric::odeint
namespace boost {
namespace numeric {
namespace odeint {
template <>
struct is_resizeable<arma::vec>
{
typedef boost::true_type type;
const static bool value = type::value;
};
template <>
struct same_size_impl<arma::vec, arma::vec>
{
static bool same_size(const arma::vec& x, const arma::vec& y)
{
return x.size() == y.size();
}
};
template<>
struct resize_impl<arma::vec, arma::vec>
{
static void resize(arma::vec& v1, const arma::vec& v2)
{
v1.resize(v2.size());
}
};
}
}
} // namespace boost::numeric::odeint
#define CTSTM_OUT(OBJECT,COUNTER) \
OBJECT.sample_[COUNTER] = s; \
OBJECT.strategy_id_[COUNTER] = k; \
OBJECT.patient_id_[COUNTER] = i; \
OBJECT.grp_id_[COUNTER] = transmod->obs_index_.get_grp_id(); \
OBJECT.patient_wt_[COUNTER] = transmod->obs_index_.get_patient_wt(); \
OBJECT.state_id_[COUNTER] = h;
/***************************************************************************//**
* @ingroup ctstm
* Simulate disease progression (i.e., a path through a multi-state model).
* This function is exported to @c R and used in @c CohortCtstmTrans$sim_stateprobs() and
* in @c CohortCtstm$sim_disase().
* @param R_CtstmTrans An R object of class @c CohortCtstmTrans.
* @param start_state The starting health state for each patient and random sample
* of the parameter sets. (TODO: allow for probabilities in the initial states.)
* @param start_age The starting age of each patient in the simulation.
* @param start_time The starting time of the simulation.
* @param max_t The maximum time to simulate the model for.
* @param max_age The maximum age that a patient can live to.
* @return An R data frame of the same format as stateprobs_out.
******************************************************************************/
// [[Rcpp::export]]
Rcpp::List C_cohort_ctstm_sim(Rcpp::Environment R_CtstmTrans,
Rcpp::List R_CostsStateVals,
Rcpp::Environment R_QALYsStateVal,
std::vector<int> start_state, // by patient
std::vector<double> start_age, // by patient
std::vector<double> times, // common
std::string clock,
std::vector<int> transition_types,
int progress = 0,
double dr_qalys = 0.0,
double dr_costs = 0.0,
std::string type="predict",
double eps = 1e-100) {
using namespace boost::numeric::odeint;
typedef arma::vec state_type;
enum TransitionType {tt_time, tt_age};
// Initialize
std::unique_ptr<hesim::ctstm::transmod> transmod = hesim::ctstm::transmod::create(R_CtstmTrans);
int n_costs = R_CostsStateVals.size();
std::vector<hesim::statevals> costs_lookup;
std::vector<hesim::statmods::obs_index> obs_index_costs;
for (int cost_=0; cost_<n_costs; ++cost_) {
costs_lookup.push_back(hesim::statevals(Rcpp::as<Rcpp::Environment>(R_CostsStateVals[cost_])));
obs_index_costs.push_back(hesim::statmods::obs_index(hesim::statmods::get_id_object(Rcpp::as<Rcpp::Environment>(R_CostsStateVals[cost_]))));
}
// use unique_ptr because statevals and obs_index can be null and do not have null constructors
std::unique_ptr<hesim::statevals> qalys_lookup;
std::unique_ptr<hesim::statmods::obs_index> obs_index_qalys;
bool valid_R_QALYsStateVal = R_QALYsStateVal.exists("method");
if (valid_R_QALYsStateVal) {
qalys_lookup = std::unique_ptr<hesim::statevals>(new hesim::statevals(R_QALYsStateVal));
auto id_object = hesim::statmods::get_id_object(R_QALYsStateVal);
obs_index_qalys =
std::unique_ptr<hesim::statmods::obs_index>(new hesim::statmods::obs_index(id_object));
}
int n_samples = transmod->get_n_samples();
int n_strategies = transmod->get_n_strategies();
int n_patients = transmod->get_n_patients(); // actually, the number of cohorts
int n_states = transmod->trans_mat_.n_states_;
int n_times = times.size();
int N = n_samples *
n_strategies *
n_patients *
n_states*
n_times;
hesim::stateprobs_out out(N);
int N2 = n_samples *
n_strategies *
n_patients *
n_states *
(1+(valid_R_QALYsStateVal ? 1 : 0)+n_costs); // LYs, QALYs, cost categories
hesim::ev_out out2(N2);
int counter = 0, counter2 = 0;
// Loop
for (int s = 0; s < n_samples; ++s){
if (progress > 0){
if ((s + 1) % progress == 0){ // R-based indexing
Rcpp::Rcout << "sample = " << s + 1 << std::endl;
}
}
for (int k = 0; k < n_strategies; ++k){
transmod->obs_index_.set_strategy_index(k);
for (int i = 0; i < n_patients; ++i){
transmod->obs_index_.set_patient_index(i);
arma::vec p0 = arma::zeros(n_states*(3+R_CostsStateVals.size()));
p0[start_state[i]] = 1.0;
int n_trans = transmod->trans_mat_.n_trans_;
arma::uvec from(n_trans);
arma::uvec to(n_trans);
for (int state_ = 0; state_<n_states; ++state_) {
std::vector<int> trans_ids = transmod->trans_mat_.trans_id(state_);
std::vector<int> tos = transmod->trans_mat_.to(state_);
int n_trans_state = trans_ids.size();
for (int trans_ = 0; trans_ < n_trans_state; ++trans_){ // NB: within a state
int trans_id = trans_ids[trans_];
from(trans_id) = state_;
to(trans_id) = tos[trans_];
}
}
arma::mat report(n_times,p0.size());
report.row(0) = p0.t();
auto stepper = make_dense_output(1.0e-10, 1.0e-10, runge_kutta_dopri5<state_type>());
for (size_t j=1; j<n_times; j++) {
size_t n = integrate_adaptive(stepper,
[&](const state_type &Y , state_type &dYdt, const double t)
{
arma::vec rates(n_trans);
arma::vec utilities(n_states);
arma::mat costs_by_state(n_states,n_costs);
arma::mat costs_by_transition(n_trans,n_costs);
double age = start_age[i]+t;
double tstar = std::max(eps,t);
// rate calculations
for (int state_ = 0; state_<n_states; ++state_) {
for (int trans_ = 0; trans_ < n_trans; ++trans_){
if (clock == "forward"){
rates[trans_] = transmod->summary(trans_, s, {tstar}, "hazard")[0];
} else { // clock == "mixt"
if (transition_types[trans_] == tt_age){
rates[trans_] = transmod->summary(trans_, s, {age}, "hazard")[0];
} else { // transition_types[trans_] == tt_time (tt_reset is *not* currently available)
rates[trans_] = transmod->summary(trans_, s, {tstar}, "hazard")[0];
}
}
} // end loop over transitions
} // end loop over states
// utility input calculations
if (valid_R_QALYsStateVal) {
int t_index_qalys = hesim::hesim_bound(t, obs_index_qalys->time_start_);
for (int state_ = 0; state_<n_states; ++state_) {
int obs_qalys = (*obs_index_qalys)(k, // strategy
i, // patient
state_,
t_index_qalys);
utilities(state_) = qalys_lookup->sim(s, obs_qalys, type);
} // end loop over states
}
// cost input calculations
for (int cost_=0; cost_<n_costs; ++cost_) {
int t_index_costs = hesim::hesim_bound(t, obs_index_costs[cost_].time_start_);
if (costs_lookup[cost_].method_ == "wlos") {
for (int state_ = 0; state_<n_states; ++state_) {
int obs_costs = obs_index_costs[cost_](k, // strategy
i, // patient
state_,
t_index_costs);
costs_by_state(state_, cost_) += costs_lookup[cost_].sim(s, obs_costs, type);
} // end loop over states
} else if (costs_lookup[cost_].method_ == "starting") {
for (int trans_=0; trans_<n_trans; ++trans_) {
int obs_costs = obs_index_costs[cost_](k, // strategy
i, // patient
to(trans_),
t_index_costs);
costs_by_state(to(trans_), cost_) = costs_lookup[cost_].sim(s, obs_costs, type);
} // end loop over transitions
} else { // costs_lookup[cost_].method_ == "transition"
for (int trans_=0; trans_<n_trans; ++trans_) {
int obs_costs = obs_index_costs[cost_](k, // strategy
i, // patient
trans_, // assumes transition
t_index_costs);
costs_by_transition(trans_, cost_) = costs_lookup[cost_].sim(s, obs_costs, type); // by trans_
} // end loop over transitions
} // end case: "transition"
} // end loop over costs
dYdt = dYdt*0.0;
arma::vec delta = Y(from) % rates;
dYdt(to) += delta;
dYdt(from) -= delta;
double drr_utilities = std::exp(-dr_qalys*t);
double drr_costs = std::exp(-dr_costs*t);
// undiscounted life-years
dYdt(arma::span(n_states,2*n_states-1)) = Y(arma::span(0,n_states-1));
// discounted utilities
dYdt(arma::span(2*n_states,3*n_states-1)) = utilities % Y(arma::span(0,n_states-1))*drr_utilities;
// discounted costs
for (int cost_=0; cost_<n_costs; ++cost_) {
if (costs_lookup[cost_].method_ == "wlos")
dYdt(arma::span((3+cost_)*n_states,(4+cost_)*n_states-1)) += costs_by_state.col(cost_) % Y(arma::span(0,n_states-1))*drr_costs;
else if (costs_lookup[cost_].method_ == "starting") {
for (int trans_=0; trans_<n_trans; ++trans_)
dYdt(to[trans_] + (3+cost_)*n_states) += costs_by_state(to[trans_],cost_) * Y(from[trans_]) * rates(trans_) *drr_costs;
} else { // costs_lookup[cost_].method_ == "transition"
for (int trans_=0; trans_<n_trans; ++trans_)
dYdt(from[trans_] + (3+cost_)*n_states) += // costs arbitrarily assigned to the state of origin (from state)
costs_by_transition(trans_,cost_) * Y(from(trans_)) * rates(trans_) * drr_costs;
}
}
},
p0,
times[j-1],
times[j],
times[j]-times[j-1]);
report.row(j) = p0.t();
} // end j times_ loop
// Output
for (int h = 0; h < n_states; ++h){
for (int ti = 0; ti < n_times; ++ti){
CTSTM_OUT(out,counter);
out.t_[counter] = times[ti];
out.prob_[counter] = report(ti, h);
++counter;
} // end cycles loop
// report for life-years
CTSTM_OUT(out2,counter2);
out2.dr_[counter2] = 0.0; // assume no discounting for life-years
out2.outcome_[counter2] = "ly";
out2.value_[counter2] = report(n_times-1,n_states+h);
++counter2;
// report for qalys
if (valid_R_QALYsStateVal) {
CTSTM_OUT(out2,counter2);
out2.dr_[counter2] = dr_qalys;
out2.outcome_[counter2] = "qaly";
out2.value_[counter2] = report(n_times-1,2*n_states+h);
++counter2;
}
// report for costs
for (int cost_=0; cost_<n_costs; ++cost_) {
CTSTM_OUT(out2,counter2);
out2.dr_[counter2] = dr_costs;
out2.outcome_[counter2] = "Category " + std::to_string(cost_+1);
out2.value_[counter2] = report(n_times-1,(3+cost_)*n_states+h);
++counter2;
} // end cost category loop
} // end state loop
} // end patient loop
} // end strategy loop
} // end parameter sampling loop
// Return
using namespace Rcpp;
return(List::create(_["stateprobs"]=out.create_R_data_frame(),
_["ev"]=out2.create_R_data_frame()));
}
#undef CTSTM_OUT
// test example that does not use hesim
class TestODE {
public:
typedef arma::vec state_type;
size_t n_states;
double discount_rate;
TestODE(double discount_rate = 0.0)
: n_states(4), discount_rate(discount_rate) { }
inline double hweibull(const double t, const double shape, const double scale) {
return R::dweibull(std::max(1e-100,t),shape,scale,0)/R::pweibull(t,shape,scale,0,0);
}
void operator() (const state_type &Y , state_type &dYdt, const double t)
{
run_step(Y, dYdt, t);
}
void run_step(const state_type &Y , state_type &dYdt, const double t)
{
using namespace arma;
// conv_to<uvec>::from(an_ivec)
// vec rates = {0.1, 0.2, 0.3};
vec rates = {hweibull(t,0.8,1.0), hweibull(t,0.8,1.0),
hweibull(t,0.8,2.0), hweibull(t,0.8,3.0)};
vec utilities = {0.9, 0.8, 0.7, 0.0};
vec costs_per_unit_time = {10000.0, 20000.0, 30000.0, 0.0}; // by state
vec costs_per_transition = {10.0, 15.0, 20.0, 30.0}; // by transition
vec starting_costs = {0.0, 2000.0, 3000.0, 0.0}; // by state
uvec from = {0, 1, 1, 2};
uvec to = {1, 0, 2, 3};
dYdt = dYdt*0.0;
// state occupation/transition probabilities
vec delta = Y(from) % rates;
dYdt(to) += delta;
dYdt(from) -= delta;
// double dr = 1.0/pow(1.0+discount_rate,t);
double dr = std::exp(-discount_rate*t);
// undiscounted life-years
dYdt(span(n_states,2*n_states-1)) = Y(span(0,n_states-1));
// discounted utilities
dYdt(span(2*n_states,3*n_states-1)) = utilities % Y(span(0,n_states-1))*dr;
// discounted costs
dYdt(span(3*n_states,4*n_states-1)) = costs_per_unit_time % Y(span(0,n_states-1))*dr;
// discounted transition costs
for (int j=0; j<4; ++j)
dYdt(from[j] + 4*n_states) += costs_per_transition[j] * Y(from[j]) * rates(j) * dr;
// discounted starting costs
for (int j=0; j<4; ++j)
dYdt(to[j] + 5*n_states) += starting_costs[to[j]] * Y(from[j]) * rates(j) * dr;
}
arma::mat run(arma::vec p0, arma::vec times) {
using namespace boost::numeric::odeint;
size_t n_times = times.size();
// combine the results
arma::mat combined(n_times,p0.size());
combined.row(0) = p0.t();
auto stepper = make_dense_output(1.0e-10, 1.0e-10, runge_kutta_dopri5<state_type>());
for (size_t i=1; i<n_times; i++) {
size_t n = integrate_adaptive(stepper,
[this](const state_type &Y , state_type &dYdt, const double t) {
this->run_step(Y,dYdt,t);
},
p0,
times[i-1],
times[i],
times[i]-times[i-1]);
combined.row(i) = p0.t();
}
return combined;
}
};
// [[Rcpp::export]]
Rcpp::List runTestODE(arma::vec p0, arma::vec times, double discount_rate = 0.0) {
using namespace Rcpp;
return List::create(_("times")=times,
_("Y")=TestODE(discount_rate).run(p0,times));
}
.cpp")
## Test code
## Simulation data
strategies <- data.frame(strategy_id = 1L)
patients <- data.frame(patient_id = 1:2, intercept=1)
## Multi-state model with transition specific models (as per testODE in the C++ code)
tmat <- rbind(c(NA, 1, NA, NA),
c(2, NA, 3, NA),
c(NA, NA, NA, 4),
c(NA, NA, NA, NA))
## short utility function
make_param <- function(shape,scale,dist="weibull")
params_surv(
coefs = list(
shape = data.frame(
intercept = log(shape)
),
scale = data.frame(
intercept = log(scale)
)
),
dist = dist)
fits = params_surv_list(make_param(0.8, 1.0), make_param(0.8, 1.0),
make_param(0.8, 2.0), make_param(0.8, 3.0))
## Simulation model
hesim_dat <- hesim_data(strategies = strategies,
patients = patients)
fits_data <- expand(hesim_dat)
## create_IndivCtstmTrans can also be used for cCTSTM:)
transmod <- create_IndivCtstmTrans(fits, input_data = fits_data,
trans_mat = tmat,
start_age=50,
start_state=1,
clock="mixt",
transition_types=c("time","time","time","time"))
## utilities
utilities_tbl = stateval_tbl(data.frame(state_id=1:4,
est=c(0.9, 0.8, 0.7, 0.0)),
dist="fixed")
utilities = create_StateVals(utilities_tbl, hesim_data=hesim_dat, n=1)
## costs
costs_wlos_tbl = stateval_tbl(data.frame(state_id=1:4,
est=c(10000.0, 20000.0, 30000.0, 0.0)),
dist="fixed")
costs_wlos = create_StateVals(costs_wlos_tbl, hesim_data=hesim_dat, n=1)
## for the transition costs, "state_id" is actually the id for the transition
costs_transition_tbl <- stateval_tbl(tbl = data.frame(state_id = 1:4,
est = c(10, 15, 20, 30)),
dist = "fixed")
costs_transition <- create_StateVals(costs_transition_tbl, n = 1,
time_reset = FALSE, method = "transition",
hesim_data = hesim_dat)
costs_starting_tbl = stateval_tbl(data.frame(state_id=1:4,
est=c(0.0, 2000.0, 3000.0, 0.0)),
dist="fixed")
costs_starting = create_StateVals(costs_starting_tbl, hesim_data=hesim_dat, n=1,
method="starting")
test = C_cohort_ctstm_sim(transmod,
list(costs_wlos,costs_transition,costs_starting),
utilities,
start_state=transmod$start_state - 1L,
start_age=transmod$start_age,
times=c(0,1,2,10),
clock=transmod$clock,
transition_types=match(transmod$transition_types, c("time","age")) - 1L,
progress=0,
dr_qalys=0.03,
dr_costs=0.03)
local({
names = c("ly","qaly","Category 1","Category 2","Category 3")
test2 = runTestODE(c(1,0,0,0, rep(0,4*5)), c(0,1,2,10), 0.03)
Y = test2$Y |> "dim<-"(c(4,4,6))
(Y[4,,-1] - xtabs(value~state_id+outcome,test$ev,patient_id==1)[,names]) |> abs() |> max()
})
## runTestODE(c(1,0,0,0, rep(0,4*5)), c(0,1,2,10), 0.03)
# CohortCtstm class -------------------------------------------------------------
#' Individual-level continuous time state transition model
#'
#' @description
#' Simulate outcomes from a cohort-level continuous time state transition
#' model (CTSTM). The class supports "clock-forward" (i.e., Markov), and mixtures of
#' clock-age and clock-forward multi-state models as described in
#' [`IndivCtstmTrans`].
#' @format An [R6::R6Class] object.
#' @seealso The [`IndivCtstmTrans`] documentation
#' describes the class for the transition model and the [`StateVals`] documentation
#' describes the class for the cost and utility models. An [`IndivCtstmTrans`]
#' object is typically created using [create_IndivCtstmTrans()].
#'
#' There are currently no relevant vignettes.
#' @name CohortCtstm
NULL
#' @param dr Discount rate.
#' @param type `"predict"` for mean values or `"random"` for random samples
#' as in `$sim()` in [`StateVals`].
#' @export
CohortCtstm <- R6::R6Class("CohortCtstm",
public = list(
#' @field trans_model The model for health state transitions. Must be an object
#' of class [`IndivCtstmTrans`].
trans_model = NULL,
#' @field utility_model The model for health state utility. Must be an object of
#' class [`StateVals`].
utility_model = NULL,
#' @field cost_models The models used to predict costs by health state.
#' Must be a list of objects of class [`StateVals`], where each element of the
#' list represents a different cost category.
cost_models = NULL,
#' @field stateprobs_ An object of class [`stateprobs`] simulated using `$sim_stateprobs()`.
stateprobs_ = NULL,
#' @field qalys_ An object of class [`qalys`] simulated using `$sim_qalys()`.
qalys_ = NULL,
#' @field costs_ An object of class [`costs`] simulated using `$sim_costs()`.
costs_ = NULL,
#' @field times_ A vector of times.
times_ = NULL,
#' @description
#' Create a new `CohortCtstm` object.
#' @param trans_model The `trans_model` field.
#' @param utility_model The `utility_model` field.
#' @param cost_models The `cost_models` field.
#' @return A new `CohortCtstm` object.
initialize = function(trans_model = NULL, utility_model = NULL, cost_models = NULL) {
self$trans_model <- trans_model
self$utility_model <- if (is.null(utility_model)) new.env() else utility_model
self$cost_models <- cost_models
},
#' @description
#' Simulate quality-adjusted life-years (QALYs) as a function of
#' `utility_model` and costs as a function of `cost_models`.
#' @param lys If `TRUE`, then life-years are simulated in addition to QALYs.
#' @return An instance of `self` with simulated output of
#' class [qalys] stored in `qalys_` and class [costs] stored in `costs_`..
sim = function(t=c(0,1), dr_qalys = .03, dr_costs=.03, type = c("predict", "random"),
lys = TRUE, progress=NULL,
by_patient = TRUE,
eps=1e-100){
if(!inherits(self$utility_model, "StateVals")){
stop("'utility_model' must be an object of class 'StateVals'",
call. = FALSE)
}
if(!is.list(self$cost_models)){
stop("'cost_models' must be a list of objects of class 'StateVals'",
call. = FALSE)
} else{
for (i in 1:length(self$cost_models)){
if (!inherits(self$cost_models[[i]], "StateVals")){
stop("'cost_models' must be a list of objects of class 'StateVals'",
call. = FALSE)
}
}
}
type <- match.arg(type)
hesim:::check_dr(dr_qalys)
hesim:::check_dr(dr_costs)
self$times_ <- t
## browser()
C_ev <- C_cohort_ctstm_sim(R_CtstmTrans=self$trans_model,
R_CostsStateVals=self$cost_models,
R_QALYsStateVal=self$utility_model,
start_state=self$trans_model$start_state - 1L,
start_age=self$trans_model$start_age,
times=self$times_,
clock=self$trans_model$clock,
transition_types=match(self$trans_model$transition_types,
c("time","age")) - 1L,
progress=if(is.null(progress)) 0 else progress,
dr_qalys=dr_qalys,
dr_costs=dr_costs,
type=match.arg(type),
eps=eps)
self$stateprobs_ <- as.data.table(C_ev$stateprobs)
ev <- as.data.table(C_ev$ev)
self$qalys_ <- ev[outcome=="qaly"][
,`:=`(category='qalys',qalys=value,value=NULL,outcome=NULL)]
if (lys)
self$qalys_$ly = ev[outcome=="ly","value"]
self$costs_ <- do.call(rbind,
lapply(unique(grep("Category",ev$outcome,value=TRUE)),
function(pattern) {
ev2 <- ev[outcome==pattern][
,`:=`(category=pattern,costs=value,value=NULL,outcome=NULL)]
ev2
}))
invisible(self)
},
#' @description
#' Summarize costs and QALYs so that cost-effectiveness analysis can be performed.
#' See [summarize_ce()].
#' @param by_grp If `TRUE`, then costs and QALYs are computed by subgroup. If
#' `FALSE`, then costs and QALYs are aggregated across all patients (and subgroups).
summarize = function(by_grp = FALSE) {
hesim:::check_summarize(self)
hesim:::summarize_ce(self$costs_, self$qalys_, by_grp)
}
) # end public
) # end class
CohortCtstm$new(trans_model=transmod,
cost_models=list(costs_wlos,costs_transition,costs_starting),
utility_model=utilities)$sim(t=c(0,1,2,10))$summarize()
1
## ## Reconstruct hesim_data from input_mats
## ## This is essentially the reverse of hesim:::expand.hesim_data
## unexpand = function(self) {
## stopifnot(inherits(self, "input_mats"))
## out = list()
## if ("strategy_id" %in% names(self))
## out$strategies = data.frame(strategy_id = unique(self$strategy_id))
## if ("patient_id" %in% names(self)) {
## patient_data = as.data.frame(as.data.table(self))
## patient_data[c("strategy_id","group_id","time_start")] = NULL
## out$patients = unique(patient_data)
## }
## if ("group_id" %in% names(self))
## out$groups = data.frame(group_id = unique(self$group_id))
## if ("time_start" %in% names(self))
## out$time_start = data.frame(time_start = unique(self$time_start))
## structure(out, class="hesim_data")
## }
## make_default_utilities = function(self, ...) {
## stopifnot(inherits(self, "IndivCtstmTrans"),
## hasName(self, "input_data"))
## n_states = nrow(self$trans_mat)
## utilities_tbl = stateval_tbl(data.frame(state_id=1:n_states,
## est=rep(1,n_states)),
## dist="fixed")
## create_StateVals(utilities_tbl, hesim_data=unexpand(self$input_data), ...)
## }
## unexpand(transmod$input_data)
## make_default_utilities(transmod, n=1)$sim(10)
test = C_cohort_ctstm_sim(transmod,
NULL,
new.env(),
start_state=rep(0L,nrow(patients)),
start_age=rep(0,nrow(patients)),
times=c(0,1,2,10),
clock="mixt",
transition_types=c(0L,0L,0L,0L),
dr_costs=0.03,
dr_qalys=0.03)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment