Response surface designs using the generalized variance inflation factors

Abstract

We study response surface designs using the generalized variance inflation factors for subsets as an extension of the variance inflation factors.

Keywords:

• variance inflation factors
• generalized variance inflation factors
• response surface designs

Public Interest Statement

Response surface designs are a mainstay in applied statistics. The variance inflation factors VIF are a measure of collinearity for a single variable in a linear regression model. The generalization to subsets of variables is the generalized variance inflation factor GVIF. This research introduces $GVIF$ as a penalty measure for extending a linear response model to a response surface with the included quadratic terms. The methodology is demonstrated with case studies, and, in particular, it is shown that using GVIF, the H310 design can be improved for the standard global optimality criteria of $A$, $D$, and $E$.

1. Introduction

We consider a linear regression $\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon }$ with $\mathbf{X}$ a full rank $n\times p$ matrix and $\mathcal{L}(\mathit{\epsilon })=N(\mathbf{0},{\mathit{\sigma }}^2{\mathbf{I}}_n)$. The variance inflation factor $VIF$, Belsley (1986), measures the penalty for adding one non-orthogonal additional explanatory variable to a linear regression model, and they can be computed as a ratio of determinants. The extension of $VIF$ to a measure of the penalty for adding a subset of variables to a model is the generalized variance inflation factor $GVIF$ of Fox and Monette (1992), which will be used to study response surface designs, in particular, as the penalty for adding the quadratic terms to the model.

2. Variance inflation factors

For our linear model $\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon }$, let ${\mathbf{D}}_X$ be the diagonal matrix with entries on the diagonal ${\mathbf{D}}_X[i,i]={({\mathbf{X}}^{\prime }\mathbf{X})}_{i,i}^{-1/2}$. When the design has been standardized $\mathbf{X}\to \mathbf{X}{D}_X$ , the $VIF$s are the diagonal entries of the inverse of ${\mathbf{S}}_X={\mathbf{D}}_X(X^{\prime }\mathbf{X})D_X$. That is, the $VIF$s are the ratios of the actual variances for the explanatory variables to the “ideal” variances had the columns of $\mathbf{X}$ been orthogonal. Note that we follow Stewart (1987) and do not necessarily center the explanatory variables.

For our linear model $\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon }$, view $\mathbf{X}=[{\mathbf{X}}_{[p]},{\mathbf{x}}_p]$ with ${\mathbf{x}}_p$ the $p^{\mathrm{th}}$ column of $\mathbf{X}$ and ${\mathbf{X}}_{[p]}$ the matrix formed by the remaining columns. The variance inflation factor $VI{F}_p$ measures the effect of adding column ${\mathbf{x}}_p$ to ${\mathbf{X}}_{[p]}$. For notational convenience, we demonstrate $VI{F}_p$ with the last column $p$. An ideal column would be orthogonal to the previous columns with the entries in the off-diagonal elements of the $p^{\mathrm{th}}$ row and $p^{\mathrm{th}}$ column of ${\mathbf{X}}^{\prime }\mathbf{X}$ all zeros. Denote by ${\mathbf{M}}_p$ the idealized moment matrix$$\begin{array}{c}\hfill {\mathbf{M}}_p=\left[\begin{array}{cc}{\mathbf{X}}_{[p]}^{\prime }{\mathbf{X}}_{[p]}& {\mathbf{0}}_{p-1}\\ {\mathbf{0}}_{p-1}^{\prime }& {\mathbf{x}}_p^{\prime }{\mathbf{x}}_p\end{array}\right].\end{array}$$

The $VIF$s are the diagonal entries of ${\mathbf{S}}_X^{-1}={\mathbf{D}}_X^{-1}{({\mathbf{X}}^{\prime }\mathbf{X})}^{-1}{\mathbf{D}}_X^{-1}$. It remains to note that the inverse, ${\mathbf{S}}_X^{-1}$, can be computed using cofactors ${\mathbf{C}}_{i,j}$. In particular,(1) $$\begin{array}{cc}\hfill VI{F}_p& ={[{\mathbf{S}}_X^{-1}]}_{p,p}={[{\mathbf{D}}_X^{-1}{({\mathbf{X}}^{\prime }\mathbf{X})}^{-1}{\mathbf{D}}_X^{-1}]}_{p,p}\hfill \\ \hfill & ={({\mathbf{x}}_p^{\prime }{\mathbf{x}}_p)}^{1/2}\frac{det({\mathbf{C}}_{p,p})}{det({\mathbf{X}}^{\prime }\mathbf{X})}{({\mathbf{x}}_p^{\prime }{\mathbf{x}}_p)}^{1/2}=\frac{det({\mathbf{M}}_p)}{det({\mathbf{X}}^{\prime }\mathbf{X})}\hfill \end{array}$$(1)

the ratio of the determinant of the idealized moment matrix ${\mathbf{M}}_p$ to the determinant of the moment matrix ${\mathbf{X}}^{\prime }\mathbf{X}$. This definition extends naturally to subsets and is discussed in the next section.

For an alternate view of the how collinearities in the explanatory variables inflate the model variances of the regression coefficients when compared to a fictitious orthogonal reference design, consider the formula for the model variance$$\begin{array}{c}\hfill Va{r}_M({\widehat{\mathit{\beta }}}_j)=\frac{{\mathit{\sigma }}^2}{{\sum }_{i=1}^n{(x_{ij}-{\overline{x}}_j)}^2}\frac{1}{1-R_j^2}\end{array}$$

where $R_j^2$ is the square of the multiple correlation from the regression of the $j^{\mathrm{th}}$ column of $\mathbf{X}=[x_{ij}]$ on the remaining columns as in Liao and Valliant (2012). The first term ${\mathit{\sigma }}^2/\sum {(x_{ij}-{\overline{x}}_j)}^2$ is the model variance for ${\widehat{\mathit{\beta }}}_j$ had the $j^{\mathrm{th}}$ explanatory variable been orthogonal to the remaining variables. The second term $1/(1-R_j^2)$ is a standard definition of the $j^{\mathrm{th}}$VIF as in Thiel (1971).

3. Generalized variance inflation factors

In this section, we introduce the GVIFs as an extension of the classical variance inflation factors $VIF$ from Equation 1. For the linear model $\mathbf{Y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon }$, view $\mathbf{X}=[{\mathbf{X}}_1,{\mathbf{X}}_2]$ partitioned with ${\mathbf{X}}_1$ of dimension $n\times r$ usually consisting of the lower order terms and ${\mathbf{X}}_2$ of dimension $n\times s$ usually consisting of the higher order terms. The idealized moment matrix for the $(r,s)$ partitioning of $\mathbf{X}$ is$$\begin{array}{c}\hfill {\mathbf{M}}_{(r,s)}=\left[\begin{array}{cc}{\mathbf{X}}_1^{\prime }{\mathbf{X}}_1& {\mathbf{0}}_{r\times s}\\ {\mathbf{0}}_{s\times r}& {\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\end{array}\right].\end{array}$$

Following Equation 1, to measure the effect of adding ${\mathbf{X}}_2$ to the design ${\mathbf{X}}_1$, that is for ${\mathbf{X}}_2|{\mathbf{X}}_1$, we define the generalized variance inflation factor as(2) $$\begin{array}{c}\hfill GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)=\frac{det({\mathbf{M}}_{(r,s)})}{det({\mathbf{X}}^{\prime }\mathbf{X})}=\frac{det({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1)det({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2)}{det({\mathbf{X}}^{\prime }\mathbf{X})}\end{array}$$(2)

as in Equation 10 of Fox and Monette (1992), who compared the sizes of the joint confidence regions for $\mathit{\beta }$ for partitioned designs and noted when $\mathbf{X}=[{\mathbf{X}}_{[p]},{\mathbf{x}}_p]$ that $GVIF[{\mathbf{x}}_p|{\mathbf{X}}_{[p]}]=VI{F}_p$. Equation 2 is in the spirit of the efficiency comparisons in linear inferences introduced in Theorems 4 and 5 of Jensen and Ramirez (1993). A similar measure of collinearity is mentioned in Note 2 in Wichers (1975), Theorem 1 of Berk (1977), and Garcia, Garcia, and Soto (2011). For the simple linear regression model with $p=2$, Equation 2 gives $VIF=\frac{1}{1-{\mathit{\rho }}^2}$ with $\mathit{\rho }$ the correlation coefficient as required. Fox and Monette (1992) suggested that ${\mathbf{X}}_1$ contains the variables which are of “simultaneous interest,” while ${\mathbf{X}}_2$ contains additional variables selected by the investigator. We will set ${\mathbf{X}}_1$ for the constant and main effects and set ${\mathbf{X}}_2$ the (optional) quadratic terms with values from ${\mathbf{X}}_1$.

Willan and Watts (1978) measured the effect of collinearity using the ratio of the volume of the actual joint confidence region for $\widehat{\mathit{\beta }}$ to the volume of the joint confidence region in the fictitious orthogonal reference design. Their ratio is in the spirit of $GVIF$ as $det({\mathbf{X}}^{\prime }\mathbf{X})$ is inversely proportional to the square of the volume of the joint confidence region for $\widehat{\mathit{\beta }}$. They also introduced a measure of relative predictability and they note: “The existence of near linear relations in the independent variables of the actual data reduces the overall predictive efficiency by this factor.” For a simple case study, consider the simple linear regression model with $n=4$, ${\mathbf{x}}_1={[-2,-1,1,2]}^{\prime }$, and $\mathbf{y}={[4,1,1,4]}^{\prime }$. The $95\%$ prediction interval for $x_1=0$ is $2.5\pm 10.20$. If the model also includes ${\mathbf{x}}_2={[-2.001,-1.001,1.001,2.001]}^{\prime }$, then the $95\%$ prediction interval for $(x_1,x_2)=(0,0)$ is $2.5\pm 46.02$ demonstrating the loss of predictive efficiency due to the collinearity introduced by ${\mathbf{x}}_2$.

For the $(r,s)$ partition of $\mathbf{X}=[{\mathbf{X}}_1,{\mathbf{X}}_2]$ with ${\mathbf{X}}_1$ of dimension $n\times r$ and ${\mathbf{X}}_2$ of dimension $n\times s$, set$$\begin{array}{c}\hfill {\mathbf{D}}_{(r,s)}=\left[\begin{array}{cc}{\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}& \mathbf{0}\\ \mathbf{0}& {\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}\end{array}\right],\end{array}$$

and denote the canonical moment matrix as(3) $$\begin{array}{cc}\hfill \mathbf{R}& ={\mathbf{D}}_{(r,s)}({\mathbf{X}}^{\prime }\mathbf{X})D_{(r,s)}\hfill \\ \hfill & =\left[\begin{array}{cc}{\mathbf{I}}_{r\times r}& {\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}({\mathbf{X}}_1^{\prime }{\mathbf{X}}_2){\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}\\ {\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}({\mathbf{X}}_2^{\prime }{\mathbf{X}}_1){\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}& {\mathbf{I}}_{s\times s}\end{array}\right];\hfill \end{array}$$(3)

with determinant$$\begin{array}{c}\hfill det(\mathbf{R})=\frac{det({\mathbf{X}}^{\prime }\mathbf{X})}{det({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1)det({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2)}=\frac{1}{GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)};\end{array}$$

equivalently,$$\begin{array}{c}\hfill det(\mathbf{R})=det({\mathbf{I}}_{r\times r}-{\mathbf{B}}_{r\times s}{\mathbf{B}}_{s\times r}^{\prime })=det({\mathbf{I}}_{s\times s}-{\mathbf{B}}_{s\times r}^{\prime }{\mathbf{B}}_{r\times s})\end{array}$$

where ${\mathbf{B}}_{r\times s}={\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}({\mathbf{X}}_1^{\prime }{\mathbf{X}}_2){\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}$.

In the case $\{r=p-1$, $s=1\}$, ${\mathbf{X}}_2={\mathbf{x}}_p$ is a $n\times 1$ vector and the partitioned design $\mathbf{X}=[{\mathbf{X}}_1,{\mathbf{x}}_p]$ has $det(\mathbf{R})=1-[{\mathbf{x}}_p^{\prime }{\mathbf{X}}_1{\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1}{\mathbf{X}}_1^{\prime }{\mathbf{x}}_p]/\left({\mathbf{x}}_p^{\prime }{\mathbf{x}}_p\right)$. From standard facts for the inverse of a partitioned matrix, for example, Myers (1990, p. 459), $VI{F}_p={[{\mathbf{R}}^{-1}]}_{p,p}={[{\mathbf{D}}_{(p-1,1)}^{-1}{({\mathbf{X}}^{\prime }\mathbf{X})}^{-1}{\mathbf{D}}_{(p-1,1)}^{-1}]}_{p,p}$ can be computed directly as$$\begin{array}{cc}\hfill {\left({\mathbf{x}}_p^{\prime }{\mathbf{x}}_p\right)}^{1/2}{({\mathbf{X}}^{\prime }\mathbf{X})}_{p,p}^{-1}{\left({\mathbf{x}}_p^{\prime }{\mathbf{x}}_p\right)}^{1/2}& =\frac{{\mathbf{x}}_p^{\prime }{\mathbf{x}}_p}{{\mathbf{x}}_p^{\prime }{\mathbf{x}}_p-{\mathbf{x}}_p^{\prime }{\mathbf{X}}_1{\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1}{\mathbf{X}}_1^{\prime }{\mathbf{x}}_p}\hfill \\ \hfill & =\frac{1}{1-[{\mathbf{x}}_p^{\prime }{\mathbf{X}}_1{\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1}{\mathbf{X}}_1^{\prime }{\mathbf{x}}_p]/\left({\mathbf{x}}_p^{\prime }{\mathbf{x}}_p\right)}\hfill \\ \hfill & =\frac{1}{det(\mathbf{R})}=GVIF({\mathbf{X}}_2|{\mathbf{X}}_1).\hfill \end{array}$$

Table 1. CCD with parameter $a$, canonical index ${\mathit{\gamma }}_X^2$, and $GVIF$

$a$	${\mathit{\gamma }}_X^2$	$GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)$
0.000	0.4444	1.800
1.000	0.8000	5.000
1.500	0.8858	8.758
$\sqrt{2}$	0.8889	9.000
1.750	0.8514	6.729
$\infty$	0.4444	1.800

We study the eigenvalue structure of ${\mathbf{M}}_{(r,s)}$ in Appendix 1. Let $\{{\mathit{\lambda }}_1\ge {\mathit{\lambda }}_2\ge \dots \ge {\mathit{\lambda }}_{min(r,s)}\ge 0\}$ be the non-negative singular values of ${\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}({\mathbf{X}}_1^{\prime }{\mathbf{X}}_2){\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}$. It is shown in Appendix 1 that an alternative formulation for $GVIF$ is(4) $$\begin{array}{c}\hfill GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)=\prod _{i=1}^{min(r,s)}{(1-{\mathit{\lambda }}_i^2)}^{-1}.\end{array}$$(4)

4. Quadratic model with $p=3$

For the partitioning $\mathbf{X}=[{\mathbf{X}}_r|{\mathbf{X}}_s]$, the canonical moment matrix, Equation 3, has the identity matrices ${\mathbf{I}}_r$, ${\mathbf{I}}_s$ down the diagonal and off-diagonal array ${\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}{\mathbf{X}}_1^{\prime }{\mathbf{X}}_2{\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}$. For the quadratic model $y={\mathit{\beta }}_0+{\mathit{\beta }}_{1}x+{\mathit{\beta }}_2{x}^2$ and partitioning $\mathbf{X}=[\mathbf{1},\mathbf{x}|{\mathbf{x}}^2]$, we have$$\begin{array}{c}\hfill \mathbf{R}=\left[\begin{array}{ccc}1& 0& {\mathit{\rho }}_1\\ 0& 1& {\mathit{\rho }}_2\\ {\mathit{\rho }}_1& {\mathit{\rho }}_2& 1\end{array}\right].\end{array}$$

From Equation 4, $GVIF({\mathbf{x}}^2|\mathbf{1},x)={(1-{\mathit{\lambda }}^2)}^{-1}$ where $\mathit{\lambda }=\sqrt{{\mathit{\rho }}_1^2+{\mathit{\rho }}_2^2}$ is the unique positive singular value of ${[{\mathit{\rho }}_1,{\mathit{\rho }}_2]}^{\prime }$. Denote$$\begin{array}{c}\hfill {\mathit{\gamma }}_X^2={\mathit{\rho }}_1^2+{\mathit{\rho }}_2^2\end{array}$$

as the canonical index with $GVIF({\mathbf{x}}^2|\mathbf{1},\mathbf{x})=\frac{1}{1-{\mathit{\gamma }}_X^2}=\frac{1}{det(\mathbf{R})}$. Surprisingly, many higher order designs also have the off-diagonal entry of the canonical moment matrix with a unique positive singular value with $GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)=\frac{1}{1-{\mathit{\gamma }}_X^2}$ with the collinearity between the lower order terms and the upper order terms as a function of the canonical index ${\mathit{\gamma }}_X^2$.

5. Central composite and factorial designs for quadratic models $(p=6)$

In this section, we compare the central composite design (CCD) $\mathbf{X}$ of Box and Wilson (1951) and the factorial design $\mathbf{Z}.$ The design points are shown in Table A1 of Appendix 2. Both designs are $9\times 6$ and use the quadratic response model$$\begin{array}{c}\hfill y={\mathit{\beta }}_0+{\mathit{\beta }}_1{x}_1+{\mathit{\beta }}_2{x}_2+{\mathit{\beta }}_{11}{x}_1^2+{\mathit{\beta }}_{22}{x}_2^2+{\mathit{\beta }}_{12}{x}_1{x}_2+\mathit{\epsilon }.\end{array}$$

The CCD traditionally uses the value $a=\sqrt{2}$ in four entries, while the factorial design uses the value $a=1$. To study the difference in the designs with these different values, we computed the GVIF to compare the “orthogonality” between the lower order terms ${\mathbf{X}}_1$ of dimension $9\times 3$ and the higher order quadratic terms ${\mathbf{X}}_2$ of dimension $9\times 3$. The off-diagonal ${\mathbf{B}}_{3\times 3}$ entry of $\mathbf{R}$ from Equation 3 in Section 3 has the form$$\begin{array}{c}\hfill {\mathbf{B}}_{3\times 3}=\left[\begin{array}{ccc}{\mathit{\rho }}_1& {\mathit{\rho }}_2& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\end{array}$$

with ${\mathit{\rho }}_1={\mathit{\rho }}_2=\frac{2}{3}\frac{2+a^2}{\sqrt{8+2a^4}}$, canonical index ${\mathit{\gamma }}_X^2={\mathit{\rho }}_1^2+{\mathit{\rho }}_2^2$ and $GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)=\frac{1}{1-{\mathit{\gamma }}_X^2}$ as in the quadratic model case with $p=3$ shown in the Section 4. For Table 1, if $a=1$, then ${\mathit{\rho }}_1={\mathit{\rho }}_2=2/\sqrt{10}$, ${\mathit{\gamma }}_X^2=8/10$, and $GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)=5$. Surprisingly, the classical choice of $a=\sqrt{2}$ gives the largest value for $GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)$, that is the worst value, indicating the greatest collinearity between the lower and higher order terms, as noted in O’Driscoll and Ramirez (in press).

6. Larger designs $(p=10)$

We consider the quadratic response surface designs for(5) $$\begin{array}{cc}\hfill y& ={\mathit{\beta }}_0+{\mathit{\beta }}_1{x}_1+{\mathit{\beta }}_2{x}_2+{\mathit{\beta }}_3{x}_3+\hfill \\ \hfill & {\mathit{\beta }}_{11}{x}_1^2+{\mathit{\beta }}_{22}{x}_2^2+{\mathit{\beta }}_{33}{x}_3^2+{\mathit{\beta }}_{12}{x}_1{x}_2+{\mathit{\beta }}_{13}{x}_1{x}_3+{\mathit{\beta }}_{23}{x}_2{x}_3+\mathit{\epsilon }\hfill \end{array}$$(5)

with $n$ responses and with $\mathbf{X}$ partitioned into $\left[{\mathbf{X}}_1|X_2\right]$ with ${\mathbf{X}}_1$ the four lower order terms $(r=4)$ and ${\mathbf{X}}_2$ the six quadratic terms $(s=6)$. Four popular designs are given in Appendix 2. They are the hybrid designs ($H310$ and $H311B$) of Roquemore (1976) Tables A2 and A3, the Box and Behnken (1960) $(BBD)$ design Table A4, and the CCD of Box and Wilson (1951) Table A5.

For each design, we compute the $10\times 10$ canonical moment matrix. It is striking that, for all these designs, the off-diagonal $4\times 6$ array in $\mathbf{R}$ has only one non-zero singular value with its square the canonical index ${\mathit{\gamma }}_{\mathbf{X}}^2$. It follows that $GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)=\frac{1}{1-{\mathit{\gamma }}_{\mathbf{X}}^2}$.

Table 2. Hybrid designs $H310$, $H311B$, Box and Behnken $BBD,and$ CCD

Design	$n$	$r$	${\mathit{\gamma }}_{\mathbf{X}}^2$	$GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)$
$H310$	11	4	0.8199	5.553
$H311B$	11	4	0.9091	11.00
BBD	13	4	0.9231	13.00
CCD	15	4	0.9333	15.00

Table 3. Singular values for off-diagonal array of $\mathbf{R}$ for BDD and SCD with $GVIF$

	BDD		SCD
	0.9017		0.9535
$\mathit{\alpha }=$	0.3424	$\mathit{\beta }=$	0.6325
	0.3424		0.6325
	0.1877		0.6325
$GVIF=$	7.114		50.93

Table 2 reports that the design $H310$ is the most conditioned with respect to the GVIF with the least amount of collinearity between the lower and higher order terms.

7. More complicated designs with ordered singular values

Let $\mathbf{X}$ be the minimal design of Box and Draper (1974) BDD with $n=11$ from Table A6, and let $\mathbf{Z}$ be the small composite design of Hartley (1959) $SCD$ with $n=11$ from Table A7 for the quadratic response surface model $(r=4$ and $s=6)$ as in Equation (5(5) $$\begin{array}{cc}\hfill y& ={\mathit{\beta }}_0+{\mathit{\beta }}_1{x}_1+{\mathit{\beta }}_2{x}_2+{\mathit{\beta }}_3{x}_3+\hfill \\ \hfill & {\mathit{\beta }}_{11}{x}_1^2+{\mathit{\beta }}_{22}{x}_2^2+{\mathit{\beta }}_{33}{x}_3^2+{\mathit{\beta }}_{12}{x}_1{x}_2+{\mathit{\beta }}_{13}{x}_1{x}_3+{\mathit{\beta }}_{23}{x}_2{x}_3+\mathit{\epsilon }\hfill \end{array}$$(5) ). Let $\mathit{\alpha }=\{{\mathit{\alpha }}_1\ge \dots \ge {\mathit{\alpha }}_r\ge 0\}$ and $\mathit{\beta }=\{{\mathit{\beta }}_1\ge \dots \ge {\mathit{\beta }}_r\ge 0\}$ be the non-negative singular values of the off-diagonal array for ${\mathbf{R}}_X$ and ${\mathbf{R}}_Z$, respectively. As ${\mathit{\alpha }}_i\le {\mathit{\beta }}_i$$(1\le i\le r)$ in Table 3, it follows that $GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)\le GVIF({\mathbf{Z}}_2|{\mathbf{Z}}_1)$ showing less collinearity between the lower and higher order terms for the BDD design.

8. An improved H310 design

When the diagonal matrix ${\mathrm{\Lambda }}_{r\times s}$ in Equation 6 in Appendix 1 has only one non-zero entry, we have denoted the square of this value the canonical index. We extend this definition to the case when ${\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}({\mathbf{X}}_1^{\prime }{\mathbf{X}}_2){\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}$ has multiple positive singular values. The Frobenious norm for a rectangular matrix ${\mathbf{A}}_{r\times s}$ is defined by ${||\mathbf{A}||}_F^2={\sum }_{i=1}^r{\sum }_{j=1}^s{a}_{ij}^2=trace({\mathbf{A}}^{\prime }\mathbf{A})$. For a design matrix $\mathbf{X}$, we extend the definition of the canonical index with ${\mathit{\gamma }}_{\mathbf{X}}^2=||{\mathrm{\Lambda }}_{r\times s}{||}_F^2$. Alternatively, ${\mathit{\gamma }}_{\mathbf{X}}^2=trace({\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1}({\mathbf{X}}_2^{\prime }{\mathbf{X}}_1){\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1}({\mathbf{X}}_1^{\prime }{\mathbf{X}}_2))$ as in Equation 7.

We examine, in detail, the $H310$ design matrix ${\mathbf{X}}_{11\times 10}$, Table A2 in Appendix 2, with our attention to the value of $-0.1360$ in row 2 for ${\mathbf{x}}_3$. In succession, we will replace the values $\{1.1736,0.6386,-0.9273,1.0000,1.2906,-0.1360\}$ by a free parameter and use ${\mathit{\gamma }}_{\mathbf{X}}^2$ to determine an optimal value. For example, replacing the four entries which are $1.1736$ with $c_1$, we calculate the minimum value for ${\mathit{\gamma }}_{\mathbf{X}}^2=0.8199$ with $c_1=1.1768$ denoted $c_{min}$ in Table 4. These values are within the four digit accuracy of the data. We performed a similar calculation with $c_2$ using the four entries which are $0.6386$; with $c_3$ with the four entries which are $-0.9273$; with $c_4$ with the eight entries which are $1$; and with $c_5$ with the single entry $1.2906$. The original design has ${\mathit{\gamma }}_{\mathbf{X}}^2=0.8199$. The entries in the $H310$ design are given to four significant digits. With this precision, the original design is nearly optimal with respect to the canonical index ${\mathit{\gamma }}_{\mathbf{X}}^2$ for the first five entries in Table 4. The sixth entry of $c_6=-0.1360$ was not optimal with ${\mathit{\gamma }}_{\mathbf{X}}^2=0.8181$ with $c_{min}=-0.01264$, a magnitude value smaller.

Table 4. Optimal values $c_{min}$ for ${\mathit{\gamma }}_{\mathbf{X}}^2$

Design entries						$c_{min}$	${\mathit{\gamma }}_{\mathbf{X}}^2$
$c_1$	0.6386	-0.9273	1	1.2906	$-0.1360$	1.1768	0.8199
1.1736	$c_2$	-0.9273	1	1.2906	$-0.1360$	0.6356	0.8199
1.1736	0.6386	$c_3$	1	1.2906	$-0.1360$	$-0.9303$	0.8199
1.1736	0.6386	$-0.9273$	$c_4$	1.2906	$-0.1360$	0.9975	0.8199
1.1736	0.6386	$-0.9273$	1	$c_5$	$-0.1360$	1.2880	0.8199
1.1736	0.6386	$-0.9273$	1	1.2906	$c_6$	$-0.0126$	0.8181

Denote the “improved” $H310$ design as the $H310$ design with the value of $c_6=-0.01264$. The “improved” $H310$ also has a unique positive singular value for the off-diagonal of $\mathbf{R}$ with its square the canonical index ${\mathit{\gamma }}_{\mathbf{X}}^2$. All of the standard design criteria favor the “improved” $H310$ design over the $H310$ design, which was originally constructed based on the rotatability criterion to maintain equal variances for predicted responses for points that have the same distance from the design center. As usual $A(\mathbf{X})=tr({({\mathbf{X}}^{\prime }\mathbf{X})}^{-1})$, $D(\mathbf{X})=det({({\mathbf{X}}^{\prime }\mathbf{X})}^{-1})$, and $E(\mathbf{X})=max\{$eigenvalues of ${({\mathbf{X}}^{\prime }\mathbf{X})}^{-1}\}$. The small relative changes $\mathrm{\Delta }$ in the design criteria are shown in Table 5 in Column 4.

Table 5. Design criteria for the “improved” $H310$ with $\mathrm{\Delta }$ the relative change

	$c_6=-0.1360$	$c_6=-0.01264$	$\mathrm{\Delta }$
$A$	2.697	2.688	0.33%
$D$	$0.1552\times {10}^{-6}$	$0.1542\times {10}^{-6}$	0.64%
$E$	0.8611	0.8512	1.15%
${\mathit{\gamma }}_{\mathbf{X}}^2$	0.8199	0.8181	0.22%

The abnormality of the second row in $H310$ has been noted in Jensen (1998) who showed that the design is least sensitive to the second row of $\mathbf{X},$ the row containing the value $c_6=-0.1360$.

9. Conclusions

The VIF measure the penalty for adding a non-orthogonal variable to a linear regression. The $VIF$ can be computed as a ratio of determinant as in Equation 1. A similar ratio criterion was studied by Fox and Monette (1992) to measure the effect of adding a subset of new variables to a design and they dubbed it the generalized variance inflation factor $GVIF$, Equation 2. We have noted the relationship between $GVIF$ and the singular values of the off-diagonal array in the canonical moment matrix and have used $GFIV$ to study standard quadratic response designs. The $H310$ design of Roquemorer (1976) was shown not to be optimal with respect to $GFIV$ and an “improved” $H310$ design was introduced which was favored over $H310$ using the standard design criteria $A$, $D$, and $E$.

Additional information

Funding

The authors received no direct funding for this research.

Notes on contributors

Diarmuid O’Driscoll

Diarmuid O’Driscoll is the head of the Mathematics and Computer Studies Department at Mary Immaculate College, Limerick. He was awarded a Travelling Studentship for his MSc at University College Cork in 1977. He has taught at University College Cork, Cork Institute of Technology, University of Virginia, and Frostburg State University. His research interests are in mathematical education, errors in variables regression, and design criteria. In 2014, he was awarded a Teaching Heroes Award by the National Forum for the Enhancement of Teaching and Learning (Ireland).

Donald E. Ramirez

Donald E Ramirez is a full professor in the Department of Mathematics at the University of Virginia in Charlottesville, Virginia. He received his PhD in Mathematics from Tulane University in New Orleans, Louisiana. His research is in harmonic analysis and mathematical statistics. His current research interests are in statistical outliers and ridge regression.

Appendix 1

We study the eigenvalue structure of ${\mathbf{M}}_{(r,s)}$. Let $\{{\mathit{\lambda }}_1\ge {\mathit{\lambda }}_2\ge \dots \ge {\mathit{\lambda }}_{min(r,s)}\ge 0\}$ be the non-negative singular values of ${\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}({\mathbf{X}}_1^{\prime }{\mathbf{X}}_2){\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}$.

As with the canonical correlation coefficients Eaton (1983), write the off-diagonal rectangular array ${\mathbf{B}}_{r\times s}$ of $\mathbf{R}$ as $\mathbf{P}\mathrm{\Lambda }{Q}^{\prime }$with $\mathbf{P}$ and $\mathbf{Q}$ orthogonal matrices and ${\mathrm{\Lambda }}_{r\times s}$ the rectangular diagonal matrix with the non-negative singular values down the diagonal. Set$$\begin{array}{c}\hfill \mathbf{L}=\left[\begin{array}{cc}{\mathbf{P}}_{r\times r}& {\mathbf{0}}_{r\times s}\\ {\mathbf{0}}_{s\times r}& {\mathbf{Q}}_{s\times s}\end{array}\right].\end{array}$$

For notational convenience, we assume $r\le s$. The matrix $\mathbf{L}$ is orthogonal and transforms $\mathbf{R}\to$${\mathbf{L}}^{\prime }\mathbf{R}L$ into diagonal matrices:(A1) $$\begin{array}{c}\hfill \left[\begin{array}{cc}{\mathbf{I}}_r& {\mathrm{\Lambda }}_{r\times s}\\ {\mathrm{\Lambda }}_{s\times r}^{\prime }& {\mathbf{I}}_s\end{array}\right]=\left[\begin{array}{cc}{\mathbf{I}}_r& [\mathbf{S}{V}_{r\times r}|0_{r\times (s-r)}]\\ {[\mathbf{S}{V}_{r\times r}|0_{r\times (s-r)}]}^{\prime }& {\mathbf{I}}_s\end{array}\right]\end{array}$$(A1)

with ${\mathrm{\Lambda }}_{r\times s}=[\mathbf{S}{V}_{r\times r}|0_{r\times (s-r)}]$ where $\mathbf{S}{V}_{r\times r}$ is the diagonal matrix of the non-negative singular values. Since $\mathbf{L}$ is orthogonal, this transformation has not changed the eigenvalues. To compute the determinant of $\mathbf{R}$, convert the matrix in Equation 6 into an upper diagonal matrix by Gauss Elimination on ${\mathrm{\Lambda }}_{s\times r}^{\prime }$. This changes $r$ of the $1^{\prime }s$ on the diagonal in rows $r+1$ to $r+r$ into $1-{\mathit{\lambda }}_i^2$, and thus $det(\mathbf{R})={\prod }_{i=1}^{min(r,s)}(1-{\mathit{\lambda }}_i^2)$ with$$\begin{array}{c}\hfill GVIF({\mathbf{X}}_2|{\mathbf{X}}_1)=\prod _{i=1}^{min(r,s)}\frac{1}{1-{\mathit{\lambda }}_i^2}.\end{array}$$

The singular values of ${\mathbf{R}}_{12}={\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1/2}({\mathbf{X}}_1^{\prime }{\mathbf{X}}_2){\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1/2}$ are the non-negative square roots of the eigenvalues of ${\mathrm{\Lambda }}^{\prime }\mathrm{\Lambda }$ denoted by(A2) $$\begin{array}{cc}\hfill eigvals({\mathrm{\Lambda }}^{\prime }\mathrm{\Lambda })& =eigvals(({\mathbf{Q}}^{\prime }{\mathbf{R}}_{12}^{\prime }\mathbf{P})(P^{\prime }{\mathbf{R}}_{12}\mathbf{Q})))\hfill \\ \hfill & =eigvals({\left({\mathbf{X}}_1^{\prime }{\mathbf{X}}_1\right)}^{-1})({\mathbf{X}}_1^{\prime }{\mathbf{X}}_2){\left({\mathbf{X}}_2^{\prime }{\mathbf{X}}_2\right)}^{-1}({\mathbf{X}}_2^{\prime }{\mathbf{X}}_1)).\hfill \end{array}$$(A2)

If the trace of the inverse of the matrix in Equation 6 is required, then we note that$$\begin{array}{c}\hfill {\left[\begin{array}{cc}{\mathbf{I}}_r& {\mathrm{\Lambda }}_{r\times s}\\ {\mathrm{\Lambda }}_{s\times r}^{\prime }& {\mathbf{I}}_s\end{array}\right]}^{-1}=\left[\begin{array}{cc}{({\mathbf{I}}_r-{\mathrm{\Lambda }}_{r\times s}{\mathrm{\Lambda }}_{s\times r}^{\prime })}^{-1}& -{\mathrm{\Lambda }}_{r\times s}{(I_s-{\mathrm{\Lambda }}_{s\times r}^{\prime }{\mathrm{\Lambda }}_{r\times s})}^{-1}\\ -{\mathrm{\Lambda }}_{s\times r}^{\prime }{(I_r-{\mathrm{\Lambda }}_{r\times s}{\mathrm{\Lambda }}_{s\times r}^{\prime })}^{-1}& {({\mathbf{I}}_s-{\mathrm{\Lambda }}_{s\times r}^{\prime }{\mathrm{\Lambda }}_{r\times s})}^{-1}\end{array}\right]\end{array}$$

with trace given by $tr({({\mathbf{L}}^{\prime }\mathbf{R}L)}^{-1})=|r-s|+2{\sum }_{i=1}^{min(r,s)}\frac{1}{1-{\mathit{\lambda }}_i^2}$.

Appendix 2

Table A1. The lower order matrix for the CCD with center run with $a=\sqrt{2}$, $n=9$ and the lower order matrix for the factorial design with center run $n=9$

Central composite design			Factorial design
$\mathbf{1}$	${\mathbf{x}}_1$	${\mathbf{x}}_2$	$\mathbf{1}$	${\mathbf{x}}_1$	${\mathbf{x}}_2$
1	1	1	1	1	1
1	1	-1	1	1	-1
1	-1	1	1	-1	1
1	-1	-1	1	-1	-1
1	$a$	0	1	1	0
1	$-a$	0	1	-1	0
1	0	$a$	1	0	1
1	0	$-a$	1	0	-1
1	0	0	1	0	0

Table A2. The lower order matrix for the hybrid $(H310)$ design of Roquemore (1976) with center run, $n=11$

$\mathbf{1}$	${\mathbf{x}}_1$	${\mathbf{x}}_2$	${\mathbf{x}}_3$
1	0	0	1.2906
1	0	0	-0.1360
1	1	1	0.6386
1	1	-1	0.6386
1	-1	1	0.6386
1	-1	-1	0.6386
1	1.1736	0	-0.9273
1	-1.1736	0	-0.9273
1	0	1.1736	-0.9273
1	0	-1.1736	-0.9273
1	0	0	0

Table A3. The lower order matrix for the hybrid $(H311B)$ design of Roquemore (1976) with center run, $n=11$

$\mathbf{1}$	${\mathbf{x}}_1$	${\mathbf{x}}_2$	${\mathbf{x}}_3$
1	0	0	2.4495
1	0	0	-2.4495
1	2.1063	0.7507	1
1	2.1063	-0.7507	-1
1	-2.1063	-0.7507	1
1	-2.1063	0.7507	-1
1	0.7507	2.1063	-1
1	0.7507	-2.1063	1
1	-0.7507	-2.1063	-1
1	-0.7507	2.1063	1
1	0	0	0

Table A4. The lower order matrix for the Box and Behnken (1960) design $(BBD)$ with center run, $n=13$

$\mathbf{1}$	${\mathbf{x}}_1$	${\mathbf{x}}_2$	${\mathbf{x}}_3$
1	1	1	0
1	1	-1	0
1	1	0	1
1	1	0	-1
1	-1	1	0
1	-1	-1	0
1	-1	0	1
1	-1	0	-1
1	0	1	1
1	0	1	-1
1	0	-1	1
1	0	-1	-1
1	0	0	0

Table A5. The lower order matrix for the Box and Wilson (1951) CCD for $\mathit{\alpha }=1.732$ with center run, $n=15$

$\mathbf{1}$	${\mathbf{x}}_1$	${\mathbf{x}}_2$	${\mathbf{x}}_3$
1	1	1	1
1	1	1	-1
1	1	-1	1
1	1	-1	-1
1	-1	1	1
1	-1	1	-1
1	-1	-1	1
1	-1	-1	-1
1	1.732	0	0
1	-1.732	0	0
1	0	1.732	0
1	0	-1.732	0
1	0	0	1.732
1	0	0	-1.732
1	0	0	0

Table A6. The lower order matrix for the Box and Draper (1974) minimal design (BDD) with center run, $n=11$

$\mathbf{1}$	${\mathbf{x}}_1$	${\mathbf{x}}_2$	${\mathbf{x}}_3$
1	1	-1	-1
1	-1	1	-1
1	-1	-1	1
1	-1	-1	-1
1	0.1925	0.1925	-1
1	0.1925	-1	0.1925
1	-1	0.1925	0.1925
1	-0.2912	1	1
1	1	-0.2912	1
1	1	1	-0.2912
1	0	0	0

Table A7. The Lower order matrix for the small composite design of Hartley (1959) $(SCD)$ for $\mathit{\alpha }=1.732$ with center run, $n=11$

$\mathbf{1}$	${\mathbf{x}}_1$	${\mathbf{x}}_2$	${\mathbf{x}}_3$
1	1	1	1
1	1	-1	-1
1	-1	1	-1
1	-1	-1	1
1	1.732	0	0
1	-1.732	0	0
1	0	1.732	0
1	0	-1.732	0
1	0	0	1.732
1	0	0	-1.732
1	0	0	0