Simulation ofmultistate models withmultiple timescales:
simLexis in the Epi package
SDCC

th July, 2024
Thursday 18

http://BendixCarstensen.com/Epi/simLexis.pdf
Version 2.6

Compiled Thursday 18th July, 2024, 15:27

from: C:\Bendix\stat\r\lib.src\epi\pkg\vignettes/simLexis.tex

Bendix Carstensen

Steno Diabetes Center Copenhagen, Gentofte, Denmark& Department of Biostatistics, University of Copenhagenb@bxc.dk
http://BendixCarstensen.com

1 Using simLexis
1.1
1.2
1.3
1.4
1.5
1.6
1.7

1.8

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
simLexis in practice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Input for the simulation . . . . . . . . . . . . . . . . . . . . . . . . .
Setting up a Lexis object . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Analysis of rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The mortality rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5.1 Proportionality of mortality rates . . . . . . . . . . . . . . . . . . . .
1.5.2 How the mortality rates look . . . . . . . . . . . . . . . . . . . . . . .
Input to the simLexis function . . . . . . . . . . . . . . . . . . . . . . . . .
1.6.1 The transition object . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.6.2 The initial cohort . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Simulation of the follow-up . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.7.1 Using other models for simulation . . . . . . . . . . . . . . . . . . . .
Proportional hazards Poisson model . . . . . . . . . . . . . . . . . . .
Proportional hazards Cox model . . . . . . . . . . . . . . . . . . . . .
Reporting the simulation results . . . . . . . . . . . . . . . . . . . . . . . . .
1.8.1 Comparing predictions from dierent models . . . . . . . . . . . . . .

2

2
3
3
4
5
8
9
10
12
13
13
14
15
15
15
16
20

2 Simulation of transitions in multistate models

23

References

25

2.1

Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ii

23

1
Start time: 2024-07-18, 15:26:45

Chapter 1
Using simLexis
1.1

Introduction

This vignette explains the machinery behind simulation of life histories through multistatemodels implemented in simLexis. In simLexis transition rates are allowed to depend on
multiple time scales, including timescales dened as time since entry to a particular state(duration). This therefore also covers the case where time at entry into a state is an
explanatory variable for the rates, since time at entry is merely (current) time minus
duration. Thus, the set-up here goes beyond Markov- and semi-Markov-models, and brings
simulation based estimation of state-occupancy probabilities into the realm of realistic
multistate models.
The basic idea is to simulate a new Lexis object [3, 1] as dened in the Epi package forR, based on 1) a multistate model dened by its states and the transition rates between
them and 2) an initial population of individuals.
Thus the output will be a Lexis object describing the transitions of a predened set of
persons through a multistate model. Therefore, if persons are dened to be identical atstart, then calculation of the probability of being in a particular state at a given time boils
down to a simple enumeration of the fraction of the persons in the particular state at thegiven time. Bar of course the (binomial) simulation error, but this can be brought down by
simulation a suciently large number of persons.
An observed Lexis object with follow-up of persons through a number of states will
normally be the basis for estimation of transition rates between states, and thus willcontain all information about covariates determining the occurrence rates, in particular thetimescales [2]. Hence, the natural input to simulation from an estimated multistate model
will typically be an object of the same structure as the originally observed. Sincetransitions and times are what is simulated, any values of lex.Xst and lex.dur in the
input object will of course be ignored.
This rst chapter of this vignette shows by an example how to use the function
simLexis and display the results. The second chapter discusses in more detail how thesimulation machinery is implemented and is not needed for the practical use of simLexis.
2

Using simLexis

1.2

1.2 simLexis in practice3

simLexis in practiceThis section is largely a commented walk-trough of the example from the help-page ofsimLexis, with a larger number of simulated persons in order to minimize the pure
simulation variation.
When we want to simulate transition times through a multistate model where transition
rates may depend on time since entry to the current or a previous state, it is essential that
we have a machinery to keep track of the transition time on all time scales, as well as amechanism that can initiate a new time scale to 0 when a transition occurs to a state
where we shall use time since entry as determinant of exit rates from that state. This isprovided by simLexis.

1.2.1 Input for the simulation
Input for simulation of a single trajectory through a multistate model requires a
representation of the current status of a person; the starting conditions. The object that
we supply to the simulation function must contain information about all covariates and alltimescales upon which transitions depend, and in particular which one(s) of the timescales
that are dened as time since entry into a particular state. Hence, starting conditions
should be represented as a Lexis object (where lex.dur and lex.Xst are ignored, since
there is no follow-up yet), where the time scale information is in the attributes
time.scales and time.since respectively.
Note that time.scales attribute is a vector of names of variables in the Lexis object,so all of these variables should be present even if they are not used in the models for thetransitions, and they should be set to 0; if they are not in the initial dataset, simLexis willcrash, if they are NA, the simLexis will produce an object with 0 rows.
Thus there are two main arguments to a function to simulate from a multistate model:
1. A Lexis object representing the initial states and covariates of the population to be
simulated. This has to have the same structure as the original Lexis object
representing the multistate model from which transition rates in the model were
estimated. As noted above, the values for lex.Xst and lex.dur are not required(since these are the quantities that will be simulated).
2. A transition object, representing the transition intensities between states, whichshould be a list of lists of intensity representations. As an intensity representation wemean a function that for a given Lexis object can be used to produce estimates of
the transition intensities at a set of supplied time points.
The names of the elements of the transition object (which are lists) will be names of
the transient states, that is the states from which a transition can occur. The namesof the elements of each of these lists are the names of states to which transitions canoccur (which may be either transient or absorbing states).
Hence, if the transition object is called Tr then TR$A$B (or Tr[["A"]][["B"]]) will
represent the transition intensity from state A to the state B.
The entries in the transition object can be either glm objects (either with poisson or
poisreg family), representing Poisson models for the transitions, coxph objects

4

1.3 Setting up a Lexis object

Multistate models with multiple timescales

representing an intensity model along one time scale, or simply a function that takesa Lexis object as input and returns an estimated intensity for each row.
In addition to these two input items, there will be a couple of tuning parameters.
The output of the function will simply be a Lexis object with simulated transitions
between states. This will be the basis for deriving sensible statistics from the Lexis object
 see next section.

1.3

Setting up a Lexis object

