SemiMarkovModel {rdecision}

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 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:

• If there are temporary states, the result will depend on cycle length.

• Transitions are specified by their conditional probability, which is a 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.

• 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.

• The cycle time cannot be changed during the simulation.

Graph theory

A Markov model is a directed multidigraph permitting loops (a loop multidigraph), optionally labelled, or a 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 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.

Transition rates and probabilities

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.

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, p = 1-e^{-rt} (Briggs 2006) applies only in a two-state, one transition model.

Uncertainty in rates

In semi-Markov models, the conditional probabilities of the transitions from each state are usually modelled by a Dirichlet distribution. In rdecision, create a Dirichlet distribution for each state and optionally create model variables for each conditional probability (\rho_{ij}) linked to an applicable Dirichlet distribution.

Super classes

rdecision::Graph -> rdecision::Digraph -> SemiMarkovModel

Methods

Public methods

• SemiMarkovModel$new()

• SemiMarkovModel$set_probabilities()

• SemiMarkovModel$transition_probabilities()

• SemiMarkovModel$transition_cost()

• SemiMarkovModel$get_statenames()

• SemiMarkovModel$reset()

• SemiMarkovModel$get_populations()

• SemiMarkovModel$get_cycle()

• SemiMarkovModel$get_tcycle()

• SemiMarkovModel$get_elapsed()

• SemiMarkovModel$tabulate_states()

• SemiMarkovModel$cycle()

• SemiMarkovModel$cycles()

• SemiMarkovModel$modvars()

• SemiMarkovModel$modvar_table()

• SemiMarkovModel$clone()

Inherited methods

• rdecision::Graph$degree()
• rdecision::Graph$edge_along()
• rdecision::Graph$edge_at()
• rdecision::Graph$edge_index()
• rdecision::Graph$edge_label()
• rdecision::Graph$edges()
• rdecision::Graph$graph_adjacency_matrix()
• rdecision::Graph$has_edge()
• rdecision::Graph$has_vertex()
• rdecision::Graph$is_simple()
• rdecision::Graph$neighbours()
• rdecision::Graph$order()
• rdecision::Graph$size()
• rdecision::Graph$vertex_along()
• rdecision::Graph$vertex_at()
• rdecision::Graph$vertex_index()
• rdecision::Graph$vertex_label()
• rdecision::Graph$vertexes()
• rdecision::Digraph$arrow_source()
• rdecision::Digraph$arrow_target()
• rdecision::Digraph$as_DOT()
• rdecision::Digraph$digraph_adjacency_matrix()
• rdecision::Digraph$digraph_incidence_matrix()
• rdecision::Digraph$direct_predecessors()
• rdecision::Digraph$direct_successors()
• rdecision::Digraph$is_acyclic()
• rdecision::Digraph$is_arborescence()
• rdecision::Digraph$is_connected()
• rdecision::Digraph$is_polytree()
• rdecision::Digraph$is_tree()
• rdecision::Digraph$is_weakly_connected()
• rdecision::Digraph$paths()
• rdecision::Digraph$topological_sort()
• rdecision::Digraph$walk()

Method new()

Creates a semi-Markov model for cohort simulation.

Usage

SemiMarkovModel$new(
  V,
  E,
  tcycle = as.difftime(365.25, units = "days"),
  discount.cost = 0,
  discount.utility = 0
)

Arguments

V

A list of nodes (MarkovStates).

E

A list of edges (Transitions).

tcycle

Cycle length, expressed as an R difftime object.

discount.cost

Annual discount rate for future costs. Note this is a rate, not a probability (i.e. use 0.035 for 3.5%).

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%).

Details

A semi-Markov model must meet the following conditions:

1. It must have at least one node and at least one edge.

2. All nodes must be of class MarkovState;

3. All edges must be of class Transition;

4. The nodes and edges must form a digraph whose underlying graph is connected;

5. Each state must have at least one outgoing transition (for absorbing states this is a self-loop);

6. 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 Pt associated with transitions that are not defined as edges in the graph are zero. Probabilities can be changed between cycles.

7. 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.

8. The node labels must be unique to the graph.

Returns

A 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.

Method set_probabilities()

Sets transition probabilities.

Usage

SemiMarkovModel$set_probabilities(Pt)

Arguments

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 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.

Returns

Updated SemiMarkovModel object

Method transition_probabilities()

Per-cycle transition probability matrix for the model.

Usage

SemiMarkovModel$transition_probabilities()

Returns

A square matrix of size equal to the number of states. If all states are labelled, the dimnames take the names of the states.

Method transition_cost()

Return the per-cycle transition costs for the model.

Usage

SemiMarkovModel$transition_cost()

Returns

A square matrix of size equal to the number of states. If all states are labelled, the dimnames take the names of the states.

Method get_statenames()

Returns a character list of state names.

Usage

SemiMarkovModel$get_statenames()

Returns

List of the names of each state.

Method reset()

Resets the model counters.

Usage

SemiMarkovModel$reset(
  populations = NULL,
  icycle = 0L,
  elapsed = as.difftime(0, units = "days")
)

