Incorporating Model Uncertainty in Market Response Models with Multiple Endogenous Variables by Bayesian Model Averaging

Abstract

We develop a method to incorporate model uncertainty by model averaging in generalized linear models subject to multiple endogeneity and instrumentation. Our approach builds on a Gibbs sampler for the instrumental variable framework that incorporates model uncertainty in both outcome and instrumentation stages. Direct evaluation of model probabilities is intractable in this setting. However, we show that by nesting model moves inside the Gibbs sampler, a model comparison can be performed via conditional Bayes factors, leading to straightforward calculations. This new Gibbs sampler is slightly more involved than the original algorithm and exhibits no evidence of mixing difficulties. We further show how the same principle may be employed to evaluate the validity of instrumentation choices. We conclude with an empirical marketing study: estimating opening box office by three endogenous regressors (prerelease advertising, opening screens, and production budget).

1. Introduction

Market response modeling focuses on estimating the effects of marketing activities on performance. However, marketing managers are often strategic in their use of marketing activities and adapt them in response to factors unobserved by the researcher [1, 2, 3]. Endogeneity arises, for example, when a firm’s marketing strategies such as advertising spending, channel selection, and pricing are nonrandom and influenced by the firm- and industry-level factors [4, 5, 6]. Strategic management decisions are endogenous to their expected effects on market performance. Therefore, empirical market response models that seek to estimate the causal effect of multiple marketing instruments need to account for such strategic planning of marketing activities, or otherwise may suffer from an endogeneity problem, leading to biased estimates of the effects of the marketing activities on performance [1, 3, 4, 7]. Dealing with endogeneity has been extensively discussed in the marketing literature, especially concerning different forms of regression and panel models [1, 5, 8, 9, 10], choice models [11, 12], endogeneity correction based on a control function approach [13, 14], as well as structural equations models [4]. However, little research addresses incorporating model uncertainty related to endogeneity in generalized linear models.

We consider the problem of incorporating instruments and covariate uncertainty into the Bayesian estimation of an instrumental variable (IV) regression system. The concepts of model uncertainty and model averaging have received widespread attention in the economics literature for the standard linear regression framework [15, 16, 17, 18] and in generalized linear models [19, 20, 21, 22]. For a good introduction to Bayesian model averaging (BMA), see [23]. Primarily, these frameworks do not directly address the case of multiple endogenous variables, and only recently has attention been paid to model uncertainty involving multiple endogenous variables. Unfortunately, the nested nature of IV estimation renders direct model comparison difficult. In the economics literature, this has led to several different approaches [24, 25]. Durlauf et al. [25] consider approximations of marginal likelihoods in a framework similar to two-stage least squares. Lenkoski et al. [16] continue this development with the two-stage Bayesian model averaging (2SBMA) methodology, which uses a framework developed by Kleibergen and Zivot [26] to propose a two-stage extension of the unit information prior [27]. Similar approaches in closely related models have been developed by [15, 28].

Koop et al. [29] developed a fully Bayesian methodology that does not utilize approximations to integrated likelihoods. They present a reversible jump Markov chain Monte Carlo (RJMCMC) algorithm [30], which extends the methodology of Holmes et al. [31]. The authors then show that the method can handle a variety of priors, including those of [32, 33], and [34]. However, the authors note that the direct application of RJMCMC leads to significant mixing difficulties and relies on a complicated model move procedure similar to simulated tempering to escape local model modes. There is a more straightforward and relatively general model search procedure. Madigan and York [35] proposed the Markov Chain Monte Carlo Model Composition (MC3) in which one applies the same idea of a Metropolis-Hastings step for model jumps from RJMCMC but in a simplified fashion.

We propose an alternative solution to this problem: Instrumental Variable Bayesian Model Averaging (IVBMA). Our method builds on a Gibbs sampler for the IV framework, extended from that discussed in Rossi et al. [36]. While direct model comparisons are intractable, we introduce the notion of a conditional Bayes factor (CBF), first discussed by Dickey and Gunel [37] and employed in a seemingly unrelated regression context by [31]. The CBF compares two models in a nested hierarchical system, conditional on parameters not influenced by the models under consideration. We show that the CBF for both outcome and instrumental equations is exceedingly straightforward to calculate and essentially reduces to the normalizing constants of a multivariate normal distribution.

Further, we note that our method can handle generalized linear mixed models with multiple endogenous variables in a straightforward fashion. This leads to a procedure in which model moves are embedded in a Gibbs sampler, which we term MC3-within-Gibbs. Based on this order of operations, IVBMA is only trivially more complicated than a Gibbs sampler that does not incorporate model uncertainty and thus appears to have limited issues regarding mixing. This feature is essential as it shows more complicated scenarios involving endogeneity, instrumentation, and model uncertainty can be handled within this framework, an important feature when constructing more involved Bayesian hierarchical models.

When working with a large system of equations subject to endogeneity and instrumentation, there is a natural concern that the instrument assumptions may not hold. A host of frequentist-type hypotheses has been proposed to examine the instrument conditions; the most familiar to applied researchers is the test of Sargan [38]. There have been, to our knowledge, no similar checks of instrument validity proposed in the Bayesian IV literature outside of the approximate method advocated in [16]. We offer a new verification of instrument validity, also based on CBFs, which appears to be the Bayesian analog of the Sargan test. This method can integrate seamlessly with the IVBMA framework and offers a check of instrument validity.

The article proceeds as follows. The basic framework we consider and the Gibbs sampler ignoring model uncertainty is discussed in Section 2. Section 3 reviews the concept of model uncertainty, introduces the notion of CBFs, and derives the conditional model probabilities used by IVBMA. In Section 4, we propose our method of assessing instrument validity. Section 5 presents empirical illustrations of the proposed model for predicting box office revenues. Lastly, we summarize and conclude with potential applications of the IVBMA approach.

2. The instrumental variable model with multiple endogenous variables

We consider the following classic linear system model with multiple endogenous variables:

$$Y_{\mathit{ir}}={\mathit{U}}_i^{{\left(r\right)}^{\prime }}{\mathit{\beta }}_r+{\epsilon }_{\mathit{ir}},$$E1

where $r\in \left\{1,\dots ,R\right\}$ denotes the R equations in the system and $i\in \left\{1,\dots ,n\right\}$ a set of iid observations. Throughout, we assume that $Y_{i1}$ represents the dependent outcome of interest and ($Y_{i2},\dots ,Y_{\mathit{iR}})$ represent endogenously determining variables for observation i. Thus, each covariate vector ${\mathit{U}}_i^{\left(r\right)}$ has length $p_r$ and is formed such that

$${\mathit{U}}_i^{\left(1\right)}={\left(Y_{i1}\dots Y_{\mathit{iR}}{W}_{i1}\dots W_{\mathit{iq}}\right)}^{\prime .},$$

while

$${\mathit{U}}_i^{\left(r\right)}={\left(Z_{i1}\dots Z_{\mathit{\text{is}}}{W}_{i1}\dots W_{\mathit{iq}}\right)}^{\prime },$$

for r > 1. Letting ${\epsilon }_i={\left({\epsilon }_{i1},\dots ,{\epsilon }_{\mathit{iR}}\right)}^{\prime }$, we assume

$${\mathit{\epsilon }}_i~{\mathcal{N}}_R\left(0,{\mathit{K}}^{-1}\right).$$E2

