# 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.95021  0.84291 -2.3137 0.020686 *
#> smoker       2.17462  0.82118  2.6482 0.008093 **
#> gender      -0.49573  0.42045 -1.1791 0.238374
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#>
#> Count regression coefficients:
#>              Estimate    StdErr z.value   p.value
#> (Intercept)  0.716757  0.179792  3.9866 6.703e-05 ***
#> smoker      -0.611758  0.183354 -3.3365 0.0008484 ***
#> gender       0.036254  0.177468  0.2043 0.8381297
#> ---
#> 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.168   0.008 0.342 -1.937 -1.354 -1.127 -0.932 -0.636 1679.041
#>             Rhat
#> (Intercept)    1
#>
#> Count regression coefficients:
#>               mean se_mean    sd   2.5%    25%    50%    75%  97.5%    n_eff
#> (Intercept)  0.718   0.003 0.145  0.433  0.619  0.718  0.818  1.003 2699.184
#> smoker      -1.078   0.003 0.145 -1.363 -1.175 -1.077 -0.983 -0.792 2377.645
#> gender       0.173   0.003 0.142 -0.107  0.077  0.173  0.274  0.448 2809.804
#>              Rhat
#> (Intercept) 1.001
#> smoker      1.000
#> gender      1.001
#> ---
#> 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)
#>
#> Computed from 4000 by 511 log-likelihood matrix.
#>
#>          Estimate   SE
#> elpd_loo   -626.8 24.8
#> p_loo         4.5  0.4
#> looic      1253.6 49.5
#> ------
#> MCSE of elpd_loo is 0.0.
#> MCSE and ESS estimates assume independent draws (r_eff=1).
#>
#> All Pareto k estimates are good (k < 0.7).
#> See help('pareto-k-diagnostic') for details.
loo::waic(log_lik)
#>
#> Computed from 4000 by 511 log-likelihood matrix.
#>
#>           Estimate   SE
#> elpd_waic   -626.8 24.8
#> p_waic         4.5  0.4
#> waic        1253.6 49.5