1
#' @title Creates a Weighted Kaplan-Meier plot - svykm.object in survey package
2
#' @description Creates a Weighted Kaplan-Meier plot - svykm.object in survey package
3
#' @param sfit a svykm object
4
#' @param xlabs x-axis label, Default: 'Time-to-event'
5
#' @param ylabs y-axis label.
6
#' @param xlims numeric: list of min and max for x-axis. Default: NULL
7
#' @param ylims numeric: list of min and max for y-axis. Default: c(0, 1)
8
#' @param surv.scale 	scale transformation of survival curves. Allowed values are "default" or "percent".
9
#' @param ystratalabs character list. A list of names for each strata. Default: NULL
10
#' @param ystrataname The legend name. Default: 'Strata'
11
#' @param timeby numeric: control the granularity along the time-axis; defaults to 7 time-points. 
12
#' @param main plot title, Default: ''
13
#' @param pval logical: add the pvalue to the plot?, Default: FALSE
14
#' @param pval.size numeric value specifying the p-value text size. Default is 5.
15
#' @param pval.coord numeric vector, of length 2, specifying the x and y coordinates of the p-value. Default values are NULL
16
#' @param pval.testname logical: add '(Log-rank)' text to p-value. Default = F
17
#' @param legend logical. should a legend be added to the plot? Default: TRUE
18
#' @param ci logical. Should confidence intervals be plotted. Default = NULL
19
#' @param legendposition numeric. x, y position of the legend if plotted. Default: c(0.85, 0.8)
20
#' @param linecols Character. Colour brewer pallettes too colour lines. Default: 'Set1', "black" for black with dashed line.
21
#' @param dashed logical. Should a variety of linetypes be used to identify lines. Default: FALSE
22
#' @param cumhaz Show cumulaive hazard function, Default: F
23
#' @param design Data design for reactive design data , Default: NULL
24
#' @param subs = NULL,
25
#' @param table logical: Create a table graphic below the K-M plot, indicating at-risk numbers?
26
#' @param label.nrisk Numbers at risk label. Default = "Numbers at risk"
27
#' @param size.label.nrisk Font size of label.nrisk. Default = 10
28
#' @param ... PARAM_DESCRIPTION
29
#' @return plot
30
#' @details DETAILS
31
#' @examples 
32
#'  library(survey)
33
#'  data(pbc, package="survival")
34
#'  pbc$randomized <- with(pbc, !is.na(trt) & trt>0)
35
#'  biasmodel <- glm(randomized~age*edema,data=pbc)
36
#'  pbc$randprob <- fitted(biasmodel)
37
#'  dpbc <- svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))
38
#'  s1 <- svykm(Surv(time,status>0)~sex, design=dpbc)
39
#'  svyjskm(s1)
40
#' @rdname svyjskm
41
#' @import ggplot2
42
#' @importFrom stats formula
43
#' @importFrom survey svyranktest
44
#' @importFrom survival Surv
45
#' @export 
46