When $K_{1r}\ne 0$ for a given r > 1, this implies a lack of conditional independence between the residuals for the response and the associated endogenous variable. This contaminates inference if unaccounted for, necessitating the existence of instruments ${\mathit{Z}}_i$ that do not appear in ${\mathit{U}}_i^{\left(1\right)}$ and joint estimation of the parameters in Eq. (1) and Eq. (2).

Generalized linear mixed models provide a unified approach that directly acknowledges multiple levels of dependency and model different data types [39, 40, 41, 42]. Extensions to generalized linear models implicitly assume a continuous response with Gaussian errors. Extending these developments to alternative sampling models is straightforward in the context of a random-effects framework. Let g a link function such that for the response $Y_i$,

$$E\left[Y_{i1}\right]=g^{-1}\left({\mathit{U}}_i^{{\left(1\right)}^{\prime }}{\mathit{\beta }}_1+{\epsilon }_{i1}\right),$$E3

while the remaining $Y_{\mathit{ir}}$ have forms given by Eq. (1), and the residual vector remains distributed according to a $\mathcal{N}\left(0,{\mathit{K}}^{-1}\right)$ distribution. Below we first develop the normal IVBMA with an identity link.

We proceed by discussing the Bayesian estimation of these parameters under standard conjugate priors, following the developments of [36]. Accordingly, with each parameter vector, we assume

$${\mathit{\beta }}_r~\mathcal{N}\left(0,{\mathbb{I}}_{p_r}\right),$$

and

$$\mathit{K}~\mathcal{W}\left(3,{\mathbb{I}}_R\right)$$

where $\mathit{K}~\mathcal{W}\left(\delta ,\mathit{D}\right)$ represents a Wishart distribution with density

$$\mathit{pr}\left(\mathit{K}|\delta ,\mathit{D}\right)\propto {\left|\mathit{K}\right|}^{\frac{\delta -2}{2}}\mathit{exp}\left(-\frac{1}{2}\mathit{tr}\left(\mathit{KD}\right)\right){\mathbf{1}}_{K\in {\mathbb{P}}_R},$$

where ${\mathbb{P}}_R$ is the cone of symmetric positive definite matrices.

Let $\mathit{\theta }=\left\{{\mathit{\beta }}_1,\dots ,{\mathit{\beta }}_R,\mathit{K}\right\}$ represent the collection of parameters to be estimated. Denote the data $\mathcal{D}=\left\{\mathit{Y},{\mathit{U}}^{\left(1\right)},\dots ,{\mathit{U}}^{\left(R\right)}\right\}$, where $\mathit{Y}$ is the $n\times R$ matrix of responses and endogenous variables, and each ${\mathit{U}}^{\left(r\right)}$ is a $n\times p_r$ matrix. Then, our goal is to determine the posterior distribution $\mathit{pr}\left(\mathit{\theta }|\mathcal{D}\right)$. Rossi et al. [36] discuss the estimation of this model for the case when R = 2 and note that it is not possible to evaluate this posterior directly. However, an approximate inference may be performed via Gibbs sampling.

Fix r and suppose that $\mathit{K}$ and all ${\mathit{\beta }}_t$ for $t\ne r$ are given. Note, by properties of standard normal variates that

$${\epsilon }_{\mathit{ir}}\mid \mathit{K},{\left\{{\mathit{\beta }}_t\right\}}_{t\ne r}~\mathcal{N}\left({\mu }_{\mathit{ir}},K_{\mathit{rr}}^{-1}\right),$$

where

$${\mu }_{\mathit{ir}}=-\sum _{t\ne r}\frac{K_{\mathit{rt}}}{K_{\mathit{rr}}}\left(Y_{\mathit{it}}-{\mathit{U}}_i^{\left(t\right)}{\mathit{\beta }}_t\right).$$

Set ${\tilde{Y}}_{\mathit{ir}}=Y_{\mathit{ir}}-{\mu }_{\mathit{ir}}$ and thus note that

$${\tilde{Y}}_{\mathit{ir}}~\mathcal{N}\left({\mathit{U}}_i^{\left(r\right)}{\mathit{\beta }}_r,K_{\mathit{rr}}^{-1}\right).$$

The act of conditioning, therefore, turns the original system into a simple linear regression problem, and via standard results, we have that

$${\mathit{\beta }}_r\mid \mathit{K},{\left\{{\mathit{\beta }}_t\right\}}_{t\ne r}~\mathcal{N}\left({\widehat{\mathit{\beta }}}_r,{\mathit{\Omega }}_r^{-1}\right)$$E4

where

$${\mathit{\Omega }}_r=K_{\mathit{rr}}{\mathit{U}}^{{\left(r\right)}^{\prime }}{\mathit{U}}^{\left(r\right)}+{\mathbb{I}}_{p_r},$$

$${\widehat{\mathit{\beta }}}_r=K_{\mathit{rr}}{\mathit{\Omega }}_r^{-1}{\mathit{U}}^{{\left(r\right)}^{\prime }}{\tilde{\mathit{Y}}}_r.$$

Finally, suppose that all ${\mathit{\beta }}_r$ are given, then

$$\mathit{K}~\mathcal{W}\left(\delta +n,\mathit{E}+{\mathbb{I}}_R\right),$$E5

where

$$\mathit{E}={\sum }_{i=1}^n{\mathit{\epsilon }}_i{\mathit{\epsilon }}_{i^{\prime }},$$

with each ${\mathit{\epsilon }}_i$ computed relative to the current state of ${\mathit{\beta }}_1,\dots ,{\mathit{\beta }}_R$.

Eq. (4) and Eq. (5) thereby give the full conditionals necessary for the Gibbs sampler. For a basic introduction to MCMC sampling with illustration, see [43]. Our approach differs slightly from that of Rossi et al. [36], in that their Gibbs sampler features a more involved manner of updating the instrumental covariates ${\mathit{\beta }}_2$. Though the two approaches evaluate the same posterior distribution, the application of [36] when $\mathrm{R}\ge 3$ is not straightforward, and it only applies to a linear regression model. Therefore, we find that the above approach leads to more coherent implementation and description, and therefore prefer it to that of [36] for the generalized linear models with multiple endogenous variables.

For a Poisson regression using a log link in Eq. (3), the term ${\epsilon }_{i1}$ is no longer observable and is often referred to as a Poisson random effect model [41, 44, 45]. However, in a Gibbs sampling framework, these factors may be incorporated into additional parameters to be determined in the posterior. Appendix 2 shows how MCMC methods can be implemented when $Y_{i1}$ in (3) has a Poisson likelihood.

3. Incorporating model uncertainty

We outline our method for incorporating model uncertainty into the framework in Eq. (1) and Eq. (2). To explain the motivation behind our CBF approach, we first review a classic Bayesian model selection method. We then show how the concept of Bayes Factors can be usefully embedded in a Gibbs sampler yielding CBFs. These CBFs are then shown to yield straightforward calculations.

3.1 Model selection and Bayes factors

In a general framework, incorporating model uncertainty involves considering a collection of candidate models $\mathcal{I}$, using the data $\mathcal{D}$. Each model $I$ consists of a collection of probability distributions for the data $\mathcal{D}$, $\left\{\mathit{pr}\left(\mathcal{D}|\psi \right),\psi \in {\mathrm{\Psi }}_I\right\}$ where ${\mathrm{\Psi }}_I$ denotes the parameter space for the parameters of model $I$ and is a subset of the full parameter space $\mathrm{\Psi }$.

