Assume the non-informative priors p ( β ( i ) ) ∝ 1 for i = 1 , 2 and p ( η ) ∝ 1. The conditional posterior distribution in (2.3) is the same as the joint distribution of the data. It can be further factorized into p ( β ( 1 ) , β ( 2 ) , η | y , z , μ 0 , σ 2 , X ) ∝ p ( η | z , X ) ∏ i = 1 2 p ( β ( i ) | y , z , μ 0 , σ 2 , X ) , with the posterior distributions (3.1) β ( i ) | y , z , μ 0 , σ 2 , X ∼ N ( F ′ V i F ) − 1 F ′ V i ( y − μ 0 1 ) , σ 2 ( F ′ V i F ) − 1 (3.2) for i = 1 , 2 , and p ( η | z , X ) ∝ exp ∑ i = 1 n z i f ( x i ) ′ η − log ( 1 + e f ( x i ) ′ η ) . Conditioning on z, the posterior distributions of ( β 1 , β 2 ) and η are independent. Note that the noninformative prior p ( η ) ∝ 1 is proper because it leads to proper posterior p ( η | z , X ). Under the noninformative priors, the Bayesian estimation is identical to the frequentist estimation.

Using the posterior distributions (3.1)–(3.2) and the criterion function ψ ( · ) = log ( · ) in the general Bayesian optimal design criterion (2.4), we obtain the Bayesian D-optimal design criterion (3.3). (3.3) Ψ ( X | μ 0 , σ 2 ) = E z , η log ( p ( η | z , X ) ) + 1 2 ∑ i = 1 2 E η E z | η log det { ( F ′ V i F ) } + constant . The derivation is in Supplement S1. The first additive term in (3.3) is exactly the Bayesian D-optimal design criterion for GLMs. Unfortunately, its exact integration is not tractable. The common approach in experimental design for GLMs is to use a normal approximation for the posterior distribution p ( η | z , X ) [3, 28]. Such an approximation leads to (3.4) E z , η log ( p ( η | z , X ) ) ≈ E η { log det I ( η | X ) } + constant , where I ( η | X ) is the Fisher information matrix. We can easily show that I ( η | X ) = − E z ∂ 2 l ( z , η | X ) ∂ η ∂ η T = ∑ i = 1 n f ( x i ) f ( x i ) ′ π ( x i , η ) ( 1 − π ( x i , η ) ) = F ′ W 0 F , where W 0 is a diagonal weight matrix W 0 = diag { π ( x 1 , η ) ( 1 − π ( x 1 , η ) ) , … , π ( x n , η ) ( 1 − π ( x n , η ) ) } . Omitting the irrelevant constant, we approximate the exact criterion Ψ ( X | μ 0 , σ 2 ) in (3.3) as follows. (3.5) Ψ ( X | μ 0 , σ 2 ) ≈ E η { log det ( F ′ W 0 F ) } + 1 2 ∑ i = 1 2 E η E z | η log det ( F ′ V i F ) .

To construct the optimal design, we consider maximizing the approximated Ψ ( X | μ 0 , σ 2 ) in (3.5). But this is not trivial, because it involves the expectation on Z i ’s in the second additive term. To overcome this challenge, [18] constructed optimal designs by simulating samples from the joint distribution of responses and the unknown parameters. But this method can be computationally expensive for even slightly larger dimensions of experimental factors. Instead of simulating Z i ’s, we derive the following Theorem 1 that gives a tractable upper bound Q ( X ). Thus we propose using the upper bound Q ( X ) as the optimal criterion.

The proof of Theorem 1 is in Supplement S1. Note that Theorem 1 requires that F ′ W i F for i = 0 , 1 , 2 are all nonsingular. Obviously W 0 = W 1 W 2 . It is easy to see that F ′ W 0 F is nonsingular if and only if both F ′ W 1 F and F ′ W 2 F are nonsingular.

The matrices F ′ V 1 F and F ′ V 2 F involve the responses Z i ’s that are not yet observed at the experimental design stage. We can only choose the experimental run size and the design points to avoid the singularity problem with a larger probability for given values of η. Once the run size is chosen, the design points can be optimally arranged by maximizing Q ( X ). The weight matrix W 1 (or W 2 ) gives more weight to the feasible design points that are more likely to lead to Z = 1 (or Z = 0) observations so that the parameters β ( 1 ) (or β ( 2 ) ) of the linear model Y | Z = 1 (or Y | Z = 0) are more likely to be estimable. Next, we introduce some regularity conditions on the run size and number of replications to alleviate the singularity problem.

Let m be the number of distinct design points in the design matrix X, n i be the number of repeated point x i in X for i = 1 , … , m. Thus n = ∑ i = 1 m n i and n i ≥ 1 for i = 1 , … , m. First, it is necessary that m ≥ q for the linear regression model to be estimable under the noninformative priors. The if-and-only-if condition for F ′ W 0 F to be nonsingular is that rank ( F ′ W 0 F ) ≥ q. If m ≥ q and π ( x i , η ) ∈ ( 0 , 1 ) for i = 1 , … , m, then F ′ W i F for i = 0 , 1 , 2 are all nonsingular and thus positive definite. To make sure π ( x i , η ) ∈ ( 0 , 1 ) for i = 1 , … , m, it is sufficient to assume that η is finitely bounded. This condition is typically used for the frequentist D-optimal design for GLMs. For instance, [44] suggested using the centroids of the finite bounded space of η to develop the local D-optimal design for GLMs. [13] clustered different local D-optimal designs with η randomly sampled from its bounded space.

The if-and-only-if condition for F ′ V 1 F and F ′ V 2 F being nonsingular is (3.7) ∑ i = 1 m I ( ∑ j = 1 n i Z i j > 0 ) ≥ q and ∑ i = 1 m I ( ∑ j = 1 n i Z i j < n i ) ≥ q , where Z i j is the jth random binary response at the unique design point x i and I ( · ) is the indicator function. In the following, we discuss how to choose sample sizes n i and n under two scenarios (i) m = q and (ii) m > q.

Proposition 1 gives a sufficient condition on the lower bound of n i when saturated design ( m = q) is used. Under the sufficient condition, the nonsingularity of F V 1 F and F V 2 F holds with a probability larger than κ m . For Example 1 in Section 1, suppose that if the possible values of x can only be − 1 , 0 , 1, then m = q = 3. Let η = ( 1 , 1 ) ′ . If κ = 0.5, then the numbers of replications for x = − 1 , 0 , 1 need to satisfy n 1 ≥ 2, n 2 ≥ 4, and n 3 ≥ 7, respectively. If κ = 0.9, then n 1 ≥ 5, n 2 ≥ 9, and n 3 ≥ 20. Proposition 1 is useful in Step 1 to construct the initial design in Algorithm 2 in Section 5.

Proposition 2 provides one sufficient condition and one necessary condition when m > q on the smallest number of replications and the overall run size. But these conditions only examine the nonsingularity of the two matrices with large probability, which is weaker than Proposition 1. For given η value, Algorithm 2 in Section 5.1 can return the local D-optimal design. Proposition 2 can be useful to check the local D-optimal design, as π min and π max depend on η. Take the artificial example in Section 6.1 for instance. The local D-optimal design for ρ = 0 (Table S1 in Supplement S2) has m = 50 unique design points and there are q = 22 effects. According to Proposition 2, the sufficient condition requires n 0 ≥ 7 and the necessary condition requires n 0 ≥ 1. The local D-optimal design in Table S1 only satisfies the necessary condition. To meet the sufficient condition n has to be much larger. For the global optimal design considering all possible η values, Proposition 2 can provide some guidelines for the design construction when η is varied in a relatively small range.

For the parameters β ( 1 ) and β ( 2 ) , the conjugate priors are normal distribution since Y | Z follows normal distribution. Thus we consider their priors as β ( 1 ) ∼ N ( 0 , τ 2 R 1 ) , β ( 2 ) ∼ N ( 0 , τ 2 R 2 ) . where τ 2 is the prior variance and R i is the prior correlation matrix of β ( i ) for i = 1 , 2. Here we use the same prior variance τ 2 only for simplicity. The matrix R i can be specified flexibly such as using ( F ′ F ) − 1 , or those in [20] for factorial designs.

For the parameters η, we choose the conjugate prior derived in [4]. It takes the form (4.1) η ∼ D ( s , b ) ∝ exp ∑ i = 1 n s b i f ( x i ) ′ η − log ( 1 + e f ( x i ) ′ η ) , where D ( s , b ) is the distribution with parameters ( s , b ). Here s is a scalar factor and b ∈ ( 0 , 1 ) n is the marginal mean of z as shown in [10]. The value of b can be interpreted as a prior prediction (or guess) for E ( Z ). Based on the priors for ( β ( 1 ) , β ( 2 ) , η ) we can derive the posteriors as follows.

