Bayesian Variable Selection for Mixture Process Variable Design Experiment

Sadiah M. A. Aljeddani

Department of Mathematics, Al-Lith College, Umm Al-Qura University, Al-Lith, Saudi Arabia.

DOI: 10.4236/ojmsi.2022.104022                

Abstract

This paper discussed Bayesian variable selection methods for models from split-plot mixture designs using samples from Metropolis-Hastings within the Gibbs sampling algorithm. Bayesian variable selection is easy to implement due to the improvement in computing via MCMC sampling. We described the Bayesian methodology by introducing the Bayesian framework, and explaining Markov Chain Monte Carlo (MCMC) sampling. The Metropolis-Hastings within Gibbs sampling was used to draw dependent samples from the full conditional distributions which were explained. In mixture experiments with process variables, the response depends not only on the proportions of the mixture components but also on the effects of the process variables. In many such mixture-process variable experiments, constraints such as time or cost prohibit the selection of treatments completely at random. In these situations, restrictions on the randomisation force the level combinations of one group of factors to be fixed and the combinations of the other group of factors are run. Then a new level of the first-factor group is set and combinations of the other factors are run. We discussed the computational algorithm for the Stochastic Search Variable Selection (SSVS) in linear mixed models. We extended the computational algorithm of SSVS to fit models from split-plot mixture design by introducing the algorithm of the Stochastic Search Variable Selection for Split-plot Design (SSVS-SPD). The motivation of this extension is that we have two different levels of the experimental units, one for the whole plots and the other for subplots in the split-plot mixture design.

Keywords

Variable Selection, Bayesian Analysis, Mixture Experiment, Split-Plot Design

1. Introduction

Bayesian methods are important approaches due to their ability to quantify uncertainty. In such an approach, prior distributions that represent subjective beliefs about parameters are assigned to the regression coefficients. By applying Bayes’ rule, prior beliefs are updated by the data and transformed into posterior distributions, on which all inference is based.

An important task in regression building is to determine which variables should be included in the model. Therefore, the principle of Bayesian variable selection is to get the large and the small effects distinctive, and effective prior essences mass around zero and distribute the remaining over the parameter space. Such a prior represents the fact that there are small coefficients close to zero on one hand and larger coefficients on the other hand. These priors can be built as a combination of two distributions, along with a narrow normal continuous distribution centred at zero with a small variance called a “spike”, and the other with a flat normal continuous distribution with a large variance to spread over a wider range of parameter values called a “slab”. This type of priors is called “Spike-and-Slab” priors (see Figure 1). Practically saying, those priors are beneficial for purposes of variable selection because they permit the classification of the regression coefficients into two groups: one group with large, significant effects, and the other group with small, negligible effects.

As reviewed by [1] [2] introduced Bayesian variable selection via “spike-and-slab” prior distributions. The spike prior that they used was a probability mass at zero to remove the non-significant variables. Their slab is the uniform distribution with a large symmetric range in order to keep the significant variables. Following their work, many priors were proposed to implement the spike-and-slab property. [3] proposed the Stochastic Search Variable Selection (SSVS) in which the coefficients are sampled from a mixture of two normal distributions with different variances. The spike part is the distribution with a small variance while the slab

Figure 1. Gaussian mixture prior for $\beta$.

part is the distribution with a much larger variance. Also, [4] proposed positive mass at zero for the spike part and a normal distribution for the slab part. In addition, [5] proposed a Bayesian approach for model selection in fractionated split-plot experiments with application to “robust-parameter design”. In their work, they extend the SSVS algorithm of [6] to account for the split-plot error structure. They derive an expression for the posterior probability of a model that requires the computation of, at most, two unidimensional integrals, and employ this quantity for model selection. They were able to integrate the coefficients and the variance components from the joint posterior distribution of all parameters because they use the conjugate normal-inverse gamma prior for these parameters. The integrals are computed with Gaussian quadrature, and Global and Local search algorithms to find models with high posterior probabilities.

The above review summarised the former research and analyzed the problems in the Bayesian variable selection field. On the other side, we notice that the frequentist methods analysis of this kind of experimental design would fail in high Type I error rate [7]. The difficulty of using the frequentist analysis to estimate the variance components is explained in [8]. To overcome this issue, we apply the Bayesian variable selection as this method introduces prior distribution about the variance component in the linear mixed models.

The novel contribution of this paper to Bayesian variable selection is motivated by a very specific experimental design of data from experiments subject to restricted randomisation. We have two different levels of the experimental units one for the whole plots and the other for the subplots in the split-plot design; see, for example, [9]. To address this issue, we adapt the SSVS algorithm in which we sample the subplot coefficients using a mixture of normal posterior distributions with a slab variance different from the slab variance which will be used in the mixture of normal posterior distributions for the whole-plot coefficients. This method reduces Type I and II error rates as well as reduces the prediction error for split-plot design rather than applying the SSVS algorithm in which all coefficients will be sampled from a mixture of normal posterior distributions with one slab variance. We called the approached method the Stochastic Search Variable Selection for Split-Plot Design (SSVS-SPD). The frequentist analysis is dependent on the estimates of the variance components, yet these estimates cannot be precisely calculated because of the deficiency of the degrees of freedom for the random effects in the split-plot design. This issue was discussed by [8]. Introducing a prior distribution for the variance components in the linear mixed model provides additional information to overcome the problem of variance estimation.

This paper differs from [10] in which the dataset used in this paper is the Vinyl Thickness experiment, which had not been used in variable selection. Also, it differs in the choice of the shape parameter of the correlation parameter as in this paper we used $a=b=2$ to fit this data. A preprint has previously been published [11].

2. Split-Plot Design and Sample Model

The model for the block experiments includes two types of errors: block error and residual error. Hence, linear mixed models (LMMs) are used to analyse responses from the blocked experiments.

Linear mixed-effects models (LMMs) introduce correlations between observations using random effects. This leads to the use of generalised least squares (GLS) estimation, combined with restricted maximum likelihood estimation (REML) of the variance components as will be discussed. This type of analysis is used by the most design of experiments textbooks that deal with blocked designs. In matrix notation, the model corresponding to a blocked design is written as

$Y=X\beta +Z\gamma +\epsilon ,$ (1)

where Y is $n\times 1$ vector of observations on the response of interest, X is the $n\times p$ model design matrix containing the polynomial expansions of the m factor levels at the n experimental runs, $\beta$ is the $p\times 1$ vector of unknown fixed parameters, Z is an $n\times b$ random design matrix which represents the allocation of the runs to blocks, and whose $\left(i\mathrm{,}j\right)$ ^th element is one where the i^th observation belongs to the j^th blocks, and zero otherwise. If the runs of the experiment are grouped per block, then Z is of the form

$Z=\text{diag}\left[1_{k_1}\mathrm{,}{1}_{k_2}\mathrm{,}\cdots \mathrm{,}{1}_{k_b}\right]\mathrm{,}$ (2)

where $1_k$ is a k vector of ones, and $k_1\mathrm{,}{k}_2\mathrm{,}\cdots \mathrm{,}{k}_b$ are the blocks sizes. The random effects of the b blocks are contained within the $b\times 1$ vector $\gamma$, and the random errors are contained within the $n\times 1$ vector $\epsilon$. It is assumed that $\gamma$ and $\epsilon$ are independent and normally distributed, i.e. $\text{Cov}\left(\gamma ,\epsilon \right)=0_{b\times n}$, where $0_{b\times n}$ is the $b\times n$ matrix of zeros. Hence, $\gamma ~\text{N}\left(0_b\mathrm{,}{\sigma }_{\gamma }^2{{\rm I}}_b\right)$, and $\epsilon ~\text{N}\left(0_n\mathrm{,}{\sigma }_{\epsilon }^2{I}_n\right)$, where $0_b$ and $0_n$ are the b and n column vectors of zeros respectively, and $I_b$ and $I_n$ are the b-dimensional and n-dimensional identity matrices respectively.

Under these assumptions, Y is a normally distributed random variable with mean $\mathbb{E}\left(Y\right)=X\beta$, and the variance-covariance matrix of the response Y can be written as

$V=\text{Var}\left(Y\right)=\text{Var}\left(X\beta +Z\gamma +\epsilon \right)$ (3)

$=\text{Var}\left(Z\gamma \right)+\text{Var}\left(\epsilon \right)$ (4)

$=Z\text{Var}\left(\gamma \right)Z^{\prime }+{\sigma }_{\epsilon }^2{I}_n$ (5)

$={\sigma }_{\gamma }^{2}Z{Z}^{\prime }+{\sigma }_{\epsilon }^2{I}_n\mathrm{.}$ (6)

V can be given as a block diagonal,

$V=\left[\begin{array}{cccc}{V}_1& 0& \cdots & 0\\ 0& V_2& \cdots & 0\\ ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & 0& V_b\end{array}\right],$

where

$V_i={\sigma }_{\epsilon }^2{I}_{k_i}+{\sigma }_{\gamma }^2{1}_{k_i}{1^{\prime }}_{k_i},$

and

$V_i=\left[\begin{array}{cccc}{\sigma }_{\epsilon }^2+{\sigma }_{\gamma }^2& {\sigma }_{\gamma }^2& \cdots & {\sigma }_{\gamma }^2\\ {\sigma }_{\gamma }^2& {\sigma }_{\epsilon }^2+{\sigma }_{\gamma }^2& \cdots & {\sigma }_{\gamma }^2\\ ⋮& \ddots & \ddots & ⋮\\ {\sigma }_{\gamma }^2& \cdots & {\sigma }_{\gamma }^2& {\sigma }_{\epsilon }^2+{\sigma }_{\gamma }^2\end{array}\right].$

As a result, the variance-covariance matrix $V_i$ of all observations within one block is compound symmetric: the main diagonal of the matrix contains the variances of the observations, while the off-diagonal elements are covariances. However, $V_i$ can be rewritten as

$V_i={\sigma }_{\epsilon }^2\left(I_{k_i\times k_i}+\frac{{\sigma }_{\gamma }^2}{{\sigma }_{\epsilon }^2}{1}_{k_i}{1^{\prime }}_{k_i}\right)$ (7)

$={\sigma }_{\epsilon }^2\left(I_n+\eta Z{Z}^{\prime }\right),$ (8)

where $\eta ={\sigma }_{\gamma }^2/{\sigma }_{\epsilon }^2$ is a measure for the extent to which observations within the same block are correlated. The larger this variance ratio, the stronger observations within the same block are correlated.

When the random error terms as well as the group effects are normally distributed, the maximum likelihood estimate of the unknown model parameter $\beta$ in Equation (1) is the generalised least squares (GLS) estimate. Detecting the estimator $\hat{\beta }$ of $\beta$, requires to minimise