By letting the model become an additional parameter to be assessed in the posterior, we aim to calculate the posterior model probabilities given the data $\mathcal{D}$. By Bayes’ rule

$$\mathit{pr}\left(I|\mathcal{D}\right)=\frac{\mathit{pr}\left(\mathcal{D}|I\right)\mathit{pr}\left(I\right)}{{\sum }_{I\prime \in \mathcal{I}}‍\mathit{pr}\left(\mathcal{D}|I^{\prime }\right)\mathit{pr}\left(I^{\prime }\right)}$$E6

where $\mathit{pr}\left(I\right)$, denotes the prior probability for model $I\in \mathcal{I}$.

The integrated likelihood $\mathit{pr}\left(\mathcal{D}|I\right)$, is defined by

$$\mathit{pr}\left(\mathcal{D}|I\right)={\int }_{{\mathrm{\Psi }}_I}‍\mathit{pr}\left(\mathcal{D}|\psi \right)\mathit{pr}\left(\psi |I\right)\mathit{d\psi }$$E7

where $\mathit{pr}\left(\psi |I\right)$ is the prior for $\psi$ under model $I$, which by definition has all its mass on ${\mathrm{\Psi }}_I$.

One possibility for pairwise comparison of models is offered by the Bayes factor (BF), which is in most cases defined together with the posterior odds [22, 46]. The posterior odds of model $I$ versus model $I^{\prime }$ are given by

$$\frac{\mathit{pr}\left(I|\mathcal{D}\right)}{\mathit{pr}\left(I^{\prime }|\mathcal{D}\right)}=\frac{\mathit{pr}\left(\mathcal{D}|I\right)}{\mathit{pr}\left(\mathcal{D}|I^{\prime }\right)}\frac{\mathit{pr}\left(I\right)}{\mathit{pr}\left(I^{\prime }\right)},$$

where

$$\frac{\mathit{pr}\left(\mathcal{D}|I\right)}{\mathit{pr}\left(\mathcal{D}|I^{\prime }\right)}\text{and}\frac{\mathit{pr}\left(I\right)}{\mathit{pr}\left(I^{\prime }\right)}$$

denote the Bayes factor and the prior odds of $I$ versus $I^{\prime }$, respectively.

When the integrated likelihood in Eq. (7) and thus the BF can be computed directly, a straightforward method for exploring the model space, Markov Chain Monte Carlo Model Composition (MC3), was developed by Madigan and York [35]. MC3 determines posterior model probabilities by generating a stochastic process that moves through the model space $\mathcal{I}$ and has equilibrium distribution $\mathit{pr}\left(I|\mathcal{D}\right)$. Given the current state $I^{\left(s\right)}$, MC3 (a) proposes a new model $I^{\prime }$ according to a proposal distribution $q\left(·|·\right)$, (b) calculates

$$\alpha =\frac{\mathit{pr}\left(\mathcal{D}|I^{\prime }\right)\mathit{pr}\left(I^{\prime }\right)q\left(I^{\left(s\right)}|I^{\prime }\right)}{\mathit{pr}\left(\mathcal{D}|I^{\left(s\right)}\right)\mathit{pr}\left(I^{\left(s\right)}\right)q\left(I^{\prime }|I^{\left(s\right)}\right)},$$

and (c) sets $I^{\left(s+1\right)}=I^{\prime }$ with probability $min\left\{\alpha ,1\right\}$, otherwise setting $I^{\left(s+1\right)}=I^{\left(s\right)}$. It is important to note that moving between models via the MC3 approach constitutes a valid MCMC transition. This feature is critical in the development below, in that MC3 moves may be nested inside larger structures in a manner similar to Gibbs updates.

3.2 Model determination

Incorporating model uncertainty into the system Eq. (1) involves considering a separate model space ${\mathcal{M}}_r$ for each equation in the system. A given model $M_r\in {\mathcal{M}}_r$ thus restricts certain elements of ${\mathit{\beta }}_r$ to zero, and we write ${\mathit{\beta }}_{M_r}$ to indicate the non-zero elements according to $M_r$. We further let ${\Lambda }_{M_r}$ be the subspace of ${ℝ}^{p_r}$ spanned by ${\mathit{\beta }}_{M_r}$. Ideally, we would be able to incorporate model uncertainty into this system in a manner analogous to that described above. Unfortunately, the following cannot be directly calculated in any discernible way.

$$\mathit{pr}\left(\mathcal{D}|M_1,\dots M_R\right)={\int }_{{\mathbb{P}}_R}{\int }_{{\Lambda }_{M_1}}\cdots {\int }_{{\Lambda }_{M_R}}\mathit{pr}\left(\mathcal{D}|{\left\{{\mathit{\beta }}_{M_r}\right\}}_{r=1}^R,\mathit{K}\right)\mathit{pr}\left(\mathit{K}\right){\prod }_{r=1}^R\mathit{pr}\left({\mathit{\beta }}_{M_r}\right)d{\mathit{\beta }}_{M_1}\cdots d{\mathit{\beta }}_{M_R}d\mathit{K}$$

Therefore, an implementation of MC3 in the product space of ${\mathcal{M}}_1\times \cdots \times {\mathcal{M}}_R$ is infeasible. What we show below, however, is that embedding MC3 within the Gibbs sampler, and therefore calculation using CBFs to move between models offers an extremely efficient solution. CBFs were initially discussed in Dickey and Gunel [37] in a different context.

Given the system Eq. (1), fix r and suppose that ${\mathit{\theta }}_{-r}=\left\{\mathit{K},{\left\{{\mathit{\beta }}_t\right\}}_{t\ne r}\right\}$ is given. Now consider comparing two models $M_r,L_r\in {\mathcal{M}}_r$. Finally, suppose that the prior over models ${\mathcal{M}}_r$ is set independent of ${\mathit{\theta }}_{-r}$. We then have

$$\frac{\mathit{pr}\left(M_r|\mathcal{D},{\mathit{\theta }}_{-r}\right)}{\mathit{pr}\left(L_r|\mathcal{D},{\mathit{\theta }}_{-r}\right)}=\frac{\mathit{pr}\left(\mathcal{D}|M_r,{\mathit{\theta }}_{-r}\right)}{\mathit{pr}\left(\mathcal{D}|L_r,{\mathit{\theta }}_{-r}\right)}\times \frac{\mathit{pr}\left(M_r\right)}{\mathit{pr}\left(L_r\right)}.$$E8

Thus, the conditional posterior odds depend on calculating a Bayes factor conditional on the current state of ${\mathit{\theta }}_{-r}$.

Calculating the relevant terms in Eq. (6) is straightforward. In particular, we note that

$$\mathit{pr}\left(\mathcal{D}|M_r,{\mathit{\theta }}_{-r}\right)={\int }_{{\Lambda }_{M_r}}\mathit{pr}\left(\mathcal{D}|{\mathit{\beta }}_{M_r},{\mathit{\theta }}_{-r}\right)\mathit{pr}\left({\mathit{\beta }}_{M_r}|M_r\right)d{\mathit{\beta }}_{M_r},$$

which is, in essence, an integrated likelihood for model $M_r$ conditional on fixed values of ${\mathit{\theta }}_{-r}$. In Appendix 1, we show that

