Bayesian $D$-Optimal Design of Experiments with Quantitative and Qualitative Responses

Assume the non-informative priors $p(\bm{\beta}^{(i)})\propto 1$ for $i=1,2$ and $p(\bm{\eta})\propto 1$. The conditional posterior distribution in (7) is the same as the joint distribution of the data. It can be further factorized into

$$p(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X})\propto p(\bm{\eta}|\bm{z},\bm{X})\prod_{i=1}^{2}p(\bm{\beta}^{(i)}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X}),$$

with the posterior distributions

		$\displaystyle\bm{\beta}^{(i)}|\bm{y},\bm{z},\mu_{0},\sigma^{2},\bm{X}\sim N\left((\bm{F}^{\prime}\bm{V}_{i}\bm{F})^{-1}\bm{F}^{\prime}\bm{V}_{i}(\bm{y}-\mu_{0}{\bf 1}),\sigma^{2}(\bm{F}^{\prime}\bm{V}_{i}\bm{F})^{-1}\right)\textrm{ for }i=1,2,$		(9)
		$\displaystyle p(\bm{\eta}|\bm{z},\bm{X})\propto\exp\left\{\sum_{i=1}^{n}\left(z_{i}\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}-\log(1+e^{\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}})\right)\right\}.$		(10)

Conditioning on $\bm{z}$, the posterior distributions of $(\bm{\beta}^{1},\bm{\beta}^{2})$ and $\bm{\eta}$ are independent. Note that the noninformative prior $p(\bm{\eta})\propto 1$ is proper because it leads to proper posterior $p(\bm{\eta}|\bm{z},\bm{X})$. Under the noninformative priors, the Bayesian estimation is identical to the frequentist estimation.

Using the posterior distributions (9)–(10) and the criterion function $\psi(\cdot)=\log(\cdot)$ in the general Bayesian optimal design criterion (8), we obtain the Bayesian $D$-optimal design criterion (11).

	$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})$	$\displaystyle=\mathbb{E}_{\bm{z},\bm{\eta}}\left\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\}$		(11)
		$\displaystyle+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\{\log\det\{(\bm{F}^{\prime}\bm{V}_{i}\bm{F})\}\right\}+\textrm{constants}.$

The derivation is in Appendix A1. The first additive term in (11) is exactly the Bayesian $D$-optimal design criterion for GLMs. Unfortunately, its exact integration is not tractable. The common approach in experimental design for GLMs is to use a normal approximation for the posterior distribution $p(\bm{\eta}|\bm{z},\bm{X})$ (Chaloner and Verdinelli, 1995; Khuri et al., 2006). Such an approximation leads to

$$\mathbb{E}_{\bm{z},\bm{\eta}}\left\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\}\approx\mathbb{E}_{\bm{\eta}}\{\log\det\bm{I}(\bm{\eta}|\bm{X})\}+\textrm{constant},$$

(12)

where $\bm{I}(\bm{\eta}|\bm{X})$ is the Fisher information matrix. We can easily show that

$$\bm{I}(\bm{\eta}|\bm{X})=-\mathbb{E}_{\bm{z}}\left(\frac{\partial^{2}l(\bm{z},\bm{\eta}|\bm{X})}{\partial\bm{\eta}\partial\bm{\eta}^{T}}\right)=\sum_{i=1}^{n}\bm{f}(\bm{x}_{i})\bm{f}(\bm{x}_{i})^{\prime}\pi(\bm{x}_{i},\bm{\eta})(1-\pi(\bm{x}_{i},\bm{\eta}))=\bm{F}^{\prime}\bm{W}_{0}\bm{F},$$

where $\bm{W}_{0}$ is a diagonal weight matrix

$$\bm{W}_{0}=\mbox{diag}\{\pi(\bm{x}_{1},\bm{\eta})(1-\pi(\bm{x}_{1},\bm{\eta})),\ldots,\pi(\bm{x}_{n},\bm{\eta})(1-\pi(\bm{x}_{n},\bm{\eta}))\}.$$

Omitting the irrelevant constant, we approximate the exact criterion $\Psi(\bm{X}|\mu_{0},\sigma^{2})$ in (11) as follows.

