search

DoseResponse / drc

Fitting dose-response models in R
https://doseresponse.github.io/drc/
21 stars 16 forks source link

predict cannot calculate confidence intervals for models with fixed parameters #40

Open FalkoHof opened 1 month ago

FalkoHof commented 1 month ago

Hey there, I was trying to calculate RNA-seq read saturation and ran into the issue:

Is this related to drm returning a model with 1 less parameter? Or somewhat intended behaviour as a fixed parameter contains no uncertainty?

When using MM.2() kinetics with fixed asymptote:

# input data
> head(sat)
   saturation total_reads
        <num>       <int>
1:  0.1685059      494784
2:  0.2575534     2974401
3:  0.3329197     5445163
4:  0.3980089     7920298
5:  0.4538969    10405968
6:  0.5019354    12884013

> model.drm <- drm(saturation ~ total_reads, data = sat, fct = MM.2(fixed = c(1,NA), names=c("d","e")))
# data to extrapolate and predict confidence intervals to
> mml <- data.frame(reads = seq(0, max(sat$total_reads)*2, length.out = 1000))
> summary(model.drm)
Model fitted: Michaelis-Menten (1 parms)
Parameter estimates:

              Estimate Std. Error t-value   p-value    
e:(Intercept) 11937044     626398  19.057 2.214e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error:

 0.03482782 (18 degrees of freedom)

# prediction fails with the error below
> pred <- predict(model.drm, newdata = mml, interval = "confidence", level = 0.95)
Error in indexMat[, groupLevels, drop = FALSE] : 
  incorrect number of dimensions

When using MM.2() kinetics with no fixed parameters:

> model.drm <- drm(sat ~ total_reads, data = sat, fct = MM.2())
> summary(model.drm)

Model fitted: Michaelis-Menten (2 parms)

Parameter estimates:

                Estimate Std. Error t-value p-value    
d:(Intercept) 9.9953e-01 1.3612e-02  73.432  <2e-16 ***
e:(Intercept) 1.1920e+07 8.0763e+05  14.759   4e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error:

 0.03583626 (17 degrees of freedom)

> pred <- predict(model.drm, newdata = mml, interval = "confidence", level = 0.95)
> head(pred)

 Prediction      Lower      Upper
[1,] 0.00000000 0.00000000 0.00000000
[2,] 0.07677375 0.06779986 0.08574764
[3,] 0.14259479 0.12724126 0.15794832
[4,] 0.19965087 0.17974278 0.21955896
[5,] 0.24958346 0.22642738 0.27273955
[6,] 0.29364824 0.26819302 0.31910345

Session Info

> sessionInfo()
R version 4.4.0 (2024-04-24)
Platform: x86_64-pc-linux-gnu
Running under: Rocky Linux 9.4 (Blue Onyx)

Matrix products: default
BLAS/LAPACK: FlexiBLAS OPENBLAS;  LAPACK version 3.11.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Berlin
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] drc_3.0-1                   MASS_7.3-60.2               pbapply_1.7-2               scales_1.3.0               
 [5] rhdf5_2.48.0                glue_1.7.0                  ggtext_0.1.2                lubridate_1.9.3            
 [9] forcats_1.0.0               stringr_1.5.1               dplyr_1.1.4                 purrr_1.0.2                
[13] readr_2.1.5                 tidyr_1.3.1                 tibble_3.2.1                ggplot2_3.5.1              
[17] tidyverse_2.0.0             DropletUtils_1.24.0         SingleCellExperiment_1.26.0 SummarizedExperiment_1.34.0
[21] Biobase_2.64.0              GenomicRanges_1.56.0        GenomeInfoDb_1.40.0         IRanges_2.38.0             
[25] S4Vectors_0.42.0            BiocGenerics_0.50.0         MatrixGenerics_1.16.0       matrixStats_1.3.0          
[29] data.table_1.15.4          

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1          farver_2.1.2              R.utils_2.12.3            TH.data_1.1-2            
 [5] timechange_0.3.0          lifecycle_1.0.4           survival_3.6-4            statmod_1.5.0            
 [9] magrittr_2.0.3            compiler_4.4.0            rlang_1.1.4               tools_4.4.0              
[13] plotrix_3.8-4             utf8_1.2.4                labeling_0.4.3            S4Arrays_1.4.0           
[17] dqrng_0.4.1               DelayedArray_0.30.1       xml2_1.3.6                pkgload_1.3.4            
[21] multcomp_1.4-25           abind_1.4-5               BiocParallel_1.38.0       HDF5Array_1.32.0         
[25] withr_3.0.0               R.oo_1.26.0               grid_4.4.0                fansi_1.0.6              
[29] beachmat_2.20.0           colorspace_2.1-0          Rhdf5lib_1.26.0           edgeR_4.2.0              
[33] gtools_3.9.5              mvtnorm_1.2-4             cli_3.6.3                 crayon_1.5.3             
[37] generics_0.1.3            rstudioapi_0.16.0         httr_1.4.7                tzdb_0.4.0               
[41] DelayedMatrixStats_1.26.0 scuttle_1.14.0            splines_4.4.0             zlibbioc_1.50.0          
[45] parallel_4.4.0            XVector_0.44.0            vctrs_0.6.5               sandwich_3.1-0           
[49] Matrix_1.7-0              jsonlite_1.8.8            carData_3.0-5             car_3.1-2                
[53] hms_1.1.3                 locfit_1.5-9.9            limma_3.60.0              codetools_0.2-20         
[57] stringi_1.8.4             gtable_0.3.5              UCSC.utils_1.0.0          munsell_0.5.1            
[61] pillar_1.9.0              rhdf5filters_1.16.0       GenomeInfoDbData_1.2.12   R6_2.5.1                 
[65] sparseMatrixStats_1.16.0  lattice_0.22-6            R.methodsS3_1.8.2         gridtext_0.1.5           
[69] Rcpp_1.0.12               SparseArray_1.4.3         zoo_1.8-12                pkgconfig_2.0.3

Best wishes, Falko