Optimal Design of Controlled Experiments for Personalized Decision Making in the Presence of Observational Covariates

Controlled experiments are widely applied in many areas such as clinical trials or user behavior studies in IT companies. Recently, it is popular to study the experimental design problem to facilitate personalized decision making. In this paper, we investigate the problem of optimal design of multiple treatment allocation for personalized decision making in the presence of observational covariates associated with experimental units (often, patients or users). We assume that the response of a subject assigned to a treatment follows a linear model which includes the interaction between covariates and treatments to facilitate precision decision making. We define the optimal objective as the maximum variance of estimated personalized treatment effects over different treatments and different covariates values. The optimal design can be obtained by minimizing this objective. Under a semi-definite program reformulation of the original optimization problem, we use a YALMIP and MOSEK based optimization solver to provide the optimal design. Numerical studies are provided to assess the quality of the optimal design.


Introduction
Optimal designs of experiments are developed to reduce the variance of estimated model parameters (Pukelsheim, 2006) by optimizing a function of the information matrix usually under a (generalized) linear model assumption.Examples include the determinant, the trace, and the minimum eigenvalue of the information matrix, which result in D-, A-, and Eoptimal designs, respectively.In many areas, such as clinical trials or user behavior studies in IT companies, the experiments are conducted to investigate the treatment effects over different users by designing controlled experiments, i.e., divide the users into different groups and assigning each group with a treatment or control.For example, optimal design of experiments is applied to design of clinical trials in the era of precision medicine, where the patient response may depend on some biomarkers (covariates).In fact, there is a significant growing evidence that the response of the patients may depend on their covariates in a wide range of applications.For instance, Schork (2015) showed that the top ten highest-grossing drugs in the US are only effective in 4% to 25% of the patient population, a shocking result for the current clinical practice.Biswas and López-Fidalgo (2013) applied a modified Doptimal design to early stage clinical trials when the covariates are uncontrollable.Lee and Wason (2019) applied a weighted L-optimal design for late stage Phase III clinical trials with observed covariates.As another emerging application, Tech Giants such as Google, Facebook, and LinkedIn frequently use A/B testing for marketing, web design, and datadriven decision making (Siroker and Koomen, 2013).In fact, applying controlled experiments such as A/B testing experiments in business has improved the performance of companies (Koning et al., 2022).Also, there are some methodological contributions for optimal design of experiments applied for A/B testing.For example, Bhat et al. (2020) studied an offline and online standard A/B testing with the goal of maximizing the precision of least square estimations.
The literature of optimal design of experiments for linear models, which we also use in this study, is vast.From the optimal design perspective, Atkinson (2015) investigated the problem with two treatments, and developed the D-optimal design for both covariates and treatment allocation under the assumption of a linear model with interactions between covariates and treatments.Wang and Ai (2016) generalized the work of Atkinson (2015) to experiments with multiple treatments.Following the same model setting, Rosa (2020) developed the Aand E-optimal designs for treatments and covariates.
In practice, the user covariates may not be controllable, due to the restriction in the user recruiting process.Thus, an alternative problem is to allocate the treatments to a given set of users with observed covariates.From the optimal design perspective, the aim is to obtain an exact design by optimizing a design objective with the treatment allocation as the design variable given fixed covariates, which often results in an optimization problem with binary decision variables.For two treatment allocation, under the linear model assumption without interaction between covariates and treatments, Hore et al. (2014) proposed computationally tractable algorithms to search D-and A-optimal designs.Also, with a particular interest of minimizing the variance of the estimated treatment effect (i.e., a D s optimal design criterion), Bhat et al. (2020) demonstrated that the problem can be solved efficiently based on a generalization of the MAX-CUT semi-definite programming (SDP) relaxation of Goemans and Williamson (1995).Given the linear model with interaction between treatments and covariates, Zhang et al. (2022) proposed an optimal design objective for personalized decision making with two treatments and a corresponding approximation solution approach to obtain near optimal allocations efficiently.
In this paper, we generalize the work in Zhang et al. (2022) to multiple treatment allocation for personalized decision making with observational covariates.Specifically, we consider the situation that there are a finite number of treatment options and subjects with given covariates.We assume that the response of a subject assigned to a treatment follows a linear model which includes the interaction between covariates and treatments.Accordingly, we can define the best treatment option for each subject groups, which can be estimated with data.This setup facilitates precision decision making.However, by including more than two treatments, the resulting design objective can not be simplified by the approximation approach developed for the two treatment cases in Zhang et al. (2022).Since the design objective matches the E-optimal design criterion, we reformulate the optimization problem as a SDP with binary decision variables.We use a YALMIP and MOSEK based optimization solver for SDP in a branch-and-bound scheme to provide the optimal design.In our numerical study, we assess the quality of the optimal designs provided by the solver.

