example
draft
Author

Kyle Baron

Published

May 10, 2024

1 Introduction

This post shows how you can fit a multi-state model with mrgsolve. The post is still in draft mode, but I’m pushing it out to the blog in order to share it with someone who was asking about it.

``````library(msm)
library(mrgsolve)
library(dplyr)
library(minqa)
options(pillar.width = Inf)``````

We’ll estimate a multi-state model using the `cav` data from the msm package using both msm and mrgsolve.

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

1.1 Fit the model with msm

We’ll use this transition matrix

``````qq <- rbind(
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``````
``````cav.msm <- msm(
state ~ years,
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 ``````
``plot(cav.msm)``

1.2 Now use mrgsolve

1.2.1 Set up

Modify the data to use with mrgsolve

``````data <- rename(
as_tibble(cav),
ID = PTNUM,
time = years
) ``````

Set up a model with the four states

``````code <- '
\$CMT A1 A2 A3 A4

\$PLUGIN autodec evtools

\$PARAM
k12 = 0.1/2
k21 = 0.1/2
k23 = 0.1/2
k32 = 0.1/2
k14 = 0.1/2
k24 = 0.1/2
k34 = 0.1/2

\$INPUT
firstobs = 0
state = 1
statemax = 1

\$GLOBAL
void resetCMT(databox& self, const double amt, const int cmt) {
evt::ev x = evt::bolus(amt, cmt);
x.evid = 8;
self.push(x);
}

\$MAIN
A1_0 = 1;
A2_0 = 0;
A3_0 = 0;
A4_0 = 0;

\$DES

dxdt_A1 = -A1 * k12 + A2 * k21                       - A1 * k14;
dxdt_A2 =  A1 * k12 - A2 * k21 - A2 * k23 + A3 * k32 - A2 * k24;
dxdt_A3 =  A2 * k23 - A3 * k32                       - A3 * k34;
dxdt_A4 =  A1 * k14 + A2 * k24 + A3 * k34;

\$ERROR

if(state==1 && EVID==0) Y = A1;
if(state==2 && EVID==0) Y = A2;
if(state==3 && EVID==0) Y = A3;
if(state==4 && EVID==0) {
Y = A1 * k14 + A2 * k24 + A3 * k34;
}

if(EVID==0) {
for(int i = 1; i <= 4; ++i) {
resetCMT(self, 0, i);
}
resetCMT(self, 1, state);
}

\$CAPTURE Y
'``````

The key to this model is the `resetCMT()` function, which will reset all compartments to `0` when there is a transition and then initialize the appropriate compartment with a `1`.

1.2.2 Fit

Load the model and set up for estimation

``mod <- mcode("msm", code)``
``Building msm ... done.``
``````theta <- as.numeric(param(mod))
theta <- theta[grep("^k", names(theta))]
tnames <- names(theta)

pred <- function(p, data) {
p <- lapply(p, exp)
names(p) <- tnames
mod <- param(mod, p)
out <- mrgsim_q(mod, data)
-2*sum(log(out\$Y))
}``````

Fit the model

``````fit <- newuoa(
par = log(theta),
fn = pred,
data = data,
control = list(iprint = 0)
)``````

1.2.3 Compare

Compare Objective function value

``fit\$fval``
``[1] 3968.798``
``cav.msm\$minus2loglik``
``[1] 3968.798``

Compare estimates

``````est <- exp(fit\$par)

names(est) <- tnames

est %>% sort()``````
``````       k24        k14        k12        k32        k21        k34        k23
0.04026564 0.04248537 0.12787425 0.13062398 0.22510175 0.30645998 0.34259602 ``````
``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 ``````