${\left(y-X\beta \right)}^{\prime }{V}^{-1}\left(y-X\beta \right)=y^{\prime }{V}^{-1}y-2{\beta }^{\prime }{X}^{\prime }{V}^{-1}y+{\beta }^{\prime }{X}^{\prime }{V}^{-1}X\beta$ (9)

with respect to $\beta$, which is tantamount to detecting $\hat{\beta }$, so that

$\left(X^{\prime }{V}^{-1}X\right)\hat{\beta }=X^{\prime }{V}^{-1}y.$ (10)

Therefore, the generalised least squares (GLS) estimator of $\beta$ is

$\hat{\beta }={\left(X^{\prime }{V}^{-1}X\right)}^{-1}{X}^{\prime }{V}^{-1}Y,$ (11)

and the variance-covariance matrix of the estimators is given by

$\begin{array}{c}\text{Var}\left(\hat{\beta }\right)=\text{Var}\left({\left(X^{\prime }{V}^{-1}X\right)}^{-1}\left(X^{\prime }{V}^{-1}Y\right)\right)\\ ={\left(X^{\prime }{V}^{-1}X\right)}^{-1}{X}^{\prime }{V}^{-1}\text{Var}\left(Y\right){\left({\left(X^{\prime }{V}^{-1}X\right)}^{-1}{X}^{\prime }{V}^{-1}\right)}^{\prime }\\ ={\left(X^{\prime }{V}^{-1}X\right)}^{-1}{X}^{\prime }{V}^{-1}V{V}^{-1}X{\left(X^{\prime }{V}^{-1}X\right)}^{-1}\\ ={\left(X^{\prime }{V}^{-1}X\right)}^{-1}\left(X^{\prime }{V}^{-1}X\right){\left(X^{\prime }{V}^{-1}X\right)}^{-1}\\ ={\left(X^{\prime }{V}^{-1}X\right)}^{-1}.\end{array}$ (12)

Often, the variances ${\sigma }_{\gamma }^2$ and ${\sigma }_{\epsilon }^2$ are not known and therefore, Equation (11) and Equation (12) cannot be used directly. Instead, the estimates of the variance components, ${\hat{\sigma }}_{\gamma }^2$ and ${\hat{\sigma }}_{\epsilon }^2$, are substituted in the GLS estimator as in Equation (11), yielding

$\hat{\beta }={\left(X^{\prime }{V}^{-1}X\right)}^{-1}{X}^{\prime }{V}^{-1}Y\mathrm{,}$ (13)

where

$V={\hat{\sigma }}_{\epsilon }^2{I}_n+{\hat{\sigma }}_{\gamma }^{2}Z{Z}^{\prime }\mathrm{.}$ (14)

In that case, the variance-covariance matrix in Equation (12) can be approximated by

$\text{Var}\left(\hat{\beta }\right)={\left(X^{\prime }{V}^{-1}X\right)}^{-1}\mathrm{.}$ (15)

The generalised least square (GLS) estimator is unbiased, meaning that $\mathbb{E}\left(\hat{\beta }\right)=\beta$, and is equal to the maximum likelihood estimator (MLE). The likelihood function defined as it is the joint probability density function for the observed data examined as a function of the parameters. Hence, the likelihood function for Y in Equation (1) is

$L\left(\beta \mathrm{<mo>\; |\; </mo>}Y\right)={\left(2\pi \right)}^{-n/2}{\left|V\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\left(Y-X\beta \right)}^{\prime }{V}^{-1}\left(Y-X\beta \right)\right]\mathrm{,}$ (16)

where $\pi$ is a constant which does not depend on $\beta$. The maximum likelihood estimator (MLE) is the estimator that maximises the likelihood function, which is tantamount to detecting the $\hat{\beta }$ as

$\frac{\partial }{\partial \beta }L\left(\hat{\beta }\mathrm{<mo>\; |\; </mo>}Y\right)=\mathrm{0,}$ (17)

which is equal to

$\frac{\partial }{\partial \beta }\ln L\left(\hat{\beta }\mathrm{<mo>\; |\; </mo>}Y\right)=\mathrm{0,}$ (18)

where $\ln L\left(\hat{\beta }\mathrm{<mo>\; |\; </mo>}Y\right)$ is the log likelihood function. As Equation (9) is proportionate to log of Equation (16), the GLS estimator in Equation (11) is the result of Equation (17) and Equation (18).

The restricted maximum likelihood (REML) used to estimate ${\sigma }_{\epsilon }^2$ and ${\sigma }_{\gamma }^2$ is

$l_{\text{REML}}\left({\sigma }_{\epsilon }^2\mathrm{,}{\sigma }_{\gamma }^2\mathrm{<mo>\; ;\; </mo>}Y\right)=-\frac{1}{2}\text{ln}\left|V\right|-\frac{1}{2}\text{ln}\left|X^{\prime }{V}^{-1}X\right|-\frac{1}{2}{\left(Y-X\hat{\beta }\right)}^{\prime }{V}^{-1}\left(Y-X\hat{\beta }\right)\mathrm{,}$ (19)

where $\hat{\beta }$ is defined in Equation (13). The restricted log-likelihood $l_{\text{REML}}\left({\sigma }_{\epsilon }^2\mathrm{,}{\sigma }_{\gamma }^2\mathrm{<mo>\; ;\; </mo>}Y\right)$ is minimised with respect to the variance components ${\sigma }_{\epsilon }^2$ and ${\sigma }_{\gamma }^2$ to obtain an unbiased estimate for the variance components.

3. The Bayesian Methodology

3.1. Bayesian Framework

The essential philosophy behind Bayesian inference is to update a prior distribution for an unidentified parameter to a posterior distribution by Bayes’ theorem. Bayes’ theorem can be used to estimate the conditional distributions. While the frequentist approach treats the parameters as unknown and fixed, the Bayesian approach regards them as random variables. We can define the prior distribution $p\left(\theta \right)$ as the probability density (or mass) function which reflects our beliefs about $\theta$ in the parameter space $\Theta$. For given data $y={\left(y_1\mathrm{,}{y}_2\mathrm{,}\cdots \mathrm{,}{y}_n\right)}^{\prime }$, the likelihood function $f\left(y\mathrm{<mo>\; |\; </mo>}\theta \right)$ can then be defined given the parameter $\theta$ for the data y. Also, we can define the posterior density (or mass) function $p\left(\theta \mathrm{<mo>\; |\; </mo>}{y}_1\mathrm{,}{y}_2\mathrm{,}\cdots \mathrm{,}{y}_n\right)$ which represents our updated belief about $\theta$ given the observed data y.

Using Bayes’ theorem, the posterior density of $\theta$ given y is:

$p\left(\theta \mathrm{<mo>\; |\; </mo>}y\right)=\frac{f\left(y\mathrm{<mo>\; |\; </mo>}\theta \right)p\left(\theta \right)}{{\displaystyle {\int }_{\Theta }f\left(y\mathrm{<mo>\; |\; </mo>}\theta \right)p\left(\theta \right)\text{d}\theta }}\mathrm{.}$ (20)

Bayesian inference continues from this distribution. The denominator of Equation (20) is the marginal likelihood of y, and it often does not need to be calculated because it is independent of $\theta$. Bayes’ rule can then be written as:

$p\left(\theta \mathrm{<mo>\; |\; </mo>}y\right)\propto f\left(y\mathrm{<mo>\; |\; </mo>}\theta \right)p\left(\theta \right)\mathrm{.}$ (21)

Equation (21) defines the unnormalised posterior density. The posterior then is proportional to the likelihood × the prior. For more details on Bayesian inference, see [12] and [13].

A prior distribution can be selected based on past information or experimental practice. It can be informative or uninformative. The informative distribution is given numerical information to estimate the parameter of concern. The uninformative reflects equilibrium among outcomes when weak information about the parameter is presented. There are two types of uninformative priors: proper prior and improper prior. The density for proper prior distribution integrates to 1 whereas the integral of the density for an improper distribution is not finite. If the prior integrates to any positive finite value, it is called an unnormalised density and can be renormalised-multiplied by a constant to integrate to 1 [12] [13].

3.2. Markov Chain Monte Carlo (MCMC) Methods

Markov Chain Monte Carlo simulation is a general method based on drawing values of the $\theta$ from approximate distributions, and then correcting those draws to better approximate the target posterior distribution $p\left(\theta \mathrm{<mo>\; |\; </mo>}y\right)$ [12] [13] [14]. A Markov chain can be defined as a sequence of random variables ${\theta }^1\mathrm{,}{\theta }^2\mathrm{,}\cdots$ for which for any iteration t, the distribution of ${\theta }^t$ depends only on the most recent value ${\theta }^{t-1}$ [12] [13] [14]. A Markov chain is generated by sampling ${\theta }^t~p\left(\theta \mathrm{<mo>\; |\; </mo>}{\theta }^{t-1}\right)$. Thisp is called the transition kernel of the Markov chain. Therefore, ${\theta }^t$ depends only on ${\theta }^{t-1}$, not on ${\theta }^0\mathrm{,}{\theta }^1\mathrm{,}\cdots \mathrm{,}{\theta }^{t-2}$.

As $t\to \infty$, the sampling from Markov chain converges to the posterior for the right choice of transition kernel $p\left(\theta \mathrm{<mo>\; |\; </mo>}y\right)$. Thus, we should run the simulation long enough so that the distribution of the current draws is close enough to $p\left(\theta \mathrm{<mo>\; |\; </mo>}y\right)$.

3.2.1. Metropolis-Hastings Sampling

Metropolis-Hastings sampling was proposed by [15] and [16]. The theory of this sampling is based on rejection sampling. The acceptance-rejection method is a technique of getting samples from a distribution with an unknown form. The Metropolis-Hastings algorithm is a common expression for a family of Markov chain simulation methods. It is worth describing the Metropolis algorithm first, then broadening it to discuss the Metropolis-Hastings algorithm. Let $p\left(\theta \mathrm{<mo>\; |\; </mo>}y\right)$ be the conditional posterior distribution where we want to sample from. Let ${\theta }^t$ be the current parameter value, and let $\pi \left(\theta \right)$ be the proposal density. The proposal density is much like a conventional transition operator for a Markov chain, the proposal distribution depends only on the previous state in the chain. However, the transition operator for the Metropolis algorithm has a additional step that assesses whether or not the target distribution has a sufficiently large density near the proposed state to warrant accepting the chain.

3.2.2. Gibbs Sampling

A Gibbs sampler is the simplest of the Markov chain simulation algorithms, and it is used to sample from the conditional conjugate models, where we can directly sample from each conditional posterior [12] [13]. It is rare to find all the conditional posteriors in a model in known forms. One may find some conditional posterior distributions that are possible to directly sample from. Furthermore, one may find some of the conditional posteriors that cannot be straightforwardly sampled from. Therefore, the procedure for this issue is to update the parameters one at a time with the Gibbs sampler used where possible, and one-dimensional Metropolis updating where necessary. This process is called the Metropolis-Hastings within Gibbs sampling and will be used in this work.