The proof of Proposition 3 can be derived following the standard Bayesian framework, thus is omitted. To derive the Bayesian D-optimal design criterion, we take the posterior distributions in Proposition 3 to (2.4) and set ψ ( · ) = log ( · ). The derivation is very similar to that in (3.3), and thus we obtain the exact design criterion as (4.2) Ψ ( X | μ 0 , σ 2 ) = E z , η log ( p ( η | z , X ) ) + 1 2 ∑ i = 1 2 E η E z | η log det ( F ′ V i F + ρ R i − 1 ) + constant . As the integration of E z , η log ( p ( η | z , X ) ) is not tractable, we adopt the same normal approximation of the posterior distribution p ( η | z , X ) as in (3.4). A straightforward calculation leads to getting Fisher information matrix I ( η | X ) = ( 1 + s ) F ′ W 0 F. Thus we have (4.3) E z , η log ( p ( η | z , X ) ) ≈ E η { log det ( F ′ W 0 F ) } + constant . Disregarding the constant, we can approximate the exact Ψ ( X | μ 0 , σ 2 ) by Ψ ( X | μ 0 , σ 2 ) ≈ E η { log det ( F ′ W 0 F ) } + 1 2 ∑ i = 1 2 E η E z | η log det ( F ′ V i F + ρ R i − 1 ) . The following Theorem 2 gives an upper bound of the approximated criterion Ψ ( X | μ 0 , σ 2 ) to avoid the integration with respect to z, which plays the same role as Theorem 1.

For the same argument as in Section 3, we use the upper bound in (4.4) as the optimal design criterion. Note that since ρ R i − 1 is added to F ′ V i F and F ′ W i F, F ′ V i F + ρ R i − 1 and F ′ W i F + ρ R i − 1 are nonsingular. The derivation of Ψ ( X | μ 0 , σ 2 ) and Q ( X ) only needs F ′ W 0 F to be nonsingular, which requires m ≥ q and η to be finitely bounded as in Section 3.

Note that the criterion in (4.4) has a similar formulation with Q ( X ) in (3.6). The only difference is that (3.6) does not involve ρ R i − 1 . For consistency, we use the formula (4.4) as the design criterion Q ( X ) for both cases. When noninformative priors for β ( 1 ) and β ( 2 ) are used, we set ρ = 0. From another point of view, as τ 2 → ∞, ρ → 0, the variances in the priors p ( β 1 ) and p ( β 2 ) diffuse and result in a noninformative priors.

The criterion Q ( X ), consisting of three additive terms, can be interpreted intuitively. The first additive term E η { log det ( F ′ W 0 F ) } is known as the Bayesian D-optimal criterion for logistic regression and E η { log det ( F ′ W i F + ρ R i − 1 ) } is the Bayesian D-optimal criterion for the linear regression model of Y. To explain the weights, we rewrite Q ( X ) as follows. Q ( X ) = 1 · E η log det ( F ′ W 0 F ) + 1 · 1 2 ∑ i = 1 2 E η ( log det ( F ′ W i F + ρ R i − 1 ) ) . Since there are equal numbers of binary and continuous response observations, the design criterion should put the same weight (equal to 1) on both design criteria for Z and Y. For the two criteria for the linear regression models, the same weight 1/2 is used. This is also reasonable because we assume π i ∈ ( 0 , 1 ). Then none of the diagonal entries of W 1 and W 2 are zero, so the two terms should split the total weight 1 assigned for the entire linear regression part. Therefore, even though Q ( X ) are derived analytically, all the additive terms and their weights make sense intuitively.

Note that the conjugate prior p ( η ) requires prior parameters ( s , b ) to be specified. Moreover, the prior distribution (4.1) contains f ( x i ), which depends on the design points. When sampling η from the prior (4.1), it does not matter whether f ( x i )’s are actually from the design points. If relevant historical data is available, we can simply sample η from the likelihood of the data. Alternatively, one can adopt the method in [4] to estimate the parameters ( s , b ). Without the relevant data, we would use the noninformative prior for η, i.e., p ( η ) ∝ 1 in the bounded region for η.

The design criterion Q ( X ) contains some unknown parameters, including the noise-to-signal ratio ρ = σ 2 / τ 2 and the correlation matrices R i ’s for i = 1 , 2. The value of ρ has to be specified either from the historical data or from the domain knowledge. Typically we would assume ρ < 1 such that the measurement error has a smaller variance than the signal variance.

