@@ -85,19 +85,65 @@
Loading
85 85
      }
86 86
  }
87 87
  
88 -
  if (class(sfit) == "svykmlist"){
89 -
    if(is.null(ystrataname)) ystrataname <- as.character(formula(sfit)[[3]])
90 -
    if (is.null(ci)){
88 +
  if (is.null(ci)){
89 +
    if (class(sfit) == "svykmlist"){
91 90
      ci <- "varlog" %in% names(sfit[[1]])
91 +
    } else if (class(sfit) == "svykm"){
92 +
      ci <- "varlog" %in% names(sfit)
92 93
    }
94 +
  }
95 +
  
96 +
  
97 +
  if (ci & !is.null(cut.landmark)){
98 +
    if (is.null(design)){
99 +
      design <- tryCatch(get(as.character(attr(sfit, "call")$design)), error = function(e) e)
100 +
      if ("error" %in% class(design)){
101 +
        stop("'pval' option requires design object. please input 'design' option") 
102 +
      }
103 +
    }
104 +
    var.time <- as.character(formula(sfit)[[2]][[2]])
105 +
    var.event <- as.character(formula(sfit)[[2]][[3]])
106 +
    if (length(var.event) > 1){
107 +
      var.event <- setdiff(var.event, as.character(as.symbol(var.event)))
108 +
      var.event <- var.event[sapply(var.event, function(x) {"warning" %in% class(tryCatch(as.numeric(x), warning = function(w) w))})]
109 +
    }
110 +
    design1 <- design
111 +
    design1$variables[[var.event]][design1$variables[[var.time]] >= cut.landmark] <- 0
112 +
    design1$variables[[var.time]][design1$variables[[var.time]] >= cut.landmark] <- cut.landmark
113 +
    
114 +
    sfit2 <- survey::svykm(formula(sfit), design = subset(design, get(var.time) >= cut.landmark), se = T)
115 +
  }
116 +
  
117 +
  
118 +
  
119 +
  
120 +
  if (class(sfit) == "svykmlist"){
121 +
    if(is.null(ystrataname)) ystrataname <- as.character(formula(sfit)[[3]])
122 +
93 123
    if (ci){
94 124
      if ("varlog" %in% names(sfit[[1]])){
95 -
        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))))}))
125 +
        df <- do.call(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))))}))
126 +
        if (!is.null(cut.landmark)){
127 +
          df2 <- do.call(rbind, lapply(names(sfit2), function(x){data.frame("strata" = x, "time" = sfit2[[x]]$time, "surv" = sfit2[[x]]$surv, "lower" = pmax(0, exp(log(sfit2[[x]]$surv) - 1.96 * sqrt(sfit2[[x]]$varlog))), "upper" = pmin(1, exp(log(sfit2[[x]]$surv) + 1.96 * sqrt(sfit2[[x]]$varlog))))}))
128 +
          df <- rbind(df[df$time < cut.landmark, ], data.frame("strata" = unique(df$strata), "time" = cut.landmark, "surv" = 1, "lower" = 1, "upper" = 1), df2)
129 +
        }
130 +
        
96 131
      } else{
97 132
        stop("No CI information in svykmlist object. please run svykm with se = T option.")
98 133
      }
99 134
    } else{
100 -
      df <- Reduce(rbind, lapply(names(sfit), function(x){data.frame("strata" = x, "time" = sfit[[x]]$time, "surv" = sfit[[x]]$surv)}))
135 +
      df <- do.call(rbind, lapply(names(sfit), function(x){data.frame("strata" = x, "time" = sfit[[x]]$time, "surv" = sfit[[x]]$surv)}))
136 +
      if (!is.null(cut.landmark)){
137 +
        for (v in unique(df$strata)){
138 +
          if (nrow(subset(df, time == cut.landmark & strata == v)) == 0){
139 +
            df <- rbind(df, data.frame(strata = v, time = cut.landmark, surv = 1))
140 +
          } else{
141 +
            df[df$time == cut.landmark & df$strata == v, "surv"] <- 1
142 +
          }
143 +
          
144 +
          df[df$time > cut.landmark & df$strata == v, "surv"] <- df[df$time > cut.landmark & df$strata == v, "surv"]/min(df[df$time < cut.landmark & df$strata == v, "surv"])
145 +
        }
146 +
      }
101 147
    }
102 148
    
103 149
    df$strata <- as.factor(df$strata)
@@ -111,17 +157,28 @@
Loading
111 157
    
112 158
  } else if(class(sfit) == "svykm"){
113 159
    if(is.null(ystrataname)) ystrataname <- "Strata"
114 -
    if (is.null(ci)){
115 -
      ci <- "varlog" %in% names(sfit)
116 -
    }
160 +
117 161
    if (ci){
118 162
      if ("varlog" %in% names(sfit)){
119 163
        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))))
164 +
        if (!is.null(cut.landmark)){
165 +
          df2 <- data.frame("strata" = "All", "time" = sfit2$time, "surv" = sfit2$surv,  "lower" = pmax(0, exp(log(sfit2$surv) - 1.96 * sqrt(sfit2$varlog))), "upper" = pmax(0, exp(log(sfit2$surv) + 1.96 * sqrt(sfit2$varlog))))
166 +
          df <- rbind(df[df$time < cut.landmark, ], data.frame("strata" = "All", "time" = cut.landmark, "surv" = 1, "lower" = 1, "upper" = 1), df2)
167 +
        }
120 168
      } else{
121 169
        stop("No CI information in svykm object. please run svykm with se = T option.")
122 170
      }
123 171
    } else{
124 172
      df <- data.frame("strata" = "All", "time" = sfit$time, "surv" = sfit$surv)
173 +
      if (!is.null(cut.landmark)){
174 +
        if (nrow(subset(df, time == cut.landmark)) == 0){
175 +
          df <- rbind(df, data.frame(strata = "All", time = cut.landmark, surv = 1))
176 +
        } else{
177 +
          df[df$time == cut.landmark, "surv"] <- 1
178 +
        }
179 +
        
180 +
        df[df$time > cut.landmark, "surv"] <- df[df$time > cut.landmark, "surv"]/min(df[df$time < cut.landmark, "surv"])
181 +
      }
125 182
    }
126 183
    
127 184
    times <- seq(0, max(sfit$time), by = timeby)
@@ -220,7 +277,7 @@
Loading
220 277
    scale_colour_brewer(name = ystrataname, palette=linecols)
221 278
  
222 279
  #Add 95% CI to plot
223 -
  if(ci == TRUE){
280 +
  if(ci){
224 281
    if (linecols2 == "black"){
225 282
      p <- p +  geom_ribbon(data=df, aes(ymin = lower, ymax = upper), alpha=0.25, colour=NA) 
226 283
    } else{
@@ -228,27 +285,67 @@
Loading
228 285
    } 
229 286
  }
230 287
  
288 +
  if (!is.null(cut.landmark)){
289 +
    p <- p + geom_vline(xintercept = cut.landmark, lty = 2)
290 +
  }
291 +
  
231 292
  ## p-value
232 293
  if(class(sfit) == "svykm") pval <- FALSE
233 294
  #if(is.null(design)) pval <- FALSE
234 295
  
235 -
  if(pval == TRUE) {
236 -
    if(is.null(design)){
237 -
      sdiff <- survey::svylogrank(formula(sfit), design = get(as.character(attr(sfit, "call")$design)))
238 -
    } else{
239 -
      sdiff <- survey::svylogrank(formula(sfit), design = design)
296 +
  if(pval) {
297 +
    if (is.null(design)){
298 +
      design <- tryCatch(get(as.character(attr(sfit, "call")$design)), error = function(e) e)
299 +
      if ("error" %in% class(design)){
300 +
        stop("'pval' option requires design object. please input 'design' option") 
301 +
      }
240 302
    }
241 -
    pvalue <- sdiff[[2]][2]
242 303
    
243 -
    pvaltxt <- ifelse(pvalue < 0.001,"p < 0.001",paste("p =", round(pvalue, 3)))
244 -
    if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)")
245 -
    
246 -
    # MOVE P-VALUE LEGEND HERE BELOW [set x and y]
247 -
    if (is.null(pval.coord)){
248 -
      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)
304 +
    if (is.null(cut.landmark)){
305 +
      sdiff <- survey::svylogrank(formula(sfit), design = design)
306 +
      pvalue <- sdiff[[2]][2]
307 +
      
308 +
      pvaltxt <- ifelse(pvalue < 0.001, "p < 0.001", paste("p =", round(pvalue, 3)))
309 +
      if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)")
310 +
      
311 +
      # MOVE P-VALUE LEGEND HERE BELOW [set x and y]
312 +
      if (is.null(pval.coord)){
313 +
        p <- p + annotate("text",x = as.integer(max(sapply(sfit, function(x){max(x$time)/10}))), y = 0.1 + ylims[1],label = pvaltxt, size  = pval.size)
314 +
      } else{
315 +
        p <- p + annotate("text",x = pval.coord[1], y = pval.coord[2], label = pvaltxt, size  = pval.size)
316 +
      } 
249 317
    } else{
250 -
      p <- p + annotate("text",x = pval.coord[1], y = pval.coord[2], label = pvaltxt, size  = pval.size)
318 +
      if (is.null(design)){
319 +
        design <- tryCatch(get(as.character(attr(sfit, "call")$design)), error = function(e) e)
320 +
        if ("error" %in% class(design)){
321 +
          stop("'pval' option requires design object. please input 'design' option") 
322 +
        }
323 +
      }
324 +
      var.time <- as.character(formula(sfit)[[2]][[2]])
325 +
      var.event <- as.character(formula(sfit)[[2]][[3]])
326 +
      if (length(var.event) > 1){
327 +
        var.event <- setdiff(var.event, as.character(as.symbol(var.event)))
328 +
        var.event <- var.event[sapply(var.event, function(x) {"warning" %in% class(tryCatch(as.numeric(x), warning = function(w) w))})]
329 +
      }
330 +
      design1 <- design
331 +
      design1$variables[[var.event]][design1$variables[[var.time]] >= cut.landmark] <- 0
332 +
      design1$variables[[var.time]][design1$variables[[var.time]] >= cut.landmark] <- cut.landmark
333 +
      
334 +
      sdiff1 <- survey::svylogrank(formula(sfit), design = design1)
335 +
      sdiff2 <- survey::svylogrank(formula(sfit), design = subset(design, get(var.time) >= cut.landmark))
336 +
      pvalue <- sapply(list(sdiff1, sdiff2), function(x){x[[2]][2]})
337 +
      
338 +
      pvaltxt <- ifelse(pvalue < 0.001, "p < 0.001", paste("p =", round(pvalue, 3)))
339 +
      if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)")
340 +
      
341 +
      if (is.null(pval.coord)){
342 +
        p <- p + annotate("text",x = c(as.integer(max(sapply(sfit, function(x){max(x$time)/10}))), as.integer(max(sapply(sfit, function(x){max(x$time)/10}))) + cut.landmark), y = 0.1 + ylims[1],label = pvaltxt, size  = pval.size)
343 +
      } else{
344 +
        p <- p + annotate("text",x = c(pval.coord[1], pval.coord[1] + cut.landmark), y = pval.coord[2], label = pvaltxt, size  = pval.size)
345 +
      }
346 +
      
251 347
    }
348 +
    
252 349
  }
253 350
  
254 351
  ## Create a blank plot for place-holding

@@ -185,22 +185,31 @@
Loading
185 185
    upper = sfit$upper[subs2],
186 186
    lower = sfit$lower[subs2]
187 187
  )
188 +
  
189 +
  form <- sfit$call$formula
188 190
189 191
  if (!is.null(cut.landmark)){
190 192
    if (is.null(data)){
191 -
      stop("Landmark analysis requires data object. please check 'data' option")
193 +
      data <- tryCatch(eval(sfit$call$data), error = function(e) e)
194 +
      if ("error" %in% class(data)){
195 +
        stop("Landmark analysis requires data object. please input 'data' option") 
196 +
      }
192 197
    }
193 198
    
194 -
    var.time <- as.character(sfit$call$formula[[2]][[2]])
195 -
    var.event <- as.character(sfit$call$formula[[2]][[3]])
199 +
    var.time <- as.character(form[[2]][[2]])
200 +
    var.event <- as.character(form[[2]][[3]])
201 +
    if (length(var.event) > 1){
202 +
      var.event <- setdiff(var.event, as.character(as.symbol(var.event)))
203 +
      var.event <- var.event[sapply(var.event, function(x) {"warning" %in% class(tryCatch(as.numeric(x), warning = function(w) w))})]
204 +
    }
196 205
    data1 <- data
197 206
    data1[[var.event]][data1[[var.time]] >= cut.landmark] <- 0
198 207
    data1[[var.time]][data1[[var.time]] >= cut.landmark] <- cut.landmark
199 208
  
200 -
    sfit1 <- survfit(as.formula(sfit$call$formula), data1)
201 -
    sfit2 <- survfit(as.formula(sfit$call$formula), data[data[[var.time]] >= cut.landmark, ])
209 +
    sfit1 <- survfit(as.formula(form), data1)
210 +
    sfit2 <- survfit(as.formula(form), data[data[[var.time]] >= cut.landmark, ])
202 211
    
203 -
    if (levels(Factor) == "All"){
212 +
    if (length(levels(Factor)) == 1){
204 213
      df2 <- merge(subset(df, time >= cut.landmark)[, c("time", "n.risk", "n.event", "n.censor", "strata")], 
205 214
                   data.frame(time = sfit2$time, surv = sfit2$surv, strata = "All", upper = sfit2$upper, lower = sfit2$lower), 
206 215
                   by = c("time", "strata"))
@@ -228,7 +237,7 @@
Loading
228 237
  
229 238
  #Final changes to data for survival plot
230 239
  levels(df$strata) <- ystratalabs
231 -
  zeros <- data.frame(time = 0, surv = 1,
240 +
  zeros <- data.frame(time = 0, n.risk = NA, n.event = NA, n.censor = NA, surv = 1,
232 241
                      strata = factor(ystratalabs, levels=levels(df$strata)),
233 242
                      upper = 1, lower = 1)
234 243
  if (cumhaz){
@@ -237,7 +246,7 @@
Loading
237 246
    zeros$upper = 0
238 247
  }
239 248
  
240 -
  df <- plyr::rbind.fill(zeros, df)
249 +
  df <- rbind(zeros, df)
241 250
  d <- length(levels(df$strata))
242 251
  
243 252
  ###################################
@@ -322,40 +331,73 @@
Loading
322 331
  #####################a
323 332
  
324 333
  if(length(levels(summary(sfit)$strata)) == 0) pval <- F
325 -
  if(!is.null(cut.landmark)) pval <- F
334 +
  #if(!is.null(cut.landmark)) pval <- F
326 335
  
327 336
  if(pval == TRUE) {
328 337
    if (is.null(data)){
329 -
      data = eval(sfit$call$data)
330 -
    }
331 -
  
332 -
    sdiff <- survival::survdiff(as.formula(sfit$call$formula), data = data)
333 -
    pvalue <- pchisq(sdiff$chisq,length(sdiff$n) - 1,lower.tail = FALSE)
334 -
    
335 -
    ## cluster option
336 -
    if (cluster.option == "cluster" & !is.null(cluster.var)){
337 -
      form.old <- as.character(sfit$call$formula)
338 -
      form.new <- paste(form.old[2], form.old[1], " + ", form.old[3], " + cluster(", cluster.var, ")", sep="")
339 -
      sdiff <- survival::coxph(as.formula(form.new), data = data, model = T, robust = T)
340 -
      pvalue <- summary(sdiff)$robscore["pvalue"]
341 -
    } else if (cluster.option == "frailty" & !is.null(cluster.var)){
342 -
      form.old <- as.character(sfit$call$formula)
343 -
      form.new <- paste(form.old[2], form.old[1], " + ", form.old[3], " + frailty(", cluster.var, ")", sep="")
344 -
      sdiff <- survival::coxph(as.formula(form.new), data =data, model = T)
345 -
      pvalue <- summary(sdiff)$logtest["pvalue"]
338 +
      data <- tryCatch(eval(sfit$call$data), error = function(e) e)
339 +
      if ("error" %in% class(data)){
340 +
        stop("'pval' option requires data object. please input 'data' option") 
341 +
      }
346 342
    }
347 343
    
348 -
    pvaltxt <- ifelse(pvalue < 0.001,"p < 0.001",paste("p =", round(pvalue, 3)))
349 -
    
350 -
    if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)")
351 -
    
352 -
    # MOVE P-VALUE LEGEND HERE BELOW [set x and y]
353 -
    if (is.null(pval.coord)){
354 -
      p <- p + annotate("text",x = (as.integer(max(sfit$time)/5)), y = 0.1 + ylims[1],label = pvaltxt, size  = pval.size)
355 -
    } else{
356 -
      p <- p + annotate("text",x = pval.coord[1], y = pval.coord[2], label = pvaltxt, size  = pval.size)
344 +
    if (is.null(cut.landmark)){
345 +
      sdiff <- survival::survdiff(as.formula(form), data = data)
346 +
      pvalue <- pchisq(sdiff$chisq,length(sdiff$n) - 1,lower.tail = FALSE)
347 +
      
348 +
      ## cluster option
349 +
      if (cluster.option == "cluster" & !is.null(cluster.var)){
350 +
        form.old <- as.character(form)
351 +
        form.new <- paste(form.old[2], form.old[1], " + ", form.old[3], " + cluster(", cluster.var, ")", sep="")
352 +
        sdiff <- survival::coxph(as.formula(form.new), data = data, model = T, robust = T)
353 +
        pvalue <- summary(sdiff)$robscore["pvalue"]
354 +
      } else if (cluster.option == "frailty" & !is.null(cluster.var)){
355 +
        form.old <- as.character(form)
356 +
        form.new <- paste(form.old[2], form.old[1], " + ", form.old[3], " + frailty(", cluster.var, ")", sep="")
357 +
        sdiff <- survival::coxph(as.formula(form.new), data =data, model = T)
358 +
        pvalue <- summary(sdiff)$logtest["pvalue"]
359 +
      }
360 +
      
361 +
      pvaltxt <- ifelse(pvalue < 0.001, "p < 0.001", paste("p =", round(pvalue, 3)))
362 +
      if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)")
363 +
      
364 +
      # MOVE P-VALUE LEGEND HERE BELOW [set x and y]
365 +
      if (is.null(pval.coord)){
366 +
        p <- p + annotate("text",x = (as.integer(max(sfit$time)/5)), y = 0.1 + ylims[1],label = pvaltxt, size  = pval.size)
367 +
      } else{
368 +
        p <- p + annotate("text",x = pval.coord[1], y = pval.coord[2], label = pvaltxt, size  = pval.size)
369 +
      }
370 +
    } else {
371 +
      sdiff1 <- survival::survdiff(as.formula(form), data1)
372 +
      sdiff2 <- survival::survdiff(as.formula(form), data[data[[var.time]] >= cut.landmark, ])
373 +
      pvalue <- sapply(list(sdiff1, sdiff2), function(x){pchisq(x$chisq,length(x$n) - 1,lower.tail = FALSE)})
374 +
      
375 +
      ## cluster option
376 +
      if (cluster.option == "cluster" & !is.null(cluster.var)){
377 +
        form.old <- as.character(form)
378 +
        form.new <- paste(form.old[2], form.old[1], " + ", form.old[3], " + cluster(", cluster.var, ")", sep="")
379 +
        sdiff1 <- survival::coxph(as.formula(form.new), data = data1, model = T, robust = T)
380 +
        sdiff2 <- survival::coxph(as.formula(form.new), data = data[data[[var.time]] >= cut.landmark, ], model = T, robust = T)
381 +
        pvalue <- sapply(list(sdiff1, sdiff2), function(x){summary(x)$robscore["pvalue"]})
382 +
      } else if (cluster.option == "frailty" & !is.null(cluster.var)){
383 +
        form.old <- as.character(form)
384 +
        form.new <- paste(form.old[2], form.old[1], " + ", form.old[3], " + frailty(", cluster.var, ")", sep="")
385 +
        sdiff1 <- survival::coxph(as.formula(form.new), data = data1, model = T)
386 +
        sdiff2 <- survival::coxph(as.formula(form.new), data = data[data[[var.time]] >= cut.landmark, ], model = T)
387 +
        pvalue <- sapply(list(sdiff1, sdiff2), function(x){summary(x)$logtest["pvalue"]})
388 +
      }
389 +
      
390 +
      pvaltxt <- ifelse(pvalue < 0.001, "p < 0.001", paste("p =", round(pvalue, 3)))
391 +
      
392 +
      if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)")
393 +
      
394 +
      if (is.null(pval.coord)){
395 +
        p <- p + annotate("text",x = c(as.integer(max(sfit$time)/10), as.integer(max(sfit$time)/10) + cut.landmark), y = 0.1 + ylims[1],label = pvaltxt, size  = pval.size)
396 +
      } else{
397 +
        p <- p + annotate("text",x = c(pval.coord[1], pval.coord[1] + cut.landmark), y = pval.coord[2], label = pvaltxt, size  = pval.size)
398 +
      }
357 399
    }
358 -
    
400 +
  
359 401
  }
360 402
  
361 403
  ###################################################
Files Coverage
R 63.99%
Project Totals (2 files) 63.99%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading