Competing risks with ci.Crisk
Bendix Carstensen Steno Diabetes Center
Gentofte, Denmark
http://BendixCarstensen.com

SDCC - R study group, 16 April 2021

From /home/bendix/teach/AdvCoh/00/slides/crisk.tex

Friday 16th April, 2021, 09:22

1/ 30

Vignette in the Epi:
Competing risks with Lexis, parametric rates and simulation based
confidence intervals
or in http://BendixCarstensen.com/EpiDM
22,920.2
7,532
2,830

2,957
(12.9)

0

OAD
2,957

689(3.0)

Ins0689

1,056(4.6)

0Dead1,056
2/ 30

Competing risks. . . is where persons are subject to a number of different causes ofdeath, “cause1”, “cause2” etc.
ICauses of death are required to be
exhaustive and mutually exclusive:III

you will eventually die from one of the designated causes
you can only die from one.

The observed data will be a survival time and an exit status
which is either “censored alive” or one of the causes.

3/ 30

Competing risks
In situations where the causes are not causes of death but
other events, it is implicit that we only consider the firstoccurrence of an event from the state “alive”, and ignore whatoccurs after.
I Suppose we consider:
IOAD initiation
I Insulin initiation
I. . . we must also always include
II

Death (without OAD or Insulin)

to make sure possible events are exhaustive

4/ 30

Cause specific rates
λc (t) = P {event c in (t, t + h)| alive at t} /hAnalysis is in principle straight forward: estimate a model foreach of the cause-specific rates, λc , c = 1, 2, . . .
I modeling for λc :I

Event is occurrence of cause cAll other events are lumped together with the censoring
I (or: Just use the persons-years in the “alive” state)IIthe cause-specific rates will together form a complete modelfor the entire competing risks problem
I If the cause-specific rates are all we want, then we are done
I

5/ 30

DM
22,920.2
7,532
2,830

2,957
(12.9)

0

OAD
2,957

689
(3.0)

Ins
0689

1,056
(4.6)

0

Dead
1,056

6/ 30

Example data
> 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 )NOTE: entry.status has been set to "DM" for all.
NOTE: Dropping 4 rows with duration of follow up < tol> Mdm <- mcutLexis(Ldm,+
wh = c('dooad','doins'),+
new.states = c('OAD','Ins'),+
seq.states = FALSE,+
ties = TRUE)
NOTE: Precursor states set to DM
NOTE: 15 records with tied events times resolved (adding 0.01 random uniform),so results are only reproducible if the random number seed was set.
7/ 30

Example data> summary(Mdm)
Transitions:
ToFromDM Dead OAD Ins Ins+OAD
DM
2830 1056 2958 688
0OAD
0 992 3327
01006
Ins0 1520 462
171
Ins+OAD
0 299
00878
Sum
2830 2499 6285 1150
2055

Records:
7532
5325
785
1177
14819

Events: Risk time:4702
22920.19
199822965.26
323
3883.08299
4504.737322
54273.27

Persons:
7532
5325
785
1177
9996

8/ 30

Cause-specific rates
689
(3.0)

DM
22,920.2
7,532
2,830

Ins
3,883.1
97462

1,056(4.6)
152(3.9)

2,957(12.9)
0Dead2,499172(4.4)

992(4.3)
OAD
22,965.2
2,367
3,327

1,005(4.4)

299(6.6)
Ins+OAD
4,504.70878
9/ 30

Split data—only time in DM> Sdm <- splitMulti(factorize(subset(Mdm, lex.Cst == "DM")),+
tfd = seq(0, 20, 1/12))
NOTE: lex.Cst and lex.Xst now have levels:
DM Dead OAD Ins> summary(Sdm)
Transitions:To
FromDM Dead OAD Ins
DM 274263 1056 2958 688

Records:278965

Events: Risk time:4702
22920.19Persons:7532

10/ 30DM
22,920.2
7,532
2,830

2,957(12.9)

0

OAD
2,957

689(3.0)

Ins0689

1,056(4.6)

0

Dead
1,056

11/ 30Poisson analysis of time-split data

> 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. . . wrapper for:> mI <- gam(cbind(lex.Xst == "Ins" & lex.Cst != lex.Xst,+lex.dur)+~ s(tfd, k = 5),+family = poisreg,+
data = subset(Sdm, lex.Cst == "DM"))
12/ 30Cause-specific rates
500
200
10050Rates per 1000 PY

20
105210.50.2OAD
DeadIns0.10.050.02024
6Time since DM diagnosis (years)

810

13/ 30Cumulative risks
In addition to the rates we are interested in:survival probability, andI the cumulative risks — probability of dying from each cause
I i.e. the probabilities of being en each of the 4 boxes at any onetimeI

14/ 30Cumulative risks
ISurvival probability!
!
Z tX
X
S (t) = exp −λc (s) ds = exp −Λc (t)
0cc—the cumulative risk of survival. . .
I Cumulative risks!
Z tZ tXΛc (u) duRc (t) = λc (u)S (u) du = λc (u)exp −00c

These are probabilities, dimensionless.
I Note that integrals are just sums, so easy to computenumerically
I15/ 30

Computing the rates
We need λD (t), λI (t) and λO (t).
So first choose the ts, then compute rates at these times> # data frame with times where we want predicted rates, space: int> int <- 1 / 100> nd <- data.frame(tfd = seq(0, 10, int))> #> # utility function that calculates the midpoints between values in a vector> mp <- function(x) x[-1] - diff(x) / 2> #> # rates at midpoints (mp) of intervals> lD <- mp(ci.pred(mD, nd)[,1])> lI <- mp(ci.pred(mI, nd)[,1])> lO <- mp(ci.pred(mO, nd)[,1])> #> head(cbind(mp(nd$tfd), lD, lI, lO))
lD
lI
lO
2 0.005 0.07008416 0.3350871 0.4586874
3 0.015 0.06987361 0.3235601 0.4541849

16/ 30Computing the survival function
XZ tS (t) = exp −λc (s) dsc!

!
X= exp −Λc (t)0c> # cumulative rates and survival function at right border of the intervals> LD <- cumsum(lD) * int> LI <- cumsum(lI) * int> LO <- cumsum(lO) * int> Sv <- exp(- LD - LI - LO )> head(cbind(Sv, LD, LI, LO))
Sv
LDLI
LO
2 0.9913986 0.0007008416 0.003350871 0.004586874
3 0.9830309 0.0013995777 0.006586472 0.009128723
4 0.9748877 0.0020962145 0.009710769 0.013625990
5 0.9669607 0.0027907585 0.012727593 0.018079112
6 0.9592415 0.0034832158 0.015640641 0.022488524
7 0.9517225 0.0041735928 0.018453485 0.026854654
17/ 30Computing the cumulative risk functions
Z t
Rc (t) = λc (u)S (u) du0

Note that the survival starts at 1, cumulative risks at 0
> # when integrating to get the cumulative risks we use the average (mp) Sv> 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)
> head(cbind(Sv, rD, rI, rO))
Sv
rDrIrO
1.0000000 0.0000000000 0.000000000 0.000000000
2 0.9913986 0.0006978275 0.003336460 0.004567147
3 0.9830309 0.0013876301 0.006530693 0.009050927
4 0.9748877 0.0020696092 0.009589253 0.013453568
5 0.9669607 0.0027439587 0.012518360 0.017777213
6 0.9592415 0.0034108651 0.015323920 0.022023922
18/ 30Stacking the probabilities withmat2pol, 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",+col = c("black","red","blue","forestgreen"))> text(9, mp(zz["9", ]), c("Dead", "Ins", "OAD"," DM"), col = "white")> box(col = "white", lwd = 3)

19/ 301.0DM
0.8

