library(bellreg)
data(faults)
# ML approach:
mle <- bellreg(nf ~ lroll, data = faults, approach = "mle")
summary(mle)
#> Call:
#> bellreg(formula = nf ~ lroll, data = faults, approach = "mle")
#>
#> Coefficients:
#> Estimate StdErr z.value p.value
#> (Intercept) 0.98525432 0.33219450 2.9659 0.003018 **
#> lroll 0.00190932 0.00049003 3.8963 9.767e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> logLik = -88.96139 AIC = 181.9228
# Bayesian approach:
bayes <- bellreg(nf ~ lroll, data = faults, approach = "bayes", refresh = FALSE)
summary(bayes)
#>
#> bellreg(formula = nf ~ lroll, data = faults, approach = "bayes",
#> refresh = FALSE)
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> (Intercept) 0.976 0.007 0.337 0.310 0.754 0.976 1.204 1.633 2076.328 1.000
#> lroll 0.002 0.000 0.000 0.001 0.002 0.002 0.002 0.003 2350.826 1.001
#>
#> Inference for Stan model: bellreg.
#> 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 32 log-likelihood matrix.
#>
#> Estimate SE
#> elpd_loo -214.0 38.1
#> p_loo 71.7 22.5
#> looic 428.0 76.1
#> ------
#> 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) 18 56.2% 527
#> (0.7, 1] (bad) 8 25.0% <NA>
#> (1, Inf) (very bad) 6 18.8% <NA>
#> See help('pareto-k-diagnostic') for details.
loo::waic(log_lik)
#> Warning:
#> 21 (65.6%) p_waic estimates greater than 0.4. We recommend trying loo instead.
#>
#> Computed from 4000 by 32 log-likelihood matrix.
#>
#> Estimate SE
#> elpd_waic -209.3 37.5
#> p_waic 67.0 22.1
#> waic 418.6 75.0
#>
#> 21 (65.6%) p_waic estimates greater than 0.4. We recommend trying loo instead.