R/tb_tools_present.R

In olssol/tburden: Tumor burden based treatment effect evaluation

#' Spider plot of tumor burden
#'
#' @export
#'
tb_plt_tb <- function(dat_tb, sel_ids = NULL, by_var = c("ARM"),
                      ...,
                      ref_line      = c("none", "mean", "median", "mean_ci"),
                      highlight_obs = FALSE,
                      col           = c("brown", "red")) {

    ref_line <- match.arg(ref_line)
    if (ref_line != "none")
        col <- c(" gray70", "gray30")

    dat_tb_sub <- tkt_subset(dat_tb, ...)
    rst <- ggplot(data = dat_tb_sub, aes(x = DAY, y = PCHG)) +
        geom_line(aes(group = SUBJID), col = col[1])

    if (TRUE == highlight_obs &
        ("day_last_tb" %in% names(dat_tb_sub))) {

        dat_extra <- dat_tb_sub %>%
            left_join(dat_tb_sub %>%
                      filter(DAY <= day_last_tb) %>%
                      group_by(SUBJID) %>%
                      summarize(cut = max(DAY)),
                      by = "SUBJID") %>%
            filter(DAY > cut)

        rst <- rst + geom_line(data = dat_extra,
                      aes(x = DAY, y = PCHG, group = SUBJID),
                      col = col[2])
    }

    if (length(by_var) > 0) {

        if (1 == length(by_var)) {
            s_fml <- paste("~", paste(by_var, collapse = "+"))
            rst   <- rst + facet_wrap(as.formula(s_fml))
        } else {
            s_fml <- paste(by_var[1], "~", paste(by_var[-1],
                                                 collapse = "+"))
            rst   <- rst + facet_grid(as.formula(s_fml))
        }
    }

    ## selected pt
    d_sel <- dat_tb %>%
        filter(SUBJID %in% sel_ids)

    if (nrow(d_sel) > 0)
        rst <- rst +
            geom_line(data = d_sel,
                      aes(x = DAY, y = PCHG, group = SUBJID),
                      col = "green",
                      lwd = 1.5)

    ## ref curve
    dat_ref <- NULL
    if (ref_line != "none") {
        f_ref <- switch(ref_line,
                        mean    = mean,
                        median  = median,
                        mean_ci = mean)

        dat_ref <- dat_tb_sub %>%
            group_by_at(c(by_var, "DAY")) %>%
            summarize(ref_y = f_ref(PCHG),
                      n     = n(),
                      lb    = mean(PCHG) - 1.96 * sd(PCHG) / sqrt(n),
                      ub    = mean(PCHG) + 1.96 * sd(PCHG) / sqrt(n))

        if (length(by_var) > 0) {
            dat_ref$Overlay_Group <- apply(dat_ref[, by_var],
                                           1,
                                           paste,
                                           collapse = "|")
        }

        rst <- rst +
            geom_line(data = dat_ref,
                      aes(x = DAY, y = ref_y), col = "brown",
                      lwd = 1.5)

        if ("mean_ci" == ref_line) {
            rst <- rst +
                geom_ribbon(data = dat_ref,
                            aes(x = DAY, y = ref_y, ymin = lb, ymax = ub),
                            col = "yellow",
                            alpha = 0.4)
        }
    }

    rst <- rst +
        theme_bw() +
        theme(legend.position = "none")


    list(plot    = rst,
         dat_ref = dat_ref)
}

#' Histogram of TB at a given day for all patients
#'
#' @export
#'
tb_plt_tb_histogram <- function(dat_tb, day, by_var = c("ARM"), ...) {

    xx     <- unique(dat_tb$DAY)
    xx_inx <- which.min(abs(xx - day))

    dat_tb_sub <- tkt_subset(dat_tb, ...) %>%
        filter(DAY == xx[xx_inx])

    dat_sum <- dat_tb_sub %>%
        group_by(across(all_of(by_var))) %>%
        summarize(mean = mean(PCHG),
                  median = median(PCHG))

    s_fml <- paste("~", paste(by_var, collapse = "+"))

    rst <- ggplot(data = dat_tb_sub, aes(x = PCHG)) +
        geom_histogram() +
        geom_density() +
        geom_vline(aes(xintercept = mean),
                   data = dat_sum, col = "red") +
        geom_vline(aes(xintercept = median),
                   data = dat_sum, col = "brown", lty = 2) +
        facet_wrap(as.formula(s_fml)) +
        theme_bw() +
        theme(legend.position = "none")
}


