Nothing
one.boot.plot <-
function(x, ci, fixed.values, fix.ind, out.ind){
if( x$model.fit$family$family == "time-to-event"){
mysurv <- with(x$derived.data, eval(x$formula[[2]]))
event <- mysurv[,2]
}else{
event <- x$derived.data[[as.character(x$formula[[2]])]]
}
myboot.sample <- x$functions$boot.sample( event,
x$derived.data[[x$treatment.name]],
rho = x$model.fit$cohort.attributes)
#this makes it work for step function
if(substr(ci, 1,1)=="h") addind = 0
else addind = 1
rho.b <- myboot.sample[1:7]
ind <- myboot.sample[-c(1:7)]
if( x$model.fit$family$family == "time-to-event"){
event.b <- 0
}else{
event.b <- x$derived.data[[as.character(x$formula[[2]])]][ind]
}
trt.b <- x$derived.data[[x$treatment.name]][ind]
if(x$model.fit$family$family == "risks_provided"){
obsrisk.t0.b <- x$derived.data$fittedrisk.t0[ind]
obsrisk.t1.b <- x$derived.data$fittedrisk.t1[ind]
linkinvfun <- NULL
marker.b <- obsrisk.t0.b - obsrisk.t1.b#
}else if(x$model.fit$family$family == "time-to-event"){
coxfit <- do.call(coxph, list(x$formula, x$derived.data[ind,]))
obsrisk.t0.b <- get.risk.t_coxph(coxfit, x$treatment.name, x$derived.data[ind,], x$prediction.time, t = 0)
obsrisk.t1.b <- get.risk.t_coxph(coxfit, x$treatment.name, x$derived.data[ind,], x$prediction.time, t = 1)
#we still need to incorporate the nelson aalen baseline haz to get absolute risk at t = 'prediction.time'
}else{
# browser()
coef <- unname(get.coef(x$formula, x$treatment.name, x$derived.data[ind,],
x$model.fit$study.design,
rho.b,
family = x$model.fit$family)[,1])
linkinvfun <- x$model.fit$family$linkinv
obsrisk.t0.b <- get.risk.t(coef, x$formula, x$treatment.name, data = x$derived.data[ind,], linkinvfun, t = 0)
obsrisk.t1.b <- get.risk.t(coef, x$formula, x$treatment.name, data = x$derived.data[ind,], linkinvfun, t = 1)
wi = 0
}
obsdelta.b <-obsrisk.t0.b - obsrisk.t1.b#
## if there is a single marker, use that to calculate F.Y, otherwise use treatment
## effect as the marker
if(length(x$model.fit$marker.names) ==1 ){
marker.b <- x$derived.data[ind, x$model.fit$marker.names]
}else{
marker.b <- NULL
}
if(is.null(marker.b)) marker.b <- obsrisk.t0.b - obsrisk.t1.b#
F.Y <- x$functions$get.F( marker.b, event.b, trt.b, rho.b)*100#
F.D <- x$functions$get.F( obsdelta.b, event.b, trt.b, rho.b)*100#
theta.c <- EventRateVec(obsrisk.t0.b, obsrisk.t1.b, F.D, rho.b, event.b, trt.b)
#all
all <- cbind( F.Y, obsrisk.t0.b, obsrisk.t1.b, F.D, obsdelta.b, theta.c)
all <- unique(all)
# browser()
if(length(fix.ind) > 1){
myorder <- apply(all[,fix.ind], 2, order)
out <- matrix(0, ncol = length(fixed.values), nrow = length(fix.ind))
for( i in 1:length(fix.ind)){
#for decreasing pred curves with horizontal bands
if(is.element(fix.ind[i], c(2,3)) & !(all.equal(order(all[,fix.ind[i]]), order(all[,1]))==TRUE)){ addind = 1}
tmpind <- sum.I(fixed.values, ">=", all[myorder[,i],fix.ind[i]])+addind
tmpind[tmpind<=0] <- NA
tmpall <- all[myorder[,i],out.ind[i]]
out[i, ] <- tmpall[tmpind]
}
}else{
myorder <- order(all[,fix.ind])
out <- numeric( length(fixed.values))
tmpind <- sum.I(fixed.values, ">=", all[myorder,fix.ind]) +addind
tmpind[tmpind==0] <- NA
tmpall <- all[myorder,out.ind]
out <- tmpall[tmpind]
}
out
}
one.boot.plot_disc <-
function(x){
if( x$model.fit$family$family == "time-to-event"){
mysurv <- with(x$derived.data, eval(x$formula[[2]]))
event <- mysurv[,2]
}else{
event <- x$derived.data[[as.character(x$formula[[2]])]]
}
myboot.sample <- x$functions$boot.sample( event = event,
trt = x$derived.data[[x$treatment.name]],
rho = x$model.fit$cohort.attributes)
rho.b <- myboot.sample[1:7]
ind <- myboot.sample[-c(1:7)]
# browser()
#event.b <- x$derived.data[[as.character(x$formula[[2]])]][ind]
# trt.b <- x$derived.data[[x$treatment.name]][ind]
marker.b <- x$derived.data[ind, x$model.fit$marker.names]
if(x$model.fit$family$family == "risks_provided")
{
provided_risk <- cbind(x$derived.data$fittedrisk.t0[ind],
x$derived.data$fittedrisk.t1[ind])
} else provided_risk = NULL
tmp.trtsel <- trtsel.boot(x$formula, x$treatment.name,
data = x$derived.data[ind, ],
d=x$model.fit$thresh,
rho = rho.b,
study.design = x$model.fit$study.design,
family = x$model.fit$family,
disc.rec.no.trt = x$model.fit$disc.rec.no.trt,
provided_risk = provided_risk,
prediction.time = x$prediction.time)
####
unique.fitted.vals <- unique(cbind(marker.b, tmp.trtsel$derived.data[,c("fittedrisk.t0", "fittedrisk.t1", "trt.effect")]))
mval <- sort(unique(marker.b))
out <- with( unique.fitted.vals,
c( event.trt0.mkr0 = fittedrisk.t0[marker.b==mval[1]],
event.trt1.mkr0 = fittedrisk.t1[marker.b==mval[1]],
event.trt0.mkr1 = fittedrisk.t0[marker.b==mval[2]],
event.trt1.mkr1 = fittedrisk.t1[marker.b==mval[2]],
trteff.mkr0 = trt.effect[marker.b==mval[1]],
trteff.mkr1 = trt.effect[marker.b==mval[2]]))
return(out)
}
one.boot.plot.compare <-
function(data1, data2, formulas, event.names, treatment.names, ci,
fixeddeltas.y1, fixeddeltas.y2,
rho = rho, study.design, obp.boot.sample, obp.get.F, fix.ind, out.ind, link,
provided_risk = NULL,
prediction.times){
if(substr(ci, 1,1)=="h") addind = 0
else addind = 1
if( link$family == "time-to-event"){
mysurv <- with(data1, eval(event.names[[1]]))
event1 <- mysurv[,2]
mysurv <- with(data2, eval(event.names[[2]]))
event2 <- mysurv[,2]
}else{
event1 <- data1[[event.names[[1]]]]
event2 <- data2[[event.names[[2]]]]
}
myboot.sample <- obp.boot.sample( event1, data1[[treatment.names[1]]], rho)
rho.b <- myboot.sample[1:7]
ind <- myboot.sample[-c(1:7)]
event.b = event1[ind]
trt.b <- data1[[treatment.names[1]]][ind]
### marker1
fixed.values <- fixeddeltas.y1
if(link$family == "risks_provided"){
obsrisk.t0.b <- provided_risk[ind,1]
obsrisk.t1.b <- provided_risk[ind,2]
}else if(link$family == "time-to-event"){
coxfit <- do.call(coxph, list(formulas[[1]], data1[ind,]))
obsrisk.t0.b <- get.risk.t_coxph(coxfit, treatment.names[1], data1[ind,], prediction.time = prediction.times[1], t = 0)
obsrisk.t1.b <- get.risk.t_coxph(coxfit, treatment.names[1], data1[ind,], prediction.time = prediction.times[1], t = 1)
#we still need to incorporate the nelson aalen baseline haz to get absolute risk at t = 'prediction.time'
coef <- summary(coxfit)$coefficients
}else{
coef <- unname(get.coef(formulas[[1]], treatment.names[1], data1[ind,], study.design, rho.b, link)[,1])
linkinvfun <- link$linkinv
obsrisk.t0.b <- get.risk.t(coef, formulas[[1]], treatment.names[1], data1[ind,], linkinvfun, t = 0)
obsrisk.t1.b <- get.risk.t(coef, formulas[[1]], treatment.names[1], data1[ind,], linkinvfun, t = 1)
#obsrisk.t0.b <- c(ifelse(coef[3]>0, 0,1), get.risk.t0(coef, marker.b, linkinvfun))
#obsrisk.t1.b <- c(ifelse(sum(coef[3:4])>0, 0,1), get.risk.t1(coef, marker.b, linkinvfun))
}
obsdelta.b <- obsrisk.t0.b- obsrisk.t1.b
F.D <- obp.get.F( obsdelta.b, event.b, trt.b, rho.b)*100
#obsdelta.b <- c(-1, obsrisk.t0.b[-1] - obsrisk.t1.b[-1])
#F.D <- c(0, obp.get.F( obsdelta.b[-1], event.b, trt.b, rho.b))*100
all <- cbind( F.D, obsdelta.b)
myorder <- order(all[,fix.ind])
out <- numeric( length(fixed.values))
tmpind <- sum.I(fixed.values, ">=", all[myorder,fix.ind])
tmpind[tmpind==0] <- NA
tmpall <- all[myorder,out.ind]
out <- tmpall[tmpind]
marker.vals.delta.Y1 <- out
##
### Y2
fixed.values <- fixeddeltas.y2
if(link$family == "risks_provided"){
obsrisk.t0.b <- provided_risk[ind,3]
obsrisk.t1.b <- provided_risk[ind,4]
}else if(link$family == "time-to-event"){
coxfit <- do.call(coxph, list(formulas[[2]], data2[ind,]))
obsrisk.t0.b <- get.risk.t_coxph(coxfit, treatment.names[2], data2[ind,], prediction.time = prediction.times[2], t = 0)
obsrisk.t1.b <- get.risk.t_coxph(coxfit, treatment.names[2], data2[ind,], prediction.time = prediction.times[2], t = 1)
#we still need to incorporate the nelson aalen baseline haz to get absolute risk at t = 'prediction.time'
coef <- summary(coxfit)$coefficients
}else{
coef <- unname(get.coef(formulas[[2]], treatment.names[2], data2[ind,],study.design, rho.b, link)[,1])
linkinvfun <- link$linkinv
obsrisk.t0.b <- get.risk.t(coef, formulas[[2]], treatment.names[2], data2[ind,], linkinvfun, t = 0)
obsrisk.t1.b <- get.risk.t(coef, formulas[[2]], treatment.names[2], data2[ind,], linkinvfun, t = 1)
#obsrisk.t0.b <- c(ifelse(coef[3]>0, 0,1), get.risk.t0(coef, marker2.b, linkinvfun))
#obsrisk.t1.b <- c(ifelse(sum(coef[3:4])>0, 0,1), get.risk.t1(coef, marker2.b, linkinvfun))
}
obsdelta.b <- obsrisk.t0.b - obsrisk.t1.b
F.D <- obp.get.F( obsdelta.b, event.b, trt.b, rho.b)*100
#obsdelta.b <- c(-1, obsrisk.t0.b[-1] - obsrisk.t1.b[-1])
#F.D <- c(0, obp.get.F( obsdelta.b[-1], event.b, trt.b, rho.b))*100
all <- cbind( F.D, obsdelta.b)
myorder <- order(all[,fix.ind])
out <- numeric( length(fixed.values))
tmpind <- sum.I(fixed.values, ">=", all[myorder,fix.ind]) + addind
tmpind[tmpind==0] <- NA
tmpall <- all[myorder,out.ind]
out <- tmpall[tmpind]
marker.vals.delta.Y2 <- out
##
return( rbind( marker.vals.delta.Y1, marker.vals.delta.Y2))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.