Nothing
#' @title A semi-Markov model for cohort simulation
#' @description An R6 class representing a semi-Markov model for
#' cohort simulation.
#' @details A class to represent a continuous time semi-Markov chain, modelled
#' using cohort simulation. As interpreted in \pkg{rdecision}, semi-Markov
#' models
#' may include temporary states and transitions are defined by per-cycle
#' probabilities. Although used widely in health economic modelling, the
#' differences between semi-Markov models and Markov processes introduce
#' some caveats for modellers:
#' \itemize{
#' \item{If there are temporary states, the result will depend on cycle length.}
#' \item{Transitions are specified by their conditional probability, which
#' is a \emph{per-cycle} probability of starting a cycle in one state and
#' ending it in another; if the cycle length changes, the probabilities should
#' change, too.}
#' \item{Probabilities and rates cannot be linked by the Kolmogorov forward
#' equation, where the per-cycle probabilities are given by the matrix
#' exponential of the transition rate matrix, because this equation does not
#' apply if there are temporary states. In creating semi-Markov models, it is
#' the modeller's task to estimate probabilities from published data on
#' event rates.}
#' \item{The cycle time cannot be changed during the simulation.}
#' }
#' @section Graph theory:
#' A Markov model is a directed multidigraph permitting loops (a loop
#' multidigraph), optionally labelled, or a \dfn{quiver}. It is a
#' multidigraph because there are potentially two edges between each pair of
#' nodes \{A,B\} representing the transition probabilities from A to B and
#' \emph{vice versa}. It is a directed graph because the transition
#' probabilities refer to transitions in one direction. Each edge can be
#' optionally labelled. It permits self-loops (edges whose source and target
#' are the same node) to represent patients that remain in the same state
#' between cycles.
#' @section Transition rates and probabilities:
#' \subsection{Why semi-Markov?}{
#' Beck and Pauker (1983) and later Sonnenberg and Beck (1993) proposed the
#' use of Markov processes to model the health economics of medical
#' interventions. Further, they introduced the additional concept of temporary
#' states, to which patients who transition remain for exactly one cycle. This
#' breaks the principle that Markov processes are memoryless
#' and thus the underlying mathematical formalism, first
#' developed by Kolmogorov, is not applicable. For example, ensuring that all
#' patients leave a temporary state requires its transition rate to be infinite.
#' Hence, such models are usually labelled as semi-Markov processes.
#' }
#' \subsection{Rates and probabilities}{
#' Miller and Homan (1994) and Fleurence & Hollenbeak (2007) provide advice
#' on estimating probabilities from rates. Jones (2017) and Welton (2005)
#' describe methods for estimating probabilities in multi-state,
#' multi-transition models, although those methods may not apply to
#' semi-Markov models with temporary states. In particular, note that the
#' "simple" equation, \eqn{p = 1-e^{-rt}} (Briggs 2006) applies only in a
#' two-state, one transition model.
#' }
#' \subsection{Uncertainty in rates}{
#' In semi-Markov models, the conditional probabilities of the transitions
#' from each state are usually modelled by a Dirichlet distribution. In
#' \pkg{rdecision}, create a Dirichlet distribution for each state and
#' optionally create model variables for each conditional probability
#' (\eqn{\rho_{ij}}) linked to an applicable Dirichlet distribution.
#' }
#' @references{
#' Beck JR and Pauker SG. The Markov Process in Medical Prognosis.
#' \emph{Med Decision Making}, 1983;\strong{3}:419–458.
#'
#' Briggs A, Claxton K, Sculpher M. Decision modelling for health economic
#' evaluation. Oxford, UK: Oxford University Press; 2006.
#'
#' Fleurence RL and Hollenbeak CS. Rates and probabilities in economic
#' modelling. \emph{PharmacoEconomics}, 2007;\strong{25}:3--6.
#'
#' Jones E, Epstein D and García-Mochón L. A procedure for deriving
#' formulas to convert transition rates to probabilities for multistate
#' Markov models. \emph{Medical Decision Making} 2017;\strong{37}:779–789.
#'
#' Miller DK and Homan SM. Determining transition probabilities: confusion
#' and suggestions. \emph{Medical Decision Making} 1994;\strong{14}:52-58.
#'
#' Sonnenberg FA, Beck JR. Markov models in medical decision making: a
#' practical guide. \emph{Medical Decision Making}, 1993:\strong{13}:322.
#'
#' Welton NJ and Ades A. Estimation of Markov chain transition probabilities
#' and rates from fully and partially observed data: uncertainty
#' propagation, evidence synthesis, and model calibration.
#' \emph{Medical Decision Making}, 2005;\strong{25}:633-645.
#' }
#' @docType class
#' @author Andrew J. Sims \email{andrew.sims@@newcastle.ac.uk}
#' @export
SemiMarkovModel <- R6::R6Class(
classname = "SemiMarkovModel",
lock_class = TRUE,
inherit = Digraph,
private = list(
# private class variables
smm.tcycle = NULL,
smm.Pt = NULL,
smm.discost = NULL,
smm.disutil = NULL,
smm.pop = NULL,
smm.icycle = NULL,
smm.elapsed = NULL,
# @description Low-level model cycling function
# @details Iterates the model for a single cycle. Updates the state
# populations, cycle number and elapsed time.
# return Matrix of transition counts between states. Row and column names
# are state names.
cycle_pop = function() {
# find the number of transitions between each pair of states by
# multiplying state populations with each column of the transition
# matrix
p_t <- self$transition_probabilities()
n_t <- private$smm.pop * p_t
# update the state populations, cycle count and elapsed time
private$smm.pop <- drop(private$smm.pop %*% p_t)
private$smm.icycle <- private$smm.icycle + 1L
private$smm.elapsed <- private$smm.elapsed + private$smm.tcycle
# return the transition matrix
return(n_t)
},
# @description Low-level model cycling function that calculates values,
# using reference costs supplied as arguments, i.e., without doing graph
# traversal to get values which tend not to change between cycles.
# @details Iterates the model for a single cycle. Updates the state
# populations, cycle number and elapsed time, and calculates values (costs
# and utilities).
# @param hcc.pop See function cycle
# @param hcc.cost See function cycle
# @param hcc.QALY See function cycle
# @param m_transition_costs Square matrix of order number of states, with
# each cell being the cost of transition between a pair of states, as
# returned by self$transition_cost().
# v_state_utilities Vector of utilities, by state.
# v_state_costs Vector of annual occupancy costs, by state.
# @return Matrix with one row per state and columns for population,
# entry cost, total cost and QALY.
cycle_pop_with_values = function(
hcc.pop, hcc.cost, hcc.QALY, m_transition_costs, v_state_utilities,
v_state_costs
) {
# get the cycle time, in years
dty <- as.numeric(private$smm.tcycle, units = "days") / 365.25
# save the state populations and elapsed time before cycling
v_pop_start <- self$get_populations()
elapsed_start <- as.numeric(private$smm.elapsed, units = "days") / 365.25
# cycle the population and get the number of transitions between states
m_nt <- private$cycle_pop()
# calculate transition costs as number of transitions between each
# pair of states, Hadamard multiplied by transition costs, summed by
# the state being entered (entry costs)
m_ct <- m_nt * m_transition_costs
v_ec <- colSums(m_ct)
# get the state populations and elapsed time after cycling
v_pop_end <- self$get_populations()
elapsed_end <- as.numeric(private$smm.elapsed, units = "days") / 365.25
# apply half cycle correction to populations
if (hcc.pop) {
v_pop <- (v_pop_start + v_pop_end) / 2.0
} else {
v_pop <- v_pop_end
}
# apply half cycle correction to costs
if (hcc.cost) {
v_cpop <- (v_pop_start + v_pop_end) / 2.0
elapsed_c <- (elapsed_start + elapsed_end) / 2.0
} else {
v_cpop <- v_pop_end
elapsed_c <- elapsed_end
}
dfc <- 1.0 / (1.0 + private$smm.discost) ^ elapsed_c
v_oc <- v_cpop * v_state_costs * dty * dfc
v_ec <- v_ec * dfc
# apply half cycle correction to QALYs
if (hcc.QALY) {
v_upop <- (v_pop_start + v_pop_end) / 2.0
elapsed_u <- (elapsed_start + elapsed_end) / 2.0
} else {
v_upop <- v_pop_end
elapsed_u <- elapsed_end
}
dfu <- 1.0 / (1.0 + private$smm.disutil) ^ elapsed_u
v_qaly <- v_upop * v_state_utilities * dty * dfu
# return the populations, costs and utilities, by state
N <- sum(v_pop_end)
rm <- matrix(
data = c(
v_pop,
v_oc / N,
v_ec / N,
(v_oc + v_ec) / N,
v_qaly / N
),
nrow = length(v_upop), ncol = 5L, byrow = FALSE,
dimnames = list(
rownames(m_transition_costs),
c("Pop", "OccCost", "EntryCost", "Cost", "QALY")
)
)
return(rm)
}
),
public = list(
#' @description Creates a semi-Markov model for cohort simulation.
#' @details A semi-Markov model must meet the following conditions:
#' \enumerate{
#' \item It must have at least one node and at least one edge.
#' \item All nodes must be of class \code{MarkovState};
#' \item All edges must be of class \code{Transition};
#' \item The nodes and edges must form a digraph whose underlying
#' graph is connected;
#' \item Each state must have at least one outgoing transition (for
#' absorbing states this is a self-loop);
#' \item For each state the sum of outgoing conditional transition
#' probabilities must be one. For convenience, one outgoing transition
#' probability from each state may be set to NA when the
#' probabilities are defined. Typically, probabilities for self
#' loops would be set to NA. Transition probabilities in \eqn{Pt}
#' associated with transitions that are not defined as edges in the
#' graph are zero. Probabilities can be changed between cycles.
#' \item No two edges may share the same source and target nodes (i.e.
#' the digraph may not have multiple edges). This is to ensure that there
#' are no more transitions than cells in the transition matrix.
#' \item The node labels must be unique to the graph.
#' }
#' @param V A list of nodes (\code{MarkovState}s).
#' @param E A list of edges (\code{Transition}s).
#' @param tcycle Cycle length, expressed as an R \code{difftime} object.
#' @param discount.cost Annual discount rate for future costs. Note this
#' is a rate, not a probability (i.e. use 0.035 for 3.5\%).
#' @param discount.utility Annual discount rate for future incremental
#' utility. Note this is a rate, not a probability (i.e. use 0.035
#' for 3.5\%).
#' @return A \code{SemiMarkovModel} object. The population of the first
#' state is set to 1000 and from each state there is an equal
#' conditional probability of each allowed transition.
initialize = function(V, E, tcycle = as.difftime(365.25, units = "days"),
discount.cost = 0.0, discount.utility = 0.0) {
# initialize the base class(es)
super$initialize(V, E)
# check minimum number of nodes and edges
abortifnot(self$order() >= 1L && self$size() >= 1L,
message = "The model must have at least 1 node and 1 edge",
class = "invalid_graph"
)
# check that all nodes inherit from MarkovState
S <- which(
vapply(V, FUN.VALUE = TRUE, FUN = inherits, what = "MarkovState")
)
abortifnot(setequal(seq_along(V), S),
message = "Each node must be a 'MarkovState'.",
class = "invalid_state"
)
# check that all edges inherit from Transition
tedges <- which(
vapply(E, FUN.VALUE = TRUE, inherits, what = "Transition")
)
abortifnot(setequal(seq_along(E), tedges),
message = "Each edge must be a 'Transition'.",
class = "invalid_transition"
)
# check that the underlying graph is connected
abortifnot(self$is_weakly_connected(),
message = "The underlying graph of {V,E} must be connected",
class = "invalid_graph"
)
# check that there are no multiple edges
A <- self$digraph_adjacency_matrix()
abortif(any(A > 1L),
message = "The digraph must not have multiple edges",
class = "multiple_edges"
)
# check that there is at least one outgoing transition from each state
abortif(any(rowSums(A) < 1L),
message = "Every state must have at least one outgoing transition",
class = "missing_transition"
)
# check that the node names are unique
abortifnot(length(unique(self$get_statenames())) == self$order(),
message = "State labels must be unique",
class = "invalid_state_names"
)
# check that the cycle time is an interval
abortifnot(inherits(tcycle, what = "difftime"),
message = "Argument 'tcycle' must be of class 'difftime'.",
class = "invalid_tcycle"
)
private$smm.tcycle <- tcycle
# check and set discounts
abortifnot(is.numeric(discount.cost),
message = "Discount rate must be numeric",
class = "invalid_discount"
)
private$smm.discost <- discount.cost
abortifnot(is.numeric(discount.utility),
message = "Discount rate must be numeric",
class = "invalid_discount"
)
private$smm.disutil <- discount.utility
# set the initial probability matrix (equal chance from each state)
Pt <- matrix(
data = 0.0,
nrow = self$order(), ncol = self$order(),
dimnames = list(
source = self$get_statenames(), target = self$get_statenames()
)
)
for (e in self$edges()) {
is <- self$vertex_index(e$source())
it <- self$vertex_index(e$target())
Pt[is, it] <- 1.0
}
Pt <- Pt / rowSums(Pt)
self$set_probabilities(Pt)
# reset the model to its ground state
self$reset()
# return a new SemiMarkovModel object
return(invisible(self))
},
#' @description Sets transition probabilities.
#' @param Pt Per-cycle transition probability matrix. The row and
#' column labels must be the state names and each row must sum to one.
#' Non-zero probabilities for undefined transitions are not allowed. At
#' most one \code{NA} may appear in each row. If an NA is present in a row,
#' it is replaced by 1 minus the sum of the defined probabilities.
#' @return Updated \code{SemiMarkovModel} object
set_probabilities = function(Pt) {
# check form of Pt
nstates <- self$order()
sn <- self$get_statenames()
abortifnot(
is.matrix(Pt),
nrow(Pt) == nstates && ncol(Pt) == nstates,
setequal(sn, dimnames(Pt)[[1L]]),
setequal(sn, dimnames(Pt)[[2L]]),
message = paste(
"'Pt' must be a matrix",
"of size", nstates, "x", nstates, ".",
"Each row and column must have a unique state name."
),
class = "invalid_Pt"
)
# check NAs and replace them
nNA <- rowSums(is.na(Pt))
abortif(any(nNA > 1L),
message = "No more than one NA per row is allowed",
class = "invalid_Pt"
)
sumP <- rowSums(Pt, na.rm = TRUE)
for (r in seq_len(nstates)) {
if (nNA[[r]] > 0L) {
Pt[[r, which(is.na(Pt[r, ]))]] <- 1.0 - sumP[[r]]
}
}
# check that all cells are a probability
abortif(
any(Pt < 0.0),
any(Pt > 1.0),
message = "All elements of 'Pt' must be probabilities",
class = "invalid_Pt"
)
# reorder rows and columns to match internal vertex order, allowing for
# the possibility of a 1 x 1 model
Pt <- Pt[sn, sn, drop = FALSE]
# check that Pt[i,j]==0 if edge i->j is not in graph
M <- self$digraph_adjacency_matrix(boolean = TRUE)
AA <- (Pt != 0.0)
abortif(sum(M | AA) > sum(M),
message = "All non-zero Pt cells must correspond to transitions",
class = "invalid_Pt"
)
# check that all rows of Pt sum to 1
sumP <- rowSums(Pt)
abortif(any(abs(sumP - 1.0) > sqrt(.Machine$double.eps)),
message = "All rows of Pt must sum to 1",
class = "invalid_Pt"
)
# set the class variable
private$smm.Pt <- Pt
# return updated model
return(invisible(self))
},
#' @description Per-cycle transition probability matrix for the model.
#' @return A square matrix of size equal to the number of states. If all
#' states are labelled, the dimnames take the names of the states.
transition_probabilities = function() {
return(private$smm.Pt)
},
#' @description Return the per-cycle transition costs for the model.
#' @return A square matrix of size equal to the number of states. If all
#' states are labelled, the dimnames take the names of the states.
transition_cost = function() {
# get the state names
state_names <- self$get_statenames()
# construct the matrix
o <- self$order()
Ic <- matrix(
data = 0.0,
nrow = o, ncol = o,
dimnames = list(source = state_names, target = state_names)
)
# populate the cells with costs
for (ie in self$edge_along()) {
e <- self$edge_at(ie)
is <- self$vertex_index(e$source())
it <- self$vertex_index(e$target())
Ic[[is, it]] <- e$cost()
}
return(Ic)
},
#' @description Returns a character list of state names.
#' @return List of the names of each state.
get_statenames = function() {
iv <- self$vertex_along()
statenames <- self$vertex_label(iv)
return(statenames)
},
#' @description Resets the model counters.
#' @details Resets the state populations, next cycle number and elapsed time
#' of the model. By default the model is returned to its ground state (zero
#' people in the all states; next cycle is labelled
#' zero; elapsed time (years) is zero). Any or all of these can be set via
#' this function. \code{icycle} is simply an integer counter label for each
#' cycle, \code{elapsed} sets the elapsed time in years from the index time
#' from which discounting is assumed to apply.
#' @param populations A named vector of populations for
#' the start of the state. The names should be the state names.
#' Due to the R implementation of matrix algebra, \code{populations}
#' must be a numeric type and is not restricted to being an integer. If
#' NULL, the population of all states is set to zero.
#' @param icycle Cycle number at which to start/restart.
#' @param elapsed Elapsed time since the index (reference) time used for
#' discounting as an R \code{difftime} object.
#' @return Updated \code{SemiMarkovModel} object.
reset = function(populations = NULL, icycle = 0L,
elapsed = as.difftime(0.0, units = "days")) {
# get the state names
sn <- self$get_statenames()
# check and set the populations
if (is.null(populations)) {
private$smm.pop <- vector(mode = "numeric", length = self$order())
names(private$smm.pop) <- sn
} else {
abortifnot(
length(populations) == self$order(),
setequal(sn, names(populations)),
is.numeric(populations),
message = paste(
"'populations' must be a numeric vector with one named element",
"per state"
),
class = "invalid_populations"
)
# set the populations
for (s in sn) {
private$smm.pop[s] <- populations[s]
}
}
# check and set the cycle number
abortifnot(
is.integer(icycle),
icycle >= 0L,
message = "'icycle' must be an integer >= 0",
class = "invalid_icycle"
)
private$smm.icycle <- icycle
# check and update the elapsed time
abortifnot(inherits(elapsed, what = "difftime"),
message = "Argument 'elapsed' must be of class 'difftime'.",
class = "invalid_elapsed"
)
private$smm.elapsed <- elapsed
# return updated object
return(invisible(self))
},
#' @description Gets the occupancy of each state.
#' @return A numeric vector of populations, named with state names.
get_populations = function() {
return(private$smm.pop)
},
#' @description Gets the current cycle number.
#' @return Current cycle count, as an integer.
get_cycle = function() {
return(private$smm.icycle)
},
#' @description Gets the cycle duration.
#' @return Current cycle duration, as a difftime object.
get_tcycle = function() {
return(private$smm.tcycle)
},
#' @description Gets the current elapsed time.
#' @details The elapsed time is defined as the difference between the
#' current time in the model and an index time used as the reference
#' time for applying discounting. By default the elapsed time starts at
#' zero. It can be set directly by calling \code{reset}. It is incremented
#' after each call to \code{cycle} by the cycle duration to the time at the
#' end of the cycle (even if half cycle correction is used). Thus, via the
#' \code{reset} and \code{cycle} methods, the time of each cycle relative
#' to the discounting index and its duration can be controlled arbitrarily.
#' @return Elapsed time as an R \code{difftime} object.
get_elapsed = function() {
return(private$smm.elapsed)
},
#' @description Tabulation of states
#' @details Creates a data frame summary of each state in the model.
#' @return A data frame with the following columns:
#' \describe{
#' \item{Name}{State name}
#' \item{Cost}{Annual cost of occupying the state}
#' \item{Utility}{Incremental utility associated with being in the state.}
#' }
tabulate_states = function() {
v <- self$vertexes()
df <- data.frame(
Name = vapply(X = v, FUN.VALUE = "x", FUN = function(x) x$name()),
Cost = vapply(X = v, FUN.VALUE = 1.0, FUN = function(x) x$cost()),
Utility = vapply(X = v, FUN.VALUE = 1.0, FUN = function(x) {
x$utility()
})
)
return(df)
},
#' @description Applies one cycle of the model.
#' @param hcc.pop Determines the state populations returned by this
#' function. If FALSE, the end of cycle populations apply; if TRUE the
#' mid-cycle populations and time apply. The mid-cycle populations are
#' taken as the mean of the start and end populations and the discount
#' time as the mid-point. The value of this parameter does not affect the
#' state populations or elapsed time passed to the next cycle or available
#' via \code{get_populations}; those are always the end cycle values.
#' @param hcc.cost Determines the state occupancy costs returned by this
#' function and the time at which the cost discount is applied to the
#' occupancy costs and the entry costs. If FALSE, the end of cycle
#' populations and time apply; if TRUE the mid-cycle populations and time
#' apply, as per \code{hcc.pop}. The value of this parameter does not affect
#' the state populations or elapsed time passed to the next cycle or
#' available via \code{get_populations}; those are always the end cycle
#' values.
#' @param hcc.QALY Determines the incremental quality adjusted life years
#' returned by this function and the time at which the utility discount is
#' applied. If FALSE, the end of cycle population and reference time are
#' applied to the utilities of each state; if TRUE the mid-cycle populations
#' and time are applied to the state utilities. The value of this parameter
#' does not affect the state populations or elapsed time passed to the next
#' cycle or available via \code{get_populations}; those are always the
#' end cycle values.
#' @return Calculated values, one row per state, as a data frame with the
#' following columns:
#' \describe{
#' \item{\code{State}}{Name of the state.}
#' \item{\code{Cycle}}{The cycle number.}
#' \item{\code{Time}}{Clock time in years of the end of the cycle.}
#' \item{\code{Population}}{Populations of the states; see \code{hcc.pop}.}
#' \item{\code{OccCost}}{Cost of the population occupying the state for
#' the cycle. Discounting and half cycle correction is applied, if those
#' options are set. The costs are normalized by the model population. The
#' cycle costs are derived from the annual occupancy costs of the
#' \code{MarkovState}s.}
#' \item{\code{EntryCost}}{Cost of the transitions \emph{into} the state
#' during the cycle. Discounting is applied, if the option is set.
#' The result is normalized by the model population. The cycle costs
#' are derived from \code{Transition} costs.}
#' \item{\code{Cost}}{Total cost, normalized by model population.}
#' \item{\code{QALY}}{Quality adjusted life years gained by occupancy of
#' the states during the cycle. Half cycle correction and discounting are
#' applied, if these options are set. Normalized by the model population.}
#' }
cycle = function(hcc.pop = TRUE, hcc.cost = TRUE, hcc.QALY = TRUE) {
# check and set half cycle correction parameters
abortifnot(
is.logical(hcc.pop),
is.logical(hcc.cost),
is.logical(hcc.QALY),
message = "Half cycle correction arguments must be logical.",
class = "invalid_hcc"
)
# get the transition costs (matrix), incremental state occupancy
# costs and per-cycle incremental state utilities (vectors)
m_transition_costs <- self$transition_cost()
nstates <- self$order()
v_state_costs <- vector(mode = "numeric", length = nstates)
v_state_utilities <- vector(mode = "numeric", length = nstates)
for (iv in self$vertex_along()) {
v <- self$vertex_at(iv)
v_state_costs[[iv]] <- v$cost()
v_state_utilities[[iv]] <- v$utility()
}
# cycle population and calculate values
m_c <- private$cycle_pop_with_values(
hcc.pop, hcc.cost, hcc.QALY, m_transition_costs, v_state_utilities,
v_state_costs
)
# create return data frame
rc <- data.frame(
State = rownames(m_c),
Cycle = rep(private$smm.icycle, times = nstates),
Time = rep(
as.numeric(private$smm.elapsed, units = "days") / 365.25,
times = nstates
),
Population = m_c[, "Pop"],
OccCost = m_c[, "OccCost"],
EntryCost = m_c[, "EntryCost"],
Cost = m_c[, "Cost"],
QALY = m_c[, "QALY"]
)
return(rc)
},
#' @description Applies multiple cycles of the model.
#' @details The starting populations are redistributed through the
#' transition probabilities and the state occupancy costs are
#' calculated, using function \code{cycle}. The end populations are
#' then fed back into the model for a further cycle and the
#' process is repeated. For each cycle, the state populations and
#' the aggregated occupancy costs are saved in one row of the
#' returned data frame, with the cycle number. If the cycle count
#' for the model is zero when called, the first cycle reported
#' will be cycle zero, i.e. the distribution of patients to starting
#' states.
#' @param ncycles Number of cycles to run; default is 2.
#' @param hcc.pop Determines the state populations returned by this
#' function. If FALSE, the end of cycle populations apply; if TRUE the
#' mid-cycle populations and time apply. The mid-cycle populations are
#' taken as the mean of the start and end populations and the discount
#' time as the mid-point. The value of this parameter does not affect the
#' state populations or elapsed time passed to the next cycle or available
#' via \code{get_populations}; those are always the end cycle values.
#' @param hcc.cost Determines the state occupancy costs returned by this
#' function and the time at which the cost discount is applied to the
#' occupancy costs and the entry costs. If FALSE, the end of cycle
#' populations and time apply; if TRUE the mid-cycle populations and time
#' apply, as per \code{hcc.pop}. The value of this parameter does not affect
#' the state populations or elapsed time passed to the next cycle or
#' available via \code{get_populations}; those are always the end cycle
#' values.
#' @param hcc.QALY Determines the incremental quality adjusted life years
#' returned by this function and the time at which the utility discount is
#' applied. If FALSE, the end of cycle population and reference time are
#' applied to the utilities of each state; if TRUE the mid-cycle populations
#' and time are applied to the state utilities. The value of this parameter
#' does not affect the state populations or elapsed time passed to the next
#' cycle or available via \code{get_populations}; those are always the
#' end cycle values.
#' @return Data frame with cycle results, with the following columns:
#' \describe{
#' \item{\code{Cycle}}{The cycle number.}
#' \item{\code{Years}}{Elapsed time at end of cycle, years}
#' \item{\code{Cost}}{Cost associated with occupancy and transitions between
#' states during the cycle.}
#' \item{\code{QALY}}{Quality adjusted life years associated with occupancy
#' of the states in the cycle.}
#' \item{\code{<name>}}{Population of state \code{<name>} at the end of
#' the cycle.}
#' }
cycles = function(
ncycles = 2L, hcc.pop = TRUE, hcc.cost = TRUE, hcc.QALY = TRUE
) {
# show zero?
nzero <- as.integer(private$smm.icycle == 0L)
# list of state names
statenames <- self$get_statenames()
# construct output matrix
rm <- matrix(
nrow = ncycles + nzero, ncol = 4L + length(statenames),
dimnames = list(NULL, c("Cycle", "Years", "Cost", "QALY", statenames))
)
# add zero
if (nzero > 0L) {
rm[1L, ] <- c(0L, 0.0, 0.0, 0.0, private$smm.pop)
}
# get the transition costs (matrix), incremental state occupancy
# costs and per-cycle incremental state utilities (vectors)
m_transition_costs <- self$transition_cost()
nstates <- self$order()
v_state_costs <- vector(mode = "numeric", length = nstates)
v_state_utilities <- vector(mode = "numeric", length = nstates)
for (iv in self$vertex_along()) {
v <- self$vertex_at(iv)
v_state_costs[[iv]] <- v$cost()
v_state_utilities[[iv]] <- v$utility()
}
# cycle the model
for (i in (1L + nzero) : (ncycles + nzero)) {
# cycle population and calculate values
m_c <- private$cycle_pop_with_values(
hcc.pop, hcc.cost, hcc.QALY,
m_transition_costs, v_state_utilities, v_state_costs
)
# add result to the return matrix (one row per state)
rm[i, ] <- c(
Cycle = private$smm.icycle,
Years = as.numeric(private$smm.elapsed, units = "days") / 365.25,
Cost = sum(m_c[, "Cost"]),
QALY = sum(m_c[, "QALY"]),
m_c[, "Pop"]
)
}
# convert to data frame (note: it is more efficient to return as a matrix)
DF <- as.data.frame(rm)
DF$Cycle <- as.integer(DF$Cycle)
# return summary data frame
return(DF)
},
#' @description Find all the model variables in the Markov model.
#' @details Returns variables of type \code{ModVar} that have been
#' specified as values associated with transition rates and costs and
#' the state occupancy costs and utilities.
#' @return A list of \code{ModVar}s.
modvars = function() {
# create list
mv <- list()
# find the ModVars in the transitions
for (e in self$edges()) {
if (inherits(e, what = "Transition")) {
mv <- c(mv, e$modvars())
}
}
# find the ModVars in the states
for (v in self$vertexes()){
if (inherits(v, what = "MarkovState")) {
mv <- c(mv, v$modvars())
}
}
# return a unique list
return(unique(mv))
},
#' @description Tabulate the model variables in the Markov model.
#' @param expressions A logical that defines whether expression model
#' variables should be included in the tabulation.
#' @return Data frame with one row per model variable, as follows:
#' \describe{
#' \item{\code{Description}}{As given at initialization.}
#' \item{\code{Units}}{Units of the variable.}
#' \item{\code{Distribution}}{Either the uncertainty distribution, if
#' it is a regular model variable, or the expression used to create it,
#' if it is an \code{ExprModVar}.}
#' \item{\code{Mean}}{Mean; calculated from means of operands if
#' an expression.}
#' \item{\code{E}}{Expectation; estimated from random sample if expression,
#' mean otherwise.}
#' \item{\code{SD}}{Standard deviation; estimated from random sample if
#' expression, exact value otherwise.}
#' \item{\code{Q2.5}}{p=0.025 quantile; estimated from random sample if
#' expression, exact value otherwise.}
#' \item{\code{Q97.5}}{p=0.975 quantile; estimated from random sample if
#' expression, exact value otherwise.}
#' \item{\code{Est}}{TRUE if the quantiles and SD have been estimated by
#' random sampling.}
#' }
modvar_table = function(expressions = TRUE) {
# create list of model variables in this decision tree, excluding
# expressions if not wanted
mvlist <- self$modvars()
if (!expressions) {
mvlist <- mvlist[vapply(mvlist, FUN.VALUE = TRUE, FUN = function(v) {
!v$is_expression()
})]
}
# create a data frame of model variables
DF <- data.frame(
Description = vapply(mvlist, FUN.VALUE = "x", FUN = function(x) {
x$description()
}),
Units = vapply(mvlist, FUN.VALUE = "x", FUN = function(x) {
x$units()
}),
Distribution = vapply(mvlist, FUN.VALUE = "x", FUN = function(x) {
x$distribution()
}),
Mean = vapply(mvlist, FUN.VALUE = 1.0, FUN = function(x) {
x$mean()
}),
E = vapply(mvlist, FUN.VALUE = 1.0, FUN = function(x) {
rv <- ifelse(x$is_expression(), x$mu_hat(), x$mean())
return(rv)
}),
SD = vapply(mvlist, FUN.VALUE = 1.0, FUN = function(x) {
rv <- ifelse(x$is_expression(), x$sigma_hat(), x$SD())
return(rv)
}),
Q2.5 = vapply(mvlist, FUN.VALUE = 1.0, FUN = function(x) {
rv <- ifelse(
x$is_expression(),
x$q_hat(probs = 0.025),
x$quantile(probs = 0.025)
)
return(rv)
}),
Q97.5 = vapply(mvlist, FUN.VALUE = 1.0, FUN = function(x) {
rv <- ifelse(
x$is_expression(),
x$q_hat(probs = 0.975),
x$quantile(probs = 0.975)
)
return(rv)
}),
Est = vapply(mvlist, FUN.VALUE = TRUE, FUN = function(exp) {
exp$is_expression()
})
)
# Return the table
return(DF)
}
)
)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.