Bayesian decision-theoretic design of experiments under an alternative model
Antony M. Overstall James M. McGree




#1 1




1
  
  
    
  Bayesian decision-theoretic design of experiments under an alternative model
  

  
    Antony M. Overstall

    Southampton Statistical Sciences Research Institute,
 University of Southampton,
 Southampton, UK
 (A.M.Overstall@soton.ac.uk) 


    James M. McGree

    School of Mathematical Sciences,
 Queensland University of Technology,
 Brisbane, Australia 
 (james.mcgree@qut.edu.au)
    

===============================================================================================================================================================================================================================================================================================================

 

0
  

    Title


   


Traditionally Bayesian decision-theoretic design of experiments proceeds by choosing a design to minimise expectation of a given loss function over the space of all designs. The loss function encapsulates the aim of the experiment, and the expectation is taken with respect to the joint distribution of all unknown quantities implied by the statistical model that will be fitted to observed responses. In this paper, an extended framework is proposed whereby the expectation of the loss is taken with respect to a joint distribution implied by an alternative statistical model. Motivation for this includes promoting robustness, ensuring computational feasibility and for allowing realistic prior specification when deriving a design. To aid in exploring the new framework, an asymptotic approximation to the expected loss under an alternative model is derived, and the properties of different loss functions are established.  The framework is then demonstrated on a linear regression versus full-treatment model scenario, on estimating parameters of a non-linear model under model discrepancy and a cubic spline model under an unknown number of basis functions.


Keywords:  cubic spline basis; expected loss function; full-treatment model; model discrepancy; non-linear model; normal linear model.



§ INTRODUCTION
 

The Bayesian decision-theoretic approach <cit.> is a natural framework to plan experiments in many fields of science and engineering. It starts with specification of a loss function representing the aim of the experiment. A Bayesian design then minimises the expectation of the loss over the space of all possible designs where expectation is with respect to a probability distribution over all unknown quantities implied by the statistical model that will be fitted upon observation of the experimental responses. 

The approach is scientifically appealing since the precise aim of the experiment is encoded by the loss function and, through specification of a probability distribution over all unknown quantities, pre-experimental (lack of) knowledge is incorporated. In practice, finding a Bayesian design by minimising the expected loss is a challenging computational problem <cit.>. The expected loss is typically analytically intractable and the space of all designs potentially high-dimensional. However, recently there has been significant progress in the development of novel computational methodology for finding Bayesian designs originating in diverse fields of study <cit.>.

In this paper, we propose an extended framework for Bayesian decision-theoretic design of experiments. This is achieved by defining the expected loss by taking expectation with respect to a probability distribution implied by an alternative statistical model (termed a designer model).  <cit.> considered the case where the fitted and designer prior distributions were different but the joint distributions for the responses were identical. We extend this by considering the case where both the prior and distribution for responses can be different. There are several reasons why a design may be sought under a different model to the fitted model. For example, the fitted model may emulate a partially observed process, but by considering the whole data-generating process, robustness is introduced into the design procedure. Assuming a simpler fitted model can also be used to induce robustness to model misspecification <cit.> or may be necessary to ensure computational feasibility of finding Bayesian designs <cit.>. Conversely, reliable prior specification can be assisted by the designer model being simpler than the fitted model.  Accordingly, we provide the framework and supporting methodology to enable Bayesian design within such settings, and apply these methods to typical design problems found in the literature.

The outline of the paper is as follows. In Section <ref>, the extended framework for Bayesian design of experiments under an alternative model is described and justified, and the choice of loss function considered. In Section <ref>, an asymptotic approximation to the expected loss under the extended framework is developed to explore the consequences of designing under an alternative model. Under the traditional approach, designs found under a loss, and its posterior expectation, are equivalent. This is no longer the case under the proposed extended framework. We consider the choice between these pairs of loss functions and show that there is essentially a bias-variance trade off to be made. To demonstrate the extended framework, three examples are considered in Section <ref> and Section SM6 in the Supplementary Material. The first considers a linear regression versus full-treatment model scenario where application of the extended framework results in objective functions sharing properties of classical design criteria that promote increased replication <cit.>. The second considers estimating parameters of non-linear models under model discrepancy, and finds that the resulting designs resemble a compromise of designs found under the traditional Bayesian decision-theoretic framework and space-filling designs. Finally, the third (in the Supplementary Material), considers estimation of the parameters of a model formed using cubic splines, with an unknown number of basis functions.

The Supplementary Material includes proofs and derivations of the results, and the cubic spline model example.
 


§ BAYESIAN DESIGN UNDER AN ALTERNATIVE MODEL
 
 




 §.§ Internal and external expected loss
 

Suppose the ultimate aim of the experiment is to learn the relationship between k controllable variables and a response. The experiment consists of n runs where, for i=1,…,n, the ith run involves specifying a treatment-level combination, i.e. values for the k controllable variables Δ_i = (Δ_i1, …, Δ_ik)^, and subsequent observation of a response y_i. Let Δ denote the design; the n × k matrix with ith row given by Δ_i^ and let = (y_1,…,y_n)^ denote the n × 1 vector of observed responses. 
 
On the completion of the experiment, a statistical model is assumed describing the relationship between the controllable variables and response. That is,  is assumed to be a realisation from a multivariate probability distribution with density/mass function π(|, F,Δ) completely known up to a p × 1 vector of unknown parameters  with prior density π(|F,Δ). We refer to π(|, F,Δ) and π(|F,Δ) as the fitted likelihood and prior, respectively, and collectively as the fitted model, where the conditioning on F is to make it clear that this is under the fitted model. Using Bayes' theorem, the fitted posterior distribution is π(|, F,Δ) ∝π(|, F,Δ) π(| F,Δ). Note that  refers to parameters of interest. There may be additional nuisance parameters present but these have been marginalised, with respect to a prior distribution, to obtain π(|, F,Δ).

Since, apart from the values of , the fitted model completely specifies the relationship between the controllable variables and response, the experimental aim reduces to determining the value of . This aim can manifest itself in various different ways. Bayesian decision-theoretic design of experiments starts by encoding the exact aim in a loss function, λ(, , F,Δ), giving the loss of estimating  via the fitted posterior distribution conditional on . The choice of loss function is considered in more detail in Section <ref>. 

Prior to experimentation, the aim is to specify Δ to best estimate  as measured by the loss function. Traditionally, a Bayesian design minimises the expected loss over the space of all designs, where expectation is with respect to the joint distribution of  and  as implied by the fitted model <cit.>. Mathematically, the expected loss is 

    L_FF(Δ) = E_, | F,Δ[ λ(, , F,Δ) ]

and referred to, in this paper, as the internal expected loss under the fitted model. Unless otherwise specified; if referring to internal expected loss, it is under the fitted model.

Instead, consider minimising the expectation of the loss with respect to an alternative model: the designer model. Here we suppose that  is a realisation of a probability distribution with density/mass π(|, D,Δ), where  is a p_θ× 1 vector of parameters with prior density π(|D,Δ). We refer to π(|, D,Δ) and π(|D,Δ) as the designer likelihood and prior, respectively, and the conditioning on D is to make it clear that this is under the designer model. Similar to the fitted model, any additional nuisance parameters have been marginalised to obtain π(|, D,Δ). 

Let = (^, ^)^, where  is a p_γ× 1 vector of parameters only present in the fitted model and  are parameters common to both models. If the fitted and designer models share no common parameters, i.e. =, then the two models are referred to as disjoint. Conversely, if all parameters of interest are common to both models, i.e. =, then the models are referred to as compatible. As such, prior to marginalisation, the fitted model can be nested within the designer model, vice-versa, or no nesting may be present.

