How to generate a design for a response surface with a discrete input random variable?

Question

I'm trying to generate a design of experiments for fitting a response surface for a quantity $Y$ such that, $Y = f(X_1, ..., X_6)$. I'm open to a factorial design as well as a random design. $X_1,..., X_5$ are continuous valued with specified upper bound and lower bound values however, $X_6$ variable is discrete valued. It can only take integer values between 1 and 30.

How do I generate a design of experiments considering the fact that $X_6$ needs to be discrete valued?

I need a second order response surface i.e. there will be interaction terms of the input random variables. Therefore, a simple full factorial wouldn't fit my requirements.

For e.g., consider the case of a central composite design (CCD) for two variables, $X_1$ and $X_6$.

array([[-0.70710678, -0.70710678],
   [ 0.70710678, -0.70710678],
   [-0.70710678,  0.70710678],
   [ 0.70710678,  0.70710678],
   [ 0.        ,  0.        ],
   [ 0.        ,  0.        ],
   [ 0.        ,  0.        ],
   [ 0.        ,  0.        ],
   [-1.        ,  0.        ],
   [ 1.        ,  0.        ],
   [ 0.        , -1.        ],
   [ 0.        ,  1.        ],
   [ 0.        ,  0.        ],
   [ 0.        ,  0.        ],
   [ 0.        ,  0.        ],
   [ 0.        ,  0.        ]])

Let lower bound of $X_1$ be $2.3$ and upper bound be 2.7. Therefore, -0.707 value for $X_1$ will be 2.358 and 0.707 value will be 2.641. These numbers make sense to my model. However, I don't get how can I get the corresponding values for $X_6$ where -1 corresponds to 1 and +1 corresponds to 30.

Alternatively, can this be done using a randomized design like Latin Hypercube Sampling?

Answer 1

You can do this with a Latin hypercube sample in R:

require(lhs)

upper <- c(1,2,3,4,2)
lower <- c(0,0,1,0,1)
stopifnot(length(upper) == length(lower))
int_upper <- 30
int_lower <- 1
N <- 100

X <- randomLHS(N, length(upper) + 1)

Y <- X
for (i in seq_along(upper))
{
  Y[,i] <- qunif(X[,i], lower[i], upper[i])
}

Y[,6] <- floor(X[,6]*(int_upper + 1 - int_lower)) + int_lower

# Note Y[,6] will only be uniform if N %% (int_upper - int_lower + 1) == 0
table(Y[,6])

y <- 3 + X[,1] + 2*X[,2] + 3*X[,3] + 4*X[,4] + X[,5] + X[,6] + rnorm(N, 0, 1.5)

colnames(X) <- paste0("X", 1:6)

lm1 <- lm(y ~ (X1 + X2 + X3 + X4 + X5 + X6)^2, data = as.data.frame(X))
summary(lm1)