4. A Hierarchical Mixture Model for Variable Selection

The linear mixed model fitted to data from a split-plot experiment with n responses is

$y~N\left({\beta }_0{1}_n+X\beta \mathrm{,}V\right)\mathrm{,}$ (22)

where y is $n\times 1$ vector of random responses, ${\beta }_0$ is the intercept, $1_n$ is a $n\times 1$ vector of ones, X is the $n\times p$ model matrix without the column of the intercept, $\beta$ is the $p\times 1$ vector of fixed effect parameters and V is

$V={\sigma }_{\epsilon }^2\left(I_n+\frac{{\sigma }_{\gamma }^2}{{\sigma }_{\epsilon }^2}Z{Z}^{\prime }\right),$

where Z is the random effect design matrix. As $\rho =\frac{{\sigma }_{\gamma }^2}{{\sigma }^2}$, and ${\sigma }^2={\sigma }_{\epsilon }^2+{\sigma }_{\gamma }^2$, then V can be written as

$V={\sigma }^2\left(1-\rho \right)\left(I_n+\frac{\rho }{1-\rho }Z{Z}^{\prime }\right)\mathrm{.}$ (23)

We need to find the highest posterior probability of an indicator vector $\nu =\left({\nu }_1\mathrm{,}{\nu }_2\mathrm{,}\cdots ,{\nu }_p\right)$ such that

${\nu }_j\mathrm{<mo>\; =\; </mo>}(\begin{array}{l}0\text{if}{\beta }_j=0\\ 1\text{if}{\beta }_j\ne 0\end{array}$

for $j=1,2,\cdots ,p$. When ${\nu }_j=1$ the term is assumed to be active and will be included in the model, and when ${\nu }_j=0$ the term is assumed to be non-active and will not be included in the model.

Following [6], and [5], we assume that $\beta \mathrm{<mo>\; |\; </mo>}{\sigma }^2\mathrm{,}\nu \mathrm{,}c~N\left(0_p\mathrm{,}{\sigma }^2{D}_{\nu \mathrm{,}c}\right)$, where $\nu$ is the indicator vector, c is the prior variance of the slab distribution, and $D_{\nu \mathrm{,}c}$ is a diagonal matrix with the j^th diagonal element $cI\left({\nu }_j=1\right)+dI\left({\nu }_j=0\right),j=1,\cdots ,p$. The parameters ${\sigma }^2\mathrm{,}\nu$ and c will be given prior distributions, and the parameter d is assumed to be a small fixed non-negative number because we want the spike distribution to have a smaller variance than the slab distribution. Formally the prior construction of $\beta$ is the following:

${\beta }_j\mathrm{<mo>\; |\; </mo>}{\sigma }^2\mathrm{,}{\nu }_j\mathrm{,}c~\left(1-{\nu }_j\right)N\left(\mathrm{0,}d{\sigma }^2\right)+{\nu }_{j}N\left(\mathrm{0,}c{\sigma }^2\right)\mathrm{.}$

For every coefficient ${\beta }_j$, a Bernoulli variable ${\nu }_j$ is defined taking values 1 and 0 with probability of inclusion $\omega$, as $p\left({\nu }_j=1\right)=\omega$ and $p\left({\nu }_j=0\right)=\left(1-\omega \right)$. Often, ${\nu }_j$ ’s are taken as independent Bernoulli ( $\omega$ ) random variables, where $0<\omega <1$. It is common to fix $\omega$ in the normal mixture, however, we shall deal with $\omega$ as a parameter to investigate different values of $\omega$, and sample it from the Beta distribution as it will be explained in Section 6.3.

5. Prior Distributions

1) Following the prior distributions used by [5], we assume that the prior distribution for the fixed effects is $\beta ~N\left(\mathrm{0,}{\sigma }^2{D}_{\nu \mathrm{,}c}\right)$

$p\left(\beta \mathrm{<mo>\; |\; </mo>}{\sigma }^2\mathrm{,}\nu \mathrm{,}c\right)\propto {\left|{\sigma }^2{D}_{\nu \mathrm{,}c}\right|}^{-1/2}\text{exp}\left(-\frac{1}{2}{\beta }^{\prime }{\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}\beta \right)\mathrm{.}$

2) The prior distribution for the total variance is ${\sigma }^2~IG\left(a\mathrm{,}b\right)$,

$p\left({\sigma }^2\right)\propto {\left({\sigma }^2\right)}^{-a-1}\text{exp}\left(-\frac{b}{{\sigma }^2}\right)\mathrm{.}$

For this work, we used $a=0$ and $b=0$ following [5] as this yields the common non-informative prior for ${\sigma }^2$. This prior is improper, however we will sample from the posterior distribution, which should be a proper gamma distribution.

3) The prior distribution for the correlation parameter is $\rho ~Beta\left(a^{\prime }\mathrm{,}{b}^{\prime }\right)$ with shape parameters $a^{\prime },b^{\prime }>0$. We consider $a^{\prime }=b^{\prime }=2$ following [8]. According to [8], “A $Beta\left(a^{\prime }\mathrm{,}{b}^{\prime }\right)$ prior distribution for a correlation parameter can be interpreted as indicating a prior point estimate of $a^{\prime }/\left(a^{\prime }+b^{\prime }\right)$, this prior information being worth $a^{\prime }+b^{\prime }$ observations”. Our prior was selected to be centred at $2/\left(2+2\right)=0.5$ and to be worth four observations in each block. For an experiment with $a^{\prime }+b^{\prime }$ observations, the posterior distribution would give equal weight to the prior and the likelihood [8]. Similar choice of the prior is in [6] as they set $a=b=\text{2}$ so that the $p\left(\rho \right)$ is symmetric with a mode of 0.5. This has the effect of pulling the posterior mode of $p\left(\rho \right)$ towards 0.5 when the data are scarce. The prior density for $\rho$ is

$p\left(\rho \right)\propto {\rho }^{\left(a^{\prime }-1\right)}{\left(1-\rho \right)}^{\left(b^{\prime }-1\right)}\mathrm{.}$

4) The prior distribution for the elements of the indicator vector is ${\nu }_j~Bernoulli\left(\omega \right)$,

$p\left({\nu }_j\right)=(\begin{array}{l}\omega \text{if}{\nu }_j=1\\ 1-\omega \text{if}{\nu }_j=0\end{array}$

where $\omega$ is the prior probability that ${\beta }_j$ is active following [5].

5) The prior distribution for the elements of the probability of inclusion is $\omega ~Beta\left(c_0\mathrm{,}{d}_0\right)$,

$p\left(\omega \right)\propto {\omega }^{\left(c_0-1\right)}{\left(1-\omega \right)}^{\left(d_0-1\right)}\mathrm{.}$

[5] set $\omega =0.25$. However, we select $c_0$ and $d_0$ such that the prior of $\omega$ has a mode = 0.25. The choice of $c_0=2$ and $d_0=4$ results in a prior with a mode = 0.25 and the upper cumulative percentile at 5% equals 0.66. Meaning a 5% chance the observations have a probability density function ≥ 0.66.

6) The prior distribution for the slab variance c is a discrete uniform prior distribution with support points T = {1/4, 9/16, 1, 4, 9, 16, 25} as given by [5]. They found that large values of c tend to favor sparse models with large effects and in this case small effects will be missed. On the other hand, small values of c tend to favor less sparse models. Moreover, very small values of c tend to favor sparse models again. They select the support points in T such that it covers small and large values of c. The prior distribution for c is

$p\left(c\right)=(\begin{array}{l}\frac{1}{7}\text{if}c\in \text{T}\\ 0\text{otherwise}\end{array}$

6. Full Conditional Distributions

We use the prior distributions presented in Section 5 to derive the full conditional distributions. The likelihood of the data depends on $\beta$, so we can derive the conditional distribution for $\beta$ using the prior distribution $\beta \mathrm{<mo>\; |\; </mo>}{\sigma }^2\mathrm{,}\nu \mathrm{,}c~N\left(0_p\mathrm{,}{\sigma }^2{D}_{\nu \mathrm{,}c}\right)$ and the likelihood for the model 22

$L\left(y\mathrm{<mo>\; |\; </mo>}\theta \right)\propto {\left|V\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\left(y-X\beta \right)}^{\prime }{V}^{-1}\left(y-X\beta \right)\right]\mathrm{.}$

Note that we standardise both X and y so the fixed effect vector $\beta$ does not include the intercept. The conditional distribution for $\beta$ can be expressed as:

$\begin{array}{l}p\left(\beta \mathrm{<mo>\; |\; </mo>}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\omega \mathrm{,}\rho \mathrm{,}y\right)\propto p\left(y\mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\omega \mathrm{,}\rho \right)p\left(\beta \mathrm{<mo>\; |\; </mo>}{\sigma }^2\mathrm{,}\nu \mathrm{,}c\right)\\ \propto {\left|V\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\left(y-X\beta \right)}^{\prime }{V}^{-1}\left(y-X\beta \right)\right]\\ \times {\left|{\sigma }^2{D}_{\nu \mathrm{,}c}\right|}^{-1/2}\text{exp}\left(-\frac{1}{2}{\beta }^{\prime }{\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}\beta \right)\\ \propto {\left|V\right|}^{-1/2}{\left|{\sigma }^2{D}_{\nu \mathrm{,}c}\right|}^{-1/2}\text{exp}[-\frac{1}{2}\left({\beta }^{\prime }{\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}\beta \right)\\ -\frac{1}{2}\left(y^{\prime }{V}^{-1}y-y^{\prime }{V}^{-1}X\beta -{\beta }^{\prime }{X}^{\prime }{V}^{-1}y+{\beta }^{\prime }{X}^{\prime }{V}^{-1}X\beta \right)]\end{array}$

$\begin{array}{l}\propto {\left|V\right|}^{-1/2}{\left|{\sigma }^2{D}_{\nu \mathrm{,}c}\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\beta }^{\prime }\left({\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right)\beta \right]\\ \times \text{exp}\left[-\frac{1}{2}\left(-y^{\prime }{V}^{-1}X\right)\beta \right]\text{exp}\left[-\frac{1}{2}{\beta }^{\prime }\left(-X^{\prime }{V}^{-1}y\right)\right]\\ \propto {\left|V\right|}^{-1/2}{\left|{\sigma }^2{D}_{\nu \mathrm{,}c}\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\beta }^{\prime }\left({\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right)\beta \right]\\ \times \text{exp}\left[\frac{1}{2}\left(y^{\prime }{V}^{-1}X\right)\beta +\frac{1}{2}{\beta }^{\prime }\left(X^{\prime }{V}^{-1}y\right)\right]\\ \propto {\left|V\right|}^{-1/2}{\left|{\sigma }^2{D}_{\nu \mathrm{,}c}\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\beta }^{\prime }\left({\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right)\beta +{\beta }^{\prime }\left(X^{\prime }{V}^{-1}y\right)\right]\mathrm{.}\end{array}$

