Skip to contents

Simulate ordinal patient trajectories from a Markov model

Source: R/simulate-markov.R
sim_markov.Rd

Generates daily ordinal states for each patient from a first-order Markov model: the state on day t depends on the state on day t - 1, the patient's covariates, and the day, through a proportional-odds or partial proportional-odds model. Choose a ready-made model of a published trial with preset, or describe your own model with preset = NULL.

The presets are:

  • "actt1": 8 states of the NIAID ordinal scale for a population like that of the ACTT-1 trial of Remdesivir in COVID-19, built from the ACTT-2 model (see actt1_preset).

  • "actt2": 8 states of the NIAID ordinal scale, from a model of the ACTT-2 trial of Baricitinib plus Remdesivir in COVID-19 (see actt2_preset).

  • "violet": 4 states (home, hospital, ventilator or ARDS, dead), from Frank Harrell's model of the VIOLET trial of vitamin D3 in critically ill adults (see violet_preset).

With a preset, n_patients sets the sample size and treatment_effect the treatment effect. Any model argument you supply, such as intercepts or baseline_data, replaces that part of the preset and keeps the rest.

Usage

sim_markov(
  preset,
  n_patients = 1000,
  treatment_prob = 0.5,
  treatment_effect = NULL,
  treatment_effect_decay = 0,
  follow_up_time = 28,
  baseline_data = NULL,
  intercepts = NULL,
  lp_function = NULL,
  extra_params = NULL,
  absorbing_states = NULL,
  covariate_names = NULL,
  seed = NULL
)

Arguments

preset

The model to simulate from: "actt1", "actt2", or "violet", or NULL to supply your own model through baseline_data, intercepts, and lp_function. There is no default.

n_patients

Number of patients to sample from the preset's population. Must be a positive integer (default: 1000). Cannot be combined with baseline_data.

treatment_prob

Probability that a sampled patient is assigned to treatment (tx = 1); the other patients are controls (tx = 0). Must be between 0 and 1 (default: 0.5). Cannot be combined with baseline_data.

treatment_effect

Log odds ratio of being in a worse state on day 1 for treated (tx = 1) versus control (tx = 0) patients. Values below 0 mean that treatment helps. NULL (the default) uses the preset's value: -0.0730 for "actt1" and "actt2", and log(0.75) for "violet". With preset = NULL it means 0, no effect.

treatment_effect_decay

Nonnegative exponential decay rate of the treatment effect per day. The effect on day t is treatment_effect * exp(-treatment_effect_decay * (t - 1)), so -log(0.01) / 14 shrinks it to 1% of its day-1 size by day 15. The default 0 gives a constant effect.

follow_up_time

Number of days to simulate (default: 28).

baseline_data

A data frame with one row per patient, replacing the patients that a preset would sample. It must contain id, yprev (the state on the day before the first simulated day), and every covariate that the model uses: tx for "actt1" and "actt2", and age, sofa, and tx for "violet". Required when preset = NULL.

intercepts

Numeric vector of thresholds, one fewer than the number of states, in increasing order. With K states numbered from 1 (best) to K (worst), intercepts[j] plus the linear predictor is the log odds of being in state K + 1 - j or worse: intercepts[1] belongs to the worst state alone and intercepts[K - 1] to states 2 and worse. Required when preset = NULL; otherwise NULL uses the preset's thresholds.

lp_function

A function that calculates the linear predictor, which is added to the thresholds. It is called with these named arguments:

  • yprev: previous state, a factor with levels 1 to K

  • t: current day

  • one argument for each column named in covariate_names, and tx when baseline_data has a tx column

  • treatment_effect: the treatment effect on day t, already reduced by treatment_effect_decay

  • extra_params: the value of extra_params

Each argument holds one patient's value. The function must return either one number, added to every threshold (a proportional-odds model), or a numeric vector with length(intercepts) elements, added to the corresponding thresholds (a partial proportional-odds model). Larger values make worse states more likely. Required when preset = NULL.

extra_params

Optional object, usually a named numeric vector of coefficients, passed unchanged to lp_function. NULL uses the preset's coefficients.

absorbing_states

Integer vector of states that patients cannot leave once they enter them. NULL (the default) means the worst state, length(intercepts) + 1, which is death in every preset. Use integer() when no state is absorbing.

covariate_names

Character vector of baseline_data columns passed to lp_function as arguments with the same names. A tx column is always passed. NULL (the default) uses the preset's covariates, or, with preset = NULL, every column except id and yprev.

seed

