jskm

Jinseob Kim

2024-12-25

Install

install.packages("devtools")
library(devtools)
install_github("jinseob2kim/jskm")
library(jskm)

Example

Survival probability

# Load dataset
library(survival)
data(colon)
#> Warning in data(colon): data set 'colon' not found
fit <- survfit(Surv(time, status) ~ rx, data = colon)

# Plot the data
jskm(fit)

jskm(fit,
  table = T, pval = T, label.nrisk = "No. at risk", size.label.nrisk = 8,
  xlabs = "Time(Day)", ylabs = "Survival", ystratalabs = c("Obs", "Lev", "Lev + 5FU"), ystrataname = "rx",
  marks = F, timeby = 365, xlims = c(0, 3000), ylims = c(0.25, 1), showpercent = T
)
#>        strata time n.risk
#> 1      rx=Obs    0    630
#> 2      rx=Obs  365    519
#> 3      rx=Obs  730    417
#> 4      rx=Obs 1095    360
#> 5      rx=Obs 1460    318
#> 6      rx=Obs 1825    288
#> 7      rx=Obs 2190    187
#> 8      rx=Obs 2555     75
#> 9      rx=Obs 2920     13
#> 10     rx=Obs 3285      0
#> 11     rx=Lev    0    620
#> 12     rx=Lev  365    502
#> 13     rx=Lev  730    406
#> 14     rx=Lev 1095    348
#> 15     rx=Lev 1460    319
#> 16     rx=Lev 1825    299
#> 17     rx=Lev 2190    201
#> 18     rx=Lev 2555     87
#> 19     rx=Lev 2920     13
#> 20     rx=Lev 3285      4
#> 21 rx=Lev+5FU    0    608
#> 22 rx=Lev+5FU  365    531
#> 23 rx=Lev+5FU  730    453
#> 24 rx=Lev+5FU 1095    420
#> 25 rx=Lev+5FU 1460    391
#> 26 rx=Lev+5FU 1825    361
#> 27 rx=Lev+5FU 2190    247
#> 28 rx=Lev+5FU 2555    102
#> 29 rx=Lev+5FU 2920     24
#> 30 rx=Lev+5FU 3285      4
#> Warning: Removed 16 rows containing missing values or values outside the scale range
#> (`geom_step()`).
#> Warning: Removed 3 rows containing missing values or values outside the scale range
#> (`geom_text()`).

Cumulative incidence: 1- Survival probability

jskm(fit, ci = T, cumhaz = T, mark = F, ylab = "Cumulative incidence (%)", surv.scale = "percent", pval = T, pval.size = 6, pval.coord = c(300, 0.7))

Landmark analysis

jskm(fit, mark = F, surv.scale = "percent", pval = T, table = T, cut.landmark = 500, showpercent = T)
#>        strata time n.risk
#> 1      rx=Obs    0    630
#> 2      rx=Obs  500    470
#> 3      rx=Obs 1000    372
#> 4      rx=Obs 1500    315
#> 5      rx=Obs 2000    256
#> 6      rx=Obs 2500     90
#> 7      rx=Obs 3000     11
#> 8      rx=Lev    0    620
#> 9      rx=Lev  500    464
#> 10     rx=Lev 1000    360
#> 11     rx=Lev 1500    318
#> 12     rx=Lev 2000    266
#> 13     rx=Lev 2500    107
#> 14     rx=Lev 3000      8
#> 15 rx=Lev+5FU    0    608
#> 16 rx=Lev+5FU  500    498
#> 17 rx=Lev+5FU 1000    425
#> 18 rx=Lev+5FU 1500    387
#> 19 rx=Lev+5FU 2000    328
#> 20 rx=Lev+5FU 2500    127
#> 21 rx=Lev+5FU 3000     14

Competing risk analysis

status2 variable: 0 - censoring, 1 - event, 2 - competing risk

## Make competing risk variable, Not real
colon$status2 <- colon$status
colon$status2[1:400] <- 2
colon$status2 <- factor(colon$status2)
fit2 <- survfit(Surv(time, status2) ~ rx, data = colon)
jskm(fit2, mark = F, surv.scale = "percent", table = T, status.cmprsk = "1")
#>        strata time n.risk
#> 1      rx=Obs    0    630
#> 2      rx=Obs  500    470
#> 3      rx=Obs 1000    372
#> 4      rx=Obs 1500    315
#> 5      rx=Obs 2000    256
#> 6      rx=Obs 2500     90
#> 7      rx=Obs 3000     11
#> 8      rx=Lev    0    620
#> 9      rx=Lev  500    464
#> 10     rx=Lev 1000    360
#> 11     rx=Lev 1500    318
#> 12     rx=Lev 2000    266
#> 13     rx=Lev 2500    107
#> 14     rx=Lev 3000      8
#> 15 rx=Lev+5FU    0    608
#> 16 rx=Lev+5FU  500    498
#> 17 rx=Lev+5FU 1000    425
#> 18 rx=Lev+5FU 1500    387
#> 19 rx=Lev+5FU 2000    328
#> 20 rx=Lev+5FU 2500    127
#> 21 rx=Lev+5FU 3000     14