The key to deriving the joint posterior distribution is to rewrite the expression in the exponential part in a more convenient form. This can happen by using the multivariate completion of squares:

$\text{<msup> <mi> U </mi> <mo> ′ </mo> </msup> <mi> AU </mi>}-\text{<mi> 2 </mi> <msup> <mi> U </mi> <mo> ′ </mo> </msup>}\alpha ={\left(\text{U}-{\text{A}}^{-1}\alpha \right)}^{\prime }\text{A}\left(\text{U}-{\text{A}}^{-1}\alpha \right)-{\alpha }^{\prime }{\text{A}}^{-1}\alpha \mathrm{,}$

where A is a symmetric positive definite (hence invertible) matrix. We assume $\text{U}=\beta$, $\text{A}={\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X$, and $\alpha =X^{\prime }{V}^{-1}y$.

The conditional distribution for $\beta$ can be written as:

$\begin{array}{l}p\left(\beta \mathrm{<mo>\; |\; </mo>}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\omega \mathrm{,}\rho \mathrm{,}y\right)\\ \propto {\left|V\right|}^{-1/2}{\left|{\sigma }^2{D}_{\nu \mathrm{,}c}\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\beta }^{\prime }\left({\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right)\beta -2{\beta }^{\prime }\left(X^{\prime }{V}^{-1}y\right)\right]\\ \propto {\left|V\right|}^{-1/2}{\left|{\sigma }^2{D}_{\nu \mathrm{,}c}\right|}^{-1/2}\text{exp}[{\left[\beta -{\left({\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right)}^{-1}\left(X^{\prime }{V}^{-1}y\right)\right]}^{\prime }\\ \times \left[{\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right]\left[\beta -{\left({\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right)}^{-1}\left(X^{\prime }{V}^{-1}y\right)\right]]\\ \propto \text{exp}\left[-\frac{1}{2}{\left(\beta -{\beta }_{\ast }\right)}^{\prime }{\text{D}}_{\ast }^{-1}\left(\beta -{\beta }_{\ast }\right)\right]\mathrm{.}\end{array}$

Thus, we can sample $\beta$ from the conditional posterior $\text{N}\left({\beta }_{\ast },{\text{D}}_{\ast }\right)$, where

${\beta }_{\ast }={\left({\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right)}^{-1}\left(X^{\prime }{V}^{-1}y\right)\text{and}{\text{D}}_{\ast }={\left({\left({\sigma }^2{D}_{\nu \mathrm{,}c}\right)}^{-1}+X^{\prime }{V}^{-1}X\right)}^{-1}\mathrm{.}$ (24)

6.1. The Conditional Distribution for ρ

The likelihood of the data depends on $\rho$, so the conditional distribution for $\rho$ can be derived by

$\begin{array}{l}p\left(\rho \mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\omega \mathrm{,}y\right)\propto p\left(y\mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\omega \mathrm{,}\rho \right)p\left(\rho \right)\\ \propto {\left|V\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\left(y-X\beta \right)}^{\prime }{V}^{-1}\left(y-X\beta \right)\right]\times {\rho }^{\left(a^{\prime }-1\right)}{\left(1-\rho \right)}^{\left(b^{\prime }-1\right)}\mathrm{.}\end{array}$

We note here that the likelihood depends on $\rho$ through $V$ as in (4.6), so we can express $V$ as a function of $\rho$,

$\text{V}={\sigma }^2\left(1-\rho \right)\left(I_n+\frac{\rho }{1-\rho }Z{Z}^{\prime }\right)\mathrm{.}$

The conditional distribution for $\rho$ is a non-standard distribution that cannot be sampled directly. Therefore, we use the Metropolis-Hastings (M-H) rejection sampling. Our correlation parameter is $\rho \in \left(\mathrm{0,1}\right]$, and has a prior $\beta \left(a^{\prime }\mathrm{,}{b}^{\prime }\right)$. We apply the Random-Walk Metropolis-Hastings algorithm, and select a proposal distribution of log-normal distribution for the variance ratio $\eta$ where $\eta =f\left(\rho \right)=\frac{\rho }{1-\rho }$ with a mean equal to the current value of ${\eta }^t$ at iteration t and variance $s^2$. The choice of $s^2$ affects the jumping rule in the random walk proposal distribution. As we have one parameter to be updated in the random walk algorithm which is $\rho$, we follow [12] and [17] to set $s^2=g^2\Sigma$. The most efficient jump has a scale $g\approx 2.4/\sqrt{h}$ where h is the number of parameters which will be updated. In this work, we set $g=2.4$ and $h=1$ following [17], and we set $\Sigma =100$ as this yields an appropriate acceptance rate associated with the independent sampler of the ACF plot. Thus, $\eta \in \left(\mathrm{0,}\infty \right)$ and

$g\left(\eta \right)=\frac{1}{\eta \sqrt{2\pi s^2}}\exp \left[-\frac{1}{2s^2}{\left(\ln \eta -{\eta }^t\right)}^2\right]$

We can use $\eta =\frac{\rho }{1-\rho }$ as a transformation function between $\eta$ and $\rho$ as $\rho =\frac{\eta }{1+\eta }$ and the Jacobian function of $\rho$ is $J\left(\rho \right)=\frac{\text{d}\eta }{\text{d}\rho }=\frac{1}{{\left(1-\rho \right)}^2}$.

We draw a proposal value ${\eta }_{\ast }$ from a log-normal( ${\eta }^t\mathrm{,}{s}^2$ ) distribution, and the probability of accepting or rejecting ${\eta }_{\ast }$ is the minimum of 1 and the ratio r where r is

$r=\frac{p\left({\rho }_{\ast }\mathrm{<mo>\; |\; </mo>}\text{all}\right)}{p\left({\rho }^t\mathrm{<mo>\; |\; </mo>}\text{all}\right)}\times \frac{q\left({\eta }^t\mathrm{<mo>\; |\; </mo>}{\eta }_{\ast }\right)}{q\left({\eta }_{\ast }\mathrm{<mo>\; |\; </mo>}{\eta }^t\right)}\mathrm{,}$

which is equivalent to

$r=\frac{p\left({\rho }_{\ast }\mathrm{<mo>\; |\; </mo>}\text{all}\right)}{p\left({\rho }^t\mathrm{<mo>\; |\; </mo>}\text{all}\right)}\times \frac{q\left({\rho }^t\mathrm{<mo>\; |\; </mo>}{\rho }_{\ast }\right)J\left({\rho }^t\right)}{q\left({\rho }_{\ast }\mathrm{<mo>\; |\; </mo>}{\rho }^t\right)J\left({\rho }_{\ast }\right)}\mathrm{.}$

Our proposal ratio is

$\frac{q\left({\eta }^t\mathrm{<mo>\; |\; </mo>}{\eta }_{\ast }\right)}{q\left({\eta }_{\ast }\mathrm{<mo>\; |\; </mo>}{\eta }^t\right)}=\frac{{\eta }_{\ast }\exp \left[-\frac{1}{2s^2}{\left(\ln {\eta }^t-{\eta }_{\ast }\right)}^2\right]}{{\eta }^t\exp \left[-\frac{1}{2s^2}{\left(\ln {\eta }_{\ast }-{\eta }^t\right)}^2\right]}\mathrm{,}$

which is equivalent to

$\frac{q\left({\rho }^t\mathrm{<mo>\; |\; </mo>}{\rho }_{\ast }\right)J\left({\rho }^t\right)}{q\left({\rho }_{\ast }\mathrm{<mo>\; |\; </mo>}{\rho }^t\right)J\left({\rho }_{\ast }\right)}=\frac{\left(\frac{{\rho }_{\ast }}{1-{\rho }_{\ast }}\right)\exp \left[-\frac{1}{s^2}{\left(\ln \left(\frac{{\rho }^t}{1-{\rho }^t}\right)-\left(\frac{{\rho }_{\ast }}{1-{\rho }_{\ast }}\right)\right)}^2\right]\times \left|\frac{1}{{\left(1-{\rho }^t\right)}^2}\right|}{\left(\frac{{\rho }^t}{1-{\rho }^t}\right)\exp \left[-\frac{1}{s^2}{\left(\ln \left(\frac{{\rho }_{\ast }}{1-{\rho }_{\ast }}\right)-\left(\frac{{\rho }^t}{1-{\rho }^t}\right)\right)}^2\right]\times \left|\frac{1}{{\left(1-{\rho }_{\ast }\right)}^2}\right|}\mathrm{.}$

The ratio r can be expressed as

$\begin{array}{c}r=\frac{{\left|V\left({\rho }_{\ast }\right)\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\left(y-X\beta \right)}^{\prime }V{\left({\rho }_{\ast }\right)}^{-1}\left(y-X\beta \right)\right]\times {\left({\rho }_{\ast }\right)}^{\left(a^{\prime }-1\right)}{\left(1-{\rho }_{\ast }\right)}^{\left(b^{\prime }-1\right)}}{{\left|V\left({\rho }^t\right)\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\left(y-X\beta \right)}^{\prime }V{\left({\rho }^t\right)}^{-1}\left(y-X\beta \right)\right]\times {\left({\rho }^t\right)}^{\left(a^{\prime }-1\right)}{\left(1-{\rho }^t\right)}^{\left(b^{\prime }-1\right)}}\\ \times \frac{{\rho }_{\ast }\left(1-{\rho }^t\right)\exp \left[-\frac{1}{s^2}{\left(\ln \left(\frac{{\rho }^t}{1-{\rho }^t}\right)-\left(\frac{{\rho }_{\ast }}{1-{\rho }_{\ast }}\right)\right)}^2\right]\left|{\left(1-{\rho }^t\right)}^{-2}\right|}{{\rho }^t\left(1-{\rho }_{\ast }\right)\exp \left[-\frac{1}{s^2}{\left(\ln \left(\frac{{\rho }_{\ast }}{1-{\rho }_{\ast }}\right)-\left(\frac{{\rho }^t}{1-{\rho }^t}\right)\right)}^2\right]\left|{\left(1-{\rho }_{\ast }\right)}^{-2}\right|}\mathrm{.}\end{array}$ (25)

where $\text{V}\left({\rho }_{\ast }\right)={\sigma }^2\left(1-{\rho }_{\ast }\right)\left(I_n+\frac{{\rho }_{\ast }}{1-{\rho }_{\ast }}Z{Z}^{\prime }\right)$, and $\text{V}\left({\rho }^t\right)={\sigma }^2\left(1-{\rho }^t\right)\left(I_n+\frac{{\rho }^t}{1-{\rho }^t}Z{Z}^{\prime }\right)$.

6.2. The Conditional Distribution for σ²

The likelihood of the data depends on ${\sigma }^2$, so we can express the conditional distribution of ${\sigma }^2$ as

$\begin{array}{l}p\left({\sigma }^2\mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}\rho \mathrm{,}\nu \mathrm{,}c\mathrm{,}\omega \mathrm{,}y\right)\propto p\left(y\mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\omega \mathrm{,}\rho \right)p\left({\sigma }^2\right)\\ \propto {\left|V\right|}^{-1/2}\text{exp}\left[-\frac{1}{2}{\left(y-X\beta \right)}^{\prime }{V}^{-1}\left(y-X\beta \right)\right]\times {\left({\sigma }^2\right)}^{-a-1}\text{exp}\left(-\frac{b}{{\sigma }^2}\right)\mathrm{.}\end{array}$

We know that $\text{V}={\sigma }^2\left(1-\rho \right)\left(I_n+\frac{\rho }{1-\rho }Z{Z}^{\prime }\right)$, so the conditional posterior for ${\sigma }^2$ can be written as

$\begin{array}{l}p\left({\sigma }^2\mathrm{<mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\\ \propto {\left|\left(1-\rho \right){\sigma }^2\left(I_n+\frac{\rho }{1-\rho }Z{Z}^{\prime }\right)\right|}^{-1/2}{\left({\sigma }^2\right)}^{-a-1}\times \text{exp}\left(-\frac{b}{{\sigma }^2}\right)\\ \times \text{exp}\left[-\frac{1}{2}{\left(y-X\beta \right)}^{\prime }{\left(\left(1-\rho \right){\sigma }^2\left(I_n+\frac{\rho }{1-\rho }Z{Z}^{\prime }\right)\right)}^{-1}\left(y-X\beta \right)\right]\\ \propto {\left({\sigma }^2\right)}^{-\left(a+\frac{n}{2}\right)-1}\text{exp}\left(-\frac{1}{{\sigma }^2}\left[\frac{{\left(y-X\beta \right)}^{\prime }{\left(\left(1-\rho \right)\left(I_n+\frac{\rho }{1-\rho }Z{Z}^{\prime }\right)\right)}^{-1}\left(y-X\beta \right)}{2}+b\right]\right)\mathrm{.}\end{array}$

This is the inverse gamma distribution with a shape parameter $a^{\ast }$ and a scale parameter $b^{\ast }$ such that

$a^{\ast }=a+\frac{n}{2}\text{and}{b}^{\ast }=\frac{{\left(y-X\beta \right)}^{\prime }{\left(\left(1-\rho \right)\left(I_n+\frac{\rho }{1-\rho }Z{Z}^{\prime }\right)\right)}^{-1}\left(y-X\beta \right)}{2}+b\mathrm{.}$ (26)

The indicator vector can be drawn conditionally on the regressor coefficient and computation of the marginal likelihood is not required. The prior probabilities for ${\nu }_j$ are

$p\left({\nu }_j\right)=(\begin{array}{l}\omega \text{if}{\nu }_j=1\\ 1-\omega \text{if}{\nu }_j=0\end{array}$

where $\omega$ is the prior probability that ${\beta }_j$ is active. The joint conditional posterior distribution for $\nu$ has mass function

$\begin{array}{c}p\left(\nu \mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}y\right)\propto p\left(y\mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}\nu \right)p\left(\beta \mathrm{,}\nu \right)\\ =p\left(y\mathrm{<mo>\; |\; </mo>}\beta \right)p\left(\beta \mathrm{,}\nu \right)\\ \propto p\left(\beta \mathrm{,}\nu \right)\\ =p\left(\beta \mathrm{<mo>\; |\; </mo>}\nu \right)p\left(\nu \right)\mathrm{.}\end{array}$

The conditional density for $\beta$ given $\nu$ is

$\begin{array}{c}p\left(\beta \mathrm{<mo>\; |\; </mo>}\nu \right)\propto {\left|{\sigma }^2\text{diag}\left[c\text{I}\left({\nu }_j=1\right)+d\text{I}\left({\nu }_j=0\right)\right]\right|}^{-1/2}\\ \times \text{exp}\left[-\frac{1}{2}{\beta }^{\prime }{\left({\sigma }^2\text{diag}\left[c\text{I}\left({\nu }_j=1\right)+d\text{I}\left({\nu }_j=0\right)\right]\right)}^{-1}\beta \right]\mathrm{.}\end{array}$

The conditional distribution for the j^th component given ${\nu }_j$ is

$\begin{array}{c}p\left({\beta }_j\mathrm{<mo>\; |\; </mo>}{\nu }_j\right)\propto {\left|{\sigma }^2\left[c\text{I}\left({\nu }_j=1\right)+d\text{I}\left({\nu }_j=0\right)\right]\right|}^{-1/2}\\ \times \text{exp}\left[-\frac{{\beta }_j^2}{2{\sigma }^2\left[c\text{I}\left({\nu }_j=1\right)+d\text{I}\left({\nu }_j=0\right)\right]}\right]\mathrm{.}\end{array}$

The conditional posterior probabilities for ${\nu }_j$ are therefore

$\begin{array}{c}p\left({\nu }_j=1|{\beta }_j,y\right)=p\left({\nu }_j=1\right)p\left({\beta }_j|{\nu }_j=1\right)\\ \propto \omega {\left|c{\sigma }^2\right|}^{-1/2}\text{exp}\left[-\frac{{\beta }_j^2}{2c{\sigma }^2}\right]\mathrm{,}\end{array}$ (27)

and

$\begin{array}{c}p\left({\nu }_j=0|{\beta }_j,y\right)=p\left({\nu }_j=0\right)p\left({\beta }_j|{\nu }_j=0\right)\\ \propto \left(1-\omega \right){\left|d{\sigma }^2\right|}^{-\mathrm{1\; <mo>\; /\; </mo>2}}\text{exp}\left[-\frac{{\beta }_j^2}{2d{\sigma }^2}\right]\mathrm{.}\end{array}$ (eq:4.10)

6.3. The Conditional Distribution for ω

The probability of inclusion $\omega$ can be drawn conditionally on the indicator and computation of the marginal likelihood is not required. Hence the conditional distribution for $\omega$ is

$\begin{array}{c}p\left(\omega \mathrm{<mo>\; |\; </mo>}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\beta \mathrm{,}\rho \mathrm{,}y\right)\propto p\left(y\mathrm{<mo>\; |\; </mo>}\omega \mathrm{,}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\beta \mathrm{,}\rho \right)p\left(\omega \mathrm{,}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\beta \mathrm{,}\rho \right)\\ =p\left(y\mathrm{<mo>\; |\; </mo>}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\mathrm{,}\beta \mathrm{,}\rho \right)p\left(\omega \mathrm{,}\nu \right)\\ \propto p\left(\omega \mathrm{,}\nu \right)\\ =p\left(\nu \mathrm{<mo>\; |\; </mo>}\omega \right)p\left(\omega \right)\\ \propto {\omega }^{{\displaystyle {\sum }_{j=1}^p{\nu }_j}}{\left(1-\omega \right)}^{p-{\displaystyle {\sum }_{j=1}^p{\nu }_j}}\times {\omega }^{\left(c_0-1\right)}{\left(1-\omega \right)}^{\left(d_0-1\right)}\\ \propto {\omega }^{c_0+{\displaystyle {\sum }_{j=1}^p{\nu }_j}-1}{\left(1-\omega \right)}^{p-{\displaystyle {\sum }_{j=1}^p{\nu }_j}+d_0-1}\mathrm{.}\end{array}$

Hence,

$\omega \mathrm{<mo>\; |\; </mo>}\nu ~\text{Beta}\left(c_0+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\nu }_j,p-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\nu }_j+d_0\right)\mathrm{,}\text{where}p-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\nu }_j={\displaystyle \underset{j=1}{\overset{p}{\sum }}}\text{I}\left({\nu }_j=0\right)\mathrm{.}$ (28)