As an example we will use the DMlate dataset from the Epi package; it is a dataset
simulated to resemble a random sample of 10,000 patients from the Danish National
Diabetes Register.
We start by loading the Epi package:> options( width=90 )> library( Epi )> print( sessionInfo(), l=F )R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)
Matrix products: default
attached base packages:[1] statsgraphics grDevices utils

datasets

methodsbaseother attached packages:[1] Epi_2.52
loaded via a namespace (and not attached):[1] etm_1.1.1
vctrs_0.6.5[5] rlang_1.1.4
cmprsk_2.2-12[9] zoo_1.8-12
glue_1.7.0[13] grid_4.4.1
tibble_3.2.1[17] lifecycle_1.0.4
compiler_4.4.1[21] pkgconfig_2.0.3
mgcv_1.9-1[25] tidyselect_1.2.1
utf8_1.2.4[29] splines_4.4.1
magrittr_2.0.3

nlme_3.1-165
generics_0.1.3
plyr_1.8.9
MASS_7.3-61
dplyr_1.1.4
lattice_0.22-6
pillar_1.9.0
Matrix_1.7-0

cli_3.6.3
data.table_1.15.4
fansi_1.0.6
numDeriv_2016.8-1.1
Rcpp_1.0.12
R6_2.5.1
parallel_4.4.1
survival_3.7-0

First we load the diabetes data and set up a simple illness-death model:> data(DMlate)> dml <- Lexis( entry = list(Per=dodm, Age=dodm-dobth, DMdur=0 ),+exit = list(Per=dox),+
exit.status = factor(!is.na(dodth),labels=c("DM","Dead")),+data = DMlate )NOTE: entry.status has been set to "DM" for all.
NOTE: Dropping 4 rows with duration of follow up < tolThis is just data for a simple survival model with states DM and Dead. Now we cut thefollow-up at insulin start, which for the majority of patients (T2D) is a clinical indicator ofdeterioration of disease regulation. We therefore also introduce a new timescale, and split

Using simLexis

1.4 Analysis of rates

the non-precursor states, so that we can address the question of ever having been on
insulin:> dmi <- cutLexis( dml, cut = dml$doins,+
pre = "DM",+
new.state = "Ins",+
new.scale = "t.Ins",+
split.states = TRUE )> summary( dmi, timeScales=T )Transitions:To
From
DM Ins Dead Dead(Ins) Records: Events: Risk time: Persons:
DM 6157 1694 2048
09899
3742 45885.49
9899
Ins
0 1340
0
451
1791
451
8387.77
1791
Sum 6157 3034 2048
451
11690
4193 54273.27
9996
Timescales:
Per Age DMdur t.Ins"""""" "Ins"

Note that we show the time scales in the Lexis object, and that it is indicated that thetime scale t.Ins is dened as time since entry into stat state Ins.
We can show how many person-years we have and show the number of transitions andtransition rates (per 1000), using the boxes.Lexis function to display the states and thenumber of transitions between them:> boxes( dmi, boxpos = list(x=c(20,20,80,80),+
y=c(80,20,80,20)),+
scale.R = 1000, show.BE = TRUE )1.4

Analysis of rates

In the Lexis object (which is just a data frame) each person is represented by one record
for each transient state occupied, thus in this case either 1 or 2 records; those who have arecorded time both without and with insulin have two records.
In order to be able to t Poisson models with occurrence rates varying by the dierent
time-scales, we split the follow-up in 3-month intervals for modeling:> Si <- splitLexis( dmi, seq(0,20,1/4), "DMdur" )> summary( Si )Transitions:To
From
DM Ins Dead Dead(Ins) Records: Events: Risk time: Persons:
DM 184986 1694 2048
01887283742
45885.49
9899
Ins
0 34707
0
451
35158
451
8387.771791
Sum 184986 36401 2048
451
2238864193 54273.27
9996> print( subset( Si, lex.id==97 )[,1:10], digits=6 )lex.idPer
Age DMdur t.Ins lex.dur lex.Cst lex.Xst sex
dobth
97 1997.55 58.93 0.00
NA
0.25
DM
DM
F 1938.62
97 1997.80 59.18 0.25
NA
0.25
DM
DM
F 1938.62
97 1998.05 59.43 0.50
NA
0.25
DM
DM
F 1938.62
97 1998.30 59.68 0.75
NA
0.25
DM
DM
F 1938.62
97 1998.55 59.93 1.00
NA
0.25
DM
DM
F 1938.62

5

6

1.4 Analysis of rates

Multistate models with multiple timescalesDM
45,885.5
9,899
6,157

2,048
(44.6)
0

Dead
2,048

1,694
(36.9)

Ins
8,387.8
97
1,340

451
(53.8)

Dead(Ins)
0451

Figure 1.1: Data overview for the dmi dataset. Numbers in the boxes are person-years andthe number of persons who begin, resp. end their follow-up in each state, and numbers on thearrows are no. of transitions and rates (transition intensities) per 1000 PY../simLexis-boxes
97 1998.80 60.18
97 1999.05 60.43
97 1999.30 60.68
97 1999.55 60.93
97 1999.72 61.11
97 1999.80 61.18
97 2000.05 61.43
97 2000.30 61.68
97 2000.55 61.93

1.25
1.50
1.75
2.00
2.18
2.25
2.50
2.75
3.00

NA
NA
NA
NA
0.00
0.07
0.32
0.57
0.82

0.25
0.25
0.25
0.18
0.07
0.25
0.25
0.25
0.01

DM
DM
DM
DM
DM
DM
DM
Ins
Ins
Ins
Ins
Ins
Ins
Ins
Ins
Ins
Ins Dead(Ins)

F 1938.62
F 1938.62
F 1938.62
F 1938.62
F 1938.62
F 1938.62
F 1938.62
F 1938.62
F 1938.62

Note that when we split the follow-up, each person's follow up now consists of manyrecords, each with the current values of the timescales at the start of the interval
represented by the record. In the modeling we shall assume that the rates are constantwithin each 6-month interval, but the size of these rates we model as smooth functions of
the time scales (that is the values at the beginning of each interval).
The approach often used in epidemiology where one parameter is attached to eachinterval of time (or age) is not feasible when more than one time scale is used, becauseintervals are not classied the same way on all timescales.
We shall use natural splines (restricted cubic splines) for the analysis of rates, and hence
we must allocate knots for the splines. This is done for each of the time-scales, andseparately for the transition out of states DM and Ins. For age, we place the knots so thatthe number of events is the same between each pair of knots, but only half of this beyond
each of the boundary knots, whereas for the timescales DMdur and tIns where we haveobservation from a well-dened 0, we put knots at 0 and place the remaining knots so thatthe number of events is the same between each pair of knots as well as outside theUsing simLexis