jskm(fit2, mark = F, surv.scale = "percent", table = T, status.cmprsk = "1", showpercent = T, cut.landmark = 500)
#>        strata time n.risk
#> 1      rx=Obs    0    630
#> 2      rx=Obs  500    470
#> 3      rx=Obs 1000    372
#> 4      rx=Obs 1500    315
#> 5      rx=Obs 2000    256
#> 6      rx=Obs 2500     90
#> 7      rx=Obs 3000     11
#> 8      rx=Lev    0    620
#> 9      rx=Lev  500    464
#> 10     rx=Lev 1000    360
#> 11     rx=Lev 1500    318
#> 12     rx=Lev 2000    266
#> 13     rx=Lev 2500    107
#> 14     rx=Lev 3000      8
#> 15 rx=Lev+5FU    0    608
#> 16 rx=Lev+5FU  500    498
#> 17 rx=Lev+5FU 1000    425
#> 18 rx=Lev+5FU 1500    387
#> 19 rx=Lev+5FU 2000    328
#> 20 rx=Lev+5FU 2500    127
#> 21 rx=Lev+5FU 3000     14

Theme Jama

jskm(fit, theme = "jama", cumhaz = T, table = T, mark = F, ylab = "Cumulative incidence (%)", surv.scale = "percent", pval = T, pval.size = 6, pval.coord = c(300, 0.7))
#>        strata time n.risk
#> 1      rx=Obs    0    630
#> 2      rx=Obs  500    470
#> 3      rx=Obs 1000    372
#> 4      rx=Obs 1500    315
#> 5      rx=Obs 2000    256
#> 6      rx=Obs 2500     90
#> 7      rx=Obs 3000     11
#> 8      rx=Lev    0    620
#> 9      rx=Lev  500    464
#> 10     rx=Lev 1000    360
#> 11     rx=Lev 1500    318
#> 12     rx=Lev 2000    266
#> 13     rx=Lev 2500    107
#> 14     rx=Lev 3000      8
#> 15 rx=Lev+5FU    0    608
#> 16 rx=Lev+5FU  500    498
#> 17 rx=Lev+5FU 1000    425
#> 18 rx=Lev+5FU 1500    387
#> 19 rx=Lev+5FU 2000    328
#> 20 rx=Lev+5FU 2500    127
#> 21 rx=Lev+5FU 3000     14

Theme Nejmoa

jskm(fit, theme = "nejm", nejm.infigure.ratiow = 0.6, nejm.infigure.ratioh = 0.4, nejm.infigure.ylim = c(0, 0.7), cumhaz = T, table = T, mark = F, ylab = "Cumulative incidence (%)", surv.scale = "percent", pval = T, pval.size = 6, pval.coord = c(300, 0.7))
#>        strata time n.risk
#> 1      rx=Obs    0    630
#> 2      rx=Obs  500    470
#> 3      rx=Obs 1000    372
#> 4      rx=Obs 1500    315
#> 5      rx=Obs 2000    256
#> 6      rx=Obs 2500     90
#> 7      rx=Obs 3000     11
#> 8      rx=Lev    0    620
#> 9      rx=Lev  500    464
#> 10     rx=Lev 1000    360
#> 11     rx=Lev 1500    318
#> 12     rx=Lev 2000    266
#> 13     rx=Lev 2500    107
#> 14     rx=Lev 3000      8
#> 15 rx=Lev+5FU    0    608
#> 16 rx=Lev+5FU  500    498
#> 17 rx=Lev+5FU 1000    425
#> 18 rx=Lev+5FU 1500    387
#> 19 rx=Lev+5FU 2000    328
#> 20 rx=Lev+5FU 2500    127
#> 21 rx=Lev+5FU 3000     14
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.

Weighted Kaplan-Meier plot - svykm.object in survey package

library(survey)
#> Loading required package: grid
#> Loading required package: Matrix
#> 
#> Attaching package: 'survey'
#> The following object is masked from 'package:graphics':
#> 
#>     dotchart
data(pbc, package = "survival")
pbc$randomized <- with(pbc, !is.na(trt) & trt > 0)
biasmodel <- glm(randomized ~ age * edema, data = pbc)
pbc$randprob <- fitted(biasmodel)

dpbc <- svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, randomized))

s1 <- svykm(Surv(time, status > 0) ~ 1, design = dpbc)
s2 <- svykm(Surv(time, status > 0) ~ sex, design = dpbc)

svyjskm(s1)

svyjskm(s2)

svyjskm(s2, cumhaz = T, ylab = "Cumulative incidence(%)", surv.scale = "percent", showpercent = T)

If you want to get confidence interval, you should apply se = T option to svykm object.

s3 <- svykm(Surv(time, status > 0) ~ sex, design = dpbc, se = T)
svyjskm(s3)

svyjskm(s3, ci = F, showpercent = T)

svyjskm(s3, ci = F, surv.scale = "percent", pval = T, table = T, cut.landmark = 1000)

Theme

JAMA

svyjskm(s2, theme = "jama", pval = T, table = T, design = dpbc)

NEJM

svyjskm(s2, theme = "nejm", nejm.infigure.ratiow = 0.4, nejm.infigure.ratioh = 0.4, nejm.infigure.ylim = c(0.2, 1), pval = T, table = T, design = dpbc)