$${\int }_{{\Lambda }_{M_r}}\mathit{pr}\left(\mathcal{D}|{\mathit{\beta }}_{M_r},{\mathit{\theta }}_{-r}\right)d{\beta }_{M_r}\propto {\left|{\mathit{\Omega }}_{M_r}\right|}^{-1/2}\mathit{exp}\left(\frac{1}{2}{{\widehat{\mathit{\beta }}}_{M_r}}^{\prime }{\mathit{\Omega }}_{M_r}{\widehat{\mathit{\beta }}}_{M_r}\right)$$E9

where ${\widehat{\mathit{\beta }}}_{M_r}$ and ${\mathit{\Omega }}_{M_r}$ relative to the subspace ${\Lambda }_{M_r}$.

The power of this result is that the model $M_r$ and the associated parameter ${\mathit{\beta }}_{M_r}$ may then be updated in a block. In particular, we note that

$$\mathit{pr}\left({\mathit{\beta }}_r,M_r|{\mathit{\theta }}_{-r},\mathcal{D}\right)=\mathit{pr}\left({\mathit{\beta }}_r|M_r,{\mathit{\theta }}_{-r},\mathcal{D}\right)\times \mathit{pr}\left(M_r|{\mathit{\theta }}_{-r},\mathcal{D}\right).$$

Since MC3 constitutes a valid MCMC transition in the model space ${\mathcal{M}}_r$, we may first attempt to update $M_r$ via Eq. (8) and then subsequently resample ${\beta }_{M_r}$ via Eq. (4). By cycling through all R equations in Eq. (1) in this manner, and then subsequently updating $\mathit{K}$, we have proposed a computationally efficient estimation strategy for incorporating model uncertainty in IV frameworks.

4. Assessing instrument validity

For the estimates ${\mathit{\beta }}_1$ to have appropriate inferential properties, it is critical that the instrumental variables $\mathit{Z}$ be valid. In other words, $E\left[{{\mathit{Z}}_i}^{\prime }{\epsilon }_{i1}|{\epsilon }_{i2},\dots ,{\epsilon }_{\mathit{iR}}\right]=\mathbf{0}$. Many tools exist for evaluating the validity of this assumption in frequentist settings, and the most popular method is the test of Sargan [38]. To our knowledge, consideration of similar assessments in a Bayesian framework has not been explored beyond the approximate analysis proposed in [16]. We offer a Bayesian evaluation of instrument validity, borrowing many of the ideas above and merging them with the idea of the Sargan test.

Suppose that all residuals were known. Let $\mathit{\varsigma }$ be such that

$${\varsigma }_i={\epsilon }_{i1}+\sum _{r=2}^R\frac{K_{1r}}{K_{11}}{\epsilon }_{\mathit{ir}}.$$

The essential notion of the Sargan test is to consider the model,

$${\varsigma }_i=\mathit{Z}{\prime }_i\mathit{\xi }+{\eta }_i,{\eta }_i~\mathcal{N}\left(0,{\tau }^{-1}\right)$$

and test whether $\mathit{\xi }\ne \mathbf{0}$. The mechanics of the Sargan test ultimately rely on asymptotic theory, and Lenkoski et al. [16] discuss its poor performance in low sample size environments.

Our approach is to model this in a Bayesian context. In particular, we consider two models: $J_0$ which states that $\mathit{\xi }=\mathbf{0}$ and $J_1$ which puts $\mathit{\xi }\in {ℝ}^q$. We then aim to determine whether $\mathit{pr}\left(J_0|\mathcal{D}\right)$ is large, indicating instrument validity. Note that this can be represented as the following marginalization

$$\mathit{pr}\left(J_0|\mathcal{D}\right)=\int \mathit{pr}\left(J_0|\mathit{\varsigma },\mathcal{D}\right)\mathit{pr}\left(\mathit{\varsigma }|\mathcal{D}\right)d\mathit{\varsigma }.$$E10

Let $\left\{{\mathit{\theta }}^{\left(1\right)},\dots ,{\mathit{\theta }}^{\left(S\right)}\right\}$ be an MCMC sample of $\mathit{pr}\left(\mathit{\theta }|\mathcal{D}\right)$ and $\left\{{\varsigma }^{\left(1\right)},\dots ,{\varsigma }^{\left(S\right)}\right\}$ be the associated realization from each MCMC draw. This draw then enables us to approximate (10) with

$$\int \mathit{pr}\left(J_0|\varsigma ,\mathcal{D}\right)\mathit{pr}\left(\varsigma |\mathcal{D}\right)\mathit{d\varsigma }=\frac{1}{S}\sum _{s=1}^S\mathit{pr}\left(J_0|{\varsigma }^{\left(s\right)},\mathcal{D}\right).$$

Note that

$$\mathit{pr}\left(J_0|{\varsigma }^{\left(s\right)},\mathcal{D}\right)=\frac{1}{1+\frac{\mathit{pr}\left(J_1|{\varsigma }^{\left(s\right)},\mathcal{D}\right)}{\mathit{pr}\left(J_0|{\varsigma }^{\left(s\right)},\mathcal{D}\right)}}$$

and therefore, we have reduced the problem of assessing $\mathit{pr}\left(J_0|\mathcal{D}\right)$ to evaluating several CBFs. At this juncture, note that

$$\mathit{pr}\left(J_0|{\varsigma }^{\left(s\right)},\mathcal{D}\right)\propto \mathit{pr}\left({\varsigma }^{\left(s\right)}|J_0,\mathcal{D}\right)\bullet \mathit{pr}\left(J_0\right)={\int }_0^{\infty }\mathit{pr}\left({\varsigma }^{\left(s\right)}|\tau ,\mathcal{D}\right)\mathit{pr}\left(\tau \right)\mathit{d\tau }\bullet \mathit{pr}\left(J_0\right),$$

while

$$\mathit{pr}\left(J_1|{\varsigma }^{\left(s\right)},\mathcal{D}\right)\propto \mathit{pr}\left({\varsigma }^{\left(s\right)}|J_1,\mathcal{D}\right)\bullet \mathit{pr}\left(J_1\right)={\int }_0^{\infty }{\int }_{{ℝ}^q}\mathit{pr}\left({\varsigma }^{\left(s\right)}|\tau ,\mathit{\xi },\mathcal{D}\right)\mathit{pr}\left(\mathit{\xi },\tau \right)d\mathit{\xi }\mathit{d\tau }\bullet \mathit{pr}\left(J_1\right).$$

Evaluation of these integrals thus requires the specification of priors $\mathit{pr}\left(\tau \right)$ under $J_0$ and $\mathit{pr}\left(\mathit{\xi },\tau \right)$ under $J_1$ . Under the model $J_0$, we propose the standard prior

$$\tau ~\Gamma \left(1/2,1/2\right)$$

which yields

$$\mathit{pr}\left(J_0|{\varsigma }^{\left(s\right)},\mathcal{D}\right)\propto {\left(\frac{1}{2}+\frac{{\mathit{\varsigma }}^{{\left(s\right)}^{\prime }}{\mathit{\varsigma }}^{\left(s\right)}}{2}\right)}^{-\left(n+1\right)/2}.$$E11

For $J_1$, we use the prior

$$\tau ~\Gamma \left(1/2,1/2\right)$$

