I’m trying to fit an item response model (Bürkner, 2020) using brms. I’d like to fit a two-parameter IRT model that estimates difficulty/easiness (i.e., eta) and discrimination (i.e., alpha). The items are count data that have zero-inflation and heavy skewness/overdispersion, such that the standard deviation is greater than the mean. A negative binomial response distribution appears to fit the data better than a poisson distribution. I’ve been trying to model the response distribution with a hurdle negative binomial model to account for the zero-inflation, but am having difficulty getting the model to converge.

Because of difficulties getting the model to converge, I tried simulating data for which the model should be able to fit. However, I’m still having difficulty with model convergence (e.g., high R-hat values) even with the simulated data. A reproducible example is below. I simulated the items from a Poisson distribution (I’m not sure how to simulate correlated data from a negative binomial distribution).

Here is the code for the model:

```
brms_formula <- bf(
response ~ exp(logalpha) * eta,
eta ~ 1 + predictor + (1 | item) + (1 | subid),
logalpha ~ 1 + predictor + (1 | item),
shape ~ 1 + predictor,
hu ~ 1 + predictor,
nl = TRUE
)
model_fit <- brm(
formula = brms_formula,
data = simData,
family = brmsfamily("hurdle_negbinomial"),
chains = 1,
iter = 10000,
init = 0,
seed = 52242,
control = list(adapt_delta = 0.85),
backend = "cmdstanr"
)
```

Here is the reprodicible example:

```
library("simstudy")
library("brms")
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.20.4). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#>
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#>
#> ar
library("cmdstanr")
#> This is cmdstanr version 0.6.1
#> - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
#> - CmdStan path: D:/Documents/.cmdstan/cmdstan-2.33.1
#> - CmdStan version: 2.33.1
library("tidyverse")
# Simulate Data
set.seed(52242)
lambdas <- c(1, 8, 10, 12, 20)
simData <- genCorGen(
1000,
nvars = length(lambdas),
params1 = lambdas,
dist = "poisson",
rho = .6,
corstr = "cs",
wide = FALSE)
simData <- simData %>%
rename(
item = period,
subid = id,
response = X
)
simData$predictor <- simData$response + rnorm(nrow(simData), sd = 20)
cor(simData)
#> subid item response predictor
#> subid 1.000000000 0.0000000 0.01872772 0.002656721
#> item 0.000000000 1.0000000 0.85882741 0.279491061
#> response 0.018727723 0.8588274 1.00000000 0.324219762
#> predictor 0.002656721 0.2794911 0.32421976 1.000000000
# Model syntax
brms_formula <- bf(
response ~ exp(logalpha) * eta,
eta ~ 1 + predictor + (1 | item) + (1 | subid),
logalpha ~ 1 + predictor + (1 | item),
shape ~ 1 + predictor,
hu ~ 1 + predictor,
nl = TRUE
)
# Model priors
get_prior(
brms_formula,
data = simData,
family = brmsfamily("hurdle_negbinomial")
)
#> prior class coef group resp dpar nlpar lb ub
#> (flat) b eta
#> (flat) b Intercept eta
#> (flat) b predictor eta
#> student_t(3, 0, 2.5) sd eta 0
#> student_t(3, 0, 2.5) sd item eta 0
#> student_t(3, 0, 2.5) sd Intercept item eta 0
#> student_t(3, 0, 2.5) sd subid eta 0
#> student_t(3, 0, 2.5) sd Intercept subid eta 0
#> (flat) b logalpha
#> (flat) b Intercept logalpha
#> (flat) b predictor logalpha
#> student_t(3, 0, 2.5) sd logalpha 0
#> student_t(3, 0, 2.5) sd item logalpha 0
#> student_t(3, 0, 2.5) sd Intercept item logalpha 0
#> (flat) b hu
#> (flat) b predictor hu
#> logistic(0, 1) Intercept hu
#> (flat) b shape
#> (flat) b predictor shape
#> student_t(3, 0, 2.5) Intercept shape
#> source
#> default
#> (vectorized)
#> (vectorized)
#> default
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> default
#> (vectorized)
#> (vectorized)
#> default
#> (vectorized)
#> (vectorized)
#> default
#> (vectorized)
#> default
#> default
#> (vectorized)
#> default
# Fit the model
model_fit <- brm(
formula = brms_formula,
data = simData,
family = brmsfamily("hurdle_negbinomial"),
chains = 1,
iter = 10000,
init = 0,
seed = 52242,
control = list(adapt_delta = 0.85),
backend = "cmdstanr"
)
#> In file included from stan/lib/stan_math/stan/math/prim/prob/von_mises_lccdf.hpp:5,
#> from stan/lib/stan_math/stan/math/prim/prob/von_mises_ccdf_log.hpp:4,
#> from stan/lib/stan_math/stan/math/prim/prob.hpp:356,
#> from stan/lib/stan_math/stan/math/prim.hpp:16,
#> from stan/lib/stan_math/stan/math/rev.hpp:14,
#> from stan/lib/stan_math/stan/math.hpp:19,
#> from stan/src/stan/model/model_header.hpp:4,
#> from C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.hpp:2:
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp: In function 'stan::return_type_t<T_x, T_sigma, T_l> stan::math::von_mises_cdf(const T_x&, const T_mu&, const T_k&)':
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: '-Wmisleading-indentation' is disabled from this point onwards, since column-tracking was disabled due to the size of the code/headers
#> 194 | if (cdf_n < 0.0)
#> |
#> stan/lib/stan_math/stan/math/prim/prob/von_mises_cdf.hpp:194: note: adding '-flarge-source-files' will allow for more column-tracking support, at the expense of compilation time and memory
#> Start sampling
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is inf, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is inf, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is inf, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is inf, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is inf, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is inf, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is 0, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is inf, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is 0, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Log location parameter is inf, but must be finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: neg_binomial_2_log_lpmf: Precision parameter is 0, but must be positive finite! (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 71, column 6 to line 73, column 53) (in 'C:/Users/ITPETE~1/AppData/Local/Temp/RtmpK82p8f/model-cec4452265f.stan', line 195, column 6 to column 82)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Warning: 5000 of 5000 (100.0%) transitions ended with a divergence.
#> See https://mc-stan.org/misc/warnings for details.
# Model summary
summary(model_fit)
#> Warning: 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.
#> Warning: There were 5000 divergent transitions after warmup. Increasing
#> adapt_delta above 0.85 may help. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> Family: hurdle_negbinomial
#> Links: mu = log; shape = log; hu = logit
#> Formula: response ~ exp(logalpha) * eta
#> eta ~ 1 + predictor + (1 | item) + (1 | subid)
#> logalpha ~ 1 + predictor + (1 | item)
#> shape ~ 1 + predictor
#> hu ~ 1 + predictor
#> Data: simData (Number of observations: 5000)
#> Draws: 1 chains, each with iter = 10000; warmup = 5000; thin = 1;
#> total post-warmup draws = 5000
#>
#> Group-Level Effects:
#> ~item (Number of levels: 5)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> sd(eta_Intercept) 1.00 0.00 1.00 1.00 1.39 2
#> sd(logalpha_Intercept) 1.00 0.00 1.00 1.00 1.97 1
#> Tail_ESS
#> sd(eta_Intercept) 20
#> sd(logalpha_Intercept) 10
#>
#> ~subid (Number of levels: 1000)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(eta_Intercept) 1.00 0.00 1.00 1.00 1.72 2 30
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> shape_Intercept -2.52 0.01 -2.53 -2.50 1.76 2 12
#> hu_Intercept 0.05 0.00 0.05 0.06 1.52 2 11
#> eta_Intercept 0.08 0.00 0.08 0.08 1.89 1 17
#> eta_predictor 0.13 0.00 0.12 0.13 1.44 4 24
#> logalpha_Intercept 0.02 0.00 0.01 0.02 1.98 1 20
#> logalpha_predictor -0.02 0.00 -0.02 -0.02 1.70 2 17
#> shape_predictor 0.23 0.00 0.23 0.23 1.76 2 11
#> hu_predictor -0.01 0.00 -0.01 -0.01 1.44 2 11
#>
#> Draws were sampled using sample(hmc). 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).
sessionInfo()
#> R version 4.3.1 (2023-06-16 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19045)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/Chicago
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.3
#> [5] purrr_1.0.2 readr_2.1.4 tidyr_1.3.0 tibble_3.2.1
#> [9] ggplot2_3.4.4 tidyverse_2.0.0 cmdstanr_0.6.1 brms_2.20.4
#> [13] Rcpp_1.0.11 simstudy_0.7.0
#>
#> loaded via a namespace (and not attached):
#> [1] gridExtra_2.3 inline_0.3.19 rlang_1.1.1
#> [4] magrittr_2.0.3 matrixStats_1.0.0 compiler_4.3.1
#> [7] loo_2.6.0 callr_3.7.3 vctrs_0.6.4
#> [10] reshape2_1.4.4 pkgconfig_2.0.3 crayon_1.5.2
#> [13] fastmap_1.1.1 backports_1.4.1 ellipsis_0.3.2
#> [16] utf8_1.2.3 threejs_0.3.3 promises_1.2.1
#> [19] rmarkdown_2.25 tzdb_0.4.0 markdown_1.11
#> [22] ps_1.7.5 xfun_0.40 reprex_2.0.2
#> [25] jsonlite_1.8.7 later_1.3.1 styler_1.10.2
#> [28] parallel_4.3.1 prettyunits_1.2.0 R6_2.5.1
#> [31] dygraphs_1.1.1.6 stringi_1.7.12 StanHeaders_2.26.28
#> [34] rstan_2.32.3 knitr_1.44 zoo_1.8-12
#> [37] base64enc_0.1-3 R.utils_2.12.2 mvnfast_0.2.8
#> [40] bayesplot_1.10.0 timechange_0.2.0 httpuv_1.6.11
#> [43] Matrix_1.6-1.1 R.cache_0.16.0 igraph_1.5.1
#> [46] tidyselect_1.2.0 rstudioapi_0.15.0 abind_1.4-5
#> [49] yaml_2.3.7 codetools_0.2-19 miniUI_0.1.1.1
#> [52] curl_5.1.0 processx_3.8.2 pkgbuild_1.4.2
#> [55] lattice_0.21-9 plyr_1.8.9 shiny_1.7.5.1
#> [58] withr_2.5.1 bridgesampling_1.1-2 posterior_1.4.1
#> [61] coda_0.19-4 evaluate_0.22 RcppParallel_5.1.7
#> [64] xts_0.13.1 pillar_1.9.0 tensorA_0.36.2
#> [67] checkmate_2.2.0 DT_0.30 stats4_4.3.1
#> [70] shinyjs_2.1.0 distributional_0.3.2 generics_0.1.3
#> [73] hms_1.1.3 rstantools_2.3.1.1 munsell_0.5.0
#> [76] scales_1.2.1 gtools_3.9.4 xtable_1.8-4
#> [79] glue_1.6.2 tools_4.3.1 shinystan_2.6.0
#> [82] data.table_1.14.8 colourpicker_1.3.0 fs_1.6.3
#> [85] mvtnorm_1.2-3 grid_4.3.1 QuickJSR_1.0.7
#> [88] crosstalk_1.2.0 colorspace_2.1-0 nlme_3.1-163
#> [91] cli_3.6.1 fansi_1.0.5 Brobdingnag_1.2-9
#> [94] V8_4.4.0 gtable_0.3.4 R.methodsS3_1.8.2
#> [97] digest_0.6.33 htmlwidgets_1.6.2 farver_2.1.1
#> [100] htmltools_0.5.6.1 R.oo_1.25.0 lifecycle_1.0.3
#> [103] mime_0.12 shinythemes_1.2.0
```

^{Created on 2023-10-21 with reprex v2.0.2}

I’m not sure if I need to make an adjustment to the model, priors, and/or response distribution. In the models with my actual (non-simulated) data, I’ve also further increased adapt_delta, but no luck. I had similar convergence difficulties when fitting a model to my actual data using other response distributions (e.g., zero-inflated negative binomial, hurdle/zero-inflated poisson). I’m hopeful that working through and solving the issue with fitting the model to these simulated data will help me make the necessary adjustments to fit the model to my actual data. I’m a bit new to the hurdle negative binomial model, so any help and feedback would be greatly appreciated.