The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.