Hi @paul.buerkner ,
I’m trying to build on the above thread to implement a “switch model” as described in the supplementary information (page 2) of this manuscript:
It’s very similar to a piecewise linear model, only it has no slopes, and it allows for the variance to change before and after the “switch” (break) point.
I think I got the model structure right (?) but the fitting is on the slow-ish end and I get some convergence issues. Moreover, the effective sample size is quite low although the mean effects are correctly estimated. I was wondering if you have any ideas on how to improve this implementation to speed up the code and eliminate the convergence issues? Reprex below.
Thanks!
library(tidyverse)
library(brms)
rstan::rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
set.seed(10)
b0 <- -1
b1 <- 3
b2 <- -4
s1 <- 0.5
s2 <- 1
omega <- 30
group_error <- rnorm(50, 0, 0.01)
names(group_error) <- as.character(1:50)
df <- data.frame(predictor = seq(0.2, 1e2, 0.2)) %>%
dplyr::mutate(group_id = rep(1:50, each = 10) %>%
sample %>%
as.factor,
mean_error = ifelse(predictor < omega,
rnorm(sum(predictor < omega), 0, s1),
rnorm(sum(predictor >= omega), 0, s2)),
response = b0 + ifelse(predictor - omega < 0, b1, b2) +
mean_error + group_error[group_id])
ggplot(data = df, mapping = aes(y = response, x = predictor)) +
geom_point()
bform_switch <- brms::bf(response ~ b0 + b1 * step(omega - predictor) +
b2 * step(predictor - omega),
# keep omega within the range of predictor
brms::nlf(omega ~ inv_logit(alpha) * 1e2),
# allow error to change with switch
brms::nlf(sigma ~ s1 * step(omega - predictor) +
s2 * step(predictor - omega)),
b0 ~ 1 + (1 | group_id),
s1 + s2 + b1 + b2 + alpha ~ 1,
nl = TRUE)
bprior <- prior(normal(0, 2), nlpar = "b0") +
prior(normal(0, 2), nlpar = "b1") +
prior(normal(0, 2), nlpar = "b2") +
prior(normal(0, 1), nlpar = "alpha") +
prior(normal(0, 1), nlpar = "s1") +
prior(normal(0, 1), nlpar = "s2")
fit_s <- brms::brm(bform_switch, data = df, prior = bprior)
These are the returned warning messages I got:
1: There were 4000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
2: Examine the pairs() plot to diagnose sampling problems
3: The largest R-hat is 1.72, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
4: 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
5: 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
The outcome seems to capture well the simulated estimates (note that values of parameters s1 and s2 are estimated on the log scale, so expect s1 \approx log(0.5) and s2 \approx log(1)).
> fit_s
Family: gaussian
Links: mu = identity; sigma = log
Formula: response ~ b0 + b1 * step(omega - predictor) + b2 * step(predictor - omega)
omega ~ inv_logit(alpha) * 100
sigma ~ s1 * step(omega - predictor) + s2 * step(predictor - omega)
b0 ~ 1 + (1 | group_id)
s1 ~ 1
s2 ~ 1
b1 ~ 1
b2 ~ 1
alpha ~ 1
Data: df (Number of observations: 500)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~group_id (Number of levels: 50)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(b0_Intercept) 0.04 0.04 0.00 0.14 1.61 7 12
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b0_Intercept -1.61 0.52 -2.40 -0.37 1.29 12 25
s1_Intercept -0.60 0.06 -0.71 -0.48 1.01 463 860
s2_Intercept -0.02 0.04 -0.09 0.06 1.01 368 905
b1_Intercept 3.62 0.52 2.36 4.41 1.30 11 25
b2_Intercept -3.43 0.52 -4.70 -2.65 1.30 11 24
alpha_Intercept -0.85 0.00 -0.86 -0.85 1.00 961 2520
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Warning message:
Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results! We recommend running more iterations and/or setting stronger priors.
Including the before-after changes in residual error:
plot(brms::conditional_effects(fit_s, method = "posterior_predict"),
points = TRUE)
For the sake of transparency, this is output of my sessionInfo():
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] brms_2.14.4 Rcpp_1.0.5 forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4 readr_1.4.0 tidyr_1.1.2 tibble_3.0.4 ggplot2_3.3.2 tidyverse_1.3.0
loaded via a namespace (and not attached):
[1] minqa_1.2.4 colorspace_2.0-0 ellipsis_0.3.1 ggridges_0.5.2 rsconnect_0.8.16 estimability_1.3 markdown_1.1 base64enc_0.1-3 fs_1.5.0
[10] rstudioapi_0.13 farver_2.0.3 rstan_2.21.3 DT_0.16 fansi_0.4.1 mvtnorm_1.1-1 lubridate_1.7.9.2 xml2_1.3.2 codetools_0.2-16
[19] bridgesampling_1.0-0 splines_4.0.2 shinythemes_1.1.2 bayesplot_1.7.2 projpred_2.0.2 jsonlite_1.7.1 nloptr_1.2.2.2 broom_0.7.2 dbplyr_2.0.0
[28] shiny_1.5.0 compiler_4.0.2 httr_1.4.2 emmeans_1.5.0 backports_1.2.0 assertthat_0.2.1 Matrix_1.2-18 fastmap_1.0.1 cli_2.1.0
[37] later_1.1.0.1 prettyunits_1.1.1 htmltools_0.5.0 tools_4.0.2 igraph_1.2.6 coda_0.19-4 gtable_0.3.0 glue_1.4.2 reshape2_1.4.4
[46] V8_3.4.0 cellranger_1.1.0 vctrs_0.3.5 nlme_3.1-148 crosstalk_1.1.0.1 ps_1.4.0 lme4_1.1-25 rvest_0.3.6 mime_0.9
[55] miniUI_0.1.1.1 lifecycle_0.2.0 gtools_3.8.2 statmod_1.4.35 MASS_7.3-51.6 zoo_1.8-8 scales_1.1.1 colourpicker_1.1.0 hms_0.5.3
[64] promises_1.1.1 Brobdingnag_1.2-6 parallel_4.0.2 inline_0.3.16 shinystan_2.5.0 curl_4.3 gamm4_0.2-6 gridExtra_2.3 StanHeaders_2.21.0-6
[73] loo_2.3.1 stringi_1.5.3 dygraphs_1.1.1.6 boot_1.3-25 pkgbuild_1.1.0 rlang_0.4.8 pkgconfig_2.0.3 matrixStats_0.57.0 lattice_0.20-41
[82] labeling_0.4.2 rstantools_2.1.1 htmlwidgets_1.5.2 processx_3.4.4 tidyselect_1.1.0 plyr_1.8.6 magrittr_2.0.1 R6_2.5.0 generics_0.1.0
[91] DBI_1.1.0 pillar_1.4.6 haven_2.3.1 withr_2.3.0 mgcv_1.8-31 xts_0.12.1 abind_1.4-5 modelr_0.1.8 crayon_1.3.4
[100] grid_4.0.2 readxl_1.3.1 callr_3.5.1 threejs_0.3.3 reprex_0.3.0 digest_0.6.27 xtable_1.8-4 httpuv_1.5.4 RcppParallel_5.0.2
[109] stats4_4.0.2 munsell_0.5.0 shinyjs_2.0.0