Blocking first order response surface designs with interaction under correlated error

Abstract

This paper deals with blocking of first order response surface model with interaction under equi- and auto-correlated error structure. The conditions for orthogonal blocking, orthogonal estimation of the parameters and constancy of the variance of the predicted response have been derived for the proposed model. D-optimality criterion for orthogonally blocked first order response surface designs with interaction has also been discussed. A method of construction of design satisfying the derived conditions is given under the proposed model along with an illustration.

Keywords

Blocking in response surface model orthogonal estimation equi-correlated errors auto-correlated errors

1. Introduction

Consider the functional form $y_{u}=f(x_{iu},\bm{\theta})+\varepsilon_{u}$ used in response surface methodology to explain the functional relationship between ‘ $N$ ’ experimental runs on ‘ $v$ ’ input variables $x_{1},x_{2},\ldots,x_{v}$ and single response $y_{u},u=1,2,\ldots,N$ . The parameter vector $\bm{\theta}$ is estimated using ordinary least square method and $\epsilon_{u}$ is the error term associated with the $u^{\rm th}$ observation assumed to follow normal distribution with mean zero and constant variance $\sigma^{2}$ (Box and Draper, (1987)).

One of the most important assumptions in response surface methodology is that the experimental units are homogeneous. If the experimental units are not homogeneous then the experiment is conducted in groups or in blocks, so that the experimental units within each group/block are homogeneous. The blocking is done by protecting the block effects in the model that are orthogonal to the parameters to be estimated. Khuri (1994) studied the effect of blocking in estimating the mean response, prediction variance and optimum value of the response surface. He found that the prediction variance increases as a result of blocking. However, for an orthogonally blocked design, the least squares estimates of the model parameters remain unchanged, except possibly for the intercept.

Blocking when observations are correlated has some real-time applications. For example, in order to find out the optimum process temperature and time of thermally processed fish product, quality characteristics resulting from containers viz: aluminium can and retort pouch can be predicted as a function of process temperature and time. The quality response variables measured during different process intervals are found to be correlated. The containers can be treated as groups/blocks in order to accommodate one type of variability.

Another example is the optimization of process conditions for the isolation of targeted microorganisms from the fish sample collected from different locations. The process parameters may be inoculation medium concentration and time. The fish samples collected from different locations may have different initial loads of microorganisms. The response measured on different combination of process conditions tend to be correlated. The response variable namely total plate count can be estimated as a function of inoculation medium concentration and time in groups/blocks of samples from different locations.

Das and Park (2009) developed a measure of robust slope rotatability for second order response surface experimental designs. They have derived conditions for slope rotatability with correlated error. There have been several studies on response surface designs under correlated error by Panda and Das (1994); Das (2003a, 2003b) and Varghese et al. (2013). They all derived the rotatability conditions for first and second order response surface models under correlated errors. Rajyalaksmi and Victorbabu (2015) studied second order slope rotatable designs under intra-class correlated error structure using balanced incomplete block designs with $v=$ 7, $b=$ 7, $r=$ 3, $k=$ 3 and $\lambda=$ 1.

Mann et al. (2010) studied the robustness aspects of the blocked designs of experiments by taking into account correlated error in the linear models and proposed a generalized least squares estimator to construct robust designed experiments with blocks. Varghese and Jaggi (2011) studied orthogonal blocking of first order response surface models, incorporating neighbor effects for overcoming the heterogeneity among experimental units and obtained conditions for orthogonal estimation of the parameters of the model. Myers and Montgomery (2002) mentioned orthogonal blocking in second order response surface designs and derived conditions for orthogonal blocking in second order response surface designs. They also pointed out the orthogonality conditions are true for central composite designs, Box-Behnken designs and two level factorial designs under orthogonal blocking. Goos and Donev (2006) studied the effect of blocking on response surface designs and derived optimality and orthogonal conditions for blocked experiments.

The present study focused on the orthogonal blocking/grouping of first order response surface models with interaction under correlated error structure. The conditions for orthogonal estimation of model parameters from the blocked experiment with correlated error structure have been derived. The effect of blocking on the prediction variance of the response variables is also estimated and the constancy of variance of estimated parameters is derived.

In the next section, statistical model with correlated error is introduced. The conditions for orthogonal estimation of the parameters, constancy of the variance of the predicted response and method of design construction and illustration for equi- and auto-correlated error are given in Sections 3 and 4, respectively. Effect of blocking on the prediction variance and D-optimality under correlated error is discussed respectively in Sections 5 and 6. In Section 7, analysis of simulated data with real time application is given. Concluding remarks are given in Section 8.

2. Response surface model with blocks under correlated error

A first order response surface model with interaction for ‘ $v$ ’ input factors and ‘ $b$ ’ blocks can be written as

$y_{lu}=\beta_{0}+\sum\nolimits_{i=1}^{v}{x_{ilu}+}\sum{\sum\nolimits_{i<j=2}^{% v}{\beta_{ij}x_{ilu}x_{jlu}}}+\sum\nolimits_{l=1}^{b}{\delta_{l}}B_{lu}+% \varepsilon_{u},1\leqslant i\leqslant v,1\leqslant u\leqslant N,$

$\bm{Y}=\beta_{0}{\bf 1}_{N}+\bm{X}_{1}\bm{\beta}+\bm{B}\bm{\delta}+\bm{e},$ (1)

where $\bm{Y}$ is a response variable of order $N\times 1$ , $\beta_{0}$ is intercept, ${\bf 1}_{N}$ is a column vector of 1’s of order $N\times 1$ , ${\bm{X}}_{1}$ is a design matrix of order $N\times(v+{v\choose 2})$ with ${\bm{X}}_{1}=\left(x_{1},x_{2},\ldots,x_{v};x_{1}\otimes x_{2},\ldots,x_{\left% (v-1\right)}\otimes x_{v}\right),x_{i}=\left(x_{i1},x_{i2},\ldots,x_{iN}\right% )^{\prime},x_{i}\otimes x_{j}=\left(x_{i1}x_{j1},x_{i2}x_{j2},\ldots,x_{iN}x_{% jN}\right)^{\prime};$ $1\leqslant i,j\leqslant v$ , and $\otimes$ denotes the Hadamard product as defined by Das and Parker (2006), $\bm{\beta}={(\beta_{1},\ldots,\beta_{v},\beta_{12},\ldots,\beta_{\left(v-1% \right)v})}^{\prime}$ is a $(v+{v\choose 2})\times 1$ vector of regression coefficients; ${\bm{\delta}}=\left(\delta_{1},\delta_{2},\ldots,\delta_{b}\right)^{\prime};$ $\delta_{l}$ is the effect of $l^{\rm th}$ block ( $l=1,2,\ldots,b$ ) and ${\bm{B}}$ is the block diagonal matrix of the form ${\bm{B}}=\text{diag}({\bf 1}_{n_{1}},{\bf 1}_{n_{2}},\ldots,{\bf 1}_{n_{b}})$ , where $n_{l}$ is the size of the $l^{\rm th}$ block such that $N=\sum\nolimits_{l=1}^{b}n_{l}$ . Further, term ${\bm{e}}$ is an $N\times 1$ vector of errors which follows N-variate normal distribution with $E({\bm{e}})=0$ and $D({\bm{e}})=\bm{V}$ with $\text{rank}(\bm{V})=N$ .

The above model in Eq. 1 is not of full rank since the sum of columns of ${\bm{B}}$ is ${\bf 1}_{N}$ , that is $\bm{B}{\bf 1}_{b}={\bf 1}_{N}$ . Therefore, the model can be re-written as

$\bm{Y}=\bm{X}_{1}\bm{\beta}+\bm{B}\bm{\tau}+\bm{e},$ (2)

where ${\bm{\tau}}=\beta_{0}{\bf 1}_{b}+{\bm{\delta}}$ . If the columns of ${\bm{B}}$ are linearly independent of those of ${\bm{X}}_{1}$ , then above model is of full rank. Thus ${\bm{\beta}}$ and ${\bm{\tau}}$ can be uniquely estimated by ordinary least squares method. It is not possible to estimate $\beta_{0}$ independent of ${\bm{\delta}}$ unless certain constraint is imposed on the element of ${\bm{\delta}}$ (Khuri and Cornell (1987)). For this purpose we may assume $\sum\nolimits_{l=1}^{b}\delta_{l}=0$ , so that $\beta_{0}=\frac{1}{b}\sum\nolimits_{l=1}^{b}\tau_{l}=\frac{1}{b}{\bf 1}_{{b}}^% {\prime}{\bm{\tau}}$ .

