myplotcompare.trtsel <-
function(x, bootstraps = 500, alpha = .05,
ci = "horizontal", marker.names = c("Marker 1", "Marker 2"),
fixeddeltas.y1 = NULL, fixeddeltas.y2,
xlab = NULL,
ylab = NULL,
xlim = NULL,
ylim = NULL,
main = NULL, offset = offset, conf.bands)
{
ts1 <- x$x
ts2 <- x$x2
if(substr(ci, 1, 4) =="hori") {
## trt selcurve plot with horizontal ci bands
fix.ind = 2 # fix delta
out.ind = 1 # vary F.marker
}else if(substr(ci, 1, 4) =="vert"){
## trt sel curve plot with vertical ci bands
fix.ind = 1 #fix F.event
out.ind = 2 #vary delta
}
fittedrisk.t0.y1 <- ts1$derived.data$fittedrisk.t0
fittedrisk.t1.y1 <- ts1$derived.data$fittedrisk.t1
delta.y1 <- ts1$derived.data$trt.effect
link <- ts1$model.fit$family
fittedrisk.t0.y2 <- ts2$derived.data$fittedrisk.t0
fittedrisk.t1.y2 <- ts2$derived.data$fittedrisk.t1
delta.y2 <- ts2$derived.data$trt.effect
rho <- ts1$model.fit$cohort.attributes
study.design <- ts1$model.fit$study.design
trt <- ts1$derived.data[[ts1$treatment.name]]
if( link$family == "time-to-event"){
event.name1 = ts1$formula[[2]]
event.name2 = ts2$formula[[2]]
}else{
event.name1 = as.character(ts1$formula[[2]])
event.name2 = as.character(ts2$formula[[2]])
}
boot.sample <- ts1$functions$boot.sample
get.F <- ts1$functions$get.F
###
#old.par <- par(no.readonly = TRUE)
if( length(fixeddeltas.y1) > 0 & bootstraps >0 ) {
#bootstrapping done by fixing response and strapping marker
# if(is.null(fixed.values) & ci =="horizontal"){
# fixeddeltas.y1 <- seq(from = min(delta.y1), to = max(delta.y1), length.out = 100)
# fixeddeltas.y2 <- seq(from = min(delta.y2), to = max(delta.y2), length.out = 100)
# }else if(is.null(fixed.values) & ci =="vertical"){
# fixeddeltas.y1 <- 1:100/100
# fixeddeltas.y2 <- 1:100/100
# }else{
#fixeddeltas.y1 <- fixed.values
#fixeddeltas.y2 <- fixed.values
# }
if(link$family == "risks_provided"){
provided_risk <- cbind(fittedrisk.t0.y1,
fittedrisk.t1.y1,
fittedrisk.t0.y2,
fittedrisk.t1.y2)
}else{
provided_risk = NULL
}
boot.dat <- replicate( bootstraps, one.boot.plot.compare(data1 = ts1$derived.data, data2 = ts2$derived.data,
formulas = list(ts1$formula, ts2$formula),
event.names = c(event.name1, event.name2),
treatment.names = c(ts1$treatment.name, ts2$treatment.name),
ci = ci,
fixeddeltas.y1 = fixeddeltas.y1, fixeddeltas.y2 = fixeddeltas.y2,
rho = rho, study.design = study.design, obp.boot.sample = boot.sample, obp.get.F = get.F, fix.ind, out.ind, link = link,
provided_risk = provided_risk,
prediction.times = c(x$x$prediction.time, x$x2$prediction.time)))
if(length(fixeddeltas.y1)==1){
bounds.delta.y1<- quantile(boot.dat[1,,], probs = c(alpha/2, 1-alpha/2), na.rm = TRUE)
bounds.delta.y2<- quantile(boot.dat[2,,], probs = c(alpha/2, 1-alpha/2), na.rm = TRUE)
bounds.delta.y1 = t(t(bounds.delta.y1))
bounds.delta.y2 = t(t(bounds.delta.y2))
}else{
bounds.delta.y1<- apply(boot.dat[1,,], 1, function(x, ...){quantile(unlist(x), ...)}, probs = c(alpha/2, 1-alpha/2), na.rm = TRUE)
bounds.delta.y2<- apply(boot.dat[2,,], 1, function(x, ...){quantile(unlist(x), ...)}, probs = c(alpha/2, 1-alpha/2), na.rm = TRUE)
}
}else{
bounds.delta.y1 <- NULL
bounds.delta.y2 <- NULL
}
## end bootstrap
#set up the ylimit of the treatment effect plot
if(is.null(ylim)){
if(substr(ci, 1,1) %in% c("v", "V")){
min.delta <- min(c(delta.y1, delta.y2, bounds.delta.y1, bounds.delta.y2), na.rm=TRUE)
max.delta <- max(c(delta.y1, delta.y2, bounds.delta.y1, bounds.delta.y2), na.rm=TRUE)
cen <- mean(c(min.delta, max.delta), na.rm=TRUE)
}else{
min.delta <- min(c(delta.y1, delta.y2), na.rm=TRUE)
max.delta <- max(c(delta.y1, delta.y2), na.rm=TRUE)
cen <- mean(c(min.delta, max.delta), na.rm=TRUE)
}
ran <- max.delta - min.delta
ran <- ran*1.1
ylim <- c(cen-ran/2, cen+ran/2)
}
ts1.curves <- trteffectPLOTcompare_gg(x1=ts1, x2 = ts2, ci = ci, ci.bounds = rbind(bounds.delta.y1, bounds.delta.y2), get.F = get.F, fixed.values = rbind(fixeddeltas.y1, fixeddeltas.y2+offset), conf.bands = conf.bands, rho=rho, xlab=xlab, ylab = ylab, xlim=xlim, ylim = ylim, main = main)
p <- ts1.curves[[1]]
p <- p + scale_linetype_manual(name = "", breaks = c("1","2", "3", "4"), values = c(1, 2, 3, 4), labels = c(marker.names, "Mean", "Zero"))
p <- p + scale_size_manual(name = "", breaks = c("1", "2", "3", "4"), values = c(1, 1, .5, .5), labels = c(marker.names, "Mean", "Zero"))
print(p)
#if(is.null(xlim)) xlim = c(0,100)
#legend(x=xlim[2]+diff(xlim)/15, y = quantile(ylim, prob = .5), legend = marker.names, lty = c(1,2), col = c("black", "black"), bty="n", cex = 1, xpd = TRUE, lwd=c(2,2))
if(bootstraps >0 & length(fixeddeltas.y1) > 0 ){
conf.ints.y1 <- as.data.frame(cbind(fixeddeltas.y1, t(bounds.delta.y1)))
names(conf.ints.y1) <- c("fixed", "lower", "upper")
conf.ints.y2 <- as.data.frame(cbind(fixeddeltas.y2, t(bounds.delta.y2)))
names(conf.ints.y2) <- c("fixed", "lower", "upper")
result <- list("plot" = p,
"x" = list( "conf.intervals" = conf.ints.y1),
"x2" = list( "conf.intervals" = conf.ints.y2))
}else{
result <- list("plot" = p)
}
invisible(result)
#par(old.par)
}
myplotcompare.trtsel_disc <-
function(x, bootstraps = 500, alpha = .05,
ci = "horizontal", marker.names = c("Marker 1", "Marker 2"),
xlab = NULL,
ylab = NULL,
xlim = NULL,
ylim = NULL,
main = NULL, offset = offset, conf.bands, annotate.plot)
{
quantile <- NULL #appease check
ts1 <- x$x
ts2 <- x$x2
fittedrisk.t0.y1 <- ts1$derived.data$fittedrisk.t0
fittedrisk.t1.y1 <- ts1$derived.data$fittedrisk.t1
marker1 = ts1$derived.data[[ts1$model.fit$marker.names]]
delta.y1 <- ts1$derived.data$trt.effect
link <- ts1$model.fit$link
fittedrisk.t0.y2 <- ts2$derived.data$fittedrisk.t0
fittedrisk.t1.y2 <- ts2$derived.data$fittedrisk.t1
marker2 <- ts2$derived.data[[ts2$model.fit$marker.names]]
delta.y2 <- ts2$derived.data$trt.effect
rho <- ts1$model.fit$cohort.attributes
study.design <- ts1$model.fit$study.design
trt <- ts1$derived.data$trt
if( ts1$model.fit$family$family == "time-to-event"){
warning("plotting comparisons of two discrete markers with a time-to-event outcome is not implemented yet.")
return(NULL)
}else{
event <- ts1$derived.data[[as.character(ts1$formula[[2]])]]
boot.sample <- ts1$functions$boot.sample
one.boot.plot_disc <-
function(event, trt, marker1, marker2, w1, w2, rho = rho, obp.boot.sample){
myboot.sample <- obp.boot.sample( event, trt, rho)
rho.b <- myboot.sample[1:7]
ind <- myboot.sample[-c(1:7)]
event.b <- event[ind]
trt.b <- trt[ind]
marker1.b <- marker1[ind]
marker2.b <- marker2[ind]
mval1 <- sort(unique(marker1))
mval2 <- sort(unique(marker2))
c( trteff.mkr10 =mean(event.b[trt.b==0 & marker1.b ==mval1[1]]) - mean(event.b[trt.b==1 & marker1.b ==mval1[1]]),
trteff.mkr11 =mean(event.b[trt.b==0 & marker1.b ==mval1[2]]) - mean(event.b[trt.b==1 & marker1.b ==mval1[2]]),
trteff.mkr20 =mean(event.b[trt.b==0 & marker2.b ==mval2[1]]) - mean(event.b[trt.b==1 & marker2.b ==mval2[1]]),
trteff.mkr21 =mean(event.b[trt.b==0 & marker2.b ==mval2[2]]) - mean(event.b[trt.b==1 & marker2.b ==mval2[2]])
)
}
if(conf.bands){
boot.data <- replicate(bootstraps, one.boot.plot_disc( event, trt, marker1, marker2, rho,obp.boot.sample = boot.sample))
mval1 <- sort(unique(marker1))
mval2 <- sort(unique(marker2))
row.names(boot.data) = c(paste("trteffect.1mkr", mval1[1], sep = ""),
paste("trteffect.1mkr", mval1[2], sep = ""),
paste("trteffect.2mkr", mval2[1], sep = ""),
paste("trteffect.2mkr", mval2[2], sep = ""))
#horizontal
if(substr(ci, 1,1) =="h") { warning("Horizontal CI bands are not allowed for treatment effect plots with a discrete marker. Vertical bands will be computed"); ci <- "vertical";}
#vertical
myconf.ints <- apply(boot.data, 1, quantile, probs = c(alpha/2, 1-alpha/2))
ci = "vertical"
}else{
myconf.ints = NULL
}
ts1.curves <- trteffectPLOTcompare_gg_disc(x1=ts1,
x2 = ts2,
ci.bounds = myconf.ints,
conf.bands = conf.bands,
offset = offset,
xlab=xlab, ylab = ylab,
xlim=xlim, ylim = ylim,
main = main,
marker.names = marker.names,
annotate.plot = annotate.plot)
p <- ts1.curves[[1]]
print(p)
#if(is.null(xlim)) xlim = c(0,100)
#legend(x=xlim[2]+diff(xlim)/15, y = quantile(ylim, prob = .5), legend = marker.names, lty = c(1,2), col = c("black", "black"), bty="n", cex = 1, xpd = TRUE, lwd=c(2,2))
result <- list("plot" = p,
"ci.bounds" = ts1.curves[[2]])
}
invisible(result)
#par(old.par)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.