Probability

0.6
OAD

0.4

Ins0.2Dead0.00246
Time since DM diagnosis (years)

81020/ 30

Confidence intervals. . . do not carry easily over through these calculations
I No analytical expression available (not even approximate)
I Thus, simulation:
I From the parameters of the model we got rates,from these we computed survival and cumulative risks.
I Instead of the parameters from the model use a random
sample from the normal distribution with:III

mean: estimates of parameters
variance: var-cov of parameters

21/ 30Confidence intervals
IInstead of the parameters from the model use a random
sample from the normal distribution with:II

mean: estimates of parameters
variance: var-cov of parameters

Compute rates, survival and cumulative risks from thisI Repeat many times (2000 for example)
I Use the 2.5 and 97.5 percentiles for the 95% c.i.
I Implemented in ci.CriskI

22/ 30Using ci.Crisk> res <- ci.Crisk(list(OAD = mO,+Ins = mI,+
Dead = mD),+
nd = data.frame(tfd = 0:1000 / 100),+nB = 500,+perm = 4:1)
NOTE: Times are assumed to be in the column tfd at equal distances of 0.01
> str(res)
List of 4$ Crisk: num [1:1001, 1:4, 1:3] 1 0.991 0.983 0.975 0.967 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfd : chr [1:1001] "0" "0.01" "0.02" "0.03" ...
.. ..$ cause: chr [1:4] "Surv" "OAD" "Ins" "Dead".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ Srisk: num [1:1001, 1:3, 1:3] 0 0.000703 0.001398 0.002084 0.002764 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfd : chr [1:1001] "0" "0.01" "0.02" "0.03" ...
.. ..$ cause: chr [1:3] "Dead" "Dead+Ins" "Dead+Ins+OAD".. ..$: chr [1:3] "50%" "2.5%" "97.5%"

23/ 30> str(res)
List of 4$ Crisk: num [1:1001, 1:4, 1:3] 1 0.991 0.983 0.975 0.967 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfd : chr [1:1001] "0" "0.01" "0.02" "0.03" ...
.. ..$ cause: chr [1:4] "Surv" "OAD" "Ins" "Dead".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ Srisk: num [1:1001, 1:3, 1:3] 0 0.000703 0.001398 0.002084 0.002764 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfd : chr [1:1001] "0" "0.01" "0.02" "0.03" ...
.. ..$ cause: chr [1:3] "Dead" "Dead+Ins" "Dead+Ins+OAD".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ Stime: num [1:1001, 1:4, 1:3] 0 0.00996 0.01983 0.02962 0.03933 ...
..- attr(*, "dimnames")=List of 3.. ..$ tfd : chr [1:1001] "0" "0.01" "0.02" "0.03" ...
.. ..$ cause: chr [1:4] "Surv" "OAD" "Ins" "Dead".. ..$: chr [1:3] "50%" "2.5%" "97.5%"$ time : num [1:1001] 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 ...
- attr(*, "int")= num 0.01
24/ 30Plotting the results ci.Crisk> matshade(res$time,+cbind(res$Crisk[,"Dead",],+res$Crisk[,"Ins" ,],+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)

25/ 300.4
OAD

Cumulative probability

DeadIns

0.3

0.2

0.1

0.0
02

4
6
Time since DM diagnosis (years)

81026/ 30

Plotting the stacked probabilities from ci.Crisk> 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.3)

27/ 301.0DM0.8

Probability0.6
OAD0.4

Ins0.2
Dead0.00246Time since DM diagnosis (years)

81028/ 30

Competing risks
Different cases of death
I Cause specific rates analyzed separately
there is no such thing as “competing risks analysis”
I With multiple causes, the survival depends on all cause specificrates.
I . . . and so does the cumulative risks
I There is no information about what happen if one cause ofdeath is eliminatedI ci.Crisk:
I Input: rate models and common prediction data frame
I Result: Arrays
I29/ 30

References