1.4 Analysis of rates

boundary knots.
> nk <- 5> ( ai.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ),+quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) )
10%
30%
50%
70%
90%
23.23751 48.82218 58.63244 67.79028 78.88542> ( ad.kn <- with( subset(Si,lex.Xst=="Dead"),+quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) )
10%
30%
50%
70%
90%
61.91951 72.52731 78.43121 83.32348 90.15195> ( di.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ),+c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) )
20%
40%
60%
80%
0.00000000 0.06570842 0.45448323 3.28761123 6.63764545> ( dd.kn <- with( subset(Si,lex.Xst=="Dead"),+c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) )
20%
40%
60%
80%
0.0000000 0.7687885 2.1327858 4.0465435 6.5232033> ( ti.kn <- with( subset(Si,lex.Xst=="Dead(Ins)"),+c(0,quantile( t.Ins+lex.dur, probs=(1:(nk-1))/nk ) )) )
20%
40%
60%
80%
0.0000000 0.3093771 1.1307324 2.5489391 4.9117043

Note that when we tease out the event records for transition to transient states (in thiscase Ins, that is lex.Xst=="Ins"), we should add lex.Cst!=lex.Xst, to include onlytransition records and avoiding including records of sojourn time in the transient state.
We then t Poisson models to transition rates, using the wrapper Ns from the Epi
package to simplify the specication of the rates:> library( splines )> DM.Ins <- glm( (lex.Xst=="Ins") ~ Ns( Age , knots=ai.kn ) ++Ns( DMdur, knots=di.kn ) ++I(Per-2000) + sex,+family=poisson, offset=log(lex.dur),+data = subset(Si,lex.Cst=="DM") )> ci.exp( DM.Ins )exp(Est.)
2.5%
97.5%
(Intercept)
1.37516629 1.21414923 1.55753699
Ns(Age, knots = ai.kn)1 0.23411761 0.19185266 0.28569349
Ns(Age, knots = ai.kn)2 0.23166177 0.19576557 0.27414000
Ns(Age, knots = ai.kn)3 0.02835009 0.02284379 0.03518363
Ns(Age, knots = ai.kn)4 0.38067392 0.32945713 0.43985278
Ns(DMdur, knots = di.kn)1 0.04462626 0.03373451 0.05903459
Ns(DMdur, knots = di.kn)2 0.22388988 0.19028663 0.26342723
Ns(DMdur, knots = di.kn)3 0.03379141 0.02574708 0.04434907
Ns(DMdur, knots = di.kn)4 0.47100646 0.40783317 0.54396529
I(Per - 2000)
0.97381513 0.96082049 0.98698552
sexF
0.73757407 0.66886782 0.81333785> class( DM.Ins )[1] "glm" "lm"7

8

1.5 The mortality rates

Multistate models with multiple timescales

We can also t this model with a slightly simpler syntax using the glm.Lexis function:> DM.Ins <- glm.Lexis( Si, from = "DM", to = "Ins",+
formula = ~ Ns( Age , knots=ai.kn ) ++Ns( DMdur, knots=di.kn ) ++I(Per-2000) + sex )stats::glm Poisson analysis of Lexis object Si with log link:
Rates for the transition:
DM->Ins> ci.exp( DM.Ins )exp(Est.)
2.5%
97.5%
(Intercept)
1.37516630 1.21415038 1.55753552
Ns(Age, knots = ai.kn)1 0.23411761 0.19185214 0.28569426
Ns(Age, knots = ai.kn)2 0.23166177 0.19576473 0.27414118
Ns(Age, knots = ai.kn)3 0.02835009 0.02284374 0.03518371
Ns(Age, knots = ai.kn)4 0.38067392 0.32945601 0.43985427
Ns(DMdur, knots = di.kn)1 0.04462625 0.03373347 0.05903639
Ns(DMdur, knots = di.kn)2 0.22388988 0.19028599 0.26342813
Ns(DMdur, knots = di.kn)3 0.03379141 0.02574705 0.04434912
Ns(DMdur, knots = di.kn)4 0.47100646 0.40783201 0.54396683
I(Per - 2000)
0.97381513 0.96082027 0.98698574
sexF
0.73757407 0.66886641 0.81333956> class( DM.Ins )[1] "glm.lex" "glm""lm"

So we have a slightly simpler syntax, and we get an informative message of whichtransition(s) we are modeling. However we do not have update method for these objects.
> DM.Dead <- glm.Lexis( Si, from = "DM", to = "Dead",+formula = ~ Ns( Age , knots=ad.kn ) ++Ns( DMdur, knots=dd.kn ) ++I(Per-2000) + sex )stats::glm Poisson analysis of Lexis object Si with log link:
Rates for the transition:
DM->Dead> Ins.Dead <- glm.Lexis( Si, from = "Ins",+
formula = ~ Ns( Age , knots=ad.kn ) ++Ns( DMdur, knots=dd.kn ) ++Ns( t.Ins, knots=ti.kn ) ++I(Per-2000) + sex )stats::glm Poisson analysis of Lexis object Si with log link:
Rates for the transition:Ins->Dead(Ins)

Note the similarity of the code used to t the three models, is is mainly redening the
response variable (to state) and the subset of the data used (from state). Also note thatthe last model need no specication of to, the default is to model all transitions from thefrom state, and his case there is only one.

1.5

The mortality rates

This section discusses in some detail how to extract ad display the mortality rates from themodels tted. But it is not necessary for understanding how to use simLexis in practice.

Using simLexis

1.5 The mortality rates

91.5.1 Proportionality of mortality rates
Note that we have tted separate models for the three transitions, there is no assumptionof proportionality between the mortality rates from DM and Ins.
However, there is nothing that prevents us from testing this assumption; we can just t amodel for the mortality rates in the entire data frame Si, and compare the deviance fromthis with the sum of the deviances from the separate models using the glm.Lexis function:> All.Dead <- glm.Lexis( Si, to = c("Dead(Ins)","Dead"),+formula = ~ Ns( Age , knots=ad.kn ) ++Ns( DMdur, knots=dd.kn ) ++lex.Cst ++I(Per-2000) + sex )stats::glm Poisson analysis of Lexis object Si with log link:
Rates for transitions:Ins->Dead(Ins)
DM->Dead> round( ci.exp( All.Dead ), 3 )(Intercept)
Ns(Age, knots = ad.kn)1
Ns(Age, knots = ad.kn)2
Ns(Age, knots = ad.kn)3
Ns(Age, knots = ad.kn)4
Ns(DMdur, knots = dd.kn)1
Ns(DMdur, knots = dd.kn)2
Ns(DMdur, knots = dd.kn)3
Ns(DMdur, knots = dd.kn)4
lex.CstIns
I(Per - 2000)
sexF