$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})\approx\mathbb{E}_{\bm{\eta}}\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\}+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\{\log\det(\bm{F}^{\prime}\bm{V}_{i}\bm{F})\right\}.$

(13)

To construct the optimal design, we consider maximizing the approximated $\Psi(\bm{X}|\mu_{0},\sigma^{2})$ in (13). But this is not trivial, because it involves the expectation on $Z_{i}$’s in the second additive term. To overcome this challenge, Hainy et al. (2013) constructed optimal designs by simulating samples from the joint distribution of responses and the unknown parameters. But this method can be computationally expensive for even slightly larger dimensions of experimental factors. Instead of simulating $Z_{i}$’s, we derive the following Theorem 1 that gives a tractable upper bound $Q(\bm{X})$. Thus we propose using the upper bound $Q(\bm{X})$ as the optimal criterion.

The proof of Theorem 1 is in Appendix A1. Note that Theorem 1 requires that $\bm{F}^{\prime}\bm{W}_{i}\bm{F}$ for $i=0,1,2$ are all nonsingular. Obviously $\bm{W}_{0}=\bm{W}_{1}\bm{W}_{2}$. It is easy to see that $\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ is nonsingular if and only if both $\bm{F}^{\prime}\bm{W}_{1}\bm{F}$ and $\bm{F}^{\prime}\bm{W}_{2}\bm{F}$ are nonsingular.

The matrices $\bm{F}^{\prime}\bm{V}_{1}\bm{F}$ and $\bm{F}^{\prime}\bm{V}_{2}\bm{F}$ involve the responses $Z_{i}$’s that are not yet observed at the experimental design stage. We can only choose the experimental run size and the design points to avoid the singularity problem with a larger probability for given values of $\bm{\eta}$. Once the run size is chosen, the design points can be optimally arranged by maximizing $Q(\bm{X})$. The weight matrix $\bm{W}_{1}$ (or $\bm{W}_{2}$) gives more weight to the feasible design points that are more likely to lead to $Z=1$ (or $Z=0$) observations so that the parameters $\bm{\beta}^{(1)}$ (or $\bm{\beta}^{(2)}$) of the linear model $Y|Z=1$ (or $Y|Z=0$) are more likely to be estimable. Next, we introduce some regularity conditions on the run size and number of replications to alleviate the singularity problem.

Let $m$ be the number of distinct design points in the design matrix $\bm{X}$, $n_{i}$ be the number of repeated point $\bm{x}_{i}$ in $\bm{X}$ for $i=1,\ldots,m$. Thus $n=\sum_{i=1}^{m}n_{i}$ and $n_{i}\geq 1$ for $i=1,\ldots,m$. First, it is necessary that $m\geq q$ for the linear regression model to be estimable under the noninformative priors. The if-and-only-if condition for $\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ to be nonsingular is that $\textrm{rank}(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\geq q$. If $m\geq q$ and $\pi(\bm{x}_{i},\bm{\eta})\in(0,1)$ for $i=1,\ldots,m$, then $\bm{F}^{\prime}\bm{W}_{i}\bm{F}$ for $i=0,1,2$ are all nonsingular and thus positive definite. To make sure $\pi(\bm{x}_{i},\bm{\eta})\in(0,1)$ for $i=1,\ldots,m$, it is sufficient to assume that $\bm{\eta}$ is finitely bounded. This condition is typically used for the frequentist $D$-optimal design for GLMs. For instance, Woods et al. (2006) suggested using the centroids of the finite bounded space of $\bm{\eta}$ to develop the local $D$-optimal design for GLMs. Dror and Steinberg (2006) clustered different local $D$-optimal designs with $\bm{\eta}$ randomly sampled from its bounded space.

The if-and-only-if condition for $\bm{F}^{\prime}\bm{V}_{1}\bm{F}$ and $\bm{F}^{\prime}\bm{V}_{2}\bm{F}$ being nonsingular is

$$\sum_{i=1}^{m}I(\sum_{j=1}^{n_{i}}Z_{ij}>0)\geq q\quad\textrm{and}\quad\sum_{i=1}^{m}I(\sum_{j=1}^{n_{i}}Z_{ij}<n_{i})\geq q,$$

(15)

where $Z_{ij}$ is the $j$th random binary response at the unique design point $\bm{x}_{i}$ and $I(\cdot)$ is the indicator function. In the following, we discuss how to choose sample sizes $n_{i}$ and $n$ under two scenarios (i) $m=q$ and (ii) $m>q$.

Proposition 1 gives a sufficient condition on the lower bound of $n_{i}$ when saturated design ($m=q$) is used. Under the sufficient condition, the nonsingularity of $\bm{F}\bm{V}_{1}\bm{F}$ and $\bm{F}\bm{V}_{2}\bm{F}$ holds with a probability larger than $\kappa^{m}$. For Example 1 in Section 1, suppose that if the possible values of $x$ can only be $-1,0,1$, then $m=q=3$. Let $\bm{\eta}=(1,1)^{\prime}$. If $\kappa=0.5$, then the numbers of replications for $x=-1,0,1$ need to satisfy $n_{1}\geq 2$, $n_{2}\geq 4$, and $n_{3}\geq 7$, respectively. If $\kappa=0.9$, then $n_{1}\geq 5$, $n_{2}\geq 9$, and $n_{3}\geq 20$. Proposition 1 is useful in Step 1 to construct the initial design in Algorithm 2 in Section 5.

Proposition 2 provides one sufficient condition and one necessary condition when $m>q$ on the smallest number of replications and the overall run size. But these conditions only examine the nonsingularity of the two matrices with large probability, which is weaker than Proposition 1. For given $\bm{\eta}$ value, Algorithm 2 in Section 5.1 can return the local $D$-optimal design. Proposition 2 can be useful to check the local $D$-optimal design, as $\pi_{\min}$ and $\pi_{\max}$ depend on $\bm{\eta}$. Take the artificial example in Section 6.1 for instance. The local $D$-optimal design for $\rho=0$ (Table A1 in Appendix A2) has $m=50$ unique design points and there are $q=22$ effects. According to Proposition 2, the sufficient condition requires $n_{0}\geq 7$ and the necessary condition requires $n_{0}\geq 1$. The local $D$-optimal design in Table A1 only satisfies the necessary condition. To meet the sufficient condition $n$ has to be much larger. For the global optimal design considering all possible $\bm{\eta}$ values, Proposition 2 can provide some guidelines for the design construction when $\bm{\eta}$ is varied in a relatively small range.

For the parameters $\bm{\beta}^{(1)}$ and $\bm{\beta}^{(2)}$, the conjugate priors are normal distribution since $Y|Z$ follows normal distribution. Thus we consider their priors as

$$\bm{\beta}^{(1)}\sim N(\bm{0},\tau^{2}\bm{R}_{1}),\quad\bm{\beta}^{(2)}\sim N(\bm{0},\tau^{2}\bm{R}_{2}).$$

where $\tau^{2}$ is the prior variance and $\bm{R}_{i}$ is the prior correlation matrix of $\bm{\beta}^{(i)}$ for $i=1,2$. Here we use the same prior variance $\tau^{2}$ only for simplicity. The matrix $\bm{R}_{i}$ can be specified flexibly such as using $(\bm{F}^{\prime}\bm{F})^{-1}$, or those in Joseph (2006) for factorial designs.

For the parameters $\bm{\eta}$, we choose the conjugate prior derived in Chen and Ibrahim (2003). It takes the form

$$\bm{\eta}\sim D(s,\bm{b})\propto\exp\left\{\sum_{i=1}^{n}s\left(b_{i}\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}-\log(1+e^{\bm{f}(\bm{x}_{i})^{\prime}\bm{\eta}})\right)\right\},$$

(22)

where $D(s,\bm{b})$ is the distribution with parameters $(s,\bm{b})$. Here $s$ is a scalar factor and $\bm{b}\in(0,1)^{n}$ is the marginal mean of $\bm{z}$ as shown in Diaconis and Ylvisaker (1979). The value of $\bm{b}$ can be interpreted as a prior prediction (or guess) for $\mathbb{E}(\bm{Z})$. Based on the priors for $(\bm{\beta}^{(1)},\bm{\beta}^{(2)},\bm{\eta})$ we can derive the posteriors as follows.