Problem Description
A decision-making problem in practice can often be simplified to choose a treatment from a finite candidate set {1, . . ., K} that has the best performance measure µ k , i.e., k * ∈ argmax k∈{1,...,K} µ k . ( The values of µ k 's are unknown, which can be estimated based on collected user responses to their treatment allocations.Particularly, given n users, the experimenter allocates one of the K treatments to each of the users.The response of the i-th user is denoted by y i .
Assume that y i is a continuous scalar that follows a model given by where ε i is an additive error term, and x ik ∈ {0, 1} with x ik = 1 indicating that the k-th treatment is allocated to the i-th user.We require that K k=1 x ik = 1 to guarantee that each user is only exposed to a single treatment.Under the model assumption in (2), the performance measure µ k represents the average treatment performance over the target user population, which can be estimated by μk given responses y i 's and treatment allocation x ik 's.
Then a data-driven solution to (1) is given by k * ∈ argmax k∈{1,...,K} μk . (3) The accuracy of this decision is associated with the accuracy of the estimates μk 's.To reduce the uncertainty of the decision, we can reduce the variances of the estimates μk 's.
Under the context of personalized decision making, a user is associated with observed covari- where Z is the covariates space of the target population.
Examples of the covariates include the demographic information, social behavior and network connections.The treatment performance can often be related to the covariates values.
Let µ k (z) be a function of the users' covariates z representing the personalized treatment performance measure.In this paper, we further assume that µ k (z) is a linear function over z.Therefore, the linear model in ( 2) is represented by where β β β k = (β 1k , . . ., β pk ) ⊤ is a vector of linear coefficients.The first feature in each z i is fixed to be one as the linear intercept.Similar to Wang and Ai (2016), this model incorporates the interaction between treatment and covariates, which enables the estimation of heterogeneity of treatment performance over covariates z.
Under this model, the problem of personalized decision making becomes With estimated parameters, a data-driven solution to (5) can be obtained for any z ∈ Z,

Design Criterion of Controlled Experiments for Personalized Decision Making
Let y = (y 1 , . . ., y n ) ⊤ be the responses of n users.Correspondingly, let Z be an n × p matrix, whose rows are the users' covariates z ⊤ 1 , . . ., z ⊤ n (the first entry is loaded by one as the intercept).The treatment allocation is recorded in an n × K matrix X, where the i-th row of X is denoted by x ⊤ i = (x i1 , . . ., x iK ) ⊤ representing the treatment allocation of the i-th user.We assume that the responses are generated under the model in (4).The model covariates matrix is then given by which we assume has a full column rank.The notation ⊗ denotes the Kronecker product.
By stacking We assume that the additive error term ε i 's in (4) are independent and identically distributed random variables with mean zero and variance σ 2 .The variance-covariance matrix of the estimated b in ( 7) can be expressed by , which is the result of that x i x ⊤ i is a K ×K diagonal matrix with diagonal entries x i1 , . . ., x iK .
Then the variance-covariance matrix of b can be simplified to For any k ∈ {1, . . ., K} and z ∈ Z, we have that As noted earlier, we assume that the experiments are conducted with n existing users.Thus, the personalized information matrix Z is observed, and the experimental design problem is to determine the treatment allocation matrix X.A necessary condition for positive definiteness of the variance-covariance matrix requires that n i=1 x ik ≥ p for any k ∈ {1, . . ., K} and n ≥ Kp.
Our goal is to minimize the variance of z ⊤ β β β k for any k and z to improve the accuracy of the data-driven decision in (6).Instead of minimizing all the variances simultaneously, we formulate the optimization problem by minimizing the worst case, i.e., the maximum variance of any k ∈ {1, . . ., K} and z ∈ Z.Therefore, the design criterion is given by The optimal design is then given by solving the problem min where n i=1 x ik ≥ p is required for the positive definiteness of the matrix n i=1 z i z ⊤ i x ik , and the constraint K k=1 x ik = 1 ensures that each user is exactly assigned to one treatment.
Although the optimal design criterion has a closed-form expression, the corresponding optimization problem is still a challenge to solve directly as a mixed-integer minimax problem.
We assume that the covariates space Z is given by a unit ball, i.e., Z = {z ∈ R p | ∥z∥ ≤ 1}.
For any X that gives a positive definite variance-covariance matrix, we have that max where λ max and λ min produce the maximum and minimum eigenvalues of a matrix, respectively.Accordingly, the optimization problem in ( 10) is equivalent to max This design objective is known as the E-optimal design in the literature (e.g., Pukelsheim (2006) and Dette and Studden (1993)), i.e., min k=1,2,...,K .
Next, we provide a SDP-based solution approach to solve (11) with observed covariates Z and a binary decision X.

A SDP-based Solution Approach
Following the reformulation technique for E-optimal design in Boyd and Vandenberghe (2004), we cast problem (11) as a SDP with binary decision variables.The reformulation is based on the following observation.Let A be a real p × p symmetric matrix.It is well known in matrix analysis that for any δ ∈ R, λ min (A − δI) = λ min (A) − δ, where I is the p × p identity matrix.Therefore, λ min (A) can be expressed as the largest value of δ such that λ min (A − δI) ≥ 0, i.e., where A − δI ⪰ 0 means that A − δI is a positive semi-definite matrix.
We apply the above observation to problem (11).For k = 1, . . ., K, Consequently, problem (11) can be equivalently written as max Note that, the reformulation technique is the same as used in the SDP reformulation proposed for the approximate E-optimal design in Section 7.5 of Boyd and Vandenberghe (2004).
However, the discrete feature in problem ( 11) cannot be averted by simply considering a continuous relaxation as in Boyd and Vandenberghe (2004), even when n is large.Thus, the implementation of a SDP problem with continuous decisions cannot solve our problem.As a mixed binary SDP, the reformulation (12) can be solved by a branch-and-bound (BnB) algorithm, where each relaxation problem in a BnB node is solved by an SDP solver (Gally et al., 2018).A built-in solver in YALMIP (Löfberg, 2004) provides such a BnB implementation for MATLAB (MATLAB, 2021), and MOSEK (ApS, 2019) can be used as the embedded SDP solver.
In practice, the experimenter often requires a (roughly) balanced allocation across the K treatments.If K can be divided by n, we can include balanced constraints which also potentially accelerates the computation time of solving this problem by reducing the number of feasible solutions.If K can not be divided by n, a set of roughly balanced constraints can be added similarly.For convenience, we set n as a multiple of K and include where I{•} denotes the indicator function, and the "Relative Improvement (RI)" is defined by The two quality measures compare the objective given by the proposed solution approach with the population and the median of the objectives given by random designs.The higher values of "Percentile" and "RI" indicate better quality of the optimal design, i.e., the value of adopting the proposed optimization approach.Next, we investigate the performance of the optimal design under two common examples of fixed covariates Z.
Example I.In this example, we generate the entries of the covariates Z independently from the standard normal distribution.We consider a batch of cases with all the combinations of n ∈ {120, 240, 360}, p ∈ {3, 5} and K ∈ {2, 3, 4}.For each combination of n, p and K, we generate a set of covariates z 1 , . . ., z n , compute the quality measures Percentile and RI, and replicate this process for 100 times.The resulting 100 copies of Percentile and RI are depicted as boxplots in Figures 2 and 3, respectively.The results give the following observations: 1.The values of Percentile and RI increase with K.For K = 2, the proposed solution approach does not always give better design than the random designs.
2. The values of Percentile increase with n, whereas the values of RI decrease with n.For n = 120, the objective values given by the proposed approach may not outperform the random designs that give median or higher quantiles of objectives.
3. The values of Percentile and RI increase slightly as we increase p from 3 to 5.
These observations indicate that the cases with K = 3 or 4 or n = 240 or 360 provide solutions of higher quality than the case with K = 2 and n = 120, which is possibly due to that the difficulty in solving the problem caused by the binary constraints is more significant when K and n are small.q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q K = 2 Example II.In this example, we generate the entries of the covariates Z independently from a 0-1 Bernoulli distribution with the proportion of ones equal to r.The aim of this example is to demonstrate the impact of imbalance in covariates to the quality of optimal design.Therefore, we fix n = 200, K = 4 and p = 3, and consider r ∈ {0.1, 0.5}, where r = 0.5 gives balanced covariates and 0.1 gives imbalanced covariates.The values of Percentile and RI are depicted as boxplots in Figure 4.The results show that, for balanced covariates, the optimal design given by the proposed solution approach can often reach to the optimal value, but the values of RI are often around 20%.However, for imbalanced covariates, the relative improvement can often reach 50%-200%, but the Percentile may stand as low as 75% for some instances.This example demonstrates that the optimization problem under balanced covariates provides higher-quality solutions than that under imbalanced covariates, whereas, the benefit of optimal design under balanced covariates may not be as significant as that under imbalanced covariates.q q q q q q q q q q q q q q q q q q q q q q q q q K = 2 Relative Improvement (%) Figure 3: The quality measure RI (given by (15)) of different cases in Example I, i.e., Relative Improvement of the optimal objective λ opt with respect to the median of objectives given by 1000 random designs.
After 300 seconds, the solver outputs the optimality gap (0%-100%), with a smaller optimality gap indicating that the output solution is closer to the true optimal solution.For both examples, the average optimality gaps of solving for optimal design is reduced to 10% or lower after 300 seconds.This shows that, for the sizes of the above examples, the proposed approach is appropriate to solve and returns a solution that is close to the exact optimal design.We also point out that for cases with K over 10 and/or n over 400, the proposed solution approach did not return a feasible solution within 300 seconds.For those cases, a heuristic and/or randomized solution approach can be a reasonable choice.q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Percentile(%) Relative_Improvement(%)

Conclusion
To improve the accuracy of personalized decision, this paper investigates the optimal design to minimize the maximum variance of estimated personalized treatment effects over different treatments and different covariates values.The resulting design objective matches the Eoptimal design criterion.To provide the optimal design of multiple treatment allocation in the presence of observed covariates, a SDP solution approach is applied to solve the optimization problem.The proposed solution approach is able to provide optimal designs efficiently for n from 100-300 as shown in our numerical results.We point out potential future directions.First, it is desired to develop more efficient reformulations or optimization algorithms to handle the proposed design objective with larger sizes.Second, in practice, it is useful to investigate online optimal allocation assuming that the users are recruited in an online dynamic fashion under the proposed design objective.
i.e., k * (z) = argmax k=1,2,...,K z ⊤ β β β k (6) Similar to (3), to improve the accuracy of personalized decisions, it is desired to reduce the variance of estimated µ k (z) = z ⊤ β β β k for each k ∈ {1, . . ., K} and each z ∈ Z.Following the optimal design perspective in Zhang et al. (2022), we can reduce the the variances by searching an optimal allocation of the treatments x ik 's in the presence of fixed covariates of recruited users as illustrated in Figure 1.Next, we investigate the design criterion that characterizes the variance of estimated µ k (z) = z ⊤ β β β k for each k ∈ {1, . . ., K} and any z ∈ Z.

Figure 1 :
Figure 1: An illustration of the motivation of the proposed approach for two treatment allocation (i.e., + and − in the figure): We aim to obtain an optimal design of the treatment allocations to reduce the variances of estimated personalized treatment effects for users with any covariates z.By reducing these variances, it is clearer how to make the decision accurately.
size Kp, we can express the least squares estimator of b by b

Figure 2 :
Figure2: The quality measure Percentile (given by (14)) of different cases in Example I, i.e., Percentile of the optimal objective λ opt in the population of objectives given by 1000 random designs.