exp(Est.) 2.5% 97.5%
0.057 0.049 0.065
4.101 3.462 4.858
4.661 4.064 5.346
15.434 13.548 17.583
7.509 6.695 8.421
0.466 0.384 0.565
0.642 0.563 0.731
0.229 0.165 0.318
0.796 0.713 0.888
2.168 1.947 2.415
0.965 0.954 0.977
0.665 0.614 0.721

Incidentally we could have dispensed with the to= argument too, because the default is totake to to be all absorbing states in the model.
From the parameter values we would in a simple setting just claim that start ofinsulin-treatment was associated with a slightly more than doubling of mortality.
The model All.dead assumes that the age- and DM-duration eects on mortality in the
DM and Ins states are the same, and moreover that there is no eect of insulin duration,
but merely a mortality that is larger by a multiplicative constant not depending on insulin
duration. The model DM.Dead has 8 parameters to describe the dependency on age andDM duration, the model Ins.Dead has 12 for the same plus the insulin duration (a natural
spline with k knots gives k − 1 parameters, and we chose k = 5 above).
We can compare the t of the simple proportional hazards model with the t of theseparate models for the two mortality rates, by adding up the deviances and d.f. fromthese:> what <- c("null.deviance","df.null","deviance","df.residual")> ( rD <- unlist( DM.Dead[what] ) )null.deviancedf.nulldeviance df.residual
22535.77
188727.00
20412.81
188717.00> ( rI <- unlist( Ins.Dead[what] ) )10

1.5 The mortality rates

null.deviance4867.127

df.null
35157.000

Multistate models with multiple timescales
deviance
4211.735

df.residual
35143.000> ( rA <- unlist( All.Dead[what] ) )null.deviance27415.21

df.null
223885.00

deviance
24705.70

df.residual
223874.00> round( c( dd <- rA-(rI+rD), "pVal"=1-pchisq(dd[3],dd[4]+1) ), 3 )null.deviance12.314

df.null
1.000

deviance
81.154

df.residual pVal.deviance14.0000.000

Thus we see there is a substantial non-proportionality of mortality rates from the twostates; but a test provides no clue whatsoever to the particular shape of the
non-proportionality.
To this end, we shall explore the predicted mortalities under the two models
quantitatively in more detail. Note that the reason that there is a dierence in the nulldeviances (and a dierence of 1 in the null d.f.) is that the null deviance of All.Dead referto a model with a single intercept, that is a model with constant and identical mortality
rates from the states DM and Ins, whereas the null models for DM.Dead and Ins.Dead haveconstant but dierent mortality rates from the states DM and Ins. This is however
irrelevant for the comparison of the residual deviances.

1.5.2 How the mortality rates look
If we want to see how the mortality rates are modelled in DM.Dead and Ins.Dead in
relation to All.Dead, we make a prediction of rates for say men diagnosed in dierent ages
and going on insulin at dierent times after this. So we consider men diagnosed in ages 40,
50, 60 and 70, and who either never enter insulin treatment or do it 0, 2 or 5 years after
diagnosis of DM.
To this end we create a prediction data frame where we have observation times fromdiagnosis and 12 years on (longer would not make sense as this is the extent of the data).
But we start by setting up an array to hold the predicted mortality rates, classied bydiabetes duration, age at diabetes onset, time of insulin onset, and of course type of model.
What we want to do is to plot the age-specic mortality rates for persons not on insulin,
and for persons starting insulin at dierent times after DM. The mortality curves start at
the age where the person gets diabetes and continues 12 years; for persons on insulin they
start at the age when they initiate insulin.
> pr.rates <- NArray( list( DMdur = seq(0,12,0.1),+
DMage = 4:7*10,+
r.Ins = c(NA,0,2,5),+
model = c("DM/Ins","All"),+what = c("rate","lo","hi") ) )> str( pr.rates )logi [1:121, 1:4, 1:4, 1:2, 1:3] NA NA NA NA NA NA ...
- attr(*, "dimnames")=List of 5..$ DMdur: chr [1:121] "0" "0.1" "0.2" "0.3" ...
..$ DMage: chr [1:4] "40" "50" "60" "70"..$ r.Ins: chr [1:4] NA "0" "2" "5"..$ model: chr [1:2] "DM/Ins" "All"..$ what : chr [1:3] "rate" "lo" "hi"Using simLexis

1.5 The mortality rates

11

