saudiwin / ordbetareg_pack

Repository for R package ordbetareg, used to fit continuous data with lower and upper bounds.
Other
17 stars 3 forks source link

manual intercepts via `intercept_prior_mean` give error #14

Closed stefgehrig closed 1 year ago

stefgehrig commented 1 year ago

Since v0.7.0, handing intercept_prior_mean = 0 and intercept_prior_SD = 5 should allow to configure a prior N(0,5) for the intercept. However, the error Error: The following priors do not correspond to any model parameter: b_Intercept ~ normal(0,5) is returned.

Reproducible example:

library(dplyr)
data("pew")

model_data <- select(pew,therm, education="F_EDUCCAT2_FINAL", region="F_CREGION_FINAL", income="F_INCOME_FINAL")

ord_fit_mean <- ordbetareg(formula=therm ~ education + income + (1|region),
    data=model_data,
    intercept_prior_mean = 0,
    intercept_prior_SD  = 5,
    cores=2,chains=2)

Session info

R version 4.2.3 (2023-03-15 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 22621)

Matrix products: default

locale: [1] LC_COLLATE=English_Germany.utf8 LC_CTYPE=English_Germany.utf8 LC_MONETARY=English_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=English_Germany.utf8

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

other attached packages: [1] dplyr_1.1.1 ordbetareg_0.7.1 brms_2.19.0 Rcpp_1.0.10

loaded via a namespace (and not attached): [1] TH.data_1.1-1 colorspace_2.1-0 ellipsis_0.3.2 class_7.3-21 estimability_1.4.1 markdown_1.6 base64enc_0.1-3
[8] rstudioapi_0.14 proxy_0.4-27 farver_2.1.1 rstan_2.26.22 DT_0.27 fansi_1.0.4 mvtnorm_1.1-3
[15] bridgesampling_1.1-2 codetools_0.2-19 splines_4.2.3 shinythemes_1.2.0 bayesplot_1.10.0 jsonlite_1.8.4 gganimate_1.0.8
[22] shiny_1.7.4 compiler_4.2.3 emmeans_1.8.5 backports_1.4.1 Matrix_1.5-4 fastmap_1.1.1 cli_3.6.1
[29] later_1.3.0 tweenr_2.0.2 htmltools_0.5.5 prettyunits_1.1.1 tools_4.2.3 igraph_1.4.2 coda_0.19-4
[36] gtable_0.3.3 glue_1.6.2 reshape2_1.4.4 posterior_1.4.1 V8_4.3.0 vctrs_0.6.1 nlme_3.1-162
[43] transformr_0.1.4 crosstalk_1.2.0 tensorA_0.36.2 stringr_1.5.0 ps_1.7.4 mime_0.12 lpSolve_5.6.18
[50] miniUI_0.1.1.1 lifecycle_1.0.3 gtools_3.9.4 MASS_7.3-58.3 zoo_1.8-12 scales_1.2.1 colourpicker_1.2.0
[57] hms_1.1.3 promises_1.2.0.1 Brobdingnag_1.2-9 parallel_4.2.3 sandwich_3.0-2 inline_0.3.19 shinystan_2.6.0
[64] curl_5.0.0 gridExtra_2.3 ggplot2_3.4.2 loo_2.6.0 StanHeaders_2.26.22 stringi_1.7.12 dygraphs_1.1.1.6
[71] e1071_1.7-13 checkmate_2.1.0 pkgbuild_1.4.0 rlang_1.1.0 pkgconfig_2.0.3 matrixStats_0.63.0 distributional_0.3.2 [78] lattice_0.20-45 purrr_1.0.1 sf_1.0-12 rstantools_2.3.1 htmlwidgets_1.6.2 processx_3.8.0 tidyselect_1.2.0
[85] plyr_1.8.8 magrittr_2.0.3 R6_2.5.1 magick_2.7.4 generics_0.1.3 multcomp_1.4-23 DBI_1.1.3
[92] withr_2.2.0 pillar_1.9.0 units_0.8-1 xts_0.13.0 survival_3.5-5 abind_1.4-5 tibble_3.2.1
[99] crayon_1.5.2 KernSmooth_2.23-20 utf8_1.2.3 progress_1.2.2 grid_4.2.3 callr_3.7.3 threejs_0.3.3
[106] digest_0.6.31 classInt_0.4-9 xtable_1.8-4 tidyr_1.3.0 httpuv_1.6.9 RcppParallel_5.1.7 stats4_4.2.3
[113] munsell_0.5.0 shinyjs_2.1.0

stefgehrig commented 1 year ago

Since 0 + Intercept was removed, the intercept is in the current package version not treated as an population-level effect anymore by default. So rather than the currently implemented

  if(!is.null(intercept_prior)) {

    priors <- priors + set_prior(paste0("normal(",intercept_prior[1],",",intercept_prior[2],")"),
                                 coef="Intercept",class="b")

we'd need to specify it via class = "Intercept".

This is assuming that what is specified via intercept_prior_mean and intercept_prior_SD is meant to be the prior on the internally centered intercept (brms: "when internally centering all population-level predictors around zero"), consistent with the default brms procedure. Probably it would be good to be explicit here in the documentation whether this option to ordbetareg() is intended for the centered or uncentered intercept.

saudiwin commented 1 year ago

Thanks much for doing the detective work on this, I should be able to get it fixed shortly.