A design found under the proposed extended framework is given by minimising the external expected loss defined as

    L_DF(Δ) = E_, | D,Δ{E_| , , F,Δ[ λ(, , F,Δ)] }.

The external expected loss is formed by assuming that once responses are observed, we proceed as if the fitted model represents the true data-generating process, i.e. under the M-closed paradigm <cit.>. Initially, expectation is taken with respect to the fitted posterior distribution of  conditional on . This is because these parameters are only present in the fitted model (if the models are non-compatible). Then expectation is taken with respect to the remaining unknowns under the designer model.

Note that if the fitted and the designer models are the same, then the external and internal expected loss functions are equal, i.e. L_DF(Δ) = L_FF(Δ).

There are several reasons to consider finding a design under an alternative model. These can be separated into two cases: when the designer model is simpler than the fitted model and the vice-versa. In Section <ref> and Section SM6 in the Supplementary Material we provide examples of both cases.

First consider the case where the fitted model actually represents the, a-priori, current best representation of the data-generating process, and the designer model is actually a proxy simplification. This can occur when the fitted model is complex with a large number of parameters and it is necessary to use a fitted prior distribution representing weak prior information. It is known <cit.> that in these cases, the vaguery of the fitted prior distribution combined with the large number of parameters leads to the internal expected loss being relatively flat and, therefore, non-trivial to (numerically) minimise. However it may be possible to elicit a prior distribution for a simpler designer model and use this to take expectation of the loss function formed under the fitted model.

Now consider the case where the designer model actually represents the best representation of the data-generating process and the fitted model is a simplification. In this case, it initially would make sense to find a design by minimising the internal expected loss under the designer model, i.e.

    L_DD(Δ) = E_, | D,Δ[ λ(, , D,Δ) ].

For certain models it may not be possible to represent the experimental aim by a suitable loss function, λ(, , D,Δ), depending on the designer posterior distribution. Instead, an alternative fitted model is used for which it is possible to form a suitable loss function. However it is prudent to take expectation under the original designer model since this is the best representation of the data-generating process. For example, suppose the designer model is a full-treatment model <cit.> with each unique treatment-level combination permitting a unique mean response. This represents a plausible data-generating process but for convenience and interpretability a regression model (positing a relationship between the numerical values of the controllable variables and mean response) is usually assumed as the fitted model.

Alternatively, suppose the designer model is based on scientific reasoning through a mechanistic or phenomenological model. Again it would make sense to find the design via the internal expected loss under the designer model (<ref>). Assuming the fitted and designer models are compatible, consider Monte Carlo approximations to the external and internal expected loss functions given by
L̂_DF(Δ) = B^-1∑_b=1^B λ(_b, _b, F,Δ) and L̂_DD(Δ) = B^-1∑_b=1^B λ(_b, _b, D,Δ), respectively, where {_b, _b }_b=1^B is a sample generated from the joint distribution of  and  under the designer model. There are models, e.g. so-called intractable likelihood models, for which it is straightforward to generate such a sample but extremely difficult to evaluate the designer posterior distribution and hence the loss function λ(_b, _b, D ,Δ), for b=1,…,B. However by forming an approximating fitted model it becomes computationally feasible to find a Bayesian design via the external expected loss <cit.>. In Section <ref>, we show, under a certain class of loss function, that the external expected loss is an upper bound on the internal expected loss under the designer model, providing formal justification for using such approaches.





































 §.§ Loss functions
 

In Section <ref>, the loss function is given informally by λ(, , F,Δ), stating that it gives the loss of estimating  via the fitted posterior (conditional on ). We now consider the specification of the loss and the implications of such a choice when designing an experiment under an alternative model. For this, we consider two exemplar loss functions as defined in Table <ref> (left-hand column); self-information and squared error. Note that ‖‖_2,A^2 = ^   A and I_p is the p × p identity matrix. These loss functions have been chosen since under a normal linear model with non-informative prior distributions, the internal expected loss reduces to the objective function under D- and A-optimality, respectively; two widely-studied classical optimal design criteria. 





  §.§.§ Optimal loss functions


First consider estimating  based on observed responses . Suppose the loss for this aim is λ_0(, ) where ∈ℬ is an action for estimating  and ℬ is the action space (e.g.  is a point estimate of  in ℬ = ℝ^p). The optimal action minimises the fitted posterior expectation of λ_0(, ), i.e. _F = argmin_∈ℬE_|, F,Δ[ λ_0(, )].

The loss function, as written in Section <ref>, is given by λ(, , F,Δ) = λ_0(, _F) and is optimal under the fitted posterior distribution. In this paper, only optimal loss functions are considered. Indeed, the self-information and squared error loss functions are optimal under logarithmic scoring and λ_0(, ) = ‖ - ‖_2,I_p^2, respectively. However, the formulation in Section <ref> allows consideration of sub-optimal loss functions, if deemed appropriate, e.g. 
λ(, , F,Δ) = ‖ - M_|, F,Δ( ) ‖_2,I_p^2, where M_|, F,Δ( ) is the fitted posterior mode of .








Under optimal loss functions, the following result giving an inequality for the external expected loss in relation to the internal expected loss under the designer model. The proof is given in Section SM1 of the Supplementary Material.

 
Suppose fitted and designer models are compatible, i.e. =, with λ(, , F,Δ) = λ_0(, _F) and λ(, , D,Δ) = λ_0(, _D), then L_DF(Δ) ≥ L_DD(Δ), where L_DD(Δ) is the internal expected loss under the designer model given by (<ref>).












In Lemma <ref>, the difference between the external expected loss and the internal expected loss under the designer model is due to the loss of information caused by the fitted model being different to the model under which expectation is taken. Note that Lemma 2.1 only requires the compatible fitted and designer models to be different, not that the designer model is more complex.

By finding the design that minimises the external expected loss, we are minimising an upper bound on the internal expected loss under the designer model. A special case of Lemma <ref> for the self-information loss is called the Barber-Agakov bound <cit.> and has been used by <cit.> to approximately find Bayesian designs. Lemma <ref> now justifies the use of this approach to other loss functions as has previously been pursued for intractable likelihood models (for example, , as discussed in Section <ref>).

Lemma <ref> is uninformative on the relationship between the external expected loss and the internal expected loss under the fitted model (<ref>).



  §.§.§ Generator and composite loss functions


There are intuitive choices of loss function given as scalar summaries of the fitted posterior distribution, independent of , e.g. λ(, F,Δ) = |var(|, F,Δ) |, where |·| denotes the determinant of a square matrix. In some cases, these are given as the fitted posterior expectation of a loss function λ_G(, , F,Δ), i.e. λ_C(, F,Δ) = E_|, F,Δ[ λ_G(, , F,Δ)]. For such pairs, we term λ_C(, F,Δ) as composite and λ_G(, , F,Δ) as generator. The right-hand column of Table <ref> shows the entropy and trace variance loss functions which are composite to the generator loss functions of self-information and squared error, respectively, in the left-hand column.

The internal expected composite loss is equal to the corresponding internal expected generator loss. Apart from when fitted and designer models are disjoint, i.e. =, this equivalence is not necessarily true when considering the external expected loss. This leads to a question on how to choose between a pair of loss functions equivalent under internal expected loss, e.g. how should we choose between self-information and entropy loss functions? In the next section, to answer this question, we develop an asymptotic approximation to the external expected loss. This is used to explore differences in behaviour between the exemplar generator and composite loss functions shown in Table <ref>. We find that a composite loss solely focuses on fitted posterior precision. On the other hand, a generator loss also incorporates point estimation of common parameters . This is essentially a bias-variance trade-off. Considering this trade-off is unnecessary under internal expected loss because the fitted model is, in a sense, “unbiased” as it is identical to the designer model.