For convenience the Epi package contains a function that computes predicted (log-)rates
with c.i.  it is merely a wrapper for predict.glm.
So we set up the prediction data frame and modify it in loops over ages at onset andinsulin onset in order to collect the predicted rates in dierent scenarios:> nd <- data.frame( DMdur = as.numeric( dimnames(pr.rates)[[1]] ),+lex.Cst = factor( 1, levels=1:4,+
labels=levels(Si$lex.Cst) ),+sex = factor( 1, levels=1:2, labels=c("M","F")) )Note that we did not insert lex.dur as covariate in the prediction frame. This would be
required if we used the poisson family with the glm, but the wrapper glm.Lexis uses thepoisreg family, so lex.dur is ignored and predictions always comes in the (inverse) unitsof lex.dur. So we get rates per 1 person-year in the predictions.
> for( ia in dimnames(pr.rates)[[2]] )+{+ dnew <- transform( nd, Age = as.numeric(ia)+DMdur,+
Per = 1998+DMdur )+ pr.rates[,ia,1,"DM/Ins",] <- ci.pred( DM.Dead, newdata = dnew )+ pr.rates[,ia,1,"All" ,] <- ci.pred( All.Dead, newdata = dnew )+ for( ii in dimnames(pr.rates)[[3]][-1] )+{+ dnew = transform( dnew, lex.Cst = factor( 2, levels=1:4,+
labels=levels(Si$lex.Cst) ),+t.Ins = ifelse( (DMdur-as.numeric(ii)) >= 0,+
DMdur-as.numeric(ii), NA ) )+ pr.rates[,ia, ii ,"DM/Ins",] <- ci.pred( Ins.Dead, newdata = dnew )+ pr.rates[,ia, ii ,"All",] <- ci.pred( All.Dead, newdata = dnew )+}+}So for each age at DM onset we make a plot of the mortality as function of current age
both for those with no insulin treatment and those that start insulin treatment 0, 2 and 5
years after diabetes diagnosis, thus 4 curves (with c.i.). These curves are replicated with adierent color for the simplied model.
> par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 )> plot( NA, xlim=c(40,82), ylim=c(5,300), bty="n",+log="y", xlab="Age", ylab="Mortality rate per 1000 PY" )> abline( v=seq(40,80,5), h=outer(1:9,10^(0:2),"*"), col=gray(0.8) )> for( aa in 4:7*10 ) for( ii in 1:4 )+matshade( aa+as.numeric(dimnames(pr.rates)[[1]]),+cbind( pr.rates[,paste(aa),ii,"DM/Ins",],+pr.rates[,paste(aa),ii,"All" ,] )*1000,+
type="l", lty=1, lwd=2,+
col=c("red","limegreen") )From gure 1.2 we see that there is a substantial insulin-duration eect which is notaccommodated by the simple model with only one time-dependent variable to describe theinsulin eect. Note that the simple model (green curves) for those on insulin does notdepend in insulin duration, and hence the mortality curves for those on insulin are justparallel to the mortality curves for those not on insulin, regardless of diabetes duration (orage) at the time of insulin initiation. This is the proportional hazards assumption. Thus
the eect of insulin initiation is under-estimated for short duration of insulin andover-estimated for long duration of insulin.

12

1.6 Input to the simLexis function

Multistate models with multiple timescales

200

Mortality rate per 1000 PY

100

50

20

10

5
40

50

60
Age

70

80

Figure 1.2: Estimated mortality rates for male diabetes patients with no insulin (lower setsof curves) and insulin (upper curves), with DM onset in age 40, 50, 60 and 70. The red
curves are from the models with separate age eects for persons with and without insulin,
and a separate eect of insulin duration. The green curves are from the model with common
age-eects and only a time-dependent eect of insulin, assuming no eect of insulin duration(the classical time-dependent variable approach). Hence the upper green curve is common
for any time of insulin inception.
./simLexis-mort-int
This is the major discrepancy between the two models, and illustrates the importance ofbeing able to accommodate dierent time scales, but there is also a declining overall insulineect by age which is not accommodated by the proportional hazards approach.
Finally, this plot illustrates an important feature in reporting models with multiple
timescales; all timescales must be represented in the predicted rates, only reporting theeect of one timescale, conditional on a xed value of other timescales is misleading since
all timescales by denition advance at the same pace. For example, the age-eect for axed value of insulin duration really is a misnomer since it does not correspond to any real
person's follow-up, but to the mortality of persons in dierent ages but with the same
duration of insulin use.

1.6

Input to the simLexis function

We want to estimate the cumulative probability of being in each of the 4 states, so that we
can assess the fraction of diabetes pateints that go on insulin

Using simLexis

1.6 Input to the simLexis function

13

In order to simulate from the multistate model with the estimated transition rates, andget the follow-up of a hypothetical cohort, we must supply both the transition rates and thestructure of the model as well as the initial cohort status to simLexis.

1.6.1 The transition object
We rst put the models into an object representing the transitions; note this is a list oflists, the latter having glm objects as elements:> Tr <- list( "DM" = list( "Ins"= DM.Ins,+"Dead"= DM.Dead ),+"Ins" = list( "Dead(Ins)" = Ins.Dead ) )Now we have the description of the rates and of the structure of the model. The Tr object
denes the states and models for all transitions between them; the object Tr$A$B is themodel for the transition intensity from state A to state B.

1.6.2 The initial cohort
We now dene an initial Lexis object of persons with all relevant covariates dened. Notethat we use NULL as row indicator in the Lexis object we used for modeling; this conserves
the time.scale and time.since attributes which are needed for the simulation:> str( ini <- Si[NULL,1:9] )Classes 'Lexis' and 'data.frame':
0 obs. of 9 variables:$ lex.id : int
$ Per
: num
$ Age
: num
$ DMdur : num$ t.Ins : num$ lex.dur: num$ lex.Cst: Factor w/ 4 levels "DM","Ins","Dead",..:$ lex.Xst: Factor w/ 4 levels "DM","Ins","Dead",..:$ sex: Factor w/ 2 levels "M","F":- attr(*, "time.scales")= chr [1:4] "Per" "Age" "DMdur" "t.Ins"- attr(*, "time.since")= chr [1:4] "" "" "" "Ins"- attr(*, "breaks")=List of 4..$ Per : NULL..$ Age : NULL..$ DMdur: num [1:81] 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 ...
..$ t.Ins: NULL

We now have an empty Lexis object with attributes reecting the timescales in themultistate model we want to simulate from. But we must enter some data to represent theinitial state of the persons whose follow-up we want to simulate through the model; so llin data for one man and one woman:> ini[1:2,"lex.id"] <- 1:2> ini[1:2,"lex.Cst"] <- "DM"> ini[1:2,"Per"] <- 1995
> ini[1:2,"Age"] <- 60> ini[1:2,"DMdur"] <- 5
> ini[1:2,"sex"] <- c("M","F")
> ini

14

1.7 Simulation of the follow-up

Multistate models with multiple timescales

lex.id Per Age DMdur t.Ins lex.dur lex.Cst lex.Xst sex
1 1995 60
5
NA
NA
DM
<NA>
M
2 1995 60
5
NA
NA
DM
<NA>
F

So the persons starts in age 60 in 1995 with 5 years of diabetes duration. Note that thet.Ins is NA, because this is a timescale that rst comes alive if a transtion to Ins issimulated.

1.7

Simulation of the follow-up

Now we simulate life-courses of a 1000 of each of these persons using the estimated model.
The t.range argument gives the times range where the integrated intensities (cumulative
rates) are evaluated and where linear interpolation is used when simulating transition
times. Note that this must be given in the same units as lex.dur in the Lexis object usedfor tting the models for the transitions. It is not a parameter that can be easily
determined from the TR object, hence it must be supplied by the user.
> set.seed(52381764)> Nsim <- 500> system.time( simL <- simLexis( Tr,+ini,+t.range = 12,+
N = Nsim ) )bruger system forløbet
1.22
0.31
1.53

