R Under development (unstable) (2024-01-23 r85822 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Testof the impact of the three information criteria
> library(DRomics)
Loading required package: limma
Loading required package: DESeq2
Loading required package: S4Vectors
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: 'BiocGenerics'

The following object is masked from 'package:limma':

    plotMA

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
    as.data.frame, basename, cbind, colnames, dirname, do.call,
    duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
    lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
    pmin.int, rank, rbind, rownames, sapply, setdiff, table, tapply,
    union, unique, unsplit, which.max, which.min


Attaching package: 'S4Vectors'

The following object is masked from 'package:utils':

    findMatches

The following objects are masked from 'package:base':

    I, expand.grid, unname

Loading required package: IRanges

Attaching package: 'IRanges'

The following object is masked from 'package:grDevices':

    windows

Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: 'MatrixGenerics'

The following objects are masked from 'package:matrixStats':

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: 'Biobase'

The following object is masked from 'package:MatrixGenerics':

    rowMedians

The following objects are masked from 'package:matrixStats':

    anyMissing, rowMedians


DRomics has been loaded.
!!!! IMPORTANT CHANGES IN DEFAULT PLOT ARGUMENTS !!!!
Now all the plot functions use by default a log10 scale
for dose and BMD values, except curvesplot(),
for which the use of a dose log scale requires the specification
by the user of a non null minimal value (xmin).
We also put the default value of the argument scaling at TRUE
in curvesplot() and bmdplotwithgradient(),
to focus on shapes of dose-responses rather than on their amplitude.

> visualize <- FALSE # put to TRUE for a manual check 
> 
> if (visualize)
+ {
+   ### test on microarray data ######################
+   datafilename <- system.file("extdata", "transcripto_very_small_sample.txt", package="DRomics")
+ 
+   (o <- microarraydata(datafilename, check = TRUE, norm.method = "cyclicloess"))
+   (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.05))
+   (fAIC <- drcfit(s_quad, information.criterion = "AIC", progressbar = TRUE))
+   (fAICc <- drcfit(s_quad, information.criterion = "AICc", progressbar = TRUE))
+   (fBIC <- drcfit(s_quad, information.criterion = "BIC", progressbar = TRUE))
+   
+   table(fAIC$fitres$model)
+   table(fAICc$fitres$model)
+   table(fBIC$fitres$model)
+   
+   head(fAIC$information.criterion.val)
+   head(fAICc$information.criterion.val)
+   head(fBIC$information.criterion.val)
+ 
+   # check of values on the linear model
+   (npts <- length(o$dose))
+   k <- 3 # mod lin
+   # correction to get AICc
+   fAIC$information.criterion.val$AIC.L + 2*k*(k+1)/(npts -k -1)
+   fAICc$information.criterion.val$AIC.L
+   # correction to get BIC
+   fAIC$information.criterion.val$AIC.L - 2*k + log(npts)*k
+   fBIC$information.criterion.val$AIC.L
+   
+   plot(fAIC)
+   plot(fAICc)
+   plot(fBIC)
+   
+   ################# on larger data sets
+   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.05))
+   (fAIC <- drcfit(s_quad, information.criterion = "AIC", progressbar = TRUE))
+   (fAICc <- drcfit(s_quad, information.criterion = "AICc", progressbar = TRUE))
+   (fBIC <- drcfit(s_quad, information.criterion = "BIC", progressbar = TRUE))
+   
+   table(fAIC$fitres$model)
+   table(fAICc$fitres$model)
+   table(fBIC$fitres$model)
+   
+ 
+   ### test on metabolo data with 4 doses ######################
+   data(Scenedesmus_metab)
+   head(Scenedesmus_metab)
+   # build of a dataset with 
+   Scenedesmus_metab2 <- Scenedesmus_metab[, c(1:8,10, 12, 14, 15, 18, 19, 22, 23)]
+   head(Scenedesmus_metab2)
+   
+   (o <- continuousomicdata(Scenedesmus_metab2))
+   (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.05))
+   (fAIC <- drcfit(s_quad, information.criterion = "AIC", progressbar = TRUE))
+   (fAICc <- drcfit(s_quad, information.criterion = "AICc", progressbar = TRUE))
+   (fBIC <- drcfit(s_quad, information.criterion = "BIC", progressbar = TRUE))
+   
+   table(fAIC$fitres$model)
+   table(fAICc$fitres$model)
+   table(fBIC$fitres$model)
+   
+   table(fAIC$fitres$nbpar)
+   table(fAICc$fitres$nbpar)
+   table(fBIC$fitres$nbpar)
+   
+   plot(fAIC)
+   plot(fAICc)
+   plot(fBIC)
+   
+   ### test on RNAseq data with 5 doses ######################
+   data(Zhou_kidney_pce)
+   head(Zhou_kidney_pce)
+ 
+   (o <- RNAseqdata(Zhou_kidney_pce))
+   (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.01))
+   (fAIC <- drcfit(s_quad, information.criterion = "AIC", progressbar = TRUE))
+   (fAICc <- drcfit(s_quad, information.criterion = "AICc", progressbar = TRUE))
+   (fBIC <- drcfit(s_quad, information.criterion = "BIC", progressbar = TRUE))
+   
+   table(fAIC$fitres$model)
+   table(fAICc$fitres$model)
+   table(fBIC$fitres$model)
+   
+   table(fAIC$fitres$nbpar)
+   table(fAICc$fitres$nbpar)
+   table(fBIC$fitres$nbpar)
+   
+   plot(fAIC, 81)
+   plot(fAICc, 81)
+   plot(fBIC, 81)
+   
+   # exploration of simplified biphasic models with f = 0
+   # f <- fAIC
+   # f <- fBIC
+   f <- AICc
+   (id2explore <- f$fitres$id[f$fitres$model %in% c("Gauss-probit", "log-Gauss-probit") & 
+                                f$fitres$f == 0])
+   f$fitres[f$fitres$id %in%  id2explore, ]
+   plot(f, items = id2explore)
+   
+   ###### test on apical data
+   data(Scenedesmus_apical)
+   head(Scenedesmus_apical)
+   (o <- continuousanchoringdata(Scenedesmus_apical, backgrounddose = 0.1))
+   (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.01))
+   (fAIC <- drcfit(s_quad, information.criterion = "AIC", progressbar = TRUE))
+   (fAICc <- drcfit(s_quad, information.criterion = "AICc", progressbar = TRUE))
+   (fBIC <- drcfit(s_quad, information.criterion = "BIC", progressbar = TRUE))
+   
+   plot(fAIC)
+   plot(fAICc)
+   plot(fBIC)
+ 
+   ###### test on in situ RNAseq data
+   datafilename <- system.file("extdata", "insitu_RNAseq_sample.txt", package="DRomics")
+   (o <- RNAseqdata(datafilename, backgrounddose = 2e-2, transfo.method = "rlog"))
+   (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.01))
+   (fAIC <- drcfit(s_quad, information.criterion = "AIC", progressbar = TRUE))
+   (fAICc <- drcfit(s_quad, information.criterion = "AICc", progressbar = TRUE))
+   (fBIC <- drcfit(s_quad, information.criterion = "BIC", progressbar = TRUE))
+ 
+   table(fAIC$fitres$model)
+   table(fAICc$fitres$model)
+   table(fBIC$fitres$model)
+   
+   table(fAIC$fitres$nbpar)
+   table(fAICc$fitres$nbpar)
+   table(fBIC$fitres$nbpar)
+ 
+   plot(fAIC, dose_log_transfo = TRUE)
+   plot(fAICc, dose_log_transfo = TRUE)
+   plot(fBIC, dose_log_transfo = TRUE)
+   
+   nrow(fAIC$fitres)
+   nrow(fAICc$fitres)
+   nrow(fBIC$fitres)
+   
+   id2compare <- fBIC$fitres$id[50:70]
+   plot(fAIC, items = id2compare, dose_log_transfo = TRUE)
+   plot(fAICc, items = id2compare, dose_log_transfo = TRUE)
+   plot(fBIC, items = id2compare, dose_log_transfo = TRUE)
+   
+   # exploration of simplified biphasic models with f = 0
+   # f <- fAIC
+   # f <- fBIC
+   f <- fAICc
+   (id2explore <- f$fitres$id[f$fitres$model %in% c("Gauss-probit", "log-Gauss-probit") & 
+                                f$fitres$f == 0])
+   f$fitres[f$fitres$id %in%  id2explore, ]
+   plot(f, items = id2explore, dose_log_transfo = TRUE)
+   
+   head(fAIC$information.criterion.val, 20)
+   head(fAICc$information.criterion.val, 20)
+   head(fBIC$information.criterion.val, 20)
+   
+   
+ }
> 
> proc.time()
   user  system elapsed 
   8.23    0.82    9.06