§ UNDERSTANDING THE EXTERNAL EXPECTED LOSS
 

In this section, an asymptotic approximation to the external expected loss is derived.  This approximation is used to provide insight into properties of exemplar loss functions, in particular, to investigate the difference between external expected generator and composite loss functions.



 §.§ Asymptotic approximation of external expected loss
 

In general, the external expected loss (<ref>) is not available in closed form and use of an approximation is therefore necessary. Recent developments in computational methodology for approximating (and minimising) the expected loss (example references in Section <ref>) can actually be used for this task with little modification. Indeed, in the examples in Section <ref> and Section SM6 in the Supplementary Material we describe Monte Carlo approximations to the external expected loss function.

However, in this section, to gain understanding of the behaviour of the external expected loss, we develop an asymptotic approximation. The resulting approximations to the external expected loss under the exemplar loss functions in Table <ref> are derived and compared to the corresponding expressions under the internal expected loss. 

First, some preliminary definitions are required. Let = ( ^,^)^ be the p × 1 vector of parameter values for  that minimise the Kullback-Leibler divergence between the designer and fitted likelihoods for a given value of  under the designer model, i.e.  maximises

    d() = E_|,D,Δ[ logπ (|, F,Δ) ].

Define ℐ and ℐ̃ to be 

    E_|, D,Δ[ - ∂^2 logπ (|, F,Δ)/∂∂^]

evaluated at = and =, respectively. Similarly, define 𝒥̃ to be 

    E_|, D,Δ[ ∂logπ (|, F,Δ)/∂∂logπ (|, F,Δ)/∂^],

evaluated at =, respectively. Let 𝒦̃ = ℐ̃^-1𝒥̃ℐ̃^-1 with ℐ̃ and 𝒦̃ partitioned as

    ℐ̃ = [ [ ℐ̃_γγ ℐ̃_γθ; ℐ̃_θγ ℐ̃_θθ ]],         𝒦̃ = [ [ 𝒦̃_γγ 𝒦̃_γθ; 𝒦̃_θγ 𝒦̃_θθ ]]

where ℐ̃_γγ is the p_γ× p_γ sub-matrix corresponding to  with similar definitions for ℐ̃_θθ, ℐ̃_θγ= ℐ̃_γθ^ and the sub-matrices of 𝒦̃. 

Under regularity conditions (see Section SM2 in the Supplementary Material), an asymptotic approximation <cit.> to the fitted posterior distribution of  is given by

    N(, ℐ̃^-1),

where = ( ^,^)^ is the maximum likelihood estimate of  under the fitted model.

The loss function is approximated by replacing fitted posterior quantities by the corresponding approximated fitted posterior quantities under (<ref>). Such an approximation is denoted λ^*(, ,Δ). For example, the self-information loss is approximated by

    λ^*_SI(, ,Δ) = N_SI - 1/2log|ℐ̃| + 1/2‖ - ‖^2_2,ℐ̃,

where N_SI = p log (2π)/2 is a constant not depending on the design. 

To approximate the external expected loss (<ref>), the parameters present only in the fitted model, , are marginalised with respect to an approximate fitted posterior conditional on , induced from (<ref>). Then expectation with respect to the distribution of  conditional on  under the designer model can be approximated using asymptotic results for inference under an alternative (or misspecified) model <cit.>.

The result is now formally stated with proof given in Section SM2 of the Supplementary Material.

 
Under conditions given in Section SM2 of the Supplementary Material, an asymptotic approximation to the external expected loss is given by

    L^*_DF(Δ) = E_| D,Δ( E_|, D,Δ{E_|, , F,Δ[ λ^*(, ,Δ) ] })

where the two inner expectations are taken with respect to

    |, , F,Δ∼N( P̃_1  - P̃_2 , ℐ̃_γγ^-1),         |, D,Δ∼N( , 𝒦̃),

with P̃_1 = (I_p_γ, P̃_2 ) and P̃_2 = ℐ̃_γγ^-1ℐ̃_γθ.


In Theorem <ref>, 𝒦̃ has the sandwich form of the “robust variance formula" <cit.>. In inferential settings, the values of ℐ̃ and 𝒥̃ are estimated from observed responses, whereas, in design settings, the values follow from the assumption of a designer model.

Due to the tractability of the normal distribution, the two inner expectations in (<ref>) are available in closed form for many common loss functions, e.g. the exemplar loss functions in Table <ref>. Then the approximate external expected loss is written L^*_DF(Δ) = E_| D,Δ[ ℓ^*_G(,Δ) ]; the designer prior expectation of a function ℓ^*_G(,Δ), typically a function of ℐ̃ and 𝒦̃.



For loss functions independent of , including composite losses, the approximate loss is denoted λ^*(,Δ) and the approximate external loss given by
L^*_DF(Δ) = E_| D,Δ[ ℓ_C^*(,Δ) ], where ℓ^*_C(,Δ) = E_|, D,Δ[ λ^*(,Δ) ] is often available in closed form.

If the fitted and designer models are identical, then ℐ_D = 𝒥_D, and =. In this case, we obtain an asymptotic approximation to the internal expected loss given by L^*_FF(Δ) = E_| F,Δ[ ℓ_U^*(,Δ) ], where ℓ_U^*(,Δ) is the expected value of λ^*(, ,Δ) with respect to (<ref>).

















































































 §.§ Exemplar loss functions


Table <ref> shows the functions ℓ^*_G(,Δ) and ℓ^*_C(,Δ) for the exemplar loss functions considered in this paper (see Table <ref>). Derivation of these expressions is given in Section SM3 of the Supplementary Material. Also shown is the function ℓ^*_U(,Δ) for the asymptotic approximation to the internal expected loss. Note that designs that minimise the fitted prior expectation of ℓ^*_U(,Δ) under self-information and squared error are commonly referred to as pseudo-Bayesian D- and A-optimal, respectively. In Table <ref>, 𝒯̃_θθ = (ℐ̃_θθ - ℐ̃_θγℐ̃_γγ^-1ℐ̃_γθ)^-1 is the approximate marginal fitted posterior variance matrix of , under (<ref>). Additionally, 𝒮̃_θθ = I_p_θ + ℐ̃_θγℐ̃_γγ^-1ℐ̃_γγ^-1ℐ̃_γθ. 



Note that ℓ^*_G, ℓ^*_C and ℓ_U^* all feature -log|ℐ| (under self-information/entropy) or tr( ℐ^-1) (under squared error/trace variance). Minimising these terms has the effect of (approximately) maximising posterior precision. For the pseudo-Bayesian objective functions these are averaged with respect to the fitted prior. Whereas for external expected loss, they are averaged with respect to the designer prior of  through ; a function of  and the best possible estimate of  under the fitted model. The expressions for ℓ_C^* (i.e. under composite loss) give the entropy and trace of variance under the asymptotic distribution of  under a misspecified model (<ref>).

The expressions for ℓ^*_G(, Δ) under the generator squared error and self-information loss functions feature two extra terms related to the quality of  as a point estimate of the common parameters .  The first of these, ‖ - ‖_2,A^2 where A is either 𝒯̃_θθ^-1 or 𝒮̃_θθ, relates to the bias. This term penalises designs that lead to  being (on average under the designer model) far from . The second term measures the difference between the variance of  under the fitted model, i.e. 𝒯̃_θθ, and under the designer model, i.e. 𝒦̃_θθ. 



§ EXAMPLES


Two examples are considered in this section where designs are found under our proposed design framework. In both cases, the fitted model is simpler than the designer model.  As will be seen, the designs found provide robustness through replicated design points within a linear regression setting, and also through providing space filling properties for estimating parameters under a non-linear model. In Section SM6 in the Supplementary Material, a further example is provided where the fitted model is more complex (as measured by number of parameters) than the designer model.



 §.§ Normal linear regression vs full-treatment model
 

In this example we explore the implications of assuming a full-treatment designer model within a fitted linear regression setting. Under the fitted model, it is assumed

    |, σ^2, F,Δ∼N(X , σ^2 I_n),

where X is an n × p model matrix (a function of Δ), and σ^2>0 is a nuisance parameter. We assume a conjugate normal-inverse-gamma fitted prior distribution for  and σ^2 such that | σ^2, F ∼N(_F, σ^2 V_F) and σ^2 | F ∼IG(a_F/2, b_F/2). For the designer model, at this stage, it is only assumed that  has finite mean E(| D,Δ) = _D and variance var(| D,Δ) = Σ_D. Note that the fitted and designer models are disjoint, sharing no common parameters.












Under the fitted and designer models, it is possible to derive closed form expressions, to some extent, for the external expected loss. The fitted posterior distribution of  is t(_F,  b̂_F V̂_F/â_F, â_F), i.e. a multivariate t-distribution with mean _F = V̂_F (X^ + V_F^-1_F), scale matrix b̂_FV̂_F/â_F, and â_F = a_F+n degrees of freedom, where 

    V̂_F = (X^ X + V_F^-1)^-1;     b̂_F = b_F + ‖ - X _F ‖^2_2,Σ_F^-1;     Σ_F = I_n + XV_FX^ .

The external expected self-information and squared error loss functions are given by

    L_DF,SI(Δ)     =     H_SI,1 + 1/2log|V̂_F| + p/2E_| D,Δ(logb̂_F),
    
    L_DF,SE(Δ)     =    E_| D,Δ(b̂_F)/â_F-2tr(V̂_F),

respectively, where E_| D,Δ(b̂_F) = b_F + ‖_D - X _F ‖^2_2,Σ_F^-1 + tr( Σ_D Σ_F^-1) and H_SI,1 is a constant not depending on the design. Derivations of (<ref>) and (<ref>) are provided in Section SM4 in the Supplementary Material. The expectation of logb̂_F with respect to the designer model, appearing in the external expected self-information loss (<ref>), is not available in closed form. In what follows, we use a delta approximation <cit.> where E_| D,Δ(logb̂_F) ≈logE_| D,Δ(b̂_F). However, note that by Jensen's inequality, i.e. E_| D,Δ(logb̂_F) ≤logE_| D,Δ(b̂_F), the approximate external expected self-information loss (which we minimise to find designs) is actually an upper bound on the true external expected loss.

Compare (<ref>) and (<ref>) to the corresponding expressions for the internal expected self-information and squared error given by

    L_FF,SI(Δ) = H_SI,2 + 1/2log|V̂_F| ,         
    L_FF,SE(Δ) = b_F/a_F-2tr(V̂_F),

respectively, where H_SI,2 is a constant not depending on the design. The difference is that the external expected loss functions feature expectation of (a function of) b̂_F under the designer model. <cit.> describe how b̂_F summarises the inadequacy of the fitted model so it is natural that this drives the difference, with the external expected losses favouring designs for which the expected value of b̂_F is small.

Now assume a specific form for the designer model as a full-treatment model <cit.>, i.e. where each unique treatment permits a unique mean response. Mathematically, the full-treatment designer model is 

    |, σ^2, D,Δ∼N( Z , σ^2 I_n );   |σ^2, D,Δ∼N( _D, σ^2 V_D );   σ^2 | D ∼IG(a_D/2,b_D/2),

where  is a q × 1 vector of mean treatment effects. The n × q designer model matrix Z is a function of Δ where q is the number of unique treatments. For i=1,…,n, the ith row of Z is the vector _j, i.e. a vector of zeros except for the jth element which is one, indicating that run i receives unique treatment j. Therefore, _D = Z_D, Σ_D = b_D (I_n + ZV_DZ^)/(a_D-2) and E_| D,Δ(b̂_F) = b_F + ‖ Z_D - X _F ‖^2_2,Σ_F^-1 + tr( Σ_D Σ_F^-1).

Interpretation is simplified by assuming non-informative fitted and designer prior distributions where, for the fitted model b_F=0, _F = _p (the p × 1 vector of zeros) and V_F^-1 = 0_p × p (the p × p matrix of zeros), and, for the designer model, _D=_q. However, V_D needs to be finite for Σ_D to be finite and for the external expected loss to exist. It is assumed the elements of  are independent with τ_j ∼N(0, κσ^2/n_j ), where κ>0 and n_j is the number of runs receiving the jth unique treatment, for j=1,…,q. Under this prior specification, the designer posterior mean of τ_j is the sample mean of the corresponding responses shrunk toward zero by a common factor κ/(1+κ), with κ controlling the amount of shrinkage. Consequently, V_D = κ(Z^  Z)^-1.

Under these prior specifications, minimising the external expected self-information and squared error loss functions is equivalent to minimising

    L̅_DF,SI(Δ)     =     p log[(1+κ)(n - p) - κ d] - log| X^ X |, 
    L̅_DF,SE(Δ)     =    [(1+κ)(n - p) - κ d)] tr[ (X^ X)^-1]

respectively, where d=n-q is known as the pure error degrees of freedom, i.e. the number of repeated treatments. Derivations of expressions (<ref>) and (<ref>) are provided in Section SM5 of in the Supplementary Material. We refer to designs that minimise (<ref>) and (<ref>) as DE- and AE-optimal, respectively. The corresponding expressions to (<ref>) and (<ref>)  under the internal expected loss are

    L̅_FF,SI(Δ) =- log| X^ X |,         L̅_FF,SE(Δ) = tr[ (X^ X)^-1],

respectively, and designs that minimise L̅_FF,SI(Δ) and L̅_FF,SE(Δ) are called D- and A-optimal, respectively. Hence it can be seen that under the external expected loss there is a trade-off between precision of estimation of , as measured by functions of X^ X, and pure error degrees of freedom, i.e. DE- and AE-optimal designs will tend to feature more replication than D- and A-optimal designs, respectively. A numerical example verifying this is now provided.


















Consider Example 1 from <cit.> involving an experiment with n=16 runs and k=3 design variables. The fitted model is a second-order model including an intercept, three first-order terms, three quadratic terms and three pairwise interactions, i.e. p = 10. We specify κ = n leading to a unit-information prior for  <cit.>.

DE- and AE-efficiency of a design Δ are defined by

    Eff_DE(Δ)    =   exp{1/p[L̅_DF,SI(Δ^*_DF,SI) - L̅_DF,SI(Δ)] }
    Eff_AE(Δ)    =   L̅_DF,SE(Δ^*_DF,SE)/L̅_DF,SE(Δ),

respectively, where Δ^*_DF,SI and Δ^*_DF,SE are the DE- and AE-optimal designs, respectively, with similar expressions for D- and A-efficiency. The D-efficiency of the DE-optimal design is 85%, whereas the DE-efficiency of the D-optimal design is only 7%. The corresponding figures for squared error are 84% and 10%. Clearly, the D- and A-optimal designs are not robust to the external expected loss, whereas the DE- and AE-optimal designs are robust to the internal expected loss. The D-optimal design has q=16 unique treatments compared to q=10 for the DE-optimal design. The equivalent values for the A- and AE-optimal designs are q=14 and q=10, respectively. The minimum number of unique treatments to be able to estimate all elements of  under the fitted model is q=p=10. Thus the DE- and AE-optimal designs feature the maximum amount of replication possible whilst still being able to estimate .

