$ D $-optimal designs for Poisson regression with synergetic interaction effect

We characterize $D$-optimal designs in the two-dimensional Poisson regression model with synergetic interaction and provide an explicit proof. The proof is based on the idea of reparameterization of the design region in terms of contours of constant intensity. This approach leads to a substantial reduction of complexity as properties of the sensitivity can be treated along and across the contours separately. Furthermore, some extensions of this result to higher dimensions are presented.


Introduction
Count data plays an important role in medical and pharmaceutical development, marketing, or psychological research. For example, Vives, Losilla, and Rodrigo [21] performed a review on articles published in psychological journals in the period from 2002 to 2006. There they found out that a substantial part of these articles dealt with count data for which the mean was quite low (for details we refer to the discussion in Graßhoff et al. [8]). In these situations, standard linear models are not applicable because they cannot account for the inherent heteroscedasticity. Instead Poisson regression models are often more appropriate to describe such data. As an early source in psychological research we may refer to the Rasch Poisson counts model introduced by Rasch [15] in 1960 to predict person ability in an item response setup. The Poisson regression model can be considered as a particular Generalized Linear Model (see McCullagh and Nelder [13]). For the analysis of count data in the Poisson regression model there is a variety of literature (see e. g. Cameron and Trivedi [3]) and the statistical analysis is implemented in main standard statistical software packages (cf. "glm" in R,"GENLIN" in SPSS, "proc genmod" in SAS), But only few work has been done to design such experiments. Ford, Torsney and Wu derived optimal designs for the one-dimensional Poisson regression model in their pioneering paper on canonical transformations [7]. Wang et al. [22] obtained numerical solutions for optimal designs in two-dimensional Poisson regression models both for the main effects only (additive) model as well as for the model with interaction term. For the main effects only model the optimality of their design was proven analytically by Russell et al. [17] even for larger dimensions. Rodríguez-Torreblanca and Rodríguez-Díaz [16] extended the result by Ford et al. for one-dimensional Poisson regression to overdispersed data specified by a negative binomial regression model, and Schmidt and Schwabe [18] generalized the result by Russell et al. for higher-dimensional Poisson regression to a much broader class of additive regression models. Graßhoff et al. [8] gave a complete characterization of optimal designs in an ANOVA-type setting for Poisson regression with binary predictors and Kahle et al. [9] indicate, how interactions could be incorporated in this particular situation.
In this paper, we find D-optimal designs for the two-dimensional Poisson regression model with synergetic interaction as before considered numerically by Wang et al. [22]. We show the D-optimality by reparametrizing the design space via hyperbolic coordinates, such that the inequalities in the Kiefer-Wolfowitz equivalence theorem only need to be checked on the boundary and the diagonal of the design region. This allows us to find an analytical proof for the D-optimality of the proposed design. Furthermore, we extend this result in various ways to higher-dimensional Poisson regression. First, we find D-optimal designs for first-order and second-order interactions, given that the prespecified interaction parameters are zero. Second, we present a D-optimal design for Poisson regression with first-order synergetic interaction where the design space is restricted to the union of the two-dimensional faces of the positive orthant.
The paper is organized as follows. In the next section we introduce the basic notations for Poisson regression models and specify the corresponding concepts of information and design in Section 3. Results for two-dimensional Poisson regression with interaction are established in Section 4. In Section 5, we present some extensions to higher-dimensional Poisson regression models. Further extensions are discussed in Section 6. Technical proofs have been deferred to an Appendix. We note that most of the inequalities there have first been detected by using the computer algebra system Mathematica [23], but analytical proofs are provided in the Appendix for the readers' convenience.

