Competing risks with
Lexis, parametric rates andsimulation based condence intervals
SDCC
April 2024

http://bendixcarstensen.com/
Version 5th July, 2024, 15:26
from: C:\Bendix\stat\r\lib.src\epi\pkg\vignettes/crisk.tex
Compiled Thursday 18

Bendix Carstensen

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

1 Competing risks
1.1

1.2

1

Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1.1

Cause specic rates and likelihood . . . . . . . . . . . . . . . . . . . .

1

1.1.2

Survival and cumulative risks. . . . . . . . . . . . . . . . . . . . . .

1

1.1.3

Sojourn times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.1.4

The time scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2Rates and cumulative risks . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2.1

Condence intervals by simulation . . . . . . . . . . . . . . . . . . . .

3

1.2.2

Subdistribution hazard . . . . . . . . . . . . . . . . . . . . . . . . . .

32 Example data4

2.1

A Lexis object. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

Models for rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.3

Cumulative rates and risks . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.4

Sojourn times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

3 Condence intervals

12

3.1

Common parameters across cause-specic rates. . . . . . . . . . . . . . . .

13

3.2

Simulation based condence intervals . . . . . . . . . . . . . . . . . . . . . .

13

3.3

Simulated condence intervals for rates . . . . . . . . . . . . . . . . . . . . .

15

3.4

Condence intervals for cumulative risks. . . . . . . . . . . . . . . . . . . .

17

3.5

Condence intervals for stacked cumulative risks . . . . . . . . . . . . . . . .

18

3.6

Sojourn times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

4 A simple illustration

20

4.1

Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

4.2

A Lexis object with 3 causes of death. . . . . . . . . . . . . . . . . . . . .

21

4.3

Models for the rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

4.4

Derived measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

4.4.1

Cumulative risks. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

4.4.2

Stacked cumulative risks . . . . . . . . . . . . . . . . . . . . . . . . .

23

4.4.3

Sojourn times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

4.4.4

Comparing groups

26. . . . . . . . . . . . . . . . . . . . . . . . . . . .

ii

Chapter 1
Competing risks
1.1 Concepts
The concept of competing risks is one where persons in a given state, alive, say, aresubject to a number of dierent causes of death, cause1, cause2 etc. Causes of death arerequired to be exhaustive and mutually exclusive. That is, you will eventually die from oneof the designated causes, and you can only die from one.
The

cause-specic rate for cause c is dened as:

λc (t) = P{death from cause c in (t, t + h] | alive at t} /h. . . formally, the limit of this quantity as h → 0.
The observed data will be a survival time and an exit status which is either censoredalive or one of the causes. In situations where the causes are not causes of death but otherevents, it is implicit that we only consider the rst occurrence of an event from the statealive, and ignore whatever occurs after that.

1.1.1 Cause specic rates and likelihood
The likelihood for observations from a competing risk scenario is a function of thecause-specic transition rates, and is a

product of the likelihoods that would emerge if we

considered each cause as being the only possible event. Thus analysis is in principle
straight forward: estimate a model for each of the cause-specic rates; these will together
form a complete model for the competing risks problem.
If the cause-specic rates are all we want to assess then we are done.

1.1.2 Survival and cumulative risks
In addition to the rates we might however also be interested in thethecumulative risks of each cause of death.

survival probability and

The survival is the probability of still being alive at a given time after entrya functionof time since entry. The cumulative risk of dying from cause c is the probability of having
died from cause c as a function of time since entry.
This means that a time of entry is required for the calculations these quantities.

1

2

1.2

Rates and cumulative risks

CmpRskParSim

1.1.3 Sojourn times
sojourn time for cause c is the time spent in the cause c state before a given time, t,say. This is also called the expected lifetime lost to cause c, truncated at the time t. For
The

the state alive it will be the expected time lived before t.

1.1.4 The time scale
The cause specic rates will be functions of come covariates, notably a time scale, be thatage or time since entry to the study or even calendar time. But the cumulative risks areprobabilities that refer to timesince some origin. Thus cumulative risks (and survival) areonly meaningful in relation to a time that begins at 0. Though not a formal mathematicalrequirement this implies that we should have data starting at time 0.
If we were to use age as timescale for cumulative risk, we would want data available sincebirth; if we only had observations where most people entered between 20 and 40 years ofage, we could mathematically compute cumulative risk to some age, but it would be
nonsense. Instead we would compute the cumulative risk

given that a person attained age

40, say. In that case the time scale would be age − 40.

1.2 Rates and cumulative risks
Each of the cumulative risks is a function of the survival function which in turn depends on

all rates. Specically, if the cause-specic rates are λc (t), c = 1, 2, . . ., then the survival
function (probability of being alive at time t) is:

!
!
Z tX
X
λc (s) ds = exp −Λc (t)
S(t) = exp −0
The quantities Λc (t) =

c

(1.1)

c

rates

Rt

λ (s) ds are called cumulative(probabilists call them
0 cintegrated intensities), although they are not rates. Cumulative rates are dimensionless,but they have no probability interpretation of any kind.
The cumulativerisk, the probability of dying from cause k, say, before time t, Rk (t) is:

!
Z t
Z t
X
Λj (u) duRk (t) = λk (u)S(u) du = λk (u) exp −
0

0

(1.2)

j

The rationale is that λk (u) du is the probability of death from cause k in the small interval(u, u + du), conditional on being alive at time u. If this is multiplied with the probabilityof being alive at u, S(u), Bayes rule tells us that we get the marginal probability of deathfrom cause k in the small interval (u, u + du), This function of u is the argument in the
integral; so integration from u = 0 to u = t, will give the probability of death from cause canywhere in (0, t)the cumulative risk of cause k at t.
Parametric models for the cause-specic rates can produce estimated transition rates λc
at closely spaced intervals, and the cumulative risks can then be estimated from these bysimple numerical integration; this is illustrated in the next chapter.

Note that at any one time every person is either alive or dead from one of the causes, so
the sum of the survival and the cumulative risks is always 1:

1 = S(t) + R1 (t) + R2 (t) + · · · ,∀tCompeting risks

1.2

Rates and cumulative risks 31.2.1 Condence intervals by simulation
But even if we from the modeling of the λc s may have standard errors of log


λc (t) , thestandard errors of the Rc s will be analytically intractable from these.
In practice, the only viable way to get condence intervals for the cumulative risks, Rc , istherefore by calculation of a set of rates λc (t) by sampling from the posterior distributionof the parameters in the models for log


λc (t) , and then compute the integrals numericallyfor each simulated sample, according to formulae 1.1 and 1.2. This will produce a so-calledparametric bootstrap sample of the cumulative risks from which we can derive condence
intervals
The simulation approach also allows calculation of condence intervals for sums of the
cumulative risks, R1 (t) + R2 (t), for example, which will be needed if we want to show
stacked cumulative risks.
Finally, it will also allow calculation of standard errors of sojourn times in each of thestates alive and cause1, cause2,. . . . While the latter two may not be of direct interest,thendierences between such sojourn times between dierent groups can be interpreted as