Arguments

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, 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.

icycle

Cycle number at which to start/restart.

elapsed

Elapsed time since the index (reference) time used for discounting as an R difftime object.

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. icycle is simply an integer counter label for each cycle, elapsed sets the elapsed time in years from the index time from which discounting is assumed to apply.

Returns

Updated SemiMarkovModel object.

Method get_populations()

Gets the occupancy of each state.

Usage

SemiMarkovModel$get_populations()

Returns

A numeric vector of populations, named with state names.

Method get_cycle()

Gets the current cycle number.

Usage

SemiMarkovModel$get_cycle()

Returns

Current cycle count, as an integer.

Method get_tcycle()

Gets the cycle duration.

Usage

SemiMarkovModel$get_tcycle()

Returns

Current cycle duration, as a difftime object.

Method get_elapsed()

Gets the current elapsed time.

Usage

SemiMarkovModel$get_elapsed()

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 reset. It is incremented after each call to cycle by the cycle duration to the time at the end of the cycle (even if half cycle correction is used). Thus, via the reset and cycle methods, the time of each cycle relative to the discounting index and its duration can be controlled arbitrarily.

Returns

Elapsed time as an R difftime object.

Method tabulate_states()

Tabulation of states

Usage

SemiMarkovModel$tabulate_states()

Details

Creates a data frame summary of each state in the model.

Returns

A data frame with the following columns:

Name

State name

Cost

Annual cost of occupying the state

Utility

Incremental utility associated with being in the state.

Method cycle()

Applies one cycle of the model.

Usage

SemiMarkovModel$cycle(hcc.pop = TRUE, hcc.cost = TRUE, hcc.QALY = TRUE)

Arguments

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 get_populations; those are always the end cycle values.

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 hcc.pop. The value of this parameter does not affect the state populations or elapsed time passed to the next cycle or available via get_populations; those are always the end cycle values.

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 get_populations; those are always the end cycle values.

Returns

Calculated values, one row per state, as a data frame with the following columns:

State

Name of the state.

Cycle

The cycle number.

Time

Clock time in years of the end of the cycle.

Population

Populations of the states; see hcc.pop.

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 MarkovStates.

EntryCost

Cost of the transitions 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 Transition costs.

Cost

Total cost, normalized by model population.

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.

Method cycles()

Applies multiple cycles of the model.

Usage

SemiMarkovModel$cycles(
  ncycles = 2L,
  hcc.pop = TRUE,
  hcc.cost = TRUE,
  hcc.QALY = TRUE
)

Arguments

ncycles

Number of cycles to run; default is 2.

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 get_populations; those are always the end cycle values.

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 hcc.pop. The value of this parameter does not affect the state populations or elapsed time passed to the next cycle or available via get_populations; those are always the end cycle values.

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 get_populations; those are always the end cycle values.

Details

The starting populations are redistributed through the transition probabilities and the state occupancy costs are calculated, using function 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.

Returns

Data frame with cycle results, with the following columns:

Cycle

The cycle number.

Years

Elapsed time at end of cycle, years

Cost

Cost associated with occupancy and transitions between states during the cycle.

QALY

Quality adjusted life years associated with occupancy of the states in the cycle.

<name>

Population of state <name> at the end of the cycle.

Method modvars()

Find all the model variables in the Markov model.

Usage

SemiMarkovModel$modvars()

Details

Returns variables of type ModVar that have been specified as values associated with transition rates and costs and the state occupancy costs and utilities.

Returns

A list of ModVars.

Method modvar_table()

Tabulate the model variables in the Markov model.

Usage

SemiMarkovModel$modvar_table(expressions = TRUE)

Arguments

expressions

A logical that defines whether expression model variables should be included in the tabulation.

Returns

Data frame with one row per model variable, as follows:

Description

As given at initialization.

Units

Units of the variable.

Distribution

Either the uncertainty distribution, if it is a regular model variable, or the expression used to create it, if it is an ExprModVar.

Mean

Mean; calculated from means of operands if an expression.

E

Expectation; estimated from random sample if expression, mean otherwise.

SD

Standard deviation; estimated from random sample if expression, exact value otherwise.

Q2.5

p=0.025 quantile; estimated from random sample if expression, exact value otherwise.

Q97.5

p=0.975 quantile; estimated from random sample if expression, exact value otherwise.

Est

TRUE if the quantiles and SD have been estimated by random sampling.

Method clone()

The objects of this class are cloneable with this method.

Usage

SemiMarkovModel$clone(deep = FALSE)

Arguments

deep

Whether to make a deep clone.

Author(s)

Andrew J. Sims [e-mail]

References

Beck JR and Pauker SG. The Markov Process in Medical Prognosis. Med Decision Making, 1983;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. PharmacoEconomics, 2007;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. Medical Decision Making 2017;37:779–789.

Miller DK and Homan SM. Determining transition probabilities: confusion and suggestions. Medical Decision Making 1994;14:52-58.

Sonnenberg FA, Beck JR. Markov models in medical decision making: a practical guide. Medical Decision Making, 1993: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. Medical Decision Making, 2005;25:633-645.