Journal of Statistical Theory and Applications

<Previous Article In Issue
Download article (PDF)
Next Article In Issue>
Volume 18, Issue 3, September 2019, Pages 212 - 221

Simultaneous Optimization of Multiple Responses That Involve Correlated Continuous and Ordinal Responses According to the Gaussian Copula Models

Authors
Fatemeh Jiryaie, Ahmad Khodadadi*
Department of Statistics, Shahid Beheshti University, 1983963113 G. C., Tehran, Iran
*Corresponding author. Email: a_khodadadi@sbu.ac.ir
Corresponding Author
Ahmad Khodadadi
Received 14 February 2017, Accepted 14 June 2017, Available Online 9 September 2019.
DOI
10.2991/jsta.d.190701.001How to use a DOI?
Keywords
Gaussian copula; Mixed outcomes; Multivariate distribution; Simultaneous optimization
Abstract

This study investigates the simultaneous optimization of multiple correlated responses that involve mixed ordinal and continuous responses. The proposed approach is applicable for responses that have either an all ordinal categorical form are continuous but have different marginal distributions, or when standard multivariate distribution of responses is not applicable or does not exist. These multiple responses have rarely been the focus of studies despite their high occurrence during experiments. The copula functions have been used to construct a multivariate model for mixed responses. To resolve the computational problems of estimation under a high dimension of responses, we have estimated parameters of the model according to a pairwise likelihood estimation method. We adapted the generalized distance approach to determine settings of the factors that simultaneously optimized the mean of continuous responses and desired cumulative categories of the ordinal responses. A simulation study was used to evaluate the performance of the estimators from the pairwise likelihood approach. Finally, we presented an application of the proposed method in a real data example of a semiconductor manufacturing process.

Copyright
© 2019 The Authors. Published by Atlantis Press SARL.
Open Access
This is an open access article distributed under the CC BY-NC 4.0 license (http://creativecommons.org/licenses/by-nc/4.0/).

1. INTRODUCTION

Improvements to product or process performance are important problems in pharmacology, agriculture, and industry, which is the focus of attention by manufacturers. Consequently, of particular interest is detection of optimal settings for the control factors at which the response presents certain desired characteristics. Several publications focus on optimizing a single response (Draper [1]; Hoerl [2]; Peace [3]; Fowlkes and Creveling [4]; Paul and Khuri [5]). However, numerous situations exist where multiple responses need to be simultaneously optimized. For example, in clinical trials it is important to determine the combination of drugs with maximum therapeutic effects and the lowest level of toxicity. In a semiconductor manufacturing process minimizing the defect count in the sensitive area and achieving the target amount of an ion implanted in a wafer may simultaneously need to be considered for an ion implantation process.

Most quality improvement studies research approaches for simultaneous optimization of multiple responses. The desirability function is one of the most popular methods to optimize a multi-response system (Harrington [6]; Derringer and Suich [7]; Kim and Lin [8]). This approach, turns the several responses into a single response, which results in a combined desirability. The desirability function approach is easy to apply and allows the user to make a subjective judgment on the importance of each response; however, this approach does not take into account the variance–covariance structure of the responses. Ignoring the possible correlations between the responses may be misleading and leads to incorrect optimization decisions. To overcome this difficulty, Elsayed and Chen [9], Ko et al. [10], Pignatiello [11], Tsui [12], and Vining [13] have proposed the use of the loss function approach to optimize multiple correlated responses. Khuri and Conlon [14] introduced an efficient optimization algorithm based on a generalized distance approach. These researchers assumed that all mean responses in the system depended on the same set of controllable variables via a polynomial regression model. The first step of their algorithm was to obtain individual optima of the estimated responses over the experimental region, after which they measured the deviation from the ideal optimum by means of a distance function expressed in terms of the estimated responses along with their variance–covariance structure. Finally, this function could be minimized to arrive at a set of suitable operating conditions.

The following regression models are used in the mentioned multiple response optimization procedure: ordinary least squares (OLS), generalized least squares (GLS), and multivariate regression (MVR), all of which are under the assumption of normal errors. OLS and GLS regressions consider modeling responses individually with the assumption of independent responses. OLS regression is under the homogeneity of error variances, whereas GLS is free from error variances. MVR models the responses simultaneously, by taking into consideration a correlation between the errors. Normality assumption of responses and homogeneity of error variances may be violated in situations where the responses are not normal, discrete, or exhibit heterogeneous variances. Such situations happen frequently in clinical and epidemiological studies. Mukhopadhyay and Khuri [15] have recently modified Khuri and Conlon's [14] algorithm to a generalized linear model (GLM), which can be used to handle multiple discrete responses and bypass the heterogeneity of error variances. They assume all margin distributions are the same and the joint distribution of responses are available and belong to the multivariate exponential family.

A number of different proposed methods optimize multiple continuous responses. Su and Tong [16] have used principle component analysis. However, Lai and Chang [17], Lu and Antony [18], and Tong and Su [19] used the fuzzy theorem whereas Wu and Chyu [20] suggested a mathematical programming method.

However, ordinal responses observed are observed in number of experiments due to the quality characteristic or the convenience of the measurement technique and cost-effectiveness. In the optimization of ordinal responses, Taguchi [21–23] primarily employed the accumulation analysis (AA) method. In this method the corresponding cumulative categories are defined, then the researcher determines the effects of the factor levels according to the probability distribution by the categories. Finally, the optimal control factor settings are obtained by the desired cumulative category and the important location effects are taken into consideration. Nair [24] has proposed two scoring schemes that separately detect the location and dispersion effects. Jeng and Guo [25] presented a weighted probability-scoring scheme (WPSS) to avoid the computational complexity of Nair's scoring scheme. Thy considered the location and dispersion effects. Chipman and Hamada [26] used a GLM with Bayesian estimation techniques to optimize these type of responses. Computational complexity was more than Nair's scoring scheme.

In case of mixed responses very few studies have been conducted that optimize the ordinal-continuous responses. Hsieh and Tong [27] have employed an artificial neural network technique to optimize the ordinal-continuous responses, a method that is employed with difficulty in industrial settings. Wu [28] has presented an approach based on the quality loss function of Taguchi [23] where ordinal responses can be treated as continuous responses and a weighted average quality loss is defined for the ordinal responses. This approach is easier than the approach by Hsieh and Tong [27], however there is no correlation between the responses.

In this article, we introduce an approach to simultaneously optimize mixed correlated continuous and ordinal responses. Our procedure can be easily applied when responses are all ordinal, all continuous with different types of marginal distributions, or in cases where standard multivariate distribution of responses is not applicable or does not exist. For example, when the entire marginal distribution of responses is gamma and the responses are correlated, it is difficult to determine the multivariate exponential distribution. In this approach we have used the Gaussian copula function. We extended the regression models for a bivariate mixed outcomes of De Leon and Wu [29] to the multivariate mixed discrete and continuous outcomes through pairwise fitting of models for the joint modeling of a multivariate mixed outcome based on the concept by Fieuws and Verbeke [30]. After specifying the effects of the factor levels on the mean continuous responses and probability distributions by the categories of the ordinal responses, we have adopted the generalized distance approach of Khuri and Conlon [14] to carry out the optimal control factor settings by mean of continuous responses and desired cumulative categories of ordinal responses. Copula-based dependencies, introduced in statistical literature by Skalar [31], allows one to model the dependence structure independently of marginal distributions. This approach provides an alternative and more useful representation of multivariate distribution compared to traditional approaches such as multivariate normality. Formally, copula can be defined as follows: Suppose that we have K marginal CDFs, FX1⋅,…,FXK⋅, where X1,…,XK are the random variables. Sklar's theorem states that every k-dimensional cumulative distribution FX1,…,XK⋅=PX1≤x1,…,XK≤xk of a random vector X1,…,XK can be expressed by involving only the marginals FX1x1,…,FXKxk as FX1,…,XKx1,…,xk=CFX1x1,…,FXKxk, where C is a copula. Gaussian copulas are an important family which has been used in a variety of applications Song [32]. The P-dimensional Gaussian copula is defined as

CFX1x1,…,FXKxK=ΦKΦ−1FX1x1,…,Φ−1FXKxK;R,
where Φ−1⋅ is the inverse function of the standard normal distribution Φ⋅ and ΦK⋅;R is the k-dimensional standard multivariate normal distribution function. A more through overview can be found in the reference works by Joe [33] or by Nelsen [34].

We considered the real data obtained from a semiconductor manufacturing process in which the defect counted on the sensitive area (an ordinal response) and the amount of ion implanted (a continuous response) require simultaneous investigation for an ion implantation process (Hsieh and Tong [27]), as discussed in Section 6.

This paper is organized as follows: We define a multivariate model for mixed responses in Section 2. In Section 3, parameters of regression models and variance–covariance of the parameters are simultaneously estimated. The confidence region of the parameters, estimated mean of continuous responses, and estimated desired cumulative categories of ordinal responses are obtained in this section. We use all of these for the optimization algorithm. Section 4 outlines the optimization algorithm according to a generalized distance approach. In Section 5, we have conducted a simulation study to compare the performance of estimators from pairwise and full likelihood estimation. An application of the proposed optimization algorithm is described in Section 6 with a real data example. Finally, concluding remarks are presented in Section 7.

2. A MULTIVARIATE MODEL FOR MIXED RESPONSES

We consider a mixed multi-response obtained from the ith run of the experiment, i=1,…,n, is Yi=Yi1,…,YiP,Zi1,…,ZiQ′, where Yips, for p=1,…,P are continuous responses and Ziqs, for q=1,…,Q are ordinal responses with Kq levels. Underling Ziq is Yiq*, a continuous latent variable, such that

Ziq=kq⇔γqkq<Yiq*≤γqkq+1,kq=0,…,Kq−1,
where γq0=−∞, γqKq=+∞ and γq=γq1,…,γqKq−1′ is an unknown vector of the threshold parameters. Let cumulative distribution function (CDF) of Yip and Yiq* be FYip.;μpxi,θp and FYiq*.;μq*xi,θq* respectively, where μpxi=EYip=hpδp′xiβp, μq*xi=EYiq*=hqδq*′xiβq*. Here, δpxi and δq*xi are known vector functions of the covariate vector xi=xi1,…,xim′, βp and βq* are vectors of regression coefficient and θp and θq* are the marginal parameters. We use the Gaussian copula function to jointly analyze these types of responses. This function has been used in a variety of applications (Song [32]) because of its flexibility and analytical tractability, and to extend the regression models for a bivariate mixed outcomes of De Leon and Wu [29] to the multivariate mixed outcomes. As it was shown in recent paper Jiryaie et al. [35] the joint CDF and density of Yi can be obtained as follows:
PYi1≤yi1,…,YiP≤yiP,Zi1=k1,…,ZiQ=kQ=PYi1≤yi1,…,YiP≤yiP,γ1k1<Yi1*≤γ1k1+1,…,γQkQ<YiQ*≤γQkQ+1;R =∑ε1=01…∑εQ=01−1Q+∑j=1QεjΦP*t1,…,tP,s1k1+ε1,…,sQkQ+εQ;R
and
fYiyi;θ=∑ε1=01…∑εQ=01−1Q+∑j=1Qεj∂P∂yP…∂y1ΦP*t1,…,tP,s1k1+ε1,…,sQkQ+εQ;R,
respectively, where tp=Φ−1FYipyip;μp,θp, p=1,…,P, sqkq+εq=Φ−1FYiq*γqkq+εq;μq*,θq*, q=1,…,Q, P*=P+Q, θ is a vector parameter that contains βp, βq*, γq, θp, θq* and R, Φ−1. is the inverse function of the univariate standard normal CDF, while ΦP*.;R is the P*-dimensional standard CDF with correlation matrix Rrij and yi=yi1,…,yiP,k1,…,kQ′.

2.1. Estimation

The maximum likelihood estimator (MLE) of θ can be found by maximizing the log-likelihood function ℓ(θ)=∑i=1nlogfYiyi;θ, with an iterative technique such as the Newton–Raphson updating scheme. Hence we should layout the density (1), let fYPyp be the density of Yp, Sq=Φ−1FYq*Yq*, Tp=Φ−1FYpYp, and R|T1:Tk=RTk+1:SQ|T1:Tk be the partial correlation matrix for, Tk+1,…,TP,S1,…,SQ, after eliminating T1,…,Tk, for k=1,…,P−1. Note that

ΦP*t1,…,tp,s1k1+ε1,…,sQkQ+εQ;R=∫−∞t1ϕtΦP*−1t2|1,…,tP|1,s1|1k1+ε1,…,sQ|1kQ+εQ;R|T1dt,
where ϕ is a standard normal density. Since ∂t1/∂y1=fY1y1/ϕt1 the differentiation of (2) with respect to y1 is
ΦP*−1t2|1,…,tP|1,s1|1k1+ε1,…,sQ|1kQ+εQ;R|T1fY1y1.

Next with applying (2), differentiation of function (3) with respect to y2 is

ΦP*−2t3|1:2,…,tP|1:2,s1|1:2k1+ε1,…,sQ|1:2kQ+εQ;R|T1:T2×ϕt2|11−rT1T22fY2y2ϕt2×fY1y1,
where ∂t2|1/∂y2=ϕt2|1∂t2/∂y2; so by repeated application of (2) we can find the densities of Y. We remark that notice the elements of R|T1:Tk can be computed recursively. For example,
rTpSq|T1:Tk=rTpSq|T1:Tk−1−rTkTp|T1:Tk−1rTkSq|T1:Tk−11−rTkTp|T1:Tk−121−rTkSq|T1:Tk−12,
for p=k+1,…,P, q=1,…,Q.

In situations where the dimension of the vector variable Y is elevated, some difficulties will arise to layout the density (1). This is particularly true for high dimensional continuous vector variables due to differentiating the nested conditional normal distribution. In addition, the dimension of vector parameter θ in the joint distribution of Y will increase, which leads to computational problems to estimate the parameters. Fortunately for every dimension of Y all pair marginal density can be easily found from (1) as follows:

fYp1,Yp2yp1,yp2;μp1x,μp2x,θp1,θp2=fYp1yp1fYp2yp2×|R*|−12exp12tp1,tp2′I2−R* −1tp1,tp2,fYp,Zqyp,kq;μpx,μq*x,θp,θq*=fYpypDrqpsqkq+1,tp−Drqpsqkq,tp,fZq1,Zq2kq1,kq2;μq1*x,μq2*x,θq1*,θq2*=∑εq1=01∑εq2=01−1εq1+εq2ΦPsq1k1+ε1,sq2k2+ε2.

Herein q1≠q2, p1≠p2, q1,q2,q=1,…,Q and p1,p2,p=1,…,P, R*=1rp1p2rp2p1 1 and Dra,b=Φa−rb1−r2.

In order to overcome these complicated problems, we can use the pairwise likelihood estimation procedure of Fieuws and Verbeke [30]. In this approach instead of maximizing the full log-likelihood, each pairwise log-likelihood is separately maximized. Let the vector parameter of all possible pair likelihoods be ψ=ψ1 ′,…,ψm ′′, ψj represents the vector of all parameters in the jth bivariate joint model, j=1,…,m, m=P*P*−1/2 is the total number of possible pairs. Note the vector parameter θ of density (2) and ψ are not equivalent, elements of correlation matrix R in θ have a single counterpart in ψ, while the marginal parameters in θ have P*−1 counterparts in ψ. Thus

θ^=Aψ^,
where A is a matrix containing the appropriate coefficients to calculate the averages and ψ^=ψ^1 ′,…,ψ^m ′′ which ψ^j, j=1,…,m is obtained from maximizing the jth pair log-likelihood, ℓjψj=∑i=1nlog fYijyij;ψj, where fYijyij;ψj is the jth pair joint density with yi1=yi1,yi2,…,ym=yiP−1,yiP. It can be shown that θ^−θ∼N0,Σ, with Σ=AJ−1KJ−1A′, where J is a block-diagonal matrix with diagonal blocks Jtt and K is a symmetric matrix containing blocks Ktr, as Jtt=∑i=1n∂ℓitψt∂ψt∂ℓitψt∂ψt′ and Ktr=∑i=1n∂lit(ψt)∂ψt ∂lir(ψr)∂ψr ′, t,r=1,…,m.

For computational convenience in the estimation step the constraints on rij∈−1,1 can be removed with the Fisher's Z-transformation rij*=12log1+rij1−rij. So by the delta method, we have SEr^ij=1+r^ij1−r^ijSEr^ij*, SE is the standard error.

Based on Wald [36] an approximate 1001−α% confidence region for β=β1′,…,βp′,β1†′,…,βQ†′′ with βq†=γqkq,βq*′′, q=1,…,Q is given by

C=β:β^−β′varβ^−1β^−β≤χα,p~2,
here p~ is the length of the vector parameter β.

In the ordinal-continuous multi-response system the goal of simultaneous optimization is determining a point, xo, in the design region, R, at which the estimated mean responses of continuous variables, μ^pxo=hpδp′xoβ^p, and cumulative probabilities of desired categories kq=kqo,q=1,…,Q, p^Zq≤kqo=FYq*δ  q ′†(xo)β^q†;0,θ^q*, δ  q†(x)=1,−δ  q ′*(x) ′, are optimal. Let μ(x)=μ1(x),…μp(x),pZ1≤k10,…,pZQ≤kQ0′, ηpx=δp′xβp and ηqo(x)=δ  q ′†(x)βq† a first-order Taylor expansion of μ^px and p^Zq≤kqo around ηpx and ηqox, respectively are

μ^px=μpx+∂hpηpx∂ηpxη^px−ηpx
and
p^Zq≤kqo=pZq≤kqo+∂FYq*ηqox;0,θq*∂ηqoxη^qox−ηqox.

Thus an approximation of varμ^x is given by

varμ^x=ABvarβ^B′A′.

Herein A=diag∂μ1x∂η1x,…,∂μPx∂ηPx,∂FY1*η1ox;0,θ1∂η1ox,…,∂FYQ*ηQox;0,θQ∂ηQox, B is a block-diagonal matrix with diagonal blocks δpx, p=1,…,P and δq†x, q=1,…,Q. So the estimation of (6) is

var^μ^x=A^Bvarβ^B′A^′,
where A^ is the MLE of A.

3. THE SIMULTANEOUS OPTIMIZATION PROCEDURE

At the outset, we individually optimize each estimated mean response of continuous variables, μ^px, p=1,…,P and cumulative probabilities of desired categories FYq*δ  q ′†(xo)β^q†;0,θ^q*, q=1,…,Q over the experimental region R. We denote the individual optimum value of μ^px, p=1,…,P and FYq*δ  q ′†(xo)β^q†;0,θ^q*, q=1,…,Q by κ^j, j=1,…,P*. If all the individual optima κ^1,…,κ^P* are attained at the same setting of x in R, then the optimization problem will be solved and no further study will be required; otherwise, it is necessary to search compromise conditions favorable for all responses. To access these compromise conditions, we follow the generalized distance approach by Khuri and Conlon [14] that Mukhopadhyay and Khuri [15] adapted for the multivariate GLM situation. This approach attempts to find the conditions on x that minimizes the distance between the estimated mean responses vector μ^x and the individual optima vector κ^x=κ^1x,…,κ^p*x′. Such a distance function is denoted by ρμ^x,κ^. A variety of choices is possible for the distance measure ρ, Khuri and Conlon [14] have suggested the two following distance functions:

ρ1μ^x,κ^=μ^x−κ^′var^μ^x−1μ^x−κ^12
and
ρ2μ^x,κ^=∑j=1mμ^jx−κ^j2κ^j212.

Since κ^1,…,κ^P* are the individual optimum values of the random variables μ^px, p=1,…,P and FYq*δ    q ′†(xo)β^q†;0,θ^q*, q=1,…,Q, they are random variables themselves. We need to incorporate the variability of κ^ into ρ. For this purpose, suppose ξjβ to be the true optimum value of the jth j=1,…,P* mean response or cumulative probability optimized individually over R, and let ξβ=ξ1β,…,ξP*β′. Minimizing ρμ^x,ξβ is impossible because ξβ is unknown. We instead minimize an upper bound of ρμ^x,ξβ to carry out the optimization step, as in Khuri and Conlon [4]. Let Dξ be a confidence region for the true optima vector ξβ. By employing the 1001−α% confidence region C on β given in (5), Mukhopadhyay and Khuri [15] show that an appropriate choice of Dξ is given by

Pξβ∈×j=1P*DjC=Pξjβ∈DjC|j=1,…,P*≥Pβ∈C≈1−α,
where × denotes the Cartesian product and DjC=minγ∈Cξjγ,maxγ∈Cξjγ is an individual confidence interval of ξjβ.

If ξβ∈Dξ, then we have

ρμ^x,ξβ≤max ξ∈Dξ ρμ^x,ξ,
and hence
min x∈R ρμ^x,ξβ≤min x∈R max ξ∈Dξ ρμ^x,ξ.

Therefore by minimizing the right-hand side of (10) over the region R, we adopt a conservative distance approach to the optimization problem.

4. SIMULATION STUDY: FULL VERSUS PAIRWISE LIKELIHOOD ESTIMATION

In order to evaluate and compare the performance of the estimators from the pairwise and full likelihood approaches for the mixed ordinal-continuous responses, we have considered a simple 3×3 design with xi=xi1,xi2′, xij∈−1,0,1, i=1,…,9 and j=1,2. Distribution of the responses have the following form:

Yi1∼Gamma μ1xi,ν,   Yi2∼normal μ2xi,σ2,
where μ1xi=exp β10+β11xi1+β12xi2, μ2xi=β20+β21xi1+β22xi2 and the density of the Gamma μ1xi,ν distribution is specified as
fYi1yi1;μ1xi,ν=1Γννμ1xiνyi1ν−1e−νμ1xiyi1,
with EYi1=μ1xi and CDF FYi1.;μ1xi,ν,
Zi1=0Yi1*+μ1*xi≤γ11Yi1*+μ1*xi>γ1,Zi2=0Yi2*+μ2*xi≤γ211γ21<Yi2*+μ2*xi≤γ222Yi2*+μ2*xi>γ22,
where Y1*,Y2*∼normal 0,1, μ1*xi=β11*xi1+β12*xi2, μ2*xi=β21*xi1+β22*xi2. Following (1) the joint density of Yi=Yi1,Yi2,Zi1,Zi2′ can be simplified as
fYiyi;θ=fYi1yi1;μ1xi,νϕt2|1σ1−rT1,T22φyi1,yi2,k1,k2,
herein φyi1,yi2,k1,k2 is evaluated for yi1≥0 and yi2∈R at k1=0,1 and k2=0,1,2 in Table 1 with Φ2⋅;r* the standard bivariate normal CDF with correlation r* and t2=yi2−μ2xi/σ.
t1=Φ−1Fyi1;μ1xi,ν,t2|1=t2−rT1T2t11−rT1T22,
s1|1=γ1−μ1*xi−rT1S1t11−rT1S12,s1|1:2=s1|1−rT2S1|T1t2|11−rT2S1|T12,
s2|1j=γ2j−μ2*xi−rT1S2t11−rT1S22,s2|1:2j=s2|1j−rT2S2|T1t2|11−rT2S2|T12,
rT2Sj|T1=rT2Sj−rT1T2rT1Sj1−rT1T221−rT1Sj2,r*=rS1S2|T1−rT2S1|T1rT2S2|T11−rT2S1|T121−rT2S2|T12,
for j=1,2. A total of R = 1000 repeated samples of size ni=10, i=1,…,9 for each run of experiments were generated with the true values of parameters given in Table 2. Estimation of parameters and calculation of the likelihood score functions at estimated parameters were implemented in R using the “optim” and “fdHess” functions, respectively. The results of this simulation study are reported in Table 2. Relative biases and mean square errors of the pairwise and full likelihood suggest that the pairwise likelihood estimation obtains suitable point estimates with small mean square errors. Furthermore relative efficiencies that are generally close to 1 in the pairwise method show that pairwise likelihood estimations have standard errors which reflect these estimates true sampling variability.

k1=0 k1=1
k2=0 Φ2s1|1:2,s2|1:21;r* Φ2s1|1:2,s2|1:22;r*−Φ2s1|1:2,s2|1:21;r*
k2=1 Φ2s1|1:2,s2|1:22;r*−Φ2s1|1:2,s2|1:21;r* Φs2|1:22−Φs2|1:21−Φ2s1|1:2,s2|1:22;r*+Φ2s1|1:2,s2|1:21;r*
k2=2 Φs1|1:2−Φ2s1|1:2,s2|1:22;r* 1−Φs2|1:22−Φs1|1:2+Φ2s1|1:2,s2|1:22;r*
Table 1

φyi1,yi2,k1,k2 evaluated at k1=0,1, and k2=0,1,2.

Parameter Truth RBp×100 RBf×100 Ave SEp SDp REp MSEp MSEf
β10 1 −0.8577 −25.0141 0.1035 0.1107 0.9354 0.0123 0.0798
β11 −3 −0.1182 −0.2141 0.1260 0.1559 0.8080 0.0243 0.0254
β12 −2 −0.2783 −1.1402 0.1259 0.1374 0.9164 0.0189 0.0217
β20 −1 3.6322 16.4251 0.1867 0.1641 1.1380 0.0282 0.0608
β21 1 0.7144 −0.1973 0.2410 0.2026 1.1893 0.0411 0.0451
β22 −1 −1.4872 −2.2351 0.2396 0.2007 1.1939 0.0405 0.0454
γ1 1 11.4930 −29.3023 0.3426 0.3477 0.9854 0.1272 0.0890
β11* −1 10.1537 −10.6593 0.3854 0.3956 0.9741 0.1513 0.0567
β12* 2 10.6494 −14.6922 0.5231 0.4963 1.0539 0.3493 0.0919
γ11 −0.25 1.6823 116.7175 0.2153 0.1948 1.1054 0.0380 0.0880
γ12 1 6.0361 −2.0552 0.2547 0.2356 1.0808 0.0592 0.0396
β21* 1 6.4798 −15.4204 0.2652 0.2426 1.0930 0.0631 0.0521
β22* −2 5.1658 −14.7951 0.3553 0.3184 1.1161 0.1120 0.0920
σ 1 −4.3821 29.3046 0.0753 0.1196 0.6298 0.0162 0.0890
ν 0.4 3.5435 −20.1926 0.0534 0.0521 1.0258 0.0029 0.0090
rT1T2 0.5 −1.5025 18.1566 0.0918 0.0851 1.0793 0.0073 0.0098
rT1S1 0.5 4.4983 −40.5024 0.1917 0.1621 1.1827 0.0268 0.0424
rT1S2 0.5 2.7490 −40.6943 0.1427 0.1278 1.1160 0.0165 0.0428
rT2S1 0.5 5.1281 −31.6022 0.2037 0.1746 1.1671 0.0311 0.0331
rT2S2 0.5 3.3858 −33.5907 0.1484 0.1312 1.1306 0.0175 0.0340
rS1S2 0.5 6.8221 18.7066 0.3100 0.2837 1.0927 0.0817 0.0099

True values are given, the relative bias (RB = Bias ÷ Parameter) under the pairwise RBp and full likelihood RBf, average estimate of standard errors Ave SEp, empirical standard division SDp and relative efficiency REp=Ave SEp÷SDp under the pairwise approach are obtained. Mean squared error under pairwise MSEp and full likelihood MSEf are listed.

Table 2

Result of simulation study for pairwise and full likelihood estimation.

5. ILLUSTRATIVE EXAMPLE

In dealing with simultaneous optimization of mixed ordered categorical and continuous responses, a case study of an ion implantation process from a Taiwanese integrated circuit (IC) fabrication manufacturer was conducted by Hsieh and Tong [27] based on artificial neural networks. This example contained two quality responses: i) the amount of ion implanted in a wafer, continuous response denoted by Y, and ii) the defect situation of a sensitive area in the wafer, an ordered response denoted by Z, which included five ordered categories: very good, good, not good not bad, bad, and very bad. The responses was listed in a progressively worse order; we denoted these categories as 1–5. Each wafer had 36 sensitive areas that were tested independently. There were six control factors denoted by X1,…,X6. Between these, X1 was discrete whereas the others were continuous. Table 3 lists the control factors with their levels and coded levels of X1,…,X6 denoted by x1,…,x6. Therefore the region of the experiment transformed to the R=x=x1,…,x6′|x1=0,1;1≤xl≤3,l=2,…6. The two mixed responses data in a L18 orthogonal array are given in Table 4. In this table mik, i=1,…,18, k=1,…5 are the number of 36 sensitive areas in the ith wafer, which was tested in the ith run (ith level of x) of the experiment which fell into the kth category.

Level X1 X2 X3 X4 X5 X6 x1 x2 x3 x4 x5 x6
Level 1 Type 1 6 50 5 4 25 1 1 1 1 1 1
Level 2 Type 2 12 100 10 8 50 0 2 2 2 2 2
Level 3 18 150 15 12 75 3 3 3 3 3
Table 3

Control factors and coded factors with their levels.

i yi1 mi1 mi2 mi3 mi4 mi5
1 741.4 33 3 0 0 0
2 972.1 24 5 6 1 0
3 796.1 6 2 20 8 0
4 797.8 0 28 4 4 0
5 796.6 2 2 4 12 16
6 802.1 4 0 20 4 8
7 908.2 0 2 6 14 14
8 645.7 10 2 8 4 12
9 650.3 0 0 0 24 12
10 1072.5 34 0 2 0 0
11 1316.1 30 2 4 0 0
12 890.5 10 10 12 0 4
13 886.6 14 8 10 4 0
14 826.5 8 16 12 0 0
15 800.1 0 8 6 4 18
16 816.1 18 12 6 0 0
17 824.2 10 6 0 4 16
18 735.6 0 4 2 6 24
Table 4

Experimental data.

The continuous response is a nominal-the-best (NTB) with a target value of 1000 (after the data was transformed). First, for the continuous response Y we fitted the normal GLM regression with link μxi=δ′xiβ and gamma GLM regression with link μxi=exp δ′xiβ where the density of the Gamma μi,ν distribution was specified as (11), δxi=1,xi1,…,xi6′ and β=β0,β1,…,β6′. Akaike information criteria (AIC) were 230.37 and 225.1, respectively. Gamma distribution appeared to have a better fit to this response compared to normal distribution. We assumed that the dispersion parameter ν was known and did not need to be estimated in the joint regression model. In our example we used its MLE in the marginal gamma GLM, which was 35.9348.

Zi, i=1,…,n, n=18 has the ordinal form with K = 5 categories. Therefore we have used the cumulative probabilities. The kth cumulative probability in the ith run of the experiment is

pZi≤k=FY*γk−δ*xi′β*=expγk−δ*xi′β*1+expγk−δ*xi′β*,
where k=0,…,K−1, δ*=xi1,…,xi6′, β*=β1*,…,β6*′, and γ=γ1,…γ4′.

Yi=Yi,Zi′ given covariate xi=xi1,…,xi6′, i=1,…,18 are independent, Yi∈R+ and Zi∈0,…,4, using (1), we can write the joint density of Yi and Zi as follows:

fYi,Ziyi,k;θ=fYiyi;μxiDrsk+1,t−Drsk,t,
where t=Φ−1FYiyi;μxi,ν, sk=Φ−1FYi*γki−μ*xi, Dra,b=Φa−rb1−r2 and θ=β′,β′*,γ′,r′.

Therefore log-likelihood for the presented data set in Table 4 is

ℓθ=∑i=1nlog fYiyi;μxi+∑i=1n∑k=0K−1mik+1logDrsk+1,t−Drsk,t.

Following this log-likelihood, we estimated the parameters and calculated the likelihood's score functions at the estimated parameters in R software using the “optim” and “fdHess” functions, respectively. Table 5 lists the estimated parameters, their standard errors and p-values. In this table r^=0.0092 with a p-value = 0.4672 shows we can accept that the amount of ion implanted and the defect situation of a sensitive area in a wafer are independent. In this example we aim to reach a point in the design region that simultaneously minimize |μ^x−1000| and maximize π^1x as cumulative probabilities of the desired category, herein π^kx=p^Z=k−1, k=1,…,5. Table 6 shows κ^1 and κ^2, the individual optimum of μ^x and π^1x over R, do not have the same location in the design region. So its needed to continue, we know the components of κ^x are random variables, using these variablity in κ^ a confidence region such as (12) can be constructed, before using the distance metric. For this purpose 5000 values were randomly selected from C, the confidence region of β′,β*′′, and individual optimum values, μx,π1x′, of each point over the experimental region R were computed. We denote these optimum by ξ1,ξ2′ and we find the minimum and maximum 5000 values of ξ1,ξ2′ to get the approximate lower and upper bound of DiC, i = 1, 2, where in Table 6 D1C=822.0784,1000 and D2C=0.8486,0.9648 are confidence intervals with a 95% convergence probability for μ and π1, respectively.

Parameter Estimate Std.Error z=θ^seθ^ p-value
β^0 7.1070 0.2268 31.3394 0.0000
β^1 −0.1316 0.0789 −1.6682 0.0476
β^2 −0.1143 0.0481 −2.3778 0.0087
β^3 −0.0557 0.0488 −1.1416 0.1268
β^4 −0.0115 0.0474 −0.2426 0.4042
β^5 −0.0281 0.0487 −0.5772 0.2819
β^6 0.0566 0.0488 1.1610 0.1228
β^1* 0.7580 0.1561 4.8551 0.0000
β^2* 1.4816 0.1073 13.8052 0.0000
β^3* 1.1839 0.1009 11.7374 0.0000
β^4* −0.2319 0.0939 −2.4699 0.0068
β^5* 0.1324 0.0970 1.3653 0.0861
β^6* 0.3071 0.0942 3.2607 0.0006
γ^1 4.9282 0.4683 10.5234 0.0000
γ^2 5.9847 0.4880 12.2635 0.0000
γ^3 7.2068 0.5166 13.9513 0.0000
γ^4 8.2651 0.5351 15.4455 0.0000
r^ 0.0092 0.1125 0.0822 0.4672
Table 5

Estimates and standard errors.

Gamma Response Ordinal Response
Location (1, 6.002, 63.848, 6.08, 4.401, 71.88) (0, 6, 50, 15, 4, 25)
Optimum mean response 1000 0.925
Confidence region (822.078, 1000) (0.849, 0.965)
Table 6

The individual optima and confidence region Dξ.

Taghuchi Neural Network Gaussian Copula
Location (1, 12, 50, 15, 8, 25) (1, 6.06, 46.32, 12.19, 11.63, 52.06) (0,7.26, 54.9, 12.55, 11.91, 25.42)
Ordinal response (0.54, 0.23, 0.15, 0.05, 0.03) (0.75, 0.14, 0.07, 0.02, 0.01) (0.85,0.094, 0.04, 0.01, 0.01)
Gamma response 778.113 912.682 946.277
Maxρ2 0.496 0.245 0.194
Table 7

Simultaneous optima (Example 1).

In this example we used the distance measure ρ2 due to two responses that had an equal importance in the ion implantation process by Hsieh and Tong [27] and were independent. For each of the 5000 values of x randomly selected from the region R, we computed the maximum distance function ρ2μ^x,ξ with respect to ξ∈Dξ, where Dξ=×j=12DjC, and μ^x=μ^x,π^1x′. Next a minimum of 5000 values of this maximum distance were obtained. This distance, the corresponding simultaneous maximum of μ^x and π^1x and their locations are given in Table 7.

For the purpose of compression we considered two points of the design region that Hsieh and Tong [27] introduced for simultaneous optimization of these responses with the Taguchi and Artificial Neural Network methods. By using these locations, estimated parameters in Table 5 and confidence region in Table 6, μ^, π^1,…,π^5′ and maximum of ρ2μ^x,ξ with respect to this confidence region were computed (Table 7). These results showed that the location founded by the Gaussian copula had better optimum values for these responses.

6. CONCLUSION

In the simultaneous optimization problem, due to the inherent nature of the data and convenience of measurements, it is not feasible to report all of the responses as continuous variables with normal distribution. One of the most popular means is to represent the data in the ordinal categorical form. Thus the outputs may involve mixed continuous and ordinal variables. The innovative use of the copula function permits a model of various types of correlated responses, such as mixed continuous and ordinal responses, those with all ordinal categorical forms, continuous responses that have different marginal distributions, or where standard multivariate distribution of the responses is not applicable or does not exist. This paper used the pairwise likelihood estimation method for a high dimension of responses, and alleviated the computational demands of estimation. The results of the simulation study showed the usefulness of this method. Adopting the generalized distance approach would allow us to simultaneously optimize such responses by considering the dependency between them. An example demonstrated the effectiveness of the proposed method. The published methods could not be directly applied to simultaneous optimization of such correlated responses.

REFERENCES

2.A.E. Hoerl, Chem. Eng. Prog. Symposium Ser., Vol. 60, 1964, pp. 67-77.
3.G.S. Peace, Taguchi Methods: A Hands-on Approach, Addison-Wesley, Reading, 1993.
4.W.Y. Fowlkes and C.M. Creveling, Engineering Methods for Robust Product Design: Using Taguchi Methods in Technology and Product Development, Addison-Wesley, Reading, 1995.
6.E.C. Harrington, Ind. Qual. Control, Vol. 21, 1965, pp. 494-498.
21.G. Taguchi, Taguchi Methods: Case Studies from the U.S. and Europe, American Supplier Institute, Michigan, 1991.
22.G. Taguchi, Taguchi Methods: Research and Development, American Supplier Institute, Michigan, 1991.
23.G. Taguchi, Taguchi Methods: Signal-to-Noise Ratio for Quality Evaluation, American Supplier Institute, Michigan, 1991.
28.F.C. Wu, Int. J. Ind. Eng., Vol. 15, 2008, pp. 231-238.
31.A. Sklar, Fonctions de Répartition à n Dimensions et Leurs Marges, Publications de l'Institut Statistique de l'Université de Paris, Paris, Vol. 8, 1959, pp. 229-231.
32.P.X.K. Song, Correlated Data Analysis: Modeling, Analytics, and Applications, Springer, New York, 2007.
<Previous Article In Issue
Download article (PDF)
Next Article In Issue>
Journal
Journal of Statistical Theory and Applications
Volume-Issue
18 - 3
Pages
212 - 221
Publication Date
2019/09/09
ISSN (Online)
2214-1766
ISSN (Print)
1538-7887
DOI
10.2991/jsta.d.190701.001How to use a DOI?
Copyright
© 2019 The Authors. Published by Atlantis Press SARL.
Open Access
This is an open access article distributed under the CC BY-NC 4.0 license (http://creativecommons.org/licenses/by-nc/4.0/).

Cite this article

risenwbib
TY  - JOUR
AU  - Fatemeh Jiryaie
AU  - Ahmad Khodadadi
PY  - 2019
DA  - 2019/09/09
TI  - Simultaneous Optimization of Multiple Responses That Involve Correlated Continuous and Ordinal Responses According to the Gaussian Copula Models
JO  - Journal of Statistical Theory and Applications
SP  - 212
EP  - 221
VL  - 18
IS  - 3
SN  - 2214-1766
UR  - https://doi.org/10.2991/jsta.d.190701.001
DO  - 10.2991/jsta.d.190701.001
ID  - Jiryaie2019
ER  -