#' Histogram of last slope
#'
#' @export
#'
tb_plt_tb_slope <- function(dat_tb, by_var = c("ARM"), ...) {

    dat_tb_sub <- tkt_subset(dat_tb, ...)
    s_fml      <- paste("~", paste(by_var, collapse = "+"))

    rst <- ggplot(data = dat_tb_sub, aes(x = slope)) +
        geom_histogram() +
        geom_density() +
        facet_wrap(as.formula(s_fml)) +
        theme_bw() +
        theme(legend.position = "none")
}


#' Spider plot of tumor burden for individual subject with multiple imputations
#'
#'
#' @export
#'
tb_plt_tb_ind <- function(ind_dat) {

    ind_mean <- ind_dat %>%
        group_by(x) %>%
        summarize(y = mean(y))

    rst <- ggplot(data = ind_dat, aes(x = x, y = y)) +
        geom_line(aes(group = imp), lty = 2, col = "gray20") +
        geom_line(data = ind_mean, aes(x = x, y = y),
                  col = "red") +
        theme_bw() +
        theme(legend.position = "none")
}


#' Survival curves
#'
#' @export
#'
tb_plt_km <- function(dat_surv, type = c("PFS", "OS"), ...) {

    type       <- match.arg(type)
    var_status <- paste(type, "_", "CNSR", sep = "")
    var_time   <- paste(type, "_", "DAYS", sep = "")

    plot_km(dat_surv, var_time, var_status, lab_y = type, ...)$plot
}


#' Survival curves for imputed survival
#'
#' @export
#'
tb_plt_km_imp <- function(imp_surv, dat_surv, inx_imp = NULL,
                          type = c("PFS", "OS"), ...) {
    type     <- match.arg(type)
    dat_surv <- imp_surv %>%
        left_join(dat_surv, by = "SUBJID") %>%
        mutate(status = 0)

    if (!is.null(inx_imp)) {
        dat_surv <- dat_surv %>%
            filter(Imp == inx_imp)
    }

    stopifnot(nrow(dat_surv) > 0)

    if ("PFS" == type) {
        dat_surv$time <- apply(dat_surv[, c("IT_PFS", "IT_OS")], 1,
                               function(x) min(x, na.rm = TRUE))
    } else {
        dat_surv$time <- dat_surv$IT_OS
    }

    plot_km(dat_surv, "time", "status", lab_y = type, ...)$plot
}



#' Plot correlation of utilities
#'
#'
#'
#'
#' @export
#'
tb_plt_estimate <- function(rst_estimate, var1 = "uti_tb", var2 = "uti_event") {

    rst_estimate$x <- rst_estimate[[var1]]
    rst_estimate$y <- rst_estimate[[var2]]

    sum_lm <- rst_estimate %>%
        group_by(imp, ARM) %>%
        summarize(R2 =  cor(x, y)) %>%
        mutate(R2 = round(R2, 3))

    ggplot(data = rst_estimate, aes(x = x, y = y)) +
        geom_point() +
        geom_smooth(method = "lm", se = FALSE) +
        theme_bw() +
        facet_grid(imp ~ ARM) +
        labs(x = var1, y = var2) +
        geom_label(data = sum_lm,
                   aes(x = -Inf, y = Inf,
                       label = paste("R2 = ", R2, sep = " ")),
                   hjust = 0, vjust = 1)
}

## -----------------------------------------------------------------
##
##                SURVIVAL  PRESENTATION
##
## -----------------------------------------------------------------