The proof of Proposition 3 can be derived following the standard Bayesian framework, thus is omitted. To derive the Bayesian $D$-optimal design criterion, we take the posterior distributions in Proposition 3 to (8) and set $\psi(\cdot)=\log(\cdot)$. The derivation is very similar to that in (11), and thus we obtain the exact design criterion as

	$\displaystyle\Psi(\bm{X}|\mu_{0},\sigma^{2})$	$\displaystyle=\mathbb{E}_{\bm{z},\bm{\eta}}\left\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\}$		(25)
		$\displaystyle+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\{\log\det(\bm{F}^{\prime}\bm{V}_{i}\bm{F}+\rho\bm{R}_{i}^{-1})\right\}+\textrm{constant}.$

As the integration of $\mathbb{E}_{\bm{z},\bm{\eta}}\left\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\}$ is not tractable, we adopt the same normal approximation of the posterior distribution $p(\bm{\eta}|\bm{z},\bm{X})$ as in (12). A straightforward calculation leads to getting Fisher information matrix $\bm{I}(\bm{\eta}|\bm{X})=(1+s)\bm{F}^{\prime}\bm{W}_{0}\bm{F}$. Thus we have

$$\mathbb{E}_{\bm{z},\bm{\eta}}\left\{\log(p(\bm{\eta}|\bm{z},\bm{X}))\right\}\approx\mathbb{E}_{\bm{\eta}}\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\}+\textrm{constant}.$$

(26)

Disregarding the constant, we can approximate the exact $\Psi(\bm{X}|\mu_{0},\sigma^{2})$ by

$$\Psi(\bm{X}|\mu_{0},\sigma^{2})\approx\mathbb{E}_{\bm{\eta}}\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\}+\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}\mathbb{E}_{\bm{z}|\bm{\eta}}\left\{\log\det(\bm{F}^{\prime}\bm{V}_{i}\bm{F}+\rho\bm{R}_{i}^{-1})\right\},$$

The following Theorem 2 gives an upper bound of the approximated criterion $\Psi(\bm{X}|\mu_{0},\sigma^{2})$ to avoid the integration with respect to $\bm{z}$, which plays the same role as Theorem 1.

For the same argument as in Section 3, we use the upper bound in (27) as the optimal design criterion. Note that since $\rho\bm{R}_{i}^{-1}$ is added to $\bm{F}^{\prime}\bm{V}_{i}\bm{F}$ and $\bm{F}^{\prime}\bm{W}_{i}\bm{F}$, $\bm{F}^{\prime}\bm{V}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}$ and $\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}$ are nonsingular. The derivation of $\Psi(\bm{X}|\mu_{0},\sigma^{2})$ and $Q(\bm{X})$ only needs $\bm{F}^{\prime}\bm{W}_{0}\bm{F}$ to be nonsingular, which requires $m\geq q$ and $\bm{\eta}$ to be finitely bounded as in Section 3.

Note that the criterion in (27) has a similar formulation with $Q(\bm{X})$ in (14). The only difference is that (14) does not involve $\rho\bm{R}_{i}^{-1}$. For consistency, we use the formula (27) as the design criterion $Q(\bm{X})$ for both cases. When noninformative priors for $\bm{\beta}^{(1)}$ and $\bm{\beta}^{(2)}$ are used, we set $\rho=0$. From another point of view, as $\tau^{2}\rightarrow\infty$, $\rho\rightarrow 0$, the variances in the priors $p(\bm{\beta}_{1})$ and $p(\bm{\beta}_{2})$ diffuse and result in a noninformative priors.

The criterion $Q(\bm{X})$, consisting of three additive terms, can be interpreted intuitively. The first additive term $\mathbb{E}_{\bm{\eta}}\{\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\}$ is known as the Bayesian $D$-optimal criterion for logistic regression and $\mathbb{E}_{\bm{\eta}}\{\log\det(\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1})\}$ is the Bayesian $D$-optimal criterion for the linear regression model of $Y$. To explain the weights, we rewrite $Q(\bm{X})$ as follows.

