Copula-based robust optimal block designs

Blocking is often used to reduce known variability in designed experiments by collecting together homogeneous experimental units. A common modelling assumption for such experiments is that responses from units within a block are dependent. Accounting for such dependencies in both the design of the experiment and the modelling of the resulting data when the response is not normally distributed can be challenging, particularly in terms of the computation required to find an optimal design. The application of copulas and marginal modelling provides a computationally efficient approach for estimating population-average treatment effects. Motivated by an experiment from materials testing, we develop and demonstrate designs with blocks of size two using copula models. Such designs are also important in applications ranging from microarray experiments to experiments on human eyes or limbs with naturally occurring blocks of size two. We present methodology for design selection, make comparisons to existing approaches in the literature and assess the robustness of the designs to modelling assumptions.


INTRODUCTION AND MOTIVATION
Statistical design of experiments underpins much quantitative work in the biological, physical, and engineering sciences, providing a principled approach to the efficient allocation of (typically sparse) experimental resources to address the aims of the study. Often, experiments aim to understand a process by modeling discrete data, for example, arising from the observation of a binary or count response. For completely randomized experiments, assuming homogeneous experimental units, a generalized linear model (GLM) may provide an appropriate description, and there has been much research into the construction of optimal and efficient designs for multifactor GLMs. [1][2][3]4 See the work of Atkinson and Woods 5 for a comprehensive review.
When heterogeneous experimental units can be grouped into more homogeneous groups, or blocks, accounting for this grouping can improve the precision of inferences made from the experimental data. Methods to find block designs for discrete data have recently been proposed by Woods and vande Ven, 6 Niaparast and Schwabe, 7 and Waite and Woods, 8 among others. Two modeling paradigms have been adopted in the design literature: conditional models where the joint distribution of the data is derived by explicitly including block-specific random effects (eg, generalized linear mixed models 9 ) and marginal models, where the dependence structure of the data is specified separately from the marginal distribution of each response (eg, with parameters estimated via generalized estimating equations [GEEs] 10 ). For the linear model, these two modeling approaches coincide. In this paper, we find optimal designs under a marginal modeling approach when the intrablock dependence structure is defined via a copula. Such models are particularly appropriate when block effects are not of interest in themselves and the aim of the experiment is to understand the effects of treatment factors averaged across blocks. Optimal designs for marginal models using alternative definitions of the dependence structure have been found by other works. [11][12][13] Defining dependence via a copula model has the advantages of providing a flexible dependence modeling separate to the marginal probability models, and a more interpretable approach to defining the degree of dependence via commonly used measures; see Section 2 and, in particular, Section 2.3.
Although our methods can be generalized to arbitrary block sizes, we focus on the important special case of experiments with blocks of size two (see the work of Godolphin 14 ). Such blocks occur routinely in microarray experiments 15,16 and in experiments on people, for example, with eyes or arms as experimental units. 17 Practical motivation for our work comes from a materials science experiment. In Section 3, we find designs appropriate for aerospace materials testing experiments similar to those performed by our collaborators at the UK Defence Science and Technology Laboratory. The aim of these experiments is to compare the thermal properties of a set of novel materials against a reference material. In particular, one aim is to assess the probability of failure due to the exposure to extreme (high) temperatures. The experiment is performed using an arc jet to heat material samples, which are held in one of six "wedges," each of which holds a pair of samples on a strut attached to a circular carousel (see Figure 1). Hence, the experiment can be considered as a block design with six blocks, each containing two units. In the particular experiment considered here, six materials were tested, a reference and five novel samples. A variety of measures are made on each tested sample, including a visual inspection of quality to assess material failure, which leads to a binary (pass/fail) response. It is this response for which we find optimal designs.
In common with most nonlinear models, the performance of a given design for a copula-based GLM model may depend on the values of the model parameters that define both the marginal model and the dependence structure. If strong prior information is available, then locally optimal designs can be sought for given values of the model parameters. Otherwise, Bayesian 18 or maximin 19 approaches can be adopted. In common with much of the recent literature on designs for GLMs, we find optimal designs robust to the values of the model parameters via a pseudo-Bayesian approach 20(ch. 18) with a classical quantity for design performance averaged with respect to a prior distribution on the parameters. Here, we adopt variants of D-optimality for design selection.
The remainder of the paper is organized as follows. In Section 2, we introduce the statistical models we employ, including copulas, and develop design methods for blocked experiments. An illustrative comparison is made to previous design approaches based on GEEs using an example from the work of Woods and vande Ven. 6 In Section 3, we demonstrate and FIGURE 1 Arc jet carousel, struts, and "wedges" (left) and schematic (right). In addition to the six wedges for holding material samples, the carousel had two further wedges used for temperature measurement [Colour figure can be viewed at wileyonlinelibrary.com] assess our methods via application to the materials testing example. In particular, we show how prior information on the parameters influences the choice of optimal design. We provide a brief discussion and some areas for future work in Section 4.

DESIGNS FOR COPULA-BASED MARGINAL MODELS
Suppose the experiment varies m treatment factors, x T = (x 1 , … , x m ), and the experiment has b blocks of size k; throughout, our examples will assume k = 2. The jth unit in the ith block receives treatment x T i = (x 1i , … , x mi ) (i = 1, … , b; j = 1, … , k) and realizes observation Y i j . The x i j are chosen from a (typically standardized) compact design space , which could also be discrete and are not necessarily distinct. Independence of observations Y i , Y i ′ ′ , for i, i ′ = 1, … , b; j, j ′ = 1, … , k, is assumed across blocks (i ≠ i ′ ), but we allow dependence within a block (i = i ′ ), which we describe via a copula model.

Statistical modeling via copulas
The problem of specifying a probability model for dependent random variables Y i1 , … , Y jk can be simplified by express- and an associated k-copula (or dependence function) C defined as follows (cf the work of Nelsen 21 ).
and if all coordinates of u are 1 except u i , then The connection between a copula and a joint probability distribution is given by Sklar's theorem, 22 which affirms that for all 1 The copula C may not be unique for discrete margins; however, the practical limitations for statistical purposes are little; cf the work of Genest and Nešlehová. 23 Owing to Sklar's theorem, parametric families of copulas represent a powerful tool to describe the joint relationship between dependent random variables. Selecting the appropriate dependence within an assumed parametric copula family reduces to the selection of copula parameters, which correspond, for example, to a specific measure of association for the modeled random variables. Assuming Y i1 , … , Y i k are continuous random variables with associated copula C(·; ), one measure of association proposed by Joe 24 is given by Equation (2) is a generalized version of Kendall's and, hence, establishes a correspondence between a scalar copula parameter and the degree of dependence. More details and properties of this quantity, and another more traditional measure of concordance can be found in the work of Genest et al. 25

