R/PlotSTEPP.R
In lengning/gClinBiomarker: gClinBiomarker: an R package for clinical biomarker analyses
#' Create a STEPP (Subpopulation Treatment Effect Pattern Plot)
#'
#' This function creates a STEPP from the given point estimates and confidence intervals at desired percentiles.
#'
#' @author Alexey Pronin \email{pronin.alexey@gene.com}, Ning Leng \email{leng.ning@gene.com}, and previous team members (see DESCRIPTION)
#'
#' @param data input data frame. Rows are patients and columns are variables (e.g. demographics variables, time to event variables,
#' biomarker variables, treatment indicator, etc.). One patient per row.
#' @param outcome.var outcome variable. In case of a 'survival', variable, it will be a vector of two variables: 1) time to event 2) censorship
#' @param outcome.class outcome class of the 'outcome' variable. Can be either "continuous", "binary", or "survival".
#' @param trt name of the treatment variable.
#' @param var name of the biomarker variable.
#' @param covariate vector specifying the covariate variables.
#' This can be added to adjust for in the analysis for survival and continuous outcome variable classes.
#' @param strata vector specifying the stratification variables. This can be added for the survival outcome variable class.
#' @param placebo.code name of the control group within the treatment variable
#' @param active.code name of the treatment/experimental group within the treatment variable
#' @param quantile.type an integer between 1 and 9 selecting one of the nine quantile algorithms. See \code{\link{quantile}}. Default is 2.
#' @param alpha confidence level (CI) for point estimate, i.e. 0.05 for 95 percent CI. Default is 0.05.
#' @param min.pt minimum center.pt. Default is NULL.
#' @param max.pt maximum center.pt. Default is NULL.
#' @param window.width width of each window. Default is 0.25.
#' @param by size of 'slide', i.e. speed of the window that moves along the x-axis. Default is 0.05.
#' @param yrange.lower value of the lower y-axis. Default is NULL.
#' @param yrange.upper value of the upper y-axis. Default is NULL.
#' @param xticks x-tick marks. Default is NULL.
#' @param show.refline if TRUE, the reference line will be displayed. Default is TRUE.
#' @param refline.color color of the reference line. Default is "grey".
#' @param show.refline.ac if TRUE, the reference line for effect in All Comers will be displayed. Default is TRUE.
#' @param refline.color.ac color of the reference line for effect in All Comers. Default is "lightblue".
#' @param estimate.color color of the estimate line. Default is "blue".
#' @param estimate.lty type of the estimate line. Default is 1.
#' @param estimate.lwd width of the estimate line. Default is 2.
#' @param surv.conf.type confidence interval type. Default is "plain". see conf.type in survfit.
#' @param ci.color color of the CI lines. Default is "black".
#' @param ci.lty type of the CI lines. Default is 2.
#' @param ci.lwd width of the CI lines. Default is 1.
#' @param ci.shade if TRUE, the area between the lower and upper CI lines will be shaded in gray. Default is TRUE.
#' @param plot.title title of the plot. Default is "STEPP: Subgroup Treatment Effect Pattern Plot".
#' @param sub.title subtitle of the plot, this will be displayed on the bottom of the plot. Default is NULL.
#' @param xlabel x-axis label. Default is "Biomarker Percentile".
#' @param ylabel y-axis label. Default is NULL.
#' @param show.legend if TRUE, a legend will be displayed. Default is TRUE.
#' @param legend.loc location of the legend to be displayed. Default is "topright".
#' @param legend.text text to be displayed in the legend.
#' Default is \code{c("Point Estimate",paste((1-alpha)*100," percent Confidence Interval",sep=""))}
#' @param legend.col colors of the legend. Default is \code{c(estimate.color, ci.color)}.
#' @param legend.lty type of the legend lines. Default is \code{c(estimate.lty, ci.lty)}.
#' @param legend.lwd width of the legend lines. Default is \code{c(estimate.lwd, ci.lwd)}.
#' @param legend.bty type of box for the legend. Default is "n".
#' @param bm.digits digits to be displayed/used for the lower and upper confidence level estimates. Default is 2.
#' @param actual.scale if TRUE, it generates the figure using the actual scale instead of the scale in percentage. Default is FALSE.
#' @param equal.in.LL,equal.in.UL if both are TRUE, window is defined using >= and <= (in legend: "\code{[ ]}").
#' if (TRUE, FALSE), window is defined using >= and < (in legend: "\code{[ )}").
#' if (FALSE, TRUE), window is defined using > and <= (in legend: "\code{( ]}").
#' if both are FALSE, window is defined using > and < (in legend: "\code{()}"). Default is (TRUE, FALSE).
#' @param csv.name csv file name (includes numbers used in the graphs). If NULL (default), the function will return the result.
#' @param pdf.name name of output pdf file. If it's NULL (default), the plots will be displayed but not saved as pdf.
#' @param pdf.param a list of parameters that define pdf graphics device. See \code{\link{pdf}}. Default is \code{list(width=11, height=8.5)}.
#' @param par.param a list of parameters that define graphcial parameters. See \code{\link{par}}. Default is \code{list(mar=c(4,4,3,2))}.
#' @param ties Default is "efron". To match internal sas results, use "exact". See parameter "ties" in coxph.
#' @note Patient data without corresponding biomarker data are automatically removed.
#' For survival data, censorship variable is 1 if an event happened, 0 if censored.
#'
#' @importFrom survival coxph Surv strata
#'
#' @examples
#'
#' data(input)
#' PlotSTEPP(data = input,
#' outcome.var = c("PFS", "PFS.event"),
#' outcome.class = "survival",
#' trt = "Arm",
#' var = "KRAS.exprs",
#' covariate = "Sex",
#' strata = "Age",
#' placebo.code = "CTRL",
#' active.code = "TRT")
#'
#' @seealso
#' \link{SummaryTwoGroups}
#'
#' @export
PlotSTEPP <- function(data,
outcome.var,
outcome.class,
trt=NULL,
var,
covariate = NULL,
strata = NULL,
placebo.code = NULL,
active.code = NULL,
quantile.type = 1,
alpha = .05,
window.width = .25,
min.pt = NULL,
max.pt = NULL,
by = .05,
yrange.lower=NULL,
yrange.upper=NULL,
xticks = NULL,
show.refline = TRUE,
refline.color = "grey",
show.refline.ac = TRUE,
refline.color.ac = "lightblue",
estimate.color="blue",
estimate.lty = 1,
estimate.lwd = 2,
surv.conf.type="plain",
ties="efron",
ci.color = "black",
ci.lty = 2,
ci.lwd = 1,
ci.shade = TRUE,
plot.title = "STEPP: Subgroup Treatment Effect Pattern Plot",
sub.title = NULL,
xlabel = "Biomarker Percentile",
ylabel = NULL,
show.legend = TRUE,
legend.loc = "topright",
legend.text = c("Point Estimate", paste((1-alpha)*100," percent Confidence Interval",sep="")),
legend.col = c(estimate.color, ci.color),
legend.lty = c(estimate.lty, ci.lty),
legend.lwd = c(estimate.lwd, ci.lwd),
legend.bty = "n",
bm.digits = 2,
actual.scale = FALSE,
equal.in.LL = TRUE,
equal.in.UL = FALSE,
pdf.name = NULL,
pdf.param = list(width= 11, height=8.5),
par.param = list(mar=c(4,4,3,2)),
csv.name = NULL) {
stopifnot(class(data)=="data.frame")
outcome.class <- match.arg(outcome.class, c("survival", "binary", "continuous"))
stopifnot(all(c(var, outcome.var, trt)%in%colnames(data)))
if(any(is.na(data[[var]])))message("some NA in var column, will ignore NA entries")
# Remove patient data without corresponding biomarker value
data.ori <- data
data <- data[which(!is.na(data[[var]])),]
# Read off the data
Outcome <- data[, outcome.var]
Biomarker <- data[, var]
if(!is.null(trt))if(length(unique(data[,trt]))==1) trt <- NULL
if(is.null(trt)) nArms <- 1
if (!is.null(trt)) { # multi-arm study
Treatment <- data[,trt]
Arms <- levels(factor(Treatment))
if (!is.null(placebo.code)) { # reference arm specified
if(length(placebo.code)!=1)stop("placebo.code should have length 1")
if(!placebo.code %in% Arms)stop("placebo.code should be an element in treatment column")
Arms <- c(placebo.code, setdiff(Arms, placebo.code)) # first placebo then others
if (!is.null(active.code)) {
if(length(active.code)!=1)stop("active.code should have length 1")
if(!all(active.code %in% Arms))stop("active.code should be elements in treatment column")
if(length(intersect(placebo.code, active.code))>1)
stop("code cannot be in both active.code and placebo.code!")
Arms <- c(Arms[1], active.code) #order by specified input
}
}
nArms <- length(Arms)
if(!nArms%in%c(2))stop("only 2-arm is allowed")
data[,trt] <- factor(data[,trt],levels=Arms)
if(is.null(placebo.code))placebo.code <- Arms[1]
if(is.null(active.code))active.code <- Arms[-1]
}
if (!is.null(covariate)) {
Covariate <- data[, covariate]
}
if (!is.null(strata)) {
Strat.factor <- data[, strata]
}
# If min.pt or max.pt are NULL
if (is.null(min.pt)) {
min.pt <- window.width/2
}
if (is.null(max.pt)) {
max.pt <- 1 - window.width/2
}
nbins <- round((max.pt - min.pt)/by) + 1
# X-axis ticks (if not specified by user)
if (is.null(xticks)) {
if (nbins > 4) {
allticks <- seq(min.pt, max.pt, by = by)
if (nbins %% 2 == 0) {
xticks <- c(min.pt, quantile(allticks, c(.33, .67), type = quantile.type), max.pt)
} else {
xticks <- c(min.pt, quantile(allticks, c(.25, .5, .75), type = quantile.type), max.pt)
}
} else {
xticks <- seq(min.pt, max.pt, by = by)
}
}
# Calculate summary statistics for each bin
cn <- c("Center.pt", "Effect.Size", "Lower", "Upper", "BMV.LL", "BMV.UL", "BMV.Center", "Left.pt", "Right.pt", "N")
if (outcome.class == "survival") {
cn <- c(cn, "Events")
}
sdata <- matrix(NA, nbins, length(cn))
colnames(sdata) <- cn
sdata[, "Center.pt"] <- seq(min.pt, max.pt, by)
# Effect size in All Comers
if (is.null(covariate) & is.null(strata)) {
effect.ac <- SummaryTwoGroups(Outcome, 1:length(Biomarker), Treatment,
placebo.code, active.code, outcome.class,
alpha, surv.conf.type=surv.conf.type, ties=ties,
covariate.var = NULL,
strat.factor.var = NULL)["Effect.Size"]
} else if (!is.null(covariate) & is.null(strata)) {
effect.ac <- SummaryTwoGroups(Outcome, 1:length(Biomarker), Treatment,
placebo.code, active.code, outcome.class,
alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = Covariate,
strat.factor.var = NULL)["Effect.Size"]
} else if (is.null(covariate) & !is.null(strata)) {
effect.ac <- SummaryTwoGroups(Outcome, 1:length(Biomarker), Treatment,
placebo.code, active.code, outcome.class,
alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = NULL,
strat.factor.var = Strat.factor)["Effect.Size"]
} else if (!is.null(covariate) & !is.null(strata)) {
effect.ac <- SummaryTwoGroups(Outcome, 1:length(Biomarker), Treatment,
placebo.code, active.code, outcome.class,
alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = Covariate,
strat.factor.var = Strat.factor)["Effect.Size"]
}
options(scipen=999)
for (i in 1:nbins) {
# Define this window
START_SEQ <- seq( 0 , 1 - min.pt*2 , by = by)
start <- START_SEQ[[i]]
end <- start + window.width
# Save in sdata
sdata[i, "Left.pt"] <- start
sdata[i, "Right.pt"] <- end
# Biomarker values corresponding to the percentiles
# 'type'=2 is used to match the result from median()
LL <- quantile(Biomarker, start, type = quantile.type)
UL <- quantile(Biomarker, end, type = quantile.type)
# Round LL and UL
LL <- round(LL, bm.digits)
UL <- round(UL, bm.digits)
if (LL == 0 & UL == 0) {
stop("No enough unique values, please consider to reduce number of windows (increase window.width) or increase number of digits when rounding (bm.digits)")
}
sdata[i, c("BMV.LL", "BMV.UL")] <- c(LL, UL)
# The subgroup index for this window
if (equal.in.LL == TRUE & equal.in.UL == FALSE) {
sindex <- which(Biomarker >= LL & Biomarker < UL)
} else if (equal.in.LL == TRUE & equal.in.UL == TRUE) {
sindex <- which(Biomarker >= LL & Biomarker <= UL)
} else if (equal.in.LL == FALSE & equal.in.UL == TRUE) {
sindex <- which(Biomarker > LL & Biomarker <= UL)
} else if (equal.in.LL == FALSE & equal.in.UL == FALSE) {
sindex <- which(Biomarker > LL & Biomarker < UL)
} else {
stop("equal.in.LL and equal.in.UL can be equal to TRUE or FALSE only!")
}
# Calculate median for each bean for Biomarker
sdata[i, "BMV.Center"] <- median(Biomarker[sindex])
# Sample size in each bin
sdata[i, "N"] <- length(sindex)
if (outcome.class == "survival") {
sdata[i, "Events"] <- sum(Outcome[sindex, 2], na.rm = TRUE)
}
if (is.null(covariate) & is.null(strata)) {
sdata[i, c("Effect.Size", "Lower", "Upper")] <-
SummaryTwoGroups(Outcome, sindex, Treatment,
placebo.code, active.code,
outcome.class, alpha,
covariate.var = NULL,
strat.factor.var = NULL)[c("Effect.Size", "Lower", "Upper")]
} else if (!is.null(covariate) & is.null(strata)) {
sdata[i, c("Effect.Size", "Lower", "Upper")] <-
SummaryTwoGroups(Outcome, sindex, Treatment,
placebo.code, active.code,
outcome.class, alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = Covariate,
strat.factor.var = NULL)[c("Effect.Size", "Lower", "Upper")]
} else if (is.null(covariate) & !is.null(strata)) {
sdata[i, c("Effect.Size", "Lower", "Upper")] <-
SummaryTwoGroups(Outcome, sindex, Treatment,
placebo.code, active.code,
outcome.class, alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = NULL,
strat.factor.var = Strat.factor)[c("Effect.Size", "Lower", "Upper")]
} else if (!is.null(covariate) & !is.null(strata)) {
sdata[i, c("Effect.Size", "Lower", "Upper")] <-
SummaryTwoGroups(Outcome, sindex, Treatment,
placebo.code, active.code,
outcome.class, alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = Covariate,
strat.factor.var = Strat.factor)[c("Effect.Size", "Lower", "Upper")]
}
} # end of for (i in 1:nbins)
##### Print out sdata
sdata2 <- sdata
colnames(sdata2)[1:9] <- c("Window Center", "Effect.Size", "CI Lower", "CI Upper",
"BM Lower", "BM Upper", "BM Center", "Window Left", "Window Right")
if (!is.null(ylabel)) {
colnames(sdata2) <- gsub("Effect.Size", ylabel, colnames(sdata2))
} else {
if (outcome.class == "survival") {
colnames(sdata2) <- gsub("Effect.Size", "Hazard Ratio", colnames(sdata2))
} else if (outcome.class == "binary") {
colnames(sdata2) <- gsub("Effect.Size", "Proportion Difference", colnames(sdata2))
} else if (outcome.class == "continuous") {
colnames(sdata2) <- gsub("Effect.Size", "Mean Difference", colnames(sdata2))
}
}
if (!is.null(csv.name)) {
write.csv(sdata2, file = csv.name, row.names = FALSE)
}
##### Draw the plot
PlotParam(pdf.name, pdf.param, par.param)
if(is.null(sub.title))sub.title <- paste("soc:", placebo.code, "; trt:", active.code)
center.pt <- sdata[, "Center.pt"]
effect.size <- sdata[, "Effect.Size"]
lower = sdata[, "Lower"]
upper = sdata[, "Upper"]
bml = sdata[, "BMV.LL"]
bmu = sdata[, "BMV.UL"]
if (outcome.class == "survival") {
num.events = sdata[, "Events"]
}
num.pats = sdata[, "N"]
if (actual.scale == FALSE) {
# If center.pt is NULL
if (is.null(center.pt)) {
center.pt <- seq(min.pt, max.pt, by = by)
}
# x values should always be between 0 and 1
if (sum(center.pt <= 0 | center.pt >= 1) > 0) {
warning("Center percentile must be between 0 and 1 !")
}
}
# Starting point should be smaller than Ending point
if (min.pt < 0 | min.pt > max.pt | max.pt > 1 ) {
warning("Inappropriate center percentile; Please check.")
print(paste("window.width", window.width))
print(paste("min.pt", min.pt))
print(paste("max.pt", max.pt))
}
if (is.null(yrange.lower)) {
yrange.lower= min(lower, na.rm=TRUE)
}
if (is.null(yrange.upper)) {
yrange.upper= max(upper, na.rm=TRUE)
}
if (outcome.class == "binary") {
if(is.null(ylabel)) {
ylabel <- "Proportion Difference (%)"
}
if (yrange.lower < -1 | yrange.upper > 1) {
warning("Input yrange.lower and yrange.upper should be between -1 and 1 !")
}
} else if (outcome.class == "survival") {
if (is.null(ylabel)) {
ylabel <- "Hazard Ratio"
}
if (yrange.lower < 0 | yrange.upper < 0) {
warning("Input yrange.lower and yrange.upper should be >0 !")
}
} else if (outcome.class == "continuous") {
if (is.null(ylabel)) {
ylabel <- "Mean Difference"
}
}
# Set an empty panel
if (outcome.class == "binary") {
plot(c(ifelse(actual.scale == FALSE, 0, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, 1, sdata[nbins, "BMV.Center"])),
c(yrange.lower, yrange.upper),
type = "n",
axes = FALSE,
xlim = c(ifelse(actual.scale == FALSE, min.pt, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, max.pt, sdata[nbins, "BMV.Center"])),
xlab = xlabel,
ylab = ylabel,
main = plot.title,
sub = sub.title)
# y-axis in % scale
axis(2,
-10:10/10,
paste(c(rep("", 11), rep("+", 10)), -10:10*10, sep = ""),
las = 2)
h_ <- 0
text_ <- paste(round(effect.ac*100), "%", sep = "")
} else if (outcome.class == "survival") {
plot(c(ifelse(actual.scale == FALSE, 0, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, 1, sdata[nbins, "BMV.Center"])),
c(yrange.lower, yrange.upper),
type = "n",
log = "y",
xaxt = "n",
xlim = c(ifelse(actual.scale == FALSE, min.pt, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, max.pt, sdata[nbins, "BMV.Center"])),
xlab = xlabel,
ylab = ylabel,
las = 2,
main = plot.title,
sub = sub.title)
h_ <- 1
text_ <- round(effect.ac, 2)
} else if (outcome.class == "continuous") {
plot(c(ifelse(actual.scale == FALSE, 0, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, 1, sdata[nbins, "BMV.Center"])),
c(yrange.lower, yrange.upper),
type = "n",
xaxt = "n",
xlim = c(ifelse(actual.scale == FALSE, min.pt, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, max.pt, sdata[nbins, "BMV.Center"])),
xlab = xlabel,
ylab = ylabel,
main = plot.title,
sub = sub.title)
h_ <- 0
text_ <- round(effect.ac, 1)
}
#title(line = 0, xlab = xlabel)
# Reference lines for no effect/all comers result
if (show.refline) {
abline(h = h_, col = refline.color)
}
if (show.refline.ac) {
abline(h = effect.ac, col = refline.color.ac)
mtext(text = text_,
side = 4,
line = 0,
at = effect.ac,
las = 1)
}
# X-axis ticks for actual scale
if (actual.scale == TRUE) {
xticks2 <- sdata[sdata[, "Center.pt"] %in% xticks, "BMV.Center"]
old.xticks = xticks
xticks = xticks2
}
# Annotate windows
abline(v = xticks, col = refline.color, lty = 2)
if (actual.scale == TRUE) {
xticks = old.xticks
}
# ticks on x-axis: percentiles
lefts <- paste((xticks - window.width/2)*100, "%", sep = "")
rights <- paste((xticks + window.width/2)*100, "%", sep = "")
# ticks on x-axis: biomarker values
bmlefts <- bml[match(as.numeric(xticks), as.numeric(center.pt))]
bmrights <- bmu[match(as.numeric(xticks), as.numeric(center.pt))]
if (actual.scale == TRUE) {
xticks = xticks2
#mtext(paste("Median value:", xticks[1], sep = ""), side = 1, line = 0.2, adj=1, at = xticks[1], cex = 0.6)
#for (i in 2:length(xticks)) {
# mtext(xticks[i], side = 1, line = 0.2, at = xticks[i], cex = 0.6)
#}
mtext("Median value: ", side = 1, line = 0.2, adj=1, at = xticks[1], cex = 0.6)
mtext(xticks, side = 1, line = 0.2, at = xticks, cex = 0.6)
}
if (equal.in.LL == TRUE & equal.in.UL == FALSE) {
for (i in 1:(length(xticks)-1)) {
axis(1, xticks[i], paste("[", lefts[i], ", ", rights[i], ")", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts[i], ", ", bmrights[i], ")", sep = ""), side = 1, line = 1.7, at = xticks[i], cex = 0.6)
}
axis(1, xticks[length(xticks)], paste("[", lefts[length(xticks)], ", ", rights[length(xticks)], "]", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts[length(xticks)], ", ", bmrights[length(xticks)], "]", sep = ""), side = 1, line = 1.7, at = xticks[length(xticks)], cex = 0.6)
} else if (equal.in.LL == TRUE & equal.in.UL == TRUE) {
axis(1, xticks, paste("[", lefts, ", ", rights, "]", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts, ", ", bmrights, "]", sep = ""), side = 1, line = 1.7, at = xticks, cex = 0.6)
} else if (equal.in.LL == FALSE & equal.in.UL == TRUE) {
axis(1, xticks[1], paste("[", lefts[1], ", ", rights[1], "]", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts[1], ", ", bmrights[1], "]", sep = ""), side = 1, line = 1.7, at = xticks[1], cex = 0.6)
for (i in 2:length(xticks)) {
axis(1, xticks[i], paste("(", lefts[i], ", ", rights[i], "]", sep = ""), cex.axis = 0.6)
mtext(paste("(", bmlefts[i], ", ", bmrights[i], "]", sep = ""), side = 1, line = 1.7, at = xticks[i], cex = 0.6)
}
} else {
axis(1, xticks[1], paste("[", lefts[1], ", ", rights[1], ")", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts[1], ", ", bmrights[1], ")", sep = ""), side = 1, line = 1.7, at = xticks[1], cex = 0.6)
for (i in 2:(length(xticks)-1)) {
axis(1, xticks[i], paste("(", lefts[i], ", ", rights[i], ")", sep = ""), cex.axis = 0.6)
mtext(paste("(", bmlefts[i], ", ", bmrights[i], ")", sep = ""), side = 1, line = 1.7, at = xticks[i], cex = 0.6)
}
axis(1, xticks[length(xticks)], paste("(", lefts[length(xticks)], ", ", rights[length(xticks)], "]", sep = ""), cex.axis = 0.6)
mtext(paste("(", bmlefts[length(xticks)], ", ", bmrights[length(xticks)], "]", sep = ""), side = 1, line = 1.7, at = xticks[length(xticks)], cex = 0.6)
}
if (actual.scale == TRUE) {
xticks = old.xticks
}
# add events and patients
num.pats.rights <- num.pats[match(as.numeric(xticks), as.numeric(center.pt))]
if (outcome.class == "survival") {
num.events.lefts <- num.events[match(as.numeric(xticks), as.numeric(center.pt))]
if (actual.scale == TRUE) {
xticks = xticks2
}
mtext("Events/N: ", side = 1, line = 2.35, at = xticks[1], adj = 1, cex = 0.6)
mtext(paste(num.events.lefts, "/", num.pats.rights, sep = ""), side = 1, line = 2.35, at = xticks, cex = 0.6)
} else {
if (actual.scale == TRUE) {
xticks = xticks2
}
mtext(paste("N: ", num.pats.rights, sep = ""), side = 1, line = 2.35, at = xticks, cex = 0.6)
}
if (actual.scale == TRUE) {
xticks = old.xticks
}
# Wrap around with a box
box()
# Define center points as median values
if (actual.scale == TRUE) {
center.pt <- sdata[, "BMV.Center"]
}
# Actual STEPP lines are drawn at last on top of all the reference lines
lines(center.pt, effect.size, lwd = 2, col = estimate.color)
lines(center.pt, lower, lty = 2, col = ci.color)
lines(center.pt, upper, lty = 2, col = ci.color)
# Add shading
if (ci.shade) {
x <- c(center.pt[1], center.pt)
xp <- c(x, rev(x))
yp <- c(upper[1], upper, rev(lower), lower[length(lower)])
polygon(xp, yp, density=NULL, col= rgb(0.6,0.6,0.6,0.1), lwd=0.5, border=NA)
}
# Add legend
if (show.legend) {
# Display All Comers reference line in legend
if (show.refline.ac) {
legend.text <- c(legend.text, "Biomarker Population")
legend.col <- c(legend.col, refline.color.ac)
legend.lty <- c(legend.lty, 1)
legend.lwd <- c(legend.lwd, 1)
}
legend(legend.loc, legend = legend.text, col = legend.col, lty = legend.lty, lwd = legend.lwd, bty = legend.bty)
}
PlotParam()
return(sdata2)
}
lengning/gClinBiomarker documentation built on May 9, 2019, 2:55 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' Create a STEPP (Subpopulation Treatment Effect Pattern Plot)
#'
#' This function creates a STEPP from the given point estimates and confidence intervals at desired percentiles.
#'
#' @author Alexey Pronin \email{pronin.alexey@gene.com}, Ning Leng \email{leng.ning@gene.com}, and previous team members (see DESCRIPTION)
#'
#' @param data input data frame. Rows are patients and columns are variables (e.g. demographics variables, time to event variables,
#' biomarker variables, treatment indicator, etc.). One patient per row.
#' @param outcome.var outcome variable. In case of a 'survival', variable, it will be a vector of two variables: 1) time to event 2) censorship
#' @param outcome.class outcome class of the 'outcome' variable. Can be either "continuous", "binary", or "survival".
#' @param trt name of the treatment variable.
#' @param var name of the biomarker variable.
#' @param covariate vector specifying the covariate variables.
#' This can be added to adjust for in the analysis for survival and continuous outcome variable classes.
#' @param strata vector specifying the stratification variables. This can be added for the survival outcome variable class.
#' @param placebo.code name of the control group within the treatment variable
#' @param active.code name of the treatment/experimental group within the treatment variable
#' @param quantile.type an integer between 1 and 9 selecting one of the nine quantile algorithms. See \code{\link{quantile}}. Default is 2.
#' @param alpha confidence level (CI) for point estimate, i.e. 0.05 for 95 percent CI. Default is 0.05.
#' @param min.pt minimum center.pt. Default is NULL.
#' @param max.pt maximum center.pt. Default is NULL.
#' @param window.width width of each window. Default is 0.25.
#' @param by size of 'slide', i.e. speed of the window that moves along the x-axis. Default is 0.05.
#' @param yrange.lower value of the lower y-axis. Default is NULL.
#' @param yrange.upper value of the upper y-axis. Default is NULL.
#' @param xticks x-tick marks. Default is NULL.
#' @param show.refline if TRUE, the reference line will be displayed. Default is TRUE.
#' @param refline.color color of the reference line. Default is "grey".
#' @param show.refline.ac if TRUE, the reference line for effect in All Comers will be displayed. Default is TRUE.
#' @param refline.color.ac color of the reference line for effect in All Comers. Default is "lightblue".
#' @param estimate.color color of the estimate line. Default is "blue".
#' @param estimate.lty type of the estimate line. Default is 1.
#' @param estimate.lwd width of the estimate line. Default is 2.
#' @param surv.conf.type confidence interval type. Default is "plain". see conf.type in survfit.
#' @param ci.color color of the CI lines. Default is "black".
#' @param ci.lty type of the CI lines. Default is 2.
#' @param ci.lwd width of the CI lines. Default is 1.
#' @param ci.shade if TRUE, the area between the lower and upper CI lines will be shaded in gray. Default is TRUE.
#' @param plot.title title of the plot. Default is "STEPP: Subgroup Treatment Effect Pattern Plot".
#' @param sub.title subtitle of the plot, this will be displayed on the bottom of the plot. Default is NULL.
#' @param xlabel x-axis label. Default is "Biomarker Percentile".
#' @param ylabel y-axis label. Default is NULL.
#' @param show.legend if TRUE, a legend will be displayed. Default is TRUE.
#' @param legend.loc location of the legend to be displayed. Default is "topright".
#' @param legend.text text to be displayed in the legend.
#' Default is \code{c("Point Estimate",paste((1-alpha)*100," percent Confidence Interval",sep=""))}
#' @param legend.col colors of the legend. Default is \code{c(estimate.color, ci.color)}.
#' @param legend.lty type of the legend lines. Default is \code{c(estimate.lty, ci.lty)}.
#' @param legend.lwd width of the legend lines. Default is \code{c(estimate.lwd, ci.lwd)}.
#' @param legend.bty type of box for the legend. Default is "n".
#' @param bm.digits digits to be displayed/used for the lower and upper confidence level estimates. Default is 2.
#' @param actual.scale if TRUE, it generates the figure using the actual scale instead of the scale in percentage. Default is FALSE.
#' @param equal.in.LL,equal.in.UL if both are TRUE, window is defined using >= and <= (in legend: "\code{[ ]}").
#' if (TRUE, FALSE), window is defined using >= and < (in legend: "\code{[ )}").
#' if (FALSE, TRUE), window is defined using > and <= (in legend: "\code{( ]}").
#' if both are FALSE, window is defined using > and < (in legend: "\code{()}"). Default is (TRUE, FALSE).
#' @param csv.name csv file name (includes numbers used in the graphs). If NULL (default), the function will return the result.
#' @param pdf.name name of output pdf file. If it's NULL (default), the plots will be displayed but not saved as pdf.
#' @param pdf.param a list of parameters that define pdf graphics device. See \code{\link{pdf}}. Default is \code{list(width=11, height=8.5)}.
#' @param par.param a list of parameters that define graphcial parameters. See \code{\link{par}}. Default is \code{list(mar=c(4,4,3,2))}.
#' @param ties Default is "efron". To match internal sas results, use "exact". See parameter "ties" in coxph.
#' @note Patient data without corresponding biomarker data are automatically removed.
#' For survival data, censorship variable is 1 if an event happened, 0 if censored.
#'
#' @importFrom survival coxph Surv strata
#'
#' @examples
#'
#' data(input)
#' PlotSTEPP(data = input,
#' outcome.var = c("PFS", "PFS.event"),
#' outcome.class = "survival",
#' trt = "Arm",
#' var = "KRAS.exprs",
#' covariate = "Sex",
#' strata = "Age",
#' placebo.code = "CTRL",
#' active.code = "TRT")
#'
#' @seealso
#' \link{SummaryTwoGroups}
#'
#' @export
PlotSTEPP <- function(data,
outcome.var,
outcome.class,
trt=NULL,
var,
covariate = NULL,
strata = NULL,
placebo.code = NULL,
active.code = NULL,
quantile.type = 1,
alpha = .05,
window.width = .25,
min.pt = NULL,
max.pt = NULL,
by = .05,
yrange.lower=NULL,
yrange.upper=NULL,
xticks = NULL,
show.refline = TRUE,
refline.color = "grey",
show.refline.ac = TRUE,
refline.color.ac = "lightblue",
estimate.color="blue",
estimate.lty = 1,
estimate.lwd = 2,
surv.conf.type="plain",
ties="efron",
ci.color = "black",
ci.lty = 2,
ci.lwd = 1,
ci.shade = TRUE,
plot.title = "STEPP: Subgroup Treatment Effect Pattern Plot",
sub.title = NULL,
xlabel = "Biomarker Percentile",
ylabel = NULL,
show.legend = TRUE,
legend.loc = "topright",
legend.text = c("Point Estimate", paste((1-alpha)*100," percent Confidence Interval",sep="")),
legend.col = c(estimate.color, ci.color),
legend.lty = c(estimate.lty, ci.lty),
legend.lwd = c(estimate.lwd, ci.lwd),
legend.bty = "n",
bm.digits = 2,
actual.scale = FALSE,
equal.in.LL = TRUE,
equal.in.UL = FALSE,
pdf.name = NULL,
pdf.param = list(width= 11, height=8.5),
par.param = list(mar=c(4,4,3,2)),
csv.name = NULL) {
stopifnot(class(data)=="data.frame")
outcome.class <- match.arg(outcome.class, c("survival", "binary", "continuous"))
stopifnot(all(c(var, outcome.var, trt)%in%colnames(data)))
if(any(is.na(data[[var]])))message("some NA in var column, will ignore NA entries")
# Remove patient data without corresponding biomarker value
data.ori <- data
data <- data[which(!is.na(data[[var]])),]
# Read off the data
Outcome <- data[, outcome.var]
Biomarker <- data[, var]
if(!is.null(trt))if(length(unique(data[,trt]))==1) trt <- NULL
if(is.null(trt)) nArms <- 1
if (!is.null(trt)) { # multi-arm study
Treatment <- data[,trt]
Arms <- levels(factor(Treatment))
if (!is.null(placebo.code)) { # reference arm specified
if(length(placebo.code)!=1)stop("placebo.code should have length 1")
if(!placebo.code %in% Arms)stop("placebo.code should be an element in treatment column")
Arms <- c(placebo.code, setdiff(Arms, placebo.code)) # first placebo then others
if (!is.null(active.code)) {
if(length(active.code)!=1)stop("active.code should have length 1")
if(!all(active.code %in% Arms))stop("active.code should be elements in treatment column")
if(length(intersect(placebo.code, active.code))>1)
stop("code cannot be in both active.code and placebo.code!")
Arms <- c(Arms[1], active.code) #order by specified input
}
}
nArms <- length(Arms)
if(!nArms%in%c(2))stop("only 2-arm is allowed")
data[,trt] <- factor(data[,trt],levels=Arms)
if(is.null(placebo.code))placebo.code <- Arms[1]
if(is.null(active.code))active.code <- Arms[-1]
}
if (!is.null(covariate)) {
Covariate <- data[, covariate]
}
if (!is.null(strata)) {
Strat.factor <- data[, strata]
}
# If min.pt or max.pt are NULL
if (is.null(min.pt)) {
min.pt <- window.width/2
}
if (is.null(max.pt)) {
max.pt <- 1 - window.width/2
}
nbins <- round((max.pt - min.pt)/by) + 1
# X-axis ticks (if not specified by user)
if (is.null(xticks)) {
if (nbins > 4) {
allticks <- seq(min.pt, max.pt, by = by)
if (nbins %% 2 == 0) {
xticks <- c(min.pt, quantile(allticks, c(.33, .67), type = quantile.type), max.pt)
} else {
xticks <- c(min.pt, quantile(allticks, c(.25, .5, .75), type = quantile.type), max.pt)
}
} else {
xticks <- seq(min.pt, max.pt, by = by)
}
}
# Calculate summary statistics for each bin
cn <- c("Center.pt", "Effect.Size", "Lower", "Upper", "BMV.LL", "BMV.UL", "BMV.Center", "Left.pt", "Right.pt", "N")
if (outcome.class == "survival") {
cn <- c(cn, "Events")
}
sdata <- matrix(NA, nbins, length(cn))
colnames(sdata) <- cn
sdata[, "Center.pt"] <- seq(min.pt, max.pt, by)
# Effect size in All Comers
if (is.null(covariate) & is.null(strata)) {
effect.ac <- SummaryTwoGroups(Outcome, 1:length(Biomarker), Treatment,
placebo.code, active.code, outcome.class,
alpha, surv.conf.type=surv.conf.type, ties=ties,
covariate.var = NULL,
strat.factor.var = NULL)["Effect.Size"]
} else if (!is.null(covariate) & is.null(strata)) {
effect.ac <- SummaryTwoGroups(Outcome, 1:length(Biomarker), Treatment,
placebo.code, active.code, outcome.class,
alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = Covariate,
strat.factor.var = NULL)["Effect.Size"]
} else if (is.null(covariate) & !is.null(strata)) {
effect.ac <- SummaryTwoGroups(Outcome, 1:length(Biomarker), Treatment,
placebo.code, active.code, outcome.class,
alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = NULL,
strat.factor.var = Strat.factor)["Effect.Size"]
} else if (!is.null(covariate) & !is.null(strata)) {
effect.ac <- SummaryTwoGroups(Outcome, 1:length(Biomarker), Treatment,
placebo.code, active.code, outcome.class,
alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = Covariate,
strat.factor.var = Strat.factor)["Effect.Size"]
}
options(scipen=999)
for (i in 1:nbins) {
# Define this window
START_SEQ <- seq( 0 , 1 - min.pt*2 , by = by)
start <- START_SEQ[[i]]
end <- start + window.width
# Save in sdata
sdata[i, "Left.pt"] <- start
sdata[i, "Right.pt"] <- end
# Biomarker values corresponding to the percentiles
# 'type'=2 is used to match the result from median()
LL <- quantile(Biomarker, start, type = quantile.type)
UL <- quantile(Biomarker, end, type = quantile.type)
# Round LL and UL
LL <- round(LL, bm.digits)
UL <- round(UL, bm.digits)
if (LL == 0 & UL == 0) {
stop("No enough unique values, please consider to reduce number of windows (increase window.width) or increase number of digits when rounding (bm.digits)")
}
sdata[i, c("BMV.LL", "BMV.UL")] <- c(LL, UL)
# The subgroup index for this window
if (equal.in.LL == TRUE & equal.in.UL == FALSE) {
sindex <- which(Biomarker >= LL & Biomarker < UL)
} else if (equal.in.LL == TRUE & equal.in.UL == TRUE) {
sindex <- which(Biomarker >= LL & Biomarker <= UL)
} else if (equal.in.LL == FALSE & equal.in.UL == TRUE) {
sindex <- which(Biomarker > LL & Biomarker <= UL)
} else if (equal.in.LL == FALSE & equal.in.UL == FALSE) {
sindex <- which(Biomarker > LL & Biomarker < UL)
} else {
stop("equal.in.LL and equal.in.UL can be equal to TRUE or FALSE only!")
}
# Calculate median for each bean for Biomarker
sdata[i, "BMV.Center"] <- median(Biomarker[sindex])
# Sample size in each bin
sdata[i, "N"] <- length(sindex)
if (outcome.class == "survival") {
sdata[i, "Events"] <- sum(Outcome[sindex, 2], na.rm = TRUE)
}
if (is.null(covariate) & is.null(strata)) {
sdata[i, c("Effect.Size", "Lower", "Upper")] <-
SummaryTwoGroups(Outcome, sindex, Treatment,
placebo.code, active.code,
outcome.class, alpha,
covariate.var = NULL,
strat.factor.var = NULL)[c("Effect.Size", "Lower", "Upper")]
} else if (!is.null(covariate) & is.null(strata)) {
sdata[i, c("Effect.Size", "Lower", "Upper")] <-
SummaryTwoGroups(Outcome, sindex, Treatment,
placebo.code, active.code,
outcome.class, alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = Covariate,
strat.factor.var = NULL)[c("Effect.Size", "Lower", "Upper")]
} else if (is.null(covariate) & !is.null(strata)) {
sdata[i, c("Effect.Size", "Lower", "Upper")] <-
SummaryTwoGroups(Outcome, sindex, Treatment,
placebo.code, active.code,
outcome.class, alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = NULL,
strat.factor.var = Strat.factor)[c("Effect.Size", "Lower", "Upper")]
} else if (!is.null(covariate) & !is.null(strata)) {
sdata[i, c("Effect.Size", "Lower", "Upper")] <-
SummaryTwoGroups(Outcome, sindex, Treatment,
placebo.code, active.code,
outcome.class, alpha,surv.conf.type=surv.conf.type, ties=ties,
covariate.var = Covariate,
strat.factor.var = Strat.factor)[c("Effect.Size", "Lower", "Upper")]
}
} # end of for (i in 1:nbins)
##### Print out sdata
sdata2 <- sdata
colnames(sdata2)[1:9] <- c("Window Center", "Effect.Size", "CI Lower", "CI Upper",
"BM Lower", "BM Upper", "BM Center", "Window Left", "Window Right")
if (!is.null(ylabel)) {
colnames(sdata2) <- gsub("Effect.Size", ylabel, colnames(sdata2))
} else {
if (outcome.class == "survival") {
colnames(sdata2) <- gsub("Effect.Size", "Hazard Ratio", colnames(sdata2))
} else if (outcome.class == "binary") {
colnames(sdata2) <- gsub("Effect.Size", "Proportion Difference", colnames(sdata2))
} else if (outcome.class == "continuous") {
colnames(sdata2) <- gsub("Effect.Size", "Mean Difference", colnames(sdata2))
}
}
if (!is.null(csv.name)) {
write.csv(sdata2, file = csv.name, row.names = FALSE)
}
##### Draw the plot
PlotParam(pdf.name, pdf.param, par.param)
if(is.null(sub.title))sub.title <- paste("soc:", placebo.code, "; trt:", active.code)
center.pt <- sdata[, "Center.pt"]
effect.size <- sdata[, "Effect.Size"]
lower = sdata[, "Lower"]
upper = sdata[, "Upper"]
bml = sdata[, "BMV.LL"]
bmu = sdata[, "BMV.UL"]
if (outcome.class == "survival") {
num.events = sdata[, "Events"]
}
num.pats = sdata[, "N"]
if (actual.scale == FALSE) {
# If center.pt is NULL
if (is.null(center.pt)) {
center.pt <- seq(min.pt, max.pt, by = by)
}
# x values should always be between 0 and 1
if (sum(center.pt <= 0 | center.pt >= 1) > 0) {
warning("Center percentile must be between 0 and 1 !")
}
}
# Starting point should be smaller than Ending point
if (min.pt < 0 | min.pt > max.pt | max.pt > 1 ) {
warning("Inappropriate center percentile; Please check.")
print(paste("window.width", window.width))
print(paste("min.pt", min.pt))
print(paste("max.pt", max.pt))
}
if (is.null(yrange.lower)) {
yrange.lower= min(lower, na.rm=TRUE)
}
if (is.null(yrange.upper)) {
yrange.upper= max(upper, na.rm=TRUE)
}
if (outcome.class == "binary") {
if(is.null(ylabel)) {
ylabel <- "Proportion Difference (%)"
}
if (yrange.lower < -1 | yrange.upper > 1) {
warning("Input yrange.lower and yrange.upper should be between -1 and 1 !")
}
} else if (outcome.class == "survival") {
if (is.null(ylabel)) {
ylabel <- "Hazard Ratio"
}
if (yrange.lower < 0 | yrange.upper < 0) {
warning("Input yrange.lower and yrange.upper should be >0 !")
}
} else if (outcome.class == "continuous") {
if (is.null(ylabel)) {
ylabel <- "Mean Difference"
}
}
# Set an empty panel
if (outcome.class == "binary") {
plot(c(ifelse(actual.scale == FALSE, 0, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, 1, sdata[nbins, "BMV.Center"])),
c(yrange.lower, yrange.upper),
type = "n",
axes = FALSE,
xlim = c(ifelse(actual.scale == FALSE, min.pt, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, max.pt, sdata[nbins, "BMV.Center"])),
xlab = xlabel,
ylab = ylabel,
main = plot.title,
sub = sub.title)
# y-axis in % scale
axis(2,
-10:10/10,
paste(c(rep("", 11), rep("+", 10)), -10:10*10, sep = ""),
las = 2)
h_ <- 0
text_ <- paste(round(effect.ac*100), "%", sep = "")
} else if (outcome.class == "survival") {
plot(c(ifelse(actual.scale == FALSE, 0, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, 1, sdata[nbins, "BMV.Center"])),
c(yrange.lower, yrange.upper),
type = "n",
log = "y",
xaxt = "n",
xlim = c(ifelse(actual.scale == FALSE, min.pt, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, max.pt, sdata[nbins, "BMV.Center"])),
xlab = xlabel,
ylab = ylabel,
las = 2,
main = plot.title,
sub = sub.title)
h_ <- 1
text_ <- round(effect.ac, 2)
} else if (outcome.class == "continuous") {
plot(c(ifelse(actual.scale == FALSE, 0, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, 1, sdata[nbins, "BMV.Center"])),
c(yrange.lower, yrange.upper),
type = "n",
xaxt = "n",
xlim = c(ifelse(actual.scale == FALSE, min.pt, sdata[1, "BMV.Center"]), ifelse(actual.scale == FALSE, max.pt, sdata[nbins, "BMV.Center"])),
xlab = xlabel,
ylab = ylabel,
main = plot.title,
sub = sub.title)
h_ <- 0
text_ <- round(effect.ac, 1)
}
#title(line = 0, xlab = xlabel)
# Reference lines for no effect/all comers result
if (show.refline) {
abline(h = h_, col = refline.color)
}
if (show.refline.ac) {
abline(h = effect.ac, col = refline.color.ac)
mtext(text = text_,
side = 4,
line = 0,
at = effect.ac,
las = 1)
}
# X-axis ticks for actual scale
if (actual.scale == TRUE) {
xticks2 <- sdata[sdata[, "Center.pt"] %in% xticks, "BMV.Center"]
old.xticks = xticks
xticks = xticks2
}
# Annotate windows
abline(v = xticks, col = refline.color, lty = 2)
if (actual.scale == TRUE) {
xticks = old.xticks
}
# ticks on x-axis: percentiles
lefts <- paste((xticks - window.width/2)*100, "%", sep = "")
rights <- paste((xticks + window.width/2)*100, "%", sep = "")
# ticks on x-axis: biomarker values
bmlefts <- bml[match(as.numeric(xticks), as.numeric(center.pt))]
bmrights <- bmu[match(as.numeric(xticks), as.numeric(center.pt))]
if (actual.scale == TRUE) {
xticks = xticks2
#mtext(paste("Median value:", xticks[1], sep = ""), side = 1, line = 0.2, adj=1, at = xticks[1], cex = 0.6)
#for (i in 2:length(xticks)) {
# mtext(xticks[i], side = 1, line = 0.2, at = xticks[i], cex = 0.6)
#}
mtext("Median value: ", side = 1, line = 0.2, adj=1, at = xticks[1], cex = 0.6)
mtext(xticks, side = 1, line = 0.2, at = xticks, cex = 0.6)
}
if (equal.in.LL == TRUE & equal.in.UL == FALSE) {
for (i in 1:(length(xticks)-1)) {
axis(1, xticks[i], paste("[", lefts[i], ", ", rights[i], ")", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts[i], ", ", bmrights[i], ")", sep = ""), side = 1, line = 1.7, at = xticks[i], cex = 0.6)
}
axis(1, xticks[length(xticks)], paste("[", lefts[length(xticks)], ", ", rights[length(xticks)], "]", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts[length(xticks)], ", ", bmrights[length(xticks)], "]", sep = ""), side = 1, line = 1.7, at = xticks[length(xticks)], cex = 0.6)
} else if (equal.in.LL == TRUE & equal.in.UL == TRUE) {
axis(1, xticks, paste("[", lefts, ", ", rights, "]", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts, ", ", bmrights, "]", sep = ""), side = 1, line = 1.7, at = xticks, cex = 0.6)
} else if (equal.in.LL == FALSE & equal.in.UL == TRUE) {
axis(1, xticks[1], paste("[", lefts[1], ", ", rights[1], "]", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts[1], ", ", bmrights[1], "]", sep = ""), side = 1, line = 1.7, at = xticks[1], cex = 0.6)
for (i in 2:length(xticks)) {
axis(1, xticks[i], paste("(", lefts[i], ", ", rights[i], "]", sep = ""), cex.axis = 0.6)
mtext(paste("(", bmlefts[i], ", ", bmrights[i], "]", sep = ""), side = 1, line = 1.7, at = xticks[i], cex = 0.6)
}
} else {
axis(1, xticks[1], paste("[", lefts[1], ", ", rights[1], ")", sep = ""), cex.axis = 0.6)
mtext(paste("[", bmlefts[1], ", ", bmrights[1], ")", sep = ""), side = 1, line = 1.7, at = xticks[1], cex = 0.6)
for (i in 2:(length(xticks)-1)) {
axis(1, xticks[i], paste("(", lefts[i], ", ", rights[i], ")", sep = ""), cex.axis = 0.6)
mtext(paste("(", bmlefts[i], ", ", bmrights[i], ")", sep = ""), side = 1, line = 1.7, at = xticks[i], cex = 0.6)
}
axis(1, xticks[length(xticks)], paste("(", lefts[length(xticks)], ", ", rights[length(xticks)], "]", sep = ""), cex.axis = 0.6)
mtext(paste("(", bmlefts[length(xticks)], ", ", bmrights[length(xticks)], "]", sep = ""), side = 1, line = 1.7, at = xticks[length(xticks)], cex = 0.6)
}
if (actual.scale == TRUE) {
xticks = old.xticks
}
# add events and patients
num.pats.rights <- num.pats[match(as.numeric(xticks), as.numeric(center.pt))]
if (outcome.class == "survival") {
num.events.lefts <- num.events[match(as.numeric(xticks), as.numeric(center.pt))]
if (actual.scale == TRUE) {
xticks = xticks2
}
mtext("Events/N: ", side = 1, line = 2.35, at = xticks[1], adj = 1, cex = 0.6)
mtext(paste(num.events.lefts, "/", num.pats.rights, sep = ""), side = 1, line = 2.35, at = xticks, cex = 0.6)
} else {
if (actual.scale == TRUE) {
xticks = xticks2
}
mtext(paste("N: ", num.pats.rights, sep = ""), side = 1, line = 2.35, at = xticks, cex = 0.6)
}
if (actual.scale == TRUE) {
xticks = old.xticks
}
# Wrap around with a box
box()
# Define center points as median values
if (actual.scale == TRUE) {
center.pt <- sdata[, "BMV.Center"]
}
# Actual STEPP lines are drawn at last on top of all the reference lines
lines(center.pt, effect.size, lwd = 2, col = estimate.color)
lines(center.pt, lower, lty = 2, col = ci.color)
lines(center.pt, upper, lty = 2, col = ci.color)
# Add shading
if (ci.shade) {
x <- c(center.pt[1], center.pt)
xp <- c(x, rev(x))
yp <- c(upper[1], upper, rev(lower), lower[length(lower)])
polygon(xp, yp, density=NULL, col= rgb(0.6,0.6,0.6,0.1), lwd=0.5, border=NA)
}
# Add legend
if (show.legend) {
# Display All Comers reference line in legend
if (show.refline.ac) {
legend.text <- c(legend.text, "Biomarker Population")
legend.col <- c(legend.col, refline.color.ac)
legend.lty <- c(legend.lty, 1)
legend.lwd <- c(legend.lwd, 1)
}
legend(legend.loc, legend = legend.text, col = legend.col, lty = legend.lty, lwd = legend.lwd, bty = legend.bty)
}
PlotParam()
return(sdata2)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.