Expressions for the external expected loss under self-information (<ref>) and squared error (<ref>) should be compared to the following objective functions under the DP- and AP-criteria <cit.> 

    L̅_DP(Δ)     =     plog F_p,d,1-α- log| X^ X |
    L̅_AP(Δ)     =     F_1,d,1-αtr[ (X^ X)^-1].

respectively. Note that L̅_DP(Δ) in (<ref>) is actually the natural logarithm of the objective function in <cit.> but the optimisation problem is equivalent due to log being a monotonically increasing function. In (<ref>) and (<ref>), F_p,d,1-α the 1-α-quantile of the F_p,d distribution with F_p,d,1-α being a decreasing function of d. The motivation behind the DP- and AP-criteria is to minimise the volume of the confidence region or sum of confidence intervals, respectively, for , where the response variance σ^2 is estimated under the full-treatment model; a model-independent estimate. In both cases, the objective functions are modified by a decreasing function of the pure error degrees of freedom, d. As shown above we are able to provide a different, Bayesian decision-theoretic, motivation for the same outcome.


































 §.§ Estimating parameters of a Michaelis-Menten model under model discrepancy
 

In this section we consider designing an experiment to estimate the parameters of a non-linear model. Finding designs for such models is widely featured in the design of experiments literature (see, e.g. ; ). Under such models, it is hypothesised that the expectation of y_i is η(, Δ_i); a mathematical model based on mechanistic or phenomenological arguments. As such the parameters  can have interpretable physical meaning. 

As a specific example, consider the Michaelis-Menten model <cit.>, commonly-used in the study of pharmacokinetics and chemical kinetics. A standard experiment has a response given by reaction velocity y_i measured for k=1 controllable variable: the scaled concentration of substrate Δ_i = x_i ∈ [0,1], for i=1,…,n. The model has 

    η(, x) = θ_1 xL/(xL+θ_2),

where θ_1>0 is the maximum velocity of the chemical reaction, θ_2>0 is the substrate concentration at which the reaction rate is half of θ_1, and L is the maximum possible concentration of substrate, which we assume is L=400.

Supposing a normal distribution, we assume the following fitted model

    |, σ^2, F, Δ∼N( (, Δ), σ^2 I_n )

where (, Δ) = (η(, x_1), …, η(, x_n))^ and σ^2 > 0 is the unknown response variance (a nuisance parameter, not of direct interest). The design problem is to choose the values of Δ = (x_1,…,x_n)^, i.e. the n concentrations.

The fitted model given by (<ref>) is assumed for pragmatic convenience. However the reasoning behind the mathematical model is usually incomplete, meaning that there can be a systematic discrepancy between the true expected value of y_i and η(, x_i). The <cit.> framework aims to model this discrepancy using a Gaussian process and we assume this for the designer model. Specifically, the designer model is assumed to be 

    |, σ^2, , D, Δ∼N( (, Δ) + , σ^2 I_n ).

where 

    |α, ρ, σ^2, D, Δ∼N( _n, σ^2 ρ C ).

In (<ref>), ρ>0, τ^2>0, and C is an n × n correlation matrix with ijth element C_ij = c(| x_i - x_j | ; α) where c( | x_i- x_j | ; α) is a correlation function depending on α>0 and the distance between x_i and x_j, with the property that c( 0 ; α) = 1. In what follows, we choose the Mátern (5/2) correlation function where

    c(d ; α) = ( 1 + d/α + d^2/3 α^2) exp( - d/α)

a commonly-used correlation function <cit.>. It follows from (<ref>) and (<ref>) that

    |, σ^2, α, ρ, D, Δ∼N[(, Δ) , σ^2 ( I_n + ρ C ) ].


We consider three designs found by minimising

    [    ;     ]

In each case, the expected loss is not available in closed form. Indeed, since the fitted posterior is not of known form, the loss function itself is analytically intractable. To approach this problem, we use a nested Monte Carlo approximation <cit.>. The following algorithm is used to approximate the external expected squared error loss.


  * Generate a sample {_b, σ^2_b, α_b, ρ_b }_b=1^B of size B from the designer prior distribution of , σ^2, α and ρ.

  * For b=1,…,B, generate _b from the designer likelihood 

    |_b, σ_b^2, α_b, ρ_b, D, Δ∼N[ (_b, Δ) , σ_b^2 ( I_n + ρ_b C_b ) ].


  * 
Approximate the external expected squared error loss by the following Monte Carlo approximation

    L̂_DF,SE(Δ) = 1/B∑_b=1^B ‖_b - Ê(|_b, F, Δ) ‖^2_2, I_2.

In (<ref>), Ê(|_b, F, Δ) is a Monte Carlo approximation to the fitted posterior mean of  conditional on _b. This approximation is formed as follows. Generate a sample {_j, σ̃_j^2 }_j=1^B̃ of size B̃ from the fitted prior distribution of  and σ^2. Then 

    Ê(|_b, F, Δ) = ∑_j=1^B̃_j π(_b |_j, σ̃_j^2, F, Δ)/∑_j=1^B̃π(_b |_j, σ̃_j^2, F, Δ),

for b=1,…,B.

Note that the sample generated in Step <ref> is reused to calculate Ê(|_b, F, Δ) for all b=1,…,B.

In a similar fashion, the external expected trace variance loss can be approximated by

    L̂_DF,TV(Δ) = 1/B∑_b=1^B tr{v̂âr̂( |_b, F, Δ) }

with

    v̂âr̂( |_b, F, Δ) = Ê(^T |_b, F, Δ) - Ê(|_b, F, Δ)Ê(|_b, F, Δ)^,

and

    Ê(^|_b, F, Δ) = ∑_j=1^B̃_j _j^π(_b |_j, σ̃_j^2, F, Δ)/∑_j=1^B̃π(_b |_j, σ̃_j^2, F, Δ).


The internal expected squared error loss is approximated by, first, generating a sample {_b, σ_b^2 }_b=1^B of size B from the fitted prior distribution of  and σ^2. For b=1, …, B, generate _b from the fitted likelihood

    |_b, σ_b^2, F, Δ∼N( (_b, Δ), σ_b^2 I_n ).

Then

    L̂_FF,SE(Δ) = 1/B∑_b=1^B ‖_b - Ê(|_b, F, Δ) ‖^2_2, I_2.


We find designs with n=20 concentrations. The fitted prior distribution is such that  and σ^2 are independent with θ_1, θ_2 iid∼U[20,200] and σ^2 ∼Exp(1); an exponential distribution with mean one. The designer prior distribution is such that , σ^2, α and ρ are independent with θ_1, θ_2 iid∼U[20,200], σ^2, ρiid∼Exp(1) and α∼Exp(5).



To minimise the nested Monte Carlo approximations to the expected loss functions, we use the approximate coordinate exchange algorithm <cit.>. We use B = B̃ = 20,000 but it is noted that if ‖ - Ê(|, F, Δ) ‖^2_2, I_2 is continuously differentiable then the convergence rate can be slightly improved <cit.> by setting B ∝B̃^2 (whilst holding BB̃ fixed).










Table <ref> shows the efficiencies of the three designs under the three expected loss functions and the number of unique concentrations. The top panel of Figure <ref> shows μ(, x) plotted against x for 100 values generated from the designer prior of  and the bottom panel shows the location of the concentrations for each of the three designs, i.e. found by minimising (i) external expected squared error loss; (ii) external expected trace variance loss; and (iii) internal expected squared error loss. The size of the plotting symbol is an increasing function of the replication of each concentration. From Figure <ref>, the three designs have common clusters of concentrations around 0.2 and 1. The designs found by minimising the external expected trace variance and internal expected squared error loss functions are broadly similar. This is confirmed in Table <ref> by the relatively high efficiencies for these designs. Compared to the other designs, the concentrations under the design that minimises the external expected squared error loss are more “spread out", i.e. a compromise between the designs found under external trace variance and internal squared error and a space-filling design.
 
