search

stan-dev / rstanarm

rstanarm R package for Bayesian applied regression modeling
https://mc-stan.org/rstanarm
GNU General Public License v3.0
388 stars 133 forks source link

feature request: implementing `se` method for `stanmvreg` objects #404

Open IndrajeetPatil opened 4 years ago

IndrajeetPatil commented 4 years ago

Running the following code gives me error:

Error in se.stanmvreg(x) : Not currently implemented

Reprex:

library(rstanarm)
#> Loading required package: Rcpp
#> Registered S3 method overwritten by 'xts':
#>   method     from
#>   as.zoo.xts zoo
#> rstanarm (Version 2.19.2, packaged: 2019-10-01 20:20:33 UTC)
#> - Do not expect the default priors to remain the same in future rstanarm versions.
#> Thus, R scripts should specify priors explicitly, even if they are just the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores())
#> - bayesplot theme set to bayesplot::theme_default()
#>    * Does _not_ affect other ggplot2 plots
#>    * See ?bayesplot_theme_set for details on theme setting

# Univariate joint model, with association structure based on the
# current value of the linear predictor
f1 <- stan_jm(formulaLong = logBili ~ year + (1 | id),
              dataLong = pbcLong,
              formulaEvent = Surv(futimeYears, death) ~ sex + trt,
              dataEvent = pbcSurv,
              time_var = "year",
              # this next line is only to keep the example small in size!
              chains = 1, cores = 1, seed = 12345, iter = 1000)
#> Loading required namespace: data.table
#> Fitting a univariate joint model.
#> 
#> Please note the warmup may be much slower than later iterations!
#> 
#> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
#> Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
#> Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 5.253 seconds (Warm-up)
#> Chain 1:                4.426 seconds (Sampling)
#> Chain 1:                9.679 seconds (Total)
#> Chain 1:
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#tail-ess

se(f1)

Created on 2019-11-29 by the reprex package (v0.3.0)

Session info ``` r devtools::session_info() #> - Session info --------------------------------------------------------------- #> setting value #> version R version 3.6.1 (2019-07-05) #> os Windows 10 x64 #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate English_United States.1252 #> ctype English_United States.1252 #> tz Europe/Berlin #> date 2019-11-29 #> #> - Packages ------------------------------------------------------------------- #> package * version date lib source #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0) #> backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.1) #> base64enc 0.1-3 2015-07-28 [1] CRAN (R 3.5.0) #> bayesplot 1.7.0 2019-05-23 [1] CRAN (R 3.6.0) #> boot 1.3-23 2019-07-05 [1] CRAN (R 3.6.1) #> callr 3.3.2 2019-09-22 [1] CRAN (R 3.6.1) #> cli 1.1.0 2019-03-19 [1] CRAN (R 3.6.0) #> codetools 0.2-16 2018-12-24 [1] CRAN (R 3.5.2) #> colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.0) #> colourpicker 1.0 2017-09-27 [1] CRAN (R 3.5.1) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.1) #> crosstalk 1.0.0 2016-12-21 [1] CRAN (R 3.5.1) #> data.table 1.12.6 2019-10-18 [1] CRAN (R 3.6.1) #> desc 1.2.0 2019-11-11 [1] Github (r-lib/desc@61205f6) #> devtools 2.2.1 2019-09-24 [1] CRAN (R 3.6.1) #> digest 0.6.23 2019-11-23 [1] CRAN (R 3.6.1) #> dplyr 0.8.3.9000 2019-10-10 [1] Github (tidyverse/dplyr@dcfc1d1) #> DT 0.10 2019-11-12 [1] CRAN (R 3.6.1) #> dygraphs 1.1.1.6 2018-07-11 [1] CRAN (R 3.5.1) #> ellipsis 0.3.0 2019-09-20 [1] CRAN (R 3.6.1) #> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0) #> fastmap 1.0.1 2019-10-08 [1] CRAN (R 3.6.1) #> fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0) #> ggplot2 3.2.1 2019-08-10 [1] CRAN (R 3.6.1) #> ggridges 0.5.1 2018-09-27 [1] CRAN (R 3.5.1) #> glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.0) #> gridExtra 2.3 2017-09-09 [1] CRAN (R 3.5.1) #> gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.0) #> gtools 3.8.1 2018-06-26 [1] CRAN (R 3.5.0) #> highr 0.8 2019-03-20 [1] CRAN (R 3.6.0) #> htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.1) #> htmlwidgets 1.5.1 2019-10-08 [1] CRAN (R 3.6.1) #> httpuv 1.5.2 2019-09-11 [1] CRAN (R 3.6.1) #> igraph 1.2.4.1 2019-04-22 [1] CRAN (R 3.5.3) #> inline 0.3.15 2018-05-18 [1] CRAN (R 3.5.1) #> knitr 1.26 2019-11-12 [1] CRAN (R 3.6.1) #> later 1.0.0 2019-10-04 [1] CRAN (R 3.6.1) #> lattice 0.20-38 2018-11-04 [2] CRAN (R 3.6.1) #> lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.6.0) #> lifecycle 0.1.0 2019-08-01 [1] CRAN (R 3.6.1) #> lme4 1.1-21 2019-03-05 [1] CRAN (R 3.6.0) #> loo 2.1.0 2019-03-13 [1] CRAN (R 3.6.0) #> magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.1) #> markdown 1.1 2019-08-07 [1] CRAN (R 3.6.1) #> MASS 7.3-51.4 2019-03-31 [1] CRAN (R 3.6.0) #> Matrix 1.2-18 2019-11-27 [1] CRAN (R 3.6.1) #> matrixStats 0.55.0 2019-09-07 [1] CRAN (R 3.6.1) #> memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.0) #> mime 0.7 2019-06-11 [1] CRAN (R 3.6.0) #> miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 3.5.1) #> minqa 1.2.4 2014-10-09 [1] CRAN (R 3.5.1) #> munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.1) #> nlme 3.1-140 2019-05-12 [2] CRAN (R 3.6.1) #> nloptr 1.2.1 2018-10-03 [1] CRAN (R 3.5.1) #> pillar 1.4.2 2019-06-29 [1] CRAN (R 3.6.0) #> pkgbuild 1.0.6 2019-10-09 [1] CRAN (R 3.6.1) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.1) #> pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.6.0) #> plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.1) #> prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.1) #> processx 3.4.1 2019-07-18 [1] CRAN (R 3.6.1) #> promises 1.1.0 2019-10-04 [1] CRAN (R 3.6.1) #> ps 1.3.0 2018-12-21 [1] CRAN (R 3.6.0) #> purrr 0.3.3 2019-10-18 [1] CRAN (R 3.6.1) #> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.1) #> Rcpp * 1.0.3 2019-11-08 [1] CRAN (R 3.6.1) #> remotes 2.1.0 2019-06-24 [1] CRAN (R 3.6.0) #> reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.5.1) #> rlang 0.4.2 2019-11-23 [1] CRAN (R 3.6.1) #> rmarkdown 1.18 2019-11-27 [1] CRAN (R 3.6.1) #> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.1) #> rsconnect 0.8.15 2019-07-22 [1] CRAN (R 3.6.1) #> rstan 2.19.2 2019-07-09 [1] CRAN (R 3.6.1) #> rstanarm * 2.19.2 2019-10-03 [1] CRAN (R 3.6.1) #> rstantools 2.0.0 2019-09-15 [1] CRAN (R 3.6.1) #> scales 1.1.0 2019-11-18 [1] CRAN (R 3.6.1) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0) #> shiny 1.4.0 2019-10-10 [1] CRAN (R 3.6.1) #> shinyjs 1.0 2018-01-08 [1] CRAN (R 3.5.1) #> shinystan 2.5.0 2018-05-01 [1] CRAN (R 3.5.1) #> shinythemes 1.1.2 2018-11-06 [1] CRAN (R 3.6.0) #> StanHeaders 2.19.0 2019-09-07 [1] CRAN (R 3.6.1) #> stringi 1.4.3 2019-03-12 [1] CRAN (R 3.6.0) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.0) #> survival 3.1-7 2019-11-09 [1] CRAN (R 3.6.1) #> testthat 2.3.0 2019-11-05 [1] CRAN (R 3.6.1) #> threejs 0.3.1 2017-08-13 [1] CRAN (R 3.5.1) #> tibble 2.1.3 2019-06-06 [1] CRAN (R 3.6.0) #> tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.1) #> usethis 1.5.1.9000 2019-11-28 [1] Github (r-lib/usethis@c7314cf) #> withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.1) #> xfun 0.11 2019-11-12 [1] CRAN (R 3.6.1) #> xtable 1.8-4 2019-04-21 [1] CRAN (R 3.5.3) #> xts 0.11-2 2018-11-05 [1] CRAN (R 3.6.0) #> yaml 2.2.0 2018-07-25 [1] CRAN (R 3.5.1) #> zoo 1.8-6 2019-05-28 [1] CRAN (R 3.6.0) #> #> [1] C:/Users/inp099/Documents/R/win-library/3.6 #> [2] C:/Program Files/R/R-3.6.1/library ```
sambrilleman commented 4 years ago

Will running:

f1$stan_summary[, "sd"]

essentially get you what you want?

@jgabry or @bgoodri might be able to explain why running:

count_data <- data.frame(
 counts = c(18,17,15,20,10,20,25,13,12),
 outcome = gl(3,1,9),
 treatment = gl(3,3)
)
fit3 <- stan_glm(
  counts ~ outcome + treatment, 
  data = count_data, 
  family = poisson(link="log"),
  prior = normal(0, 2, autoscale = FALSE),
  # for speed of example only
  chains = 2, iter = 250 
)
se(fit3)
fit3$stan_summary[, "sd"]

doesn't return the same thing, because I can't quite remember how/why the se.stanreg method differs from the posterior SD.