$$Q(\bm{X})=1\cdot\mathbb{E}_{\bm{\eta}}\left(\log\det(\bm{F}^{\prime}\bm{W}_{0}\bm{F})\right)+1\cdot\left(\frac{1}{2}\sum_{i=1}^{2}\mathbb{E}_{\bm{\eta}}(\log\det(\bm{F}^{\prime}\bm{W}_{i}\bm{F}+\rho\bm{R}_{i}^{-1}))\right)$$

Since there are equal numbers of binary and continuous response observations, the design criterion should put the same weight (equal to 1) on both design criteria for $Z$ and $Y$. For the two criteria for the linear regression models, the same weight 1/2 is used. This is also reasonable because we assume $\pi_{i}\in(0,1)$. Then none of the diagonal entries of $\bm{W}_{1}$ and $\bm{W}_{2}$ are zero, so the two terms should split the total weight 1 assigned for the entire linear regression part. Therefore, even though $Q(\bm{X})$ are derived analytically, all the additive terms and their weights make sense intuitively.

Note that the conjugate prior $p(\bm{\eta})$ requires prior parameters $(s,\bm{b})$ to be specified. Moreover, the prior distribution (22) contains $f(\bm{x}_{i})$, which depends on the design points. When sampling $\bm{\eta}$ from the prior (22), it does not matter whether $f(\bm{x}_{i})$’s are actually from the design points. If relevant historical data is available, we can simply sample $\bm{\eta}$ from the likelihood of the data. Alternatively, one can adopt the method in Chen and Ibrahim (2003) to estimate the parameters $(s,\bm{b})$. Without the relevant data, we would use the noninformative prior for $\bm{\eta}$, i.e., $p(\bm{\eta})\propto 1$ in the bounded region for $\bm{\eta}$.

The design criterion $Q(\bm{X})$ contains some unknown parameters, including the noise-to-signal ratio $\rho=\sigma^{2}/\tau^{2}$ and the correlation matrices $\bm{R}_{i}$’s for $i=1,2$. The value of $\rho$ has to be specified either from the historical data or from the domain knowledge. Typically we would assume $\rho<1$ such that the measurement error has a smaller variance than the signal variance.

The setting of $\bm{R}_{i}$ can also be specified flexibly. If historical data are available, $\bm{R}_{i}$ can be set as the estimated correlation matrix of $\bm{\beta}^{(i)}$. Otherwise, we can use the correlation matrix in Joseph (2006) and Kang and Joseph (2009), which is targeted for factorial design. Specifically, let $\bm{\beta}$ be the unknown coefficients of the linear regression model and the prior distribution is $\bm{\beta}\sim N(\bm{0},\tau^{2}\bm{R})$. For 2-level factor coded in $-1$ and 1, Joseph (2006) suggests that $\bm{R}$ is a diagonal matrix and the priors for individual $\beta_{j}$ is

		$\displaystyle\beta_{0}\sim N(0,\tau^{2}),$		(28)
		$\displaystyle\beta_{j}\sim N(0,\tau^{2}r),\qquad i=1,\ldots,p,$
		$\displaystyle\beta_{j}\sim N(0,\tau^{2}r^{2}),\qquad i=p+1,\ldots,p+\binom{p}{2},$
		$\displaystyle\vdots$
		$\displaystyle\beta_{2^{p}-1}\sim N(0,\tau^{2}r^{p}),$

where $\beta_{j}$ $i=1,\ldots,p$ are main effects, $\beta_{j}$ $j=p+1,\ldots,p+\binom{p}{2}$ are 2-factor-interactions and up to the $p$-factor-interaction $\beta_{2^{p}-1}$. The variance of $\beta_{j}$ decreases exponentially with the order of their corresponding effects by $r\in(0,1)$, thus it incorporates the effects hierarchy principle (Wu and Hamada, 2011). Joseph (2006) showed that if $\bm{f}(\bm{x})$ contains all the $2^{p}$ effects of all $p$ orders, $\tau^{2}\bm{R}$ can be represented alternatively by Kronecker product as $\tau^{2}\bm{R}=\varsigma^{2}\bigotimes_{j=1}^{p}\bm{F}_{j}(x_{j})^{-1}\bm{\Psi}_{j}(x_{j})(\bm{F}_{j}(x_{j}))^{-1}$. The model matrix for the 2-level factor and the correlation matrix are