The similarity of the designs found under the external expected trace variance (ii) and the internal expected squared error (iii) loss functions agrees with the asymptotic results from Section <ref>. In this example, the function d() given by (<ref>) is

    d(_θ, t_σ^2)     =    E_σ^2, ρ, α| D, Δ{E_|, σ^2, ρ, α, D, Δ[ logπ(|_θ, t_σ^2, F, Δ) ] }
        =     - n/2log (2π) - n/2log t_σ^2 - 1/2 t_σ^2‖ (, Δ) -  (_θ, Δ) ‖^2_2,I_n - n σ_0^2 ( 1 + ρ_0 )/2 t_σ^2

where _θ and t_σ^2 are the elements of  corresponding to  and σ^2, respectively, and σ^2_0 = E_σ^2 | D, Δ (σ^2) and ρ_0 = E_ρ | D, Δ (ρ) are the designer prior means of σ^2 and ρ, respectively. Clearly, d(_θ, t_σ^2) is maximised by = and σ̃^2 = σ_0^2 ( 1 + ρ_0 ). It then follows from (<ref>) that the fitted marginal posterior distribution of  is approximated by 

    N(, ℐ̃_θθ^-1),

where

    ℐ̃_θθ^-1 = σ_0^2( 1 + ρ_0 ) [ ∑_i=1^n . ∂η(, Δ_i)/∂|_ = . ∂η(, Δ_i)/∂^|_ = ]^-1.

Therefore, from Table <ref>, the external expected trace variance loss is approximated by 

    L^*_DF,TV(Δ) = E_| D, Δ[ tr( ℐ̃_θθ^-1) ].

By contrast, from Table <ref>, the internal expected squared error loss is approximated by

    L^*_FF,SE(Δ) = E_| F, Δ[ tr( ℐ_θθ^-1) ],

where

    ℐ_θθ^-1 = E_σ^2 | F, Δ (σ^2)  [ ∑_i=1^n . ∂η(, Δ_i)/∂|_ = . ∂η(, Δ_i)/∂^|_ = ]^-1.

Since σ_0^2( 1 + ρ_0 ) and E_σ^2 | F, Δ(σ^2) do not depend on the design (for this example), the approximations (<ref>) and (<ref>) are minimised by the same design. 




































































§ DISCUSSION


An extended framework is proposed for Bayesian design whereby the expectation of the loss is taken with respect to the probability distribution implied by an alternative model (the designer model) to the one being fitted to the observed responses. Properties of the framework are explored. The key findings are that the external expected loss is an upper bound on the corresponding internal expected loss under the designer model, providing justification for using a simpler fitted model for the purposes of computing a design. We also observed that designs found under composite and generator loss functions are no longer equivalent. By developing an asymptotic approximation to the external expected loss we were able to show that the composite loss functions focus solely on maximising posterior precision, whereas generator loss functions also incorporate point estimation.







We encourage researchers consider the extended framework when designing experiments under the Bayesian decision-theoretic approach, if not to specify the design to be used in practice, but at least as a benchmark to assess the performance of that design. Designs can be found under the proposed framework with only minor modification of existing computational methods (see Section <ref> and Section SM6 in the Supplementary Material).

In two of the examples, the designer model is more complex than the fitted model, and we are, in a sense, attempting to insert robustness into the design procedure. Robust Bayesian inference (i.e. inference under the M-open paradigm of , page 385) has become an important research field in recent years <cit.>. Related to this is a body of current work on the estimation of the parameters of a non-linear model (see Section <ref>) taking account of the incompleteness of the underlying mathematical model <cit.>. Considering Bayesian decision-theoretic design of experiments under these approaches will be the focal point of future work.

Lastly, we have focused on an experimental aim of parameter estimation. In principle, the extended framework can be generalised to model selection and prediction aims. Exploring this will be a further avenue for future work.



§ ACKNOWLEDGEMENT

The authors would like to thank the associate editor and two anonymous reviewers for providing comments which led to substantial improvement in the paper. They would also like to thank members of the Design Study Group at the University of Southampton for initial discussions and feedback.  JM was supported by an Australian Research Council Discovery Project (DP200101263).



  
 
 
  








§ PROOF OF LEMMA 2.1
 

Recall that the fitted and designer models are compatible, so that =. Therefore, the internal expected loss under the designer model can be written as

    L_DD(Δ)     =    E_, | D,Δ[ λ(, , D,Δ) ] 
        =    E_| D,Δ{E_|, D,Δ[ λ(, , D,Δ)]}
        =    E_| D,Δ{E_|, D,Δ[ λ_0(, _D)]},

where _D = min_∈ℬE_|, D,Δ[ λ_0(, )]. Since = _D minimises E_|, D,Δ[ λ_0(, )] over ℬ, it follows that

    E_|, D,Δ[ λ_0(, _D)] ≤E_|, D,Δ[ λ_0(, _F)],

where _F = min_∈ℬE_|, F,Δ[ λ_0(, )]. Since λ(, , F,Δ) = λ_0(, _F), it follows that

    E_|, D,Δ[ λ_0(, _D)] ≤E_|, D,Δ[ λ(, , F,Δ)],

and from (<ref>), that

    L_DD(Δ) ≤E_| D,Δ{E_|, D,Δ[ λ(, , F,Δ)]} = L_DF(Δ),

for compatible models.



§ PROOF OF THEOREM 3.1
 

The following regularity conditions are sufficient to prove Theorem 3.1.



  A1
∈ℬ where ℬ is an open, convex and bounded subset of ℝ^p.

  A2
The fitted likelihood satisfies local asymptotic normality (LAN; see, for example, , Chapter 7).

  A3
The fitted log-likelihood log(|, F, Δ) converges almost surely to d() for all ∈ℬ.

  A4
∈ℬ uniquely maximises d().

  A5
For 𝒜_ϵ = {∈ℬ : d() > d() - ϵ},

    ∫_𝒜_ϵπ(| F, Δ) d > 0,

for all ϵ >0.

  A6
The fitted prior density π(| F, Δ) is continuous and positive at .

  A7
- . ∂^2 logπ (|, F, Δ)/ ∂∂^|_ = converges in probability to ℐ̃.

  A8
. ∂logπ (|, F, Δ)/ ∂|_ = converges in distribution to N( _p, 𝒥̃).


By assumption of conditions A3-A5, Theorem 2.2 of <cit.> establishes a Bayesian consistency for the fitted posterior distribution, i.e. 

    ∫_𝔸_ϵπ(|, F, Δ) d→ 1

as n →∞. 

Then, by assumption of conditions A1 and A2, Theorem 2.1 of <cit.> shows that the fitted posterior distribution converges to (5) in total variation, as n →∞.

The external expected loss given by (2) is approximated by

    L^*_DF(Δ)    =   E_, | D,Δ{E_| , , F,Δ[ λ^*(, , Δ)] }
        =    E_| D,Δ( E_|, D,Δ{E_| , , F,Δ[ λ^*(, , Δ)] }),

by replacing the loss by the approximate loss under (5). 

Consider the inner expectation of λ^*(, , Δ) with respect to | , , F,Δ. Using properties of the multivariate normal and block matrices <cit.>, it is straightforward to derive

    |, , F,Δ∼N( P̃_1  - P̃_2 , ℐ̃_γγ^-1),

as an approximation to | , , F,Δ.

Now consider taking expectation with respect to |, D,Δ. Note that E_|, , F,Δ[ λ^*(, , Δ)] only depends on  through , thus it suffices to take expectation with respect to |, D,Δ. 