$$\mathit{\xi }\mid \tau ~\mathcal{N}\left(0,{\tau }^{-1}{\mathbb{I}}_q\right)$$

which yields

$$\mathit{pr}\left(J_1|{\varsigma }^{\left(s\right)},\mathcal{D}\right)\propto {\left|\mathit{\Xi }\right|}^{-\frac{1}{2}}{\left(\frac{1}{2}+\frac{{\left({\varsigma }^{\left(s\right)}-\mathit{Z}{\widehat{\mathit{\xi }}}^{\left(s\right)}\right)}^{\prime }\left({\varsigma }^{\left(s\right)}-\mathit{Z}{\widehat{\mathit{\xi }}}^{\left(s\right)}\right)}{2}\right)}^{-\frac{n+1}{2}},$$E12

where

$$\mathit{\Xi }=\tau \left({\mathit{Z}}^{\prime }\mathit{Z}+{\mathbb{I}}_q\right),$$

$$\widehat{\mathit{\xi }}=\tau {\mathit{\Xi }}^{-1}\mathit{Z}{\mathit{\varsigma }}^{\left(s\right)}.$$

This approach offers similar performance to the Sargan test, which has the desirable feature that it is a fully Bayesian approach, as opposed to the approximate test of [16], and it can be directly embedded in the Gibbs sampling procedures outlined above. We emphasize in the discussion section that further work can be done on this diagnostic.

5. Empirical study: determinants of opening box office

In this section, we consider a generalized linear model with an identity link in the presence of multiple endogenous variables and covariates based on the IVBMA framework incorporating model uncertainty. Based on previous studies of box office revenues, we estimate the effects of three endogenous predictors, prelaunch advertising spending, the number of screens, and production budget with other covariates on opening box office.

Several studies have established a significant link between advertising expenditures and box-office grosses [47, 48, 49, 50]. Almost 90% of a movie’s advertising budget is allocated in the weeks leading up to the theatrical launch [49] shows the importance of prerelease advertising. The number of screens on which a movie is released has been recognized as one of the most significant factors related to the box office [51, 52, 53]. Prerelease advertising spending and the number of opening screens need to be considered endogenous because it is plausible for movies that are expected to generate high box office gross to receive more advertising and distribution. That is, advertising spending and distribution are more likely to be determined by expected box office revenues.

Major studios dominate the movie marketplace regarding film production and distribution. The production budget is an essential predictor because big budgets translate into the casting of top actors and directors, lavish sets and costumes, special effects, and expensive digital manipulations, leading to heightened audience attractiveness [54, 55]. Previous studies [55, 56, 57] used production budget as a direct influencer or moderating variable, but it is also the studio’s strategic decision using knowledge about viewers and competitors’ actions, that is, the data reflect firm’s strategic behavior [58]. While researchers examined endogeneity in advertising responsiveness using a control function approach [14] or price endogeneity using Gaussian copula [9], they did not simultaneously control for multiple endogenous variables or incorporate model uncertainty. The proposed approach can test the effects of three endogenous variables in a generalized framework.

5.1 Description of the data

Starting from all movies released by major studios from 2006 to 2007, we analyzed 130 movies, including 16 animation and 50 R-rated movies, based on the IMDb database. We have excluded films without the complete prerelease advertising information from TNS Media Intelligence. Advertising data include the total dollar value of prerelease media expenditure across 17 different media. The number of opening screens, production budget, and opening box office gross are obtained from IMDb.com and BoxOfficeMojo.com. Table 1 shows the summary statistics of the dependent variable and three endogenous variables. Opening box office gross varies from less than a million to over 100 million dollars. The production budget represents the most significant expense for movie studios [49]. For movies in our sample, they are about $52 million on average and vary from $4 million to $210 million. It becomes crucial for films with high production costs to succeed at the box office to recover their costs, resulting in higher advertising spending and showing at more theaters.

	Mean	Median	SD	Range
Opening box office	20.48	14.32	17.90	(0.72, 102.75)
Prerelease advertising	4.39	4.16	2.17	(0.69, 10.79)
Number of opening screens	2729	2692	707	(825, 4054)
Production budget	52.15	35.0	44.21	(4, 210)

Table 1.

Summary statistics.

The three endogenous predictors were regressed on eleven potential instruments and thirteen additional covariates, summarized in Table 2. Covariates such as genre, MPAA rating, animation, sequels, and release date are publicly available on IMDb and The Numbers. The genre is classified into seven categories (action, comedy, drama, horror, Sci-Fi, mystery/suspense, and romance), and the MPAA rating into two dummy variables (R, PG-13, and others).

Instrumental variables (Z)	Release time	Period indicator based on 10-year box office gross
		(1 = May, June, July, December/0 = other months)
	Expert	Marketability ratings of industry experts
	Direct	Production and distribution by the same company
	Distributor	Production studio dummy variables (SD1–SD6)
		(1 = FOX, 2 = COLUMBIA, 3 = WARNER BROTHERS, 4 = UNIVERSAL, 5 = PARAMOUNT, 6 = DREAMWORKS, 7 = Others)
Covariates (W)	Seasonality	Seasonal index by decomposition model
	Sequels	Dummy variable
	Animation	Animation movies
	Critics review	Movie ratings from Rotten Tomatoes (0–100 points)
	GD1–GD7	Genre dummy variables
		(1 = action/adventure, 2 = comedy, 3 = drama, 4 = horror, 5 = Sci-Fi, 6 = mystery/suspense, 7 = romance, 8 = others)
	RD1–RD2	MPAA rating dummy variables
		(1 = R, 2 = PG13, 3 = Others)

Table 2.

Description of the instruments (Z) and covariates (W).

The MPAA rating is related to the potential size of viewers. Not R-rated movies are open to more moviegoers from the outset, making it necessary to have wider releases and intensive advertising. Critics’ ratings are obtained from the Rotten Tomatoes website, which gives a composite score of 1–100 based on evaluations from movie critics. A monthly seasonality index was obtained by estimating a decomposition model using a time series of monthly box office gross. The seasonal parameter was optimized at 0.56 with the mean absolute percentage error of 10.5%.

For the two endogenous variables, prerelease advertising and opening number of screens, we have used four common instruments of the 11 variables: (a) movie distributors, (b) release time, (c) average marketability ratings by three industry experts in one of the major studios, and (d) whether the same studio did production and distribution. Studios have considerable discretion over the amount and schedule of prelaunch advertising they allocate to each movie [51]. Because advertising elasticities for motion pictures are significantly higher compared to other industries [52], studios’ decisions on prerelease advertising spending and opening screens would have a significant impact on the success at the box office. We have included eight major studios to examine any studio-specific effects on advertising and distribution. Release time is another critical characteristic since movie advertising is seasonal, as heavily supported movies are usually released in peak seasons [51]. Based on the monthly box office gross from 2001 to 2010, we have found a substantial increase in box office gross in May–July and December. A dummy variable is used to indicate those months. For the third endogenous variable, production budget, we exclude release time and expert ratings since they are unavailable at the time of budget decision. Similarly, the seasonal index and critics review were also excluded from the regression of the production budget. Some major studios like 20th Century Fox and Paramount are vertically integrated, having their distribution division. A dummy variable Direct indicates whether both production and distribution divisions finance a movie. For the common instruments on each endogenous regressor, the proposed IVBMA approach has a built-in capability of variable selection using the posterior inclusion probability.