The Eq. 2 can be written as

$\bm{Y}=\bm{X}\bm{\theta}+\bm{e},\text{ where }\bm{X}=[\bm{X}_{1}:\bm{B}]\text{% and }{\bm{\theta}}=\left({\bm{\beta}}^{\prime}:{\bm{\tau}}^{\prime}\right).$ (3)

Suppose the elements of ${\bm{e}}$ are correlated and ‘ ${\bm{V}}$ ’ is known, an extension to usual error structure, the procedure for estimation of parameters will be different. The dependence of error components might have occurred due to some ordering of experimental units in the block. The usual estimation procedure for model Eq. 1 mentioned above will not be valid in the presence of correlated errors.

The best linear unbiased estimator of ${\bm{\theta}}$ under correlated error is $\hat{\bm{\theta}}$ , given by

$\hat{\bm{\theta}}=\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right)^{-1}{\bm{X% }}^{\prime}\bm{V}^{-1}{\bm{Y}}\text{ and }\bm{D}(\hat{\bm{\theta}})=\left({\bm% {X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right)^{-1}\sigma^{2}.$ (4)

We partition $\bm{D}(\hat{\bm{\theta}})$ as follows, in order to simplify the computations under different error structures

$\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right)^{-1}=\begin{bmatrix}{\bm{M}}% _{1}&{\bm{M}}_{2}&{\bm{M}}_{3}\\ {\bm{M}}_{2}^{\prime}&{\bm{M}}_{4}&{\bm{M}}_{5}\\ {\bm{M}}_{3}^{\prime}&{\bm{M}}_{5}^{\prime}&{\bm{M}}_{6}\\ \end{bmatrix}_{\left(v+{v\choose 2}+b\right)\times(v+{v\choose 2}+b)},$

where ${\bm{M}}_{1},{\bm{M}}_{4}$ and ${{\bm{M}}}_{6}$ are symmetric matrices and all the matrices are defined below:

$\displaystyle{\bm{M}}_{1}={{\bm{x}}}_{i}^{\prime}\bm{V}^{-1}{\bm{x}}_{j}=\left% (v_{i.j}\right)_{v\times v};1\leqslant i,j\leqslant v,$ $\displaystyle{\bm{M}}_{2}={{\bm{x}}}_{i}^{\prime}\bm{V}^{-1}\left({\bm{x}}_{j}% \otimes{\bm{x}}_{s}\right)=\left(v_{i.js}\right)_{v\times{v\choose 2}};1% \leqslant i,j<s\leqslant v,$ $\displaystyle{\bm{M}}_{3}={\bm{x}}_{i}^{\prime}\bm{V}^{-1}\bm{B}={\left(v_{i.b% }\right)}_{v\times b};1\leqslant i\leqslant v,$ $\displaystyle{\bm{M}}_{4}=\left({\bm{x}}_{i}\otimes{\bm{x}}_{j}\right)^{\prime% }\bm{V}^{-1}\left({\bm{x}}_{s}\otimes{\bm{x}}_{t}\right)=\left(v_{ij.st}\right% )_{{v\choose 2}\times{v\choose 2}};1\leqslant i,s<j,t\leqslant v,$ $\displaystyle{\bm{M}}_{5}=\left({\bm{x}}_{i}\otimes{\bm{x}}_{j}\right)^{\prime% }\bm{V}^{-1}\bm{B}={\left(v_{ij.b}\right)}_{{v\choose 2}\times b};1\leqslant i% ,j\leqslant v,$ $\displaystyle{\bm{M}}_{6}={\bm{B}}^{\prime}\bm{V}^{-1}\bm{B}={\left(v_{b.b}% \right)}_{b\times b}.$

Thus, the estimated value of ${\bm{Y}}$ at the point ${\bm{g}(\bm{x})}^{\prime}$ is obtained as

$\hat{y}(x)={{\bm{g}(\bm{x})}}^{\prime}\hat{\theta},\text{∼{}where∼{}}{{\bm{g}(% \bm{x})}}^{\prime}=[\left(x_{{10}},x_{20},\ldots,x_{v0},{x}_{10}x_{20},\ldots,% x_{(v-1)0}x_{v0}\right)^{\prime}:\frac{1}{b}{\bf 1}_{b}^{\prime})]$ (5)

and the prediction variance of $\hat{y}(x)$ at a point ${{\bm{g}(\bm{x})}}^{\prime}$ is given by

$V\left(\hat{y}(x)\right)={{\bm{g}(\bm{x})}}^{\prime}\left({\bm{X}}^{\prime}\bm% {V}^{-1}{\bm{X}}\right)^{-1}{\bm{g}}\left({\bm{x}}\right)\sigma^{2}.$ (6)

In the next section, we simplify the discussion when error vector has a specific correlation structure.

3. Blocking of response surface models under equi-correlated error structure

The $N\times N$ matrix of equi-correlation structure ( $\bm{V}_{{\bm{E}}}$ ) of errors ‘ ${\bm{e}}$ ’ given in the model Eq. 3 is given by

$\bm{V}_{{\bm{E}}}=\left(1-\rho\right){\bm{I}}_{N}+\rho{\bm{J}}_{N\times N}% \text{ and }\bm{V}_{{\bm{E}}}^{-1}=\left[\left(a-b\right){\bm{I}}_{N}+b{\bm{J}% }_{N\times N}\right],$

where $\rho$ is correlation coefficient, ${\bm{I}}$ is an identity matrix of order $N\times N$ and ${\bm{J}}$ is a matrix of order $N\times N$ , which contains ‘ $b$ ’ block diagonal matrices of order $n\times n$ of 1’s and all other elements are zeros. Further

$a=\frac{1+(N-2)\rho}{(1-\rho)(1+\left(N-1\right)\rho)},b=-\frac{\rho}{\left(1-% \rho\right)\left(1+\left(N-1\right)\rho\right)}\text{∼{}and∼{}}\rho>-\frac{1}{% (N-1)}.$

The generalized least square estimator of ${\bm{\theta}}$ is $\hat{\bm{\theta}}=\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right)^{-1}{\bm{X% }}^{\prime}\bm{V}^{-1}{\bm{Y}}$ , here

${\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}=\frac{1}{1-\rho}\left[\left({\bm{X}}^{% \prime}{\bm{I}\bm{X}}\right)-\frac{\rho}{1+(n-1)\rho}{\bm{X}}^{\prime}{\bm{J}% \bm{X}}\right].$ (7)

3.1 Conditions for orthogonal blocking and parameter estimation

The moment matrix $\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right)^{-1}$ should be diagonal for the orthogonal estimation of parameter ${\bm{\theta}}$ . We attain this objective by imposing the following restrictions;

1.
$\displaystyle v_{i.b}=\sum\nolimits_{u=1}^{n}x_{iu}=0,1\leqslant i\leqslant v$ and $\displaystyle v_{ij.b}=\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}}=0,1\leqslant i% \neq j\leqslant v$ .
2.
$\displaystyle v_{i.j}=b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}}-\frac{\rho}{1+% \left(n-1\right)\rho}b\sum\nolimits_{u=1}^{n}x_{iu}\sum\nolimits_{u=1}^{n}x_{% ju}=0,1\leqslant i\neq j\leqslant v$ .
3.
$\displaystyle v_{ii.j}=b\sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}}-\frac{\rho}{% 1+\left(n-1\right)\rho}b\sum\nolimits_{u=1}^{n}x_{iu}\sum\nolimits_{u=1}^{n}{x% _{iu}x_{ju}}=0,1\leqslant i\neq j\leqslant v$ .
4.
$\displaystyle v_{ij.s}=b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}x_{su}}-\frac{\rho% }{1+\left(n-1\right)\rho}b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}}\sum\nolimits_{% u=1}^{n}x_{su}=0,1\leqslant i<j,s\leqslant v$ .
5.
$\displaystyle v_{ii.js}=b\sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}x_{su}}-\frac% {\rho}{1+\left(n-1\right)\rho}b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}}\sum% \nolimits_{u=1}^{n}{x_{iu}x_{su}}=0,1\leqslant i,j<s\leqslant v$ .
6.
$\displaystyle v_{ij.st}=b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}x_{su}x_{tu}}-% \frac{\rho}{1+\left(n-1\right)\rho}b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}}\sum% \nolimits_{u=1}^{n}{x_{su}x_{tu}}=0,1\leqslant i,s<j,t\leqslant v$ .