47
svyjskm <- function(sfit,
48
                    xlabs = "Time-to-event",
49
                    ylabs = "Survival probability",
50
                    xlims = NULL,
51
                    ylims = c(0,1),
52
                    ystratalabs = NULL,
53
                    ystrataname = NULL,
54
                    surv.scale = c("default", "percent"),
55
                    timeby = NULL,
56
                    main = "",
57
                    pval = FALSE,
58
                    pval.size = 5, 
59
                    pval.coord = c(NULL, NULL),
60
                    pval.testname = F,
61
                    legend = TRUE,
62
                    legendposition=c(0.85,0.8),
63
                    ci = NULL,
64
                    linecols="Set1",
65
                    dashed= FALSE,
66
                    cumhaz = F,
67
                    design = NULL,
68
                    subs = NULL,
69
                    table = F,
70
                    label.nrisk = "Numbers at risk",
71
                    size.label.nrisk = 10,
72
                    ...) {
73
  
74 1
  surv <- strata <- lower <- upper <- NULL
75
  
76
  
77 1
  if (is.null(timeby)){
78 1
    if (class(sfit) == "svykmlist"){
79 1
      timeby <- signif(max(sapply(sfit, function(x){max(x$time)}))/7, 1)
80 1
    } else if(class(sfit) == "svykm"){
81
      
82 1
      timeby <- signif(max(sfit$time)/7, 1)
83
      }
84
  }
85
  
86 1
  if (class(sfit) == "svykmlist"){
87 1
    if(is.null(ystrataname)) ystrataname <- as.character(formula(sfit)[[3]])
88 1
    if (is.null(ci)){
89 1
      ci <- "varlog" %in% names(sfit[[1]])
90
    }
91 1
    if (ci){
92 1
      if ("varlog" %in% names(sfit[[1]])){
93 1
        df <- Reduce(rbind, lapply(names(sfit), function(x){data.frame("strata" = x, "time" = sfit[[x]]$time, "surv" = sfit[[x]]$surv, "lower" = pmax(0, exp(log(sfit[[x]]$surv) - 1.96 * sqrt(sfit[[x]]$varlog))), "upper" = pmin(1, exp(log(sfit[[x]]$surv) + 1.96 * sqrt(sfit[[x]]$varlog))))}))
94
      } else{
95 1
        stop("No CI information in svykmlist object. please run svykm with se = T option.")
96
      }
97
    } else{
98 1
      df <- Reduce(rbind, lapply(names(sfit), function(x){data.frame("strata" = x, "time" = sfit[[x]]$time, "surv" = sfit[[x]]$surv)}))
99
    }
100
    
101 1
    df$strata <- as.factor(df$strata)
102 1
    times <- seq(0, max(sapply(sfit, function(x){max(x$time)})), by = timeby)
103 1
    if (is.null(ystratalabs)){
104 1
      ystratalabs <- names(sfit)
105
    }
106 1
    if (is.null(xlims)){
107 1
      xlims <- c(0,max(sapply(sfit, function(x){max(x$time)})))
108
    }
109
    
110 1
  } else if(class(sfit) == "svykm"){
111 1
    if(is.null(ystrataname)) ystrataname <- "Strata"
112 1
    if (is.null(ci)){
113 0
      ci <- "varlog" %in% names(sfit)
114
    }
115 1
    if (ci){
116 1
      if ("varlog" %in% names(sfit)){
117 1
        df <- data.frame("strata" = "All", "time" = sfit$time, "surv" = sfit$surv,  "lower" = pmax(0, exp(log(sfit$surv) - 1.96 * sqrt(sfit$varlog))), "upper" = pmax(0, exp(log(sfit$surv) + 1.96 * sqrt(sfit$varlog))))
118
      } else{
119 1
        stop("No CI information in svykm object. please run svykm with se = T option.")
120
      }
121
    } else{
122 1
      df <- data.frame("strata" = "All", "time" = sfit$time, "surv" = sfit$surv)
123
    }
124
    
125 1
    times <- seq(0, max(sfit$time), by = timeby)
126 1
    if (is.null(ystratalabs)){
127 1
      ystratalabs <- "All"
128
    }
129 1
    if (is.null(xlims)){
130 1
      xlims <- c(0,max(sfit$time))
131
    }
132
  }
133
  
134 1
  m <- max(nchar(ystratalabs))
135
  
136
  
137
  
138
  
139
  
140 1
  if (cumhaz){
141 1
    df$surv <- 1 - df$surv
142 1
    if (ci){
143 1
      df$lower <- 1 - df$upper
144 1
      df$upper <- 1 - df$lower
145
    }
146
    }
147
  
148
  #Final changes to data for survival plot
149 1
  levels(df$strata) <- ystratalabs
150 1
  zeros <- data.frame("strata" = factor(ystratalabs, levels=levels(df$strata)), "time" = 0, "surv" = 1)
151 1
  if (ci){
152 1
    zeros$upper <- 1
153 1
    zeros$lower <- 1
154
  }
155
  
156 1
  if (cumhaz){
157 1
    zeros$surv <- 0
158 1
    if (ci){
159 1
      zeros$lower <- 0
160 1
      zeros$upper <- 0
161
    }
162
  }
163
  
164 1
  df <- rbind(zeros, df)
165 1
  d <- length(levels(df$strata))
166
  
167
  ###################################
168
  # specifying axis parameteres etc #
169
  ###################################
170
  
171 1
  if(dashed == TRUE | linecols == "black"){
172 0
    linetype=c("solid", "dashed", "dotted", "dotdash", "longdash", "twodash", "1F", "F1", "4C88C488", "12345678")
173
  } else {
174 1
    linetype=c("solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid")
175
  }
176
  
177
  # Scale transformation
178
  #::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
179 1
  surv.scale <- match.arg(surv.scale)
180 1
  scale_labels <-  ggplot2::waiver()
181 1
  if (surv.scale == "percent") scale_labels <- scales::percent
182
  
183
  
184 1
  p <- ggplot2::ggplot( df, aes(x=time, y=surv, colour=strata, linetype=strata)) + ggtitle(main)
185 1
  linecols2 <- linecols
186 1
  if (linecols == "black"){
187 0
    linecols <- "Set1"
188 0
    p <- ggplot2::ggplot( df, aes(x=time, y=surv, linetype=strata)) + ggtitle(main)
189
  }
190
  
191
  #Set up theme elements
192 1
  p <- p + theme_bw() +
193 1
    theme(axis.title.x = element_text(vjust = 0.7),
194 1
          panel.grid.minor = element_blank(),
195 1
          axis.line = element_line(size =0.5, colour = "black"),
196 1
          legend.position = legendposition,
197 1
          legend.background = element_rect(fill = NULL),
198 1
          legend.key = element_rect(colour = NA),
199 1
          panel.border = element_blank(),
200 1
          plot.margin = unit(c(0, 1, .5,ifelse(m < 10, 1.5, 2.5)),"lines"),
201 1
          panel.grid.major = element_blank(),
202 1
          axis.line.x = element_line(size = 0.5, linetype = "solid", colour = "black"),
203 1
          axis.line.y = element_line(size = 0.5, linetype = "solid", colour = "black")) +
204 1
    scale_x_continuous(xlabs, breaks = times, limits = xlims) +
205 1
    scale_y_continuous(ylabs, limits = ylims, labels = scale_labels)
206
  
207
  
208
 
209
  #Removes the legend:
210 1
  if(legend == FALSE)
211 0
    p <- p + theme(legend.position="none")
212
  
213
  #Add lines too plot
214 1
  p <- p + geom_step(size = 0.75) +
215 1
    scale_linetype_manual(name = ystrataname, values=linetype) +
216 1
    scale_colour_brewer(name = ystrataname, palette=linecols)
217
  
218
  #Add 95% CI to plot
219 1
  if(ci == TRUE){
220 1
    if (linecols2 == "black"){
221 0
      p <- p +  geom_ribbon(data=df, aes(ymin = lower, ymax = upper), alpha=0.25, colour=NA) 
222
    } else{
223 1
      p <- p +  geom_ribbon(data=df, aes(ymin = lower, ymax = upper, fill = strata), alpha=0.25, colour=NA) + scale_fill_brewer(name = ystrataname, palette=linecols)
224
    } 
225
  }
226
  
227
  ## p-value
228 1
  if(class(sfit) == "svykm") pval <- FALSE
229
  #if(is.null(design)) pval <- FALSE
230
  
231 1
  if(pval == TRUE) {
232 1
    if(is.null(design)){
233 0
      sdiff <- survey::svylogrank(formula(sfit), design = get(as.character(attr(sfit, "call")$design)))
234
    } else{
235 1
      sdiff <- survey::svylogrank(formula(sfit), design = design)
236
    }
237 1
    pvalue <- sdiff[[2]][2]
238
    
239 1
    pvaltxt <- ifelse(pvalue < 0.0001,"p < 0.0001",paste("p =", round(pvalue, 4)))
240 1
    if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)")
241
    
242
    # MOVE P-VALUE LEGEND HERE BELOW [set x and y]
243 1
    if (is.null(pval.coord)){
244 1
      p <- p + annotate("text",x = (as.integer(max(sapply(sfit, function(x){max(x$time)/5})))), y = 0.1 + ylims[1],label = pvaltxt, size  = pval.size)
245
    } else{
246 1
      p <- p + annotate("text",x = pval.coord[1], y = pval.coord[2], label = pvaltxt, size  = pval.size)
247
    }
248
  }