Integer. Random seed for reproducibility (default: NULL, no seed set). The seed is set once, before patients and treatment are sampled, and the trajectories continue from the same random-number stream.

Value

A data frame with one row per patient and day and columns:

  • id: patient identifier

  • time: day, from 1 to follow_up_time

  • y: state on this day

  • yprev: state on the previous day

  • all other columns of baseline_data, such as tx

Details

On each day the function:

  1. calculates the linear predictor for each patient not in an absorbing state,

  2. adds it to the thresholds to get the probability of each state or worse,

  3. converts these to the probability of each state,

  4. samples the next state for each patient, and

  5. keeps patients in absorbing states where they are.

States must be consecutive integers from 1 (best) to K (worst). If lp_function returns threshold-specific linear predictors, the probability of being in a state or worse must not increase from one state to the next worse state; otherwise the function stops because the implied state probabilities would be negative.

Each preset's thresholds follow the order described under intercepts. To change one of them, start from the preset's help page, which lists them.

See also

actt1_preset, actt2_preset, and violet_preset for the presets' models.

Examples

# 50 patients from the ACTT-2 model
actt2 <- sim_markov("actt2", n_patients = 50, seed = 123)
head(actt2)
#>   id time y yprev tx
#> 1  1    1 5     5  0
#> 2  1    2 5     5  0
#> 3  1    3 5     5  0
#> 4  1    4 5     5  0
#> 5  1    5 4     5  0
#> 6  1    6 4     4  0

# A treatment effect that fades to 1% of its day-1 size by day 15
decaying <- sim_markov(
  "actt2",
  n_patients = 50,
  treatment_effect = log(0.8),
  treatment_effect_decay = -log(0.01) / 14,
  seed = 123
)

# Your own patients with the ACTT-1 model
patients <- data.frame(id = 1:4, yprev = c(5, 5, 7, 7), tx = c(0, 1, 0, 1))
sim_markov("actt1", baseline_data = patients, follow_up_time = 7, seed = 1)
#>    id time y yprev tx
#> 1   1    1 5     5  0
#> 2   1    2 5     5  0
#> 3   1    3 5     5  0
#> 4   1    4 5     5  0
#> 5   1    5 5     5  0
#> 6   1    6 4     5  0
#> 7   1    7 4     4  0
#> 8   2    1 5     5  1
#> 9   2    2 4     5  1
#> 10  2    3 5     4  1
#> 11  2    4 5     5  1
#> 12  2    5 2     5  1
#> 13  2    6 1     2  1
#> 14  2    7 1     1  1
#> 15  3    1 7     7  0
#> 16  3    2 7     7  0
#> 17  3    3 7     7  0
#> 18  3    4 7     7  0
#> 19  3    5 7     7  0
#> 20  3    6 7     7  0
#> 21  3    7 8     7  0
#> 22  4    1 7     7  1
#> 23  4    2 7     7  1
#> 24  4    3 7     7  1
#> 25  4    4 7     7  1
#> 26  4    5 7     7  1
#> 27  4    6 7     7  1
#> 28  4    7 7     7  1

# Your own model with three states: 1 = home, 2 = in hospital, 3 = dead
baseline <- data.frame(id = 1:20, yprev = 2, tx = rep(0:1, 10))
hospital_lp <- function(yprev, t, tx, treatment_effect, extra_params) {
  # Hospitalized patients are more likely to stay in hospital or die
  5 * (yprev == "2") - 0.1 * t + treatment_effect * tx
}
trajectories <- sim_markov(
  NULL,
  baseline_data = baseline,
  follow_up_time = 10,
  intercepts = c(-9, -3),
  lp_function = hospital_lp,
  treatment_effect = log(0.8),
  seed = 1
)
head(trajectories)
#>   id time y yprev tx
#> 1  1    1 2     2  0
#> 2  1    2 1     2  0
#> 3  1    3 1     1  0
#> 4  1    4 1     1  0
#> 5  1    5 1     1  0
#> 6  1    6 1     1  0

# Partial proportional odds: treatment affects only the odds of death
death_only_lp <- function(yprev, t, tx, treatment_effect, extra_params) {
  hospital_lp(yprev, t, tx = 0, treatment_effect, extra_params) +
    c(treatment_effect * tx, 0)
}
trajectories_ppo <- sim_markov(
  NULL,
  baseline_data = baseline,
  follow_up_time = 10,
  intercepts = c(-9, -3),
  lp_function = death_only_lp,
  treatment_effect = log(0.8),
  seed = 1
)