These restrictions ensure the rotatability of the experimental design ${\bm{X}}_{1}$ under orthogonal blocks. Now, ${\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}$ can be written as

${\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}=\frac{1}{1-\rho}\text{diag}\begin{pmatrix% }\displaystyle b\sum\nolimits_{u=1}^{n}x_{1u}^{2},\ldots,b\sum\nolimits_{u=1}^% {n}x_{iu}^{2},\ldots,b\sum\nolimits_{u=1}^{n}x_{vu}^{2},\\ \displaystyle b\sum\nolimits_{u=1}^{n}{x_{1u}^{2}x_{2u}^{2}},\ldots,b\sum% \nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}^{2}},\ldots,b\sum\nolimits_{u=1}^{n}{x_{% \left(v-1\right)u}^{2}x_{vu}^{2}},\\ \displaystyle\frac{n\left(1-\rho\right)}{1+\left(n-1\right)\rho},\ldots,\frac{% n\left(1-\rho\right)}{1+\left(n-1\right)\rho}\\ \end{pmatrix}.$ (8)

Then, the parameters ( ${\bm{\theta}}$ ’s) are estimated from the normal equation $\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right){\bm{\theta}}={{\bm{X}}}^{% \prime}\bm{V}^{-1}{\bm{Y}}$

$\displaystyle\text{i.e., ∼{}∼{}}\begin{bmatrix}{\bm{S}}_{1}&0^{\prime}\\ 0&{\bm{S}}_{2}\\ \end{bmatrix}\begin{bmatrix}{\bm{\beta}}\\ {\bm{\tau}}\\ \end{bmatrix}=\frac{1}{1-\rho}\begin{bmatrix}{\bm{T}}\\ {\bm{B}}\\ \end{bmatrix},$ (9)

where ${\bm{T}}=[T_{1},\ldots,T_{i},\ldots,T_{v},T_{12},\ldots,T_{ij},\ldots,T_{\left% (v-1\right)v}]^{\prime}$ and ${\bm{B}}={[B_{1},\ldots,B_{l},\ldots,B_{b},]}^{\prime}$ are the vectors of treatment combination totals and block totals, respectively;

$\displaystyle T_{i}=\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=1}^{n}{x_{iu}y_{lu% }};T_{ij}=\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}y_{lu}}% \text{∼{}and∼{}}$ $\displaystyle B_{l}=\sum\nolimits_{u=1}^{n}y_{lu}-\frac{\rho}{1+\left(n-1% \right)\rho}n\sum\nolimits_{u=1}^{n}y_{lu}.$ $\displaystyle{\bm{S}}_{1}=\frac{1}{1-\rho}\text{diag}\begin{Bmatrix}% \displaystyle b\sum\nolimits_{u=1}^{n}x_{1u}^{2},\ldots,b\sum\nolimits_{u=1}^{% n}x_{iu}^{2},\ldots,b\sum\nolimits_{u=1}^{n}x_{vu}^{2},\\ \displaystyle b\sum\nolimits_{u=1}^{n}{x_{1u}^{2}x_{2u}^{2}},\ldots,b\sum% \nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}^{2}},\ldots,b\sum\nolimits_{u=1}^{n}{x_{% \left(v-1\right)u}^{2}x_{vu}^{2}}\\ \end{Bmatrix},$ $\displaystyle{\bm{S}}_{2}=\frac{1}{1-\rho}\text{diag}\begin{Bmatrix}% \displaystyle\frac{n\left(1-\rho\right)}{1+\left(n-1\right)\rho},\ldots,\frac{% n\left(1-\rho\right)}{1+\left(n-1\right)\rho}\end{Bmatrix}.$

Hence,

$\displaystyle\hat{\bm{\theta}}=\begin{bmatrix}\hat{\bm{\beta}}\\ \hat{\bm{\tau}}\\ \end{bmatrix}=\begin{bmatrix}{\bm{S}}_{1}^{-1}{\bm{T}}\\ {\bm{S}}_{2}^{-1}{\bm{B}}\\ \end{bmatrix}\text{ and }{\bm{D}}\begin{pmatrix}\hat{\bm{\beta}}\\ \hat{\bm{\tau}}\\ \end{pmatrix}=\sigma^{2}\begin{bmatrix}{\bm{S}}_{1}^{-1}&0^{\prime}\\ 0&{\bm{S}}_{2}^{-1}\\ \end{bmatrix}.$ (10)

The variance of the estimated parameters is obtained as

$\displaystyle{V(\hat{\beta}}_{i})=\sigma^{2}\left(1-\rho\right)\frac{1}{% \displaystyle b\sum\nolimits_{u=1}^{n}x_{iu}^{2}-\frac{\rho}{1+\left(n-1\right% )\rho}b\left(\sum\nolimits_{u=1}^{n}x_{iu}\right)^{2}},$ (11) $\displaystyle{V(\hat{\beta}}_{ij})={\sigma}^{2}\left(1-\rho\right)\frac{1}{% \displaystyle b\sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}^{2}}-\frac{\rho}{% \displaystyle 1+\left(n-1\right)\rho}{b\left(\sum\nolimits_{u=1}^{n}x_{iu}x_{% ju}\right)}^{2}}$ (12)

and

$\displaystyle{V(\hat{\tau}}_{l})=\sigma^{2}\left(1-\rho\right)\frac{1+\left(n-% 1\right)\rho}{n\left(1-\rho\right)}={\sigma}^{2}\frac{1+\left(n-1\right)\rho}{% n}.$ (13)

The estimated response at the point ${{\bm{g}(\bm{x})}}^{\prime}$ and its prediction variance via Eq. 6 is given by

$V\left(\hat{y}(x)\right)=(1-\rho)\sigma^{2}\begin{Bmatrix}\displaystyle\frac{x% _{{10}}^{2}}{\displaystyle b\sum\nolimits_{u=1}^{n}x_{1u}^{2}}+\ldots+\frac{x_% {i0}^{2}}{\displaystyle b\sum\nolimits_{u=1}^{n}x_{iu}^{2}}+\ldots+\frac{x_{v0% }^{2}}{\displaystyle b\sum\nolimits_{u=1}^{n}x_{vu}^{2}}+\\ \displaystyle\frac{x_{{10}}^{2}x_{{20}}^{2}}{\displaystyle b\sum\nolimits_{u=1% }^{n}{x_{1u}^{2}x_{2u}^{2}}}+\ldots+\frac{x_{i0}^{2}x_{j0}^{2}}{\displaystyle b% \sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}^{2}}}+\ldots+\frac{x_{\left(v-1\right% )0}^{2}x_{v0}^{2}}{\displaystyle b\sum\nolimits_{u=1}^{n}{x_{\left(v-1\right)u% }^{2}x_{vu}^{2}}}+\\ \displaystyle\frac{1+\left(n-1\right)\rho}{bn\left(1-\rho\right)}\\ \end{Bmatrix}.$ (14)

The constancy of the variance of estimated parameters is ensured by putting the following conditions

1.
$\displaystyle\sum\nolimits_{u=1}^{n}x_{iu}^{2}=g$ ,
2.
$\displaystyle\sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}^{2}}=h,$

where, $g$ and $h$ are constants for $\forall i\neq j=1,2,\ldots,v$ for all blocks of size ‘ $n$ ’. Therefore,