Design of experiments for copula models
In common with most work on optimal design of experiments, we base our criterion on the Fisher information matrix (FIM), the inverse of which provides an asymptotic approximation to the variance-covariance matrix of the maximum likelihood estimators of the model parameters.
We will work within a class of normalized block designs defined as with n ≤ b distinct (support) blocks. As defined, bw i must be an integer and it represents the replication of the ith support block (i = 1, … , n). Without loss of generality, we assume the first n blocks in the design corresponding to 1 , … , n , with the remaining b − n blocks being replicates. We relax the assumption that bw i is the integer to find the so-called approximate or continuous designs; see also the works of Cheng 26 and Waite and Woods. 8 Let Ξ denote the space of all possible designs of this form. Denote the vector of responses from the ith block as with corresponding expectation vector where (·; ·) is a known function and = ( 1 , … , r ) T is a vector of unknown parameters requiring estimation. Denote the marginal distribution function for the jth entry in the block as F Y i ( i ; x i , ), j = 1, … , k, and denote the joint distribution, derived via a copula transformation, for the k responses in the ith block as C( is the joint density function represented through a copula C in accordance with Equation (1). The FIM for an approximate block design is then given by An optimal design ⋆ maximizes a scalar function {M( ; )} of the information matrix. Previous work on optimal designs for copulas has focused on finding completely randomized locally optimal designs for multivariate responses, which can be considered as a block design where every unit within a block must receive the same treatment. We generalize these methods to allow different treatments for each unit within each block. Denman et al 27 found D-optimal designs for a bivariate response (k = 2) that maximized D {M( ; )} = det M( ; ), and Perrone and Müller 28 developed a corresponding equivalence theorem. These methods were extended to the local D A -criterion and, as a special case, for the D s -criterion in the work of Perrone et al. 29 Other relevant uses of design of experiments in copula models are those by Deldossi et al 30 and Durante and Perrone, 31 but until now, all relied on the availability of a single "best guess" vector of parameter values. It is well recognized that, for many nonlinear models, optimal designs for particular values of the model parameters may be very inefficient under different values; see the work of Woods et al 1 for the case of scalar response GLMs.
To overcome this dependence on assumed parameter values, here we adopt a pseudo-Bayesian approach for constructing block designs. Adopting this approach provides a more robust approach to design than the localized methods provided by Perrone et al. 29 Furthermore, our primary interest is typically in s meaningful linear combination of the parameters. Such combinations can be defined as elements of the vector A T , where A T is an s × (r + l) matrix of rank s < (r + l).
If M( ; ) is nonsingular, the variance-covariance matrix of the maximum likelihood estimator of A T is proportional to A T {M( ; )} −1 A. Hence, we define a robust D A -optimal block design ⋆ as the design that maximizes where G( ) is a proper prior distribution function for and Γ ⊂ R r+l is the support of G. See also the work of Woods and vande Ven. 6 Most often the main interest is in an s < (r + l)-dimensional subset of the parameters. In such a case, a robust D s -optimal block design can be found by maximizing following the partition of the information matrix as ) .
Here, M 11 is the (s × s) partition related to the parameters of interest. This criterion follows as a special case of the D A -criterion with A T = (I s 0 s×(r+l−s) ), where I s is the s × s identity matrix and 0 s×(r+l−s) is the s × (r + l − s) zero matrix. We evaluate a design via its Bayesian efficiencies under a given criterion, relative to an appropriate reference design * (see, for example, the work of Waite 32 ). Under robust D s -optimality, this efficiency is given by .
We find designs that maximize (4) and (5) numerically using a version of the Fedorov-Wynn algorithm, 33,34 as implemented in R package docopulae. 35 The optimality of a block design ⋆ under the robust D A -criterion, regardless of how it was found, can be assessed via application of the following Kiefer-Wolfowitz-type equivalence theorem. The proof is similar to that for completely randomized experiments with multivariate response; see the work of Perrone et al 29 for the locally optimal design case. Theorem 1. The following properties are equivalent: 3. over all ∈ Ξ, the design ⋆ minimizes the function where Ξ is the set of all possible block designs.

