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 cut.landmark cut-off for landmark analysis, Default = NULL
29
#' @param ... PARAM_DESCRIPTION
30
#' @return plot
31
#' @details DETAILS
32
#' @examples 
33
#'  library(survey)
34
#'  data(pbc, package="survival")
35
#'  pbc$randomized <- with(pbc, !is.na(trt) & trt>0)
36
#'  biasmodel <- glm(randomized~age*edema,data=pbc)
37
#'  pbc$randprob <- fitted(biasmodel)
38
#'  dpbc <- svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))
39
#'  s1 <- svykm(Surv(time,status>0)~sex, design=dpbc)
40
#'  svyjskm(s1)
41
#' @rdname svyjskm
42
#' @import ggplot2
43
#' @importFrom stats formula
44
#' @importFrom survey svyranktest
45
#' @importFrom survival Surv
46
#' @export 
47

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

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

356
  
357
}
358

359

360

Read our documentation on viewing source code .

Loading