years of life lost to each cause between groups.

1.2.2 Subdistribution hazard
A common concept seen in competing risks is the subdistribution hazard, and proportional
hazards models for this (the Fine-Gray model).
Suppose we only consider all-cause mortality, λ(t). Then the cumulative risk of death is: Z t
R(t) = 1 − S(t) = 1 − exp − λ(s) ds
0Solving this for λ will yield:

d log 1 − R(t)λ(t) = −dt
In the case of multiple causes of death we can dene thesubdistribution hazard for cause cby using the same transformation of the cumulative risk for cause c:d log 1 − Rc (t)λ̃c (t) = −dtor, for the cause-specic risk: Z tRc (t) = 1 − exp − λ̃c (s) ds
0Thus, the subdistribution hazard λ̃c (t) is a function which, when subjected to the(hazard→risk) function from the all-cause mortality case, yields the cause-specic risk.
The Fine-Gray model is a model for the subdistribution hazard λ̃c . It is only a model for

one cause-specic hazard. Of course it can be applied to all available causes in turn, but
the sum of the cumulative risks derived from the models may exceed 1. . .
Unlike a cause-specic hazard, which can depend on multiple time scales, thesubdistribution hazard can only depend on one, since it requires an originjust likecumulative risks.

But the interpretation of a subdistribution hazard is dicult. I have yet to see one that
goes beyond the mathematical formalism above. Therefore the subdistribution is notincluded in this vignette.

Chapter 2
Example data
2.1 A Lexis object
As an illustrative data example we use the (fake) diabetes register data; we set up the
Lexis object, an then cut the follow-up time at dates of OAD and Ins:

> library(Epi)> data(DMlate)> Ldm <- Lexis(entry = list(per = dodm,+
age = dodm-dobth,+
tfd = 0),+exit = list(per = dox),+
exit.status = factor(!is.na(dodth), labels = c("DM", "Dead")),+data = DMlate[sample(1:nrow(DMlate), 1000),])
NOTE: entry.status has been set to "DM" for all.
> summary(Ldm, t = T)
Transitions:To
From DM Dead Records: Events: Risk time: Persons:
DM 749 2511000
2515221.16
1000
Timescales:
per age tfd"" "" ""> set.seed(1952)> Mdm <- mcutLexis(Ldm,+wh = c('dooad','doins'),+
new.states = c('OAD','Ins'),+
seq.states = FALSE,+
ties = TRUE)
NOTE: Precursor states set to DM
NOTE: 3 records with tied events times resolved (adding 0.01 random uniform),so results are only reproducible if the random number seed was set.
> summary(Mdm)
Transitions:To
FromDM Dead OAD Ins Ins+OAD Records: Events: Risk time: Persons:DM307
94 290 74
0765
4582281.28765
OAD
0 103 311 0
95
509
198
2154.83
509
Ins
0 16
0 52
22
90
38
414.91
90

4

Example data

2.1

Ins+OAD 0 38 00Sum
307 251 601 126

79
196

117
1481

38
732

A Lexis object 5

370.145221.16117
1000

We initially split the FU before drug inception in intervals of 1/12 year, creating a Lexis
object for a competing risks situation with three possible event types:> Sdm <- splitLexis(factorize(subset(Mdm,+
lex.Cst == "DM")),+
time.scale = "tfd",+
breaks = seq(0, 20, 1/12))> summary(Sdm)
Transitions:To
FromDM Dead OAD Ins Records: Events: Risk time: Persons:
DM 27318 94 290 74
27776
458
2281.28765
We can illustrate the follow-up in the full data frame and in the restricted data frame> boxes(Mdm, boxpos = list(x = c(15, 50, 15, 85, 85),+
y = c(85, 50, 15, 85, 15)),+
scale.R = 100,+
show.BE = TRUE)

74
(3.2)

DM
2,281.3
765
307

Ins
414.9
16

52

94
(4.1)
16
(3.9)

290
(12.7)
0

Dead
251

22
(5.3)

103
(4.8)
OAD
2,154.8
219
311

95
(4.4)

38
(10.3)
Ins+OAD370.1
079

Figure 2.1: The transitions in the multistate model, where follow-up is extended also after
beginning of rst drug exposure. Rates in brackets are per 100 PY.
./crisk-boxes5> boxes(Relevel(Sdm, c(1, 4, 2, 3)),+
boxpos = list(x = c(15, 85, 75, 15),+y = c(85, 85, 30, 15)),+
scale.R = 100,+
show.BE = TRUE )

6

2.2

Models for rates

DM
2,281.3
765
307

290
(12.7)

CmpRskParSim

74
(3.2)

Ins
0

94
(4.1)

0

0

74

Dead
94

OAD
290

Figure 2.2: The transitions in the competing risks model, where follow-up is stopped at drug
exposure. By that token only the DM state has person-years; a characteristic of a competingrisks situation.
./crisk-boxes4

2.2 Models for rates
Now that we have set up a dataset with three competing events, we can model thecause-specic rates separately by time from diagnosis as the only underlying time scale.
This is done by Poisson-regression on the time-split data set; since the dataset is in

Lexis format we can use the convenience wrapper gam.Lexis to model rates a smoothfunctions of time (tfd). Note that we only need to specify the to= argument because there
is only one possible from for each to (incidentally the same for all to states, namely DM):
> mD <- gam.Lexis(Sdm, ~ s(tfd, k = 5), to = 'Dead')
mgcv::gam Poisson analysis of Lexis object Sdm with log link:
Rates for the transition:
DM->Dead> mO <- gam.Lexis(Sdm, ~ s(tfd, k = 5), to = 'OAD' )mgcv::gam Poisson analysis of Lexis object Sdm with log link:
Rates for the transition:
DM->OAD> mI <- gam.Lexis(Sdm, ~ s(tfd, k = 5), to = 'Ins' )mgcv::gam Poisson analysis of Lexis object Sdm with log link:
Rates for the transition:
DM->Ins
We see that gam.Lexis (just like glm.Lexis would) tells us what transitions are modeled.
With these models tted we can compute the rates, cumulative rates and the cumulative
risks and sojourn times in states using the usual formulae. First we compute the rates inExample data

2.3

