Zero-inflated Bell model

library(bellreg)

data(cells)

# ML approach:
mle <- zibellreg(cells ~ smoker+gender|smoker+gender, data = cells, approach = "mle")
summary(mle)
#> Call:
#> zibellreg(formula = cells ~ smoker + gender | smoker + gender, 
#>     data = cells, approach = "mle")
#> 
#> Zero-inflated regression coefficients:
#>             Estimate   StdErr z.value  p.value   
#> (Intercept) -1.95238  0.84517 -2.3100 0.020886 * 
#> smoker       2.17674  0.82339  2.6436 0.008202 **
#> gender      -0.49592  0.42059 -1.1791 0.238354   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> 
#> Count regression coefficients:
#>              Estimate    StdErr z.value   p.value    
#> (Intercept)  0.716517  0.179858  3.9838 6.783e-05 ***
#> smoker      -0.611658  0.183407 -3.3350 0.0008531 ***
#> gender       0.036358  0.177476  0.2049 0.8376820    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> --- 
#> logLik = -610.3234   AIC = 1232.647

# Bayesian approach:
bayes <- zibellreg(cells ~ 1|smoker+gender, data = cells, approach = "bayes", refresh = FALSE)
summary(bayes)
#> Call:
#> zibellreg(formula = cells ~ 1 | smoker + gender, data = cells, 
#>     approach = "bayes", refresh = FALSE)
#> 
#> Zero-inflated regression coefficients:
#>               mean se_mean    sd   2.5%    25%    50%    75%  97.5%    n_eff
#> (Intercept) -1.162   0.008 0.336 -1.934 -1.344 -1.122 -0.932 -0.637 1785.709
#>              Rhat
#> (Intercept) 1.003
#> 
#> Count regression coefficients:
#>               mean se_mean    sd   2.5%    25%    50%    75%  97.5%    n_eff
#> (Intercept)  0.717   0.003 0.151  0.423  0.614  0.715  0.818  1.021 2920.121
#> smoker      -1.072   0.003 0.148 -1.355 -1.171 -1.074 -0.976 -0.775 2482.000
#> gender       0.170   0.003 0.146 -0.124  0.077  0.170  0.270  0.456 2801.897
#>              Rhat
#> (Intercept) 1.000
#> smoker      1.000
#> gender      0.999
#> --- 
#> Inference for Stan model: zibellreg.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
#> post-warmup draws per chain=1000, total post-warmup draws=4000.

log_lik <- loo::extract_log_lik(bayes$fit)
loo::loo(log_lik)
#> Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
#> 
#> Computed from 4000 by 511 log-likelihood matrix.
#> 
#>          Estimate    SE
#> elpd_loo  -1064.4  53.3
#> p_loo       220.5  22.0
#> looic      2128.8 106.6
#> ------
#> MCSE of elpd_loo is NA.
#> MCSE and ESS estimates assume independent draws (r_eff=1).
#> 
#> Pareto k diagnostic values:
#>                          Count Pct.    Min. ESS
#> (-Inf, 0.7]   (good)     432   84.5%   165     
#>    (0.7, 1]   (bad)       57   11.2%   <NA>    
#>    (1, Inf)   (very bad)  22    4.3%   <NA>    
#> See help('pareto-k-diagnostic') for details.
loo::waic(log_lik)
#> Warning: 
#> 93 (18.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.
#> 
#> Computed from 4000 by 511 log-likelihood matrix.
#> 
#>           Estimate   SE
#> elpd_waic  -1008.7 48.3
#> p_waic       164.8 15.3
#> waic        2017.4 96.6
#> 
#> 93 (18.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.