$V\left(\hat{y}(x)\right)=(1-\rho)\sigma^{2}\left\{\frac{(1+\left(n-1\right)% \rho)}{bn(1-\rho)}+\frac{\displaystyle\sum\nolimits_{i=1}^{v}x_{i0}^{2}}{bg}+% \frac{\displaystyle\sum\nolimits_{i\neq j=1}^{v}{x_{i0}^{2}x_{j0}^{2}}}{bh}% \right\}.$ (15)

This shows that the variances of estimators of linear regression coefficients $\hat{\beta}_{i}$ ’s and interaction regression coefficients $\hat{\beta}_{ij}$ ’s $(i\neq j=1,2,\ldots,v)$ are same and the variance of the estimated response is a function of $\sum\nolimits_{i=1}^{v}x_{i0}^{2}$ and $\sum\nolimits_{i\neq j=1}^{v}{x_{i0}^{2}x_{j0}^{2}}$ for $\hat{\beta}_{i}$ ’s and $\hat{\beta}_{ij}$ ’s, respectively. For a given $\rho$ , the estimated response will have the same variance under equi-correlated error structure.
3.2 Method of construction and illustration

We considered a method of construction given by Varghese and Jaggi (2011) for a 2 ${}^{v}$ full factorial with two factor interaction for $v$ factors each at 2 levels and arrange the combinations in lexicographic order. This arrangement contains $v$ columns and 2 ${}^{v}$ runs. Arrange all these runs in first block, and the second block is obtained by rotating the columns of first block in circular fashion. Similarly, rotating the $v$ columns of 2 ${}^{v}$ factorial points $(v-1)$ times we get $v\times 2^{v}$ design points in $v$ blocks each of size 2 ${}^{v}$ and accommodate two factor interactions in the design matrix. The design so obtained satisfies all the conditions mentioned in Section 3.1 under equi-correlated error structure.

We illustrate the method by taking $v=$ 3 (say $x_{1}$ , $x_{2}$ and $x_{3}$ ) with each factor at two levels. The full factorial with two factor interaction contains 8 runs which constitute the first block and the other blocks are obtained by rotating the columns of the first block in a circular fashion. The design matrix ${\bm{X}}$ is obtained by incorporating ${v\choose 2}$ ways of possible interactions among input variables. The ${{\bm{X}}}^{\prime}\bm{V}^{-1}{\bm{X}},$ where $\bm{X}=\left[{\bm{X}}_{1}:{\bm{B}}\right]$ matrix is obtained under equi-correlated error structure for $\rho=0.4$ as

${{\bm{X}}}^{\prime}\bm{V}^{-1}{\bm{X}}=\begin{bmatrix}40&0&0&0&0&0&0&0&0\\ 0&40&0&0&0&0&0&0&0\\ 0&0&40&0&0&0&0&0&0\\ 0&0&0&40&0&0&0&0&0\\ 0&0&0&0&40&0&0&0&0\\ 0&0&0&0&0&40&0&0&0\\ 0&0&0&0&0&0&2.105&0&0\\ 0&0&0&0&0&0&0&2.105&0\\ 0&0&0&0&0&0&0&0&{2.105}\\ \end{bmatrix}$

and $\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right)^{-1}$ is obtained by taking the reciprocal of the diagonal elements.

Thus, ${V(\hat{\beta}}_{i})={0.025\sigma}^{2},{V(\hat{\beta}}_{ij})={0.025\sigma}^{2}$ for $i=1,2,3$ , ${V(\hat{\tau}}_{l})=0.475\sigma^{2}$ for $l=1,2,3,$ and $V\left(\hat{y}(x)\right)={0.308\sigma}^{2}$ and the variance of the estimated responses at all points are the same. The D-optimality of the design obtained from $\det\left({{\bm{X}}_{1}}^{\prime}\bm{V}^{-1}{\bm{X}}_{1}\right)$ is 4096000000 for $\rho=0.4$ . The computational details of D-optimality are given in Section 6.

4. Blocking of response surface models under auto-correlated error structure

The $N\times N$ matrix of auto-correlation structure of order 1 ( $\bm{V}_{{E}}$ ) of errors ‘ $\bm{e}$ ’ in the model Eq. 3 is given by $\bm{V}_{{E}}=\rho^{\left|i-j\right|},1\leqslant i,j\leqslant N$ , then

$\bm{V}_{{E}}^{-1}=\frac{1}{(1-\rho^{2})}\left[\left(1+\rho^{2}\right){\bm{I}}_% {N}-\rho^{2}{\bm{A}}_{N\times N}-\rho{\bm{C}}_{N\times N}\right],$

where $\rho$ is correlation coefficient, ${\bm{I}}_{N}$ is an $N\times N$ identity matrix, ${\bm{A}}$ is an $N\times N$ matrix with elements $a_{11}=$ $a_{NN}=$ 1 and all other elements zeros, and ${\bm{C}}$ is an $N\times N$ matrix with $c_{ij}=$ 1 for $|i-j|=1$ and all other elements are zeros. Then, the generalized least square estimator of ${\bm{\theta}}$ is $\hat{\bm{\theta}}=\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right)^{-1}{\bm{X% }}^{\prime}\bm{V}^{-1}{\bm{Y}}$ , here

${\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}=\frac{1}{(1-\rho^{2})}\left[\left(1+\rho^% {2}\right){\bm{X}}^{\prime}{\bm{I}}_{{N}}{\bm{X}}-\rho^{2}{\bm{X}}^{\prime}{% \bm{A}\bm{X}-\rho}{\bm{X}}^{\prime}{\bm{C}\bm{X}}\right],$ (16)

where ${\bm{I}}_{{\bm{N}}}$ is an $N\times N$ identity matrix, ${\bm{A}}$ is a matrix of order $N\times N$ , which contains ‘ $b$ ’ block diagonal matrices of order $n\times n$ with elements $a_{111}=a_{1nn}=1$ and rest are zeros. Further ${\bm{C}}$ is a matrix of order $N\times N$ , which contains ‘ $b$ ’ block diagonal matrices of order $n\times n$ with $c_{1ij}=1$ for $|i-j|=1$ on the diagonal and all other elements are zeros.

4.1 Conditions for orthogonal blocking and parameter estimation

1.
$\displaystyle v_{i.b}=b\sum\nolimits_{u=1}^{n}x_{iu}+b\rho^{2}\sum\nolimits_{u% =2}^{n-1}x_{iu}-b\rho\left(\sum\nolimits_{u=1}^{n}x_{iu}+\sum\nolimits_{u=2}^{% n-1}x_{iu}\right)=0,1\leqslant i\leqslant v$ .
2.
$\displaystyle v_{ij.b}=b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}}+b\rho^{2}\sum% \nolimits_{u=2}^{n-1}{x_{iu}x_{ju}}-b\rho\left(\sum\nolimits_{u=1}^{n}{x_{iu}x% _{ju}}+\sum\nolimits_{u=2}^{n-1}{x_{iu}x_{ju}}\right)=0$ ,

$1\leqslant i\neq j\leqslant v$ .
3.
$\displaystyle v_{i.j}=b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}}+b\rho^{2}\sum% \nolimits_{u=2}^{n-1}{x_{iu}x_{ju}}-b\rho\left(\sum\nolimits_{u=1}^{n-1}{x_{iu% }x_{j(u+1)}}+\sum\nolimits_{u=1}^{n-1}{x_{i(u+1)}x_{ju}}\right)=0$ ,

$1\leqslant i\neq j\leqslant v$ .
4.
$\displaystyle v_{ii.j}=b\sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}}+b\rho^{2}% \sum\nolimits_{u=2}^{n-1}{x_{iu}^{2}x_{ju}}-b\rho\left({\begin{array}[]{l}% \displaystyle\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i(u+1)}x_{ju}}+\\ \displaystyle\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i(u+1)}x_{j(u+1)}}\\ \end{array}}\right)=0,1\leqslant i\neq j\leqslant v$ .
5.
$\displaystyle v_{ij.s}=b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}x_{su}}+b\rho^{2}% \sum\nolimits_{u=2}^{n-1}{x_{iu}x_{ju}x_{su}}-b\rho\left({\begin{array}[]{l}% \displaystyle\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{j\left(u+1\right)}x_{s\left(u+% 1\right)}}+\\ \displaystyle\sum\nolimits_{u=1}^{n-1}{x_{i\left(u+1\right)}x_{ju}x_{su}}\\ \end{array}}\right)=0$ ,

$1\leqslant i<j,s\leqslant v$ .
6.
$\displaystyle v_{ii.js}=b\sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}x_{su}}+b\rho% ^{2}\sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}x_{su}}-b\rho\left({\begin{array}[% ]{l}\displaystyle\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u+1\right)}x_{j% \left(u+1\right)}x_{su}}+\\ \displaystyle\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u+1\right)}x_{ju}x_{s% \left(u+1\right)}}\\ \end{array}}\right)=0$ ,

$1\leqslant i,j<s\leqslant v$ .
7.
$\displaystyle v_{ij.st}=b\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}x_{su}x_{tu}}+b% \rho^{2}\sum\nolimits_{u=1}^{n}{x_{iu}x_{ju}x_{su}x_{tu}}-b\rho\left({\begin{% array}[]{l}\displaystyle\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{j(u+1)}x_{s(u+1)}x_% {tu}}+\\ \displaystyle\sum\nolimits_{u=1}^{n-1}{x_{i(u+1)}x_{ju}x_{su}x_{t(u+1)}}\\ \end{array}}\right)=0$ ,

$1\leqslant i,s<j,t\leqslant v$ .

These restrictions ensure the rotatability of the experimental design ${\bm{X}}_{1}$ under orthogonal blocks.

Now, ${\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}$ can be written as

$\displaystyle{\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}=\frac{1}{\left(1-\rho^{2}% \right)}\text{diag}$ $\displaystyle\left\{\begin{array}[]{@{}l@{}}\displaystyle b\sum\nolimits_{u=1}% ^{n}x_{1u}^{2}+b\rho^{2}\sum\nolimits_{u=2}^{n-1}x_{1u}^{2}-2b\rho\sum% \nolimits_{u=1}^{n-1}{x_{1u}x_{1\left(u+1\right)}},\ldots,\\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{iu}^{2}+{b\rho}^{2}\sum\nolimits_{u=2% }^{n-1}x_{iu}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u+1\right)}}% ,\ldots,\\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{vu}^{2}+b\rho^{2}\sum\nolimits_{u=2}^% {n-1}x_{vu}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{vu}x_{v\left(u+1\right)}},% \\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{1u}^{2}x_{2u}^{2}+b\rho^{2}\sum% \nolimits_{u=2}^{n-1}x_{1u}^{2}x_{2u}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{1% u}x_{1\left(u+1\right)}x_{2u}x_{2\left(u+1\right)}},\ldots,\\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{iu}^{2}x_{ju}^{2}+b\rho^{2}\sum% \nolimits_{u=2}^{n-1}x_{iu}^{2}x_{ju}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{% iu}x_{i\left(u+1\right)}x_{ju}x_{j\left(u+1\right)}},\ldots,\\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{\left(v-1\right)u}^{2}x_{vu}^{2}+b% \rho^{2}\sum\nolimits_{u=2}^{n-1}x_{\left(v-1\right)u}^{2}x_{vu}^{2}-2b\rho% \sum\nolimits_{u=1}^{n-1}{x_{\left(v-1\right)u}x_{\left(v-1\right)\left(u+1% \right)}x_{vu}x_{v\left(u+1\right)}},\\ \displaystyle\left(1-\rho\right)\left(n-\left(n-2\right)\rho\right),\ldots,% \left(1-\rho\right)\left(n-\left(n-2\right)\rho\right)\\ \end{array}\right\}.$ (17)

Then, the parameters ( ${\bm{\theta}}$ ’s) are estimated from the normal equation $\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right){\bm{\theta}=}{{\bm{X}}}^{% \prime}\bm{V}^{-1}{\bm{Y}}$

$\text{i.e.,∼{}}\begin{bmatrix}{\bm{S}}_{1}&0^{\prime}\\ 0&{\bm{S}}_{2}\\ \end{bmatrix}\begin{bmatrix}{\bm{\beta}}\\ {\bm{\tau}}\\ \end{bmatrix}=\frac{1}{1-\rho}\begin{bmatrix}{\bm{T}}\\ {\bm{B}}\\ \end{bmatrix},$

where ${\bm{T}}={[T_{1},\ldots,T_{i},\ldots,T_{v},T_{12},\ldots,T_{ij},\ldots,T_{% \left(v-1\right)v}]}^{\prime}$ and ${\bm{B}}={[B_{1},\ldots,B_{l},\ldots,B_{b},]}^{\prime}$ are the vector of treatment combination totals and block totals, respectively;

$\displaystyle T_{i}=\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=1}^{n}{x_{iu}y_{lu% }}+\rho^{2}\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=2}^{n-1}{x_{iu}y_{lu}}$ $\displaystyle\qquad-\rho\left(\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=2}^{n}{x% _{iu}y_{l(u-1)}}+\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=1}^{n-1}{x_{iu}y_{l(u% +1)}}\right);$ $\displaystyle T_{ij}=\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=1}^{n}{x_{iu}x_{% ju}y_{lu}}+\rho^{2}\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=2}^{n-1}{x_{iu}x_{% ju}y_{lu}}-\rho\left({\begin{array}[]{l}\displaystyle\sum\nolimits_{l=1}^{b}% \sum\nolimits_{u=2}^{n}{x_{iu}x_{ju}y_{l\left(u-1\right)}}+\\ \displaystyle\sum\nolimits_{l=1}^{b}\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{ju}y_{l% \left(u+1\right)}}\\ \end{array}}\right)$ $\displaystyle\text{and∼{}}B_{l}=\sum\nolimits_{u=1}^{n}y_{lu}+\rho^{2}\sum% \nolimits_{u=2}^{n-1}y_{lu}-\rho\left(\sum\nolimits_{u=1}^{n}y_{lu}+\sum% \nolimits_{u=2}^{n-1}y_{lu}\right)$ $\displaystyle{\bm{S}}_{1}=\frac{1}{\left(1-\rho^{2}\right)}\text{diag}$ $\displaystyle\left\{\begin{array}[]{@{}l@{}}\displaystyle b\sum\nolimits_{u=1}% ^{n}x_{1u}^{2}+b\rho^{2}\sum\nolimits_{u=2}^{n-1}x_{1u}^{2}-2b\rho\sum% \nolimits_{u=1}^{n-1}{x_{1u}x_{1\left(u+1\right)}},\ldots,\\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{vu}^{2}+b\rho^{2}\sum\nolimits_{u=2}^% {n-1}x_{vu}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{vu}x_{v\left(u+1\right)}},% \\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{1u}^{2}x_{2u}^{2}+b\rho^{2}\sum% \nolimits_{u=2}^{n-1}x_{1u}^{2}x_{2u}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{1% u}x_{1\left(u+1\right)}x_{2u}x_{2\left(u+1\right)}},\ldots,\\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{iu}^{2}x_{ju}^{2}+b\rho^{2}\sum% \nolimits_{u=2}^{n-1}x_{iu}^{2}x_{ju}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{% iu}x_{i\left(u+1\right)}x_{ju}x_{j\left(u+1\right)}},\ldots,\\ \displaystyle b\sum\nolimits_{u=1}^{n}x_{\left(v-1\right)u}^{2}x_{vu}^{2}+b% \rho^{2}\sum\nolimits_{u=2}^{n-1}x_{\left(v-1\right)u}^{2}x_{vu}^{2}-2b\rho% \sum\nolimits_{u=1}^{n-1}{x_{\left(v-1\right)u}x_{\left(v-1\right)\left(u+1% \right)}x_{vu}x_{v\left(u+1\right)}},\\ \end{array}\right\}$ $\displaystyle{\bm{S}}_{2}=\frac{1}{\left(1-\rho^{2}\right)}\text{diag}\left[% \left(1-\rho\right)\left(n-\left(n-2\right)\rho\right),\ldots,\left(1-\rho% \right)\left(n-\left(n-2\right)\rho\right)\right].$

Hence,

$\hat{\bm{\theta}}=\begin{bmatrix}\hat{\bm{\beta}}\\ \hat{\bm{\tau}}\\ \end{bmatrix}=\begin{bmatrix}{\bm{S}}_{1}^{-1}{\bm{T}}\\ {\bm{S}}_{2}^{-1}{\bm{B}}\\ \end{bmatrix}\text{ and }{\bm{D}}\begin{bmatrix}\hat{\bm{\beta}}\\ \hat{\bm{\tau}}\\ \end{bmatrix}=\sigma^{2}\begin{bmatrix}{\bm{S}}_{1}^{-1}&0^{\prime}\\ 0&{\bm{S}}_{2}^{-1}\\ \end{bmatrix}.$ (18)

The variance of the estimated parameters is obtained as

$\displaystyle V(\hat{\beta}_{i})={\sigma}^{2}\left(1-\rho^{2}\right)\frac{1}{% \displaystyle b\sum\nolimits_{u=1}^{n}x_{iu}^{2}+{b\rho}^{2}\sum\nolimits_{u=2% }^{n-1}x_{iu}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u+1\right)}}},$ (19) $\displaystyle V(\hat{\beta}_{ij})={\sigma}^{2}\left(1-\rho^{2}\right)\frac{1}{% \displaystyle b\sum\nolimits_{u=1}^{n}x_{iu}^{2}x_{ju}^{2}+b\rho^{2}\sum% \nolimits_{u=2}^{n-1}x_{iu}^{2}x_{ju}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{% iu}x_{i\left(u+1\right)}x_{ju}x_{j\left(u+1\right)}}},$ (20)

and

$\displaystyle V(\hat{\tau}_{l})=\sigma^{2}\left(1-\rho^{2}\right)\frac{1}{% \displaystyle\left(1-\rho\right)\left(n-\left(n-2\right)\rho\right)}.$ (21)

The estimated response at the point ${{\bm{g}(\bm{x})}}^{\prime}$ and its prediction variance via Eq. 6 is given by

$\displaystyle V\left(\hat{y}(x)\right)=\left(1-\rho^{2}\right)\sigma^{2}$ $\displaystyle\begin{Bmatrix}\displaystyle\frac{x_{10}^{2}}{\displaystyle b\sum% \nolimits_{u=1}^{n}x_{1u}^{2}+b\rho^{2}\sum\nolimits_{u=2}^{n-1}x_{1u}^{2}-2b% \rho\sum\nolimits_{u=1}^{n-1}{x_{1u}x_{1\left(u+1\right)}}}+\ldots+\\ \displaystyle\frac{x_{i0}^{2}}{\displaystyle b\sum\nolimits_{u=1}^{n}x_{iu}^{2% }+{b\rho}^{2}\sum\nolimits_{u=2}^{n-1}x_{iu}^{2}-2b\rho\sum\nolimits_{u=1}^{n-% 1}{x_{iu}x_{i\left(u+1\right)}}}+\ldots+\\ \displaystyle\frac{x_{v0}^{2}}{\displaystyle b\sum\nolimits_{u=1}^{n}x_{vu}^{2% }+b\rho^{2}\sum\nolimits_{u=2}^{n-1}x_{vu}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}% {x_{vu}x_{v\left(u+1\right)}}}+\\ \displaystyle\frac{x_{10}^{2}x_{20}^{2}}{\displaystyle b\sum\nolimits_{u=1}^{n% }x_{1u}^{2}x_{2u}^{2}+b\rho^{2}\sum\nolimits_{u=2}^{n-1}x_{1u}^{2}x_{2u}^{2}-2% b\rho\sum\nolimits_{u=1}^{n-1}{x_{1u}x_{1\left(u+1\right)}x_{2u}x_{2\left(u+1% \right)}}}+\ldots+\\ \displaystyle\frac{x_{i0}^{2}x_{j0}^{2}}{\displaystyle b\sum\nolimits_{u=1}^{n% }x_{iu}^{2}x_{ju}^{2}+b\rho^{2}\sum\nolimits_{u=2}^{n-1}x_{iu}^{2}x_{ju}^{2}-2% b\rho\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u+1\right)}x_{ju}x_{j\left(u+1% \right)}}}+\ldots+\\ \displaystyle\frac{x_{\left(v-1\right)0}^{2}x_{v0}^{2}}{\displaystyle b\sum% \nolimits_{u=1}^{n}x_{\left(v-1\right)u}^{2}x_{vu}^{2}+b\rho^{2}\sum\nolimits_% {u=2}^{n-1}x_{\left(v-1\right)u}^{2}x_{vu}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}% {x_{\left(v-1\right)u}x_{\left(v-1\right)\left(u+1\right)}x_{vu}x_{v\left(u+1% \right)}}}+\\ \displaystyle\frac{1}{b\left(1-\rho\right)\left(n-\left(n-2\right)\rho\right)}% \\ \end{Bmatrix}.$ (22)

The constancy of the variance of estimated parameters is ensured by putting the following conditions

1.
$\displaystyle\sum\nolimits_{u=1}^{n}x_{iu}^{2}=g_{1}$ and $\displaystyle\sum\nolimits_{u=1}^{n}x_{iu}^{2}x_{ju}^{2}=h_{1}$ ,
2.
$\displaystyle\sum\nolimits_{u=2}^{n-1}x_{iu}^{2}=g_{2}$ and $\displaystyle\sum\nolimits_{u=2}^{n-1}{x_{iu}^{2}x_{ju}^{2}}=h_{2}$ ,
3.
$\displaystyle\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u+1\right)}}=g_{3}$ and $\displaystyle\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u+1\right)}}x_{ju}x_{j% \left(u+1\right)}=h_{3}$ ,

where $g_{1},g_{2}$ and $g_{3}$ ; $h_{1},h_{2}$ and $h_{3}$ are constants for $\forall i\neq j=1,2,\ldots,v$ for all blocks of size ‘ $n$ ’. Therefore,

$V\left(\hat{y}(x)\right)=\left(1-\rho^{2}\right)\sigma^{2}\begin{Bmatrix}% \displaystyle\frac{\displaystyle\sum\nolimits_{i=1}^{v}x_{i0}^{2}}{b\left(g_{1% }+\rho^{2}g_{2}-2\rho g_{3}\right)}+\frac{\displaystyle\sum\nolimits_{i\neq j=% 1}^{v}{x_{i0}^{2}x_{j0}^{2}}}{b\left(h_{1}+{\rho}^{2}h_{2}-2\rho h_{3}\right)}% \\ \displaystyle+\frac{1}{b\left(1-\rho\right)\left(n-\left(n-2\right)\rho\right)% }\\ \end{Bmatrix}.$ (23)

This illustrates that the variances of estimators of linear regression coefficients $\hat{\beta}_{i}$ ’s and interaction regression coefficients $\hat{\beta}_{ij}$ ’s $(i\neq j=1,2,\ldots,v)$ are same and the variance of the estimated response is a function of $\sum\nolimits_{i=1}^{v}x_{i0}^{2}$ and $\sum\nolimits_{i\neq j=1}^{v}{x_{i0}^{2}x_{j0}^{2}}$ for $\hat{\beta}_{i}$ ’s and $\hat{\beta}_{ij}$ ’s, respectively. For a given $\rho$ , the estimated response will have the same variance under auto-correlated error structure.
4.2 Method of construction and illustration

We considered the method of construction given in the Section 3.2. This set of $v\times 2^{v}$ design points is extended to $\{v[(2\times 2^{v})+1]\}$ points by adding $(2^{v}+1)$ zero points in each block. One zero point is added before and after the design points and one zero point is added in between each set of design points in the sequence of each block. The designs so obtained has all the factors each at three levels in $\{v[(2\times 2^{v})+1]\}$ runs and satisfies all the conditions mentioned in Section 4.1. Thus, the obtained design is rotatable under auto-correlated error structure.

We illustrate the method by taking $v=3$ (say $x_{1}$ , $x_{2}$ and $x_{3}$ ) with each factor at two levels, the full factorial with two factor interaction contain 8 runs. Arrange these runs in v blocks as mentioned above and add one zero point before and after the design points and one zero point is added in between each set of design points in the sequence. Thus each block contains 17 runs, where each column is a run and the different block matrix is obtained as

Block 1	$x_{1}$	1	1	1	1	$-$ 1	$-$ 1	$-$ 1	$-$ 1
	$x_{2}$	1	1	$-$ 1	$-$ 1	1	1	$-$ 1	$-$ 1
	$x_{3}$	1	$-$ 1	1	$-$ 1	1	$-$ 1	1	$-$ 1
Block 2	$x_{1}$	1	1	$-$ 1	$-$ 1	1	1	$-$ 1	$-$ 1
	$x_{2}$	1	$-$ 1	1	$-$ 1	1	$-$ 1	1	$-$ 1
	$x_{3}$	1	1	1	1	$-$ 1	$-$ 1	$-$ 1	$-$ 1
Block 3	$x_{1}$	1	$-$ 1	1	$-$ 1	1	$-$ 1	1	$-$ 1
	$x_{2}$	1	1	1	1	$-$ 1	$-$ 1	$-$ 1	$-$ 1
	$x_{3}$	1	1	$-$ 1	$-$ 1	1	1	$-$ 1	$-$ 1

The ${{\bm{X}}}^{\prime}\bm{V}^{-1}\bm{X}$ , where $\bm{X}=\left[{\bm{X}}_{1}:{\bm{B}}\right]$ matrix is obtained under auto correlated error structure for $\rho=0.4$ as

${{\bm{X}}}^{\prime}\bm{V}^{-1}{\bm{X}}=\begin{bmatrix}{33.143}&0&0&0&0&0&0&0&0% \\ 0&{33.143}&0&0&0&0&0&0&0\\ 0&0&{33.143}&0&0&0&0&0&0\\ 0&0&0&{33.143}&0&0&0&0&0\\ 0&0&0&0&{33.143}&0&0&0&0\\ 0&0&0&0&0&{33.143}&0&0&0\\ 0&0&0&0&0&0&{7.86}&0&0\\ 0&0&0&0&0&0&0&{7.86}&0\\ 0&0&0&0&0&0&0&0&{7.86}\\ \end{bmatrix}$

and $\left({\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right)^{-1}$ is obtained by taking the reciprocal of the diagonal elements.

Thus, ${V(\hat{\beta}}_{i})={0.0302\sigma}^{2},{V(\hat{\beta}}_{ij})={0.0302\sigma}^{2}$ for $i=1,2,3,{V(\hat{\tau}}_{l})=0.127\sigma^{2}$ for $l=1,2,3$ and $V(\hat{y}(x))$ $=$ $0.221\sigma^{2}$ the variance of the estimated responses at all points are the same. The D-optimality of the design obtained from $\det\left({\bm{X}}_{1}^{\prime}\bm{V}^{-1}{\bm{X}}_{1}\right)$ is 1336209499 for $\rho=0.4$ . The computational details of D-optimality are given in Section 6.

5. Effect of blocking on the prediction variance under correlated error

As described by Khuri (1994), the prediction variance, $V\left(\hat{y}(x)\right)$ at a point ${{\bm{g}(\bm{x})}}^{\prime}$ for the interaction model under equi- and auto-correlated error can be represented respectively as

$V\left(\hat{y}(x)\right)=V(\hat{y}_{0}(x))+(1-\rho)\sigma^{2}\left\{\frac{1+% \left(n-1\right)\rho}{bn\left(1-\rho\right)}-\frac{1+\left(N-1\right)\rho}{N% \left(1-\rho\right)}\right\}$

and

$V\left(\hat{y}(x)\right)=V(\hat{y}_{0}(x))+(1-\rho^{2})\sigma^{2}\left\{\frac{% 1}{b\left(1-\rho\right)\left(n-\left(n-2\right)\rho\right)}-\frac{1}{\left(1-% \rho\right)\left(N-\left(N-2\right)\rho\right)}\right\},$

where $V(\hat{y}_{0}(x))$ is the prediction variance obtained for interaction model at the point $\bm{g}(\bm{x})^{\prime}$ under correlated error structure when block effects are zero, that is $\bm{\delta}=0$ .

Thus, $V\left(\hat{y}(x)\right)\geqslant V(\hat{y}_{0}(x))$ since $V(\hat{y}_{0}(x))\geqslant 0$ . Thus, if the design blocks are orthogonal under correlated error structure, then the prediction variance is influenced by block size ‘ $n$ ’ and ‘ $\rho$ ’.

6. D-optimality under correlated error

Goos and Donev (2006) derived D-optimality criteria for an extended design matrix ${\bm{X}}_{1}$ and a given block structure ${\bm{B}}$ when the block effect is fixed and additive by maximizing $D=\left|{\bm{X}}^{\prime}{\bm{X}}\right|$ , where ${\bm{X}}=\left[{\bm{X}}_{1}:{\bm{B}}\right]$ . In this section, D-optimality of orthogonally blocked first order response surface designs with interaction under equi- and auto-correlated error is discussed. A D-optimal orthogonal blocked design under correlated error is obtained by maximizing $D=\left|{\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}\right|$ , where ${\bm{X}}=\left[{\bm{X}}_{1}:{\bm{B}}\right]$ , as described above.

The information matrix for a given design ${\bm{X}}_{1}$ and a given block structure $\bm{B}$ of a fixed block effect model under correlated error structure can be written as

$\sigma^{-2}{\bm{X}}^{\prime}\bm{V}^{-1}{\bm{X}}={\sigma}^{-2}\begin{bmatrix}{% \bm{X}}_{1}^{\prime}\bm{V}^{-1}{\bm{X}}_{1}&{\bm{X}}_{1}^{\prime}\bm{V}^{-1}{% \bm{B}}\\ {\bm{B}}^{\prime}\bm{V}^{-1}{\bm{X}}_{1}&{\bm{B}}^{\prime}\bm{V}^{-1}{\bm{B}}% \\ \end{bmatrix}$

and its determinant is given by

$\det{({\bm{B}}^{\prime}\bm{V}^{-1}{\bm{B}})}\det\left({\bm{X}}_{1}^{\prime}\bm% {V}^{-1}{\bm{X}}_{1}-{{\bm{X}}}_{1}^{\prime}\bm{V}^{-1}{\bm{B}}\left({\bm{B}}^% {\prime}\bm{V}^{-1}{\bm{B}}\right)^{-1}{\bm{B}}^{\prime}\bm{V}^{-1}{\bm{X}}_{1% }\right)\sigma^{-2}.$

The first factor $\det({\bm{B}}^{\prime}\bm{V}^{-1}{\bm{B}})$ is a constant because the blocking arrangement is predefined and the D-optimal design is obtained by maximizing the second factor. To maximize the above determinant, we impose restrictions described in Sections 3.1 and 4.1 for equi- and auto-correlated errors respectively.

Hence,

$\det\left({\bm{X}}_{1}^{\prime}\bm{V}^{-1}{\bm{X}}_{1}-{\bm{X}}_{1}^{\prime}% \bm{V}^{-1}{\bm{B}}\left({\bm{B}}^{\prime}\bm{V}^{-1}{\bm{B}}\right)^{-1}{\bm{% B}}^{\prime}\bm{V}^{-1}{\bm{X}}_{1}\right)=\det{({\bm{X}}_{1}^{\prime}\bm{V}^{% -1}{\bm{X}}_{1}).}$

Thus, if $\det{({\bm{X}}_{1}^{\prime}\bm{V}^{-1}{\bm{X}}_{1})}$ is maximum for a chosen design ${\bm{X}}_{1}$ , then the design is D-optimal first order response surface design with interaction under correlated error structure.

The D-optimality of equi-correlated error structure is obtained from

$D=\frac{1}{(1-\rho)}\det\left\{\text{diag}(D_{1},\ldots,D_{i},\ldots,D_{v},D_{% 12},\ldots,D_{ij},\ldots,D_{\left(v-1\right)v})\right\},$

where,

$D_{i}=b\sum\nolimits_{u=1}^{n}x_{iu}^{2},\forall i=1,2,\ldots,v\text{∼{}and∼{}% }D_{ij}=b\sum\nolimits_{u=1}^{n}{x_{iu}^{2}x_{ju}^{2}},\forall i\neq j=1,2,% \ldots,v.$

The D-optimality of auto-correlated error structure is obtained from

$D=\frac{1}{(1-\rho^{2})}\mathrm{det}\left\{\text{diag}(D_{1},\ldots,D_{i},% \ldots D_{v},D_{12},\ldots,D_{ij},\ldots,D_{\left(v-1\right)v})\right\},$

where,

$D_{i}=b\sum\nolimits_{u=1}^{n}x_{iu}^{2}+{b\rho}^{2}\sum\nolimits_{u=2}^{n-1}x% _{iu}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u+1\right)}},\forall i% =1,2,\ldots v\text{∼{}and}$