Comparative example
We demonstrate robust optimal block designs for copula models using a simple example from the work of Woods and vande Ven, 6 which allows comparison to the designs found by those authors for a GEE model. We again stress that the copula approach allows us to explicitly specify this dependence in terms of interpretable quantities (eg, generalized Kendall's ), unlike the GEE model under which the dependence is only specified implicitly. We find robust designs for a single-factor log-linear regression model assuming Poisson marginal distributions and quadratic linear predictor, implying log{ (x; )} = 0 + 1 x + 2 x 2 . The prior distribution G is uniform on the parameter space [−1, 1] × [4,5] × [0.5, 1.5]. In line with our motivating example, we assume blocks of size k = 2 and intrablock dependence defined according to one of the following bivariate copula functions.
The first copula is chosen for reference purposes; the latter two represent opposing dependencies in the tails (lower-tail dependence for the Clayton versus upper-tail dependence for the Gumbel). To isolate the effect of the copula structure from the strength of the dependence, we set for each copula such that the values for Kendall's coincide at three levels 2 = , 1∕10, 1∕3, respectively. Here, = 10 −9 > 0 is a small number to approximate the zero case but avoids singularity issues. Note that, as blocks always consists of pairs of points, here the design space for each block is [−1, 1] 2 .
To find robust D-optimal designs, objective function (4) was evaluated using quadrature. 36 Optimal designs under the Clayton and Gumbel copulas are shown in Figure 2 and demonstrate that increasing the generalized dependence (ie, increasing 2 ) leads to designs placing more weight on support blocks with points on the edge of the design space. All the designs display a "mirror-image" structure, with all design points having x > 0. These features are common in designs for Poisson regression (see the work of Russell et al 4 ). The designs found under the Gumbel copula tend to include more support blocks, but the pattern in the changes to these blocks as 2 is increased is similar for both copulas.  For reference purposes, the optimal design using the independence copula, ie, an optimal design assuming no block effect, was evaluated. It showed little difference to setting the nominal level for 2 = 0 for a particular copula. In particular, the D-efficiencies (with respect to the reference design assuming no block effect) for the Clayton and Gumbel model were 96.3% and 99.7%, respectively. This efficiency expectedly decreases as the association within the block increases, for 2 = 1∕3; for instance, it is already down to 65.0% and 61.3%, respectively. In the work of Woods and vande Ven, 6 robust D-optimal designs were found under the same Poisson marginal models and prior distribution but with the dependence described using a GEE approach with an exchangeable correlation matrix and pairwise working correlation of 0.5. The optimal design found was given by That is, for example, the first support block is 1 = (0.03, 1). This design is somewhat different in structure to the copula designs, without the same mirror structure. Quantitatively, the comparison shows the efficiencies under various scenarios given in Table 1. Surprisingly, the design from the work of Woods and vande Ven 6 seems to be most compatible with an independence assumption.

APPLICATION TO THE MATERIALS EXAMPLE
In this section, we return to the materials testing example to find and assess designs for comparing six materials in block of size two under a variety of modeling assumptions. The measured response is binary, with each material sample either passing or failing a visual check. We label the five novel materials as "treatments," with the reference material considered as a control. Marginally, we assume a logistic regression to model the differences between materials set up as where expit(u) = 1∕{1 + exp(−u)}, Y i j is the binary response from the ith unit in the jth block (i = 1, 2; j = 1, … , b), (x i j ; ) is the associated probability of success, x ijl is an indicator variable taking the value 1 if the ith unit in the jth block was assigned treatment l (l = 1, … , 5) and 0 otherwise, and 0 , … , 5 are unknown parameters to be estimated. Here, 0 is the logit for the reference material, with l being the difference in expected response, on the logit scale, between the reference material and the lth novel material or treatment. The choice of copula and the strength of intrablock association makes little difference to the design selected. However, assuming different marginal models and adopting a local or pseudo-Bayesian approach has a strong impact on the designs. The impact of the marginal model here is not surprising, as the degree of dependence between binary random variables is also strongly determined by their marginal distributions. 37(ch. 7) Example designs for the Gumbel copula are shown in Figure 3. Note that the numbers 1, … , 6 must be understood as nominal labels.
With a null marginal model, ie, T = (0, 0, 0, 0, 0, 0), when the response variance is constant, the locally D-optimal design contains all material combinations, excluding those blocks containing replicates of a single treatment. This design would also be optimal under a linear model with constant error variance. For different assumed parameter vectors, for example T = (0, −1, 2, −3, 4, −5), the optimal design contains only a few distinct treatment and treatment control combinations, with differing weights; here, (1,2), (3,4), (4,5), and (5,6)  , adjusts the weighting of the support blocks to give more emphasis on comparing treatments 2 and 4 and treatments 3 and 5. These pairs of treatments have coefficients with differences to the control with the same sign. This last prior distribution is representative of the type of materials study that might be conducted in practice, with differing prior beliefs about the size and direction of the difference between the expected response from each treatment and the reference material.

DISCUSSION
The modeling of block effects by copulas seems to be a natural choice and allows for elegant separation of the block and the marginal effects. Experimental designs for such models are now readily calculable. The pseudo-Bayesian D A -optimality criterion was added to the R package docopulae version 0.4 (see the work of Rappold 35 ) with the functions wDsensitivity and wDefficiency, both relying on a prespecified quadrature scheme for evaluation of the integrals. In this paper, we have concentrated on finding designs to estimate the complete parameter vector but the implementation provides flexibility for checking for symmetry, model discrimination, etc, as investigated in the work of Perrone et al. 29 Our examples are confined to the case k = 2. While there is no theoretical necessity for that, it is difficult to specify high-dimensional parametric copulas with a sufficient range of dependence, for details, see the excellent survey of Nikoloulopoulos. 38 However, work on this issue would go well beyond the scope of this paper. It might also be interesting to contrast our findings with some known analytic results for blocks of size two as, for example, given in the work of Cheng, 26 where a Gaussian copula is implicitly assumed.
Science Fund (FWF). D.C. Woods was partially supported through Fellowship EP/J018317/1 from the UK Engineering and Physical Sciences Research Council. We are also grateful to the associate editor and two referees, whose comments have lead to a considerable improvement of this paper.