5.2 Results

Table 3 shows the IVBMA posterior estimates of the first stage. The sum of the models’ posterior probabilities containing the variable is called the inclusion probability [16, 23]. In Table 3, column IncProb shows posterior inclusion probabilities in the first stage, which provide a direct interpretation of the efficacy of an instrumentation strategy. Related to prerelease advertising spending, we find a robust movie-type effect for animation, sequels, and PG-13. Animated family films have performed consistently well at the box office, and Pixar and DreamWorks Animation are the most represented studios. Movie sequels build on the original movies’ commercial success and can be considered a brand extension of the experiential product [59]. Given the original movie’s brand power, a sequel usually achieves box office success [60]. The negative coefficient of Sequels results from relatively low advertising costs, which is one of the benefits of brand extensions [61]. The posterior inclusion probabilities of Animation and Sequels are 0.9, which shows generous production budgets for those movies. The marketability ratings by industry experts are significant predictors for prerelease advertising and opening screens. Considering that the ratings are based on the feedback from advance movie screenings, they are reliable indicators of box office performance accompanied by heavy advertising and broader release.

	Prerelease advertising			Opening screens			Production budget
	IncProb	Mean	Quantile	IncProb	Mean	Quantile	IncProb	Mean	Quantile
Intercept	0.436	−0.130	(−1.437, 0.872)	0.542	0.356	(−0.908, 2.189)	1	16.496	(16.037, 16.909)
Sequels	1	−0.415	(−0.480, −0.348)	0.201	0.025	(0, 0.211)	0.909	0.491	(0, 0.914)
Animation	1	0.197	(0.133, 0.266)	0.365	0.061	(0, 0.296)	0.901	0.562	(0, 1.083)
Seasonal index	1	1.418	(1.004, 1.836)	0.817	0.745	(0, 1.727)
Critics review	1	1.111	(1.037, 1.184)	0.070	0	(−0.025, 0.032)
R	0.089	0.004	(0, 0.068)	0.853	−0.139	(−0.288, 0)	0.242	0.041	(−0.123, 0.445)
PG-13	1	0.511	(0.467, 0.564)	0.189	−0.017	(−0.174, 0)	0.964	0.436	(0, 0.791)
Action/adventure	0.026	0	(0, 0)	0.055	0.001	(0, 0.029)	0.987	0.494	(0.183, 0.771)
Comedy	0.048	−0.002	(−0.033, 0)	0.062	0.002	(0, 0.038)	0.539	0.181	(0, 0.631)
Drama	0.048	0.001	(0, 0.030)	0.903	−0.149	(−0.259, 0)	0.168	−0.009	(−0.225, 0.124)
Horror	1	−0.229	(−0.289, −0.169)	0.081	−0.004	(−0.084, 0)	0.446	−0.139	(−0.595, 0)
Sci-Fi	0.062	0.003	(0, 0.051)	0.102	−0.006	(−0.112, 0)	0.202	0.017	(−0.178, 0.331)
Mystery/suspense	0.057	0.002	(0, 0.035)	0.046	0	(0, 0)	0.212	0.036	(−0.004, 0.358)
Romance	0.095	−0.005	(−0.068, 0)	0.125	−0.011	(−0.150, 0)	0.316	−0.085	(−0.600, 0.046)
Direct	0.028	0	(0, 0)	0.093	0.005	(0, 0.087)	0.724	0.222	(0, 0.543)
Fox	0.999	−0.147	(−0.203, −0.090)	0.065	0	(−0.025, 0.015)	0.219	0.021	(−0.196, 0.383)
Columbia	0.028	0	(0, 0)	0.074	−0.002	(−0.057, 0.003)	0.724	0.331	(0, 0.838)
Paramount	0.030	0	(0, 0)	0.114	0.009	(0, 0.134)	0.994	0.725	(0.311, 1.124)
Universal	0.034	0	(0, 0)	0.072	0.003	(0, 0.077)	0.885	0.488	(0, 0.943)
Warner Brothers	0.053	0.002	(0, 0.036)	0.071	−0.001	(−0.041, 0.020)	0.987	0.712	(0.253, 1.144)
MGM	0.038	0.001	(0, 0)	0.089	0.005	(0, 0.118)	0.249	0.053	(−0.088, 0.504)
Lions Gate	1	−0.421	(−0.498, −0.343)	0.404	−0.073	(−0.314, 0)	0.948	−0.676	(−1.161, 0)
Buena Vista	0.040	0	(0, 0)	0.092	0	(−0.051, 0.057)	0.208	0.008	(−0.262, 0.348)
Expert	1	2.134	(1.886, 2.449)	1	1.596	(1.139, 1.947)
Release Time	0.026	0	(0, 0)	0.068	0.002	(0, 0.054)

Table 3.

IVBMA results (first stage).

As expected, a seasonal index shows a high inclusion probability for both endogenous variables, which aligns with the common belief that movies with high expected gross are carefully scheduled to be released in peak seasons. Release time, however, shows no impact, and the result is mainly due to the sample characteristic that more than 65% of the movies in the sample were released in historically no peak months. Note that a seasonal index is calculated for the duration under investigation (2006–2007) while Release time is based on a 10-year window. Therefore, a seasonal index captures short-term fluctuations more accurately.

For prerelease advertising, the PG-13 rating is included with probability one. It concerns the size of potential viewers since non-R ratings imply greater reach among moviegoers, which may result in a higher level of advertising. There is empirical evidence from more than one systematic investigation to show that R-rated movies generate smaller revenues than those with less restrictive ratings [47, 62]. We also find that a dummy variable GD5 for Horror films is a significant predictor of prerelease advertising. This result may reflect the popular trend at that time. There are 15 horror movies in the sample including I am Legend, Silent Hills, and Saw III, which have been very successful at the box office. Consistent with previous literature on critics’ reviews [49, 63], we find a significant impact of reviews on movie advertising. The industry practice of using critics’ quotations in film advertisements supports the continuing authority of film critics. The use of critics’ reviews in movie advertisements indicates distributors’ beliefs and the significance of critics as a cultural intermediary for audiences [64].

In contrast, critical reviews were not included in explaining opening screens. It is consistent with the findings that the relationship between reviews and distributor’s decision is spurious [65], and there is only a positivity bias of exhibitors such that an excellent review allows a movie to stay longer on-screen while negative reviews do not shorten a film’s run [66]. That is, critical reviews do not influence an exhibitor’s decision to keep or withdraw a movie from a theater.

As shown in Table 3, regarding production budget, distributor effects are evident from the high inclusion probabilities of the studios besides movie characteristics such as Sequels and Animation. Though 21st Century Fox and Columbia have released more movies than other studios (37 in the sample), Paramount, Universal, and Warner Brothers had a higher average production budget per movie among major studios, and Lions Gate was the leading independent producer/distributor from 2006 to 2007. PG-13 rating, combined with the Action/Adventure genre consistently performs better than others at the box office by broadening its audience appeal [47]. Interestingly, the instrument, Direct, has a high inclusion probability only for the production budget. It is the case that the deals struck between distributors and exhibitors when they are separately owned are different as vertically integrated studios that are keen to get more movies through their theaters at all times because this maximizes returns from ticket sales and ancillary items such as food and drink. When the audiences start to fall, an exhibitor will prefer to end its run and show another new movie that will boost attendance figures again. Exhibitors favor signing short-run contracts for movies, but signing can be avoided if the same studio controls production, distribution, and exhibitions [67].

Table 4 shows the IVBMA posterior estimates of the second-stage regression. As discussed in section 4, we have tested instrumental validity based on a Bayesian approach. As mentioned in Section 4, the validity score represents the probability that the instrument condition is not satisfied. All instruments used in the study are essentially zero, which strongly supports the validity of the instrumentation choices. In the second stage, several variables are essential predictors of opening box office revenues. As expected, the number of opening screens and prerelease advertising are significant determinants of opening box office gross with high inclusion probabilities. Though it is difficult to disentangle the causal effect of advertising on sales using data on actual box office receipts, it is consistent with previous findings that prerelease advertising has a positive and statistically significant impact on public awareness of a movie and its box office performance [47, 49, 50, 68]. While Elberse and Eliashberg [52] argue that movie attributes and advertising expenditures mostly influence revenues indirectly through their impact on exhibitors’ screen allocations, this result supports a significant direct effect of advertising. The number of opening screens is the most important predictor, with an inclusion probability of one, which is also consistent with previous findings [53, 69, 70]. It seems to be the case that the more screens on which new movies were released, the bigger their initial audiences. The higher the audience for a movie in the opening weekend, the higher would be its audience the following week. While audiences inevitably drop off over time, a movie’s cinema run would be longer if it got off to a good start. Considering a typically high correlation between opening screens and prerelease advertising, studios’ advertising and distribution approaches may be very similar. Other than these two factors, Sequels and Drama show high inclusion probabilities, which may only reflect the characteristics of successful movies in the sample. Though we initially expected a significant effect of seasonality, it turns out to have a weak influence, though it remains relevant. Production budget has low inclusion probability, and it suggests that a movie’s production cost is an indicator of the creative talent involved or the extent to which the movie incorporates expensive special effects or uses elaborate set designs [49], but not a good indicator of success. For about 90 films released in the United States from 2008 to 2012 with budgets of more than $100 million, most of them failed to generate enough revenues at the box office to cover their costs [71]. After all, big budgets do not guarantee success, and the only way to know how audiences react to a movie is to wait until it has been released and moviegoers have had the opportunity to see it.

	IncProb	Mean	SD	Quantile	Conditional
					Mean	SD
Constant	0.529	−0.284	0.725	(−2.078, 0.989)	−0.525	0.926
Sequels	0.963	0.564	0.203	(0, 0.915)	0.585	0.174
Animation	0.147	−0.001	0.065	(−0.166, 0.163)	−0.001	0.170
R	0.103	0.002	0.041	(−0.071, 0.096)	0.016	0.126
PG-13	0.201	0.034	0.096	(−0.001, 0.345)	0.172	0.149
Action/adventure	0.098	0.002	0.034	(−0.041, 0.097)	0.031	0.104
Comedy	0.142	−0.015	0.056	(−0.205, 0.004)	−0.102	0.116
Drama	0.797	−0.238	0.157	(−0.509, 0)	−0.298	0.113
Horror	0.264	0.051	0.113	(0, 0.383)	0.193	0.143
Sci-Fi	0.150	0.007	0.070	(−0.117, 0.215)	0.054	0.173
Mystery/suspense	0.136	−0.012	0.050	(−0.183, 0.007)	−0.089	0.108
Romance	0.129	0.003	0.052	(−0.104, 0.138)	0.019	0.143
Seasonal index	0.569	−0.446	0.668	(−2.013, 0.445)	−0.781	0.723
Critics review	0.291	0.038	0.216	(−0.371, 0.670)	0.130	0.384
Prerelease advertising	0.918	0.452	0.200	(0, 0.758)	0.492	0.153
Opening screens	1	1.287	0.306	(0.756, 1.963)	1.285	0.306
Production budget	0.120	0.008	0.046	(−0.019, 0.164)	0.071	0.115

Table 4.

IVBMA results (second stage).

6. Conclusion

Market response models often use endogenous regressors since marketing activities are nonrandom and reflect the firm’s strategic behavior. Thus, ignoring the endogeneity of marketing actions will lead to incorrect estimates of response parameters and, consequently, to biased inferences [4, 58]. While researchers have developed various approaches to dealing with endogeneity, including the control function approach, Gaussian copula, or instrument-free approaches, the IV approach remains the technique of choice when dealing with endogeneity in econometrics and other areas of applied research. Almost invariably, empirical work in economics and marketing will be subject to much uncertainty about model specifications. This may be the consequence of the existence of different theories or different ways in which theories can be implemented in empirical models or other aspects such as assumptions about heterogeneity or independence of the observables [72]. It is important to realize that this uncertainty is an inherent part of the marketing response modeling.

We have proposed a computationally efficient solution to the problem of incorporating model uncertainty into IV estimation. The IVBMA method leverages an existing Gibbs sampler and shows that by nesting model moves inside this framework, model averaging can be performed with minimal additional effort. In contrast to the approximate solution proposed by [16], our method yields a theoretically justified, fully Bayesian procedure. The applied examples show this method’s benefit, by enabling additional factors to be entertained by the researcher, which are either incorporated where appropriate or promptly dropped.

The CBF approach is only one manner of incorporating model uncertainty in the framework considered. Two other options would be reversible jump schemes [29, 30] or specify a spike and slab prior [73]. We have chosen our approach because it fits nicely into the Gibbs sampling framework, unlike the reversible jump procedure of Koop et al. [29], and still explicitly incorporates uncertainty at the model level, unlike spike and slab type priors at the variable level. However, additional research is needed to explore the tradeoffs between these alternative methods of incorporating model uncertainty.

One assumption crucial to the Gibbs sampler’s functioning is the multivariate normality of the residuals in Eq. (2). Conley et al. [74] discuss a Bayesian approach that allows nonparametric estimation of the distribution of error terms in a set of simultaneous equations using a Dirichlet process mixture (DPM). We note that the IVBMA methodology can readily incorporate the DPM framework by simply replacing the IV kernel distributions of [36] with IVBMA kernel distributions. A nonparametric IVBMA approach based on non-normal errors will be one of the model extensions in the future. Another critical issue is assessing instruments’ validity in implementing IV methods. The Bayesian version of the Sargan test that we have proposed serves as a natural starting point for more involved methodologies, including latent factors though many features still need to be investigated on this front compared to other strategies.

IVBMA has the potential to be extended to more complicated likelihood frameworks. The proposed model can be extended to latent constructs in the context of structural equations modeling with latent Gaussian factors and, at the same time, selecting the suitable path model [75]. Survival analysis is another area that can benefit from the IVBMA approach in dealing with multiple endogenous regressors and implementing more flexible hazard specifications beyond the proportional hazard model [76]. Since the entire method uses a Gibbs framework, it can be incorporated in any setting where endogeneity, model uncertainty, and latent normality are present. In particular, the linear specification can be relaxed using semiparametric methods such as splines or more flexible approaches involving Gaussian processes. While the algorithms involved would understandably become more complex, the central concept involving using CBFs to assess model uncertainty would remain pertinent.

Here we outline the calculation of $\mathit{pr}\left(\mathcal{D}|M_r,{\mathit{\beta }}_{-r},\mathit{K}\right)$ . Note that