6.4. The Conditional Distribution for c

The prior distribution for c is a discrete uniform distribution with support points T = {1/4, 9/16, 1, 4, 9, 16, 25}, and it can be drawn conditionally on the regressor coefficient. The computation of the marginal likelihood is not required. Hence, the conditional distribution for c

$\begin{array}{c}p\left(c\mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}{\sigma }^2\mathrm{,}\nu \mathrm{,}\omega \mathrm{,}\rho \mathrm{,}y\right)\propto p\left(y\mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}{\sigma }^2\mathrm{,}\nu \mathrm{,}c\right)p\left(\beta \mathrm{,}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\right)\\ =p\left(y\mathrm{<mo>\; |\; </mo>}\beta \mathrm{,}{\sigma }^2\right)p\left(\beta \mathrm{,}\nu \mathrm{,}{\sigma }^2\mathrm{,}c\right)\\ \propto p\left(\beta \mathrm{,}{\sigma }^2\mathrm{,}\nu \mathrm{,}c\right)\\ =p\left(\beta \mathrm{<mo>\; |\; </mo>}{\sigma }^2\mathrm{,}\nu \mathrm{,}c\right)p(\; c\; )\end{array}$

$\begin{array}{l}\propto \frac{1}{7}{\left|{\sigma }^2\text{diag}\left[c\text{I}\left({\nu }_j=1\right)+d\text{I}\left({\nu }_j=0\right)\right]\right|}^{-1/2}\\ \times \text{exp}\left[-\frac{1}{2}{\beta }^{\prime }{\left({\sigma }^2\text{diag}\left[c\text{I}\left({\nu }_j=1\right)+d\text{I}\left({\nu }_j=0\right)\right]\right)}^{-1}\beta \right]\\ \propto \frac{1}{7}{\left[{\displaystyle \underset{j=1}{\overset{p}{\prod }}}\left[c\text{I}\left({\nu }_j=1\right)+d\text{I}\left({\nu }_j=0\right)\right]\right]}^{-1/2}\\ \times \text{exp}\left[-\frac{1}{2c}\left({\beta }^{\prime }{\text{I}}_{{\displaystyle {\sum }_j{\nu }_j}=1}\beta \right)\right]\times \text{exp}\left[-\frac{1}{2d}\left({\beta }^{\prime }{\text{I}}_{{\displaystyle {\sum }_j{\nu }_j}=0}\beta \right)\right]\end{array}$