#' Plot survival curve with area under the curve
#'
#' @export
#'
tb_plt_surv <- function(surv_f, t_dur = NULL,
                        type = c("rmf", "rmst", "none"),
                        y_lim = c(0, 1), x_lim = NULL) {

    type     <- match.arg(type)
    surv_dur <- tb_surv_cut(surv_f, t_dur)$surv_f_dur
    surv_km  <- tb_surv_cut(surv_f, x_lim)$surv_f_dur

    ## survival curves
    rst    <- ggplot(data = data.frame(Time = surv_km[, 1],
                                       Y    = surv_km[, 2]),
                     aes(x = Time, y = Y)) +
        labs(x = "Time", y = "Survival Probability") +
        ylim(y_lim) +
        theme_bw() +
        geom_step()

    if (is.null(x_lim)) {
        rst <- rst + xlim(c(0, x_lim))
    }

    if (type == "none")
        return(rst)

    ## ploygon
    if (!is.null(t_dur)) {
        y_dur    <- switch(type,
                           rmst = rbind(c(surv_dur[nrow(surv_dur), 1],
                                          0),
                                        c(0, 0)),
                           rmf  = c(surv_dur[nrow(surv_dur), 1], 1))

        surv_poly <- NULL
        for (i in 1:(nrow(surv_dur) - 1)) {
            surv_poly <- rbind(surv_poly,
                               surv_dur[i, ],
                               c(surv_dur[i + 1, 1], surv_dur[i, 2]))
        }

        surv_poly <- rbind(surv_poly, y_dur)
        rst <- rst + geom_polygon(data = data.frame(x = surv_poly[, 1],
                                                    y = surv_poly[, 2]),
                                  aes(x = x, y = y),
                                  alpha = 0.2)
    }

    rst
}

#' Plot patients by enrollment
#'
#' @export
#'
tb_plt_onstudy <- function(t_enroll, t_time, event, t_dur,
                           add_auc = FALSE, auc_k = 1.2,
                           add_lab = TRUE, size_lab = 8, hjust_lab = -1,
                           h = 0.4) {

    lab_e <- c("Censored", "Event", "Enrolled")
    dat <- data.frame(t_enroll = t_enroll,
                      time     = t_time,
                      event    = event) %>%
        arrange(t_enroll) %>%
        mutate(y     = row_number(),
               event = factor(event, 0:2, lab_e))

    rst <- ggplot(data = dat, aes(x = time, y = y)) +
        geom_point(aes(pch = event, color = event)) +
        geom_vline(xintercept = t_dur, lty = 2) +
        geom_vline(xintercept = 0, lty = 2) +
        geom_text(aes(x = 0, y = 5, label = "Study Started"),
                  angle = 90, vjust = -0.5) +
        geom_text(aes(x = t_dur, y = 5, label = "Study Finished"),
                  angle = 90, vjust = 1) +
        labs(xlim = c(-0.1, t_dur * 1.05), lty = 2) +
        theme_bw() +
        theme(
            axis.line    = element_blank(),
            axis.text    = element_blank(),
            axis.ticks   = element_blank(),
            axis.title   = element_blank(),
            panel.grid   = element_blank(),
            legend.title = element_blank(),
            legend.position = "bottom")

    ## geom_point(data = data.frame(x     = dat$t_enroll,
    ##                              y     = dat$y,
    ##                              event = factor(2, 0:2, lab_e)),
    ##            aes(x     = x,
    ##                y     = y,
    ##                pch   = event,
    ##                color = event)) +

    ## add line
    for (i in 1:nrow(dat)) {
        pt    <- dat[i, "time"]
        t_enr <- dat[i, "t_enroll"]
        pe    <- "Event" == dat[i, "event"]
        pd    <- data.frame(x = c(t_enr, pt),
                            y = c(i, i))
        rst <- rst +
            geom_line(data = pd, aes(x = x, y = y))

        if (add_auc) {
            if (pe) {
                rst <- rst +
                    geom_line(data = data.frame(x = c(pt, pt, t_dur),
                                               y = c(i, i + h, i + h)),
                              aes(x = x, y = y))
                pt_imp <- pt
            } else {
                pt_imp <- runif(1, pt, pt + auc_k * (t_dur - pt))
                pt_imp <- min(pt_imp, t_dur)
                rst    <- rst +
                    geom_line(data = data.frame(
                                  x = c(pt, pt_imp, pt_imp, t_dur),
                                  y = c(i,  i, i + h, i + h)),
                              aes(x = x, y = y), lty = 2)
            }

            rst <- rst +
                geom_polygon(data = data.frame(
                                 x = c(pt_imp, pt_imp, t_dur, t_dur),
                                 y = c(i, i + h, i + h, i)),
                             aes(x = x, y = y),
                             fill = "gray30",
                             alpha = 0.2)
        }
    }

    ## add pt label
    if (add_lab) {
        dat_lab <- data.frame(x    = dat$t_enroll,
                              y    = dat$y,
                              labs = paste("P", dat$y, sep = ""))

        rst <- rst +
            geom_label(data = dat_lab,
                       aes(x = x, y = y, label = labs),
                       hjust = hjust_lab,
                       size  = size_lab)
    }

    ## return
    rst
}