$$\mathit{pr}\left(\mathcal{D}|M_r,{\mathit{\beta }}_{-r},\mathit{K}\right)={\int }_{{\Lambda }_{M_r}}\mathit{pr}\left(\mathcal{D}|{\mathit{\beta }}_r,{\mathit{\beta }}_{-r},\mathit{K}\right)\mathit{pr}\left({\mathit{\beta }}_r|M_r\right)d{\mathit{\beta }}_r.$$

Let ${\mathit{U}}_{M_r}^{\left(r\right)}$ be the submatrix of ${\mathit{U}}^{\left(r\right)}$ associated with the variables in $M_r$ and set ${\tilde{\mathit{Y}}}_r$ as above. Then

$${\int }_{{\Lambda }_{M_r}}\mathit{pr}\left(\mathcal{D}|{\mathit{\beta }}_r,{\mathit{\beta }}_{-r},\mathit{K}\right)\mathit{pr}\left({\mathit{\beta }}_r|M_r\right)d{\mathit{\beta }}_r\propto {\int }_{{\Lambda }_{M_r}}{\left(2\pi \right)}^{-\frac{\left|M_r\right|}{2}}\mathit{exp}\left(-\frac{1}{2}\left[-2{\widehat{\mathit{\beta }}}_{M_r}{\mathit{\Omega }}_{M_r}{\mathit{\beta }}_r+{\mathit{\beta }}_{r^{\prime }}{\mathit{\Omega }}_{M_r}{\mathit{\beta }}_r\right]\right)d{\mathit{\beta }}_r,$$

where

$${\Omega }_{M_r}=K_{\mathit{rr}}{\mathit{U}}_{M_r}^{{\left(r\right)}^{\prime }}{\mathit{U}}_{M_r}^{\left(r\right)}+{\mathbb{I}}_{\mid M_r\mid },$$

$${\widehat{\mathit{\beta }}}_{M_r}=K_{\mathit{rr}}{\mathit{\Omega }}_{M_r}^{-1}{\mathit{U}}_{M_r}^{{\left(r\right)}^{\prime }}{\tilde{\mathit{Y}}}_r.$$

We can now see that the term in the integral is the canonical form of a Gaussian distribution. Appropriate completion therefore yields

$$\mathit{pr}\left(\mathcal{D}|M_r,{\mathit{\beta }}_{-r},\mathit{K}\right)\propto {\left|{\Omega }_{M_r}\right|}^{-1/2}\mathit{exp}\left(-\frac{1}{2}{{\widehat{\mathit{\beta }}}_{M_r}}^{\prime }{\mathit{\Omega }}_{M_r}{\widehat{\mathit{\beta }}}_{M_r}\right)$$

Let

$$Y_{i1}~\mathcal{P}\left({\mathit{U}}_i^{{\left(r\right)}^{\prime }}{\mathit{\beta }}_i+{\epsilon }_{i1}\right),$$

and for $r>1$,

$$Y_{\mathit{ir}}={\mathit{U}}_i^{{\left(r\right)}^{\prime }}{\mathit{\beta }}_r+{\epsilon }_{\mathit{ir}},$$

where

$${\mathit{\epsilon }}_i~\mathcal{N}\left(0,{\mathit{K}}^{-1}\right).$$

The MCMC for this model roughly follows the algorithm mentioned above, but with the additional handling of the random effect ${\epsilon }_{i1}$ and the subsequent updating of ${\mathit{\beta }}_1$. Note that

$$\mathit{pr}\left({\epsilon }_{i1}|·\right)\propto \mathit{pr}\left(Y_i|{\mathit{U}}_i^{\left(1\right)},{\mathit{\beta }}_1,{\epsilon }_{i1}\right)\mathit{pr}\left({\epsilon }_{i1}|{\mathit{\epsilon }}_i\backslash {\epsilon }_{i1},\mathit{K}\right)$$

where

$$\mathit{pr}\left({\epsilon }_{i1}|{\mathit{\epsilon }}_i\backslash {\epsilon }_{i1},\mathit{K}\right)=\mathcal{N}\left({\eta }_i,{\kappa }_i^{-1}\right)$$

with

$${\eta }_i=-\sum _{r=2}^R\frac{K_{1r}}{K_{11}}{\epsilon }_{\mathit{ir}}$$

$${\kappa }_i=\frac{1}{K_{11}}$$

Further, denote ${\mu }_i=\underset{i}{\overset{{\left(1\right)}^{\prime }}{\mathit{U}}}{\mathit{\beta }}_1$. Then

$$\mathit{pr}\left({\epsilon }_{i1}|·\right)\propto \mathit{exp}\left(-\mathit{exp}\left({\mu }_i+{\epsilon }_{i1}\right)+\left({\mu }_i+{\epsilon }_{i1}\right)Y_{i1}\right)\mathit{exp}\left(-\frac{1}{2}{\kappa }_i{\left({\epsilon }_{i1}-{\eta }_i\right)}^2\right).$$

Writing

$$f\left({\epsilon }_{i1}\right)=-\mathit{exp}\left({\mu }_i+{\epsilon }_{i1}\right)+\left({\mu }_i+{\epsilon }_{i1}\right)Y_{i1}-\frac{1}{2}{\kappa }_i{\left({\epsilon }_{i1}-{\eta }_i\right)}^2$$

we have

$$f^{\prime }\left({\epsilon }_{i1}\right)=-\mathit{exp}\left({\mu }_i+{\epsilon }_{i1}\right)+Y_{i1}-{\kappa }_i\left({\epsilon }_{i1}-{\eta }_i\right)$$

$$f^{″}\left({\epsilon }_{i1}\right)=-\mathit{exp}\left({\mu }_i+{\epsilon }_{i1}\right)-{\kappa }_i$$

Hence, by setting

$$b\left({\epsilon }_{i1}\right)=f^{\prime }\left({\epsilon }_{i1}\right)-f^{″}\left({\epsilon }_{i1}\right){\epsilon }_{i1}$$

$$c\left({\epsilon }_{i1}\right)=-f^{″}\left({\epsilon }_{i1}\right)$$

we may sample ${{\epsilon }_{i1}}^{\prime }~\mathcal{N}\left(b\left({\epsilon }_{i1}\right)/c\left({\epsilon }_{i1}\right),1/c\left({\epsilon }_{i1}\right)\right)$ and accept this update with probability $\mathit{min}\left\{\alpha ,1\right\}$ where

$$\alpha =\frac{\mathit{pr}\left(Y_{i1}|{\mu }_i,{{\epsilon }_{i1}}^{\prime }\right)\mathit{pr}\left({{\epsilon }_{i1}}^{\prime }|{\eta }_i,{\kappa }_i\right)\mathit{pr}\left({\epsilon }_{i1}|b\left(\underset{i1}{\overset{\prime }{\epsilon }}\right),c\left(\underset{i1}{\overset{\prime }{\epsilon }}\right)\right)}{\mathit{pr}\left(Y_{i1}|{\mu }_i,{\epsilon }_{i1}\right)\mathit{pr}\left({\epsilon }_{i1}|{\eta }_i,{\kappa }_i\right)\mathit{pr}\left(\underset{i1}{\overset{\prime }{\epsilon }}|b\left({\epsilon }_{i1}\right),c\left({\epsilon }_{i1}\right)\right)}.$$

Once all ${\epsilon }_{i1}$ are updated, other updates mostly follow the steps above.