The result is a Lexis object  a data frame representing the simulated follow-up of 1000
persons (500 identical men and 500 identical women) according to the rates we estimatedfrom the original dataset.
> summary( simL, by="sex" )$M
Transitions:To
From DM Ins Dead Dead(Ins) Records: Events: Risk time: Persons:
DM 148 185 167
0
500
352
3697.59500
Ins 0 136
0
49
185
49
1023.96185
Sum 148 321 167
49
685
401
4721.55
500$FTransitions:To
From DM Ins Dead Dead(Ins) Records: Events: Risk time: Persons:
DM 214 183 103
0
500
286
4223.33500
Ins 0 144
0
39
183
39
1020.36183
Sum 214 327 103
39
683
325
5243.69500

Using simLexis

1.7 Simulation of the follow-up

15

1.7.1 Using other models for simulation
Proportional hazards Poisson modelWe tted a proportional mortality model All.Dead (which tted worse than the othertwo), this is a model for both the transition from DM to Death and from Ins to Dead(Ins),
assuming that they are proportional. But it can easily be used in the simulation set-up,
because the state is embedded in the model via the term lex.Cst, which is updated during
the simulation.
Simulation using this instead just requires that we supply a dierent transition object:> Tr.p <- list( "DM" = list( "Ins"= DM.Ins,+"Dead"= All.Dead ),+"Ins" = list( "Dead(Ins)" = All.Dead ) )> system.time( simP <- simLexis( Tr.p,+ini,+t.range = 12,+
N = Nsim ) )bruger system forløbet
1.33
0.10
1.41
> summary( simP, by="sex" )$M
Transitions:To
From DM Ins Dead Dead(Ins) Records: Events: Risk time: Persons:
DM 148 222 130
0
500
352
3703.19500
Ins 0 127
0
95
222
95
1049.95
222
Sum 148 349 130
95
722
447
4753.14500$FTransitions:To
From DM Ins Dead Dead(Ins) Records: Events: Risk time: Persons:
DM 224 171 105
0
500
276
4323.34500
Ins 0 118
0
53
171
53
861.93
171
Sum 224 289 105
53
671
329
5185.27500

Proportional hazards Cox modelA third possibility would be to replace the two-time scale proportional mortality model by
a one-time-scale Cox-model, using diabetes duration as time scale, and age at diagnosis of
DM as (xed) covariate:> library( survival )> Cox.Dead <- coxph( Surv( DMdur, DMdur+lex.dur,+
lex.Xst %in% c("Dead(Ins)","Dead")) ~+Ns( Age-DMdur, knots=ad.kn ) ++I(lex.Cst=="Ins") ++I(Per-2000) + sex,+
data = Si )> round( ci.exp( Cox.Dead ), 3 )16

1.8 Reporting the simulation results

Ns(Age - DMdur, knots = ad.kn)1
Ns(Age - DMdur, knots = ad.kn)2
Ns(Age - DMdur, knots = ad.kn)3
Ns(Age - DMdur, knots = ad.kn)4I(lex.Cst == "Ins")TRUE
I(Per - 2000)
sexF

Multistate models with multiple timescales

exp(Est.) 2.5% 97.5%
4.172 3.535 4.923
4.502 3.824 5.301
16.077 14.087 18.348
7.479 6.501 8.605
2.170 1.948 2.418
0.965 0.954 0.977
0.667 0.616 0.723

Note that in order for this model to be usable for simulation, it is necessary that we use
the components of the Lexis object to specify the survival. Each record in the data frame
Si represents follow up from DMdur to DMdur+lex.dur, so the model is a Cox model withdiabetes duration as underlying timescale and age at diagnosis, Age-DMdur, as covariate.
Also note that we used I(lex.Cst=="Ins") instead of just lex.Cst, because coxph
assigns design matrix columns to all levels of lex.Cst, also those not present in data,
which would give NAs among the parameter estimates and NAs as mortality outcomes.
We see that the eect of insulin and the other covariates are pretty much the same as in
the two-time-scale model. We can simulate from this model too; there is no restrictions onwhat type of model can be used for dierent transitions> Tr.c <- list( "DM" = list( "Ins"= Tr$DM$Ins,+"Dead"= Cox.Dead ),+"Ins" = list( "Dead(Ins)" = Cox.Dead ) )> system.time( simC <- simLexis( Tr.c,+ini,+t.range = 12,+
N = Nsim ) )bruger system forløbet
2.46
0.34
2.81
> summary( simC, by="sex" )$M
Transitions:To
From DM Ins Dead Dead(Ins) Records: Events: Risk time: Persons:
DM 171 192 137
0
500
329
3645.74500
Ins 0 123
0
69
192
69
1007.35
192
Sum 171 315 137
69
692
398
4653.09500$F
Transitions:To
From DM Ins Dead Dead(Ins) Records: Events: Risk time: Persons:
DM 227 185 88
0
500
273
4332.82500
Ins 0 138
0
47
185
47
821.79185
Sum 227 323 88
47
685
320
5154.61500

1.8Reporting the simulation results

We can now tabulate the number of persons in each state at a predened set of times on agiven time scale. Note that in order for this to be sensible, the from argument would
normally be equal to the starting time for the simulated individuals.

Using simLexis

1.8 Reporting the simulation results

17> system.time(+ nSt <- nState( subset(simL,sex=="M"),+
at=seq(0,11,0.2), from=1995, time.scale="Per" ) )bruger system forløbet
0.07
0.00
0.06
> nSt[1:10,]
State
when
DM Ins Dead Dead(Ins)
1995 500 0001995.2 497
2
1
0
1995.4 486
8
60
1995.6 477 13 1001995.8 468 16 1601996 459 20
21
01996.2 448 25 27
01996.4 442 29 28
1
1996.6 438 31 30
1
1996.8 429 37 33
1

