Error in rcpp_mSpline with Cox model

[mbg-unsw]

I get an error when setting up a Cox model, apparently in setting up the spline baseline hazard, e.g.

> library(brms)
> dt<-data.frame(y=c(rep(10,100),3,5,7), c=c(rep(1,100),0,0,0))
> brm(y | cens(c) ~ 1, data = dt, family = brmsfamily("cox"))
Error in rcpp_mSpline(x = xx, df = df, degree = degree, internal_knots = knots,  : 
  Internal knots must be set inside of boundary knots.
> sessionInfo()
R version 4.0.4 (2021-02-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux bullseye/sid

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

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

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

other attached packages:
[1] brms_2.14.4 Rcpp_1.0.5 

loaded via a namespace (and not attached):
  [1] nlme_3.1-152         matrixStats_0.58.0   xts_0.12.1          
  [4] threejs_0.3.3        splines2_0.4.1       rstan_2.21.2        
  [7] backports_1.2.1      tools_4.0.4          R6_2.5.0            
 [10] DT_0.17              DBI_1.1.1            mgcv_1.8-34         
 [13] projpred_2.0.2       colorspace_2.0-0     withr_2.4.1         
 [16] tidyselect_1.1.0     gridExtra_2.3        prettyunits_1.1.1   
 [19] processx_3.4.5       Brobdingnag_1.2-6    emmeans_1.5.4       
 [22] curl_4.3             compiler_4.0.4       cli_2.3.0           
 [25] shinyjs_2.0.0        colourpicker_1.1.0   scales_1.1.1        
 [28] dygraphs_1.1.1.6     mvtnorm_1.1-1        ggridges_0.5.3      
 [31] callr_3.5.1          stringr_1.4.0        digest_0.6.27       
 [34] StanHeaders_2.21.0-7 minqa_1.2.4          base64enc_0.1-3     
 [37] pkgconfig_2.0.3      htmltools_0.5.1.1    lme4_1.1-26         
 [40] fastmap_1.1.0        htmlwidgets_1.5.3    rlang_0.4.10        
 [43] shiny_1.5.0          generics_0.1.0       zoo_1.8-8           
 [46] jsonlite_1.7.2       crosstalk_1.1.1      gtools_3.8.2        
 [49] dplyr_1.0.4          inline_0.3.17        magrittr_2.0.1      
 [52] loo_2.4.1            bayesplot_1.8.0      Matrix_1.3-2        
 [55] munsell_0.5.0        abind_1.4-5          lifecycle_0.2.0     
 [58] stringi_1.5.3        MASS_7.3-53.1        pkgbuild_1.2.0      
 [61] plyr_1.8.6           grid_4.0.4           parallel_4.0.4      
 [64] promises_1.1.1       crayon_1.4.0         miniUI_0.1.1.1      
 [67] lattice_0.20-41      splines_4.0.4        ps_1.5.0            
 [70] pillar_1.4.7         igraph_1.2.6         boot_1.3-27         
 [73] estimability_1.3     markdown_1.1         shinystan_2.5.0     
 [76] reshape2_1.4.4       codetools_0.2-18     stats4_4.0.4        
 [79] rstantools_2.1.1     glue_1.4.2           V8_3.4.0            
 [82] RcppParallel_5.0.2   vctrs_0.3.6          nloptr_1.2.2.2      
 [85] httpuv_1.5.5         gtable_0.3.0         purrr_0.3.4         
 [88] assertthat_0.2.1     ggplot2_3.3.3        mime_0.9            
 [91] xtable_1.8-4         coda_0.19-4          later_1.1.0.1       
 [94] rsconnect_0.8.16     tibble_3.0.6         shinythemes_1.2.0   
 [97] gamm4_0.2-6          statmod_1.4.35       ellipsis_0.3.1      
[100] bridgesampling_1.0-0
> 

Thanks for your help and thank you for brms!

[paul-buerkner]

Thanks. Will investigate.

[mbg-unsw]

Hi Paul, can you give me a rough idea when you're likely to look at this?

I am planning a project and if a fix is a few months off, then I'll consider another approach.

Thank you again.

[andrjohns]

@paul-buerkner I can start digging into it this afternoon and let you know what I find

[paul-buerkner]

thank you! Andrew Johnson ***@***.***> schrieb am Fr., 7. Mai 2021, 02:38:

…

@paul-buerkner <https://github.com/paul-buerkner> I can start digging into it this afternoon and let you know what I find — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1143 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCW2AG3K6B2YRY2CFE6GRLTMMZAHANCNFSM43EYGA4A> .

[paul-buerkner]

So the problem is that the brms-set default boundary knots don't match the automatically chosen internal knots by splines2::mSpline. This is strange since the default boundary knots are simply c(min(y), max(y)) which is supposably the default of splines2::mSpline anyway according to the documentation. And indeed, when I turn off the default computation of boundary knots, the error still occurs. So there might be something off with the default settings of (boundary and internal) knots in splines2::mSpline that I don't know hot to fix from the brms side.

Fortunately, you have full control over the knots and boundary knots via the bhaz argument. For example:

library(brms)
dt<-data.frame(y=c(rep(10,100),3,5,7), c=c(rep(1,100),0,0,0))
bhaz <- list(Boundary.knots = c(0, 150))
brm(y | cens(c) ~ 1, data = dt, family = brmsfamily("cox", bhaz = bhaz))

where bhaz is a list of arguments direclty passed to splines2::mSpline

[paul-buerkner]

I think the error occurs if the variation in the response y is so small that certain tail-quantiles are exactly equal to min(y) or max(y). I will push a fix to brms master that slightly adjusts the boundary knots but be a little more more extreme than min(y) and max(y), which should fix the problem.

@andrjohns thank you for offering you help on this. I appreciate it a lot. I had scheduled checking out some brms bugs anyway and was reminded by this issue, so I jumped right into it before other TODOs take over again.

[mbg-unsw]

Thanks Paul! That explanation makes sense to me.

I will watch out for the commits and test when ready.

[paul-buerkner]

should be fixed now. :-)