The problem of optimizing several response variables

When we perform an experiment in the laboratory it is usual to quantify more than one response variable in each of our experimental units. Depending on our objectives we might be interested in maximizing some responses while the rest would only be measured in order to characterize the process. But what if what we are looking for is that each response reaches a maximum, a minimum, or stays within a defined limit?

A feasible solution is, once we fit the models for each variable, to use the desirability function method to obtain a global optimum with all response variables.

This method consists of defining a function that covers the entire experimental region and estimates the global desirability (GD). In turn, the GD is defined taking into account all the responses, which allows to reduce a multivariate optimization problem to a single variable one. The objective becomes to maximize GD.

To define GD, first transform each predicted response ypi(x) into an individual desirability value di(x) that falls in the interval [0, 1]. The transformation di(x) of each response is performed according to the target values established by the researcher and measures the desirability at the experimental point or treatment x = (x1, x2, x3, …, xk), where each “x” corresponds to each of the experimental factors considered in the design.

If the variable ypi has an upper specification (USi), a lower specification (LEi) and its target value is Ti, the transformation di(x) is defined as:

The exponents s and t are used to choose the form of transformation. If large values are assigned fir these exponents, greater than or equal to 5, for example, desirability values close to 1 will be obtained only when ypi(x) is sufficiently close to the target value. On the other hand, selecting small values of s and t (≤ 0.10) will result in desirability values close to 1 for a wide range of values within the interval [LUi, USi]. If both exponents are selected as equal to 1, a linear increase in desirability will be obtained as the ypi(x) values approach the target value.

Once the m desirabilities have been calculated for the m responses at the point within the experimental region x, the GD at x is defined by a weighted geometric mean:

In this case the weights w are constants that are assigned to balance the relative importance of each response variable compared to the others. The larger the value of w the greater the requirement for the GD result to benefit that response. If we select a value of 1 for each w, the above equation is reduced to:

The global optimum experimental point x0 is the point where the value of GD(x0) is maximum. To find this maximum value we applied a numerical method.

With this brief explanation of the method, it is now time to perform a practical example with R. For this I will use the data published in George Derringer & Ronald Suich (1980).

1. Importing data

The experimental data are stored in a CSV file in the usual structure of a design matrix.To import the data I used the read_csv() function of the readr package.

library(readr)

data_tires <- read_csv("data/data_tires.csv")

head(data_tires)
## # A tibble: 6 × 7
##      x1    x2    x3    Y1    Y2    Y3    Y4
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1    -1    -1     1   102   900   470  67.5
## 2     1    -1    -1   120   860   410  65  
## 3    -1     1    -1   117   800   570  77.5
## 4     1     1     1   198  2294   240  74.5
## 5    -1    -1    -1   103   490   640  62.5
## 6     1    -1     1   132  1289   270  67

The data correspond to an experiment where the aim is to improve the formulation of belts for tire manufacturing, in such a way that certain quality parameters are met. The x’s correspond to the factors and the Y’s correspond to the responses.

Column information is as follows:

• x1: hydrated silica level
• x2: silane coupling agent level
• x3: sulfur level
• Y1: abrasion index PICO
• Y2: 200% module
• Y3: elongation at break
• Y4: hardness

The values of each response variable must comply with the following restrictions:

• Y1 > 120
• Y2 > 1000
• 400 < Y3 < 600
• 60 < Y4 < 75

2. Fitting second order models

We fit second-order models for each response with the help of the rsm package.

library(rsm)

y1_m  <- rsm(Y1 ~ SO(x1, x2, x3), data = data_tires)
y2_m  <- rsm(Y2 ~ SO(x1, x2, x3), data = data_tires)
y3_m  <- rsm(Y3 ~ SO(x1, x2, x3), data = data_tires)
y4_m  <- rsm(Y4 ~ SO(x1, x2, x3), data = data_tires)

The models are as follows:

To get an idea of the variability explained by each model I defined a small function that returns the coefficient of determination and adjusted coefficient of determination.

# Get a R^2 and adjusted R^2 from a linear model
get_r2 <- function(model){
  sum_model <- summary(model)
  round(
    c(R2 = sum_model$r.squared, adj.R2 = sum_model$adj.r.squared),
    3
  )
}

# Table with R-squared for each model
y1_r2 <- get_r2(y1_m)
y2_r2 <- get_r2(y2_m)
y3_r2 <- get_r2(y3_m) 
y4_r2 <- get_r2(y4_m)

r2t <- t(data.frame(y1_m = y1_r2, y2_m = y2_r2, y3_m = y3_r2, y4_m = y4_r2))
r2t
##         R2 adj.R2
## y1_m 0.972  0.947
## y2_m 0.742  0.510
## y3_m 0.981  0.965
## y4_m 0.958  0.921

As seen in the table, the coefficients of determination for the Y2 model are lower compared to the rest. This should be taken into account as there could be another model that fits the Y2 data better. However, for learning purposes I used all quadratic models as done in the previously mentioned paper.

3. Simultaneous optimization

# Prediction function
rsm_opt <- function(x, dObject, space = "square"){
  
  df <- data.frame(x1 = x[1], x2 = x[2], x3 = x[3])
  
  y1 <- predict.lm(y1_m, df)
  y2 <- predict.lm(y2_m, df)
  y3 <- predict.lm(y3_m, df)
  y4 <- predict.lm(y4_m, df)
  
  out <- predict(dObject, data.frame(y1 = y1, y2 = y2, y3 = y3, y4 = y4))
  
  if(space == "circular" & sqrt(sum(x^2)) > 1.63) out <- 0
  else if(space == "square" & any(abs(x) > 1.63)) out <- 0 
  
  out
}

# Optimization function
maximize_overall <- function(int_1 = c(-1.63, 1.63),
                             int_2 = c(-1.63, 1.63),
                             int_3 = c(-1.63, 1.63),
                             dObject = NULL, 
                             space = "square"){
  
  searchGrid <- expand.grid(
    x1 = seq(int_1[1], int_1[2], length.out = 5),
    x2 = seq(int_2[1], int_2[2], length.out = 5),
    x3 = seq(int_3[1], int_3[2], length.out = 5)
  )
  
  for(i in 1:dim(searchGrid)[1]){
    
    tmp <- optim(as.vector(searchGrid[i,]),
                 rsm_opt,
                 dObject = dObject,
                 space = space, 
                 control = list(fnscale = -1))
    
    if(i == 1) best <- tmp
    if(tmp$value > best$value) best <- tmp
  }
  
  best
}

library(desirability)

D_y1 <- dMax(120, 170)
D_y2 <- dMax(1000, 1300)
D_y3 <- dTarget(400, 500, 600)
D_y4 <- dTarget(60, 67.5, 75)

D_overall <- dOverall(D_y1, D_y2, D_y3, D_y4)

overall_opt <- maximize_overall(dObject = D_overall)

overall_opt
## $par
##          x1          x2          x3 
## -0.05345894  0.14718876 -0.86635592 
## 
## $value
## [1] 0.5833527
## 
## $counts
## function gradient 
##      502       NA 
## 
## $convergence
## [1] 1
## 
## $message
## NULL

4. Obtaining predictions of each response in the global optimum.

With the optimal experimental point we predict each response.

# Optimal point as a data frame
data_opt <- data.frame(
  x1 = overall_opt$par[1],
  x2 = overall_opt$par[2],
  x3 = overall_opt$par[3]
  )

# Predict each response at overall optimization point
y1_opt <- predict.lm(y1_m, data_opt)
y2_opt <- predict.lm(y2_m, data_opt)
y3_opt <- predict.lm(y3_m, data_opt)
y4_opt <- predict.lm(y4_m, data_opt)

# Response predictions in a data frame
res_opt <- data.frame(y1_opt, y2_opt, y3_opt, y4_opt)
rownames(res_opt) <- "Optimal responses"
round(res_opt, 2)
##                   y1_opt y2_opt y3_opt y4_opt
## Optimal responses 129.43   1300 465.97  68.02

It may also be useful to obtain the individual desirabilities for each response.

# Optimal desirability values 
d1_opt <- predict(D_y1, y1_opt)
d2_opt <- predict(D_y2, y2_opt)
d3_opt <- predict(D_y3, y3_opt)
d4_opt <- predict(D_y4, y4_opt)

d_opts <- data.frame(d1_opt, d2_opt, d3_opt, d4_opt)
rownames(d_opts) <- "Optimal desiabilities"
round(d_opts, 3)
##                       d1_opt d2_opt d3_opt d4_opt
## Optimal desiabilities  0.189      1   0.66  0.931

All these values accomplish with the restrictions established at the beginning.

It should be considered whether it is possible to set the level of each factor to the optimal obtained by the simultaneous optimization. In addition, all predictions must be contrasted experimentally.

If it is not possible to carry out the optimal treatment, an alternative could be to evaluate an experimental point that is possible to fix in practice and that is close to the optimum. In these cases it would be useful to know how global desirability behaves within our experimental region.

5. Visualization of the desirability function

To perform the desirability function visualizations I also defined a couple of functions.

# Generate a matrix with desirability predictions within the experimental region
d_matrix <- function(model_1, model_2, model_3, model_4, 
                     dObject,  c3 = -1, l_x = c(-1, 1), by = 0.1){
  
  x <- seq(l_x[1], l_x[2], by = by)
  lx <- length(x)
  
  data_x <- expand.grid(x1 = x, x2 = x, x3 = c3)
  
  y_1 <- predict(y1_m, data_x)
  y_2 <- predict(y2_m, data_x)
  y_3 <- predict(y3_m, data_x)
  y_4 <- predict(y4_m, data_x)
  
  d_m <- predict(D_overall, data.frame(Y1 = y_1, Y2 = y_2, Y3 = y_3, Y4 = y_4))
  dim(d_m) <- c(lx, lx)
  
  list(d_m = d_m, x = x)
}

# Deploys a contour plot for desirability within the experimental region
contour_d <- function(data = NULL, main = " ",  xlab = "x1", ylab = "x2"){
  filled.contour(
    z = data$d_m, x = data$x, y = data$x, 
    color.palette = colorRamps::matlab.like,
    plot.title = title(main = main, xlab = xlab, ylab = ylab, cex.lab = 1.5, 
                       cex.main = 1.5),
    plot.axes = {
      axis(1, cex.axis = 1.5)
      axis(2, cex.axis = 1.5)
    }
  )
}

Since we are dealing with three factors, to make the contour plots we set the value of x3 as a constant to values of -1, 0 and 1.

# Desirability at x3 = -1
dpx3_1 <- d_matrix(y1_m, y2_m, y3_m, y4_m, D_overall, 
                   c3 = -1, l_x = c(-1.6, 1.6))
contour_d(dpx3_1, main = "x3 = -1")

# Desirability at x3 = 0
dpx3_2 <- d_matrix(y1_m, y2_m, y3_m, y4_m, D_overall, 
                   c3 = 0, l_x = c(-1.6, 1.6))
contour_d(dpx3_2, main = "x3 = 0")

# Desirability at x3 = 1
dpx3_3 <- d_matrix(y1_m, y2_m, y3_m, y4_m, D_overall, 
                   c3 = 1, l_x = c(-1.6, 1.6))
contour_d(dpx3_3, main = "x3 = 1")

Note that all the functions I defined throughout this publication can be adapted for any number of response variables and for any number of factors. If you have any questions please send me an email at the address below.

Thank you so much for your time and for visit this site.

Juan Pablo Carreón Hidalgo

[email protected]

This work is licensed under a Creative Commons Attribution 4.0 International License.