Model Specification
We consider the Poisson regression model where observations Y are Poisson distributed with intensity E(Y ) = λ(x) which depends on one or more explanatory variables x = (x 1 , ..., x k ) in terms of a generalized linear model. In particular, we assume a log-link which relates the mean λ(x) to a linear component .., f p (x)) is a vector of p known regression functions and β is a p-dimensional vector of unknown parameters. For example, if x = x is one-dimensional (k = 1), then simple Poisson regression is given by f (x) = (1, x) with p = 2, β = (β 0 , β 1 ) and intensity λ(x) = exp(β 0 + β 1 x). For two explanatory variables x = (x 1 , x 2 ) (k = 2) multiple Poisson regression without interaction is given by f (x) = (1, x 1 , x 2 ) with p = 3, β = (β 0 , β 1 , β 2 ) and intensity In what follows we will focus on the two-dimensional multiple regression (x = (x 1 , x 2 ), k = 2) with interaction term, where p = 4, f (x) = (1, x 1 , x 2 , x 1 x 2 ) , β = (β 0 , β 1 , β 2 , β 12 ) and intensity Here β 0 is an intercept term such that the mean is exp(β 0 ) when the explanatory variables are equal to 0. The quantities β 1 and β 2 denote the direct effects of each single explanatory variable, and β 12 describes the amount of the interaction effect when both explanatory variables are active (non-zero).
Typically the explanatory variables describe non-negative quantities (x 1 , x 2 ≥ 0) like doses of some chemical or pharmaceutical agents -or difficulties of tasks in item response experiments in psychology. In particular, in the latter case the expected number of counts (correct answers) decreases with increasing difficulty. Then it is reasonable to assume that the direct effects are negative (β 1 , β 2 < 0), and that the interaction effect tends into the same direction if present (β 12 ≤ 0). In the case that β 12 < 0 this will be called a synergy effect because it describes a strengthening of the effect if both components are used simultaneously.

