The possible transitions to and from each state are summarized with a
matrix in an individual patient simulation.
tmat <- rbind(c(NA, 1, NA, NA, 2),
c(NA, NA, 3, NA, 4),
c(NA, NA, NA, 5, 6),
c(NA, NA, 7, NA, 8),
c(NA, NA, NA, NA, NA))
colnames(tmat) <- rownames(tmat) <- names(labs$state_id)
tmat %>%
kable() %>%
kable_styling()
Estimates from
literature
Transition 1
The model for the first transition determines the time from the
original surgery until it is deemed “successful”. In the Markov model,
we used model cycles that were 1-year long meaning that this process
took 1-year. Since there are no model cycles in a CTSTM, we can be more
flexible. Let’s assume that time to recovery (TTR) from surgery takes on
average 6 months and follows an exponential distribution.
ttrrPTHR <- 2 # Time to recovery rate implies mean time of 1/2 years
Transition 2
The second transition depends on the operative mortality rate
following primary THR (omrPTHR). As in the cohort model we will
assume the probability of death (within 1-year) is drawn from a beta
distributions.
# 2 out of 100 patients receiving primary THR died
omrPTHR_shape1 <- 2
omrPTHR_shape2 <- 98
Below we will convert this probability to a rate so that is can be
simulated using an exponential distribution. We write a simple function
to do this.
prob_to_rate <- function(p, t = 1){
(-log(1 - p))/t
}
Transition 3
Transition 3 depends on the revision rate (rr) for a
prosthesis, which will again be modeled using a proportional hazards
Weibull model. The coefficients and variance-covariance matrix are
displayed below.
## lngamma cons age male np1
## 0.3740968 -5.4909350 -0.0367022 0.7685360 -1.3444740
## lngamma cons age male np1
## lngamma 0.002251512 -0.00569100 2.800000e-08 0.00000510 0.0002590
## cons -0.005691000 0.04321908 -7.830000e-04 -0.00724700 -0.0006420
## age 0.000000028 -0.00078300 2.715661e-05 0.00003300 -0.0001110
## male 0.000005100 -0.00724700 3.300000e-05 0.01189539 0.0001840
## np1 0.000259000 -0.00064200 -1.110000e-04 0.00018400 0.1463686
Transition 4
The fourth transition is modeled using the death rate following
recovery from the THR (i.e., while in the successful primary state). We
will again use age and sex specific mortality rates.
Since we are running the simulation for a 60-year old female, the
mortality rate at time \(0\) is \(0.0067\). We are using a clock-forward
model, so the hazard depends on time since the start of the model. We
can therefore model mortality over time with a piecewise exponential
distribution with rates that vary over time. For a 60-year female, the
rates will change at 5, 15, and 25 years.
mr <- c(.0067, .0193, .0535, .1548)
mr_times <- c(0, 5, 15, 25)
Transition 5
Transition 5 is similar to transition 1 in that it is a function of
time to recovery from a THR. While we could again model this process
with an exponential distribution, we will instead assume this time is
“fixed” for illustration purposes. Specifically, let’s assume it takes
one year for every patient.
ttrRTHR <- 1 # There is no rate, the time is fixed
Transition 6
Following the prior example, the sixth transition depends on the
operative mortality rate following revision THR (omrRTHR) and
the overall mortality rate (mr). The omrRTHR is
assumed to follow the same distribution as the omrPTHR.
These probabilities will be converted to rates below and then added
to overall mortality rate (mr) at times 0, 5, 15, and 25 years.
Transition 6 will consequently be modeled using a piecewise exponential
distribution.
Transition 7
The re-revision rate (rrr) is used to model the seventh
transition. Like the omrPTHR, the probability is modeled with a
beta distribution and will be converted to a rate below.
# 4 out of 100 patients with a successful revision needed another procedure r
omrRTHR_shape1 <- 4
omrRTHR_shape2 <- 96
Transition 8
The final transition is time to death following a successful
revision. We model it with the same model used for transition 4; that
is, we assume mortality rates following a successful revision THR are
the same as mortality rates following a successful primary THR.
Multi-state
model
The transition-specific estimates can be combined to create a
multi-state model, which is comprised of survival models for each
transition. Since we will perform a probabilistic sensitivity analysis
(PSA), we will sample the parameters from the distributions described
above. 500 iterations will be used.
Survival models are regression models, so the parameters are
regression coefficients. The coefficients of each model must be a matrix
or data frame (rows for parameter samples and columns for covariates).
Since we often use intercept only models, let’s write a simple function
to convert a vector to a data.table with a single column
for the intercept.
vec_to_dt <- function(v, n = NULL){
if (length(v) == 1) v <- rep(v, n_samples)
dt <- data.table(v)
colnames(dt) <- "cons"
return(dt)
}
Our new helper function can then be used within
define_rng() to sample distributions of the
coefficients.
transmod_coef_def <- define_rng({
omrPTHR <- prob_to_rate(beta_rng(shape1 = omrPTHR_shape1,
shape2 = omrPTHR_shape2))
mr <- fixed(mr)
mr_omrPTHR <- omrPTHR + mr
colnames(mr) <- colnames(mr_omrPTHR) <- paste0("time_", mr_times)
rr <- multi_normal_rng(mu = rr_coef, Sigma = rr_vcov)
rrr <- prob_to_rate(beta_rng(shape1 = 4, shape2 = 96))
list(
log_omrPTHR = vec_to_dt(log(omrPTHR)),
log_mr = log(mr),
log_ttrrPTHR = vec_to_dt(log(ttrrPTHR)),
log_mr_omrPTHR = log(mr_omrPTHR),
rr_shape = vec_to_dt(rr$lngamma),
rr_scale = rr[, -1,],
log_rrr = vec_to_dt(log(rrr))
)
}, n = n_samples)
transmod_coef <- eval_rng(transmod_coef_def)
To get a better understanding of the output, lets summarize the
coefficients.
## param mean sd 2.5% 97.5%
## <char> <num> <num> <num> <num>
## 1: log_omrPTHR_cons -4.10477726 0.779672371 -5.76310845 -2.88705363
## 2: log_mr_time_0 -5.00564775 0.000000000 -5.00564775 -5.00564775
## 3: log_mr_time_5 -3.94765018 0.000000000 -3.94765018 -3.94765018
## 4: log_mr_time_15 -2.92807363 0.000000000 -2.92807363 -2.92807363
## 5: log_mr_time_25 -1.86562132 0.000000000 -1.86562132 -1.86562132
## 6: log_ttrrPTHR_cons 0.69314718 0.000000000 0.69314718 0.69314718
## 7: log_mr_omrPTHR_time_0 -3.70117065 0.510505368 -4.62115646 -2.77354455
## 8: log_mr_omrPTHR_time_5 -3.26042174 0.337918149 -3.79684654 -2.58972850
## 9: log_mr_omrPTHR_time_15 -2.60918341 0.185042100 -2.87101444 -2.21420249
## 10: log_mr_omrPTHR_time_25 -1.73939478 0.081331045 -1.84553093 -1.55807582
## 11: rr_shape_cons 0.37255535 0.046646488 0.28768668 0.46684078
## 12: rr_scale_cons -5.48054126 0.206073052 -5.86057729 -5.10467604
## 13: rr_scale_age -0.03689076 0.005386752 -0.04699827 -0.02647981
## 14: rr_scale_male 0.76739930 0.107506239 0.56978939 0.97699170
## 15: rr_scale_np1 -1.33841288 0.383069512 -2.12174005 -0.62716192
## 16: log_rrr_cons -3.30694248 0.540256920 -4.48356733 -2.42015321
It is also helpful to take a look at a few of the sampled coefficient
data tables.
head(transmod_coef$log_omrPTHR)
## cons
## <num>
## 1: -3.225382
## 2: -4.166867
## 3: -4.217561
## 4: -4.518022
## 5: -3.915191
## 6: -4.319199
head(transmod_coef$rr_scale)
## cons age male np1
## <num> <num> <num> <num>
## 1: -5.370253 -0.04343824 0.5988311 -2.0592318
## 2: -5.710184 -0.02507501 0.6306379 -0.3491912
## 3: -5.789766 -0.03138916 0.5626969 -1.6720223
## 4: -5.447601 -0.03706895 0.6568232 -1.5348107
## 5: -4.978150 -0.05180358 0.8589875 -1.7743936
## 6: -5.627855 -0.02914147 0.7767687 -1.2065555
In addition to the coefficients, a complete parameterization of a
transition model requires specification of the survival distributions
and potentially auxiliary information (e.g., the times at which rates
change in a piecewise exponential model). The parameters are stored in a
params_surv_list() object, which is a list of
params_surv() objects.
To create the params_surv() objects, it’s helpful to
first write another short helper function that can convert a single data
table storing the rates for each time in a piecewise exponential to a
list of regression coefficients data tables.
as_dt_list <- function(x) {
lapply(as.list(x), vec_to_dt)
}
The params_surv_list() object is then created as
follows.
transmod_params <- params_surv_list(
# 1. Primary THR:Successful primary (1:2)
params_surv(coefs = list(rate = transmod_coef$log_ttrrPTHR),
dist = "exp"),
# 2. Primary THR:Death (1:5)
params_surv(coefs = list(rate = transmod_coef$log_omrPTHR),
dist = "exp"),
# 3. Successful primary:Revision THR (2:3)
params_surv(coefs = list(shape = transmod_coef$rr_shape,
scale = transmod_coef$rr_scale),
dist = "weibullPH"),
# 4. Successful primary:Death (2:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp"),
# 5. Revision THR:Successful revision (3:4)
params_surv(coefs = list(est = vec_to_dt(ttrRTHR)),
dist = "fixed"),
# 6. Revision THR:Death (3:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr_omrPTHR),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp"),
# 7. Successful revision:Revision THR (4:3)
params_surv(coefs = list(rate = transmod_coef$log_rrr),
dist = "exp"),
# 8. Successful revision:Death (4:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp")
)