Cumulative rates and risks 7intervals of length 1/100 years. Note that these models only have time since diagnosis ascovariates, so they are the counterpart of Nelson-Aalen estimates, albeit in a biologicallymore meaningful guise.
The points where we compute the predicted rates are midpoints of intervals of length
1/100 year. These points are unrelated to the follow-up intervals in which we split the data
for analysisthey were 1 month intervals, here we use 1/10 year (about 1 month, in reallife a smaller interval should be used, say 1/50 year):> nd <- data.frame(tfd = seq(0, 10, 1/10))> rownames(nd) <- nd$tfd> str(nd)'data.frame':
101 obs. of 1 variable:$ tfd: num 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
With this we can show the rates as a function of the time since entry (diagnosis of
diabetes):> matshade(nd$tfd, cbind(ci.pred(mD, nd),+ci.pred(mI, nd),+ci.pred(mO, nd)) * 1000,+
col = c("black", "red", "blue"), log = "y", lwd = 3, plot = TRUE,+
xlab = "Time since DM diagnosis (years)",+ylab = "Rates per 1000 PY", ylim = c(0.05,500), yaxt = "n")> axis(side = 2, at = ll<-outer(c(1,2,5),-2:3,function(x,y) x*10^y),+labels = formatC(ll,digits = 4), las = 1)
> axis(side = 2, at = ll<-outer(c(1.5,2:9),-2:3,function(x,y) x*10^y),+labels = NA, tcl = -0.3)> text(0, 0.5*0.6^c(1,2,0),+c("Dead","Ins","OAD"),+
col = c("black","red","blue"), adj = 0)Note that the graph in gure 2.3 is not normally shown in analyses of competing risks;
the competing cause-specicrates are hardly ever shown. I suspect that this is frequently

because they are often modeled by a Cox model and so are buried in the model and hardto get at.
Since we will be integrating the rates, it would be relevant to show the rates on a linear
scale instead, de-emphasizing the very small uctuations of the Ins rates:> matshade(nd$tfd, cbind(ci.pred(mD, nd),+ci.pred(mI, nd),+ci.pred(mO, nd))*1000,+
col = c("black", "red", "blue"), lwd = 3, plot = TRUE,+
xlab = "Time since DM diagnosis (years)",+ylab = "Rates per 1000 PY", ylim = c(0,500), yaxs = "i")> text(8, 500 - c(2, 3, 1) * 20,+c("Dead","Ins","OAD"),+
col = c("black","red","blue"), adj = 0)

2.3 Cumulative rates and risks
For the calculation of the cumulative rates and state probabilities, we need just theestimated rates (without CIs). The formulae 1.1 and 1.2 on page 2 are transformed to

R-code; starting with the rates, λD as lD etc:

8

2.3

Cumulative rates and risks

CmpRskParSim

500
200
100

Rates per 1000 PY

50
20
10
5
2
1
0.5

OAD
Dead

0.2

Ins

0.1
0.05
0

2

4

6

8

10

Time since DM diagnosis (years)

Figure 2.3: Estimated rates from the DM state, estimates are from gam models tted to data
split in 1 month intervals (1/12 year, that is). Rates of OAD is in the vicinity of 0.1/year,
and mortality about half of this. Rates of insulin start among persons on no other drug are
beginning high, then decreasing with a nadir at about 4 years and then increase to a peak at
8 years.
./crisk-rates> # utility function that calculates the midpoints between sucessive> # values in a vector> mp <- function(x) x[-1] - diff(x) / 2> #> int <- 1/10
> # rates at midpoints of intervals> lD <- mp(ci.pred(mD, nd)[,1])> lI <- mp(ci.pred(mI, nd)[,1])> lO <- mp(ci.pred(mO, nd)[,1])> #> # cumulative rates and survival function at right border of the intervals> LD <- cumsum(lD) * int> LI <- cumsum(lI) * int> LO <- cumsum(lO) * int> # survival function, formula (1.1)> Sv <- exp(- LD - LI - LO )> #> # when integrating to get the cumulative *risks* we use the average> # of the survival function at the two endpoints> # (adding 1 as the first), formula (1.2)> Sv <- c(1, Sv)> rD <- c(0, cumsum(lD * mp(Sv)) * int)
> rI <- c(0, cumsum(lI * mp(Sv)) * int)
> rO <- c(0, cumsum(lO * mp(Sv)) * int)

Example data

2.3

Cumulative rates and risks 9

500
OAD
Dead
Ins

Rates per 1000 PY

400

300

200

100

0
0

2

4

6

8

10

Time since DM diagnosis (years)

Figure 2.4: Estimated rates from the DM state, estimates are from gam models tted to data
split in 1 month intervals (1/12 year, that is). Rates of OAD is in the vicinity of 0.1/year,
and mortality about half of this. .
./crisk-rates-l

Now we have the cumulative risks for the three causes and the survival, computed at theend of each of the intervals. At any time point the sum of the 3 cumulative risks and thesurvival should be 1:> summary(rD + rI + rO + Sv)
Min. 1st Qu. Median
Mean 3rd Qu.
1
1
1
1
1
> oo <- options(digits = 20)> cbind(summary(Sv + rD + rI + rO))
[,1]
Min.
1.0000000000000000000
1st Qu. 1.0000900529180785092
Median 1.0000906409036778300
Mean
1.0000881929293825578
3rd Qu. 1.0000915436942121595
Max.
1.0000920704233506875> options(oo)Max.
1. . . and bar a rounding error, they are.
We can then plot the 3 cumulative risk functions stacked together using mat2pol (matrix

to polygons):> zz <- mat2pol(cbind(rD, rI, rO, Sv), x = nd$tfd,+
xlim = c(0,10), xaxs = "i", yaxs = "i", las = 1,+
xlab = "Time since DM diagnosis (years)",+ylab = "Probability",

10

2.4

Sojourn times

CmpRskParSim+
col = c("black","red","blue","forestgreen"))> text(9, mp(zz["9", ]), c("Dead", "Ins", "OAD"," DM"), col = "white")> box(col = "white", lwd = 3)

1.0

DM
0.8

Probability

0.6
OAD
0.4

Ins

0.2

Dead
0.0
0

2

4
6
Time since DM diagnosis (years)

810

Figure 2.5: Probabilities of being in the 4 dierent states as a function of time since diagnosis. Note that OAD means that OAD was initiated rst, and similarly for Ins. We are notconcerned about what occur after these events. Dead means dead without being on any drug.

./crisk-stack

2.4 Sojourn times
The sojourn times in each of the states is just the area of each of the coloured parts ofgure 2.5. Since the y -dimension of the plot is probability (dimensionless) and the x-axis
has dimension time, the computed areas will have dimension time.
Normally we will not report the sojourn times as functions of (truncation) time, but only
report them at a few select truncation points, such as 5 or 10 years. Calculation of the 10
year sojourn times would be straight-forward as integrals from 0 to 10these calculations
rely on the predicted rates from nd being for the rst 10 years:> Sj <- c(sjA = sum(Sv * int),+
sjD = sum(rD * int),+
sjI = sum(rI * int),+
sjO = sum(rO * int))> c(Sj, sum(Sj))

Example data
sjA
sjD
sjI
4.5091300 1.1125393 0.9006293

2.4