Assumption of conditions A1, A3, A4, A6-A8, and by, for example, Theorem 4.1.3 of <cit.>,

    |, D, Δ∼N( , 𝒦̃),

thus completing the proof.

LAN is a general property. However it may be non-trivial to verify for a given fitted likelihood. Instead, <cit.> provides the following conditions which are more restrictive but may be easier to verify. 



  B1
∈ℬ uniquely maximises logπ (|, F, Δ).

  B2
logπ (|, F, Δ) has continuous third derivatives, uniformly bounded on ℬ.

  B3
ℐ̃ is positive-definite.


Assumption of B1-B3, along with A1, A3, A4 and A6, and by Theorem 3.2 of <cit.>, the fitted posterior converges to (5) in total variation.





















§ DERIVATION OF EXPRESSIONS IN TABLE 2
 



 §.§ Self-information loss


For the generator self information loss,

    ℓ^*_G,SI(, Δ)     =    E_|, D, Δ[  E_|, , F, Δ[ N_SI - 1/2log|ℐ̃| + 1/2‖ - ‖^2_2,ℐ̃] ] 
        =     N_SI - 1/2log|ℐ̃|
               + 1/2E_|, D, Δ[  E_|, , F, Δ[‖ - ‖^2_2,ℐ̃] ].

The inner expectation on the right hand side of (<ref>) is

    E_|, , F, Δ[‖ - ‖^2_2,ℐ̃]     =    E_|, , F, Δ[ ‖ - ‖^2_2,ℐ̃_γγ + 2 ( - )^ℐ̃_γθ( - ) . 
              . + ‖ - ‖^2_2,ℐ̃_θθ] 
        =    ‖P̃_1  - P_2  - ‖^2_2, ℐ̃_γγ + p_γ
           + 2 (P̃_1  - P_2  - )^ℐ̃_γθ( - ) + ‖ - ‖^2_2,ℐ̃_θθ
        =    ‖ - ‖_2,𝒯̃^-1_θθ^2 + p_γ.

Substituting (<ref>) into (<ref>) we obtain

    ℓ^*_G,SI(, Δ)     =     N_SI - 1/2log|ℐ̃| + 1/2E_|, D, Δ[ ‖ - ‖^2_2,𝒯̃^-1_θθ + p_γ] 
        =     N_SI- 1/2log|ℐ̃| + 1/2‖ - ‖_2,𝒯̃^-1_θθ^2 + 1/2tr( 𝒯̃^-1_θθ𝒦̃_θθ) +p_γ/2

as required.



 §.§ Squared error loss


For the generator squared error loss

    λ^*_G,SE(, , Δ) = ‖ - ‖^2_2,I_p,

and

    ℓ^*_G,SE(, Δ) = E_|, D, Δ[  E_|, , F, Δ[ ‖ - ‖^2_2,I_p] ].

The inner expectation is

    E_|, , F, Δ[ ‖ - ‖^2_2,I_p]     =    E_|, , F, Δ[ ‖ - ‖^2_2,I_p_γ + ‖ - ‖^2_2,I_p_θ] 
        =    ‖P̃_1 - P̃_2  - ‖^2_2,I_p_γ + tr( ℐ̃_γγ^-1) + ‖ - ‖^2_2,I_p_θ
        =    ‖ - ‖^2_2,𝒮̃_θθ + tr( ℐ̃_γγ^-1).

Substituting (<ref>) into (<ref>) we obtain

    ℓ^*_G,SE(, Δ)    =   E_|, D, Δ[‖ - ‖^2_2,𝒮̃_θθ + tr( ℐ̃_γγ^-1) ] 
        =    ‖ - ‖^2_2,𝒮̃_θθ + tr( ℐ̃_γγ^-1) + tr(𝒮̃_θθ𝒦̃_θθ)
        =    ‖ - ‖^2_2,𝒮̃_θθ + tr( ℐ̃^-1) + tr[𝒮̃_θθ(𝒦̃_θθ-𝒯̃_θθ)]

as required. The last line follows from

    tr( ℐ̃^-1) = tr( ℐ̃_γγ^-1 + ℐ̃_γγ^-1ℐ̃_γθ𝒯̃_θθℐ̃_θγℐ̃_γγ^-1) + tr(𝒯̃_θθ).




 §.§ Entropy loss


The composite entropy loss is

    ℓ^*_C,E(, Δ)     =     - E_|, F, Δ[ logπ (|, F, Δ) ]
        =     N_SI - 1/2log|ℐ̃| + p/2

under (5).



 §.§ Trace variance loss


The composite trace variance loss is

    ℓ^*_C,TV(, Δ)     =    tr( var( |, F, Δ) ) 
        =    tr(  ℐ̃^-1)

under (5).



§ DERIVATION OF EXPRESSIONS (7) AND (8) IN SECTION 4.1
 

The self-information loss function is given by

    λ_SI(, , F, Δ)     =     - logπ(|, F, Δ)
        =    logΓ( â_F/2) + p/2logπ - logΓ( â_F+p/2) + p/2logb̂_F + 1/2log|V̂_F |
           + â_F + p/2log( 1 + 1/â_F(  - _F )^( b̂_F/â_FV̂_F )^-1(  - _F ) ),

where Γ(·) is the Gamma function. The expected value of λ_SI(, , F, Δ) with respect to the fitted posterior distribution of  is

    E_|, F, Δ[ λ_SI(, , F, Δ)]     =    logΓ( â_F/2) + p/2logπ - logΓ( â_F+p/2) 
           + p/2logb̂_F + 1/2log|V̂_F |
           +â_F + p/2( ψ( â_F+p/2) - ψ( p/2) ) 
        =     H_SI,1 + 1/2log|V̂_F | + p/2logb̂_F,

where (<ref>) follows from, for example, <cit.> with ψ(·) the digamma function. The expression (7) is given by taking expectation of (<ref>) with respect to the marginal distribution of  under the designer model.

The squared error loss is

    λ_SI(, , F, Δ) = ‖ - _F ‖^2_2,I_p.

From properties of the multivariate t-distribution <cit.>, the expected value of λ_SE(, , F, Δ) with respect to the fitted posterior distribution of  is

    E_|, F, Δ[ λ_SE(, , F, Δ)]     =    tr( var( |, F, Δ)) 
        =    b̂_F/â_F-2tr( V̂_F ).

The expression (8) is given by taking expectation of (<ref>) with respect to the marginal distribution of  under the designer model where

    E_| D, Δ[b̂_F]     =    E_| D, Δ[ b_F + ‖ - X _F ‖^2_2,Σ_F^-1]
        =     b_F + ‖_D - X_F ‖^2_2,Σ_F^-1 + tr(Σ_D Σ_F^-1).




§ DERIVATION OF EXPRESSIONS (9) AND (10) IN SECTION 4.1
 

Under the fitted prior distribution with V_F^-1 = 0_p × p,

    Σ_F^-1    =    (I_n + X V_F X^)^-1
        =     I_n - X(X^ X )^-1 X^
        =     I_n - H_X

where H_X = X(X^ X )^-1 X^. Under the designer prior distribution with V_D = κ(Z^ Z)^-1,

    Σ_D = b_D/a_D-2(I_n + κ H_Z ),

where H_Z = Z(Z^ Z )^-1 Z^. Then with b_F=0, _F = _p and _D = _q

    E_| D, Δ( b̂_F )     =     b_F + tr( Σ^-1_FΣ_D )
        =    b_D/a_D-2tr[ (I_n - H_X )(I_n + κ H_Z)]
        =    b_D/a_D-2tr( I_n + κ H_Z - H_X - κ H_XH_Z )
        =    b_D/a_D-2[ n + κ q - p - κtr(H_XH_Z )]
        =    b_D/a_D-2[ (1+κ)n - p - κtr(H_XH_Z ) - κ d]

