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.95263 0.84547 -2.3095 0.020914 *
#> smoker 2.17680 0.82368 2.6428 0.008223 **
#> gender -0.49596 0.42065 -1.1790 0.238393
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#>
#> Count regression coefficients:
#> Estimate StdErr z.value p.value
#> (Intercept) 0.716488 0.179872 3.9833 6.796e-05 ***
#> smoker -0.611678 0.183421 -3.3348 0.0008535 ***
#> gender 0.036384 0.177484 0.2050 0.8375724
#> ---
#> 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.158 0.007 0.319 -1.885 -1.337 -1.123 -0.938 -0.629 1901.386
#> Rhat
#> (Intercept) 1.002
#>
#> Count regression coefficients:
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff
#> (Intercept) 0.721 0.003 0.149 0.427 0.624 0.72 0.823 1.009 2935.825
#> smoker -1.079 0.003 0.145 -1.363 -1.176 -1.08 -0.979 -0.795 2437.865
#> gender 0.170 0.003 0.141 -0.110 0.078 0.17 0.264 0.441 2832.203
#> Rhat
#> (Intercept) 1.001
#> smoker 1.002
#> 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)
#> Warning: Relative effective sample sizes ('r_eff' argument) not specified.
#> For models fit with MCMC, the reported PSIS effective sample sizes and
#> MCSE estimates will be over-optimistic.
#>
#> Computed from 4000 by 511 log-likelihood matrix
#>
#> Estimate SE
#> elpd_loo -626.7 24.8
#> p_loo 4.4 0.4
#> looic 1253.4 49.6
#> ------
#> Monte Carlo SE of elpd_loo is 0.0.
#>
#> All Pareto k estimates are good (k < 0.5).
#> See help('pareto-k-diagnostic') for details.
loo::waic(log_lik)
#>
#> Computed from 4000 by 511 log-likelihood matrix
#>
#> Estimate SE
#> elpd_waic -626.7 24.8
#> p_waic 4.4 0.4
#> waic 1253.4 49.6