$\begin{array}{l}\propto \frac{1}{7}{c}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}+\frac{1}{7}{d}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p\left(1-{\nu }_j\right)}}{2}}\times \text{exp}\left[-\frac{1}{2c}\left({\beta }^{\prime }{\text{I}}_{{\displaystyle {\sum }_j{\nu }_j}=1}\beta \right)\right]\\ \times \text{exp}\left[-\frac{1}{2d}\left({\beta }^{\prime }{\text{I}}_{{\displaystyle {\sum }_j{\nu }_j}=0}\beta \right)\right]\\ \propto \frac{1}{7}{c}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2c}\left({\beta }^{\prime }{\text{I}}_{{\displaystyle {\sum }_j{\nu }_j}=1}\beta \right)\right]\\ \propto \frac{1}{7}{c}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2c}\left({\beta }^{\prime }\beta \right)\right]\mathrm{.}\end{array}$

Therefore,

$p\left(c\mathrm{<mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\propto (\begin{array}{l}\frac{1}{7}{c}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2c}\left({\beta }^{\prime }\beta \right)\right]\text{if}c\in \text{T}\\ 0\text{otherwise}\end{array}$

Then, the posterior probabilities of the conditional distribution $p\left(c\mathrm{<mo>\; |\; </mo>}y\mathrm{,}\cdots \right)$ are

$p\left(c=\frac{1}{4}\mathrm{<mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\propto \frac{1}{7}{\left(\frac{1}{4}\right)}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2\left(\frac{1}{4}\right)}\left({\beta }^{\prime }\beta \right)\right]\mathrm{,}$

$p\left(c=\frac{9}{16}\mathrm{<mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\propto \frac{1}{7}{\left(\frac{9}{16}\right)}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2\left(\frac{9}{16}\right)}\left({\beta }^{\prime }\beta \right)\right]\mathrm{,}$

$p\left(c=\mathrm{1\; <mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\propto \frac{1}{7}{\left(1\right)}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2\left(1\right)}\left({\beta }^{\prime }\beta \right)\right]\mathrm{,}$

$p\left(c=\mathrm{4\; <mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\propto \frac{1}{7}{\left(4\right)}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2\left(4\right)}\left({\beta }^{\prime }\beta \right)\right]\mathrm{,}$

$p\left(c=\mathrm{9\; <mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\propto \frac{1}{7}{\left(9\right)}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2\left(9\right)}\left({\beta }^{\prime }\beta \right)\right]\mathrm{,}$

$p\left(c=\mathrm{16\; <mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\propto \frac{1}{7}{\left(16\right)}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2\left(16\right)}\left({\beta }^{\prime }\beta \right)\right]\mathrm{,}$

and

$p\left(c=\mathrm{25\; <mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\propto \frac{1}{7}{\left(25\right)}^{\frac{-{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\times \text{exp}\left[-\frac{1}{2\left(25\right)}\left({\beta }^{\prime }\beta \right)\right]\mathrm{.}$

The conditional posterior $p\left(c\mathrm{<mo>\; |\; </mo>}y\mathrm{,}\cdots \right)$ can be written as

$p\left(c\mathrm{<mo>\; |\; </mo>}y\mathrm{,}\cdots \right)\mathrm{<mo>\; =\; </mo>}(\begin{array}{l}\frac{c^{-\frac{{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}\text{exp}\left[-\frac{{\beta }^{\prime }\beta }{2c}\right]}{{\displaystyle {\sum }_{c\in T}{c}^{-\frac{{\displaystyle {\sum }_{j\mathrm{<mo>\; =\; </mo>1}}^p{\nu }_j}}{2}}}\text{exp}\left[-\frac{{\beta }^{\prime }\beta }{2c}\right]}\text{if}c\in \text{T}\\ 0\text{otherwise}\end{array}$ (29)

7. Bayesian Variable Selection Algorithms

In this section, we introduce the computational algorithms which we use in the Bayesian analysis for variable selection using the SSVS, and the SSVS-SPD. In this work, we choose an asymmetric proposal distribution, the log-normal density. We apply the Metropolis-Hastings algorithms to sample the variance ratio $\eta$ where $\eta =f\left(\rho \right)=\frac{\rho }{1-\rho }$, and $\rho$ is the correlation parameter. This is because of the fact that in our experiments, observations from different subplots within the same wholeplot are positively correlated as $\rho =\frac{{\sigma }_{\gamma }^2}{{\sigma }_{\epsilon }^2+{\sigma }_{\gamma }^2}$ ; also observations from different wholeplots are independent [5].

7.1. The Stochastic Search Variable Selection (SSVS) Algorithm

We process the MCMC estimation of the parameters $\beta \mathrm{,}\rho \mathrm{,}{\sigma }^2\mathrm{,}\nu \mathrm{,}\omega$ and c. We use the priors of all these parameters as in Section 5. The following Metropolis-Hastings within Gibbs sampling algorithm can be implemented. Let y be the $n\times 1$ vector of random responses, X is the $n\times p$ model matrix without the column of the intercept, $\beta$ is the $p\times 1$ vector of fixed effect parameters, where p is the number of fixed effect parameters that need to be estimated. We set initial values for the parameters as ${\beta }^{\left(0\right)}=1_p$, ${\nu }^{\left(0\right)}=1_p$, ${\rho }^{\left(0\right)}=0.5$, ${\sigma }^2{}^{\left(0\right)}=10$, $c^{\left(0\right)}=1$, ${\omega }^{\left(0\right)}=0.5$, $d=0.001$. Starting at the tth iteration such that $t=1,2,\cdots ,its$ where $its=10000$, and setting $j=1,2,\cdots ,p$, the sampling algorithm is:

1) For $j=1,2,\cdots ,p$, sample ${\nu }_j^{\left(t\right)}$ of the indicator vector ${\nu }^{\left(t\right)}$ using (27) for ${\beta }_j^{\left(t-1\right)}$, $c^{\left(t-1\right)}$, ${\sigma }^2{}^{\left(t-1\right)}$, and ${\omega }^{\left(t-1\right)}$.

2) Sample the mixture weight ${\omega }^{\left(t\right)}$ using (28) for ${\nu }^{\left(t\right)}$.

3) Sample the regressor coefficients ${\beta }^{\left(t\right)}$ using (24) for X, y, $V^{\left(t-1\right)}$, $D^{\left(t-1\right)}$, $c^{\left(t-1\right)}$, ${\nu }^{\left(t\right)}$, and ${\sigma }^2{}^{\left(t-1\right)}$.

4) Sample the total variance ${\sigma }^2{}^{\left(t\right)}$ using (26) for X, y, Z, ${\beta }^{\left(t\right)}$, and ${\rho }^{\left(t-1\right)}$.

5) a) Sample ${\rho }_{\ast }^{\left(t\right)}$ from $\beta \left(a^{\prime }\mathrm{,}{b}^{\prime }\right)$ ;

b) Calculate ${\alpha }^{\left(t\right)}$ using (25) for X, y, $V\left({\rho }_{\ast }^{\left(t\right)}\right)$, $V\left({\rho }^{\left(t-1\right)}\right)$, and ${\beta }^{\left(t\right)}$ ;

c) Sample $u^{\left(t\right)}$ from $U\left(\mathrm{0,1}\right)$ ;

d) If ${\alpha }^{\left(t\right)}>u^{\left(t\right)}$, then set ${\rho }^{\left(t\right)}={\rho }_{\ast }^{\left(t\right)}$, otherwise set ${\rho }^{\left(t\right)}={\rho }^{\left(t-1\right)}$.

6) Sample $c^{\left(t\right)}$ from the set T with probability mass function given in (29) for ${\beta }^{\left(t\right)}$, and ${\nu }^{\left(t\right)}$.

7.2. The Stochastic Search Variable Selection for Split-Plot Design (SSVS-SPD) Algorithm

We adapt the SSVS for the analysis of data from split-plot designs by taking into account the two types of factors, i.e. the whole-plot factors and the subplot factors which expected to have different effect sizes for the two strata in split-plot design [9]. This approach can be reported as the Stochastic Search Variable Selection for Split-Plot Design (SSVS-SPD). While the SSVS samples all parameters from one slab variance posterior distribution, the SSVS-SPD samples the whole-plot parameters and the subplot parameters from two different slab variance posterior distributions given that the whole-plot and the subplot effects might have different sizes. While the SSVS samples all parameters from one slab variance posterior distribution, the SSVS-SPD samples the whole-plot parameters and the subplot parameters from two different slab variance posterior distributions. We use the same priors as in the SSVS for all the parameters of interest as in Section 5. Basically, the SSVS-SPD can be seen as running the SSVS twice in one process, one for subplot factors and the other one for whole-plot factors. The algorithm can be explained as follows:

We process the MCMC estimation of the parameters $\beta \mathrm{,}\rho \mathrm{,}{\sigma }^2\mathrm{,}\nu \mathrm{,}\omega$ and c. The following Metropolis-Hastings within Gibbs sampling algorithm can be implemented. Let y be the $n\times 1$ vector of random responses, X is the $n\times p$ model matrix without the column of the intercept, X.S is the $n\times p_s$ model matrix for subplot factors where $p_s$ is the number of subplot fixed effect parameters.

Also, X.W is the $n\times p_w$ model matrix for whole-plot factors where $p_w$ is the number of whole-plot fixed effect parameters. The $\beta =\left({\beta }_s\mathrm{,}{\beta }_w\right)$ is the $p\times 1$ vector of fixed effect parameters, where p is the number of fixed effect parameters that need to be estimated, ${\beta }_s$ is the $p_s\times 1$ subplot effect parameters, and ${\beta }_w$ is the $p_w\times 1$ whole-plot effect parameters.

We set initial values for the parameters as ${\beta }_s^{\left(0\right)}=1_{p_s}$, ${\beta }_w^{\left(0\right)}=1_{p_w}$. The initial values for the indicator vectors for the subplot factor ${\nu }_s$ and the whole-plot factor ${\nu }_w$ are ${\nu }_s^{\left(0\right)}=1_{ps}$ and ${\nu }_w^{\left(0\right)}=1_{pw}$. Also, ${\rho }^{\left(0\right)}=0.5$, ${\sigma }^2{}^{\left(0\right)}=10$, and $d=0.001$. The initial values for the slab variance for the subplot factors $c_s$ and for the slab variance for the whole-plot factors $c_w$ are $c_s^{\left(0\right)}=c_w^{\left(0\right)}=1$. Finally, the initial weights for the subplot factors ${\omega }_s$ and for the whole-plot factors ${\omega }_w$ are ${\omega }_s^{\left(0\right)}={\omega }_w^{\left(0\right)}=0.5$.

Starting at the tth iteration such that $t=1,2,\cdots ,its$ where $its=10000$, and setting $j=1,2,\cdots ,p_s$ and $k=1,2,\cdots ,p_w$, the sampling algorithm is:

1) For $j=1,2,\cdots ,p_s$, and $k=1,2,\cdots ,p_w$ sample ${\nu }_{sj}^{\left(t\right)}$ and ${\nu }_{wk}^{\left(t\right)}$ of the indicator vectors ${\nu }_s^{\left(t\right)}$ and ${\nu }_w^{\left(t\right)}$ using (27) for ${\beta }_{sj}^{\left(t-1\right)}$, ${\beta }_{wk}^{\left(t-1\right)}$, $c_s^{\left(t-1\right)}$, $c_w^{\left(t-1\right)}$, ${\sigma }^2{}^{\left(t-1\right)}$, ${\omega }_s^{\left(t-1\right)}$, and ${\omega }_w^{\left(t-1\right)}$.

2) Allocate ${\nu }^{\left(t\right)}=\left({\nu }_s^{\left(t\right)}\mathrm{,}{\nu }_w^{\left(t\right)}\right)$.

3) Sample the mixture weights ${\omega }_s^{\left(t\right)}$ and ${\omega }_w^{\left(t\right)}$ using (28) for ${\nu }_s^{\left(t\right)}$ and ${\nu }_w^{\left(t\right)}$.

4) Allocate ${\omega }^{\left(t\right)}=\left({\omega }_s^{\left(t\right)}\mathrm{,}{\omega }_w^{\left(t\right)}\right)$.

5) Sample the regressor coefficients ${\beta }_s^{\left(t\right)}$ and ${\beta }_w^{\left(t\right)}$ using (24) for X, y, $V^{\left(t-1\right)}$, $D_s^{\left(t-1\right)}$, $D_w^{\left(t-1\right)}$, $c_s^{\left(t-1\right)}$, $c_w^{\left(t-1\right)}$, ${\nu }_s^{\left(t\right)}$, ${\nu }_w^{\left(t\right)}$, and ${\sigma }^2{}^{\left(t-1\right)}$. Where the $D_s$ is a diagonal matrix with the j^th diagonal element $c_s^{\left(t-1\right)}I\left({\nu }_{sj}=1\right)+dI\left({\nu }_{sj}=0\right)$, and $D_w$ is a diagonal matrix with the k^th diagonal element $c_w^{\left(t-1\right)}I\left({\nu }_{wk}=1\right)+dI\left({\nu }_{wk}=0\right)$.

6) Allocate ${\nu }^{\left(t\right)}=\left({\nu }_s^{\left(t\right)}\mathrm{,}{\nu }_w^{\left(t\right)}\right)$ and $D^{\left(t\right)}=\text{diag}\left(D_s^{\left(t\right)}\mathrm{,}{D}_w^{\left(t\right)}\right)$.

7) Sample the total variance ${\sigma }^2{}^{\left(t\right)}$ using (26) for X, y, Z, ${\beta }^{\left(t\right)}$, and ${\rho }^{\left(t-1\right)}$.

8) a) Sample ${\rho }_{\ast }^{\left(t\right)}$ from $\beta \left(a\mathrm{,}b\right)$ ;

b) Calculate ${\alpha }^{\left(t\right)}$ using (25) for X, y, $V\left({\rho }_{\ast }^{\left(t\right)}\right)$, $V\left({\rho }^{\left(t-1\right)}\right)$, and ${\beta }^{\left(t\right)}$ ;

c) Sample $u^{\left(t\right)}$ from $U\left(\mathrm{0,1}\right)$ ;

d) If ${\alpha }^{\left(t\right)}>u^{\left(t\right)}$, then set ${\rho }^{\left(t\right)}={\rho }_{\ast }^{\left(t\right)}$, otherwise set ${\rho }^{\left(t\right)}={\rho }^{\left(t-1\right)}$.

9) Sample $c_s^{\left(t\right)}$ and $c_w^{\left(t\right)}$ from the set T with probability mass function given in (29) for ${\beta }_s^{\left(t\right)}$, ${\beta }_w^{\left(t\right)}$, ${\nu }_s^{\left(t\right)}$ and ${\nu }_w^{\left(t\right)}$.

10) Allocate $c^{\left(t\right)}=\left(c_s^{\left(t\right)}\mathrm{,}{c}_w^{\left(t\right)}\right)$.

8. Real Data Application

[18] and [19] described an experiment in the production of vinyl for automobile seat covers. The experiment has 28 runs and is a modified version of an example in [19]. It involves the production of vinyl for automobile seat covers. In the experiment, the effects of five factors on the thickness of the vinyl are investigated. Three of the factors are mixture components and two of them are so-called process variables. As in ordinary mixture experiments, the component proportions sum to one. In this example, the response of interest does not only depend on these proportions, but also on the effects of the process variables. The mixture components in the experiment are three plasticizers whose proportions are represented by $s_1\mathrm{,}{s}_2$ and $s_3$. The two process variables studied are rate of extrusion ( $w_1$ ) and temperature of drying ( $w_2$ ). The experiment was conducted in a split-plot format. The process variables are the whole plot factors of the experiment, whereas the mixture components are the sub-plot factors. The data are shown in Table 1.

Table 1. Data for the Vinyl thickness for three-component mixture experiment with two process variable [20].

A main effects plus two factor interactions model was assumed for the process variable $w_1$ and $w_2$. For the mixture components, the quadratic mixture model was used. The main effects of the process variables were crossed with the linear blending terms only, so the model estimated in [18] is given by

$y_{ij}={\displaystyle \underset{i=1}{\overset{3}{\sum }}}{\beta }_i{s}_i+{\displaystyle \underset{i=1}{\overset{m-1}{\sum }}}{\displaystyle \underset{j=i+1}{\overset{m}{\sum }}}{\beta }_{ij}{s}_i{s}_j+{\displaystyle \underset{i=1}{\overset{2}{\sum }}}{\alpha }_i{w}_i+\alpha w_1{w}_2+{\displaystyle \underset{i=1}{\overset{2}{\sum }}}{\displaystyle \underset{j=1}{\overset{2}{\sum }}}{\delta }_{ij}{s}_i{w}_j+{\gamma }_i+{\epsilon }_{ij}$

They computed the variance components ${\sigma }^2={\sigma }_{\epsilon }^2+{\sigma }_{\gamma }^2$ by R. The ${\hat{\sigma }}_{\epsilon }^2=6.764$ and the ${\hat{\sigma }}_{\gamma }^2=0.00001$. The factor effect estimates will be displayed in the next section to compare it with the proposed methods in this paper.

9. Analysis of the Vinyl Thickness Experiment

The real dataset for the vinyl thickness experiment has been used to apply the Bayesian variable selection approach. The model involves 5 main variables ( $w_1\mathrm{,}{w}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\mathrm{,}{s}_3$ ) and two factor interaction variables ( $w_1{s}_1\mathrm{,}{w}_1{s}_2\mathrm{,}{s}_1{s}_2\mathrm{,}{s}_1{s}_3\mathrm{,}{s}_2{s}_3\mathrm{,}{w}_2{s}_1\mathrm{,}{w}_2{s}_2\mathrm{,}{w}_1{w}_2$ ). We used the prior distributions above. Following [21], in a Bayesian framework the final model could be the median probability model consisting of those variables whose posterior inclusion probability $p\left({\nu }_j=1|y\right)\ge 0.5$. The posterior probability of parameter ${\beta }_j,j=1,2,\cdots ,p$ being active is approximated by

${\displaystyle \underset{q=1}{\overset{its}{\sum }}}\frac{{\nu }_j^{\left(q\right)}}{its},$ (30)

where ${\nu }_j^{\left(q\right)}$ is ${\nu }_j$ sampled at iteration $q=1,\cdots ,its$ of the Metropolis-Hastings within Gibbs sampling algorithm.

We summarise the results of applying the Bayesian variable selection for the data from thickness vinyl experiment. The model which will be used is:

$y_{ij}={\displaystyle \underset{i=1}{\overset{3}{\sum }}}{\beta }_i{s}_i+{\displaystyle \underset{i=1}{\overset{m-1}{\sum }}}{\displaystyle \underset{j=i+1}{\overset{m}{\sum }}}{\beta }_{ij}{s}_i{s}_j+{\displaystyle \underset{i=1}{\overset{2}{\sum }}}{\alpha }_i{w}_i+\alpha w_1{w}_2+{\displaystyle \underset{i=1}{\overset{2}{\sum }}}{\displaystyle \underset{j=1}{\overset{2}{\sum }}}{\delta }_{ij}{s}_i{w}_j+{\gamma }_i+{\epsilon }_{ij}$

The estimates of the 13 parameters of the model have been reported using (SSVS) and (SSVS-SPD) for the response y which is displayed in Table 1 as well as the estimates by the generalised least estimator (GLS) for comparison purpose.

Figure 2 shows a comparison between SSVS and SSVS-SPD with respect to the resulting approximate posterior probability for the thickness vinyl experiment. The parameters ${\beta }_{w_1}\mathrm{,}{\beta }_{s_1}\mathrm{,}{\beta }_{s_2}\mathrm{,}{\beta }_{s_3}\mathrm{,}{\beta }_{w_1{s}_1}\mathrm{,}{\beta }_{w_1{s}_2}$ have the highest posterior probability of being active by both SSVS and SSVS-SPD. This indicates that the six associated variables to these terms play a significant role in this experiment. Followed by these terms, we find the parameters ${\beta }_{w_1{w}_2}\mathrm{,}{\beta }_{s_1{s}_2}\mathrm{,}{\beta }_{s_1{s}_3}\mathrm{,}{\beta }_{s_2{s}_3}$ have an approximate posterior probability of about 0.5 and 0.6 by both SSVS and SSVS-SPD. We note that SSVS and SSVS-SPD tend to consider ${\beta }_{s_1{s}_2}\mathrm{,}{\beta }_{s_1{s}_3}\mathrm{,}{\beta }_{s_2{s}_3}$ to be significant at an approximate posterior probability of 0.5 and 0.6 while they are not significant by the GLS method. All methods consider ${\beta }_{w_2}\mathrm{,}{\beta }_{w_2{s}_1}\mathrm{,}{\beta }_{w_2{s}_2}$ to be non significant with low approximate posterior probability in this experiment. The Bayesian analysis for the real data of thickness vinyl experiment yielded 10 significant variables. In contrast, the p value for the GLS estimates in Table 2 shows there are 7 significant variables as it excludes the variables associated to the coefficients ${\beta }_{s_1{s}_2}\mathrm{,}{\beta }_{s_1{s}_3}\mathrm{,}{\beta }_{s_2{s}_3}$. This means that the SSVS and SSVS-SPD yielded in extra 3 significant variables to the model. Table 2 represents the posterior means of the coefficients and standard deviation by both the SSVS and the SSVS-SPD methods and the GLS estimates with the p values for the thickness vinyl experiment. Table 3 shows the posterior mean of the correlation $\hat{\rho }$ and the posterior mean of the total variance ${\hat{\sigma }}^2$ for both SSVS and SSVS-SPD methods. Figure 3, shows the ACF plots for Markov Chain using Metropolis Hastings within Gibbs sampling algorithm used by SSVS and SSVS-SPD to sample the correlation $\rho$ and the variance ${\sigma }^2$.