We see that as time goes by, the 500 men slowly move away from the starting state (DM).
Based on this table (nSt is a table) we can now compute the fractions in each state, or,rather more relevant, the cumulative fraction across the states in some specied order, sothat a plot of the stacked probabilities can be made, using either the default rather colorfullayout, or a more minimalist version (both in gure 1.3):> pM <- pState( nSt, perm=c(1,2,4,3) )> head( pM )
State
when
DM Ins Dead(Ins) Dead
1995 1.000 1.000
1.000
11995.2 0.994 0.998
0.998
1
1995.4 0.972 0.988
0.988
11995.6 0.954 0.980
0.980
11995.8 0.936 0.968
0.968
11996 0.918 0.958
0.958
1> par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )> plot( pM )> plot( pM, border="black", col="transparent", lwd=3 )> text( rep(as.numeric(rownames(pM)[nrow(pM)-1]),ncol(pM)),+
pM[nrow(pM),]-diff(c(0,pM[nrow(pM),]))/5,+
colnames( pM ), adj=1 )> box( col="white", lwd=3 )> box()A more useful set-up of the graph would include a more through annotation and sensible
choice of colors, as seen in gure 1.4:
> clr <- c("limegreen","orange")> # expand with a lighter version of the two chosen colors> clx <- c( clr, rgb( t( col2rgb( clr[2:1] )*2 + rep(255,3) ) / 3, max=255 ) )> par( mfrow=c(1,2), las=1, mar=c(3,3,4,2), mgp=c(3,1,0)/1.6 )> # Men> plot( pM, col=clx, xlab="Date of FU" )> lines( as.numeric(rownames(pM)), pM[,2], lwd=3 )> mtext( "60 year old male, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6)

18

1.8 Reporting the simulation results

Multistate models with multiple timescales

1.0

1.0
Dead

0.8

0.8

Dead(Ins)
0.6

Ins

Probability

Probability

0.6

0.4

0.4

DM
0.2

0.2

0.0

0.0
1996

1998

2000
2002
Time

2004

2006

1996

1998

2000
2002
Time

2004

2006

Figure 1.3: Default layout of the plot.pState graph (left), and a version with the state

probabilities as lines and annotation of states.

./simLexis-pstate0> mtext( "Survival curve", side=3, line=1.5, adj=0 )> mtext( "DM, no insulin
DM, Insulin", side=3, line=0.5, adj=0, col=clr[2] )> mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[1] )> axis( side=4 )> axis( side=4, at=1:19/20, labels=FALSE )> axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 )> # Women> pF <- pState( nState( subset(simL,sex=="F"),+at=seq(0,11,0.2),+
from=1995,+
time.scale="Per" ),+perm=c(1,2,4,3) )> plot( pF, col=clx, xlab="Date of FU" )> lines( as.numeric(rownames(pF)), pF[,2], lwd=3 )> mtext( "60 year old female, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6> mtext( "Survival curve", side=3, line=1.5, adj=0 )> mtext( "DM, no insulin
DM, Insulin", side=3, line=0.5, adj=0, col=clr[2] )> mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[1] )> axis( side=4 )> axis( side=4, at=1:19/20, labels=FALSE )> axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 )If we instead wanted to show the results on the age-scale, we would use age as timescale
when constructing the probabilities; otherwise the code is pretty much the same as before(Figure 1.5):

Using simLexis

1.8 Reporting the simulation results
60 year old female, diagnosed 1990, aged 55
Survival curveDM, no insulin DM, Insulin

1.0

1.0

1.0

1.0

0.8

0.8

0.8

0.8

0.6

0.6

0.6

0.6

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0.2

0.0
2006

0.0

0.0
1996

1998

2000
2002
Date of FU

2004

Probability

Probability

60 year old male, diagnosed 1990, aged 55
Survival curveDM, no insulin DM, Insulin

19

1996

1998

2000
2002
Date of FU

2004

0.0
2006

Figure 1.4: plot.pState graphs where persons ever on insulin are given in orange andpersons never on insulin in green, and the overall survival (dead over the line) as a blackline.
./simLexis-pstatex> par( mfrow=c(1,2), las=1, mar=c(3,3,4,2), mgp=c(3,1,0)/1.6 )> # Men> pM <- pState( nState( subset(simL,sex=="M"),+
at=seq(0,11,0.2),+
from=60,+
time.scale="Age" ),+perm=c(1,2,4,3) )> plot( pM, col=clx, xlab="Age" )> lines( as.numeric(rownames(pM)), pM[,2], lwd=3 )> mtext( "60 year old male, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6)> mtext( "Survival curve", side=3, line=1.5, adj=0 )> mtext( "DM, no insulin
DM, Insulin", side=3, line=0.5, adj=0, col=clr[2] )> mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[1] )> axis( side=4 )> axis( side=4, at=1:19/20, labels=FALSE )> axis( side=4, at=1:19/20, labels=FALSE, tcl=-0.4 )> axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 )> # Women> pF <- pState( nState( subset(simL,sex=="F"),+at=seq(0,11,0.2),+
from=60,+
time.scale="Age" ),+perm=c(1,2,4,3) )> plot( pF, col=clx, xlab="Age" )> lines( as.numeric(rownames(pF)), pF[,2], lwd=3 )20

1.8 Reporting the simulation results

