# NOT RUN {
# (1) a toy example (a very small subsample of a microarray data set)
#
datafilename <- system.file("extdata", "transcripto_very_small_sample.txt",
package="DRomics")
# to test the package on a small but not very small data set
# use the following commented line
# datafilename <- system.file("extdata", "transcripto_sample.txt", package="DRomics")
o <- microarraydata(datafilename, check = TRUE, norm.method = "cyclicloess")
s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.001)
f <- drcfit(s_quad, progressbar = TRUE)
r <- bmdcalc(f)
set.seed(1234) # to get reproducible results with a so small number of iterations
(b <- bmdboot(r, niter = 5)) # with a non reasonable value for niter
# !!!! TO GET CORRECT RESULTS
# !!!! niter SHOULD BE FIXED FAR LARGER , e.g. to 1000
# !!!! but the run will be longer
b$res
plot(b) # plot of BMD.zSD after removing of BMDs with infinite upper bounds
# }
# NOT RUN {
plot(b, remove.infinite = FALSE) # plot of BMD.zSD without removing of BMDs
# with infinite upper bounds
# }
# NOT RUN {
# bootstrap on only a subsample of items, here those best fitted by the linear model
# with a greater number of iterations
# }
# NOT RUN {
(b.lin.95 <- bmdboot(r, items = r$res$id[r$res$model == "Gauss-probit"],
niter = 1000, progressbar = TRUE))
b.lin.95$res
# same bootstrap but changing the default confidence level (0.95) to 0.90
(b.lin.90 <- bmdboot(r, items = r$res$id[r$res$model == "Gauss-probit"],
niter = 1000, conf.level = 0.9, progressbar = TRUE))
b.lin.90$res
# }
# NOT RUN {
# (2) an example on a microarray data set (a subsample of a greater data set)
#
# }
# NOT RUN {
datafilename <- system.file("extdata", "transcripto_sample.txt", package="DRomics")
(o <- microarraydata(datafilename, check = TRUE, norm.method = "cyclicloess"))
(s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.001))
(f <- drcfit(s_quad, progressbar = TRUE))
(r <- bmdcalc(f))
(b <- bmdboot(r, niter = 100)) # niter to put at 1000 for a better precision
# different plots of BMD-zSD
plot(b)
if (require(ggplot2)) plot(b) + scale_x_log10() # in BMD log10 scale
plot(b, by = "trend")
plot(b, by = "model")
plot(b, by = "typology")
# a plot of BMD-xfold (by default BMD-zSD is plotted)
plot(b, BMDtype = "xfold")
# }
# NOT RUN {
# (3) Comparison of parallel and non parallel implementations
#
# }
# NOT RUN {
# to be tested with a greater number of iterations
system.time(b1 <- bmdboot(r, niter = 100, progressbar = TRUE))
system.time(b2 <- bmdboot(r, niter = 100, progressbar = FALSE, parallel = "snow", ncpus = 2))
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab