Automated Analysis of Experiments using Hierarchical Garrote

Consider approximate the latent function $f$ by a main-effect-only linear model,

Let $\bm{\beta}=(\beta_{1},\beta_{2},\dots,\beta_{p})^{\prime}$. The least squares estimator $\hat{\bm{\beta}}^{LS}=(\mathbf{X}^{\prime}\mathbf{X})^{-1}\mathbf{X}^{\prime}\mathbf{Y}$. The original nonnegative garrote (NG) (Breiman,, 1995) works by scaling the least squares estimate with shrinkage factors $\bm{\theta}(M)$ controlled by a tuning parameter $M$. The NG estimate is defined as $\hat{\bm{\beta}}^{NG}(M)=(\hat{\beta}_{1}^{NG}(M),\dots,\hat{\beta}_{p}^{NG}(M))^{\prime}$, where $\hat{\beta}_{j}^{NG}(M)=\hat{\theta}_{j}(M)\hat{\beta}_{j}^{LS}$. For an $M\geq 0$, the shrinkage factors are obtained by solving a quadratic programming problem:

By selecting $M$ suitably, some elements of $\hat{\bm{\beta}}^{NG}(M)$ can be exactly zero. Since $\hat{\beta}_{j}^{LS}\neq 0$ with a probability of 1, $x_{j}$ is considered active if and only if $\hat{\theta}_{j}(M)>0$. In other words, $\hat{\theta}_{j}(M)$ serves as an indicator of including $x_{j}$ in the selected model.

In several problems, main effects models are inadequate to approximate the underlying function $f$ in (1). To improve the approximation, we can include the quadratic effects and two-factor interactions:

The shrinkage factors in the nonnegative garrote are obtained by

In this optimization, it is possible to select a model that includes an interaction term while excluding both of its main effects, that is, $\hat{\theta}_{ij}(M)>0$ and $\hat{\theta}_{i}(M)=\hat{\theta}_{j}(M)=0$. To preserve the effect heredity principle in nonnegative garrote, Yuan et al., (2009) proposed to add constraints that enforce strong or weak heredity among the effects.

Under the strong heredity principle, we must have $\theta_{i}>0$ and $\theta_{j}>0$ if $\theta_{ij}>0$. This can be enforced by requiring $\theta_{ij}\leq\text{min}(\theta_{i},\theta_{j}),$ which is equivalent to two convex constraints

Thus, if $\hat{\theta}_{ij}(M)>0$, (4) will guarantee that $\hat{\theta}_{i}(M)>0$ and $\hat{\theta}_{j}(M)>0$, ensuring that $x_{i}$ and $x_{j}$ are part of the chosen model.

Under the weak heredity principle, we must have at least one of $\theta_{i}$ or $\theta_{j}$ to be larger than 0 if $\theta_{ij}>0$. This can be enforced by requiring $\theta_{ij}\leq\text{max}(\theta_{i},\theta_{j}).$ However, this constraint is nonconvex, and the optimization in (3) cannot be solved using quadratic programming. To simplify the optimization, Yuan et al., (2009) proposed a convex relaxation of the constraint:

Thus, if $\hat{\theta}_{ij}(M)>0$, (5) will ensure that at least one of $\hat{\theta}_{i}(M)>0$ or $\hat{\theta}_{j}(M)>0$ holds, thereby guaranteeing that at least one of $x_{i}$ or $x_{j}$ is included in the selected model.

It is more common in the analysis of experiments (Wu and Hamada,, 2021) to approximate $f$ using a linear model with both main effects (linear and quadratic) and their cross product terms (two-factor interactions):

where $u_{1}=x_{1},\dots,u_{p}=x_{p},u_{p+1}=x_{1}^{2},\dots,u_{2p}=x_{p}^{2},\ldots,u_{P}=x_{p-1}^{2}x_{p}^{2}$, and $P=2p+\binom{2p}{2}$. For an experiment, the data size $n$ is typically much smaller than the number of parameters $P$, making it impossible to obtain the least squares estimate required in the original nonnegative garrote. Therefore, Yuan and Lin, (2006) suggested using the ridge regression to obtain an initial estimate in NG.

Let $\mathbf{U}_{P}$ be an $n\times P$ model matrix. The ridge regression estimate is given by

where $\lambda$ is a tuning parameter. Then, the NG estimate is defined as $\hat{\beta}_{j}^{NG}=\hat{\theta}_{j}(M)\hat{\beta}_{j}^{R},j=1,\dots,P$, where the shrinkage factors are obtained as

To tune the hyperparameter $M$, we employ the generalized cross-validation (GCV) criterion introduced by Golub et al., (1979), which was extended by Xiong, (2010) to handle NG with ridge estimate. The GCV criterion is defined as

where $d(\hat{\bm{\beta}}^{NG})$ is the degrees of freedom of $\hat{\bm{\beta}}^{NG}$. The degrees of freedom is computed as

where $w_{i}$ is the $(i,i)$ entry of the $P\times P$ matrix $(\mathbf{U}_{P}^{\prime}\mathbf{U}_{P}+\lambda\mathbf{I}_{P})^{-1}\mathbf{U}_{P}^{\prime}\mathbf{U}_{P}$. After constructing the solution path across a series of $M\in[0.1,n-1]$, the optimal $M$ is selected by minimizing the GCV criterion,

The upper bound of $M$ is set to $n-1$ to ensure that the number of selected terms is not too high compared to the degrees of freedom in the data.

However, we will show that ridge regression estimate is inadequate in experimental data analysis. Consider a toy example with the 12-run Plackett-Burman (PB) design. The response is generated by the true model $y=20A+10AB+5AC$. See Table 1 for the design matrix and data. We use the R package glmnet to perform the ridge regression and apply the weak heredity constraint when solving the quadratic programming problem (8). This approach identifies only the effect $A$ with $\hat{\beta}_{A}^{NG}=14.938$ and misses the two interactions $AB$ and $AC$ that are present in the true model. The issue arises from the underestimation of $\hat{\bm{\beta}}^{R}_{P}$. Due to the presence of numerous independent and irrelevant effects in the model matrix, we found that $\lambda$ is often overestimated, masking the effects of important terms. In the next section, we will propose an improved initial estimate for NG to overcome the foregoing issue.

Let the factor $x_{j}$ have $m_{j}$ levels in the experiment, $j=1,\ldots,p$. Consider the model in (6). Although it has only the main effects (linear and quadratic) and their cross product terms (two-factor interactions, 2fi), it is easier to construct the prior for all possible effects. Thus, we approximate $f$ by a linear model including main effects and all the interactions, i.e.,

where $q=\prod_{i=1}^{p}m_{i}$ and $u_{0}=1.$ To avoid confusion, we let $\bm{\beta}_{q}=(\beta_{0},\beta_{1},\dots,\beta_{q-1})^{\prime}$.

The main effects of the factor $x_{j}$ can be described by $m_{j}-1$ dummy variables, and the interactions can be expressed as the products of these dummy variables. For example, consider an experiment with two factors each at three levels. Let $u_{1}$ and $u_{2}$ be the main effects of the first factor and $u_{3}$ and $u_{4}$ be the main effects of the second factor. The interaction terms are $u_{5}=u_{1}u_{3}$, $u_{6}=u_{1}u_{4}$, $u_{7}=u_{2}u_{3}$, and $u_{8}=u_{2}u_{4}$. Many popular coding systems can be used to define dummy variables, such as treatment coding, Helmert coding, and orthogonal polynomial coding (Wu and Hamada,, 2021; Harville,, 1998).

The idea of functionally induced priors (Joseph,, 2006) is to postulate a GP prior for $f$ and use that to induce the prior for $\bm{\beta}_{q}$. Thus, assume

where $\nu^{2}\psi$ is the covariance function of the Gaussian process defined as $\text{cov}\{f(\mathbf{x}),f(\mathbf{x+h})\}=\nu^{2}\psi(\mathbf{h})$. Thus, we have $\mathbf{f}=\mathbf{U}_{q}\bm{\beta}_{q},$ where $\mathbf{f}$ is the vector of function values for the full factorial design and $\mathbf{U}_{q}$ the $q\times q$ model matrix corresponding to (14). Note that $\mathbf{f}$ has a multivariate normal distribution with $E(\mathbf{f})=\mathbf{0}$ and $\text{var}(\mathbf{f})=\nu^{2}\bm{\Psi}$, where $\bm{\Psi}$ is the corresponding correlation matrix. Since $\bm{\beta}_{q}=\mathbf{U}_{q}^{-1}\mathbf{f}$, we obtain

Assume that the correlation function $\psi$ has a product correlation structure, i.e.,

Under the product correlation structure, Joseph and Delaney, (2007) have shown that

where $\bigotimes$ denotes the Kronecker product, $\mathbf{U}_{j}$ is the model matrix for factor $x_{j}$, and $\bm{\Psi}_{j}$ the corresponding correlation matrix. For example, suppose that factor $x_{j}$ is a three-level factor with levels 1, 2, and 3, the model matrix using orthogonal polynomial coding is

and the correlation matrix is

The benefit of this structure is that we can select the most appropriate coding system and correlation function for each factor according to our modeling requirements. Due to its popularity, we use the Gaussian correlation function for each factor (Santner et al.,, 2018): $\psi_{j}(h_{j})=\text{exp}(-h_{j}^{2}/\theta_{j}^{2})$, $0<\theta_{j}<\infty.$ Let $\rho_{j}=\exp(-1/\theta_{j}^{2})$. Then, the correlation function can be written as

Although the form of the correlation function is the same for all factors, the definition of $h_{j}$ should depend on the type of factors.

For quantitative factors, we define $h_{j}$ as $|\mathbf{x}_{ij}-\mathbf{x}_{kj}|$, for two runs $i$ and $k$. Suppose that there are $p$ quantitative three-level factors and their model matrices are encoded using orthogonal polynomial coding. Then, (16) simplifies to (Joseph and Delaney,, 2007):

where

$l_{ij}=1$ if $\beta_{i}$ includes the linear effect of factor $j$ and 0 otherwise, and $q_{ij}=1$ if $\beta_{i}$ includes the quadratic effect of factor $j$ and 0 otherwise. Although the prior for $\bm{\beta}_{q}$ does includes covariance terms, we do not need them to construct the diagonal matrix $\mathbf{R}$ in (13).

Since $r_{j_{l}},\;r_{j_{q}}\in(0,1)$, the variance of an interaction effect will be smaller than its parent effects. In other words, the lower-order effects are likely to be more important than the higher-order effects, and thus the priors satisfy the effect hierarchy principle. Moreover, since the variance of an interaction effect is proportional to the product of the variances of its parent effects, the priors also satisfy the effect heredity principle.

For qualitative factors where the relative distances between levels cannot be quantified, Joseph and Delaney, (2007) defines $h_{j}=0$ if $\mathbf{x}_{ij}=\mathbf{x}_{kj}$ and 1 otherwise, i.e., $h_{j}=H(\mathbf{x}_{ij},\mathbf{x}_{kj})$, where $H(\cdot,\cdot)$ is the Hamming distance. However, this definition will assign equal importance to all the main effects of $x_{j}$ even if not all of them are active. Since the main effects are represented by dummy variables, we define the relative distance for each dummy variable of $x_{j}$ as $h_{j_{d}}=H(\mathbf{x}_{ij_{d}},\mathbf{x}_{kj_{d}}),d=1,\dots,m_{j}-1$, where the subscript $d$ stands for the $d$th dummy variable. Then, according to (18), the form of the corresponding correlation function would be $\psi_{j_{d}}(h_{j_{d}})=\rho_{j_{d}}^{h_{j_{d}}^{2}},\;0<\rho_{j_{d}}<1.$

Suppose that there are $p$ qualitative three-level factors and their model matrices are encoded according to Helmert coding. Then, as shown in the Appendix, the marginals of the prior in (16) will become

where

$m_{ij_{1}}=1$ if $\beta_{i}$ includes the first main effect of factor $j$ and 0 otherwise, and $m_{ij_{2}}=1$ if $\beta_{i}$ includes the second main effect of factor $j$ and 0 otherwise. As in the case of quantitative factors, effect hierarchy and heredity are built-in in these priors.

Our aim is to obtain an initial estimate of $\bm{\beta}_{P}$ from $\mathbf{Y}=\mathbf{U}_{P}\bm{\beta}_{P}+\epsilon$, where $\epsilon\sim\mathcal{N}(\bm{0},\sigma^{2}\mathbf{I})$ and $\bm{\beta}_{P}\sim\mathcal{N}(\mathbf{0},\tau^{2}\mathbf{R})$. Here $\mathbf{R}$ is a diagonal matrix, which is constructed as in (19) and (20). The posterior mean of $\bm{\beta}_{P}$ is given by

To obtain the initial estimate of $\bm{\beta}_{P}$ using (21), several parameters need to be specified, which include $\sigma^{2}$, $\tau^{2}$, and $\bm{\rho}$ that define $\mathbf{R}$. Let $\bm{\rho}=(\bm{\rho}_{1}^{\prime},\dots,\bm{\rho}_{p}^{\prime})^{\prime}$ where $\bm{\rho}_{j}$ can be a vector. A good estimate of $\sigma^{2}$ can be obtained if the experiment contains replicates. Suppose that there are $m$ replicates in each run and let $s_{i}^{2}$ be the sample variance for the $i$th run, $i=1,\dots,n$. Then $\hat{\sigma}^{2}=\sum_{i=1}^{n}s_{i}^{2}/n$ can be a good estimate of $\sigma^{2}$. However, if the experiment is unreplicated, there will be no information to estimate $\sigma^{2}$. Instead of assuming that there is no noise, we assume that the noise would contribute at least $1\%$ of the variation in the response to prevent overfitting. Therefore, we reparameterize (21) to

where $\lambda=\frac{\sigma^{2}}{\sigma^{2}+\nu^{2}}\geq 0.01$. This form can also be used for replicated experiments as long as we replace $\mathbf{Y}$ with the sample means and divide $\mathbf{I}_{n}$ by $m$.

Joseph and Delaney, (2007) have shown that if we use orthogonal coding for each factor,

where $\text{sum}(\mathbf{\Psi}_{j})$ is the sum of all elements in $\mathbf{\Psi}_{j}$. Therefore, we only need to focus on the specifications of $\bm{\rho}$ and $\lambda$. These hyperparameters can be estimated by minimizing the negative log-likelihood of $\mathbf{Y}$ (Santner et al.,, 2018). Thus,

where $\bm{\Psi}_{n}$ is the correlation matrix of the $n$ runs.

Different from the usual applications of GP modeling, we are not only interested in the accurate prediction of the output but also the selection of important effects. Therefore, it is very important for us to obtain the global optimum of (23). This is not easy because when $n$ is small, several models can fit the data well and therefore, the likelihood can be multi-modal with several local optima. While numerous global optimization algorithms are available for this purpose, we need one that requires only a few function evaluations and is robust to the choice of starting points. Therefore, we adopt a gradient-based multi-start global optimization strategy.

For numerical stability, we put mild constraints on the hyperparameters: $\rho_{i}\in[10^{-15},.999]$ for all $i=1,\ldots,k$, and $\lambda\in[.01,.99]$. First, we generate a space-filling design with $k+1$ starting points in the feasible region. Here, we use the MaxPro design in Joseph et al., (2015). Second, for each starting point, we search for the local optimum using a local optimization algorithm. Here we use the Method of Moving Asymptotes in Svanberg, (2002). Among those local optima, the one with the lowest negative log-likelihood value will be used as the estimate of the hyperparameters. Importantly, we choose a derivative-based approach not only because gradients help us reduce the number of function evaluations but also because it gives us a more robust result than other derivative-free algorithms.

Once the hyperparameters are estimated from (23), the shrinkage factors and the regression coefficients can be obtained as in (8) using the initial estimates from (22). The tuning parameter $M$ in (8) can be selected as described in Section 2.3. A solution path across a range of $M\in[0.1,0.3(n-1)]$ is first constructed. The optimal $M$ is then selected by minimizing the GCV criterion (9), where the degrees of freedom is calculated as in (10) with

The upper bound of $M$ used in the optimization is based on the findings of Box and Meyer, (1986) and Wu and Hamada, (2021, p. 424), which indicate that the proportion of active effects (with intercept term included) typically ranges between 0.13 and 0.27 of the design’s run size.

For replicated experiments where a reliable estimate of $\sigma^{2}$ can be obtained, the shrinkage factors can be directly computed by solving (Xiong,, 2010):

Algorithm 1 represents the proposed Hierarchical Garrote (HiGarrote) method. Most of the computational overhead of the method is in step 1 of the algorithm, i.e., the optimization in (23).

To demonstrate the effectiveness of HiGarrote, a simulation study is conducted using the toy example described in Section 2.3, based on a 12-run PB design. The response is generated from the model

The simulation is repeated 100 times, with new data generated each time. For each simulated dataset, both HiGarrote and the nonnegative garrote (NG) with ridge estimate (Yuan et al.,, 2009) are fitted using weak heredity.

Figure 1 presents the boxplots of the estimates of selected effects using NG with ridge estimate and HiGarrote.

As shown in Figure 1, the NG with ridge estimate consistently identifies only the main effect A and occasionally AB. However, the estimate for A shows large variability, and the estimate for AB remains close to 0. Importantly, the crucial interaction AC is never selected, indicating its limitation in identifying weaker yet significant interactions. In contrast, HiGarrote successfully identifies the important effects A, AB, and AC, with estimates close to their true values of 20, 10, and 5. While HiGarrote occasionally selects additional effects such as B, C, and BC, their estimates are generally very small (near 0) with low variability.