with the last line following from d=n-q. Now

    tr(H_XH_Z ) = tr[H_X Z (Z^ Z)^-1 Z^] = tr[Z^ H_X Z (Z^ Z)^-1].

The q × q matrix Z^ Z is diagonal with the lth diagonal element, denoted r_l, giving the number of times treatment l=1,...,q is replicated. Note that ∑_l=1^q r_l = n. Since Z^ Z is diagonal,

    tr(H_XH_Z ) = ∑_l=1^q 1/r_l[ Z^ H_X Z]_ll,

where [ Z^ H_X Z]_ll is the lth diagonal element of Z^ H_X Z. Now

    [ Z^ H_X Z]_ll = ∑_i=1^n ∑_j=1^n z_ilz_jl h_Xij,

with h_Xij the ijth element of H_X and z_il the ilth element of Z. Note, if runs i and j receive the same treatment (say l), then h_Xij = h_Xii = h_Xjj and z_il = z_jl = 1. Then

    tr(H_XH_Z )    =   ∑_l=1^q 1/r_l∑_i=1^n ∑_j=1^n z_ilz_jl h_Xij
        =    ∑_l=1^q 1/r_l∑_i=1^n z_il∑_j=1^n z_jl h_Xij
        =    ∑_l=1^q 1/r_l∑_i=1^n z_il r_l h_Xii
        =    ∑_l=1^q ∑_i=1^n z_il h_Xii
        =    ∑_i=1^n h_Xii∑_l=1^q z_il
        =    ∑_i=1^n h_Xii
        =     p.

Then

    E_| D, Δ( b̂_F )    =   b_D/a_D-2[ (1+κ)n - p - κ p - κ d] 
        =    b_D/a_D-2[ (1+κ)(n - p) - κ d]

Lastly (<ref>) is substituted into expressions (7) and (8) ignoring constants that do not depend on Δ.



§ ESTIMATING PARAMETERS OF CUBIC SPLINE MODELS
 

Suppose the aim is to learn the relationship between a scalar controllable variable x ∈𝒳 = [0,1] and a response y. For i=1,…,n, it is assumed that

    y_i ∼N( μ(x_i), σ^2 ),

where μ(x) and σ^2 are unknown. The function μ(x) can be approximated by the Michaelis-Menten model given by (13). However, this model is known to be an incomplete representation of reality. In light of this, the fitted model is assumed to have

    μ(x) = ∑_j=1^m γ^(m)_j g^(m)_j(x),

where {g^(m)_1(x),…, g^(m)_m(x) } are a set of m basis functions and ^(m) = (γ_1^(m),…,γ_m^(m))^ are unknown parameters. We let {g^(m)_1(x),…, g^(m)_m(x) } be cubic spline basis functions <cit.>. Instead of fixing m we allow m to be unknown with the constraint 4 ≤ m ≤ n (the lower bound following from the assumption of cubic splines).

The design problem is to choose the values of Δ = (x_1,…,x_n)^.

The fitted prior distributions for ^(m) and σ^2 are ^(m)|σ^2, F, Δ∼N( _m, σ^2 κ I_m ) and σ^2 | F, Δ∼IG(a_F/2, b_F/2 ), where κ is known. The fitted prior distribution for the discrete number of basis functions, m, is such that π(m) = 1/(n-5), i.e. a uniform prior.

The model-averaged fitted posterior mean of μ(x) is

    E( μ(x) |, F, Δ) = ∑_m=4^n [π(m |, F, Δ) ∑_j=1^m g_j(x) E( γ_j^(m)|, m, F, Δ)].

In (<ref>), E( γ_j^(m)|, m, F, Δ) is the jth element of the fitted posterior mean of ^(m) (conditional on m basis functions, for j=1,…,m) given by

    E( ^(m)|, m, F, Δ) = ( κ^-1 I_m + G_m^ G_m )^-1 G_m^,

where G_m is the n × m matrix with jth column (g^(m)_j(x_1), …, g^(m)_j(x_n) )^, for j=1,…,m. Additionally, π(m |, F, Δ) is the fitted posterior probability of m given by

    π(m |, F, Δ) ∝| R_m |^-1/2( b_F + ‖‖^2_2, R_m)^- a_F + n/2,

where R_m = I_n + κ G_mG_m^. The expressions (<ref>) and (<ref>) follow from the fact that the fitted model is a normal linear model and can be found in, for example, <cit.>.

The chosen loss function is predictive squared error 

    λ(, μ(·) , F, Δ) = ∫_𝒳( μ(x) - E( μ(x) |, F, Δ) )^2 dx,

a measure of the accuracy of the model-averaged fitted posterior mean of μ(x).

We consider finding designs by minimising


  (a)
internal expected predictive squared error loss; and

  (b)
external expected predictive squared error loss.

For (b) external expected predictive squared error loss, the disjoint designer model has μ(x) = η(, x) where η(·, x) is given by (13) with independent designer prior distributions of ξ_1 ∼U[20,200], ξ_2 ∼U[20,200] and σ^2 ∼IG(a_D/2, b_D/2). 







Neither of the expected loss functions are available in closed form and require approximation for which we use Monte Carlo. The internal expected predictive squared error loss can be approximated using the following steps.


  * For b=1, …, B, generate m_b uniformly from {4, …, n }, σ_b^2 from IG(a_F/2,b_F/2), ^(m_b) from N(_m_b, σ_b^2 κ I_m_b). Then generate _b from N(G_m_b^(m_b), σ_b^2 I_n).

  * Approximate the internal expected predictive squared error loss by

    L̂_FF,PSE(Δ) = 1/B∑_b=1^B λ[_b, ∑_j=1^m_b g_j^(m_b)(·) γ_bj^(m_b) , F, Δ]

where γ_bj^(m_b) is the jth element of ^(m_b).

The external expected predictive squared error loss can be approximated using the following steps.


  * For b=1, …, B, generate _b from the designer prior distribution and σ_b^2 from IG(a_D/2, b_D/2).Then generate _b from N((_b, Δ), σ_b^2 I_n).

  * Approximate the internal expected predictive squared error loss by

    L̂_DF,PSE(Δ) = 1/B∑_b=1^B λ[_b, η(_b, ·) , F, Δ].


For both approximations, the integration in (<ref>) is evaluated using quadrature.

The hyperparameters of the fitted and designer prior distributions of σ^2 are set as a_F=a_D = 6 and b_F = b_D = 4, ensuring that the prior mean and variance are both one. The value of κ is chosen to be 10^6, indicating a very vague fitted prior distribution for ^(m).
Designs are found with n=10 using the approximate coordinate exchange algorithm.

Table <ref> shows the efficiencies of the two designs under internal and external predictive squared error loss. Figure <ref> shows the values of Δ = (x_1,…,x_n)^ for the two designs (where the size of the plotting symbol is an increasing function of the replication.). The design found by minimising internal expected loss is similar to an evenly-spaced design whereas the design found by minimising external expected loss appears to be a compromise between the evenly-spaced design and the structure of points found for the Michaelis-Menten model in Section 4.2. Further, it was found that the internal expected loss was flat near the design minimising its value, evidenced by the high efficiency under the internal expected loss for the design minimising the external expected loss. This agrees with the observation of <cit.> that expected loss functions become flat under vague prior distributions (i.e. in this example κ = 10^6). By contrast, the design found by minimising the external expected loss provides a design of high efficiency for the fitted model but reflects best current knowledge of the underlying process.





rss



