library(msm)
library(mrgsolve)
library(dplyr)
library(minqa)
options(pillar.width = Inf)
1 Introduction
This post shows how you can fit a simple multi-state model with mrgsolve. We’ll fit a model to the cav
data from the msm package using both msm and mrgsolve and then compare the results.
Some references for the data and analysis include
Sharples, L. D., Jackson, C. H., Parameshwar, J., Wallwork, J. & Large, S. R. Diagnostic accuracy of coronary angiography and risk factors for post-heart-transplant cardiac allograft vasculopathy. Transplantation 76, 679–682 (2003). [Link]
Jackson, C. H. Multi-State Models for Panel Data: The msm Package for R. J. Stat. Softw. 38, 1–28 (2011). [Link]
Multi-state modelling in R: the
msm
package [Link]
This post updated July 2024 to use mrgsolve 1.5.1 which introduces a replace()
function in the evt
namespace (see the evtools plugin). The evt::replace()
syntax is just like evt::bolus()
, but we replace the the amount in the indicated compartment, rather than adding to it. This utilizes EVID = 8
, a long-standing feature in mrgsolve, which can be conveniently called from within your model starting with 1.5.1.
The cav
data is a “series of approximately yearly angiographic examinations of heart transplant recipients. The state at each time is a grade of cardiac allograft vasculopathy (CAV), a deterioration of the arterial walls.” (Link).
There are four states in the CAV data set
- No CAV
- Mild or moderate CAV
- Severe CAV
- Death
The data set includes 2846 observations across 622 subjects.
2 Fit the model with msm
2.1 Transition matrix
The transitions in the data set look like
statetable.msm(state, PTNUM, data=cav)
to
from 1 2 3 4
1 1367 204 44 148
2 46 134 54 48
3 4 13 107 55
We’ll follow the analysis presented by the msm package, using this transition matrix
<- rbind(
qq c(0, 0.25, 0, 0.25),
c(0.166, 0, 0.166, 0.166),
c(0, 0.25, 0, 0.25),
c(0, 0, 0, 0)
) qq
[,1] [,2] [,3] [,4]
[1,] 0.000 0.25 0.000 0.250
[2,] 0.166 0.00 0.166 0.166
[3,] 0.000 0.25 0.000 0.250
[4,] 0.000 0.00 0.000 0.000
This assumes that transitions between State 1 and State 3 must pass through State 2.
2.2 Call to fit the model
Fit the model with deathexact = 4
, indicating that 4 (death) is an absorbing state whose time of entry is known exactly, but with unknown transient state just prior to entering.
<- msm(
cav.msm ~ years,
state subject = PTNUM,
data = cav,
qmatrix = qq,
deathexact = 4
)
cav.msm
Call:
msm(formula = state ~ years, subject = PTNUM, data = cav, qmatrix = qq, deathexact = 4)
Maximum likelihood estimates
Transition intensities
Baseline
State 1 - State 1 -0.17037 (-0.19027,-0.15255)
State 1 - State 2 0.12787 ( 0.11135, 0.14684)
State 1 - State 4 0.04250 ( 0.03412, 0.05294)
State 2 - State 1 0.22512 ( 0.16755, 0.30247)
State 2 - State 2 -0.60794 (-0.70880,-0.52143)
State 2 - State 3 0.34261 ( 0.27317, 0.42970)
State 2 - State 4 0.04021 ( 0.01129, 0.14324)
State 3 - State 2 0.13062 ( 0.07952, 0.21457)
State 3 - State 3 -0.43710 (-0.55292,-0.34554)
State 3 - State 4 0.30648 ( 0.23822, 0.39429)
-2 * log-likelihood: 3968.798
The msm package provides some visualization of the result
plot(cav.msm)
3 Now use mrgsolve
3.1 Data
Modify the data to use with mrgsolve, adding columns for ID
and TIME
<- rename(
data as_tibble(cav),
ID = PTNUM,
time = years
)
The data set is really fairly simple
head(data, n = 4)
# A tibble: 4 × 10
ID age time dage sex pdiag cumrej state firstobs statemax
<int> <dbl> <dbl> <int> <int> <fct> <int> <int> <int> <dbl>
1 100002 52.5 0 21 0 IHD 0 1 1 1
2 100002 53.5 1.00 21 0 IHD 2 1 0 1
3 100002 54.5 2.00 21 0 IHD 2 2 0 2
4 100002 55.6 3.09 21 0 IHD 2 2 0 2
- The
state
column marks the state at the examination, ranging from 1 (no CAV) to 3 (severe CAV) as well as State 4 (death). - Every subject starts in
state
1 - There are also some covariates that we won’t use in this example
You can also verify that there are no records in the data set to reset or dose into any compartment; we’ll handle all of that from inside the model.
3.2 Model
Set up a model with the four states
msm-cav.mod
$PLUGIN autodec evtools
$CMT @annotated: 0 : No CAV
A1 : 0 : Mild or moderate CAV
A2 : 0 : Severe CAV
A3 : 0 : Death
A4
$PARAM @annotated: 0.1 : None to mild or moderate
k12 : 0.1 : Mild or moderate to none
k21 : 0.1 : Mild or moderate to severe
k23 : 0.1 : Severe to mild or moderate
k32 : 0.1 : None to death
k14 : 0.1 : Mild or moderate to death
k24 : 0.1 : Severe to death
k34
$INPUT= 1
state
$PK= 1;
A1_0
$DES= -A1 * k12 + A2 * k21 - A1 * k14;
dxdt_A1 = A1 * k12 - A2 * k21 - A2 * k23 + A3 * k32 - A2 * k24;
dxdt_A2 = A2 * k23 - A3 * k32 - A3 * k34;
dxdt_A3 = A1 * k14 + A2 * k24 + A3 * k34;
dxdt_A4
$ERRORif(EVID != 0) return;
if(state==1) Y = A1;
if(state==2) Y = A2;
if(state==3) Y = A3;
if(state==4) Y = A1 * k14 + A2 * k24 + A3 * k34;
for(int st = 1; st <= 4; ++st) {
::replace(self, 0, st);
evt}
::bolus(self, 1, state);
evt
$CAPTURE Y
The key to this model is the evt::replace()
function, which will reset all compartments to 0
when there is an examination and then initialize the appropriate compartment with a 1
based on the value of state
in the data set. We do this with a simple for
loop in C++. The replace()
functionality is available in the evtools
plugin starting with mrgsolve 1.5.1. More on evtools in the mrgsolve user guide.
3.3 Call to fit the model
Load the model and set up for estimation
<- mread("msm-cav.mod") mod
The initial estimates will be whatever we wrote into the model file
<- as.numeric(param(mod))
theta <- theta[grep("^k", names(theta))]
theta <- names(theta) tnames
This function takes in a set of parameters and returns the -2 log-likelihood returned from the model.
<- function(p, data) {
ofv <- lapply(p, exp)
p names(p) <- tnames
<- param(mod, p)
mod <- mrgsim_q(mod, data)
out -2*sum(log(out$Y))
}
Use minqa::newuoa()
to fit the model
<- newuoa(
fit par = log(theta),
fn = ofv,
data = data,
control = list(iprint = 1)
)
start par. = -2.302585 -2.302585 -2.302585 -2.302585 -2.302585 -2.302585 -2.302585 fn = 4208.405
At return
eval: 267 fn: 3968.7979 par: -2.05671 -1.49120 -1.07120 -2.03543 -3.15860 -3.21224 -1.18267
This takes 4 to 4.5 seconds on my MacBook M1 Pro.
4 Compare msm and mrgsolve results
The final objective function values for the msm and mrgsolve fits are similar
$fval fit
[1] 3968.79787942
$minus2loglik cav.msm
[1] 3968.79789305
Compare estimated transition intensities (there might be a naming issue on the cav.msm
estimates, but I think the values match up).
<- exp(fit$par)
est
names(est) <- tnames
%>% sort() est
k24 k14 k12 k32 k21 k34 k23
0.04026619 0.04248536 0.12787425 0.13062368 0.22510165 0.30645943 0.34259544
exp(cav.msm$estimates) %>% sort()
qbase qbase qbase qbase qbase qbase qbase
0.04021023 0.04250042 0.12787033 0.13062235 0.22511913 0.30647512 0.34261129
```