Multistate models with multiple timescales> mtext( "60 year old female, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6> mtext( "Survival curve", side=3, line=1.5, adj=0 )> mtext( "DM, no insulin
DM, Insulin", side=3, line=0.5, adj=0, col=clr[2] )> mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[1] )> axis( side=4 )> axis( side=4, at=1:9/10, labels=FALSE )> axis( side=4, at=1:19/20, labels=FALSE, tcl=-0.4 )> axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 )Note the several statements with axis(side=4,...; they are necessary to get the netick-marks in the right hand side of the plots that you will need in order to read o theprobabilities at 2006 (or 71 years).
60 year old female, diagnosed 1990, aged 55
Survival curveDM, no insulin DM, Insulin

1.0

1.0

1.0

1.0

0.8

0.8

0.8

0.8

0.6

0.6

0.6

0.6

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0.2

0.0

0.0

0.0
60

62

64

66
Age

68

70

Probability

Probability

60 year old male, diagnosed 1990, aged 55
Survival curveDM, no insulin DM, Insulin

0.0
60

62

64

66
Age

68

70

Figure 1.5: plot.pState graphs where persons ever on insulin are given in orange andpersons never on insulin in green, and the overall survival (dead over the line) as a blackline.
./simLexis-pstatey

1.8.1 Comparing predictions from dierent models
We have actually tted dierent models for the transitions, and we have simulated Lexis
objects from all three approaches, so we can plot these predictions on top of each other:> PrM <- pState( nState( subset(simP,sex=="M"),+
at=seq(0,11,0.2),+
from=60,+
time.scale="Age" ),+perm=c(1,2,4,3) )Using simLexis

1.8 Reporting the simulation results

21> PrF <- pState( nState( subset(simP,sex=="F"),+at=seq(0,11,0.2),+
from=60,+
time.scale="Age" ),+perm=c(1,2,4,3) )> CoxM <- pState( nState( subset(simC,sex=="M"),+
at=seq(0,11,0.2),+
from=60,+
time.scale="Age" ),+perm=c(1,2,4,3) )> CoxF <- pState( nState( subset(simC,sex=="F"),+at=seq(0,11,0.2),+
from=60,+
time.scale="Age" ),+perm=c(1,2,4,3) )> par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )> plot( pM, border="black", col="transparent", lwd=3 )> lines( PrM, border="blue" , col="transparent", lwd=3 )> lines( CoxM, border="red" , col="transparent", lwd=3 )> text( 60.5, 0.05, "M" )> box( lwd=5, col="white" ) ; box( lwd=2, col="black" )> plot( pF, border="black", col="transparent", lwd=3 )> lines( PrF, border="blue" , col="transparent", lwd=3 )> lines( CoxF, border="red" , col="transparent", lwd=3 )> text( 60.5, 0.05, "F" )> box( lwd=5, col="white" ) ; box( lwd=2, col="black" )From gure 1.6 it is clear that the two proportional hazards models (blue and red
curves) produce pretty much the same estimates of the state occupancy probabilities overtime, but also that they relative to the model with separately estimated transition
intensities overestimates the probability of being alive without insulin and underestimates
the probabilities of being dead without insulin. However both the overall survival, and thefraction of persons on insulin are quite well in agreement with the more elaborate model.
Thus the proportional hazards models overestimate the relative mortality of the insulin
treated diabetes patients relative to the non-insulin treated.
Interestingly, we also see a bump in the estimated probabilities from the Cox-model
based model, but this is entirely an artifact that comes from the estimation method for thebaseline hazard of the Cox-model that lets the (cumulative) hazard have large jumps at
event times where the risk set is small. So also here it shows up that an assumption ofcontinuous underlying hazards leads to more credible estimates.

22

1.8 Reporting the simulation results

Multistate models with multiple timescales

0.8

0.8

0.6

0.6
Probability

1.0

Probability

1.0

0.4

0.4

0.2

0.2

M

F

0.0

0.0
60

62

64

66
Time

68

70

60

62

64

66
Time

68

70

Figure 1.6: Comparison of the simulated state occupancy probabilities using separate Poissonmodels for the mortality rates with and without insulin (black) and using proportional hazards
Poisson models (blue) and Cox-models with diabetes duration as timescale and age at diabetes
diagnosis as covariate (red).
./simLexis-comp-0

Chapter 2
Simulation of transitions in multistatemodels2.1

TheorySuppose that the rate functions for the transitions out of the current state to, say, 3dierent states are λ1 , λ2 and λ3 , and the corresponding cumulative rates are Λ1 , Λ2 andΛ3 , and we want to simulate an exit time and an exit state (that is either 1, 2 or 3). Thiscan be done in two slightly dierent ways:
1. First time, then state:(a) Compute the survival function, S(t) = exp −Λ1 (t) − Λ2 (t) − Λ3 (t)(b) Simulate a random U(0,1) variate, u, say.
(c) The simulated exitPtime is then the solution tu to the equation
S(tu ) = u ⇔j Λj (tu ) = − log(u).
(d) A simulated transition at tu is then found by simulating a randomdraw fromP
the multinomial distribution with probabilities pi = λi (tu )/ j λj (tu ).
2. Separate cumulative incidences:(a) Simulate 3 independent U(0,1) random variates u1 , u2 and u3 .
(b) Solve the equations Λi (ti ) = − log(ui ), i = 1, 2, 3 and get (t1 , t2 , t3 ).
(c) The simulated survival time is then min(t1 , t2 , t3 ), and the simulated transition isto the state corresponding to this, that is k ∈ {1, 2, 3}, where tk = min(t1 , t2 , t3 )The intuitive argument is that with three possible transition there are 3 independent
processes running, but only the rst transition is observed. The latter approach is used in
the implementation in simLexis.
The formal argument for the equality of the two approaches goes as follows:
1. Equality of the transition times:23

24

2.1 Theory

Multistate models with multiple timescales(a) In the rst approach we simulate from a distribution with cumulative rate
Λ1 (t) + Λ2 (t) + Λ3 (t), hence from a distribution with survival function:S(t) = exp −(Λ1 (t) + Λ2 (t) + Λ3 (t))

= exp −Λ1 (t) × exp −Λ2 (t) × exp −Λ3 (t)(b) In the second approach we choose the smallest of three independent survival
times, with survival functions exp(−Λi ), i = 1, 2, 3. Now, the survival function
for the minimum of three independent survival times is:

Smin (t) = P{min(t1 , t2 , t3 ) > t}
= P{t1 > t} × P{t2 > t} × P{t3 > t}
= exp −Λ1 (t) × exp −Λ2 (t) × exp −Λ3 (t)which is the same survival function as derived above.
2. Type of transition:(a) In the rst instance the probability of a transition to state i, conditional on thetransition time being t, is as known from standard probability theory:
λi (t)/ λ1 (t) + λ2 (t) + λ3 (t) .
(b) In the second approach we choose the transition corresponding to the thesmallest of the transition times. So when we condition on the event that atransition takes place at time t, we have to show that the conditional probability
that the smallest of the three simulated transition times was actually the ith, is
as above.
But conditional on survival till t, the probabilities that events of type 1, 2, 3
takes place in the interval (t, t + dt) are λ1 (t) dt, λ2 (t) dt and λ1 (t) dt,
respectively (assuming that the probability of more than one event in theinterval of length dt is 0). Hence the conditional probabilities given a transition
time in (t, t + dt) is:λi (t)λi (t) dt=
λ1 (t) dt + λ2 (t) dt + λ3 (t) dt
λ1 (t) + λ2 (t) + λ3 (t) exactly as above.
Start time: 2024-07-18, 15:26:45
End time: 2024-07-18, 15:27:10
Elapsed time: 0.41 minutes

Bibliography[1] B Carstensen and M Plummer. Using Lexis objects for multi-state models in R.
Journal of Statistical Software, 38(6):118, 1 2011.
[2] S Iacobelli and B Carstensen. Multiple time scales in multi-state models. Stat Med,
32(30):53155327, Dec 2013.
[3] M Plummer and B Carstensen. Lexis: An R class for epidemiological studies withlong-term follow-up. Journal of Statistical Software, 38(5):112, 1 2011.

25