Sojourn times 11

sjO
3.5785921 10.1008907

We see that there is a some rounding error in the calculations; the sum should really be
exactly 10.
This was a demonstration on how to compute the rates, cumulative risks and sojourntimes. But no condence intervals.

Chapter 3Condence intervals
Besides condence intervals for each of the 4 cumulative risks, we will also be interested incondence intervals forsums of any subset of the cumulative risks, corresponding to theborders between the colours in gure 2.5. If we only had two competing risks (and hence
three states) the latter would not be an issue, because the sum of any two cumulative risks
will be 1 minus the cumulative risk of the remainder, so we could get away with thecondence intervals for the single cumulative risks. This is the reason we have chosen an
example with 3 competing risks and not just 2; we then have 4 probabilities to sum indierent order.
A short look at the formulae for cumulative risks will reveal that analytic approximation
to the standard error of these probabilities (or some transform of them) is not really a
viable way to go. Particularly if we also want condence intervals of sums of the state
probabilities as those shown in stacked plots.
So in practice, if we want condence intervals not only for the state probabilities, butalso for any sum of subsets of them we would want a large number of simulated copies of
the cumulative risks, each copy being of the same structure as the one we just extracted
from the models.
Condence intervals for sojourn times (i.e. time spent) in each state up to a given time,
would come almost for free from the simulation approach.
This means that we must devise a method to make a prediction not from the estimatedmodel, but where we instead of the model parameters use a sample from the posterior
distribution of the estimated parameters. Here, the posterior distribution of the parameters
will be taken to be the multivariate normal distribution with mean equal to the vector ofparameter estimates and variance-covariance matrix equal to the estimatedvariance-covariance matrix of the parameters.
Precisely this approach is implemented in ci.lin via the sample argument; we can get a

predicted value from a given prediction data frame just as from ci.pred resp. ci.exp;
here is shown two dierent ways of getting predicted values of the cause-specic rates:> head(cbind(ci.pred(mI, nd),+ci.exp (mI, nd)))
Estimate
2.5%
97.5% exp(Est.)
2.5%
97.5%
0 0.33135605 0.24687209 0.44475191 0.33135605 0.24687209 0.44475191
0.1 0.24070162 0.18632747 0.31094326 0.24070162 0.18632747 0.31094326
0.2 0.17487157 0.13821245 0.22125406 0.17487157 0.13821245 0.22125406
0.3 0.12709035 0.10031831 0.16100706 0.12709035 0.10031831 0.16100706
0.4 0.09242716 0.07130434 0.11980729 0.09242716 0.07130434 0.11980729

12

Condence intervals

3.1

Common parameters across cause-specic rates 13