249
  
250
  ## Create a blank plot for place-holding
251 1
  blank.pic <- ggplot(df, aes(time, surv)) +
252 1
    geom_blank() + theme_void() +             ## Remove gray color
253 1
    theme(axis.text.x = element_blank(),axis.text.y = element_blank(),
254 1
          axis.title.x = element_blank(),axis.title.y = element_blank(),
255 1
          axis.ticks = element_blank(),
256 1
          panel.grid.major = element_blank(),panel.border = element_blank())
257
  
258
  ###################################################
259
  # Create table graphic to include at-risk numbers #
260
  ###################################################
261
  
262 1
  n.risk <- NULL
263 1
  if(table == TRUE) {
264
    
265 1
    if(is.null(design)){
266 0
      sfit2 <- survival::survfit(formula(sfit), data = get(as.character(attr(sfit, "call")$design))$variables)
267
    } else{
268 1
      sfit2 <- survival::survfit(formula(sfit), data = design$variables)
269
    }
270
    
271
    #times <- seq(0, max(sfit2$time), by = timeby)
272
    
273 1
    if(is.null(subs)){
274 1
      if(length(levels(summary(sfit2)$strata)) == 0) {
275 0
        subs1 <- 1
276 0
        subs2 <- 1:length(summary(sfit2,censored=T)$time)
277 0
        subs3 <- 1:length(summary(sfit2,times = times,extend = TRUE)$time)
278
      } else {
279 1
        subs1 <- 1:length(levels(summary(sfit2)$strata))
280 1
        subs2 <- 1:length(summary(sfit2,censored=T)$strata)
281 1
        subs3 <- 1:length(summary(sfit2,times = times,extend = TRUE)$strata)
282
      }
283
    } else{
284 0
      for(i in 1:length(subs)){
285 0
        if(i==1){
286 0
          ssvar <- paste("(?=.*\\b=",subs[i],sep="")
287
        }
288 0
        if(i==length(subs)){
289 0
          ssvar <- paste(ssvar,"\\b)(?=.*\\b=",subs[i],"\\b)",sep="")
290
        }
291 0
        if(!i %in% c(1, length(subs))){
292 0
          ssvar <- paste(ssvar,"\\b)(?=.*\\b=",subs[i],sep="")
293
        }
294 0
        if(i==1 & i==length(subs)){
295 0
          ssvar <- paste("(?=.*\\b=",subs[i],"\\b)",sep="")
296
        }
297
      }
298 0
      subs1 <- which(regexpr(ssvar,levels(summary(sfit2)$strata), perl=T)!=-1)
299 0
      subs2 <- which(regexpr(ssvar,summary(sfit2,censored=T)$strata, perl=T)!=-1)
300 0
      subs3 <- which(regexpr(ssvar,summary(sfit2,times = times,extend = TRUE)$strata, perl=T)!=-1)
301
    }
302
    
303 0
    if(!is.null(subs)) pval <- FALSE
304
    
305

306
    
307 1
    if(length(levels(summary(sfit2)$strata)) == 0) {
308 0
      Factor <- factor(rep("All",length(subs3)))
309
    } else {
310 1
      Factor <- factor(summary(sfit2,times = times,extend = TRUE)$strata[subs3])
311
    }
312
    
313
    
314 1
    risk.data <- data.frame(
315 1
      strata = Factor,
316 1
      time = summary(sfit2,times = times,extend = TRUE)$time[subs3],
317 1
      n.risk = summary(sfit2,times = times,extend = TRUE)$n.risk[subs3]
318
    )
319
    
320
    
321 1
    risk.data$strata <- factor(risk.data$strata, levels=rev(levels(risk.data$strata)))
322
    
323 1
    data.table <- ggplot(risk.data,aes(x = time, y = strata, label = format(n.risk, nsmall = 0))) + 
324 1
      geom_text(size = 3.5) + theme_bw() +
325 1
      scale_y_discrete(breaks = as.character(levels(risk.data$strata)),
326 1
                       labels = rev(ystratalabs)) +
327 1
      scale_x_continuous(label.nrisk, limits = xlims) +
328 1
      theme(axis.title.x = element_text(size = size.label.nrisk, vjust = 1),
329 1
            panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
330 1
            panel.border = element_blank(),axis.text.x = element_blank(),
331 1
            axis.ticks = element_blank(),axis.text.y = element_text(face = "bold",hjust = 1)) 
332 1
    data.table <- data.table +
333 1
      theme(legend.position = "none") + xlab(NULL) + ylab(NULL)
334
    
335
    
336
    # ADJUST POSITION OF TABLE FOR AT RISK
337 1
    data.table <- data.table +
338 1
      theme(plot.margin = unit(c(-1.5, 1, 0.1, ifelse(m < 10, 2.5, 3.5) - 0.15 * m), "lines"))
339
  }
340
  
341
  #######################
342
  # Plotting the graphs #
343
  #######################
344
  
345 1
  if(table == TRUE){
346 1
    grid.arrange(p, blank.pic, data.table, clip = FALSE, nrow = 3,
347 1
                 ncol = 1, heights = unit(c(2, .1, .25),c("null", "null", "null")))
348
  } else {
349 1
    p
350
  }
351

352
  
353
}
354

355

356

Read our documentation on viewing source code .

Loading