Information and Design
In experimental situations the setting x of the explanatory variables may be chosen by the experimenter from some experimental region X . As the explanatory variables describe non-negative quantities, and if there are no further restrictions on these quantities, it is natural to assume that the design region X is the non-negative half-axis [0, ∞) or the closure of quadrant I in the Cartesian plane, [0, ∞) 2 , in one-or twodimensional Poisson regression, respectively.
To measure the contribution of an observation Y at setting x the corresponding information can be used: With the log-link the Poisson regression model constitutes a generalized linear model with canonical link [13]. Furthermore for Poisson distributed observations Y the variance and the mean coincide, Var(Y ) = E(Y ) = λ(x). Hence, according to [2] the elemental (Fisher) information for an observation Y at a setting x is a p × p matrix given by Note that on the right-hand side the intensity λ(x) = exp(f (x) β) depends on the linear component f (x) β and, hence, on the parameter vector β. Consequently also the information depends on β as indicated by the notation M β .
For N independent observations Y 1 , ..., Y N at settings x 1 , ..., x N the joint Fisher information matrix is obtained as the sum of the elemental information matrices, The collection x 1 , ..., x N of settings is called an exact design, and the aim of design optimization is to choose these settings such that the statistical analysis is improved. The quality of a design can be measured in terms of the information matrix because its inverse is proportional to the asymptotic covariance matrix of the maximum-likelihood estimator of β, see Fahrmeir and Kaufmann [4]. Hence, larger information means higher precision. However, matrices are not comparable in general. Therefore one has to confine oneself to some real valued criterion function applied to the information matrix. In accordance with the literature we will use the most popular D-criterion which aims at maximizing the determinant of the information matrix. This criterion has nice analytical properties and can be interpreted in terms of minimization of the volume of the asymptotic confidence ellipsoid for β based on the maximum-likelihood estimator. The optimal design will depend on the parameter vector β and is, hence, only locally optimal.
Finding an optimal exact design is a discrete optimization problem which is often too hard for analytical solutions. Therefore we adopt the concept of approximate designs in the spirit of Kiefer [10]. An approximate design ξ is defined as a collection x 0 , ..., x n−1 of n mutually distinct settings in the design region X with corresponding weights w 0 , ..., w n−1 ≥ 0 satisfying n−1 i=0 w i = 1. Then an exact design can be written as an approximate design, where x 0 , ..., x n−1 are the mutually distinct settings in the exact design with corresponding numbers N 0 , ..., N n−1 of replications, n−1 i=0 N i = N , and frequencies w i = N i /N , i = 0, ..., n − 1. However, in an approximate design the weights are relaxed from multiples of 1/N to non-negative real numbers which allow for continuous optimization.
For an approximate design ξ the information matrix is defined as which therefore coincides with the standardized (per observation) information matrix 1 N M β (x 1 , ..., x N ). An approximate design ξ * will be called locally D-optimal at β if it maximizes the determinant of the information matrix M β (ξ).
In the case of two-dimensional Poisson regression without interaction the design ξ * β 1 ,β 2 which assigns equal weights w * 0 = w * 1 = w * 2 = 1/3 to the three settings x * 0 = (0, 0), x * 1 = (2/|β 1 |, 0), and x * 2 = (0, 2/|β 2 |) is locally D-optimal at β on X = [0, ∞) 2 for β 1 , β 2 < 0, see Russell et al. [17]. Note that the optimal coordinates on the axes coincide with the optimal values in the one-dimensional case, see Schmidt and Schwabe [18]. In both cases the optimal design is minimally supported, i.e. the number n of support points of the design is equal to the number p of parameters. It is well-known that for D-optimal minimally supported designs the optimal weights are all equal, w * i = 1/p, see Silvey [20]. Such optimal designs are attractive as they can be realized as exact designs when the sample size N is a multiple of the number of parameters p.
Further note that these optimal designs always include the setting x 0 = 0 or x 0 = (0, 0), respectively, where the intensity λ attains its largest value.
The above findings coincide with the numerical results obtained by Wang et al. [22] who also numerically found minimally supported D-optimal designs for the case of twodimensional Poisson regression with interaction. In what follows we will give explicit formulae for these designs and establish rigorous analytical proofs of their optimality.
We start with the special situation of vanishing interaction (β 12 = 0). In this case standard methods of factorization can be applied to establish the optimal design, see Schwabe [19], section 4.
In contrast to the result of Theorem 4.1 the intensity fails to factorize in the case of a non-vanishing interaction (β 12 = 0). Thus a different approach has to be chosen. As a prerequisite we mention that in the above cases the optimal designs can be derived from those for standard parameter values β 0 = 0 and β 1 = −1 in one dimension or β 1 = β 2 = −1 in two dimensions by canonical transformations, see Ford et al. [7], or, more generally, by equivariance considerations, see Radloff and Schwabe [14]. We will adopt this approach also to the two-dimensional Poisson regression model with interaction and consider the case β 0 = 0 and β 1 = β 2 = −1 first. There the interaction effect remains a free parameter, and we denote the strength of the synergy effect by ρ = −β 12 ≥ 0.  [22] we consider a class Ξ 0 of minimally supported designs as potential candidates for being optimal. In the class Ξ 0 the designs have one setting at the origin x 0 = (0, 0), where the intensity is highest, one setting x 1 = (x 1 , 0) and x 2 = (0, x 2 ) on each of the bounding axes of the design region as for the optimal design in the model without interaction, and an additional setting x 3 = (t, t) on the diagonal of the design region, where the effects of the two components are equal. The following result is due to Könner [12].
Note that t = 2 for ρ = 0 is in accordance with the optimal product-type design in Theorem 4.1, t is continuously decreasing in ρ, and t tends to 0 when the strength of synergy ρ gets arbitrarily large. Figure 1 shows the value of t in dependence on ρ.
To establish that ξ t is locally D-optimal within the class of all designs on X we will make use of the Kiefer-Wolfowitz equivalence theorem [11] in its extended version incorporating intensities, see Fedorov [6]. For this we introduce the sensitivity function where we suppress the dependence on β in the notation. Then by the equivalence theorem a design ξ * is (locally) D-optimal if (and only if) the sensitivity function ψ(x; ξ * ) does not exceed the number p of parameters uniformly on the design region X . Equivalently we may consider the deduced sensitivity function To establish this condition we need some preparatory results on the shape of the (deduced) sensitivity function. Figure 2 shows d(x; ξ t ) for t = 2 for ρ = 0, i.e. for the standardized setting in Theorem 4.1.
Note that ξ t is invariant with respect to the permutation of x 1 and x 2 . Then, combining Lemmas 4.3 to 4.5, we obtain d(x; ξ t ) ≤ 0 for all x ∈ X which establishes the D-optimality of ξ t in view of the equivalence theorem.
Theorem 4.6. In the two-dimensional Poisson regression model with interaction the design ξ t is locally D-optimal at β = (0, −1, −1, −ρ) on X = [0, ∞) 2 which assigns equal weights 1/4 to the 4 settings x 0 = (0, 0), For the general situation of decreasing intensities (β 1 , β 2 < 0) and a synergy effect (β 12 < 0) the optimal design can be obtained by simultaneous scaling of the settings Schwabe [14]. This simultaneous scaling leaves the linear component and, hence, the intensity unchanged, f (x) β = f (x) β. If the scaling of x is applied to the settings in ξ t of Theorem 4.6, then the resulting rescaled design will be locally D-optimal at β on X as the design region is invariant with respect to scaling. Furthermore, the design optimization is not affected by the value β 0 of the intercept term because this term contributes to the intensity and, hence, to the information matrix only by a multiplicative factor, λ( We thus obtain the following result from Theorem 4.6.
Theorem 4.7. Assume the two-dimensional Poisson regression model with interaction Then the design which assigns equal weights 1/4 to the 4 settings x 0 = (0, 0), Note that the settings x 0 , x 1 , and x 2 of the locally D-optimal design ξ t in the model with interaction coincide with those of the optimal design for the model without interaction. Only a fourth setting x 3 = (t/|β 1 |, t/|β 2 |) has been added in the interior of the design region.

Higher-dimensional Models
In the present section on k-dimensional Poisson regression with k explanatory variables (x = (x 1 , x 2 , ..., x k ), k ≥ 3) we restrict to the standardized case with zero intercept (β 0 = 0) and all main effects β 1 = ... = β k equal to −1 for simplicity of notation. Extensions to the case of general β 0 and β 1 , ..., β k < 0 can be obtained by the scaling method used for Theorem 4.7.
Schmidt and Schwabe [18] more generally proved that in models without interactions the locally D-optimal design points coincide with their counterparts in the marginal one-dimensional models. This approach will be extended in Theorems 5.2 and 5.4 to two-and three-dimensional marginals with interactions.
In what follows we mainly consider the particular situation that all interactions occurring in the models have values equal to 0 and that the design region is the full orthant X = [0, ∞) k . Setting the interactions to zero does not mean that we presume to know that there are no interactions in the model. Instead we are going to determine locally optimal designs in models with interactions which are locally optimal at such β for which all interaction terms attain the value 0.
We start with a generalization of Theorem 4.1 to a k-dimensional Poisson regression model with complete interactions where the number of parameters is p = 2 k .
where the number of parameters is p = 1 + k + k(k − 1)/2. , and For illustrative purposes we specify this result for k = 3 components.  The optimal design points of Corollary 5.3 are visualized in Figure 3. Note that in in the Poisson regression model with first-order interactions the locally D-optimal design has only support points on the axes and on the diagonals of the faces, but none in the interior of the design region, and that the support points on each face coincide with the optimal settings for the corresponding two-dimensional marginal model. Thus only those settings are included from the full factorial {0, 2} k of the complete interaction case (Theorem 5.1) which have, at most, two non-zero components, and the locally D-optimal design concentrates on settings with higher intensity. This is in accordance with the findings for the Poisson regression model without interactions, where only those settings will be used which have, at most, one non-zero component, and carries over to higher-order interactions. In particular, for the Poisson regression model with second-order interactions where the number of parameters is p = 1 + k + k(k − 1)/2 + k(k − 1)(k − 2)/6, we obtain a similar result.
Theorem 5.4. In the k-dimensional Poisson regression model with second-order interactions the minimally supported design which assigns equal weights 1/p to the p = 1 + k + k(k − 1)/2 + k(k − 1)(k − 2)/6 settings x 0 = (0, 0, ..., 0), x 1 = (2, 0, ..., 0), x 2 = (0, 2, ..., 0), ..., x k = (0, ..., 0, 2), The proofs of Theorems 5.2 and 5.4 are based on symmetry properties which get lost if one or more of the interaction terms are non-zero. However, if only few components of x may be active (non-zero), then locally D-optimal designs may be obtained in the spirit of the proof of Lemma 4.4 for synergetic interaction effects. We demonstrate this in the setting of first-order interactions ρ ij = −β ij ≥ 0, when the design region X consists of the union of the two-dimensional faces of the orthant, i. e. when, at most, two components of x can be active.
This result follows as in the proof of Lemma 4.4. We believe that the D-optimality of the design in Theorem 5.5 could also hold on the whole positive orthant if we assume that the prespecified interaction parameters are identical and non-positive. A proof of this statement should follow in the spirit of Farrell et al. [5], similar to the constructions in the Lemmas 4.3 and 4.5 and the proof of Theorem 5.2.
However, in the situation of general synergy effects an analogon to Lemma 4.3 cannot be established because of the lacking symmetry. Hence, it remains open whether the design of Theorem 5.5 retains its optimality in the general setting.

Discussion
The main purpose of the present paper is to characterize locally D-optimal designs explicitly for the two-dimensional Poisson regression model with interaction on the unbounded design region of quadrant I when both main effects as well as the interaction effect are negative, and to present a rigorous proof for their optimality. Obviously the designs specified in Theorem 4.7 remain optimal on design regions which are subsets of quadrant I and cover the support points of the respective design. For example, if the design region is a rectangle, X = [0, b 1 ] × [0, b 2 ], then the design of Theorem 4.7 is optimal as long as b 1 ≥ 2/|β 1 | and b 2 ≥ 2/|β 2 | for the two components. Furthermore, if the design region is shifted, X = [a 1 , ∞) × [a 2 , ∞) or a sufficiently large subregion of that, then also the locally D-optimal design is shifted accordingly and assigns equal  1 , a 2 ), x 1 = (a 1 + 2/|β 1 |, a 2 ), x 2 = (a 1 , a 2 + 2/|β 2 |), and x 3 = (a 1 + t/|β 1 |, a 2 + t/|β 2 |) where t is defined as in Theorem 4.7.
Although the locally D-optimal designs only differ in the location of the support point on the diagonal, if the main effects are kept fixed, they are quite sensitive with respect to the strength ρ of the synergy parameter in their performance. The quality of their performance can be measured in terms of the local D-efficiency which is defined for a design ξ, where ξ * β denotes the locally D-optimal design at β. This efficiency can be interpreted as the asymptotic proportion of observations required for the locally D-optimal ξ * β to obtain the same precision as for the competing design ξ of interest. For example, in the standardized case of Subsection 4.1 the design ξ x would be locally D-optimal when the strength of synergy would be (2 − x)/x 2 . Its local D-efficiency can be calculated as eff D (ξ, β) = (x/t) exp((2t + ρt 2 − 2x − ρx 2 )/4) when ρ is the true strength of synergy and t is the corresponding optimal coordinate on the diagonal (t = ( √ 1 + 8ρ − 1)/(2ρ) for ρ > 0 and t = 2 for ρ = 0). For selected values of x the local D-efficiencies are depicted in Figure 4. The appealing product-type design ξ 2 of Theorem 4.1 rapidly loses efficiency if the strength ρ of synergy substantially increase. The triangular design ξ 1 seems to be rather robust over a wide range of strength parameters, while for smaller x the design ξ x loses efficiency when there is no synergy effect (ρ = 0). Hence, it would be desirable to determine robust designs like maximin D-efficient or weighted ("Bayesian") optimal designs (see e. g. Atkinson et al. [1]), but this would go beyond the scope of the present paper.
If in contrast to the situation of Theorems 4.6 and 4.7 there is an antagonistic interaction effect which means that β 12 is positive (ρ < 0), no optimal design will exist on quadrant I because the determinant of the information matrix becomes unbounded. However, if we restrict the design region to a rectangle one may be tempted to extend the above results. For example, in the standardized case (β 1 = β 2 = −1) on a square design region Lemma 4.2 may be extended as follows Moreover, Lemma 4.4 does not depend on ρ and, if, additionally, b ≤ 1/|ρ|, then the argumentation in the proof of Lemma 4.3 can be adopted, where now the hyperbolic coordinate system is centered at (1/|ρ|, 1/|ρ|) and v is negative (cf. the proof below). However, the inequalities of Lemma 4.5 are no longer valid, in general. In particular, for ρ less than, but close to −1/8 the (deduced) sensitivity function of the design ξ t shows a local minimum at t rather than a maximum which disproves the optimality of ξ t within the class of all designs on X = [0, b] 2 . In that case an additional fifth support point is required on the diagonal, and also the weights have to be optimized. So, in the case of an antagonistic interaction effect no general analytic solution can be expected and the numerically obtained optimal designs may become difficult to be realized as exact designs.
For even smaller design regions (b < 2) design points on the adverse boundaries (x 1 = b or x 2 = b) may occur in the optimal designs, but not in the interior besides the diagonal, both in the synergetic as well as in the antagonistic case.
It seems more promising to extend the present results to negative binomial (Poisson-Gamma) regression which is a popular generalization of Poisson regression which can cope with overdispersion as in Rodríguez-Torreblanca and Rodríguez-Díaz [16] for onedimensional regression or in Schmidt and Schwabe [18] for multidimensional regression without interaction. This will be object of further investigation. matrix and by the (n × n)-dimensional diagonal matrices Λ = diag(λ(x 0 ), ..., λ(x n−1 )) and W = diag(w 0 , ..., w n−1 ) the intensity and the weight matrix, respectively. Then the information matrix can be written as For minimally supported designs the matrices F, W and Λ are quadratic (p × p) and the determinant of the information matrix factorizes, det(M(ξ)) = det(W) det(Λ) det(F) 2 .
As W and Λ are diagonal and is a triangular matrix for ξ ∈ Ξ 0 , the determinants of these matrices are the products of their entries on the diagonal. Hence, ) and the weights as well as the single settings can be optimized separately. As for all minimally supported designs the optimal weights are all equal to 1/p which is here 1/4. The contribution x 2 j exp(−x j ) of the axial points is the same as in the corresponding marginal one-dimensional Poisson regression model with β j = (0, −1) and is optimized by x j = 2, j = 1, 2. Finally, t 4 exp(−2t − ρt 2 ) is maximized by t = ( √ 1 + 8ρ − 1)/(2ρ) for ρ > 0 and t = 2 for ρ = 0.
Proof of Lemma 4.3. The main idea behind this proof is to consider the deduced sensitivity function on contours of equal intensities. For this we reparametrize the design region and use shifted and rescaled hyperbolic coordinates, where v = (1 + ρx 1 )(1 + ρx 2 ) is the (shifted and scaled) hyperbolic distance and u = log( (1 + ρx 1 )/(1 + ρx 2 )) is the (shifted and scaled) hyperbolic angle in the case ρ > 0. The design region X = [0, ∞) 2 is covered by v ≥ 1 and |u| ≤ log(v). With these coordinates, fixing v > 1 returns a path parametrized in u which intersects the diagonal at u = 0. On each of these paths the intensity function λ(x) is constant.
Because ξ t is invariant under permutation of x 1 and x 2 , i. e. sign change of u, the deduced sensitivity function d(x; ξ t ) is symmetric in u, and we only have to consider the non-negative branch, 0 ≤ u ≤ log(v). Using cosh(2u) = 2 cosh 2 (u) − 1, we observe that d(x; ξ t ) is a quadratic polynomial in cosh(u) = (exp(u) + exp(−u))/2 on each path. Further, by the invariance of ξ t , the information matrix and, hence, its inverse i. e. x 1 = 0 or x 2 = 0). As the paths cover the whole design region, the statement of the Lemma follows for ρ > 0.
In the case ρ = 0 the contours of equal intensities degenerate to straight lines, where x 1 + x 2 is constant. Then the design region can be reparametrized by x 1 = v + u and Using similar arguments as for the case ρ > 0 we can show that the sensitivity function restricted to each of these line segments for v fixed is a symmetric polynomial in u of degree 4 with positive leading term. Hence, also in the case ρ = 0 the maximum of the sensitivity function can only be attained on the diagonal (u = 0) or on the boundary (|u| = v) which completes the proof.
Proof of Lemma 4.4. With the notation in the Proof of Lemma 4.2 the deduced sensitivity function can be written as and similarly for the deduced sensitivity function d 1 (x; ξ * −1 ) of the locally D-optimal design ξ * −1 in the one-dimensional marginal model when β 1 = (0, −1) . For settings x = (x 1 , 0) we then obtain d(x; ξ t ) = d 1 (x 1 ; ξ * −1 ) by the relation between the quantities and matrices in both models and their special structure. As ξ * −1 is D-optimal in the marginal model, its deduced sensitivity d 1 is bounded by zero by the equivalence theorem. Hence, we obtain d((x 1 , 0); ξ t ) ≤ 0 for all x 1 ≥ 0.
For reasons of symmetry we also get d((0, x 2 ); ξ t ) ≤ 0 for all x 2 ≥ 0 which completes the proof.
Proof of Lemma 4.5. First note that the relation between ρ and t = ( √ 1 + 8ρ−1)/(2ρ) is one-to-one such that conversely ρ = (2 − t)/t 2 . Then, with the transformation q = x/t, the inequality to show in Lemma 4.5 can be equivalently reformulated to by using (A.1). To prove the Lemma it is then sufficient to show that the inequality (A.2) holds for all 0 ≤ t ≤ 2 and all q ≥ 0. The idea behind the proof is to split the above function into a polynomial h 0 (q, t) = 1 2 exp(2)t 2 (q − 1) 2 q 2 + (q − 1) 2 (q(t − 1) − 1) 2 in t and q and a function involving the exponential terms such that d(x; ξ t ) = h 0 (q, t) − h 1 (q, t) and to find a suitable separating function h 2 (q, t) such that the inequalities h 0 (q, t) ≤ h 2 (q, t) and h 2 (q, t) ≤ h 1 (q, t) are easier to handle, where essentially methods for polynomials can be used for the former inequality while in the latter properties of exponential functions can be employed. This function h 2 (q, t) will be defined piecewise in q by h 2 (q, t) = 1 for q ≤ q 0 exp(t + 2)(q − 1) 2 q 2 for q > q 0 ,  Figure 6. The functions h 0 (q, t) (blue), h 1 (q, t) (orange) and h 2 (q, t) (green) for t = 0, 1/2, 1 and 2 where q 0 = 3/5, and the proof will be performed case-by-case. Figure 6 visualizes this approach for selected values of t.
Next we consider h 1 (q, t) as a function of t. Its partial derivative with respect to t is given by If we compare the exponential terms, we see that for all 0 ≤ t ≤ 2 uniformly in q. Hence, the partial derivative (A.3) is non-negative if To see this we notice ∂ ∂q q −3 (2 − q) exp(4(q − 1)) = −2q −4 (2q 2 − 5q + 3) exp(4(q − 1)) ≤ 0 for q ≤ 1 such that the expression on the left hand side of (A.5) attains its minimum at q = 1, where it is equal to 1. Combining the above results we obtain that h 1 (q, t) attains its minimum at t = 0 for all q ≤ 1. It remains to show that h 1 (q, 0) = exp(2q 2 ) − exp(2)q 4 ≥ 1 for all q ≤ q 0 . For this we check the derivative with respect to q which is positive for 0 < q < q 2 and negative for q 2 < q ≤ q 0 . where q 2 ≈ 0.451. Hence, evaluating h 1 (q, 0) a the end-points of the relevant interval, h 1 (0, 0) = 1 and h 1 (q 0 , 0) ≈ 1.097, we get h 1 (q, 0) ≥ 1 which finally implies h 0 (q, t) ≤ 1 ≤ h 1 (q, t) for all q ≤ q 0 and all 0 ≤ t ≤ 2. For the case q > q 0 the condition h 0 (q, t) ≤ h 2 (q, t) is equivalent to By the exponential series expansion, exp(t) ≥ 1 + t + t 2 /2 for t ≥ 0, the right hand side is bounded from below by (t + 1)q 2 exp(2), and for (A.6) to hold it is sufficient to show The derivative of this expression with respect to q equals for q ≥ 1/2 and all 0 ≤ t ≤ 2. Hence, the expression in (A.7) itself is bounded from below by its value at q 0 = 3/5, which is approximately 0.1001.
Hence, for q ≥ 0 the expression in (A.8) attains its maximum at q = 1, where it is equal to 1. This implies h 2 (q, t) ≤ h 1 (q, t) for all q > q 0 and all 0 ≤ t ≤ 2 which completes the proof.
Proof of Theorem 5.2. Here we only give a sketch of the proof. As in the Proof of Lemma 4.3 we see that the paths of equal intensity constitute hyper-planes intersecting the design region at equilateral simplices. On each straight line within these simplices the sensitivity function is a polynomial of degree four with positive leading term. Hence, following the idea of the proofs in Farrell et al. [5] we can conclude by symmetry considerations with respect to permutation of the entries in x we can conclude that the sensitivity function may attain a maximum in the interior of the design region only at the diagonal, where all entries in x are equal (x 1 = x 2 = ... = x k = x) and in the relative interior of each j-dimensional face of the design region on the respective diagonal, where all the j non-zero entries of x are equal to some x, 2 ≤ j ≤ k.
Similar to the Proof of Lemma 4.4 on each face the deduced sensitivity function is equal to its counterpart for the D-optimal design in the two-dimensional marginal model on that face and is, thus, bounded by 0.
Finally, to derive the deduced sensitivity function on the diagonals we specify the essential design matrix F and its inverse where A = diag(1, 2 1 k , 4 1 C(k,2) ) is a diagonal matrix related to the product of the non-zero coordinates of the design points, 1 m is a m-dimensional vector with all entries equal to 1, I m is the m × m identity matrix, C(m, n) denotes binomial coefficient m n , and S 2 is the incidence matrix of a balanced incomplete block design (BIBD) for k varieties and all C(k, 2) blocks of size 2. Then by (A.1) the deduced sensitivity function equals (C(j, 2)q 2 − jq + 1) 2 + j exp(2)((j − 1)q 2 − q) 2 + C(j, 2) exp(4)q 4 − exp(2jq) on the diagonals of all j-dimensional faces, j < k, and the interior diagonal for j = k, where q = x/2 as in the Proof of Lemma 4.5. By using Mathematica and a power series expansion of order 5 for the term exp(2kq) the above expression can be seen not to exceed 0 for all q ≥ 0 which establishes the local D-optimality in view of the equivalence theorem.
Proof of Theorem 5.4. The proof goes along the lines of the Proof of Theorem 5.2. The essential design matrix F and its inverse are specified as where now A = diag(1, 2 1 k , 4 1 C(k,2) , 8 1 C(k,3) ), S 3 is the incidence matrix of a BIBD for k varieties and all C(k, 3) blocks of size 3, and S 23 is the (generalized) C(k, 3) × C(k, 2) incidence matrix which relates all blocks of size 2 to those blocks of size 3 ln which their components are included. Then the deduced sensitivity function equals (C(j, 3)q 3 − C(j, 2)q 2 + jq − 1) 2 + j exp(2)((C(j, 2) − j + 1)q 3 − (j − 1)q 2 + q) 2 + C(j, 2) exp(4)((j − 2)q 3 − q 2 ) 2 + C(j, 3) exp(6)q 6 − exp(2jq) on the diagonals, where q = x/2. By using Mathematica and a power series expansion of order 9 for the term exp(2kq) the above expression can be seen not to exceed 0 for all q ≥ 0 which establishes the local D-optimality.