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.