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

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?

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)