0.5 0.06729135 0.04990789 0.09072967 0.06729135 0.04990789 0.09072967
And here is an illustration of the prediction with model based condence intervals for the
rates, alongside predictions based on samples from the posterior distribution of theparameters in the model:> str(ci.lin(mI, nd, sample = 4))
num [1:101, 1:4] -1.22 -1.51 -1.79 -2.08 -2.36 ...
- attr(*, "dimnames")=List of 2..$ : chr [1:101] "0" "0.1" "0.2" "0.3" ...
..$ : NULL
> head(cbind(ci.pred(mI,nd), exp(ci.lin(mI, nd, sample = 4))))
Estimate
2.5%
97.5%
0 0.33135605 0.24687209 0.44475191 0.32655476 0.31871501 0.31257509 0.26618487
0.1 0.24070162 0.18632747 0.31094326 0.24019812 0.24053668 0.22792707 0.20306301
0.2 0.17487157 0.13821245 0.22125406 0.17669693 0.18154911 0.16622527 0.15492093
0.3 0.12709035 0.10031831 0.16100706 0.13002103 0.13705671 0.12127214 0.11821613
0.4 0.09242716 0.07130434 0.11980729 0.09572787 0.10351097 0.08853954 0.09024259
0.5 0.06729135 0.04990789 0.09072967 0.07054236 0.07822822 0.06471662 0.06893161
Note that we use exp(ci.lin(...this is because the sample= argument does not workwith ci.exp.

The simulation (parametric bootstrapping) is taking place at the parameter level and thetransformation to survival and cumulative risks is simply by a function applied to eachsimulated set of rates.

3.1 Common parameters across cause-specic rates
Note that we have implicitly been assuming that the transitions are being modeled
separately. If some transitions are modeled jointlyfor example assuming that the rates of

OAD and Ins are proportional as functions of time since entry, we are in trouble, because wethen need one sample from the posterior generating two dierent predictions, one for eachof the transitions modeled together. Moreover the model will have to be a model tted to astack.Lexis object, so a little more complicated to work with.
A simple way to program this would be to reset the seed to the same value before
simulating with dierent values of nd, this is what is intended to be implemented, but is
not yet. This is mainly the complication of having dierent prediction frames for dierent
risks in this case.
However, this is not a very urgent need, since the situation where you want common
parameters for dierent rates out of a common state is quite rare. By that token thefacility is not likely to be implemented anytime soon, if ever.

3.2 Simulation based condence intervals
The parametric bootstrap is implemented in the function ci.Crisk (condence intervals
for Cumulative risks) in the Epi package:

We can now run the function using the model objects for the three competing events,using a common prediction data frame, nd for the rates. The time points in the frame must

14

3.2

Simulation based condence intervals

CmpRskParSim

be so closely spaced that it makes sense to assume the rates constant in each interval; here
we use intervals of length 1/10 years, approximately 1 monthtoo large for realapplications:> system.time(+ res <- ci.Crisk(list(OAD = mO,+
Ins = mI,+
Dead = mD),+
nd = data.frame(tfd = seq(0, 10, 1/10)),+nB = 100,+perm = 4:1))
NOTE: Times are assumed to be in the column tfd at equal distances of 0.1
bruger system forløbet
0.37
0.02
0.39
> str(res)
List of 4$ Crisk: num [1:101, 1:4, 1:3] 1 0.927 0.869 0.822 0.784 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfd : chr [1:101] "0" "0.1" "0.2" "0.3" ...
.. ..$ cause: chr [1:4] "Surv" "OAD" "Ins" "Dead".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ Srisk: num [1:101, 1:3, 1:3] 0 0.00501 0.00964 0.01398 0.0181 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfd : chr [1:101] "0" "0.1" "0.2" "0.3" ...
.. ..$ cause: chr [1:3] "Dead" "Dead+Ins" "Dead+Ins+OAD".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ Stime: num [1:101, 1:4, 1:3] 0 0.0963 0.1861 0.2707 0.3509 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfd : chr [1:101] "0" "0.1" "0.2" "0.3" ...
.. ..$ cause: chr [1:4] "Surv" "OAD" "Ins" "Dead".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ time : num [1:101] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
- attr(*, "int")= num 0.1
As we see, the returned object (res) is a list of length 4, the rst 3 components are 3-way
arrays, and the last the vector of times of the rst dimension of the arrays. The latter ismainly for convenience in further processingit is easier to write res$as.numeric(dimnames(res$ Crisk)[[1]]).
The three rst components of res represent:

time than Crisk Cumulative risks for each state Srisk Stacked cumulative risks across states
 Stime Sojourn times in each state, truncated at each point of the time dimension.
The rst dimension of each array is time corresponding to endpoints of intervals oflength int, (normally assumed starting at 0, but not necessarily). The second dimension isstates (or combinations thereof ). The last dimension of the arrays is the type of statistic;

50% is the median of the samples, and the bootstrap condence intervals as indicated;taken from the alpha argument to ci.Crisk (defaults to 0.05).
The argument perm governs in which order the state probabilities are stacked in the
Srisk element of the returned list, the default is the states in the order given in the list ofmodels in the rst argument to ci.Crisk followed by the survival.

Condence intervals

3.3

Simulated condence intervals for rates 15

If we want the bootstrap samples to make other calculations we can ask the function to
return the bootstrap samples of the rates by using the argument sim.res(defaults to 'none'):= 'rates'> system.time(+ rsm <- ci.Crisk(list(OAD = mO,+
Ins = mI,+
Dead = mD),+
nd = data.frame(tfd = seq(0, 10, 1/10)),+nB = 100,+
sim.res = 'rates'))
NOTE: Times are assumed to be in the column tfd at equal distances of 0.1
bruger system forløbet000> str(rsm)num [1:101, 1:3, 1:100] 0.473 0.43 0.391 0.356 0.324 ...
- attr(*, "dimnames")=List of 3..$ time: chr [1:101] "0" "0.1" "0.2" "0.3" ...
..$ mod : chr [1:3] "OAD" "Ins" "Dead"..$ sim : chr [1:100] "1" "2" "3" "4" ...
- attr(*, "int")= num 0.1- attr(*, "time")= num [1:101] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
This is 100 bootstrap samples of the rates evaluated at the 501 endpoints of the intervals.
Alternatively we can get the bootstrap samples of the cumulative risks by setting

sim.res = 'crisk':> system.time(+ csm <- ci.Crisk(list(OAD = mO,+
Ins = mI,+
Dead = mD),+
nd = data.frame(tfd = seq(0, 10, 1/10)),+nB = 100,+
sim.res = 'crisk'))
NOTE: Times are assumed to be in the column tfd at equal distances of 0.1
bruger system forløbet
0.11
0.02
0.13
> str(csm)num [1:101, 1:4, 1:100] 1 0.931 0.876 0.83 0.792 ...
- attr(*, "dimnames")=List of 3..$ tfd : chr [1:101] "0" "0.1" "0.2" "0.3" ...
..$ cause: chr [1:4] "Surv" "OAD" "Ins" "Dead"..$ sim : chr [1:100] "1" "2" "3" "4" ...
- attr(*, "int")= num 0.1
These are 100 simulated samples of the cumulative risks evaluated at the 501 endpoints ofnd dimension

the intervals, and also includes the survival probability in the rst slot of the 2of csm.

3.3 Simulated condence intervals for rates
In gure 2.3 we showed the rates with condence intervals from the model. But in rsm wehave 100 parametric bootstrap samples of the occurrence rates, so we can derive thebootstrap medians and the bootstrap c.i.s:

16

3.3

Simulated condence intervals for rates

CmpRskParSim> Brates <- aperm(apply(rsm,+
1:2,+
quantile,+
probs = c(.5, .025, .975)),+c(2, 3, 1))> str(Brates)
num [1:101, 1:3, 1:3] 0.452 0.412 0.374 0.341 0.31 ...
- attr(*, "dimnames")=List of 3..$ time: chr [1:101] "0" "0.1" "0.2" "0.3" ...
..$ mod : chr [1:3] "OAD" "Ins" "Dead"..$: chr [1:3] "50%" "2.5%" "97.5%"(aperm permutes the dimensions of the array). Then we can plot the bootstrap estimateson top of the estimates based on the normal approximation to distribution of theparameters. They arenot surprisinglyin close agreement since they are both based on
an assumption of normality of the parameters on the log-rate scale:

500
200
100

Rates per 1000 PY

50
20
10
5
2
1
0.5

OAD
Dead

0.2

Ins

0.1
0

2

4
6
Time since DM diagnosis (years)

8

10

Estimated rates from the DM state, estimates are from gam models tted to data
split in 1 month intervals (1/12 year, that is). The white dotted curves are the bootstrap
medians, black dotted curves are the bootstrap 95% c.i.s.
./crisk-rates-ci
Figure 3.1:

> matshade(nd$tfd, cbind(ci.pred(mD, nd),+ci.pred(mI, nd),+ci.pred(mO, nd)) * 1000,+
ylim = c(0.1,500), yaxt = "n",+ylab = "Rates per 1000 PY",+
xlab = "Time since DM diagnosis (years)",

Condence intervals

3.4

Condence intervals for cumulative risks 17+
col = c("black","red","blue"), log = "y", lwd = 3, plot = TRUE)> matlines(nd$tfd,+
cbind(Brates[,"Dead",],+Brates[,"Ins" ,],+Brates[,"OAD" ,]) * 1000,+
col = c("white", "black", "black"), lty = 3, lwd = c(3,1,1))> axis(side = 2, at = ll<-outer(c(1,2,5),-2:3,function(x,y) x*10^y),+labels = formatC(ll,digits = 4), las = 1)
> axis(side = 2, at = ll<-outer(c(1.5,2:9),-2:3,function(x,y) x*10^y),+labels = NA, tcl = -0.3)> text(0, 0.5*0.6^c(1,2,0),+c("Dead", "Ins", "OAD"),+
col = c("black", "red", "blue"), adj = 0)

3.4 Condence intervals for cumulative risks
In the Crisk component of res we have the cumulative risks as functions of of time, withbootstrap condence intervals, so we can easily plot the three cumulative risks:

0.5

0.4

Cumulative probability

OAD
Dead
Ins

0.3

0.2

0.1

0.0
0

2

4
6
Time since DM diagnosis (years)

8

10

Cumulative risks for the three types of events, with 95% bootstrap-based condence
intervals as shades.
./crisk-crates

Figure 3.2:

> matshade(res$time,+
cbind(res$Crisk[,"Dead",],+res$Crisk[,"Ins" ,],

18

3.5

Condence intervals for stacked cumulative risks

CmpRskParSim+
res$Crisk[,"OAD" ,]), plot = TRUE,+
xlim = c(0,10), xaxs = "i", yaxs = "i", las = 1,+
xlab = "Time since DM diagnosis (years)",+ylab = "Cumulative probability",+col = c("black","red","blue"))> text(8, 0.3 + c(1, 0, 2) / 25,+c("Dead", "Ins", "OAD"),+
col = c("black", "red", "blue"), adj = 0)

3.5 Condence intervals for stacked cumulative risks
Unlike the single cumulative risks where we have a condence interval for each cumulativerisk, when we want to show the stacked probabilities we must deliver the condence
intervals for the relevant sums, they are in the Srisk component of res.

> str(res$Crisk)num [1:101, 1:4, 1:3] 1 0.927 0.869 0.822 0.784 ...
- attr(*, "dimnames")=List of 3..$ tfd : chr [1:101] "0" "0.1" "0.2" "0.3" ...
..$ cause: chr [1:4] "Surv" "OAD" "Ins" "Dead"..$: chr [1:3] "50%" "2.5%" "97.5%"> str(res$Srisk)num [1:101, 1:3, 1:3] 0 0.00501 0.00964 0.01398 0.0181 ...
- attr(*, "dimnames")=List of 3..$ tfd : chr [1:101] "0" "0.1" "0.2" "0.3" ...
..$ cause: chr [1:3] "Dead" "Dead+Ins" "Dead+Ins+OAD"..$: chr [1:3] "50%" "2.5%" "97.5%"
But we start out by plotting the stacked probabilities using mat2pol (matrix to polygon),the input required is the single components from the Crisk component. Then we add thecondence intervals as white shades (using matshade):> zz <- mat2pol(res$Crisk[,c("Dead", "Ins", "OAD", "Surv"),1],+x = res$time,+
xlim = c(0,10), xaxs = "i", yaxs = "i", las = 1,+
xlab = "Time since DM diagnosis (years)",+ylab = "Probability",+col = c("black","red","blue","forestgreen") )> text(9, mp(zz["9",]), c("Dead", "Ins", "OAD", "DM"), col = "white" )> matshade(res$time,+
cbind(res$Srisk[,1,],+
res$Srisk[,2,],+
res$Srisk[,3,]),+
col = 'transparent', col.shade = "white", alpha = 0.4)

3.6 Sojourn times
From the Stime component of the res we can derive the estimated time spent in each state
during the rst, say, 5 or 10 years:
When referring to the times, we use

character values5 and 10 years are not necessarily

th and 10th positions of the rst dimension of the Stime array:at the 5

Condence intervals

3.6

Sojourn times 19

1.0

DM
0.8

Probability

0.6
OAD
0.4

Ins

0.2

Dead
0.0
0

2

4
6
Time since DM diagnosis (years)

810

Figure 3.3: Probabilities of being in the 4 dierent states as a function of time since diagnosis. Note that OAD means that OAD was initiated rst, and similarly for Ins. We are notconcerned about what occurs after these events. Dead means dead without being on any drug.
The white shadings around the borders between coloured areas represent the 95% condence
intervals for the (sum of) probabilities.
./crisk-stack-ci> s510 <- res$Stime[c("5", "10"),,]> dimnames(s510)[[1]] <- c(" 5 yr","10 yr")
> round(ftable(s510, row.vars=1:2), 2)
50% 2.5% 97.5%
tfd cause
5 yr Surv 2.80 2.63 2.92
OAD
1.42 1.30 1.57
Ins
0.43 0.34 0.52
Dead 0.36 0.28 0.45
10 yr Surv 4.42 4.09 4.69
OAD
3.54 3.27 3.90
Ins0.93 0.73 1.14
Dead 1.12 0.91 1.36
So we see that the expected life lived without pharmaceutical treatment during the rst 10
years after DM diagnosis is 4.31 years with a 95% CI of (4.21;4.42), and during the rst 5
years 2.77 (2.72;2.82).

Chapter 4
A simple illustration
The following is a terse cook-book illustration of how to use the ci.Crisk function.

4.1 Data
For illustration we simulate some causes of death in the DMlate data set; rst we sample

numbers 1, 2, 3 representing 3 dierent causes of death in DMlate:> data(DMlate)> set.seed(7465)> wh <- sample(1:3, nrow(DMlate), r=T, prob = c(4, 2, 6))Those not dead are changed to 0:> wh[is.na(DMlate$dodth)] <- 0Dene a factor in DMlate dening exit status as either alive or one of the three causes ofdeath, and check by a table that all dead have a cause:> DMlate$codth <- factor(wh, labels = c("Alive", "CVD", "Can", "Oth"))> with(DMlate, table(codth, isDead = !is.na(dodth)))isDeadcodth FALSE TRUE
Alive 7497
0
CVD
0 815
Can
0 401
Oth
0 1287

DMlate now looks like a typical data set with cause of death in a separate variable; in thiscase we also added a state, Alive, for those without a recorded death:> str(DMlate)'data.frame':
10000 obs. of 8 variables:$ sex : Factor w/ 2 levels "M","F": 2 1 2 2 1 2 1 1 2 1 ...
$ dobth: num 1940 1939 1918 1965 1933 ...
$ dodm : num 1999 2003 2005 2009 2009 ...
$ dodth: num NA NA NA NA NA ...
$ dooad: num NA 2007 NA NA NA ...
$ doins: num NA NA NA NA NA NA NA NA NA NA ...
$ dox : num 2010 2010 2010 2010 2010 ...
$ codth: Factor w/ 4 levels "Alive","CVD",..: 1 1 1 1 1 4 1 1 4 1 ...

20

A simple illustration

4.2

A Lexis object with 3 causes of death 21> head(DMlate, 12)
sex
dobth
dodm
dodth
dooad doins
dox codth
50185
F 1940.256 1998.917
NANA
NA 2009.997 Alive
307563 M 1939.218 2003.309
NA 2007.446
NA 2009.997 Alive
294104 F 1918.301 2004.552
NANA
NA 2009.997 Alive336439 F 1965.225 2009.261
NANA
NA 2009.997 Alive
245651 M 1932.877 2008.653
NANA
NA 2009.997 Alive
216824 F 1927.870 2007.886 2009.923
NA
NA 2009.923 Oth
24969
M 1946.498 2007.216
NA 2007.248
NA 2009.997 Alive
325465 M 1940.475 2006.674
NA 2006.674
NA 2009.997 Alive
430048 F 1937.083 1998.956 2008.464
NA
NA 2008.464 Oth362980 M 1933.154 2009.784
NANA
NA 2009.997 Alive
279976 F 1928.300 2000.949
NANA
NA 2009.997 Alive
279271 M 1951.213 2008.114
NA 2009.743
NA 2009.997 Alive

4.2 A Lexis object with 3 causes of deathWith cause of death in a separate variable it is easy to set up a Lexis object:> dmL <- Lexis(entry = list(per = dodm,+
age = dodm - dobth,+
tfD = 0),+exit = list(per = dox),+
exit.status = codth,+data = DMlate)
NOTE: entry.status has been set to "Alive" for all.
NOTE: Dropping 4 rows with duration of follow up < tol> summary(dmL, t = T)
Transitions:To
From
Alive CVD Can Oth Records: Events: Risk time:
Alive 7497 814 400 1285
9996
2499 54273.27

Persons:
9996

Timescales:
per age tfD"" "" ""
We can show the overall rates (the default boxes isvery primitive):> boxes(dmL, boxpos = TRUE)

4.3 Models for the rates
In order to model the cause-specic mortality rates by sex and time from diagnosis (tfD),
we rst split the data in 6-month intervals> sL <- splitLexis(dmL, time.scale="age", breaks = seq(0, 120, 1/2))> summary(sL)Transitions:To
From
Alive CVD Can Oth Records: Events: Risk time: Persons:
Alive 115974 814 400 1285
1184732499 54273.27
9996

22

4.4

Derived measures

CmpRskParSim

Alive
54,273.3
1,285
(0.0)

814
(0.0)

CVD

400
(0.0)

Oth

Can

Transitions from live to dierent causes of death. You probably want to explore
the other arguments to boxes.
./crisk-boxes
Figure 4.1:

> mCVD <- gam.Lexis(sL, ~ s(tfD, by=sex), to = "CVD")
mgcv::gam Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->CVD> mCan <- gam.Lexis(sL, ~ s(tfD, by=sex), to = "Can")
mgcv::gam Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->Can> mOth <- gam.Lexis(sL, ~ s(tfD, by=sex), to = "Oth")
mgcv::gam Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->Oth

4.4 Derived measures
With these three models for the occurrence rates we can compute the cumulative risks ofdying from each of the causes. We need a prediction data frame that gives the rates atclosely spaced times, in this case for men. For women the code would be practically thesame:> nm <- data.frame(tfD = seq(0, 15, 0.1), sex = "M")Note that we can rename the states as we please by naming the entries in the list of models
we supply to ci.Crisk:A simple illustration

4.4

Derived measures 23> cR <- ci.Crisk(list(CVD = mCVD,+
Can = mCan,+
Other = mOth),+nB = 100,+nd = nm)
NOTE: Times are assumed to be in the column tfD at equal distances of 0.1> str(cR)
List of 4$ Crisk: num [1:151, 1:4, 1:3] 1 0.993 0.987 0.981 0.975 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfD : chr [1:151] "0" "0.1" "0.2" "0.3" ...
.. ..$ cause: chr [1:4] "Surv" "CVD" "Can" "Other".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ Srisk: num [1:151, 1:3, 1:3] 0 0.00351 0.00684 0.00997 0.01297 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfD : chr [1:151] "0" "0.1" "0.2" "0.3" ...
.. ..$ cause: chr [1:3] "Other" "Other+Can" "Other+Can+CVD".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ Stime: num [1:151, 1:4, 1:3] 0 0.0997 0.1987 0.2971 0.3949 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfD : chr [1:151] "0" "0.1" "0.2" "0.3" ...
.. ..$ cause: chr [1:4] "Surv" "CVD" "Can" "Other".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ time : num [1:151] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
- attr(*, "int")= num 0.1
Note that we get three arrays: Crisk, the cumulative risks; Srisk, the stacked risks andStime, the sojourn times in each state. Finally, for convenience we also have thecomponent time, the times at which the cumulative risks are computed. It is also available
as the clumpy expression as.numeric(dimnames(cR$Crisk)[[1]]), but cR$time is easier.

4.4.1 Cumulative risks
We can plot the cumulative risks for death from each of the three causes, note we use thecolors from last. Note that the time points are in the time component of the Crisk object:> clr <- c("black", "orange", "limegreen")> matshade(cR$time, cbind(cR$Crisk[, "CVD" , ],+cR$Crisk[, "Can" , ],+cR$Crisk[, "Other", ]),+col = clr, lty = 1, lwd = 2,+
plot = TRUE, ylim = c(0, 1/3), yaxs = "i")> text(0, 1/3 - 1:3/30, c("CVD", "Can", "Oth"),+col = clr, adj = 0)
We also have the stacked probabilities so we can show how the population is distributed
across the states at any one time:

4.4.2 Stacked cumulative risks
We also get the stacked probabilities in the order that we supplied the models, so that if weplot these we get the probabilities of being dead from each cause as thedierence between

24

4.4

Derived measures

0.30

CmpRskParSim

CVD
Can

0.25
Oth

y

0.20

0.15

0.10

0.05

0.00
0

5

10

15

x
Figure 4.2:

rates.

Cumulative risks of each cause of death based on gam models for the cause-specic./crisk-cR

the curves. And the condence intervals are condence intervals for the cumulative sums ofprobabilities.

> matshade(cR$time, cbind(cR$Srisk[,1,],+cR$Srisk[,2,],+
cR$Srisk[,3,]),+
col = "black", lty = 1, lwd = 2,+
plot = TRUE, ylim = c(0,1), xaxs = "i", yaxs = "i")> text(14, mp(c(0, cR$Srisk["14", , 1], 1)),+rev(c(dimnames(cR$Crisk)[[2]])))> box(bty = "o")It is not a good idea to color the curves, they do not refer to the causes of death, it is theareas

between the curves that refer to causes.

It would be more logical to color the areas

between the curves. which can be done by

mat2pol (matrix to polygons) using the Crisk component. We can then superpose thecondence intervals for the sum of the state probabilities using matshade by adding white
shades:> zz <- mat2pol(cR$Crisk[, c("Other", "Can", "CVD", "Surv"), "50%"],+x = cR$time,+
xlim = c(0,15), xaxs = "i", yaxs = "i", las = 1,+
xlab = "Time since DM diagnosis (years)",+ylab = "Probability",+col = c("gray", "red", "blue", "limegreen") )A simple illustration

4.4

Derived measures 25

1.0

0.8
Surv

y

0.6

CVD

0.4

Can
0.2
Other

0.0
0

2

4

6

8

10

12

14

x
Figure 4.3:

Stacked cumulative risks.

./crisk-Sr1> matshade(cR$time, cbind(cR$Srisk[,1,],+cR$Srisk[,2,],+
cR$Srisk[,3,]),+
col = "transparent", col.shade = "white", alpha = 0.4)
> text(14, mp(c(0, cR$Srisk["14", , 1], 1)),+rev(c(dimnames(cR$Crisk)[[2]])), col = "white")

4.4.3 Sojourn times
The third component of the result, Stime is an array of sojourn times over intervals
starting at 0 and ending at the time indicated by the rst dimension:> ftable(round(cR$Stime[paste(1:5 * 3), , ], 1), row.vars = 1)
cause Surv
CVD
Can
Other
50% 2.5% 97.5% 50% 2.5% 97.5% 50% 2.5% 97.5%
50% 2.5% 97.5%
tfD
3
2.8 2.8 2.8 0.1 0.1 0.1 0.0 0.0 0.0 0.1 0.1 0.1
6
5.3 5.2 5.3 0.2 0.2 0.3 0.1 0.1 0.1 0.4 0.4 0.49
7.4 7.4 7.5 0.5 0.5 0.5 0.3 0.2 0.3 0.8 0.8 0.9
12
9.3 9.2 9.4 0.9 0.8 0.9 0.4 0.4 0.5 1.4 1.3
1.5
15
10.9 10.7 11.0 1.3 1.3 1.4 0.7 0.6
0.8 2.1 2.0 2.2
The sojourn times in the three dead states can be taken as the years of life lost to each of
the causes, the sum of the medians for the three causes equals the time frame (5, 10, 15)minus the Surv component.

26

4.4

Derived measures

CmpRskParSim

1.0

0.8
Surv

Probability

0.6

CVD

0.4

Can
0.2
Other

0.0
0

2

4

6
8
10
Time since DM diagnosis (years)

12

14

Figure 4.4: Stacked cumulative risks with coloring of states and overlaid with condence
intervals for the probabilities shown; that is the relevant sums.
./crisk-Sr2

So we see that during the rst 15 years after diagnosis of diabetes, the expected years
alive is 10.9 years. The distribution of lifetime lost between the causes is bogus in this caseas the causes of death were randomly generated.

4.4.4 Comparing groups
Finally, we may want to see thedierence (or ratio) of survival probabilities between men

and women, say. This can be derived from two bootstrap samples using dierent prediction
frames (the argument nd= to ci.Crisk). But the two bootstrap samples must come from thesame stream of samples from the multivariate normal. This can be obtained by explicitly
setting the seed for the random number generator to the same value before calling

ci.Crisk with each of the two dierent prediction frames as nd argument:> nm <- data.frame(tfD = seq(0, 15, 0.1), sex = "M")> nw <- data.frame(tfD = seq(0, 15, 0.1), sex = "F")> set.seed(1952)> mR <- ci.Crisk(list(CVD = mCVD,+
Can = mCan,+
Other = mOth),+nd = nm,+nB = 100,+
sim.res = "crisk" )NOTE: Times are assumed to be in the column tfD at equal distances of 0.1

A simple illustration

4.4

Derived measures 27> set.seed(1952)> wR <- ci.Crisk(list(CVD = mCVD,+
Can = mCan,+
Other = mOth),+nd = nw,+
nB = 100,+
sim.res = "crisk" )NOTE: Times are assumed to be in the column tfD at equal distances of 0.1> str(wR)num [1:151, 1:4, 1:100] 1 0.994 0.989 0.983 0.978 ...
- attr(*, "dimnames")=List of 3..$ tfD : chr [1:151] "0" "0.1" "0.2" "0.3" ...
..$ cause: chr [1:4] "Surv" "CVD" "Can" "Other"..$ sim : chr [1:100] "1" "2" "3" "4" ...
- attr(*, "int")= num 0.1
The two samples are now from identical streams of random numbers, so we can getdierences and ratios of the survival curves between men and women:> dS <- mR[,"Surv",] - wR[,"Surv",]> dS <- apply(dS, 1, quantile, probs = c(.5, .025, .975)) * 100> str(dS)num [1:3, 1:151] 0 0 0 -0.0681 -0.1771 ...
- attr(*, "dimnames")=List of 2..$
: chr [1:3] "50%" "2.5%" "97.5%"..$ tfD: chr [1:151] "0" "0.1" "0.2" "0.3" ...
> rS <- mR[,"Surv",] / wR[,"Surv",]> rS <- apply(rS, 1, quantile, probs = c(.5, .025, .975))
We can then plot the dierences and the ratios of the probabilitiesnote that thedimension of the function applyed becomes the rst dimension of the result:> par(mfrow = c(1,2))> matshade(as.numeric(colnames(dS)), t(dS), plot = TRUE,+
lwd = 3, ylim = c(-5, 5),+xlab = "Time since DM diagnosis (years)",+ylab = "Men - Women survival difference (%)")> abline(h = 0)
> matshade(as.numeric(colnames(rS)), t(rS), plot = TRUE,+
lwd = 3, ylim = c(1/1.2, 1.2), log ="y",+xlab = "Time since DM diagnosis (years)",+ylab = "Men - Women survival ratio")> abline(h = 1)
To illustrate the eect of not pairing the random samples we can generate a fresh
sample for women from a dierent stream and do the calculations to illustrate the excess
we get from not aligning samples.

> fR <- ci.Crisk(list(CVD = mCVD,+
Can = mCan,+
Other = mOth),+nd = nw,+
nB = 100,+
sim.res = "crisk" )NOTE: Times are assumed to be in the column tfD at equal distances of 0.1

28

4.4

Derived measures

CmpRskParSim

1.2

1.1
Men − Women survival ratio

Men − Women survival difference (%)

4

2

0

−2

1.0

0.9
−4

0

5
10
Time since DM diagnosis (years)

15

0

5
10
Time since DM diagnosis (years)

15

Figure 4.5: Dierences and ratios of survival between men and women, derived from thesame set of bootstrap samples from the parameter vector.
./crisk-difrat> dxS <- mR[,"Surv",] - fR[,"Surv",]> dxS <- apply(dxS, 1, quantile, probs = c(.5, .025, .975)) * 100> rxS <- mR[,"Surv",] / fR[,"Surv",]> rxS <- apply(rxS, 1, quantile, probs = c(.5, .025, .975))> par(mfrow = c(1,2))> matshade(as.numeric(colnames(dS)), t(dS), plot = TRUE,+
lwd = 3, ylim = c(-5, 5),+xlab = "Time since DM diagnosis (years)",+ylab = "Men - Women survival difference (%)")> matshade(as.numeric(colnames(dxS)), t(dxS), lty = 3, col = "forestgreen")> abline(h = 0)
> matshade(as.numeric(colnames(rS)), t(rS), plot = TRUE,+
lwd = 3, ylim = c(1/1.2, 1.2), log ="y",+xlab = "Time since DM diagnosis (years)",+ylab = "Men - Women survival ratio")> matshade(as.numeric(colnames(rxS)), t(rxS), lty = 3, col = "forestgreen")> abline(h = 1)
Start time: 2024-07-18, 15:25:10
End time: 2024-07-18, 15:25:53
Elapsed time: 0.72 minutes

A simple illustration

4.4

Derived measures 29

1.2

1.1
Men − Women survival ratio

Men − Women survival difference (%)

4

2

0

−2

1.0

0.9
−4

0

5
10
Time since DM diagnosis (years)

15

0

5
10
Time since DM diagnosis (years)

15Dierences and ratios of survival between men and women, derived from separate
bootstrap samples. The outer condence bands are from bootstrap samples not properly paired
between men end women.
./crisk-difratx

Figure 4.6: