From 0a59f2cc9749ab8141ee7b7756edf4147eb20928 Mon Sep 17 00:00:00 2001 From: hyojong myung Date: Tue, 22 Oct 2024 16:10:28 +0900 Subject: [PATCH 1/2] add median expression to the graph --- NEWS.md | 8 + R/jskm.R | 320 +++++++++++++++++++++++------------ R/svyjskm.R | 356 ++++++++++++++++++++++++++------------- README.Rmd | 2 +- README.md | 12 +- Rplot.pdf | Bin 0 -> 25891 bytes tests/testthat.R | 2 +- tests/testthat/test-km.R | 3 +- 8 files changed, 476 insertions(+), 227 deletions(-) create mode 100644 Rplot.pdf diff --git a/NEWS.md b/NEWS.md index 83f5d95..83f4f4e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,11 @@ +# jskm 0.5.6 + +* Update: Add median expression to the graph in `jskm` and `svyjskm`. + +# jskm 0.5.5 + +* Update: Align color of censoring marks with color of lines in `jskm`. + # jskm 0.5.5 * Update: Align color of censoring marks with color of lines in `jskm`. diff --git a/R/jskm.R b/R/jskm.R index fe47145..e5600dd 100644 --- a/R/jskm.R +++ b/R/jskm.R @@ -17,6 +17,7 @@ #' @param pval.testname logical: add '(Log-rank)' text to p-value. Default = F #' @param marks logical: should censoring marks be added? #' @param shape what shape should the censoring marks be, default is a vertical line +#' @param med should a median line be added to the plot? Default = F #' @param legend logical. should a legend be added to the plot? #' @param legendposition numeric. x, y position of the legend if plotted. Default=c(0.85,0.8) #' @param ci logical. Should confidence intervals be plotted. Default = FALSE @@ -89,7 +90,7 @@ jskm <- function(sfit, xlims = c(0, max(sfit$time)), ylims = c(0, 1), surv.scale = c("default", "percent"), - ystratalabs = names(sfit$strata), + ystratalabs = NULL, ystrataname = "Strata", timeby = signif(max(sfit$time) / 7, 1), main = "", @@ -99,6 +100,7 @@ jskm <- function(sfit, pval.testname = F, marks = TRUE, shape = 3, + med = FALSE, legend = TRUE, legendposition = c(0.85, 0.8), ci = FALSE, @@ -123,11 +125,11 @@ jskm <- function(sfit, ################################# # sorting the use of subsetting # ################################# - + n.risk <- n.censor <- surv <- strata <- lower <- upper <- NULL - + times <- seq(0, max(sfit$time), by = timeby) - if (!is.null(theme) && theme == "nejm") legendposition <- "right" + if (!is.null(theme) && theme == "nejm") legendposition <- c(1, 0.5) if (is.null(subs)) { if (length(levels(summary(sfit)$strata)) == 0) { subs1 <- 1 @@ -157,13 +159,13 @@ jskm <- function(sfit, subs2 <- which(regexpr(ssvar, summary(sfit, censored = T)$strata, perl = T) != -1) subs3 <- which(regexpr(ssvar, summary(sfit, times = times, extend = TRUE)$strata, perl = T) != -1) } - + if (!is.null(subs) | !is.null(sfit$states)) pval <- FALSE - + ################################## # data manipulation pre-plotting # ################################## - + if (is.null(ylabs)) { if (cumhaz | !is.null(sfit$states)) { ylabs <- "Cumulative incidence" @@ -171,29 +173,50 @@ jskm <- function(sfit, ylabs <- "Survival probability" } } - - - if (length(levels(summary(sfit)$strata)) == 0) { - # [subs1] - if (is.null(ystratalabs)) ystratalabs <- as.character(sub("group=*", "", "All")) + + if (!is.null(status.cmprsk)){ + if (length(levels(summary(sfit)$strata)) == 0) { + # [subs1] + if (is.null(ystratalabs)) ystratalabs <- as.character("All") + } else { + # [subs1] + if (is.null(ystratalabs)) ystratalabs <- as.character(names(sfit$strata)) + } } else { - # [subs1] - if (is.null(ystratalabs)) ystratalabs <- as.character(sub("group=*", "", names(sfit$strata))) + if (length(levels(summary(sfit)$strata)) == 0) { + # [subs1] + if (is.null(ystratalabs)) + ystratalabs <- as.character("All") + L <- summary(sfit)$table["0.95LCL"][[1]] + U <- summary(sfit)$table["0.95UCL"][[1]] + median_time <- summary(sfit)$table["median"][[1]] + ystratalabs2 <- paste0(ystratalabs, " (median : ", median_time, ", 0.95 CI : ",L, "-", U, ")") + } else { + # [subs1] + if (is.null(ystratalabs)) + ystratalabs <- as.character(names(sfit$strata)) + ystratalabs2 <- NULL + for (i in 1:length(levels(summary(sfit)$strata))) { + L <- summary(sfit)$table[,"0.95LCL"][[i]] + U <- summary(sfit)$table[,"0.95UCL"][[i]] + median_time <- summary(sfit)$table[,"median"][[i]] + ystratalabs2 <- c(ystratalabs2, paste0(ystratalabs[[i]], " (median : ", median_time, ", 0.95 CI : ",L, "-", U, ")")) + } + } } - if (is.null(ystrataname)) ystrataname <- "Strata" m <- max(nchar(ystratalabs)) times <- seq(0, max(sfit$time), by = timeby) - + if (length(levels(summary(sfit)$strata)) == 0) { Factor <- factor(rep("All", length(subs2))) } else { Factor <- factor(summary(sfit, censored = T)$strata[subs2], levels = names(sfit$strata)) } - + # Data to be used in the survival plot - - + + if (is.null(sfit$state)) { # no cmprsk df <- data.frame( time = sfit$time[subs2], @@ -221,9 +244,9 @@ jskm <- function(sfit, lower = sfit$lower[, col.cmprsk][subs2] ) } - + form <- sfit$call$formula - + if (!is.null(cut.landmark)) { if (is.null(data)) { data <- tryCatch(eval(sfit$call$data), error = function(e) e) @@ -231,7 +254,7 @@ jskm <- function(sfit, stop("Landmark analysis requires data object. please input 'data' option") } } - + var.time <- as.character(form[[2]][[2]]) var.event <- as.character(form[[2]][[3]]) if (length(var.event) > 1) { @@ -243,23 +266,23 @@ jskm <- function(sfit, data1 <- data data1[[var.event]][data1[[var.time]] >= cut.landmark] <- 0 data1[[var.time]][data1[[var.time]] >= cut.landmark] <- cut.landmark - + sfit1 <- survfit(as.formula(form), data1) sfit2 <- survfit(as.formula(form), data[data[[var.time]] >= cut.landmark, ]) - + if (is.null(sfit$states)) { if (length(levels(Factor)) == 1) { df2 <- merge(subset(df, time >= cut.landmark)[, c("time", "n.risk", "n.event", "n.censor", "strata")], - data.frame(time = sfit2$time, surv = sfit2$surv, strata = "All", upper = sfit2$upper, lower = sfit2$lower), - by = c("time", "strata") + data.frame(time = sfit2$time, surv = sfit2$surv, strata = "All", upper = sfit2$upper, lower = sfit2$lower), + by = c("time", "strata") ) } else { df2 <- merge(subset(df, time >= cut.landmark)[, c("time", "n.risk", "n.event", "n.censor", "strata")], - data.frame(time = sfit2$time, surv = sfit2$surv, strata = rep(names(sfit2$strata), sfit2$strata), upper = sfit2$upper, lower = sfit2$lower), - by = c("time", "strata") + data.frame(time = sfit2$time, surv = sfit2$surv, strata = rep(names(sfit2$strata), sfit2$strata), upper = sfit2$upper, lower = sfit2$lower), + by = c("time", "strata") ) } - + df11 <- rbind(subset(df, time < cut.landmark), df2[, names(df)]) df <- rbind(df11, data.frame(time = cut.landmark, n.risk = summary(sfit, times = cut.landmark)$n.risk[[1]], n.event = 0, n.censor = 0, surv = 1, strata = levels(df$strata), upper = 1, lower = 1)) } else { @@ -267,24 +290,24 @@ jskm <- function(sfit, status.cmprsk <- sfit$states[2] } col.cmprsk <- which(sfit$state == status.cmprsk) - + if (length(levels(Factor)) == 1) { df2 <- merge(subset(df, time >= cut.landmark)[, c("time", "n.risk", "n.event", "n.censor", "strata")], - data.frame(time = sfit2$time, surv = sfit2$pstate[, col.cmprsk], strata = "All", upper = sfit2$upper[, col.cmprsk], lower = sfit2$lower[, col.cmprsk]), - by = c("time", "strata") + data.frame(time = sfit2$time, surv = sfit2$pstate[, col.cmprsk], strata = "All", upper = sfit2$upper[, col.cmprsk], lower = sfit2$lower[, col.cmprsk]), + by = c("time", "strata") ) } else { df2 <- merge(subset(df, time >= cut.landmark)[, c("time", "n.risk", "n.event", "n.censor", "strata")], - data.frame(time = sfit2$time, surv = sfit2$pstate[, col.cmprsk], strata = rep(names(sfit2$strata), sfit2$strata), upper = sfit2$upper[, col.cmprsk], lower = sfit2$lower[, col.cmprsk]), - by = c("time", "strata") + data.frame(time = sfit2$time, surv = sfit2$pstate[, col.cmprsk], strata = rep(names(sfit2$strata), sfit2$strata), upper = sfit2$upper[, col.cmprsk], lower = sfit2$lower[, col.cmprsk]), + by = c("time", "strata") ) } df11 <- rbind(subset(df, time < cut.landmark), df2[, names(df)]) df <- rbind(df11, data.frame(time = cut.landmark, n.risk = summary(sfit, times = cut.landmark)$n.risk[[1]], n.event = 0, n.censor = 0, surv = 0, strata = levels(df$strata), upper = 0, lower = 0)) } } - - + + if (cumhaz & is.null(sfit$states)) { upper.new <- 1 - df$lower lower.new <- 1 - df$upper @@ -292,7 +315,7 @@ jskm <- function(sfit, df$lower <- lower.new df$upper <- upper.new } - + # Final changes to data for survival plot levels(df$strata) <- ystratalabs zeros <- data.frame( @@ -305,55 +328,55 @@ jskm <- function(sfit, zeros$lower <- 0 zeros$upper <- 0 } - + df <- rbind(zeros, df) d <- length(levels(df$strata)) - + ################################### # specifying axis parameteres etc # ################################### - + if (dashed == TRUE | all(linecols == "black")) { linetype <- c("solid", "dashed", "dotted", "dotdash", "longdash", "twodash", "1F", "F1", "4C88C488", "12345678") } else { linetype <- c("solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid") } - + # Scale transformation # :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: surv.scale <- match.arg(surv.scale) scale_labels <- ggplot2::waiver() if (surv.scale == "percent") scale_labels <- scales::percent - + p <- ggplot2::ggplot(df, aes(x = time, y = surv, colour = strata, linetype = strata)) + ggtitle(main) - - + + linecols2 <- linecols if (all(linecols == "black")) { linecols <- "Set1" p <- ggplot2::ggplot(df, aes(x = time, y = surv, linetype = strata)) + ggtitle(main) } - - + # Set up theme elements p <- p + theme_bw() + theme( axis.title.x = element_text(vjust = 0.7), panel.grid.minor = element_blank(), axis.line = element_line(linewidth = 0.5, colour = "black"), - legend.position = legendposition, + legend.position = "inside", + legend.position.inside = legendposition, legend.background = element_rect(fill = NULL), legend.key = element_rect(colour = NA), panel.border = element_blank(), - plot.margin = unit(c(0, 1, .5, ifelse(m < 10, 1.5, 2.5)), "lines"), + #plot.margin = unit(c(0, 1, .5, ifelse(m < 10, 1.5, 2.5)), "lines"), axis.line.x = element_line(linewidth = 0.5, linetype = "solid", colour = "black"), axis.line.y = element_line(linewidth = 0.5, linetype = "solid", colour = "black") ) + scale_x_continuous(xlabs, breaks = times, limits = xlims) + scale_y_continuous(ylabs, limits = ylims, labels = scale_labels) - + if (!is.null(theme) && theme == "jama") { p <- p + theme( panel.grid.major.x = element_blank() @@ -363,29 +386,40 @@ jskm <- function(sfit, panel.grid.major = element_blank() ) } - - + + # Removes the legend: if (legend == FALSE) { - p <- p + theme(legend.position = "none") + p <- p + guides(colour = "none", linetype = "none") } - + # Add lines too plot if (is.null(cut.landmark)) { - p <- p + geom_step(linewidth = linewidth) + - scale_linetype_manual(name = ystrataname, values = linetype) + if(med == T & is.null(status.cmprsk)){ + p <- p + geom_step(linewidth = linewidth) + + scale_linetype_manual(name = ystrataname, values = linetype, labels = ystratalabs2) + } else { + p <- p + geom_step(linewidth = linewidth) + + scale_linetype_manual(name = ystrataname, values = linetype, labels = ystratalabs) + } } else { - p <- p + - scale_linetype_manual(name = ystrataname, values = linetype) + - geom_step(data = subset(df, time >= cut.landmark), linewidth = linewidth) + geom_step(data = subset(df, time < cut.landmark), linewidth = linewidth) + if(med == T & is.null(status.cmprsk)){ + p <- p + + scale_linetype_manual(name = ystrataname, values = linetype, labels = ystratalabs2) + + geom_step(data = subset(df, time >= cut.landmark), linewidth = linewidth) + geom_step(data = subset(df, time < cut.landmark), linewidth = linewidth) + } else { + p <- p + + scale_linetype_manual(name = ystrataname, values = linetype, labels = ystratalabs) + + geom_step(data = subset(df, time >= cut.landmark), linewidth = linewidth) + geom_step(data = subset(df, time < cut.landmark), linewidth = linewidth) + } } - + brewer.palette <- c( "BrBG", "PiYG", "PRGn", "PuOr", "RdBu", "RdGy", "RdYlBu", "RdYlGn", "Spectral", "Accent", "Dark2", "Paired", "Pastel1", "Pastel2", "Set1", "Set2", "Set3", "Blues", "BuGn", "BuPu", "GnBu", "Greens", "Greys", "Oranges", "OrRd", "PuBu", "PuBuGn", "PuRd", "Purples", "RdPu", "Reds", "YlGn", "YlGnBu", "YlOrBr", "YlOrRd" ) - + if (!is.null(theme) && theme == "jama") { col.pal <- c("#00AFBB", "#E7B800", "#FC4E07") col.pal <- rep(col.pal, ceiling(length(ystratalabs) / 3)) @@ -395,29 +429,106 @@ jskm <- function(sfit, col.pal <- linecols col.pal <- rep(col.pal, ceiling(length(ystratalabs) / length(linecols))) } - - if (is.null(col.pal)) { - p <- p + scale_colour_brewer(name = ystrataname, palette = linecols) + + + if (is.null(cut.landmark)) { + if (med == T & is.null(status.cmprsk)){ + if (is.null(col.pal)) { + p <- p + scale_colour_brewer(name = ystrataname, palette = linecols, labels = ystratalabs2) + } else { + p <- p + scale_color_manual(name = ystrataname, values = col.pal, labels = ystratalabs2) + } + } else { + if (is.null(col.pal)) { + p <- p + scale_colour_brewer(name = ystrataname, palette = linecols, labels = ystratalabs) + } else { + p <- p + scale_color_manual(name = ystrataname, values = col.pal, labels = ystratalabs) + } + } } else { - p <- p + scale_color_manual(name = ystrataname, values = col.pal) + if (med == T & is.null(status.cmprsk)){ + if (is.null(col.pal)) { + p <- p + scale_colour_brewer(name = ystrataname, palette = linecols, labels = ystratalabs2) + } else { + p <- p + scale_color_manual(name = ystrataname, values = col.pal, labels = ystratalabs2) + } + } else { + if (is.null(col.pal)) { + p <- p + scale_colour_brewer(name = ystrataname, palette = linecols, labels = ystratalabs) + } else { + p <- p + scale_color_manual(name = ystrataname, values = col.pal, labels = ystratalabs) + } + } } - + # Add censoring marks to the line: if (marks == TRUE) { p <- p + geom_point(data = subset(df, n.censor >= 1), aes(x = time, y = surv, colour = strata), shape = shape) } - + + # Add median value + if (med == TRUE & is.null(cut.landmark) & is.null(status.cmprsk)){ + if (length(levels(summary(sfit)$strata)) == 0) { + median_time <- summary(sfit)$table["median"][[1]] + if(!is.na(median_time)) + p <- p + annotate("segment", x = xlims[1], xend = median_time, y = 0.5, yend = 0.5, linewidth = 0.3, linetype = "dashed") + + annotate("segment", x = median_time, xend = median_time, y = ylims[1], yend = 0.5, linewidth = 0.3, linetype = "dashed") + } else { + for (i in 1:length(levels(summary(sfit)$strata))){ + median_time <- summary(sfit)$table[,"median"][[i]] + if(!is.na(median_time)){ + p <- p + + annotate("segment", x = xlims[1], xend = median_time, y = 0.5, yend = 0.5, linewidth = 0.3, linetype = "dashed") + + annotate("segment", x = median_time, xend = median_time, y = ylims[1], yend = 0.5, linewidth = 0.3, linetype = "dashed") + } + } + } + } + + if (med == TRUE & !is.null(cut.landmark) & is.null(status.cmprsk)){ + if (length(levels(summary(sfit)$strata)) == 0) { + median_time <- summary(sfit1)$table[,"median"][[1]] + if(!is.na(median_time)) {p <- p + annotate("segment", x = xlims[1], xend = median_time, y = 0.5, yend = 0.5, linewidth = 0.3, linetype = "dashed") + annotate("segment", x = median_time, xend = median_time, y = ylims[1] , yend = 0.5, linewidth = 0.3, linetype = "dashed") } + median_time <- summary(sfit2)$table[,"median"][[1]] + if(!is.na(median_time)) {p <- p + annotate("segment", x = xlims[1], xend = median_time, y = 0.5, yend = 0.5, linewidth = 0.3, linetype = "dashed") + annotate("segment", x = median_time, xend = median_time, y = ylims[1], yend = 0.5, linewidth = 0.3, linetype = "dashed") } + + } else { + for (i in 1:length(levels(summary(sfit)$strata))){ + median_time <- summary(sfit1)$table[,"median"][[i]] + if(!is.na(median_time)) {p <- p + annotate("segment", x = xlims[1], xend = median_time, y = 0.5, yend = 0.5, linewidth = 0.3, linetype = "dashed") + annotate("segment", x = median_time, xend = median_time, y = ylims[1], yend = 0.5, linewidth = 0.3, linetype = "dashed") } + median_time <- summary(sfit2)$table[,"median"][[i]] + if(!is.na(median_time)) {p <- p + annotate("segment", x = xlims[1], xend = median_time, y = 0.5, yend = 0.5, linewidth = 0.3, linetype = "dashed") + annotate("segment", x = median_time, xend = median_time, y = ylims[1], yend = 0.5, linewidth = 0.3, linetype = "dashed") } + } + } + } + + + # Add 95% CI to plot if (ci == TRUE) { - if (all(linecols2 == "black")) { - p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper), alpha = 0.25, colour = NA) - } else if (is.null(col.pal)) { - 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) + if (med == FALSE | !is.null(status.cmprsk)) { + if (all(linecols2 == "black")) { + p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper), alpha = 0.25, colour = NA) + } else if (is.null(col.pal)) { + 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) + } else { + p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper, fill = strata), alpha = 0.25, colour = NA) + scale_fill_manual(name = ystrataname, values = col.pal) + } } else { - p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper, fill = strata), alpha = 0.25, colour = NA) + scale_fill_manual(name = ystrataname, values = col.pal) + if (all(linecols2 == "black")) { + p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper), alpha = 0.25, colour = NA) + } else if (is.null(col.pal)) { + 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, labels = ystratalabs2) + } else { + p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper, fill = strata), alpha = 0.25, colour = NA) + scale_fill_manual(name = ystrataname, values = col.pal, labels = ystratalabs2) + } } } - + + + + + if (!is.null(cut.landmark)) { p <- p + geom_vline(xintercept = cut.landmark, lty = 2) } @@ -452,8 +563,8 @@ jskm <- function(sfit, } } } - - + + ## Create a blank plot for place-holding blank.pic <- ggplot(df, aes(time, surv)) + geom_blank() + @@ -464,14 +575,14 @@ jskm <- function(sfit, axis.ticks = element_blank(), panel.grid.major = element_blank(), panel.border = element_blank() ) - + ##################### # p-value placement # ##################### a - + if (length(levels(summary(sfit)$strata)) == 0) pval <- F # if(!is.null(cut.landmark)) pval <- F - + if (pval == TRUE) { if (is.null(data)) { data <- tryCatch(eval(sfit$call$data), error = function(e) e) @@ -479,11 +590,11 @@ jskm <- function(sfit, stop("'pval' option requires data object. please input 'data' option") } } - + if (is.null(cut.landmark)) { sdiff <- survival::survdiff(as.formula(form), data = data) pvalue <- pchisq(sdiff$chisq, length(sdiff$n) - 1, lower.tail = FALSE) - + ## cluster option if (cluster.option == "cluster" & !is.null(cluster.var)) { form.old <- as.character(form) @@ -496,10 +607,10 @@ jskm <- function(sfit, sdiff <- survival::coxph(as.formula(form.new), data = data, model = T) pvalue <- summary(sdiff)$logtest["pvalue"] } - + pvaltxt <- ifelse(pvalue < 0.001, "p < 0.001", paste("p =", round(pvalue, 3))) if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)") - + # MOVE P-VALUE LEGEND HERE BELOW [set x and y] if (is.null(pval.coord)) { p <- p + annotate("text", x = (as.integer(max(sfit$time) / 5)), y = 0.1 + ylims[1], label = pvaltxt, size = pval.size) @@ -512,7 +623,7 @@ jskm <- function(sfit, pvalue <- sapply(list(sdiff1, sdiff2), function(x) { pchisq(x$chisq, length(x$n) - 1, lower.tail = FALSE) }) - + ## cluster option if (cluster.option == "cluster" & !is.null(cluster.var)) { form.old <- as.character(form) @@ -531,11 +642,11 @@ jskm <- function(sfit, summary(x)$logtest["pvalue"] }) } - + pvaltxt <- ifelse(pvalue < 0.001, "p < 0.001", paste("p =", round(pvalue, 3))) - + if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)") - + if (is.null(pval.coord)) { 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) } else { @@ -543,27 +654,27 @@ jskm <- function(sfit, } } } - + ################################################### # Create table graphic to include at-risk numbers # ################################################### - + n.risk <- NULL if (length(levels(summary(sfit)$strata)) == 0) { Factor <- factor(rep("All", length(subs3))) } else { Factor <- factor(summary(sfit, times = times, extend = TRUE)$strata[subs3]) } - + if (table == TRUE) { risk.data <- data.frame( strata = Factor, time = summary(sfit, times = times, extend = TRUE)$time[subs3], n.risk = summary(sfit, times = times, extend = TRUE)$n.risk[subs3] ) - + risk.data$strata <- factor(risk.data$strata, levels = rev(levels(risk.data$strata))) - + data.table <- ggplot(risk.data, aes(x = time, y = strata, label = format(n.risk, nsmall = 0))) + geom_text(size = 3.5) + theme_bw() + @@ -579,34 +690,33 @@ jskm <- function(sfit, axis.ticks = element_blank(), axis.text.y = element_text(face = "bold", hjust = 1) ) data.table <- data.table + - theme(legend.position = "none") + xlab(NULL) + ylab(NULL) - - + guides(colour = "none", linetype = "none")+ xlab(NULL) + ylab(NULL) + + # ADJUST POSITION OF TABLE FOR AT RISK data.table <- data.table + theme(plot.margin = unit(c(-1.5, 1, 0.1, ifelse(m < 10, 3.1, 4.3) - 0.38 * m), "lines")) } - - + + ####################### # Plotting the graphs # ####################### - + if (!is.null(theme) && theme == "nejm") { - p2 <- p1 + coord_cartesian(ylim = nejm.infigure.ylim) + theme( - legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), - axis.text = element_text(size = 10 * nejm.infigure.ratiow) - ) + p2 <- p1 + coord_cartesian(ylim = nejm.infigure.ylim) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), + axis.text = element_text(size = 10 * nejm.infigure.ratiow) + ) + guides(colour = "none", linetype = "none") p <- p + patchwork::inset_element(p2, 1 - nejm.infigure.ratiow, 1 - nejm.infigure.ratioh, 1, 1, align_to = "panel") } - + if (table == TRUE) { ggpubr::ggarrange(p, blank.pic, data.table, - nrow = 3, - # align = "v", - heights = c(2, .1, .25) + nrow = 3, + # align = "v", + heights = c(2, .1, .25) ) } else { p } -} +} \ No newline at end of file diff --git a/R/svyjskm.R b/R/svyjskm.R index eee69dc..759140a 100644 --- a/R/svyjskm.R +++ b/R/svyjskm.R @@ -14,9 +14,10 @@ #' @param pval.size numeric value specifying the p-value text size. Default is 5. #' @param pval.coord numeric vector, of length 2, specifying the x and y coordinates of the p-value. Default values are NULL #' @param pval.testname logical: add '(Log-rank)' text to p-value. Default = F -#' @param legend logical. should a legend be added to the plot? Default: TRUE +#' @param med should a median line be added to the plot? Default = F +#' @param legend logical. should a legend be added to the plot? +#' @param legendposition numeric. x, y position of the legend if plotted. Default=c(0.85,0.8) #' @param ci logical. Should confidence intervals be plotted. Default = NULL -#' @param legendposition numeric. x, y position of the legend if plotted. Default: c(0.85, 0.8) #' @param linecols Character or Character vector. Colour brewer pallettes too colour lines. Default ="Set1", "black" for black with dashed line, character vector for the customization of line colors. #' @param dashed logical. Should a variety of linetypes be used to identify lines. Default: FALSE #' @param cumhaz Show cumulaive incidence function, Default: F @@ -36,13 +37,13 @@ #' @return plot #' @details DETAILS #' @examples -#' library(survey) -#' 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) ~ sex, design = dpbc) +# library(survey) +# 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) ~ sex, design = dpbc) #' svyjskm(s1) #' @rdname svyjskm #' @import ggplot2 @@ -68,6 +69,7 @@ svyjskm <- function(sfit, pval.size = 5, pval.coord = c(NULL, NULL), pval.testname = F, + med = FALSE, legend = TRUE, legendposition = c(0.85, 0.8), ci = NULL, @@ -87,8 +89,8 @@ svyjskm <- function(sfit, nejm.infigure.ylim = c(0, 1), ...) { surv <- strata <- lower <- upper <- NULL - - if (!is.null(theme) && theme == "nejm") legendposition <- "right" + + if (!is.null(theme) && theme == "nejm") legendposition <- c(1, 0.5) if (is.null(timeby)) { if (inherits(sfit, "svykmlist")) { timeby <- signif(max(sapply(sfit, function(x) { @@ -98,7 +100,7 @@ svyjskm <- function(sfit, timeby <- signif(max(sfit$time) / 7, 1) } } - + if (is.null(ci)) { if (inherits(sfit, "svykmlist")) { ci <- "varlog" %in% names(sfit[[1]]) @@ -106,8 +108,8 @@ svyjskm <- function(sfit, ci <- "varlog" %in% names(sfit) } } - - + + if (ci & !is.null(cut.landmark)) { if (is.null(design)) { design <- tryCatch(get(as.character(attr(sfit, "call")$design)), error = function(e) e) @@ -123,56 +125,92 @@ svyjskm <- function(sfit, "warning" %in% class(tryCatch(as.numeric(x), warning = function(w) w)) })] } + design1 <- design design1$variables[[var.event]][design1$variables[[var.time]] >= cut.landmark] <- 0 design1$variables[[var.time]][design1$variables[[var.time]] >= cut.landmark] <- cut.landmark - + sfit2 <- survey::svykm(formula(sfit), design = subset(design, get(var.time) >= cut.landmark), se = T) } - - - - + + + + if (inherits(sfit, "svykmlist")) { if (is.null(ystrataname)) ystrataname <- as.character(formula(sfit)[[3]]) - + if (ci) { if ("varlog" %in% names(sfit[[1]])) { - 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)))) - })) if (!is.null(cut.landmark)) { + 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)))) + })) 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)))) })) df <- rbind(df[df$time < cut.landmark, ], data.frame("strata" = unique(df$strata), "time" = cut.landmark, "surv" = 1, "lower" = 1, "upper" = 1), df2) + } else { + if (med == TRUE){ + df <- do.call(rbind, lapply(names(sfit), function(x) { + df2 <- 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)))) + valid_indices <- which(df2$surv <= 0.5) + closest_index <- valid_indices[which.min(abs(df2$surv[valid_indices] - 0.5))] + df2$med <- df2[closest_index,"time"] + return(df2) + })) + } else { + 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)))) + })) + } } } else { stop("No CI information in svykmlist object. please run svykm with se = T option.") - } - } else { - df <- do.call(rbind, lapply(names(sfit), function(x) { - data.frame("strata" = x, "time" = sfit[[x]]$time, "surv" = sfit[[x]]$surv) - })) + } + }else { if (!is.null(cut.landmark)) { + df <- do.call(rbind, lapply(names(sfit), function(x) { + data.frame("strata" = x, "time" = sfit[[x]]$time, "surv" = sfit[[x]]$surv) + })) for (v in unique(df$strata)) { if (nrow(subset(df, time == cut.landmark & strata == v)) == 0) { df <- rbind(df, data.frame(strata = v, time = cut.landmark, surv = 1)) } else { df[df$time == cut.landmark & df$strata == v, "surv"] <- 1 } - + 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"]) } + } else { + if (med == TRUE){ + df <- do.call(rbind, lapply(names(sfit), function(x) { + df2 <- data.frame("strata" = x, "time" = sfit[[x]]$time, "surv" = sfit[[x]]$surv) + valid_indices <- which(df2$surv <= 0.5) + closest_index <- valid_indices[which.min(abs(df2$surv[valid_indices] - 0.5))] + df2$med <- df2[closest_index,"time"] + return(df2) + }))} else { + df <- do.call(rbind, lapply(names(sfit), function(x) { + data.frame("strata" = x, "time" = sfit[[x]]$time, "surv" = sfit[[x]]$surv) + })) + } } } - + df$strata <- factor(df$strata, levels = names(sfit)) times <- seq(0, max(sapply(sfit, function(x) { max(x$time) })), by = timeby) if (is.null(ystratalabs)) { + df3 <- df[-c(1,2),] ystratalabs <- levels(df$strata) + if (med == TRUE){ + ystratalabs2 <-NULL + for (i in 1:length(names(sfit))) { + median_time <- df3[df3$strata == names(sfit)[[i]], "med"] %>% unique + ystratalabs2 <- c(ystratalabs2, paste0(ystratalabs[[i]], " (median : ", median_time, ")")) + } + } } if (is.null(xlims)) { xlims <- c(0, max(sapply(sfit, function(x) { @@ -181,45 +219,64 @@ svyjskm <- function(sfit, } } else if (inherits(sfit, "svykm")) { if (is.null(ystrataname)) ystrataname <- "Strata" - + if (ci) { if ("varlog" %in% names(sfit)) { - 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)))) if (!is.null(cut.landmark)) { + 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)))) 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)))) df <- rbind(df[df$time < cut.landmark, ], data.frame("strata" = "All", "time" = cut.landmark, "surv" = 1, "lower" = 1, "upper" = 1), df2) + } else { + if (med == T){ + 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)))) + valid_indices <- which(df$surv <= 0.5) + closest_index <- valid_indices[which.min(abs(df$surv[valid_indices] - 0.5))] + df$med <- df[closest_index,"time"] + } else { + 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)))) + } } } else { stop("No CI information in svykm object. please run svykm with se = T option.") } } else { - df <- data.frame("strata" = "All", "time" = sfit$time, "surv" = sfit$surv) if (!is.null(cut.landmark)) { + df <- data.frame("strata" = "All", "time" = sfit$time, "surv" = sfit$surv) if (nrow(subset(df, time == cut.landmark)) == 0) { df <- rbind(df, data.frame(strata = "All", time = cut.landmark, surv = 1)) } else { df[df$time == cut.landmark, "surv"] <- 1 } - + df[df$time > cut.landmark, "surv"] <- df[df$time > cut.landmark, "surv"] / min(df[df$time < cut.landmark, "surv"]) + } else{ + if (med == T){ + df <- data.frame("strata" = "All", "time" = sfit$time, "surv" = sfit$surv) + valid_indices <- which(df$surv <= 0.5) + closest_index <- valid_indices[which.min(abs(df$surv[valid_indices] - 0.5))] + df$med <- df[closest_index,"time"] + } else { + df <- data.frame("strata" = "All", "time" = sfit$time, "surv" = sfit$surv) + } } } - + times <- seq(0, max(sfit$time), by = timeby) if (is.null(ystratalabs)) { ystratalabs <- "All" + ystratalabs2 <- paste0(ystratalabs, " (median : ", unique(df3$med), ")") } if (is.null(xlims)) { xlims <- c(0, max(sfit$time)) } } - + m <- max(nchar(ystratalabs)) - - - - - + + + + + if (cumhaz) { df$surv <- 1 - df$surv if (ci) { @@ -229,15 +286,26 @@ svyjskm <- function(sfit, df$upper <- upper.new } } - + # Final changes to data for survival plot levels(df$strata) <- ystratalabs - zeros <- data.frame("strata" = factor(ystratalabs, levels = levels(df$strata)), "time" = 0, "surv" = 1) + zeros <- if (med == T){ + data.frame("strata" = factor(ystratalabs, levels = levels(df$strata)), "time" = 0, "surv" = 1, "med" = 0.5) + } else { + data.frame("strata" = factor(ystratalabs, levels = levels(df$strata)), "time" = 0, "surv" = 1) + } + if (ci) { - zeros$upper <- 1 - zeros$lower <- 1 + if (med == T){ + zeros$med <- NULL + zeros$upper <- 1 + zeros$lower <- 1 + zeros$med <- 0.5 + } else { + zeros$upper <- 1 + zeros$lower <- 1 + } } - if (cumhaz) { zeros$surv <- 0 if (ci) { @@ -245,27 +313,27 @@ svyjskm <- function(sfit, zeros$upper <- 0 } } - + df <- rbind(zeros, df) d <- length(levels(df$strata)) - + ################################### # specifying axis parameteres etc # ################################### - + if (dashed == TRUE | all(linecols == "black")) { linetype <- c("solid", "dashed", "dotted", "dotdash", "longdash", "twodash", "1F", "F1", "4C88C488", "12345678") } else { linetype <- c("solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid", "solid") } - + # Scale transformation # :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: surv.scale <- match.arg(surv.scale) scale_labels <- ggplot2::waiver() if (surv.scale == "percent") scale_labels <- scales::percent - - + + p <- ggplot2::ggplot(df, aes(x = time, y = surv, colour = strata, linetype = strata)) + ggtitle(main) linecols2 <- linecols @@ -274,14 +342,15 @@ svyjskm <- function(sfit, p <- ggplot2::ggplot(df, aes(x = time, y = surv, linetype = strata)) + ggtitle(main) } - + # Set up theme elements p <- p + theme_bw() + theme( axis.title.x = element_text(vjust = 0.7), panel.grid.minor = element_blank(), axis.line = element_line(linewidth = 0.5, colour = "black"), - legend.position = legendposition, + legend.position = "inside", + legend.position.inside = legendposition, legend.background = element_rect(fill = NULL), legend.key = element_rect(colour = NA), panel.border = element_blank(), @@ -291,7 +360,7 @@ svyjskm <- function(sfit, ) + scale_x_continuous(xlabs, breaks = times, limits = xlims) + scale_y_continuous(ylabs, limits = ylims, labels = scale_labels) - + if (!is.null(theme) && theme == "jama") { p <- p + theme( panel.grid.major.x = element_blank() @@ -301,28 +370,51 @@ svyjskm <- function(sfit, panel.grid.major = element_blank() ) } - - + + # Removes the legend: if (legend == FALSE) { - p <- p + theme(legend.position = "none") + p <- p + guides(colour = "none", linetype = "none") } - + # Add lines too plot if (is.null(cut.landmark)) { - p <- p + geom_step(linewidth = linewidth) + - scale_linetype_manual(name = ystrataname, values = linetype) + if (med == T){ + p <- p + geom_step(linewidth = linewidth) + + scale_linetype_manual(name = ystrataname, values = linetype, labels = ystratalabs2) + } else { + p <- p + geom_step(linewidth = linewidth) + + scale_linetype_manual(name = ystrataname, values = linetype, labels = ystratalabs) + } } else { p <- p + geom_step(data = subset(df, time < cut.landmark), linewidth = linewidth) + geom_step(data = subset(df, time >= cut.landmark), linewidth = linewidth) + - scale_linetype_manual(name = ystrataname, values = linetype) + scale_linetype_manual(name = ystrataname, values = linetype, labels = ystratalabs) } - + + + # Add median value + if (med == TRUE & is.null(cut.landmark)){ + df3 <- df[-c(1,2),] + if (inherits(sfit, "svykm")) { + median_time <- df3$med %>% unique + if(!is.na(median_time)) + p <- p + annotate("segment", x = xlims[1], xend = median_time, y = 0.5, yend = 0.5, linewidth = 0.3, linetype = "dashed") + annotate("segment", x = median_time, xend = median_time, y = ylims[1], yend = 0.5, linewidth = 0.3, linetype = "dashed") + } else { + for (i in 1:length(names(sfit))){ + median_time <- df3[df3$strata == names(sfit)[[i]], "med"] %>% unique + if(!is.na(median_time)){ + p <- p + annotate("segment", x = xlims[1], xend = median_time, y = 0.5, yend = 0.5, linewidth = 0.3, linetype = "dashed") + annotate("segment", x = median_time, xend = median_time, y = ylims[1], yend = 0.5, linewidth = 0.3, linetype = "dashed") + } + } + } + } + brewer.palette <- c( "BrBG", "PiYG", "PRGn", "PuOr", "RdBu", "RdGy", "RdYlBu", "RdYlGn", "Spectral", "Accent", "Dark2", "Paired", "Pastel1", "Pastel2", "Set1", "Set2", "Set3", "Blues", "BuGn", "BuPu", "GnBu", "Greens", "Greys", "Oranges", "OrRd", "PuBu", "PuBuGn", "PuRd", "Purples", "RdPu", "Reds", "YlGn", "YlGnBu", "YlOrBr", "YlOrRd" ) - + if (!is.null(theme) && theme == "jama") { col.pal <- c("#00AFBB", "#E7B800", "#FC4E07") col.pal <- rep(col.pal, ceiling(length(ystratalabs) / 3)) @@ -332,28 +424,67 @@ svyjskm <- function(sfit, col.pal <- linecols col.pal <- rep(col.pal, ceiling(length(ystratalabs) / length(linecols))) } - - if (is.null(col.pal)) { - p <- p + scale_colour_brewer(name = ystrataname, palette = linecols) + # + # if (is.null(cut.landmark)) { + # if (is.null(col.pal)) { + # p <- p + scale_colour_brewer(name = ystrataname, palette = linecols, labels = ystratalabs2) + # } else { + # p <- p + scale_color_manual(name = ystrataname, values = col.pal, labels = ystratalabs2) + # }} else {if (is.null(col.pal)) { + # p <- p + scale_colour_brewer(name = ystrataname, palette = linecols) + # } else { + # p <- p + scale_color_manual(name = ystrataname, values = col.pal) + # } } + + + if (is.null(cut.landmark)) { + if (med == T){ + if (is.null(col.pal)) { + p <- p + scale_colour_brewer(name = ystrataname, palette = linecols, labels = ystratalabs2) + } else { + p <- p + scale_color_manual(name = ystrataname, values = col.pal, labels = ystratalabs2) + } + } else { + if (is.null(col.pal)) { + p <- p + scale_colour_brewer(name = ystrataname, palette = linecols, labels = ystratalabs) + } else { + p <- p + scale_color_manual(name = ystrataname, values = col.pal, labels = ystratalabs) + } + } } else { - p <- p + scale_color_manual(name = ystrataname, values = col.pal) + if (is.null(col.pal)) { + p <- p + scale_colour_brewer(name = ystrataname, palette = linecols, labels = ystratalabs) + } else { + p <- p + scale_color_manual(name = ystrataname, values = col.pal, labels = ystratalabs) + } } - # Add 95% CI to plot + + if (ci == TRUE) { - if (all(linecols2 == "black")) { - p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper), alpha = 0.25, colour = NA) - } else if (is.null(col.pal)) { - 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) + if (med == TRUE & is.null(cut.landmark)) { + if (all(linecols2 == "black")) { + p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper), alpha = 0.25, colour = NA) + } else if (is.null(col.pal)) { + 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, labels = ystratalabs2) + } else { + p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper, fill = strata), alpha = 0.25, colour = NA) + scale_fill_manual(name = ystrataname, values = col.pal, labels = ystratalabs2) + } } else { - p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper, fill = strata), alpha = 0.25, colour = NA) + scale_fill_manual(name = ystrataname, values = col.pal) + if (all(linecols2 == "black")) { + p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper), alpha = 0.25, colour = NA) + } else if (is.null(col.pal)) { + 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, labels = ystratalabs) + } else { + p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper, fill = strata), alpha = 0.25, colour = NA) + scale_fill_manual(name = ystrataname, values = col.pal, labels = ystratalabs) + } } } - + if (!is.null(cut.landmark)) { p <- p + geom_vline(xintercept = cut.landmark, lty = 2) } - + p1 <- p ## p-value if (inherits(sfit, "svykm")) pval <- FALSE @@ -384,14 +515,14 @@ svyjskm <- function(sfit, stop("'pval' option requires design object. please input 'design' option") } } - + if (is.null(cut.landmark)) { sdiff <- survey::svylogrank(formula(sfit), design = design) pvalue <- sdiff[[2]][2] - + pvaltxt <- ifelse(pvalue < 0.001, "p < 0.001", paste("p =", round(pvalue, 3))) if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)") - + # MOVE P-VALUE LEGEND HERE BELOW [set x and y] if (is.null(pval.coord)) { p <- p + annotate("text", x = as.integer(max(sapply(sfit, function(x) { @@ -418,16 +549,16 @@ svyjskm <- function(sfit, design1 <- design design1$variables[[var.event]][design1$variables[[var.time]] >= cut.landmark] <- 0 design1$variables[[var.time]][design1$variables[[var.time]] >= cut.landmark] <- cut.landmark - + sdiff1 <- survey::svylogrank(formula(sfit), design = design1) sdiff2 <- survey::svylogrank(formula(sfit), design = subset(design, get(var.time) >= cut.landmark)) pvalue <- sapply(list(sdiff1, sdiff2), function(x) { x[[2]][2] }) - + pvaltxt <- ifelse(pvalue < 0.001, "p < 0.001", paste("p =", round(pvalue, 3))) if (pval.testname) pvaltxt <- paste0(pvaltxt, " (Log-rank)") - + if (is.null(pval.coord)) { p <- p + annotate("text", x = c(as.integer(max(sapply(sfit, function(x) { max(x$time) / 10 @@ -439,11 +570,11 @@ svyjskm <- function(sfit, } } } - - - - - + + + + + ## Create a blank plot for place-holding blank.pic <- ggplot(df, aes(time, surv)) + geom_blank() + @@ -454,11 +585,11 @@ svyjskm <- function(sfit, axis.ticks = element_blank(), panel.grid.major = element_blank(), panel.border = element_blank() ) - + ################################################### # Create table graphic to include at-risk numbers # ################################################### - + n.risk <- NULL if (table == TRUE) { if (is.null(design)) { @@ -466,9 +597,9 @@ svyjskm <- function(sfit, } else { sfit2 <- survival::survfit(formula(sfit), data = design$variables) } - + # times <- seq(0, max(sfit2$time), by = timeby) - + if (is.null(subs)) { if (length(levels(summary(sfit2)$strata)) == 0) { subs1 <- 1 @@ -498,27 +629,27 @@ svyjskm <- function(sfit, subs2 <- which(regexpr(ssvar, summary(sfit2, censored = T)$strata, perl = T) != -1) subs3 <- which(regexpr(ssvar, summary(sfit2, times = times, extend = TRUE)$strata, perl = T) != -1) } - + if (!is.null(subs)) pval <- FALSE - - - + + + if (length(levels(summary(sfit2)$strata)) == 0) { Factor <- factor(rep("All", length(subs3))) } else { Factor <- factor(summary(sfit2, times = times, extend = TRUE)$strata[subs3]) } - - + + risk.data <- data.frame( strata = Factor, time = summary(sfit2, times = times, extend = TRUE)$time[subs3], n.risk = summary(sfit2, times = times, extend = TRUE)$n.risk[subs3] ) - - + + risk.data$strata <- factor(risk.data$strata, levels = rev(levels(risk.data$strata))) - + data.table <- ggplot(risk.data, aes(x = time, y = strata, label = format(n.risk, nsmall = 0))) + geom_text(size = 3.5) + theme_bw() + @@ -534,32 +665,31 @@ svyjskm <- function(sfit, axis.ticks = element_blank(), axis.text.y = element_text(face = "bold", hjust = 1) ) data.table <- data.table + - theme(legend.position = "none") + xlab(NULL) + ylab(NULL) - - + guides(colour = "none", linetype = "none")+ xlab(NULL) + ylab(NULL) + + # ADJUST POSITION OF TABLE FOR AT RISK data.table <- data.table + theme(plot.margin = unit(c(-1.5, 1, 0.1, ifelse(m < 10, 3.1, 4.3) - 0.38 * m), "lines")) } - + ####################### # Plotting the graphs # ####################### if (!is.null(theme) && theme == "nejm") { - p2 <- p1 + coord_cartesian(ylim = nejm.infigure.ylim) + theme( - legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), - axis.text = element_text(size = 10 * nejm.infigure.ratiow) - ) + p2 <- p1 + coord_cartesian(ylim = nejm.infigure.ylim) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), + axis.text = element_text(size = 10 * nejm.infigure.ratiow) + ) + guides(colour = "none", linetype = "none") p <- p + patchwork::inset_element(p2, 1 - nejm.infigure.ratiow, 1 - nejm.infigure.ratioh, 1, 1, align_to = "panel") } - + if (table == TRUE) { ggpubr::ggarrange(p, blank.pic, data.table, - nrow = 3, - # align = "v", - heights = c(2, .1, .25) + nrow = 3, + # align = "v", + heights = c(2, .1, .25) ) } else { p } -} +} \ No newline at end of file diff --git a/README.Rmd b/README.Rmd index a9f2b14..4bf1d46 100644 --- a/README.Rmd +++ b/README.Rmd @@ -49,7 +49,7 @@ 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, + table = T, pval = T, med = 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 ) diff --git a/README.md b/README.md index 18ea6d7..a355197 100644 --- a/README.md +++ b/README.md @@ -5,6 +5,7 @@ Kaplan-Meier Plot with ‘ggplot2’: ‘survfit’ and ‘svykm’ objects from ‘survival’ and ‘survey’ packages. + [![AppVeyor build status](https://ci.appveyor.com/api/projects/status/github/jinseob2kim/jskm?branch=master&svg=true)](https://ci.appveyor.com/project/jinseob2kim/jskm) [![Github @@ -20,7 +21,6 @@ stars](https://img.shields.io/github/stars/jinseob2kim/jskm.svg)](https://github license](https://img.shields.io/github/license/jinseob2kim/jskm.svg)](https://github.com/jinseob2kim/jskm/blob/master/LICENSE) - ## Install ``` r @@ -51,7 +51,7 @@ jskm(fit) ``` r jskm(fit, - table = T, pval = T, label.nrisk = "No. at risk", size.label.nrisk = 8, + table = T, pval = T, med = 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 ) @@ -107,7 +107,7 @@ jskm(fit2, mark = F, surv.scale = "percent", table = T, status.cmprsk = "1", sho #### JAMA ``` r -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)) +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)) ``` ![](man/figures/README-unnamed-chunk-5-1.png) @@ -115,7 +115,7 @@ jskm(fit, theme='jama', cumhaz = T, table=T, mark = F, ylab = "Cumulative incide #### NEJM ``` r -jskm(fit, theme='nejm', nejm.infigure.ratiow = 0.7, 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)) +jskm(fit, theme = "nejm", nejm.infigure.ratiow = 0.7, 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)) ``` ![](man/figures/README-unnamed-chunk-6-1.png) @@ -178,7 +178,7 @@ svyjskm(s3, ci = F, surv.scale = "percent", pval = T, table = T, cut.landmark = #### JAMA ``` r -svyjskm(s2, theme='jama', pval = T, table = T, design = dpbc) +svyjskm(s2, theme = "jama", pval = T, table = T, design = dpbc) ``` ![](man/figures/README-unnamed-chunk-9-1.png) @@ -186,7 +186,7 @@ svyjskm(s2, theme='jama', pval = T, table = T, design = dpbc) #### NEJM ``` r -svyjskm(s2, theme='nejm', nejm.infigure.ratiow = 0.45, nejm.infigure.ratioh = 0.4, nejm.infigure.ylim = c(0.2,1), pval = T, table = T, design = dpbc) +svyjskm(s2, theme = "nejm", nejm.infigure.ratiow = 0.45, nejm.infigure.ratioh = 0.4, nejm.infigure.ylim = c(0.2, 1), pval = T, table = T, design = dpbc) ``` ![](man/figures/README-unnamed-chunk-10-1.png) diff --git a/Rplot.pdf b/Rplot.pdf new file mode 100644 index 0000000000000000000000000000000000000000..0fffcd0d615e7751a5fac6a8424fca4513856582 GIT binary patch literal 25891 zcmZ^}cT|&I&@cL0@g*WEN-ru3O6V=rAfTX#^b)!OLhqf>L=1{Ruu(!4P-+Mr=}5qY z5_<232mwMdgaDy=c)xYdxoh3K?qAP-X3xyt^P8F7?-@MO5|a{_Wxf}^7)^>Eh#s-_ zW0txldCS}0jafy7Swhnf1`Tla_Id;jfZe+PNJ&~!T2@L@3M3^jEhQ)YfLY?Xx8whB zor!CJ$ElRjEp0!jkBh5=|1JH%0H45sQ!xWSZ^u9fnBRZ3vf?u0Qm3^4Y0LXR(gS8N zm>K5f_+Kxi|8H+@NtlHCz-~!sLIa>4-p;oq44}?1|64Mr3`V#99|6$+N8mrc|G52N zUIG7`T|(2_E8z6k|CYl4h)Wp3{JjJH98UGg{=fA7$MApP(1m$92e{mlmI6u2+>+38 z^$0jMK|;&p6rD#f2X9B%e;xM^I0em-SmyoMGLa-dN*|NU z)lO)vIe&LH7Pi%*PToZAnMW}gV*;%V3G%G-$sTLlfojwdZA)J5WRu+X{_pq5lhv7p zqYZpS9CvcR)cK8rw}`x7i;mm7_2gA@V3N#tyqwCUXJu~VHdj`(@Uq|L9%_^+SN z2Jm&V+r@A2{KeMmtM#0qZ96e$dim9@8?Dl-3h#IxGt9m%ny;OFA<6%Q7v62l9i7gq znq7-&ww~ijE+1CUon7U;!*b&=YdTeMslur23kf+_6m_qmQLKxvX;?~HWwo9JsJ&cr ztcUl$o4cQk)lq$w-ei0vawg!Y4WwgdxVrBhRVZq%_h2Hq#CA?}U_Y_P+(^U@vrydA zrh{LxDEf%3ZXJ@kPVbp0pWhm7Yqn(w1$*Ucxdkg}Js$@yz;4WT^(K4`NO;AADF65n zf4|=F!0(#LwvH*rOD+3WRNuK<-69ox;x#2wu$N?K-Z#;Jw2tMwV?N?=8k zL8%zA(IuUo`+px2eKza!kj0V7T=ppSjNmlp8E5&J35B(1}Qj!$EK( zu~zr?>WOcknqSP`;?djw@ydAG4}ZL_l+E6smI?dyt&zBWGZEG-YAqGr`{S~2uTos^ z4;81<pAn2Z$h;=iwE142WOjer{ zG6Ep5R|9wNv9*F9Se4%+)J`;9#R=!v2AWM+ij=T+Sli! z=8zF=kC)7%T3g=@hOycQp=4VRf-v(ROZ-*)g z*-k}^&F-tf64T3O87fN&G@Ax)a#fW_zp@4H5L=h6B(&@RZF3+@;5$VrLx^)^cimZ_ z@aFd_EtMeaL>Smgfff;{JB$ATY4B-5Rq{{y$cRHbk}0{b=snm(N^VRNVh)B(SSm}i6I34w#M&Q4zE@W&@klhVtd zNI!kh-_69s?F87`dU-IJ`d3n@bx?xvbCMhov?+X-d0u{&YfDimvc(RR;`tU_{v0%a z>ZFF%xU%(z1y!M-(0TT8PwD;rb|kf#H(Z4is125QKB58?220@Mb5q3xZ3lOLHu)__ znh2fL5nr8;uTWtHm?N1&rS`>05i{p7+N$Nu5jSicIj*BYtLCWJHZ*?=jc_~0!KgVsQ)d$c1us9n9be{Dla?4v^t)Gai z^lFHa$`P~nZ95IH!%Cly@LQs^b}ZkXm5_L-fCq$p4wbeO?w!fI`Tfsnpc3e@70nXz zXzE4OZQ7wtsLxWQ2ZVHH`^9b} z@nuN1O4@BN1$)(dtSJe277ARazKKP+vvjHb<6*6jmWl~K+JRK>r~XlUpHo43{O(&7 zh_o#S<}H~LV=g~ik@G%UU>AL8$~6gB+nlY;PKg$ylqo-J)+nFwCH;=5k(`Lf;(tt6 zVmJ^KQT}b(B z2h{-Gbs!VkTw^587lzsEsGL_2O5t;8RT@9Hts<6Lk$p3TTH~YDe5(s54hCChBUNVqS<1S{|A(NywI`!( zk=9a`u;<(_?PFWt6%SoX2!C{0q->@V2Ux7kb~#Brw}}TZwGm&1S5=^$*)eamJ2|DY zw`Q+;#~XsY%0;(Ioy~3z1Id5-t8yZo2*KW>ly=Z2F)?FvZ=t>rw8`l$!P^n5JRgje zUI1@b@{dsebn)-eaC6YMySy~zUib+r1|vMbuTaT>Y|IW_{z}aR?YegGQ_Y}U5(Alm z;Lznf$_A)^DeGAnsf0oV?H*?(B{Bf+_LO9Ya1w_44y5PVF$#QSF;ECHySD=DjoeFg z4|R@XQt|&05(Hv&MRsR+!du_opFcziT+Wf66$AXuMqyN%!c_7i8gbR zh>udGlAcQ!Rii8B%Q?mm>79A-AaCA3Dqz1(dN^@y-Xf<^8#6(Z+Z^n8d>5J109P4E zO4vK*!SD26FAsBwn^5&1KK>=cVQALox96R(N0r%|q1RV-&e);0;yGG*w_-ITNO_0F zr?Aap_WX&+J!*13U{*pl%iE3-GH=J2OjE1A|7lK{To1)nFG^LuKI!jk@lBx*0?qeu zN_MRk7m)&rar5JCm2H^=JxT&C6<(KZA%Q&-vl0{T&?)_078zUdaJ0Fd!joS+2-EUI zti?qf@Xa=^aXyOpj&j_^GI97pPU=rjLyM1r>Q5|Q-fE$6byB~Fp;FziErEvjbC{`i z$UJLSYOB1+ZVTbNBQfjU7Df8XljPhl;sGb#ml)^RoS?N-pFcU`S?DL@?DYYy;g2Hk z2-U(pvOBZuI0unaq0(Reg+0nF5&76+H;$RHzO?2sTS&OVU%)j@{&IPnq!%+-b(z#RdMU^n29QtUARcW&) z@vNtT)7^}?wN2>VrjJ_LNr!|JNy?utAEq)m2Ue&?pSD%}?Q>PFdgtPIY7tIHww z!rVTkur!ITMfejC=r>z*tLgS)3|v_2q)9~9tF!1>{ae{rF2)El$SByrDpRN$;979} z#X)xV{+(x$M|)Bl?8BEs_5G5)KfmWwJ2?!&dVl`F#Svy<(~WkJ{L*J4=@4XkWwz+{ z%!DZ+EMG6)>O+pOhDLCO$d9QqBS=ygJ zfi7WNU4xbUo(s}nqn-Ggss{$}pYpovmFLhBx`6n7- z5Kulw)7(kyde(kVtW4RzOypLghI{2#3@HU_lPF?i*I)z}kkhX6%iozlMWSOmKo?!> zS)W0F_KgBCcRuyi;d-)a8PIS7s~j7kk|*1a{cPLa2^^qdomvzi*%K_Kvf9E($O}0@a^-81$;QWwSq73H zt~OfqZu0gi6$~Hqbxgmnu(ehgQcC{sGjfnFi>o9dDbe&_Sj703uftQxa(6;~#Z0X0|J(V7ycE7NN zHvuV?f35=>0xFc-5 z*C?raMjEYVx_}32TSomCu~3$g=r+FsdJR^o68hP;_$Jw(tl|Dm_XyT;GndfK9OS(2 ze2eI7W;?8c!6(N2d}jHQ?g~TW)32PmZh?tq_z&XT|9n`2tFz;FV9;^Z z5f05cpT)z6%pS=s&aJHJUMj%}d+NrNogOFERE}qJK1+vRRu8@K5N{o8MbCQQmMJKU zMZJH^R3Nm)CqPBOa)y0Zv}Lr#q66bWMaMk0Y=v45$nsZI{cJ*VyqWP`K3oftdMgX< z2(hEPou&(i^1PE-IWU9X=FXDMmoKX^vG^7^Vp%MKtUvW#oH2!l%}aYNcV&E;#T+to z7Ht~53)w2wR`gI9R3g>H;*iiRMIkQ%d0yA74R=&-tDr%End9v<_tDQXGt^45SBe9AL-JLu^1csiF z{0Dh&&hf|9{c|W21iHkh`Pl zwbJN7JnlI<2_S6A)CB*|%O!IwhH7m@MJ%k55&n_G*?$5d)=&0v(4@yu)5~0MIrwb69n%yt9uaD%sR>)K!%p-DDh`El++SY zP@8Zk%N(+NPL|k65C?^Y!r0;YZY34Tj)4kW}5hpmeKsa>DTh#hP`E@hOvLF-Wke_QjAHkDT z7~>id=4iRVuaFVA@gXe4p=^lJ)7s)6uqAp(2Z?VaLT%tqkWblMMhxwO+Lc5KQhD4; zN*iKq2O*QR_3E=vt-x}e^2uN(03%BKKP(-7AxbU&6|mT~8o2B}Boj<{qi<~j|BJQk zAR*4-4h7v^|8VG~mD?)pfyO&efO{vn@SIkkwJmEJ&%_Mfw&5!W#!V~)>ODTL-D*2&^WhkdFBic2(rM>sHyX67=a&P@3qjD;NrY5o`(@P7Zhx_@K{$!YK}rG zsDFae;?rVNp6os?`3w|K*U_<4sf3i@<=?@&4A;ABk z=}T+LV+Z?PN_K62UJ|@Da#%nZgu<_AJrGO`P|ZP5HbT63$4FG{fl6DIqTGI+>>slv##7De3%SBRYk|!Ux!hP+QnYe9`&G}f%}K()Umr}}%=Hg&<7#zq z*5rQ>LfM?;lT-S`S#}8_#H$s&Gb+$-7l40}ST|}DfOnbaqO+GxqMw?^0g&vLC4d$O z`H7^SRP)ybeXe{J%Ixl5-^2);UBe$eR0Hus)hB!7oNkZnvWCBIu_(U$(cn30>~_bL zWjN(Fo%;TxxQCH{IVfc5$pf~sM?6}l;l+&{@|l0qrGfoel#y%@!m_>s82{$8f5-Z} z^AE__WTS_cPc-;R|M;6vL`Hm%DjR~nhcW)&6& z*D5OlYrl+<`n$g#)%baE(*A8?7F(Z^lP)fBhSgcfmEy4H>rX`m@oiF`^jkB3@ zxnQB-4~};7##dB*^B8Xja`X#7ZPGw{Z21lCFj!`0REkZ4VP<1vtee1Rra|HSjr?MS zlFQSP{oW`u_U0`PC)1$nbWvW{pSfd>9q-efJ)^4psxaCWXC+h`4<5A$0IEw`k85NNL=2wo>oFw7o+I-FJ+@Tv z3lxo70eRRQ>%^LR!yxy5g<7G?<`A%EacjkeAMYw0AQr%;<+%o5a7hjKe1ku@BrOSa zx%jA(l?Ar6rXXOtQPx~lXsWwg64?<=>DB_y_f)uj0rvI8(E{^tZuX-l0A|sAB)FWX zL{(K_tw$Ow2ftA`A(*FN?4g4>kQIms`exBUi(o1~#`dAb-%oFGZYdYT$spb0yc#BB z;Cxo!yVu(vfc}T#XneKwx80Hk8%1B$p}@l?>9WYrJzJbMz*}RsC?gH}`eznnk^4p) zCZO5-7Z9Jr@I^{nZ2jc^@#KGOjM+D9uyuZ-nRbS1Zo`Tmhvp zzrx&AjvbX(f~;;5FwNL?rQjuO$w}2b-O3AvwmBv4#XTt-dl%xIa`rvuu5^z0Y#HV4 z?wbXvuRYxz??PZ)PHx|>4`_Jt>{h&U4dZEBLJqi#Ai|b0#r24(Kr`9Dk@%s;*e}2? zY@AVOTmv4YWx9QOkfjkenOyy}ERa|uy{DNM%2B+79czqq%-rg!uz%7}SwYk)4tqTD z5}`Z~%6owHg;g`Zb5eLM=Aa@UaGFwsB7>P0|CSTzb%Do|KAIC(pFHQtZj5|;6GkLZ zBUL8xoz>D!w$J*_%3^_~$p$^bZ`3@=ITa*!S{iGKV_5%A0oiwGJd;_`R*CV^w$E#{)?y)y76u>95NJmx~u< zoFV@Ta2eG*Z+<5=_ztEkNF*?Vi2nc^0{X5#XqU+c;$$4aktEv%h0>cZpq6sFZYQ`q zV{?*?*w?QxQ4(kp@ideEOwUsh89Jd0kS1wx6&M1nHu>bUEus71NJi6b57*G@58zpb z&wbq);Pq8IZLJxU!q`r%jw0n5m4N!ot%u}jh-VQFPtI@@+|c^EedlMWqt0B8g(_Ww z*H_gJQcea?!fB^-nu_$-bpYtI&eh^Zi7oWRA%Y$n0&NjUvNmx$#vOn|PZu!|_p@aL z^j9B@H+mz(&f=??_bdfa4aD&0@sJ{p2$r>flgsf< ze)`^Mc#N}_>#Die-^+Ukxc0_}YJUHE^OlobJMMX=eS1_^+RY=p*Dorz{lnrJwAokB z`!aBfC*JfI3VI>LDZnO_WHQfbRt>^X5lxy8K537qN0hIVl;2bl6GGCTP#q$jKX|EC zK3S!Qh7qmj6)aZMn8uCXm5Ul3J7ZF5bCp4+QNLpuD+vE^2bKhI^+e}1Am5!NHRCr&?AS`8 zLn!gBk>3EX)I7hU&44p>jIE_&E!wtz9=X@gLZTwvBTvr>>$UZ+P>Qw3JLv%TzCI2! zU^~G5`GRvlKLw$EWrFHQUYu_zkz7F-mTY}U+OyHyvhIl9~e}QSEQ( z7PgZ0FKg7kkevFrze(mXv1KQ~i=~gCxgSd| za0yw(1k6LPL3XF()~}wo(&UZ^d1v`A`WZ4{Fc!hVK@GQkCR`SD{nt&>GaHhhC6pZo@0_ceUhIiMRV zWd1Jh{KlGJP?M`k?w7Xbw5}<6G52woKQD{(P{|Lo-imEQOqwUJ*Q7r$%){jktQe1^ z2_wpKlc7)EXGd#_I-`}cm3vkORZ>Ur~P({^qu{<5Mx@GCuVpPFa zq+|WpOyA}+Haqr}3VnZNFuDF>7Gwf;ok#sX?j$z^zVlc;BX;Q0&*~$k^-CfREDsGG z4_*P^>Q-K|m4jZ4dewt;HGjLg9&qN(t0JphJ($C%Uq#CZtA~_T64}M$y2_XlQz`fO zXsCzE^*du4dKu=1vWLemD%al~N7iVu9Yt*Tl4FW1CtdE!Dz~~M3WkKc%D$6SZIh4m z8Sgsh)-uG2U0r+1qwZLB;~*8ASOShNu%o1Qq<+x}7qHmZg8gzi`lVYe;g*z(jyKOFO`|{v@qL`5a zkLe>>vzLu;*7jvbcU??A#5VG;?R_55V!N=WRYYmlE1JG^W`R|1)+PMbP3v5zgo2b* zt+U;msd7lp#QOnj*8_893Z~yr~<-ifsQh1^EYjT|Tx;b7LL-;U8uu zrZ>1VvSAkamw=i0D`I)$4*CCPg(D2>+T31y(vm+!U>&b-xn%4tjzC>Jh1r{Ol;@WEB*i5Mgy#fggwfvKywONU&>{>`wI!xI@QT5P zwAkwwV>07ic`FauaUo(N4s5mBTDN_(>>g!BQw^?HWh^gW=kPSch1}k5i- zA`fwXU3qvdcDnk8qp>uZX7pb5LG%3BZ63(}(YS0?VHWTE#LgolHgHbFCDjdO-xr0@8e^k;rJpOc!b9dFT5m0K`a<6G#_SpR*`$q^@%Pf&ob<;3BDM zHAa}OFSStj+9jo-BfaGC@#{Bqb>Q0q+13Rs^_^O3#NfM){ z>QTg|ZK@^6%3h*3o>u7YAo0M27k#x>JJ&v;++p!O>G`a3NZQ?Re<)L2uo7{AE-q6=Nax=g9^>tjyDzu)_*==ty9Ulvu#|PqIvs>5SfX!3;`Li ztnNhIM}#pLoZgC4WD4S@7SW(EF-&8b(G@de7&4KeWY(5L78ITOKvJ>K zFeCp@ga~W;qrdUR&(ATpaA1;(y~YkI_9JgFi}OZp-wM3(5y35k#8zHdGfqe%RgN+G zm@k;?2mYbH=wg`a%jHadG}t$tlR5YjGS{Uk#%&sO&7CZe%#q-l^aCe* z>2h^hu6#zEgWc+;nM@URpX!izx}Jh{_~_W<7$-C(J~5i#M?I{XZu}}K+$urQ(1ws% zxgZx1ij0V~D9f5H;Qdb68u1l&Ad#~4QmEsbkE{RJxb{VBJ|$Vv$aBi ztSaYHUWj^b^A)**Sm8$|#kBkj4oW;vyZ`#!KO*VcpbYo_b^CQiZpc7??=4n3bL+9# zgQ1mvoIds@!q33P65e=-%EhQ!lnQj^sH=+1tOdYA|8mB|4FQLjwQet?BGS1K?E0}v zmHtOhZd;4hXwqwC0TWo0n5Nb!cnfb@{Z-jBWMTUbu8qx60KsP#F-&1kug0a5q zZm7L!hK)kkfS^4K)&?wLXmej;{?CV*fH&@v4Mk!1-4Od7bk0=o z0B|`N+5(P|!S%inlZ6wlG1lU06N&wZh&I1?8IW~S9BGN1?%Ck;abasR%78_#GWM}V z=5aD*gIqSi>k}b#h2EpspGgm=PM$u&z;+&}*eUFxOg&|NHHm_oTM&D)X*xYy^GZb3!Nw15u&4@%IQJj@yYGZ zW4~~tOXXeUo9lZBSfr|j=La_T@L>BJ?=lpvBaw{1zu%DrW_^xgkP6>sdVE%&inkg}l z7bZoh`aaHGYCmY5IQ~<^Utb|_#+5(~kRI)espN0)*)2ulqOR}KyeuW2 zs(eYGw!-`+VEQtBo7tMAWpH*dLpjWzxaeW`CqvGVWlPLypmi*f4YpWu@LlY7^LW>> z2w-(!eVTB;Q>}V-3uS|9fRM7<8N2>bjD}4nq48rZBBHTMML>>Ivk&g~<2VztocbxS z%^w&n?UUvD1T8Ciz~0SHR_f8vbA@W_wom7{ zeOY3|e4Ts&SqV7B7_Nc$dAFA_pZAE}A=&y`&0p`bTRN*OId)pq52*#Coj;w?bBkEfTtL9?KHZN)sDk}$-i0fkCN>}^_)K)u1DEFvy(NkEE zALOr<59)N_8tCO!W8dmE+lZ$A&<#%e3Gi13V;AF4rU0?|G?}^n1nuR{eMM5TqI7G2 zNWl)}p~?+vzwUALPtEE#O4Ezqrb^T$vqi?@uSPv?mNQg`Q6V?@6iBrz*-73l+#grmiw6DlddjPC9^dZwe(xmiu{hr z%0D^Bu=o7zHI`YEuA3hTU^sbpB98)NPhQc4X=8bX4^Q|R9?;)2AGX+!FZ2sAe9FXQ z1|~Yo&GOa3ZC@-elU51edp$Mf^nR0;>eRhohqufM&3gZ)Ey>G`oO|0^BGvj-_2rr= z)_LVZO_aAmOhOJ$hXB-57H zvu5%d!3A6Dt=w=xOW&%@4yRv6F4Id}O)p^SPq`7+hNoM(aKV-5CaihzuRC1uYhq{R z&BICy{zlI~^5Q8iE9QUf?S68(<=w76e!+cFcgKV?XRzrK(+IvRqzacl+Us&ncBJ=- zT+wi7RMGcqEgU6^8cD__rVBB)`0e%r@`IlxOd4A;C)x|Fss>L-XAC7whKoVhRSn}t zV+_1abXlLO8vQ#eY*223Hk7C9Y6@t#7(=DJMl_cUZLY(MJ#Kz0F07)zTP^b=(w?ZA z$*W$ZTl(=_x7e{(oZ8adN*-Yx<}vjSMGC+xJ;uz6g?rmii!!bz2@7xd#jG5H*T$+$ zf%(~Dfrk+*b07PJM?1dxMh^B)9dzyl6hi_{S2m26>3RGT_iXLf=thgjU$P}0&$F%_ zf1FcWJM(F^w}?a6)c8YdJtm|c)3UqAiBtt`g{kg!SmP&ozoKeBB z42dH1U9-eAxRW>h@mf|(OHT2dUC}CUxa*CE!^oB57G?U_m9a|S$WaI7wM5X2RV;{ii3j- zoD2BpEduS=ma_Vi21lOt-9Lf^#8!E% zdc%(k%1ttZ56^}@CM`0}18P#_9kfXC7Yh=#KFM9?>mz<6Y_fgs`f=%OY9s2A$afCO z_dm{fFMsUGRXeahe|^8!no!%|uWdZ_{KtuPZrLsTSv?i3WlQG6<{S*8~7U zcij;SWzj@J5_XIOt|T`yQ604SuO*C(`|Fw@{5UyZGJ~QdnUY_B114zRTV%$XnA~rQ zAr=_~u*w)g)jf8dpumf1X^y6a3twiwYE)k$=r@ptqNPNtJ#On7qh`RIC4L=wTu!CB zPk9m`s8@DrHGb`se=lMPMun2WFt4mhL1C8ov^4B~dM=c3$w`_*li})X{ktom78O-l zz^7cBw4U7g#Vkp1%5nq=Ob4*G5G#C(1^R{HH(8F#475K;S_2Zx!hzhz;G5Skf&`*Vc|3 z#68onrD765Z?*A?k1c*M&YyA8mJ8F5h{A8rJBJ-Bi?V82-IyG>Yi%DSU6~OI*2kXd z(E{^A9fZk+0uSk=Bdb1+od#WUuvSL}vm<&r;kSA1#!ulp_0;ZAD103wl&oJ*eUENbGN#OBzHLFA|$;ItLOK~zH>Zv)&i%U z8{VE?c%z!>g4`#{?##MVm5MmKmt+_!apH7?pi{SzD}ML2Z`9Q6BQw@In@<}$mDKj- zyFcqnJ+zEGB6>C)qz<1B7MozXY6e)ExH8&hfd%sfCq6=^aAb zv1znKuJ^~$l&Cc-6LFZ1u&aU8Ohr!~b^IodujeN%UztX9Wc-{?kLa@M$7&j7wq~?+ zuDh}$r4CW55*~F!Vn|B_1y?4(Hg&_b^9>>`0{S;3u5Z49jyu{!J6s|Jl+?!t6z}H zZmC`;U=ZSL=`f9buGP;kf<|5iJd->$TD^~R-u?D!k-I}2>BROi61!h?paPJYWeu~e z33xKbq(5v`A9-`NXob;QG&T&qOb))`9~W+yy(`^pmoe}1yzRE?qT#HFzX~YHJ`}On zbm2yM$z=gqx?!vfjU!#{3H&NCA$G@Y-I*4fuqfg&{&L5*YMH(2ycYc9Dyqb<=7xeD zhz0cK8u)YEH@5>Bs{`##^D3Wu(d8!1|D9hpnOq` zksJxF2zHug)apcz1t<+u(?VW<=;z2eKfG08rI-lNY2;rGF0`q6g}}q|x}NL_My$(v ziTa(N?9L40`y{G3Ow|h%{P;K+De70F5$5@jS1j!~v+wR*;FcI)`fa?+xM-gGNiO9v z>WO1>H~gQR!nXU{H)I7Q53M?m1g&p5N4Ind&YAIL*o4qknS-CPnss3P9lMKZoGSj+nn^6pk?|4v<3vQ7rqPY3;2-et$>CdjlFdLo`W79+E+b!AWg_lw(ovwh%kUA=;!oVbiOKW#gzQumO! zXb>yMeX?}e#4nGas`hE7NXZQww!Pl1kyg2SVu(e*dOZV3pfd(edJ6)+5VJT3J4+@1&&s-dt*Xnn&~8c7*;t$^O#L5) zX+dS}6-}4rmB%}82Fdq81=ruPm8yFCV5xfOpys&VB{CNtdC!L7*R1-ZU(vViD55!ApG=eS z%}y#LQ>c4zW$>JW%fQWwATDmvpbvwXQ*+vDc&_xP4m{I4hiCgfLv9BL`Asme{%XQO zoDFI+BlAyoot*EEL@M}}n429EUnJ`Ex|>hfVeRp%TwCR+rEt9hRnLk78*x$EVK?W3 z44JVPBNQdHThYOCWnznU^&p}sKFDwMMptM3f=qAtC*ECiqgi8cl#dakr}OCsT)6aT z>33dWmmUW0pufo-l|E^wCLl_I&arOq>Ep;*30umjpKAJg6;>?Bi2*!A&l;+o4(XLp zK!ZY6^>OrU$QHiyE$iuaeYpERyVB&kD_bGsP%>LQoB^)RVbPn|a`s*P6xeJB>9x5W zt6E%@$&ynD=sh`cOgbKa!PL~<+4xrmM*~Z(dnO+m_H(dN!lk9&JIlR1E;m;DAZZ5I++Fvo9`60 z``qWnK_k@6uoSd-NpASFk^PodFzeoJTjC3;Z(0bG>#$$4>khbXo2f$7xpnKTrMnfU zSFDNVW=+PVlgz$SJ*^{{6PI0#>vv_Yphnbbn!3t#*gF4w-;mSh(|W{v8X0p-d(;yp zuAnKDYtl%_UD8bsuZqIktuZ0^?zeAr6UcN&**0aazFlUT6?@Wy8UTwj9{^;fLUx>9 zO5kEmSFA%_(-)DClAaHc&MyCStz3zVipGYlRp2ce1aKXCTFqz$TRYr{DAv+CsLh1%w4aC48PBX#e4Ij~V4F;0*&))4i7{MBvKP{R0`Em8?F%oWP)l&J zc~t(?go`=A$xgFbmdSnoU5MNA!Qf&UieVRzFX;F4N|?|Ym^n}+H56B>@Q_-%*Pkaa z6`npY>WpWS$Lk7q{u$x|Wng7=Sc`<#F*`18;QxChTdijfe_d?Ul28vGhqS=W6k4Rk$YauJ~_3 z`x>XoF7qze=nn&Lqp6|G%6XtBuI<}HsJY9Gj|$sTbLsg*7({ryHM>L7(H)vZ*Tn00 z`{*TqjrPNyZPn-!*9zp>sVOBg zc!K1X*;>eZG7BrksG9Jv)kk8qWK0g@8JC?Ve}8ykso5qQ!@wM^BE0>rfBhk$kbi`K z!@28wrN)t(^%{=G$=ZUw7;hR_+bHvw#+V z=o8(Vm$hV6%j*M+Tv&HbW?1}gdYp77;>46IVj2pSNc7J$W7G*+h$Z+ zHB7PR3QHj?3W*^fw)V_geI~n$;03l1IS=@%Q%$SC4}s0lM`df*Ubjr)Yj^|d{Mt&K z2e0}2s>RHs6{~aFv(I>L4EkTv95~86g{Z-!*+k@Jcc8;v+@F?~RmXWEZFLeTED7mS^ynq|L>ho~^*2tQA)!I~2@| zaIAS2;d{ngU&6~&dU7=fl2J*rEvcWY5n1uORb7lYc%@2Ds-}L5>v#DU)T|jVYc`p% zpS?QvF;Q)8>SL|{mwBcs(5t`!1q4JudJBXmf;8zOO*)}I z!E?{K=id9i`LMHR?U`9KlP~K(d-7Y_4OSW>>5p(JIpZ~12g1jnRxVLZakl=kRQxq7 zL(sI&tFi=j7U`pGh#_p9(B!Th6v3nVc0&fFJ9_ZxK3yioo7q<2Ze{EiWgXBfpl-7pVuY%B!jIRdr z@@sH0=Co662kTzT+Br$!Ji$(?0k|J7%1v z?+;a>hLt|cihcH~l5rp0+EX!WCET>B0+`FhDLE6mpi=&Ps?=O%ZB<$$cp`# zqqS@dg;|6!4q`4@^P6bmtZnAe8_`=4&lMW3Ng39gGT6xH`Xa}0?$54I1eyri>&Im% z+-qT`r)5Yfa~1Bo|Eh`P7=Wca;&Ddxzix|B`7x7V&#Ax(#7LVn5;Ul~sN7ggT8lr=hCd8%rmrTdx|(it!3|85(7TJ};;szZoSOGJH!*JYOO2no^AVXM02Gn2HhG z?UiqUn+5kP3sk<1whgtFNnsW4D-Y`i>PJ7Gw&mdo9KDB8y2-265@m>GyxIz|CKnVM z?`4b>`kd$n&=Ans4;Jvbx&fG(Y0p)v2UKS#2B zES+q9!8I%CsYVu~Y`C8cCIsjCnA<&ju=r_XHg%NWzFrCa;9A0{y}9r)N{(dBf2m?U z>seH!o(W4yp(I~#*B^l&RQB2<8RJ-Y4lzBft@MAOrbb2FEKLg|AWboad5)PP;^UHc zt)6Kcb0yu?^5lVd+%;?Gl5;+ijy4ds3#1M0E?Kp>e6!RnpFtb+V{yp4@1v)Y;~IR( zL*VJ!nL;Jsi-!a1y_~d_7KIfl1y;I|AnfPTXG#;UZT6dHs6z7qa+5x(3EmqVqGa&& z(b&)_U#BIVjOjbxwl@4LRj+kwddPIX-k3_9vrF_eOdEb4QiI*)xZB~-xbm7R0@1S4 zSClO3!&XgVJ}S-|y(|(l(bMd3&EHS>u9t}8{PLjp^xB~fv_^5T^@mG{09wxBsS{#> zc>YCIQb?(mp9p87-TTw9;4uyMXPirW4-=Da3pAAFl~Tk5PD>^r6Gl zgu7Q`#HgPpg#L$2i-vU7@Cuvrs%ma(7J|xSXm*E+p#nhO01@8!~G1JZq)GtmPke`UlYpW5vtR5c^-gaA_-UzVH?S7f(tUZ&XW*$P(#S zQZ5y*pSfE!7!UfUhrEZ5L|tfACg}i+^SqD(Ht@@3Dwh4?ifC2hFAt5o+|J&_sX$+b zc$Mk}2cOKAuF^gnR1cCg1(g4JtamsHhcHSwf8{J@q_!sT_+C@`q=JA+8OIhB8 zG$;@9*9&?|vF8p@Ce5#@bH)oJbQq&|P4iOXyE@oI(Wz2XPRm9V&u=>R=k~rcATVZsk{w}LL8>!BYC&(;SFmrLcnZR7HVekg!OK0OH>T5 zhGqllx-#$$^KzKmVl9 z>{~{j+V$;_f7_Q8{antJH8)t*DcLmPPW)G_YPDi4O~0hC_0XMgsHzD)CF%)pTWdrT zl~-RIzGEI|7OyKT)U}(f2g?F$f&&f9$CIYxH?7Fm3$wN<<|y)z{t+E>Z-ma>#)+ir)v0DugH__*BQZG8TjjZh zJOCzKR8i)g*)-?t;YDa)GzKkt27Xp+#hi$eu6w%V^Te=+o{uJHI7!oVyz)-d>6wlH zboXInvrI;_%yhytxuEBbobx{Ds8tW}AIsgf(?QSWbvjPdB>OWvkNvPKXR9 zYRhxme(m;>H)b4{iQDoAXN00aGo8yTP#HuJUZ#I-0!w(&FJOcZiNb$}V(En8UW;Lh z2kvvWa#DbuFWsp)Vq~K^p}QgUWdqKND>aE;7_2~1L4Az<8?nFBv$yi z_h90?AnjCEVGI+iFCQSEIHXWl{6+IbZyXPdIq;!D+V%ITk&cq=W&i zMx%u1gzT}sdp+29{cEM~nr$%595KR3Sr2gMahuOM$9~GpvNtn!OXgXvFVkqGFk#bU zwXfMCwx7@eL96fiG}JO>jN@iH_jbEY5IA&C66kB}Dn6lg@1&ObmlpkmzpWh|yu%-2 zCu#mMht0^MQ$69;{Bf(aj-&~0>E>ZVl=EG>9BlT z!6=x&G+--?pvCF2-W}axUg#+AFvrrW7VqxvutMybK$JHY>+^|9T3mQemT)q{9XzKZ z3^#mbz#Y||uG8Twl)Cj6wTPDlYzY_x&TVDZpl!W`0~pSB;Wc&L<7`{Cs71fIB8zJ8 z>F%%tNsTx@CW34a;bi@--Ua~S8X~V@MQ0@ z$gUYP885r=1H#h^?kzqmnM$5h(;w4weZs2;w^cuVGhGCARO&7rx{r48Dy-nlhN4m& z)2?S0Zn+~oLI4P3&X*m%6tM?W+jL=HHP1-!SK4+Q4(Yf2MSX$tSvzHU`FEPs8Xt_M zS8D8%awym~2A+nf$|=(8@GtOI=|}rBZ3Dg)sm5N%lZkq3;qUY`I=2CM zE?YRB*V%B`6n_iPr)}O69^zXic`cQE_=Hd5!Aj-|R8dYFbJ0x1N zIhChP-||B_RML#d45x#s$iAZ7K&*yK7v0k|x_)jiXK0XbZcWGQ|KFDnP~?@USRBR>0j#~2k4 zndd`_yimHhC6OnMQ|PGVSe!2tzjda>8DC~5!yY;IfV_}sk7$nir5uq@AT+6PyknN-lZInW+Ro?CJ%(Q7N@8uR;FdNFs!1c@9%5EOya8_ zAe0+aMd0*%jPaJdTALJ+E81Dt151}iBGp;bcP01O1{sv}5tNeFX-w>95bB?gu_$|= zmZP2 zZsf|f5_gRTJ;76wJ>_QfqmjGUkQ(~2)Dq|ER^929(LNlA`B)!L>%Q^JvWX%ewxDHf z)>y)Q5^l*UaGLb>>&;N6-U(lYhS7m}ecU1nB;Jrxqv|vgWJ=*X9QoYqOM;~1q*4P) zc$sd`{*Q0CY-6O+ILb#%d>1}YU*|M)T^3zhXa|qD+et2qAzT_T!hW~+y!&^K^vcC~ zxD9a2lPa6y3#@Th$e8bsy!vpMm0BU*|7Cz-U*Fv-KdeEon>hao2ra<^T^r}=FX`8Z zf)j4N=PW^#m*>hUhw?@+`!XKrK++%g)$@4GV-)2q?@8Y)uKFgm z_P~=r?W4gP%Od(@-?1BK}=x?8DDDWf5_qbW&)Rdo}xQ^p`?hY}ETT z^Gq0%Nmqe+kM=0DMqfo)TGdsVGX@}flszxIiP(J>o3$JFU6OX8XWE@m3%%#ZHm(x- z)U0b)*b)`@3%<`9_({I@o0Vm}4&E5xg{OV&Q9IZjV^(Y1#B?;~Z0Mpra z$}braZ>C8h_oaaLonOo6N2!tWO((~YW_&yIcVkFCi?V7#Nfby@f5{~?wvvnbr)A3U) z#D?vFuH^aYr31YksihCIvJ2l&W8`7d+Q(Zz{&0{_@iFOSd?L@8ArBXv^w*II5P&>L<#3dXq<1+7B2@GfM7I$j$*dxXKeg^!yLam2#Jb`)2pCI_f zp#d*seB}>gMeiW-^UPa)lza8Q$U}-5zV>;68Aks$am_Q6`Q`o2XGcrS+thMX+Xpq6 zriS%T;bbrmj?J-oy|X3Grjzwe|D8(tvl;p0KHb>t1AO=S@i6XP)p<9<|GYEvY<)Q; zK$dwnU{GJn`m}N-%;tUat@d|6M&}Rl0cTUXvHWezoagxCrsH|Ux}aRs>1nzD@#@UU znQ;^TPxE~ zAVuN%6!SUZ>J@NOzsMHQ*8M#w{^0=wqPgiSa$KuMdxF{rBNsD8UN`~K{UtHARCGq! z&B1tcbLJjGY@JW_>}&_W;d6PikO@b zKl{ZyC(7h+Me!?OiY6Yaa6)B8v2v4zpC|xJQzsn}@QY^9%xV48{e|g%?3uht@n6Wq+6fwZ~C9(GuaOBop;7=LYb)~9&Aa&^Gh}k!PUF;oIkO;P+sF4@ z>Fn!yr*Bl#VI=-Wz!shSal7@2Fx6QsW}%vm01?WbZ5m`OU%Wir^qCQ~0U%*$l5JYu zanD>_UVdl@oz<0W^-khH?P8r2?W<J>R_vg@CKR-5}t{R2m?ih#Vr4U zY`s9P{)Jd2kYOPPFn6Fih`^>Lcz}7kTatoQ3Ak3Ej0{K>3>3QPm9UnL2k(BGt3!qA)$ZU(NlkUv`@AE_wUOUbPCcT0B1pe zVPi!l#r~hv*o&q3n;MHr_A)dWrwzq5h&SlHxf78%;AsEovij`?^>_guA(A^9pW}$` zizyqrD(lA=+yHZZf0?Tk|L%T7$R+)Et5*TaBhOr^?Tq0pJ6`mcg$?S zO~E()Mn4k8odjKL(J~=7_<3Zy%0Z0O-0UWR216HM-rYA4IhK1vUEb zo4x`YB_Ox?%zgJ!FmlNqaIg%Z7ui1t`?WC1IUg-a*%%dX1E5>Jaaad##Pw!UUfG-;R3(npMOsDUA^1-*AU1Ai=*l#Z|}#N2~e^Xh5bkP({w-aN5g}?N~uF@o3gLp)nM|lyTOL$PRa6O8wj!@czE^n1rN+pfkU7BiBfP z?Zta|(?)Y*ZC$a3|PUj1okpACp8;(enBncp*} zkG@uq+UEG+yd{&97Mhnb9l1e2HjgsN0+g(MagQ3OHWzKZcM>GtxC<#+iM!bO6ezcK z_0PHDZ!Kkp9x`7)^I)Q-IALoH10|{`_ZDaYwUVpNyQsphp%h;~ILu5tByT;-8?Q?a z|8__*kZf`}pM~LpKNArb4+W6q1(bJS)9+@7+7*0Xi4guotRw~0TQV^xbpMC-5*z|K41OtnWC16{QH&53S41VpVkiGM zPlAAkR4ZwtMRVcZejH%YNq6<#{*caZnL>(HQymMRjDlSa`&5 zWO$^e*sFME_<=*JL$Kq*)3rf~k}Y>K=a^o`4fGqy(B;X|V9`DiRuND^?aG80F!oe=gOZoNA}JLjLVe0 zU+3#5RQ3k;!;qblG8F&U4;}m8EVeCX>PnV9aegjvm(H3Ah}Uq&Q|2d4&5d8rf`W2` zy0ZIhhXy|j4Za)QB981Jg*e+`ycDIMOO|;`E$xnuH=C~cj^XL(3N+I+BYSh7ckCSQ zFf2#c{j77dShdV|xmBZF)8O*R%hD_6MbPmTazFA6wS>;#!eaQ_Wd$7EYSij}Ia&Ev zC-}a`J~yIm0$1PmVjy&mYSCaFIb)0?oL>Fn)q|Txgv@XUR0YJHi=E%SWOixt<=wrfWa0e&_b`&1J8UR|78x*0T>E9V_dil|&`&II_^7{?5qm&=!_~g>$Qst!QT6cPPCU&Yr%AOc^ zSVatnQMGw=LgsVl)8`85j`$|)TxFW&vL$AZ?$4_{j}+h=H9qV!mTZplf zQi{Tq#v|X?w~LJ9MqailzqbOX%_kcq>~M%n9jIrxWe6f{>c+RVdpWoLcABgL&6+I- z&9lraKG?YDFZUch9(bEy;$ys8E?_3x*n#K^tPZYra17w865s>*)(O~I>c$%lrscfb z^ItWOq=`^z7Z#$MqR-0MVo1IRxOp{k*2C?MQc$m~@oMi3y?M)S| z0-NoGGr}6pj;>Bj*4*PV7Qlg+Cz#$hAg@Ezx?@tLQ`SD|*&6HDL{@0N*E9oIj9bF= zs;htcd1&mZM;mLuD6?BH`1ZB_d42iX*4VI7;bUX71DA4}%|4%|qkiJt*AafHev-#b zVyfc%0ewdw>mKG7SFN=A&-$+fhF$UNi|Om*PveJ`xt8R4Y}>Rn_s)0+3{Q19l^_ah zQkX7>vV}4Qq#v&DwJLgJ@`E2*I$y8O&7=BhF)jHeX^4G2W|c{f!Mm!aGrw^D zMepVT^RTQD(qhqjoTMc@hbmKNDHAQ1cF|A(9($HOx|6z1gwp@dB!A+VX>HKOIGKLx=j?)hA`iq>t zv-}B^mHR0x^~}ak440={7d4HY4{g9U-N$xMjREC$Rwo(tQ^2(?B|G?AlZfz7s>1O(0ISAf4{-2Os^6xN|^j|p=R9pzz2pn7pNfIu? z)_<3b9b6P#kPd$tNI}XDc6KmAVuXv#2j)N#QINBPiw6=YDF%AxfpCX8!R*{gg+xR^ zPy`{D0}|-~1i{^GJqa>K|Gx8oM7HSkh_BulwfkUg}5PHfS@OCu#1X> zmbMTGA=d&a2*myK>{uacPPvb2y**9NC-h8I62w=!*1{40=ISZaCWlwa0h}A zaD)rYk>JhtchAU+$ldZENg{r0g1-XsKeI&qRtz8mk0;EVbH zw@vcDZD3JB;s5a=3YPevv4{#w5VY8T^amCd5hm!y|F%i|V|?yzwhm4(H$oy4AjrVM v2S!*cAn2is9fBZAf9t|W2!uQEPmTCnf!Vsb|6WxvVHLro?CdJKs-*u1iNbmt literal 0 HcmV?d00001 diff --git a/tests/testthat.R b/tests/testthat.R index b1a0167..572f3af 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -1,4 +1,4 @@ library(testthat) library(jskm) -test_check("jskm") +test_check("jskm") \ No newline at end of file diff --git a/tests/testthat/test-km.R b/tests/testthat/test-km.R index f3cbcd3..324ccf1 100644 --- a/tests/testthat/test-km.R +++ b/tests/testthat/test-km.R @@ -1,6 +1,7 @@ -context("Kaplan-meier plot") +#context("Kaplan-meier plot") library(survival) +library(ggplot2) # data(colon);data(pbc) test_that("Run jskm", { From 0ad71672a8a2c8f594928c5a768c0aaea7f75d17 Mon Sep 17 00:00:00 2001 From: hyojong myung Date: Tue, 22 Oct 2024 16:58:55 +0900 Subject: [PATCH 2/2] fixed error in df3 --- R/jskm.R | 20 ++++++++++++++------ R/svyjskm.R | 8 ++++---- tests/testthat/test-km.R | 4 ++-- 3 files changed, 20 insertions(+), 12 deletions(-) diff --git a/R/jskm.R b/R/jskm.R index e5600dd..0c2b155 100644 --- a/R/jskm.R +++ b/R/jskm.R @@ -129,7 +129,7 @@ jskm <- function(sfit, n.risk <- n.censor <- surv <- strata <- lower <- upper <- NULL times <- seq(0, max(sfit$time), by = timeby) - if (!is.null(theme) && theme == "nejm") legendposition <- c(1, 0.5) + if (!is.null(theme) && theme == "nejm") legendposition <- legendposition if (is.null(subs)) { if (length(levels(summary(sfit)$strata)) == 0) { subs1 <- 1 @@ -190,7 +190,7 @@ jskm <- function(sfit, L <- summary(sfit)$table["0.95LCL"][[1]] U <- summary(sfit)$table["0.95UCL"][[1]] median_time <- summary(sfit)$table["median"][[1]] - ystratalabs2 <- paste0(ystratalabs, " (median : ", median_time, ", 0.95 CI : ",L, "-", U, ")") + ystratalabs2 <- paste0(ystratalabs, " (median : ", median_time, ", 95%CI : ",L, "-", U, ")") } else { # [subs1] if (is.null(ystratalabs)) @@ -200,7 +200,7 @@ jskm <- function(sfit, L <- summary(sfit)$table[,"0.95LCL"][[i]] U <- summary(sfit)$table[,"0.95UCL"][[i]] median_time <- summary(sfit)$table[,"median"][[i]] - ystratalabs2 <- c(ystratalabs2, paste0(ystratalabs[[i]], " (median : ", median_time, ", 0.95 CI : ",L, "-", U, ")")) + ystratalabs2 <- c(ystratalabs2, paste0(ystratalabs[[i]], " (median : ", median_time, ", 95%CI : ",L, "-", U, ")")) } } } @@ -377,6 +377,10 @@ jskm <- function(sfit, scale_x_continuous(xlabs, breaks = times, limits = xlims) + scale_y_continuous(ylabs, limits = ylims, labels = scale_labels) + + + + if (!is.null(theme) && theme == "jama") { p <- p + theme( panel.grid.major.x = element_blank() @@ -467,6 +471,8 @@ jskm <- function(sfit, } # Add median value + + if (med == TRUE & is.null(cut.landmark) & is.null(status.cmprsk)){ if (length(levels(summary(sfit)$strata)) == 0) { median_time <- summary(sfit)$table["median"][[1]] @@ -506,7 +512,7 @@ jskm <- function(sfit, # Add 95% CI to plot if (ci == TRUE) { - if (med == FALSE | !is.null(status.cmprsk)) { + if (med == FALSE | !is.null(status.cmprsk) | (!is.null(theme) && theme == "nejm")) { if (all(linecols2 == "black")) { p <- p + geom_ribbon(data = df, aes(ymin = lower, ymax = upper), alpha = 0.25, colour = NA) } else if (is.null(col.pal)) { @@ -705,7 +711,8 @@ jskm <- function(sfit, if (!is.null(theme) && theme == "nejm") { p2 <- p1 + coord_cartesian(ylim = nejm.infigure.ylim) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), - axis.text = element_text(size = 10 * nejm.infigure.ratiow) + axis.text = element_text(size = 10 * nejm.infigure.ratiow), + ) + guides(colour = "none", linetype = "none") p <- p + patchwork::inset_element(p2, 1 - nejm.infigure.ratiow, 1 - nejm.infigure.ratioh, 1, 1, align_to = "panel") } @@ -719,4 +726,5 @@ jskm <- function(sfit, } else { p } -} \ No newline at end of file +} + diff --git a/R/svyjskm.R b/R/svyjskm.R index 759140a..ca97a19 100644 --- a/R/svyjskm.R +++ b/R/svyjskm.R @@ -90,7 +90,7 @@ svyjskm <- function(sfit, ...) { surv <- strata <- lower <- upper <- NULL - if (!is.null(theme) && theme == "nejm") legendposition <- c(1, 0.5) + if (!is.null(theme) && theme == "nejm") legendposition <- legendposition if (is.null(timeby)) { if (inherits(sfit, "svykmlist")) { timeby <- signif(max(sapply(sfit, function(x) { @@ -264,7 +264,7 @@ svyjskm <- function(sfit, times <- seq(0, max(sfit$time), by = timeby) if (is.null(ystratalabs)) { ystratalabs <- "All" - ystratalabs2 <- paste0(ystratalabs, " (median : ", unique(df3$med), ")") + ystratalabs2 <- paste0(ystratalabs, " (median : ", unique(df$med), ")") } if (is.null(xlims)) { xlims <- c(0, max(sfit$time)) @@ -289,14 +289,14 @@ svyjskm <- function(sfit, # Final changes to data for survival plot levels(df$strata) <- ystratalabs - zeros <- if (med == T){ + zeros <- if (med == T & is.null(cut.landmark)){ data.frame("strata" = factor(ystratalabs, levels = levels(df$strata)), "time" = 0, "surv" = 1, "med" = 0.5) } else { data.frame("strata" = factor(ystratalabs, levels = levels(df$strata)), "time" = 0, "surv" = 1) } if (ci) { - if (med == T){ + if (med == T & is.null(cut.landmark)){ zeros$med <- NULL zeros$upper <- 1 zeros$lower <- 1 diff --git a/tests/testthat/test-km.R b/tests/testthat/test-km.R index 324ccf1..d40338c 100644 --- a/tests/testthat/test-km.R +++ b/tests/testthat/test-km.R @@ -7,7 +7,7 @@ library(ggplot2) test_that("Run jskm", { fit <- survfit(Surv(time, status) ~ rx, data = colon) expect_is(jskm(fit, timeby = 500, table = T, pval = T), "gg") - expect_is(jskm(fit, + expect_is(jskm(fit, med = T, table = T, pval = T, label.nrisk = "No. at risk", timeby = 365, xlabs = "Time(Day)", ylabs = "Survival", marks = F, xlims = c(0, 3500), ylims = c(0.25, 1), ystratalabs = c("Obs", "Lev", "Lev + 5FU"), ystrataname = "rx" ), "gg") @@ -37,7 +37,7 @@ test_that("Run svyjskm", { s1 <- survey::svykm(Surv(time, status > 0) ~ sex, design = dpbc, se = T) expect_is(svyjskm(s1, ci = T), "gg") expect_is(svyjskm(s1, ci = F, ylabs = "Suvrival (%)", surv.scale = "percent"), "gg") - expect_is(svyjskm(s1, table = T, pval = T, design = dpbc), "gg") + expect_is(svyjskm(s1, table = T, pval = T, design = dpbc, med = T), "gg") expect_is(svyjskm(s1, ci = T, cumhaz = T), "gg") expect_is(svyjskm(s1, ci = F, cumhaz = T), "gg") s2 <- survey::svykm(Surv(time, status > 0) ~ sex, design = dpbc, se = F)