drc
drc copied to clipboard
DoseResponse
predict cannot calculate confidence intervals for models with fixed parameters
Hey there, I was trying to calculate RNA-seq read saturation and ran into the issue:
predictseems to be unable to calculate confidence intervals when using fixed parameters.
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
