Model selection

npar

"npar" yields the number of model parameters, which is computed by the npar() method:

npar(model_train, model_train_sparse)
#> [1] 4 1

Here, model_train has 4 parameters (a coefficient for price, time, change, and comfort, respectively), while model_train_sparse has only a single price coefficient.

LL

If "LL" is included in criteria, model_selection() returns the model’s log-likelihood values. They can also be directly accessed via the logLik() method:³

logLik(model_train)
#> 'log Lik.' -1727.71 (df=4)
logLik(model_train_sparse)
#> 'log Lik.' -1865.86 (df=1)

AIC

Including "AIC" yields the Akaike’s Information Criterion (Akaike 1974), which is computed as −2 ⋅ LL + k ⋅ npar, where LL is the model’s log-likelihood value, k is the penalty per parameter with k = 2 per default for the classical AIC, and npar is the number of parameters in the fitted model.

Alternatively, the AIC() method also returns the AIC values:

AIC(model_train, model_train_sparse, k = 2)
#>                    df      AIC
#> model_train         4 3463.419
#> model_train_sparse  1 3733.720

The AIC quantifies the trade-off between over- and under-fitting, where smaller values are preferred. Here, the increase in goodness of fit justifies the additional 3 parameters of model_train.

BIC

Similar to the AIC, "BIC" yields the Bayesian Information Criterion (Schwarz 1978), which is defined as −2 ⋅ LL + log (nobs) ⋅ npar, where LL is the model’s log-likelihood value, nobs is the number of data points (which can be accessed via the nobs() method), and npar is the number of parameters in the fitted model. The interpretation is analogue to AIC.

{RprobitB} also provided a method for the BIC value:

BIC(model_train, model_train_sparse)
#>                    df      BIC
#> model_train         4 3487.349
#> model_train_sparse  1 3739.702

WAIC (with se(WAIC) and pWAIC)

WAIC is short for Widely Applicable (or Watanabe-Akaike) Information Criterion (Watanabe and Opper 2010). As for AIC and BIC, the smaller the WAIC value the better the model. Including "WAIC" in criteria yields the WAIC value, its standard error se(WAIC), and the effective number of parameters pWAIC, see below.

The WAIC is defined as

−2 ⋅ lppd + 2 ⋅ p_WAIC,

where lppd stands for log pointwise predictive density and p_WAIC is a penalty term proportional to the variance in the posterior distribution that is sometimes called effective number of parameters, see McElreath (2020) p. 220 for a reference.

The lppd is approximated as follows. Let p_si = Pr (y_i ∣ θ_s) be the probability of observation y_i (here the single choices) given the s-th set θ_s of parameter samples from the posterior. Then

lppd = ∑_ilog (S⁻¹∑_sp_si). The penalty term is computed as the sum over the variances in log-probability for each observation: p_WAIC = ∑_i𝕍_θlog (p_si). The WAIC has a standard error of $$\sqrt{n \cdot \mathbb{V}_i \left[-2 \left(\text{lppd} - \mathbb{V}_{\theta} \log (p_{si}) \right)\right]},$$ where n is the number of choices.

Before computing the WAIC of an object, the probabilities p_si must be computed via the compute_p_si() function: ⁴

model_train <- compute_p_si(model_train)

Afterwards, the WAIC can be accessed as follows:

WAIC(model_train)
WAIC(model_train_sparse)

MMLL

"MMLL" in criteria stands for marginal model log-likelihood. The model’s marginal likelihood Pr (y ∣ M) for a model M and data y is required for the computation of Bayes factors, see below. In general, the term has no closed form and must be approximated numerically.

{RprobitB} uses the posterior Gibbs samples derived from the mcmc() function to approximate the model’s marginal likelihood via the posterior harmonic mean estimator (Newton and Raftery 1994): Let S denote the number of available posterior samples θ₁, …, θ_S. Then, $$\Pr(y\mid M) = \left(\mathbb{E}_\text{posterior} 1/\Pr(y\mid \theta,M) \right)^{-1} \approx \left( \frac{1}{S} \sum_s 1/\Pr(y\mid \theta_s,M) \right) ^{-1} = \tilde{\Pr}(y\mid M).$$

By the law of large numbers, $\tilde{\Pr}(y\mid M) \to \Pr(y\mid M)$ almost surely as S → ∞.

As for the WAIC, computing the MMLL relies on the probabilities p_si = Pr (y_i ∣ θ_s), which must first be computed via the compute_p_si() function. Afterwards, the mml() function can be called with an RprobitB_fit object as input. The function returns the RprobitB_fit object, where the marginal likelihood value is saved as the entry "mml" and the marginal log-likelihood value as the attribute "mmll":⁵

model_train <- mml(model_train)
model_train$mml
attr(model_train$mml, "mmll")

There are two options for improving the approximation: As for the WAIC, you can use more posterior samples. Alternatively, you can combine the posterior harmonic mean estimate with the prior arithmetic mean estimator (Hammersley and Handscomb 1964): For this approach, S samples θ₁, …, θ_S are drawn from the model’s prior distribution. Then,

$$\Pr(y\mid M) = \mathbb{E}_\text{prior} \Pr(y\mid \theta,M) \approx \frac{1}{S} \sum_s \Pr(y\mid \theta_s,M) = \tilde{\Pr}(y\mid M).$$

Again, it hols by the law of large numbers, that $\tilde{\Pr}(y\mid M) \to \Pr(y\mid M)$ almost surely as S → ∞. The final approximation of the model’s marginal likelihood than is a weighted sum of the posterior harmonic mean estimate and the prior arithmetic mean estimate, where the weights are determined by the sample sizes.

To use the prior arithmetic mean estimator, call the mml() function with a specification of the number of prior draws S and set recompute = TRUE:⁶

model_train <- mml(model_train, S = 1000, recompute = TRUE)

Note that the prior arithmetic mean estimator works well if the prior and posterior distribution have a similar shape and strong overlap, as Gronau et al. (2017) points out. Otherwise, most of the sampled prior values result in a likelihood value close to zero, thereby contributing only marginally to the approximation. In this case, a very large number S of prior samples is required.

BF

The Bayes factor is an index of relative posterior model plausibility of one model over another (Marin and Robert 2014). Given data y and two models mod0 and mod1, it is defined as

$$ BF(\texttt{mod0},\texttt{mod1}) = \frac{\Pr(\texttt{mod0} \mid \texttt{y})}{\Pr(\texttt{mod1} \mid \texttt{y})} = \frac{\Pr(\texttt{y} \mid \texttt{mod0} )}{\Pr(\texttt{y} \mid \texttt{mod1})} / \frac{\Pr(\texttt{mod0})}{\Pr(\texttt{mod1})}. $$

The ratio Pr (mod0)/Pr (mod1) expresses the factor by which mod0 a priori is assumed to be the correct model. Per default, {RprobitB} sets Pr (mod0) = Pr (mod1) = 0.5. The front part Pr (y ∣ mod0)/Pr (y ∣ mod1) is the ratio of the marginal model likelihoods. A value of BF(mod0, mod1) > 1 means that the model mod0 is more strongly supported by the data under consideration than mod1.

Adding "BF" to the criteria argument of model_selection yields the Bayes factors. We find a decisive evidence (Jeffreys 1998) in favor of model_train.

model_selection(model_train, model_train_sparse, criteria = c("BF"))

pred_acc

Finally, adding "pred_acc" to the criteria argument for the model_selection() function returns the share of correctly predicted choices. From the output of model_selection() above (or alternatively the one in the following) we deduce that model_train correctly predicts about 6% of the choices more than model_train_sparse:⁷

pred_acc(model_train, model_train_sparse)
#> [1] 0.6951178 0.6340048