# NOT RUN {
#Simulating data from four gaussian distributions,
#with mean risks equal to 0.5 in each arm:
x0E <- rnorm(250) # reference arm, event
x0Eb <- rnorm(250, 2) # reference arm, no event
x1E <- rnorm(250, 2) # innovative arm, event
x1Eb <- rnorm(250) # innovative arm, no event
#When working with real data. You can check the randomization constraint using the
#densCurves function:
densCurves(x0 = c(x0E, x0Eb), x1 = c(x1E, x1Eb), type = "treatment selection")
#You can also use the riskCurves function to know if low values of the marker are associated
#with a better response under the reference treatment or not:
library(mgcv)
riskCurves(x0E, x0Eb, x1E, x1Eb)
#Fit normal distributions on three groups. And let the last one (1E) be undefined (derived
#indirectly using the randomization constraint):
fit0E <- fit(x0E, "norm")
fit0Eb <- fit(x0Eb, "norm")
fit1E <- fit(x1E, "undefined")
fit1Eb <- fit(x1Eb, "norm")
#Apply the main function to estimate the optimal threshold:
# first case: the mean risks of event in the two treatment arms are left unspecified (are
# determined by the number of marker measurements in the fit0E, fi0Eb, fit1E, fit1Eb)
# }
# NOT RUN {
res <- trtSelThresh(fit0E, fit0Eb, fit1E, fit1Eb,
lowRef = FALSE, toxRef = FALSE, r = 0.02, le.MCMC = 5000, plot = TRUE,
progress.bar = "text")
# second case: the mean risks of event in the two treatment arms are given through mcmc.lists
# that correspond to their posterior distributions (see the fit man page for examples on how
# to generate posterior distributions manually)
#You can summarize the results using the summary() function:
summary(res, method = "median")
#You can extract the estimates and CI bounds of each indicator presented in the summary:
estimates(res, method = "median")
credibleIntervals(res)
#Plot the decision curves (this function is time-consuming):
dCres <- decisionCurve(res, r = seq(0, 0.2, length.out = 6))
#You can plot again the decision curves by applying the plot method to dCres,
#this function is a lot faster than the previous one. It also has more options
#to customize the plot:
plot(dCres)
plot(dCres, which = 1)
plot(dCres, which = 2)
# }
Run the code above in your browser using DataLab