The setting of R i can also be specified flexibly. If historical data are available, R i can be set as the estimated correlation matrix of β ( i ) . Otherwise, we can use the correlation matrix in [20] and [24], which is targeted for factorial design. Specifically, let β be the unknown coefficients of the linear regression model and the prior distribution is β ∼ N ( 0 , τ 2 R ). For 2-level factor coded in − 1 and 1, [20] suggests that R is a diagonal matrix and the priors for individual β j is (4.5) β 0 ∼ N ( 0 , τ 2 ) , β j ∼ N ( 0 , τ 2 r ) , i = 1 , … , p , β j ∼ N ( 0 , τ 2 r 2 ) , i = p + 1 , … , p + p 2 , ⋮ β 2 p − 1 ∼ N ( 0 , τ 2 r p ) , where β j i = 1 , … , p are main effects, β j j = p + 1 , … , p + p 2 are 2-factor-interactions and up to the p-factor-interaction β 2 p − 1 . The variance of β j decreases exponentially with the order of their corresponding effects by r ∈ ( 0 , 1 ), thus it incorporates the effects hierarchy principle [46]. [20] showed that if f ( x ) contains all the 2 p effects of all p orders, τ 2 R can be represented alternatively by Kronecker product as τ 2 R = ς 2 ⨂ j = 1 p F j ( x j ) − 1 Ψ j ( x j ) ( F j ( x j ) ) − 1 . The model matrix for the 2-level factor and the correlation matrix are (4.6) F j ( x j ) = 1 − 1 1 1 and Ψ j ( x j ) = 1 ζ ζ 1 . To keep the two different presentations equivalent, let ζ = 1 − r 1 + r and τ 2 = ( 1 + ζ 2 ) p ς 2 . For the mixed-level of 2- and 3-level experiments, [24] have extended the 2-level case to the format (4.7) τ 2 R = ς 2 ⨂ j = 1 p 2 + p 3 , c + p 3 , q F j ( x j ) − 1 Ψ j ( x j ) ( F j ( x j ) − 1 ) ′ , where p 2 is the number of 2-level factors, p 3 , c is the number of 3-level qualitative (categorical) factors, and p 3 , q is the number of 3-level quantitative factors. For all the 3-level factors, the model matrix is (4.8) F j ( x j ) = 1 − 3 2 1 2 1 0 − 2 1 3 2 1 2 . An isotropic correlation function is recommended for the 3-level qualitative factors and a Gaussian correlation function for quantitative factors. Thus, the correlation matrices for the 3-level qualitative and quantitative factors are (4.9) Ψ j ( x j ) = 1 ζ ζ ζ 1 ζ ζ ζ 1 and Ψ j ( x j ) = 1 ζ ζ 4 ζ 1 ζ ζ 4 ζ 1 , respectively. To keep the covariance (4.7) consistent with the 2-level case we still set ζ = 1 − r 1 + r . To keep the variance of the intercept equal to τ 2 [24], we set τ 2 = ς 2 1 + ζ 2 p 2 1 + 2 ζ 3 p 3 , c 3 + 4 ζ + 2 ζ 4 9 p 3 , q , and thus R = 1 + ζ 2 p 2 1 + 2 ζ 3 p 3 , c 3 + 4 ζ + 2 ζ 4 9 p 3 , q − 1 × ⨂ j = 1 p 2 + p 3 , c + p 3 , q F j ( x j ) − 1 Ψ j ( x j ) ( F j ( x j ) − 1 ) ′ . It is straightforward to prove that R is a diagonal matrix if only 2-level and 3-level qualitative factors are involved, but not so if any 3-level quantitative factors are involved, and the first diagonal entry of R is always 1.

To specify different prior distributions for β ( 1 ) and β ( 2 ) , we only need to use different values r 1 (or ζ 1 ) and r 2 (or ζ 2 ) to construct the prior correlation matrix. If the prior knowledge assumes that the two responses Z and Y are independent, one can set r 1 = r 2 = r so that the two correlation matrices R i ’s are the same, denoted as R. [24] has used r = 1 / 3 (equivalently ζ = 1 / 2) according to a meta-analysis of 113 data sets from published experiments [30]. Thus we also use r = 1 / 3 in all the examples. The readers can specify different values for r 1 and r 2 if needed.