$D_{ij}=b\sum\nolimits_{u=1}^{n}x_{iu}^{2}x_{ju}^{2}+b\rho^{2}\sum\nolimits_{u=% 2}^{n-1}x_{iu}^{2}x_{ju}^{2}-2b\rho\sum\nolimits_{u=1}^{n-1}{x_{iu}x_{i\left(u% +1\right)}x_{ju}x_{j\left(u+1\right)}},\forall i\neq j=1,2,\ldots v.$

Specific examples are discussed in Sections 3.2 and 4.2 respectively for equi- and auto-correlated errors.

Table 1
Experimental lay out and response

	Block/Group 1				Block/Group 2				Block/Group 3
	CAP	CAP	CAP	CAP	CVP	CVP	CVP	CVP	MAP	MAP	MAP	MAP
$X_{1}$	1	1	$-$ 1	$-$ 1	1	1	$-$ 1	$-$ 1	1	1	$-$ 1	$-$ 1
$X_{2}$	$-$ 1	1	$-$ 1	1	$-$ 1	1	$-$ 1	1	$-$ 1	1	$-$ 1	1
$Y$	4.25	21.00	2.29	1.80	2.90	14.00	1.89	1.48	0.25	9.24	0.03	0.58

$X_{1}=$ storage day, $X_{2}=$ storage temperature and $Y=$ histamine.

Table 2

MSE as a function of $\rho$

${\bm{\rho}}$	$-$ 0.9	$-$ 0.8	$-$ 0.7	$-$ 0.6	$-$ 0.5	$-$ 0.4	$-$ 0.3	$-$ 0.2	$-$ 0.1	0
MSE	6.10	6.08	6.05	6.03	6.01	5.99	5.97	5.95	5.94	5.94
${\bm{\rho}}$	0.9	0.8	0.7	0.6	0.5	0.4	0.3	0.2	0.1	0
MSE	7.05	6.68	6.42	6.25	6.13	6.04	5.99	5.96	5.94	5.94

7. Analysis of simulated data with real-time application

Consider an experimental condition where the Histamine formulation in milkfish (Chanos chanos) during storage under different packing conditions and temperatures. Histamine is one of the important quality indices for ascertaining the quality of the fish. The treatment combinations are obtained for different levels of storage day and temperature in 2 ${}^{2}$ factorial setup. The milkfish steaks were stored under air (CAP), vacuum (CVP) and modified atmosphere (MAP 60% CO ${}_{2}$ and 40% O ${}_{2}$ ) packing conditions at 0 ${}^{\circ}$ C and 8 ${}^{\circ}$ C. The fish sample kept at 0 ${}^{\circ}$ C and 8 ${}^{\circ}$ C in different packing conditions is drawn in different time interval (3 ${}^{\rm rd}$ and 15 ${}^{\rm th}$ day) and measured Histamine content. The experimental set up with hypothetical values of response variable is given in Table 1. The packing conditions are considered as blocks/groups as different blocks/groups are significantly affected the response.

Let us assume that ${\bm{V}}$ is the variance-covariance matrix of ${\bm{e}}$ with $\rho=0.4$ as the correlation between successive elements. The fitted model under the assumed error structure for $\rho=0.4$ is

$\begin{array}[]{l}\hat{{Y}}=4.97+3.63X_{1}+3.04X_{2}+3.10X_{1}X_{2}\text{ with% ∼{}MSE}=6.04\\ \qquad(0.7)∼{}∼{}∼{}∼{}(0.7)∼{}∼{}∼{}∼{}∼{}∼{}∼{}∼{}∼{}(0.7)∼{}∼{}∼{}∼{}∼{}∼{}% ∼{}∼{}∼{}(0.7)\end{array}$

The $R^{2}$ value of the fitted model was 0.92 and all the regression coefficients was found significant at 5% level of significance ( $p<$ 0.05). The standard error of the estimated regression coefficients are given in the parenthesis.

A SAS programme was developed to see the effect of different values of $\rho$ on the mean square error (MSE) of the estimated response. The MSE values for different values of $\rho$ are given in Table 2. and it is observed that when $|\rho|>0.4$ , the MSE value slightly increases with $|\rho|$ . The MSE values at $\rho=0.4$ is compared with $\rho=0$ using F-test and found no significant difference at 5% level of significance.

8. Conclusion

The effect of orthogonal blocking in first order response surface model with interaction was studied under equi-and auto-correlated error structure. The generalized least square estimates of linear ( $\beta_{i}$ ’s) and interaction ( $\beta_{ij}$ ’s) regression coefficients under equi-and auto-correlated error structure were not affected by the orthogonal blocking, but the mean response was affected by block size ‘ $n$ ’ and ‘ $\rho$ ’. The developed model found to be rotatable as the prediction variance was found to be constant at all points that are equi-distant from the design center. The D-optimality criterion was derived for the construction of orthogonally blocked first order response surface models with interaction under equi-and auto-correlated error structure.

Footnotes

Acknowledgments

The authors are grateful to Dr. K.V. Lalitha, Principal Scientist, ICAR-Central Institute of Fisheries Technology, Cochin for providing real-time data for generating hypothetical data. The authors also thank the referees for their constructive suggestions on the earlier draft of the paper.

References

Box

G. E. P.

, & Draper

N. R.

(1987). Empirical model building and response surfaces. John wiley and sons, New York.

Das

R. N.

(2003a). Slope-rotatabilty with correlated errors. Calcutta Statistical Association Bulletin, 54, 57-70.

Das

R. N.

(2003). Robust second order rotatable designs: Part III (RSORD). Journal of Indian Society of Agricultural Statistics, 56(2), 117-130.

Das

R. N.

, & Park

S. H.

(2009). A measure of robust slope-rotatability for second order response surface experimental designs. Journal of Applied Statistics, 36, 755-767.

Goos

, & Doonev

A. N.

(2006). Blocking response surface designs. Computational Statistics & Data Analysis, 51, 1075-1088.

Khuri

A. I.

, & Cornel

J. A.

(1987). Response surfaces, designs and analysis, Marcel Dekker in New York.

Khuri

A. I.

(1992). Response surface models with random block effects. Technometrics, 34, 26-37.

Khuri

A. I.

(1994). Effect of blocking on the estimation of a response surface. Journal of Applied Statistics, 21, 305-3016.

Mann

R. K.

Edwards

, & Zhou

(2010). Robust designs for experiments with blocks. Communication in Statistics – Theory and Methods, 45(18), 35-45.

10.

Myers

R. H.

, & Mongomery

D. C.

(2002). Response surface methodology – Process and product optimization using designed experiments. John wiley and sons, New York.

11.

Panda

R. N.

, & Das

R. N.

(1994). First order rotatable designs with correlated errors (FORDEWCE). Calcutta Statistical Association Bulletin, 44, 173-174.

12.

Rajyalakshmi

, & Victorbabu

B. R.

(2015). Second order slope rotatable designs under intra-class correlated error structure using balanced incomplete block designs. Thailand Statisticain, 13(2), 169-183.

13.

Varghese

, & Jaggi

(2011). Blocking response surface designs incorporating neighbor effects. Open Journal of Statistics, 1, 199-204.

14.

Varghese

, & Jaggi

(2013). Response surface model with neighbour effects and correlated observations, Model Assisted Statistics and Applications, 8(1), 41-49.

Blocking first order response surface designs with interaction under correlated error

Abstract

Keywords

1. Introduction

2. Response surface model with blocks under correlated error

4. Blocking of response surface models under auto-correlated error structure

5. Effect of blocking on the prediction variance under correlated error

6. D-optimality under correlated error

Table 1 Experimental lay out and response

8. Conclusion

Footnotes

Acknowledgments

References

Table 1
Experimental lay out and response