#' Summarize utility by censoring
#'
#'
#'
#' @export
#'
tb_est_by_censor <- function(rst_orig, dat_surv = NULL) {

    if (is.null(dat_surv))
        dat_surv <- rst_orig$params$dat_surv

    dat_est  <- rst_orig$estimate_sub$estimate %>%
        left_join(dat_surv %>% select(SUBJID, PFS_CNSR)) %>%
        mutate(PFS = 0 == PFS_CNSR) %>%
        group_by(ARM, PFS) %>%
        summarize(adj_utility = mean(adj_utility))

    dat_est
}


#' Plot survival curve with area under the curve
#'
#' @export
#'
tb_plt_fu <- function(dat_surv, by_var = c("ARM"), event = "PFS",
                      date_dbl = NULL, id = NULL,
                      ...) {

    dat_surv$v_days   <- dat_surv[[paste(event, "_DAYS", sep = "")]]
    dat_surv$v_status <- 0 == dat_surv[[paste(event, "_CNSR", sep = "")]]

    if (!is.null(id)) {
        tmp <- dat_surv %>%
            filter(SUBJID == id)

        if (1 == nrow(tmp)) {
            id_days <- tmp[1, "v_days"]
            id_arm  <- tmp[1, "ARM"]

            dat_surv <- dat_surv %>%
                filter((v_days > id_days  &
                        1      == v_status &
                        ARM    == id_arm) |
                       SUBJID == id)
        }
    }

    dat_surv <- dat_surv %>%
        group_by(across(all_of(by_var))) %>%
        arrange(v_days) %>%
        mutate(ID = row_number())

    if (!is.null(date_dbl)) {
        t_ana          <- as.Date(date_dbl) - as.Date(dat_surv$RANDT)
        dat_surv$t_ana <- as.numeric(t_ana)

        dat_fu <- dat_surv %>%
            mutate(v_days   = t_ana - v_days,
                   v_status = -1) %>%
            filter(v_days > 0)

        dat_surv <- rbind(dat_surv, dat_fu)
    }

    g_colors <- c("Until DBL"  = "gray80",
                  "Censored"   = "#CC6600",
                  "With Event" = "blue")

    dat_plt <- dat_surv %>%
        mutate(Type = factor(v_status,
                             c(-1, 0, 1),
                             c("Until DBL", "Censored", "With Event")))

    s_fml  <- paste("~", paste(by_var, collapse = "+"))
    ggplot(data = dat_plt, aes(x = ID, v_days)) +
        geom_bar(stat = "identity", aes(fill = Type)) +
        coord_flip() +
        theme_bw() +
        facet_wrap(as.formula(s_fml)) +
        labs(y = "Follow-up Days", x = "Patients") +
        scale_fill_manual(values = g_colors)
}