Correction of the BUG in the definition of the starting values for parameter e in the fit of the Gaussprobit model.

Author: mldelignette

NEWS.md

@@ -2,6 +2,14 @@
 
 NEW FEATURES
 
+BUG FIXES
+
+- Fix a bug in the definition of starting values of the parameter e of the Gaussprobit model. This could help the fit in few cases.
+
+# DRomics 2.6-2
+
+NEW FEATURES
+
 - Add of a component in the output of RNAseqdata, continuousomicdata(),  continuousanchoringdata(), microarraydata() and (data.sd, which gives, for each item, the sd of the response per condition - NA if no replicate for a condition).
 - Add a new argument in drcfit(), named deltaAICminfromnullmodel, in order to relax requirements on the information criterion to keep the best fitted model (see ? drcfit()).
 - Modification of curvesplot() to be able to put its argument dose_log_transfo by default at TRUE as in other functions. Curvesplot now use minimum and maximum values of the chosen BMD to fix the rage on which the theoretical curve is calculated (and plotted) ad so the chosen BMD is required in the input of the function.

R/util-basicandfitfunc.R

@@ -163,7 +163,7 @@ startvalGauss5pnls <- function(xm, ym, Ushape) {
   b <- max(xm) / 4
   # initial value of e (dose corresponding to the maximal (or minimal) signal)
   xextremum <- stats::median(xm[which(ym == yextremum)]) # just in case there is more than one dose at which ym == yextremum
-  e <- min(xextremum - (c - d) * b / (f * sqrt(2 * pi)), 1e-6)
+  e <- max(xextremum - (c - d) * b / (f * sqrt(2 * pi)), 1e-6)
   startval <- list(b = b, c = c, d = d, e = e, f = f)
 }
 

share/techdoc/impact_of_GP_startingvalues_BUG.R

@@ -0,0 +1,37 @@
+Results before the correction of the BUG in starting values of e for the GP model.extract(_______________________________________________________________________________________)
+> datafilename <- system.file("extdata", "transcripto_sample.txt", package = "DRomics")
+> 
+  > (o <- microarraydata(datafilename, check = TRUE, norm.method = "cyclicloess"))
+Just wait, the normalization using cyclicloess may take a few minutes.
+Elements of the experimental design in order to check the coding of the data:
+  Tested doses and number of replicates for each dose:
+  
+  0  0.69 1.223 2.148 3.774 6.631 
+5     5     5     5     5     5 
+Number of items: 1000 
+Identifiers of the first 20 items:
+  [1] "1"    "2"    "3"    "4"    "5.1"  "5.2"  "6.1"  "6.2"  "7.1"  "7.2"  "8.1"  "8.2"  "9.1" 
+[14] "9.2"  "10.1" "10.2" "11.1" "11.2" "12.1" "12.2"
+Data were normalized between arrays using the following method: cyclicloess 
+> (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.05))
+Removing intercept from test coefficients
+Number of selected items using a quadratic trend test with an FDR of 0.05: 318
+Identifiers of the first 20 most responsive items:
+  [1] "384.2" "383.1" "383.2" "384.1" "301.1" "363.1" "300.2" "364.2" "364.1" "363.2" "301.2" "300.1"
+[13] "351.1" "350.2" "239.1" "240.1" "240.2" "370"   "15"    "350.1"
+> (f <- drcfit(s_quad, progressbar = TRUE))
+The fitting may be long if the number of selected items is high.
+|===========================================================================================| 100%
+Results of the fitting using the AICc to select the best fit model
+60 dose-response curves out of 318 previously selected were removed because no model 
+could be fitted reliably.
+Distribution of the chosen models among the 258 fitted dose-response curves:
+  
+  Hill           linear      exponential     Gauss-probit log-Gauss-probit 
+2               85               64               89               18 
+Distribution of the trends (curve shapes) among the 258 fitted dose-response curves:
+  
+  bell  dec  inc    U 
+52   59   92   55 
+
+Results after the correction of the BUG in starting values of e for the GP model.extract(_______________________________________________________________________________________)

tests/testthat/test_windows_drcfit.R

@@ -17,7 +17,7 @@ test_that("drcfit works as expected on the model results",
     expect_equal(round(mean(f$fitres$c, na.rm = TRUE), 4), 5.3354)
     expect_equal(round(mean(f$fitres$d, na.rm = TRUE), 4), 9.2984)
     expect_equal(round(mean(f$fitres$e, na.rm = TRUE), 4), 1.608)
-    expect_equal(round(mean(f$fitres$f, na.rm = TRUE), 4), 2.5344)
+    expect_equal(round(mean(f$fitres$f, na.rm = TRUE), 4), 2.5343)
     #Select model with BIC
     f.BIC <- drcfit(s_quad, progressbar = TRUE, information.criterion = "BIC")
     tmodel.BIC <- table(f.BIC$fitres$model)
@@ -35,9 +35,9 @@ test_that("drcfit works as expected on the model results",
     expect_equal(as.numeric(tmodel.AICc["Gauss-probit"]), 23)
     expect_equal(as.numeric(tmodel.AICc["log-Gauss-probit"]), 3)
     expect_equal(round(mean(f.AICc$fitres$b, na.rm = TRUE), 4), 1.2561)
-    expect_equal(round(mean(f.AICc$fitres$c, na.rm = TRUE), 4), 4.5015)
+    expect_equal(round(mean(f.AICc$fitres$c, na.rm = TRUE), 4), 4.5016)
     expect_equal(round(mean(f.AICc$fitres$d, na.rm = TRUE), 4), 9.2571)
     expect_equal(round(mean(f.AICc$fitres$e, na.rm = TRUE), 4), 1.5162)
-    expect_equal(round(mean(f.AICc$fitres$f, na.rm = TRUE), 4), 2.8576)
+    expect_equal(round(mean(f.AICc$fitres$f, na.rm = TRUE), 4), 2.8575)
 
   })