$${\bm{F}_{j}(x_{j})}=\left(\begin{array}[]{cc}1&-1\\ 1&1\end{array}\right)\;\;\textrm{and}\;\;\bm{\Psi}_{j}(x_{j})=\left(\begin{array}[]{cc}1&\zeta\\ \zeta&1\end{array}\right).$$

(29)

To keep the two different presentations equivalent, let $\zeta=\frac{1-r}{1+r}$ and $\tau^{2}=(\frac{1+\zeta}{2})^{p}\varsigma^{2}$. For the mixed-level of 2- and 3-level experiments, Kang and Joseph (2009) have extended the 2-level case to the format

$$\tau^{2}\bm{R}=\varsigma^{2}\bigotimes_{j=1}^{p_{2}+p_{3,c}+p_{3,q}}\bm{F}_{j}(x_{j})^{-1}\bm{\Psi}_{j}(x_{j})(\bm{F}_{j}(x_{j})^{-1})^{\prime},$$

(30)

where $p_{2}$ is the number of 2-level factors, $p_{3,c}$ is the number of 3-level qualitative (categorical) factors, and $p_{3,q}$ is the number of 3-level quantitative factors. For all the 3-level factors, the model matrix is

$${\bm{F}_{j}(x_{j})}=\left(\begin{array}[]{ccc}1&-\sqrt{\frac{3}{2}}&\sqrt{\frac{1}{2}}\\ 1&0&-\sqrt{2}\\ 1&\sqrt{\frac{3}{2}}&\sqrt{\frac{1}{2}}\end{array}\right).$$

(31)

An isotropic correlation function is recommended for the 3-level qualitative factors and a Gaussian correlation function for quantitative factors. Thus, the correlation matrices for the 3-level qualitative and quantitative factors are

$${\bm{\Psi}_{j}(x_{j})}=\left(\begin{array}[]{ccc}1&\zeta&\zeta\\ \zeta&1&\zeta\\ \zeta&\zeta&1\end{array}\right)\;\;\textrm{and}\;\;\bm{\Psi}_{j}(x_{j})=\left(\begin{array}[]{ccc}1&\zeta&\zeta^{4}\\ \zeta&1&\zeta\\ \zeta^{4}&\zeta&1\end{array}\right),$$

(32)

respectively. To keep the covariance (30) consistent with the 2-level case we still set $\zeta=\frac{1-r}{1+r}$. To keep the variance of the intercept equal to $\tau^{2}$ (Kang and Joseph, 2009), we set

$$\tau^{2}=\varsigma^{2}\left(\frac{1+\zeta}{2}\right)^{p_{2}}\left(\frac{1+2\zeta}{3}\right)^{p_{3,c}}\left(\frac{3+4\zeta+2\zeta^{4}}{9}\right)^{p_{3,q}},$$

and thus

$$\bm{R}=\left(\left(\frac{1+\zeta}{2}\right)^{p_{2}}\left(\frac{1+2\zeta}{3}\right)^{p_{3,c}}\left(\frac{3+4\zeta+2\zeta^{4}}{9}\right)^{p_{3,q}}\right)^{-1}\times\bigotimes_{j=1}^{p_{2}+p_{3,c}+p_{3,q}}\bm{F}_{j}(x_{j})^{-1}\bm{\Psi}_{j}(x_{j})(\bm{F}_{j}(x_{j})^{-1})^{\prime}.$$

It is straightforward to prove that $\bm{R}$ is a diagonal matrix if only 2-level and 3-level qualitative factors are involved, but not so if any 3-level quantitative factors are involved, and the first diagonal entry of $\bm{R}$ is always 1.