Figure 2. Approximate posterior probability for the thickness vinyl experiment.

Table 2. The estimated coefficients and standard deviations (in parenthesis) for the thickness vinyl experiment by SSVS and SSVS-SPD. The first row is the estimates coefficeients and p-values (in parenthesis) for the GLS.

Table 3. Posterior means of the ${\hat{\sigma }}^2$ and the $\hat{\rho }$ by the SSVS and SSVS-SPD for the thickness vinyl experiment.

Figure 3. ACF plot for the Markov chain formed by sampling the total variance ${\sigma }^2$ and the correlation $\rho$ by SSVS and SSVS-SPD for the thickness vinyl experiment.

10. Simulation Study Using the Design of the Vinyl Thickness Experiment

We performed a simulation study by generating 1000 datasets, where each dataset would ran for 10000 iterations using MCMC. The SSVS and the SSVS-SPD would be applied at two levels of $\eta =1$ and $\eta =10$. We assume that ${\sigma }_{\epsilon }^2+{\sigma }_{\epsilon }^2=10$. Thus, the true value of the total variance ${\sigma }^2=10$. Also, the true value of $\rho$ is 0.5. We will calculate Type I and II error rates using the indicator vector $\nu$ and the approximate posterior probability in 30. If the true variable is active but the algorithm yielded a corresponding approximation posterior probability of less than 0.5, this variable would then have Type II error rate. Also, if the true variable is non-active but the algorithm yielded a corresponding approximation posterior probability larger than or equal to 0.5, this variable would then have Type I error rate. We also will calculate the precision of the point estimates by SSVS and SSVS-SPD by counting the median relative model error (MRME) for the estimates of the SSVS and SSVS-SPD. We focus on the properties of the estimated models by investigating the following properties:

1) consistency in variable selection (frequency in selecting the active/non-active variable), and

2) prediction performance.

For point 1, at 5% significant level, we report Type I error rate (an effect that is truly not significant but the corresponding procedure estimate indicates that it is significant). We also report Type II error rate (an effect that is truly present but the corresponding procedure estimate indicates that it is not significant).

For point 2, following [22] and [23], prediction accuracy is measured by computing the mean-squared error for each penalised estimate ${\hat{\beta }}_{\lambda }$ as,

$\text{ME}\left({\hat{\beta }}_{\lambda }\right)={\left(X{\hat{\beta }}_{\lambda }-X\beta \right)}^{\prime }\left(X{\hat{\beta }}_{\lambda }-X\beta \right)\mathrm{.}$

The relative model error (RME) is the ratio of the model error of the penalised estimates to the model error for the GLS estimates of the fixed effects,

$RME=\frac{ME\left({\hat{\beta }}_{\lambda }\right)}{ME\left({\hat{\beta }}_{GLS}\right)}\mathrm{,}$

where ${\hat{\beta }}_{GLS}$ in Equation (13) is the generalised least squares estimator of $\beta$. The median of the relative model error (MRME) over 1000 simulated data sets were reported. MRME values greater than one indicate that the methods estimates perform worse than the GLS estimates, values near to one indicate that the methods estimates performs in a similar way to the GLS estimates, values less than one indicate that the methods estimates performs better than the GLS estimates.

We perform a simulation study to examine the performance of the SSVS and SSVS-SPD methods. Using the design of the thickness vinyl we generate the response given the true model as:

$\mathbb{E}\left(Y\right)=4w_1-3s_1+s_3+4w_1{w}_2+2s_1{s}_2-s_1{s}_3+2w_1{s}_1+3w_2{s}_1$ (31)

Type I and II error rates are displayed in Table 4 and Table 5, for two setting of $\eta =1$ and $\eta =10$. Also, Table 6 represents the estimated Posterior means of ${\hat{\sigma }}^2$ and $\hat{\rho }$ by the SSVS and SSVS-SPD from the simulation by using the design of the thickness vinyl experiment. Figure 4 shows the MRME values at $\eta =1$ and $\eta =10$. We notes that the SSVS at $\eta =1$ have MRME greater than one which indicates that the estimates by the SSVS is worse than the GLS estimates. While at $\eta =10$, the SSVS perform better than the GLS estimates. The SSVS-SPD have similar performance with the GLS estimates at both level of $\eta$.

Type II error rates at both level of $\eta$ generally are low indicating that the active values are easy to detect by both SSVS and SSVS-SPD. Furthermore, Type II error rates by the SSVS-SPD are lower than Type II error rates by the SSVS. With regard to Type I error rate, for both level of $\eta$, the SSVS-SPD have lower Type I error rates than the SSVS. Detecting active variables by both SSVS and SSVS-SPD is better than detecting the non active variables. In summary, we found that the SSVS-SPD detects the active and non-active effect factors with a lower error rate than the SSVS. The analysis of the Median Relative Model Error (MRME) for both SSVS and SSVS-SPD using the thickness vinyl design provides similar behaviour to the GLS except at $\eta =10$ for the SSVS as detecting active effect factors was harder than the SSVS-SPD.

Table 4. Type I error rate for the simulation by using the design of thickness vinyl experiment.

Table 5. Type II error rate for the simulation by using the design of thickness vinyl experiment.

Figure 4. Median relative model error (MRME) for the simulation by using the design of thickness vinyl experiment.

Table 6. Posterior means of the ${\hat{\sigma }}^2$ and the $\hat{\rho }$ by the SSVS and SSVS-SPD for the simulation by using the design of thickness vinyl experiment.

11. Conclusion

This paper provided an analysis of data from split-plot mixture process experiments using a motivating example from the industrial environment. Specifically, we recommend the use of the SSVS-SPD method for Bayesian variable selection. In our results, we observed that the SSVS-SPD can identify the active variables (linear and two-factors interaction), much better than the SSVS and the traditional used GLS method. However, as expected this comes with the expense of slightly higher Type I error rates. We also observed a better prediction performance for the models chosen by the SSVS-SPD compared to the models chosen by the SSVS method.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Aljeddani, S. (2022) Bayesian Variable Selection for Mixture Process Variable Design Experiment. https://doi.org/10.21203/rs.3.rs-1906056/v1
[2]	Tang, X., Xu, X., Ghosh, M. and Ghosh, P. (2016) Bayesian Variable Selection and Estimation Based on Global-Local Shrinkage Priors. Sankhya A, arXiv: 1605.07981.
[3]	Mitchell, T.J. and Beauchamp, J.J. (1988) Bayesian Variable Selection in Linear Regression. Journal of the American Statistical Association, 83, 1023-1032. https://doi.org/10.1080/01621459.1988.10478694
[4]	George, E.I. and McCulloch, R.E. (1993) Variable Selection via Gibbs Sampling. Journal of the American Statistical Association, 88, 881-889. https://doi.org/10.1080/01621459.1993.10476353
[5]	Geweke. J. (1996) Variable Selection and Model Comparison in Regression. Bayesian Statistics 5, Valencia, 5-9 June 1994.
[6]	Tan, M.H.Y. and Wu, C.F.J. (2013) A Bayesian Approach for Model Selection in Fractionated Split Plot Experiments with Applications in Robust Parameter Design. Technometrics, 55, 359-372. https://doi.org/10.1080/00401706.2013.778790
[7]	George, E.I. and McCulloch, R.E. (1997) Approaches for Bayesian Variable Selection. Statistica Sinica, 7, 339-373.
[8]	Aljeddani, S. (2022) Variable Selection in Randomized Block Design Experiment. Journal of Computational Mathematics, 12, 216-231. https://doi.org/10.4236/ajcm.2022.122013
[9]	Gilmour, S.G. and Goos, P. (2009) Analysis of Data from Non-Orthogonal Multistratum Designs in Industrial Experiments. Journal of the Royal Statistical Society Series C: Applied Statistics, 58, 467-484. https://doi.org/10.1111/j.1467-9876.2009.00662.x
[10]	Jones, B. and Nachtsheim, C.J. (2009) Split-Plot Designs: What, Why, and How. Journal of Quality Technology, 41, 340-361. https://doi.org/10.1080/00224065.2009.11917790
[11]	Aljeddani, S. (2019) Statistical Analysis of Data from Experiments Subject to Restricted Randomisation. University of Southampton, Southampton.
[12]	Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. and Rubin, D.B. (2014) Bayesian Data Analysis. Vol. 2, CRC Press, Boca Raton. https://doi.org/10.1201/b16018
[13]	O’Hagan, A. and Forster, J. (2004) Kendall’s Advanced Theory of Statistics: Bayesian Statistics. Vol. 2B, Arnold, London.
[14]	Gilks, W.R., Richardson, S. and Spiegelhalter, D. (1995) Markov Chain Monte Carlo in Practice. CRC Press, Boca Raton. https://doi.org/10.1201/b14835
[15]	Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H. and Teller, E. (1953) Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21, 1087-1092. https://doi.org/10.1063/1.1699114
[16]	Hastings, W.K. (1970) Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika, 57, 97-109. https://doi.org/10.1093/biomet/57.1.97
[17]	Gelman, A., Roberts, G.O., Gilks, W.R., et al. (1996) Efficient Metropolis Jumping Rules. Bayesian Statistics, 5, 599-607.
[18]	Cornell, J.A., Kowalski, S.M. and Vining, G.G. (2002) Split-Plot Designs and Estimation Methods for Mixture Experiments with Process Variables. Technometrics, 44, 72-79. https://doi.org/10.1198/004017002753398344
[19]	Cornell, J.A. (2011) Experiments with Mixtures: Designs, Models, and the Analysis of Mixture Data. John Wiley & Sons, Hoboken.
[20]	Goos, P. and Donev, A.N. (2007) Tailor-Made Split-Plot Designs for Mixture and Process Variables. Journal of Quality Technology, 39, 326-339. https://doi.org/10.1080/00224065.2007.11917699
[21]	Barbieri, M.M., Berger, J.O. et al. (2004) Optimal Predictive Model Selection. The Annals of Statistics, 32, 870-897. https://doi.org/10.1214/009053604000000238
[22]	Fan, J. and Li, R. (2001) Variable Selection via Nonconcave Penalized Likelihood and Its Oracle Properties. Journal of the American Statistical Association, 96, 1348-1360. https://doi.org/10.1198/016214501753382273
[23]	Ibrahim, J.G., Zhu, H., Garcia, R.I. and Guo, R. (2011) Fixed and Random Effects Selection in Mixed Effects Models. Biometrics, 67, 495-503. https://doi.org/10.1111/j.1541-0420.2010.01463.x