absoluteRisk error when log of time interacts with another variable
Opened this issue · 4 comments
sahirbhatnagar commented
I'm not really sure how to debug this. Any thoughts @turgeonmaxime ?
It really seems to just be when log of time is interacting with another variable, because a spline interaction works.
library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
library(splines)
data("ERSPC")
set.seed(12)
ERSPC$ScrArm <- factor(ERSPC$ScrArm,
levels = c(0,1),
labels = c("Control group", "Screening group"))
fit1 <- casebase::fitSmoothHazard(DeadOfPrCa ~ Follow.Up.Time * ScrArm, data = ERSPC)
#> 'Follow.Up.Time' will be used as the time variable
fit2 <- casebase::fitSmoothHazard(DeadOfPrCa ~ log(Follow.Up.Time) + ScrArm, data = ERSPC)
#> 'Follow.Up.Time' will be used as the time variable
fit3 <- casebase::fitSmoothHazard(DeadOfPrCa ~ log(Follow.Up.Time) * ScrArm, data = ERSPC)
#> 'Follow.Up.Time' will be used as the time variable
fit4 <- casebase::fitSmoothHazard(DeadOfPrCa ~ bs(Follow.Up.Time) * ScrArm, data = ERSPC)
#> 'Follow.Up.Time' will be used as the time variable
new_data <- data.frame(ScrArm = c("Control group", "Screening group"))
new_time <- seq(1,12,0.1)
head(casebase::absoluteRisk(object = fit1, time = new_time, newdata = new_data))
#> time
#> 0.0 0.0000000000 0.0000000000
#> 1.0 0.0001190331 0.0001278290
#> 1.1 0.0001325949 0.0001419974
#> 1.2 0.0001464883 0.0001564372
#> 1.3 0.0001607213 0.0001711537
#> 1.4 0.0001753022 0.0001861522
head(casebase::absoluteRisk(object = fit2, time = new_time, newdata = new_data))
#> time
#> 0.0 0.000000e+00 0.000000e+00
#> 1.0 3.818107e-05 3.069906e-05
#> 1.1 4.651259e-05 3.739795e-05
#> 1.2 5.569660e-05 4.478230e-05
#> 1.3 6.573835e-05 5.285632e-05
#> 1.4 7.664273e-05 6.162395e-05
head(casebase::absoluteRisk(object = fit3, time = new_time, newdata = new_data))
#> time
#> 0.0 0 0
#> 1.0 NaN NaN
#> 1.1 NaN NaN
#> 1.2 NaN NaN
#> 1.3 NaN NaN
#> 1.4 NaN NaN
head(casebase::absoluteRisk(object = fit4, time = new_time, newdata = new_data))
#> time
#> 0.0 0.000000e+00 0.000000e+00
#> 1.0 4.999085e-05 4.507590e-05
#> 1.1 5.738346e-05 5.184092e-05
#> 1.2 6.530634e-05 5.911011e-05
#> 1.3 7.378433e-05 6.690794e-05
#> 1.4 8.284244e-05 7.525902e-05Created on 2020-08-20 by the reprex package (v0.3.0)
Session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 3.6.2 (2019-12-12)
#> os Pop!_OS 19.10
#> system x86_64, linux-gnu
#> ui X11
#> language en_US:en
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz America/Toronto
#> date 2020-08-20
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date lib source
#> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.2)
#> backports 1.1.8 2020-06-17 [1] CRAN (R 3.6.2)
#> callr 3.4.3 2020-03-28 [1] CRAN (R 3.6.2)
#> casebase * 0.9.0 2020-07-03 [1] CRAN (R 3.6.2)
#> cli 2.0.2 2020-02-28 [1] CRAN (R 3.6.2)
#> colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.2)
#> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.2)
#> data.table 1.12.8 2019-12-09 [1] CRAN (R 3.6.2)
#> desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.2)
#> devtools 2.2.2 2020-02-17 [1] CRAN (R 3.6.2)
#> digest 0.6.25 2020-02-23 [1] CRAN (R 3.6.2)
#> dplyr 1.0.0 2020-05-29 [1] CRAN (R 3.6.2)
#> ellipsis 0.3.1 2020-05-15 [1] CRAN (R 3.6.2)
#> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.2)
#> fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.2)
#> fs 1.3.2 2020-03-05 [1] CRAN (R 3.6.2)
#> generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.2)
#> ggplot2 3.3.2.9000 2020-07-02 [1] Github (tidyverse/ggplot2@3aa2937)
#> glue 1.4.1 2020-05-13 [1] CRAN (R 3.6.2)
#> gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.2)
#> highr 0.8 2019-03-20 [1] CRAN (R 3.6.2)
#> htmltools 0.5.0 2020-06-16 [1] CRAN (R 3.6.2)
#> knitr 1.29 2020-06-23 [1] CRAN (R 3.6.2)
#> lattice 0.20-38 2018-11-04 [4] CRAN (R 3.6.0)
#> lifecycle 0.2.0 2020-03-06 [1] CRAN (R 3.6.2)
#> magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.2)
#> Matrix 1.2-18 2019-11-27 [4] CRAN (R 3.6.1)
#> memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.2)
#> mgcv 1.8-31 2019-11-09 [4] CRAN (R 3.6.1)
#> munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.2)
#> nlme 3.1-143 2019-12-10 [4] CRAN (R 3.6.2)
#> pillar 1.4.4 2020-05-05 [1] CRAN (R 3.6.2)
#> pkgbuild 1.0.8 2020-05-07 [1] CRAN (R 3.6.2)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.2)
#> pkgload 1.1.0 2020-05-29 [1] CRAN (R 3.6.2)
#> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 3.6.2)
#> processx 3.4.2 2020-02-09 [1] CRAN (R 3.6.2)
#> ps 1.3.3 2020-05-08 [1] CRAN (R 3.6.2)
#> purrr 0.3.4 2020-04-17 [1] CRAN (R 3.6.2)
#> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.2)
#> Rcpp 1.0.4.6 2020-04-09 [1] CRAN (R 3.6.2)
#> remotes 2.1.1 2020-02-15 [1] CRAN (R 3.6.2)
#> rlang 0.4.6 2020-05-02 [1] CRAN (R 3.6.2)
#> rmarkdown 2.3 2020-06-18 [1] CRAN (R 3.6.2)
#> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.2)
#> scales 1.1.1 2020-05-11 [1] CRAN (R 3.6.2)
#> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.2)
#> stringi 1.4.6 2020-02-17 [1] CRAN (R 3.6.2)
#> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.2)
#> survival 3.1-8 2019-12-03 [4] CRAN (R 3.6.2)
#> testthat 2.3.2 2020-03-02 [1] CRAN (R 3.6.2)
#> tibble 3.0.1 2020-04-20 [1] CRAN (R 3.6.2)
#> tidyselect 1.1.0 2020-05-11 [1] CRAN (R 3.6.2)
#> usethis 1.6.1 2020-04-29 [1] CRAN (R 3.6.2)
#> vctrs 0.3.1 2020-06-05 [1] CRAN (R 3.6.2)
#> VGAM 1.1-3 2020-04-28 [1] CRAN (R 3.6.2)
#> withr 2.2.0 2020-04-20 [1] CRAN (R 3.6.2)
#> xfun 0.15 2020-06-21 [1] CRAN (R 3.6.2)
#> yaml 2.2.1 2020-02-01 [1] CRAN (R 3.6.2)
#>
#> [1] /home/sahir/R/x86_64-pc-linux-gnu-library/3.6
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/librarysahirbhatnagar commented
I've narrowed down the error to method="numerical" since, Monte Carlo integration works:
library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
data("ERSPC")
set.seed(12)
ERSPC$ScrArm <- factor(ERSPC$ScrArm,
levels = c(0,1),
labels = c("Control group", "Screening group"))
fit3 <- casebase::fitSmoothHazard(DeadOfPrCa ~ log(Follow.Up.Time) * ScrArm, data = ERSPC)
#> 'Follow.Up.Time' will be used as the time variable
new_data <- data.frame(ScrArm = c("Control group", "Screening group"))
new_time <- seq(1,12,0.1)
head(casebase::absoluteRisk(object = fit3, time = new_time, newdata = new_data, method = "mont"))
#> time
#> 0.0 0.000000e+00 0.000000e+00
#> 1.0 2.768649e-05 4.375383e-05
#> 1.1 3.431584e-05 5.260060e-05
#> 1.2 4.174256e-05 6.222870e-05
#> 1.3 4.991928e-05 7.257229e-05
#> 1.4 5.894599e-05 8.370706e-05Created on 2020-08-20 by the reprex package (v0.3.0)
Session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 3.6.2 (2019-12-12)
#> os Pop!_OS 19.10
#> system x86_64, linux-gnu
#> ui X11
#> language en_US:en
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz America/Toronto
#> date 2020-08-20
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date lib source
#> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.2)
#> backports 1.1.8 2020-06-17 [1] CRAN (R 3.6.2)
#> callr 3.4.3 2020-03-28 [1] CRAN (R 3.6.2)
#> casebase * 0.9.0 2020-07-03 [1] CRAN (R 3.6.2)
#> cli 2.0.2 2020-02-28 [1] CRAN (R 3.6.2)
#> colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.2)
#> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.2)
#> data.table 1.12.8 2019-12-09 [1] CRAN (R 3.6.2)
#> desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.2)
#> devtools 2.2.2 2020-02-17 [1] CRAN (R 3.6.2)
#> digest 0.6.25 2020-02-23 [1] CRAN (R 3.6.2)
#> dplyr 1.0.0 2020-05-29 [1] CRAN (R 3.6.2)
#> ellipsis 0.3.1 2020-05-15 [1] CRAN (R 3.6.2)
#> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.2)
#> fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.2)
#> fs 1.3.2 2020-03-05 [1] CRAN (R 3.6.2)
#> generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.2)
#> ggplot2 3.3.2.9000 2020-07-02 [1] Github (tidyverse/ggplot2@3aa2937)
#> glue 1.4.1 2020-05-13 [1] CRAN (R 3.6.2)
#> gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.2)
#> highr 0.8 2019-03-20 [1] CRAN (R 3.6.2)
#> htmltools 0.5.0 2020-06-16 [1] CRAN (R 3.6.2)
#> knitr 1.29 2020-06-23 [1] CRAN (R 3.6.2)
#> lattice 0.20-38 2018-11-04 [4] CRAN (R 3.6.0)
#> lifecycle 0.2.0 2020-03-06 [1] CRAN (R 3.6.2)
#> magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.2)
#> Matrix 1.2-18 2019-11-27 [4] CRAN (R 3.6.1)
#> memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.2)
#> mgcv 1.8-31 2019-11-09 [4] CRAN (R 3.6.1)
#> munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.2)
#> nlme 3.1-143 2019-12-10 [4] CRAN (R 3.6.2)
#> pillar 1.4.4 2020-05-05 [1] CRAN (R 3.6.2)
#> pkgbuild 1.0.8 2020-05-07 [1] CRAN (R 3.6.2)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.2)
#> pkgload 1.1.0 2020-05-29 [1] CRAN (R 3.6.2)
#> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 3.6.2)
#> processx 3.4.2 2020-02-09 [1] CRAN (R 3.6.2)
#> ps 1.3.3 2020-05-08 [1] CRAN (R 3.6.2)
#> purrr 0.3.4 2020-04-17 [1] CRAN (R 3.6.2)
#> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.2)
#> Rcpp 1.0.4.6 2020-04-09 [1] CRAN (R 3.6.2)
#> remotes 2.1.1 2020-02-15 [1] CRAN (R 3.6.2)
#> rlang 0.4.6 2020-05-02 [1] CRAN (R 3.6.2)
#> rmarkdown 2.3 2020-06-18 [1] CRAN (R 3.6.2)
#> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.2)
#> scales 1.1.1 2020-05-11 [1] CRAN (R 3.6.2)
#> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.2)
#> stringi 1.4.6 2020-02-17 [1] CRAN (R 3.6.2)
#> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.2)
#> survival 3.1-8 2019-12-03 [4] CRAN (R 3.6.2)
#> testthat 2.3.2 2020-03-02 [1] CRAN (R 3.6.2)
#> tibble 3.0.1 2020-04-20 [1] CRAN (R 3.6.2)
#> tidyselect 1.1.0 2020-05-11 [1] CRAN (R 3.6.2)
#> usethis 1.6.1 2020-04-29 [1] CRAN (R 3.6.2)
#> vctrs 0.3.1 2020-06-05 [1] CRAN (R 3.6.2)
#> VGAM 1.1-3 2020-04-28 [1] CRAN (R 3.6.2)
#> withr 2.2.0 2020-04-20 [1] CRAN (R 3.6.2)
#> xfun 0.15 2020-06-21 [1] CRAN (R 3.6.2)
#> yaml 2.2.1 2020-02-01 [1] CRAN (R 3.6.2)
#>
#> [1] /home/sahir/R/x86_64-pc-linux-gnu-library/3.6
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/libraryturgeonmaxime commented
My first guess is that it may have something to do with log(0) being equal to -Inf, but I thought I had taken care of that... I'll have a look
sahirbhatnagar commented
I had specified newtime > 0 for that reason. Also, I would have expected the error to occur in the additive model as well log(time) + ScrArm, but that was not the case
turgeonmaxime commented
I had specified newtime > 0 for that reason.
It doesn't matter, because absoluteRisk adds 0 back anyway:
Line 234 in 3403e04