D-optimal Designs for Multinomial Logistic Models

arXiv:1707.03063v1 [math.ST] 10 Jul 2017

Xianwei Bu1 , Dibyen Majumdar2 and Jie Yang2
1
AbbVie Inc. and 2 University of Illinois at Chicago
July 12, 2017
Abstract
We consider optimal designs for general multinomial logistic models, which cover
baseline-category, cumulative, adjacent-categories, and continuation-ratio logit models, with proportional odds, non-proportional odds, or partial proportional odds
assumption. We derive the corresponding Fisher information matrices in three different forms to facilitate their calculations, the determination of conditions for their
positive definiteness, and the search for optimal designs, respectively. We conclude
that, unlike the designs for binary responses, a feasible design for a multinomial
logistic model may contain less experimental settings than parameters, which is
of practical significance. We also conclude that even for a minimally supported
design, a uniform allocation, which is typically used in practice, is not optimal
in general for a multinomial logistic model. We develop efficient algorithms for
searching D-optimal designs. We show using examples based on real experiments
that the efficiency of an experiment can be significantly improved if our designs are
adopted.

Key words and phrases: Approximate design, exact design, Fisher information matrix,
multinomial response, minimally supported design, lift-one algorithm.

1

Introduction

Design of experiments with categorical responses is becoming increasingly popular in a
rich variety of scientific disciplines (Christensen, 2015). Examples include a wine bitterness study (Randall, 1989), trauma clinical trial (Chuang-Stein and Agresti, 1997),
emergence of house flies (Zocchi and Atkinson, 1999), a polysilicon deposition study (Wu,
2008), a toxicity study (Agresti, 2013), and an odor removal study (Yang et al., 2017).
When the response is binary, generalized linear models have been widely used (McCullagh and Nelder, 1989; Dobson and Barnett, 2008). For optimal designs with generalized
linear models, there is a growing body of literature (Khuri et al., 2006; Atkinson et al.,
2007; Stufken and Yang, 2012). In this case, the minimum number of distinct experimental settings required by a nondegenerate Fisher information matrix is equal to the number
of parameters (Fedorov, 1972; Yang and Mandal, 2015). It is also known that the experimental units should be uniformly allocated when a minimally supported design, that is,
1

a design with the least number of experimental settings, is adopted (Yang and Mandal,
2015; Yang et al., 2016).
For many applications, the responses may have three or more categories, which are
more natural or more informative than binary responses. As examples, the response of
odor removal study (Yang et al., 2017) has three levels, serious odor, medium odor,
almost no odor for manufactured bio-plastics; the outcome of trauma clinical trial
(Chuang-Stein and Agresti, 1997) has five categories known as the Glasgow Outcome
Scale (Jennett and Bond, 1975), death, vegetative state, major disability, minor
disability, and good recovery. Some responses are even hierarchical, such as the
outcome of emergence of house flies (Zocchi and Atkinson, 1999) with categories, died
before opening of the pupae, died before complete emergence, and completely
emerged.
For responses with three or more categories, the models used in the literature are special cases of the multivariate generalized linear model (McCullagh, 1980). According to
the relationship among categories, the categorical responses can be of three types: nominal, ordinal, and hierarchical (Zocchi and Atkinson, 1999). In practice, a multinomial
distribution is typically used to model the responses for a specified experimental setting.
In the statistical literature, four kinds of logit models have been commonly used to link
the categorical probabilities to experimental settings or values of covariates, the baselinecategory logit model for nominal responses (Agresti, 2013; Zocchi and Atkinson, 1999),
the cumulative logit model for ordinal responses (McCullagh, 1980; Christensen, 2015),
the adjacent-categories logit model for ordinal responses (Liu and Agresti, 2005; Agresti,
2013), and the continuation-ratio logit model for hierarchical responses (Agresti, 2013;
Zocchi and Atkinson, 1999). In the literature, three different assumptions on parameter
structures have been proposed for the four logit models, proportional odds (po) assuming
the same parameters for different categories (McCullagh, 1980), non-proportional odds
(npo) allowing parameters to change across categories (Agresti, 2013), and partial proportional odds (ppo) incorporating both po and npo components (Peterson and Harrell,
1990). The four logit models and three odds assumptions generate 12 different models
for multinomial responses; for example, we can have, a cumulative logit model with proportional odds. For practitioners, all the 12 models can be fitted using SAS (Stokes et al.,
2012) or R (Yee, 2015), and AIC or BIC criterion may be used for model selection.
Despite this rich complexity of models, the relevant results in the design literature
for multinomial responses are limited to special classes. Zocchi and Atkinson (1999) constructed a general framework of optimal designs for multinomial logistic models with
non-proportional odds. Perevozskaya et al. (2003) discussed a special class of cumulative
logit models with proportional odds. Yang et al. (2017) obtained results for cumulative
link models, which is an extension of cumulative logit models with proportional odds.
Although the cumulative logit model with proportional odds is the most popular model
in practice for ordinal responses, Agresti (2010) found strong evidence against the assumption of proportional odds for the trauma clinical trial data. For the emergence
of house flies data (Zocchi and Atkinson, 1999) a continuation-ratio logit model with
non-proportional odds fits better. It is also known that for certain applications, partial
proportional odds may have a better performance than po or npo models (Lall et al.,
2002).
2

The goal of this work is a comprehensive study of optimal designs for multinomial
responses in a general framework, which covers each of the 12 models. We derive explicit
representations of the Fisher information matrix, derive conditions for positive definiteness, which is necessary for the estimation of the parameters, and study locally D-optimal
designs, theoretically and numerically. We also study EW D-optimal designs that are
good surrogates of Bayesian D-optimal designs with substantially reduced computational
cost.
Our work shows that the optimal designs for multinomial responses with three or more
categories are surprisingly different from the ones for binary responses, and in fact defies
conventional wisdom of optimal design theory, in two major aspects: (i) the required
minimum number of experimental settings is less than the number of parameters; and,
(ii) even for a minimally supported design, the commonly used uniform allocation is not
optimal, except for regular npo models. Indeed, uniform designs can be quite inefficient.
Theoretically this work reveals significant new features of optimal designs for general
multinomial models and shows that, as we go to more and more complex models, some
conventional beliefs of design theory may not hold. For applications, this work provides
a way to accurately compute efficiency of designs of experiments based on multinomial
models, and shows that widely used designs, like uniform designs, may not be efficient
in these complex models.
Our results confirm the findings obtained by Yang et al. (2017) for cumulative link
models and proportional odds.
The rest of this paper is organized as follows: In Section 2, we formulate the general
multinomial logistic model and its Fisher information matrix; in Section 3, we derive a
necessary and sufficient condition for the Fisher information matrix to be positive definite,
which determines the minimal number of required experimental settings; in Section 4,
we formulate the determinant of the Fisher information matrix as a homogeneous polynomial of allocations of experimental units and simplify its structure for D-optimality;
in Section 5, we develop numerical algorithms for searching D-optimal designs, as well
as some analytical results for minimally supported designs; we conclude with remarks in
Section 6. Following Yang et al. (2017), we focus mainly on D-optimality. Nevertheless,
our results on Fisher information matrix are useful for other criteria as well, such as
A-optimality, E-optimality, etc.

2

Multinomial Logistic Model and Its Fisher Information Matrix

We consider an experiment with d ≥ 1 factors and m ≥ 2 distinct experimental setting
xi = (xi1 , . . . , xid )T , i = 1, . . . , m. For the ith experimental setting, ni ≥ 0 categorical
responses are collected i.i.d. from a discrete distribution with J ≥ 2 categories, with
ni = 0 indicating no experimental unit assigned to this experimental setting. When
ni > 0, the responses associated with the ith experimental setting are summarized into
a multinomial response Yi = (Yi1 , · · · , YiJ )T ∼ Multinomial(ni ; πi1 , · · · , πiJ ), where πij
is the probability that the response falls into the jth category at the ith experimental
setting.
3

Throughout this paper, we assume πij > 0 for all i = 1, . . . , m and j = 1, . . . , J,
which is necessary for a multinomial logistic model.
The general structure of models that are considered is a linear regression of log odds
on two components: one specific to the category and the other common to all categories.
Since ppo (partial proportional odds) models include po (proportional odds) and npo
(nonproportional odds) models as special cases, we write the four logit models in the
literature (baseline-category, cumulative, adjacent-categories, and continuation-ratio) in
terms of ppo structure as follows:
 
πij
log
= hTj (xi )β j + hTc (xi )ζ , baseline
πiJ


πi1 + · · · + πij
log
= hTj (xi )β j + hTc (xi )ζ , cumulative
πi,j+1 + · · · + πiJ


πij
= hTj (xi )β j + hTc (xi )ζ , adjacent
log
πi,j+1


πij
= hTj (xi )β j + hTc (xi )ζ , continuation
log
πi,j+1 + · · · + πiJ
where i = 1, . . . , m, j = 1, . . . , J − 1, hTj (·) = (hj1(·), . . . , hjpj (·)) are known functions to
determine the pj predictors for the jth response category, β j = (βj1 , . . . , βjpj )T consists
of pj unknown parameters for the jth response category, hTc (·) = (h1 (·), . . . , hpc (·)) are
known functions to determine common predictors for all categories, ζ = (ζ1 , . . . , ζpc )T
consists of pc unknown parameters. As special cases of ppo, hTj (xi ) ≡ 1 leads to po
models, and hTc (xi ) ≡ 0 leads to npo models.
Following Glonek and McCullagh (1995) and Zocchi and Atkinson (1999), we rewrite
these four logit models into a unified form
CT log(Lπ i ) = η i = Xi θ,
i = 1, · · · , m
P
where π i = (πi1 , . . . , πiJ )T satisfying Jj=1 πij = 1, η i = (ηi1 , . . . , ηiJ )T ,


IJ−1 −IJ−1 0J−1
T
C =
0TJ−1 0TJ−1
1

(1)

is a J × (2J − 1) constant matrix with Ik as the identity matrix of order k and 0k as a
vector of k zeros, L is a (2J − 1) × J constant matrix taking different forms across the
four logit models (see the Appendix for details), the model matrix

 D
h (xi ) hTrow (xi )
(2)
Xi =
0Tpc
0Tp−pc
J×p
with hD (xi ) = diag{hT1 (xi ), . . . , hTJ−1 (xi )} and hTrow (xi ) = [hc (xi ), . . . , hc (xi )]T , and
the parameter vector θ = (β T1 , · · · , β TJ−1 , ζ T )T consists of p = p1 + · · · + pJ−1 + pc
unknown parameters in total. Note that npo models are degenerate cases of ppo models
after removing hTrow (xi ) columns from Xi and ζ from θ, that is, pc = 0; po models are
degenerate cases of ppo models with hD (xi ) replaced by the identify matrix IJ−1 and thus
4

p1 = · · · = pJ−1 = 1. Note that πi1 + · · · + πiJ = 1 implies that ηiJ = 0 and thus the last
row of Xi is all 0’s. We keep ηiJ and the last row of Xi in (1) for convenience following
Glonek and McCullagh (1995). In summary, model (1) covers all four logit models and
all three odds structures (po, npo, and ppo).
Example 2.1. Yee (2015) considered an artificial ppo model with d = 4 factors, J = 3
response categories, and parameter vector θ = (β11 , β12 , β13 , β21 , β22 , β23 , ζ1 , ζ2 )T . At the
ith experimental setting xi = (xi2 , xi3 , xi4 , xi5 )T ,


1 xi3 xi5 0 0
0 xi2 xi4
Xi = 0 0
0 1 xi3 xi5 xi2 xi4  .
0 0
0 0 0
0
0
0



We will revisit this example in Section 4.

Example 2.2. Zocchi and Atkinson (1999) adopted a npo dose-response model for the
study of emergence of house flies with d = 1 factor, J = 3 response categories, and
parameters θ = (β11 , β12 , β13 , β21 , β22 )T . At a dose level xi ,


1 xi x2i 0 0
Xi = 0 0 0 1 xi 
0 0 0 0 0

which is the same for all four logit models. More details about this study will be provided
in Example 5.2 with a continuation-ratio logit link.


Example 2.3. Perevozskaya et al. (2003) studied a po dose-response model with d = 1
factor and J ≥ 3 response categories. At a dose level xi with J = 3,

 
1 0 xi
β1
Xi θ = 0 1 xi  β2 
0 0 0
ζ1



We will revisit this example in Section 5.2.

Using matrix differentiation formulas (see, for example, Seber (2008, Chapter 17)),
we obtain the Fisher information matrix for model (1) as follows:
Theorem 2.1. Consider the multinomial logistic model (1) with independent observations. The Fisher information matrix
F=

m
X

ni Fi

(3)

i=1

where
Fi = (

∂π i
∂π i T
) diag(π i )−1 T
T
∂θ
∂θ

−1
with ∂π i /∂θ T = (CT D−1
i L) Xi and Di = diag(Lπ i ).

5

(4)

Theorem 2.1 is a special case of Glonek and McCullagh (1995) who built a more
general framework for multiple categorical responses. We provide independent proofs in
the Supplementary Materials, as well as Lemma S.1 for later on use. Our results apply
to more general models than Zocchi and Atkinson (1999) who focused on npo models.
Remark 2.1. Given the experimental settings x1 , . . . , xm , one can calculate F1 , . . . , Fm
using (4). Then Theorem 2.1 provides a convenient way for calculating F as a function
of the allocation (n1 , . . . , nm ). The nonsingularity of the key matrix CT D−1
i L was guaranteed by Glonek and McCullagh (1995, Theorem 1). To facilitate the calculation, we
−1
derive explicit forms of (CT D−1
for all the four logit models in the Supplementary
i L)
Materials.

3

Positive Definiteness of the Fisher Information Matrix

A typical requirement for a design is that its information matrix is positive definite,
as this is a necessary condition for the existence of unbiased estimators of parameters
with finite variance (Stoica and Marzetta, 2001). The minimal number of experimental
settings required to keep the Fisher information matrix nonsingular is a basic question
as it deals with the cost of experimentation. In this section, we investigate when the
Fisher information matrix is nonsingular, or equivalently, positive definite, for general
multinomial logistic models.

3.1

Reformulation of the Fisher information matrix

−1
In order to investigate the nonsingularity of F, we denote (CT D−1
= (ci1 , . . . , ciJ ),
i L)
with cij standing for a J × 1 column vector. Then we can rewrite Fi into a simpler form
as a corollary of Theorem 2.1:

Corollary 3.1. Under the setup of Theorem 2.1, the Fisher information at the ith experimental setting
(5)
Fi = XTi Ui Xi
where





Ui = 


u11 (π i )
..
.

· · · u1,J−1 (π i )
..
..
.
.
uJ−1,1 (π i ) · · · uJ−1,J−1 (π i )
0
···
0

and ust (π i ) = cTis diag(π i )−1 cit , for s, t = 1, . . . , J − 1.


0
.. 
. 

0 
1

(6)

With the aid of Corollary 3.1, we can reformulate the Fisher information matrix F into
a form that facilitates the discussion on positive definiteness of F. Its proof is relegated
to the Supplementary Materials.

6

Theorem 3.1. Consider the multinomial logistic model (1) with independent observations. The Fisher information matrix
F = HUHT
where H is



HD
Hrow



D

, H or



(7)
1D
Hrow



for ppo, npo, and po models, respectively, HD = diag{H1 , . . . , HJ−1}, Hrow = [Hc , . . . , Hc ],
1D = diag{1T , . . . , 1T }, Hj = (hj (x1 ), · · · , hj (xm )), j = 1, . . . , J − 1, Hc = (hc (x1 ),
· · · , hc (xm )), and U = (Ust )s,t=1,...,J−1 with Ust = diag{n1 ust (π 1 ), . . . , nm ust(π m )}.

3.2

Positive definiteness of U

In order to determine the positive definiteness of F, we first investigate the m(J − 1) ×
m(J − 1) matrix U defined in Theorem 3.1, which is symmetric since ust (π i ) = uts (π i )
and thus Ust = Uts .
Theorem 3.2. If ni > 0 for all i = 1, . . . , m, then U is positive definite.
QJ
Q
Q
−1 −2
−1
T
J−1
.
· m
Theorem 3.3. |U| = ( m
i=1 ( j=1 πij ) |C Di L|
i=1 ni )

The proofs of Theorems 3.2 and 3.3 are relegated to the Supplementary Materials.
Note that Theorems 3.2 is not a corollary of 3.3 since nonsingularity itself does not
mean positive definiteness. Theorem 3.3 implies that U is singular if ni = 0 for some
i = 1, . . . , m. Note that F can still be positive definite even if U is singular, as long as H
is of full row rank.
units (n1 , . . . , nm )
PmIn general, given an allocation of n experimental
∗
with ni ≥ 0 and i=1 ni = n, if we denote k = #{i : ni > 0} and Ust = diag{ni ust (π i ) :
ni > 0}, then U∗ = (U∗st )s,t=1,...,J−1 is a k(J − 1) × k(J − 1) matrix. After removing all
columns of H associated with ni = 0, we denote the leftover as H∗ , which is a p×k(J −1)
matrix. It can be verified that
Lemma 3.1. HUHT = H∗ U∗ (H∗ )T .
Q
Q
Q
−2
Lemma 3.2. |U∗ | = ( i:ni >0 ni )J−1 · i:ni >0 ( Jj=1 πij )−1 |CT D−1
i L| .

Since U∗ is simply U if all ni > 0, as a corollary of Theorem 3.2, we have

Corollary 3.2. U∗ is positive definite.
As a direct conclusion of Lemmas 3.1, 3.2 and Theorem 3.2, we derive a necessary
and sufficient condition for F to be positive definite:
Theorem 3.4. The Fisher information matrix F is positive definite if and only if H∗ is
of full row rank, where H∗ is the leftover of H after removing all columns associated with
ni = 0. Furthermore, if ni > 0 for all i = 1, . . . , m, then F is positive definite if and only
if H is of full row rank.

7

3.3

Row rank of H matrix

According to Theorem 3.4, the positive definiteness of the Fisher information matrix F
depends on the row rank of H or H∗ . To simplify the notations, we assume ni > 0, i =
1, . . . , m throughout this section. In this case, H = H∗ and U = U∗ . We also assume
that
m ≥ pj , j = 1, . . . , J − 1 and m ≥ pc if applicable
(8)
since H is of full row rank only if rank(Hj ) = pj , j = 1, . . . , J − 1 and rank(Hc ) = pc if
applicable.
Since H takes different forms for ppo, npo, and po models, we investigate its row rank
case by case.
Theorem 3.5. Consider the p × m(J − 1) matrix H in Theorem 3.1.
(1) For npo models, rank(H) = rank(H1 ) + · · · + rank(HJ−1 ).

(2) For po models, rank(H) = rank( 1, HTc ) + J − 2, where 1 is a vector of all 1’s.

(3) For ppo models, rank(H) = rank(H1 )+· · ·+rank(HJ−1 )+rank(Hc )−dim[M(HTc )
T
T
T
∩ (∩J−1
j=1 M(Hj ))], where M(Hc ) stands for the column space of Hc or the row
space of Hc .

The proof of Theorem 3.5 is relegated to the Supplementary Materials. In order to apT
ply it to ppo models, we need an efficient way
dim[M(HTc ) ∩ (∩J−1
j=1 M(Hj ))].
T to calculate
We provide a formula for calculating dim( j M(HTj )) for general matrices, Theorem A.1
in the Appendix, and relegated its proof to the Supplementary Materials.
As a direct conclusion of Theorem 3.5, we have
T
Corollary 3.3. For ppo models, |F | > 0 only if m > pc +pH , where pH = dim(∩J−1
j=1 M(Hj )).

Combining Theorems 3.4, 3.5, and Corollary 3.3, we obtain the major theorem for
positive definiteness of the Fisher information matrix F:
Theorem 3.6. Consider the multinomial logistic model (1) with m distinct experimental
settings xi and the corresponding numbers of observations ni , i = 1, . . . , m. The Fisher
information matrix F is positive definite if and only if
(1) For non-proportional odds (npo) models, m ≥ max{p1 , . . . , pJ−1} and xi ’s keep Hj
of full row rank pj , j = 1, . . . , J − 1.
(2) For proportional odds (po) models, m ≥ pc + 1 and the extended matrix (1, HTc ) is
of full rank pc + 1.
(3) For partial proportional odds (ppo) models, m ≥ max{p1 , . . . , pJ−1, pc + pH } and
xi ’s keep Hj of full row rank
 pj , j = 1, . . . , J − 1; Hc of full row rank pc ; as well
T
as M(HTc ) ∩ ∩J−1
M(H
)
= {0}.
j=1
j
(4) For a special case of ppo models with H1 = · · · = HJ−1, m ≥ pc + p1 and the
extended matrix (HT1 , HTc ) is of full rank pc + p1 .
8

Actually Case (3) in Theorem 3.6 covers all of the three odds structures. For npo
models, pc = 0, pH ≤ min{p1 , . . . , pJ−1 }, then Case (3) leads to Case (1); for po models,
pH = p1 = · · · = pJ−1 = 1, Case (3) is the same as Case (2); and Case (4) is a special of
Case (3).
Theorem 3.6 implies that the number m of distinct experimental settings for a ppo
model could be as low as max{p1 , . . . , pJ−1 , pc + pH }, which is strictly less than the
number of parameters p = p1 + · · · + pJ−1 + pc . The following artificial example serves
as an illustration.
Example 3.1. Consider an experiment with four factors (d = 4), three response categories (J = 3), and four distinct experimental settings (m = 4). Then the experimental
settings are xi = (xi1 , xi2 , xi3 , xi4 )T , i = 1, 2, 3, 4. Consider a multinomial logistic model
with ppo such that






x14
1 x11
1 x11 x12 x13
 x24 
 1 x21 
 1 x21 x22 x23 
T
T





HT1 = 
 1 x31 x32 x33  , H2 =  1 x31  , Hc =  x34 
x44
1 x41
1 x41 x42 x43
That is, p1 = 4, p2 = 2, pc = 1, pH = 2,
p1 + p2 + pc = 7 parameters. In this case,

1
1
1
 x11 x21 x31

 x12 x22 x32

H=
 x13 x23 x33
 0
0
0

 0
0
0
x14 x24 x34

max{p1 , p2 , pc + pH } = 4, and there are p =
1
0
0
0
0
x41 0
0
0
0
x42 0
0
0
0
x43 0
0
0
0
0
1
1
1
1
0 x11 x21 x31 x41
x44 x14 x24 x34 x44












is 7 × 8. It can be verified that rank(H) = 7 for general x1 , x2 , x3 , x4 according to
Theorem 3.5. That is, the minimal number of experimental settings in this case is
m = max{p1 , . . . , pJ−1 , pc + pH } = 4.


4

Determinant of the Fisher Information Matrix

The D-criterion maximizes the determinant of the Fisher information matrix F. In order
to characterize theoretical properties of D-optimal designs, we need to investigate the
structure of |F|. In this section, we reformulate F into the form of GT WG with a
diagonal matrix W so that |F| becomes a homogeneous polynomial of allocations.
Theorem 4.1. Consider the multinomial logistic model (1) with independent observations. The Fisher information matrix
F = nGT WG

9

(9)

where n is the total number of observations with ni of them assigned to the ith experimental setting xi , W = diag{w1 diag(π 1 )−1 , . . . , wm diag(π m )−1 } is an mJ ×mJ diagonal
matrix with proportions wi = ni /n, and G is an mJ × p matrix which takes the forms of


PJ−1
T
c
·
h
(x
)
c11 hT1 (x1 ) · · · c1,J−1 hTJ−1 (x1 )
1j
1
c
Pj=1


J−1
c2j · hTc (x2 ) 
 c21 hT1 (x2 ) · · · c2,J−1 hTJ−1 (x2 )
j=1
Gppo = 

···
···
···
···


P
J−1
T
c
·
h
(x
)
cm1 hT1 (xm ) · · · cm,J−1 hTJ−1 (xm )
m
c
j=1 mj


T
T
c11 h1 (x1 ) · · · c1,J−1 hJ−1 (x1 )
 c21 hT1 (x2 ) · · · c2,J−1 hTJ−1 (x2 ) 

Gnpo = 


···
···
···
T
T
cm1 h1 (xm ) · · · cm,J−1 hJ−1 (xm )


PJ−1
T
c11 · · · c1,J−1
c
·
h
(x
)
1
c
j=1 1j
PJ−1


T
c
·
h
(x
 c21 · · · c2,J−1
2) 
c
j=1 2j
Gpo = 

···
···

 ··· ···
PJ−1
T
cm1 · · · cm,J−1
j=1 cmj · hc (xm )
for ppo, npo, and po models, respectively.

In order to find D-optimal designs, we need to maximize |GT WG|. Since W is
diagonal, we obtain the following theorem as a direct conclusion of Theorem 1.1.2 of
Fedorov (1972) or Lemma 3.1 of Yang and Mandal (2015):
Theorem 4.2. The determinant of the Fisher information matrix
X
αm
cα1 ,...,αm · w1α1 · · · wm
|F| = np |GT WG| = np
α1 ≥0,...,αm ≥0 :

Pm

i=1

(10)

αi =p

with cα1 ,...,αm =
X

(i1 ,...,ip )∈Λ(α1 ,...,αm )

|G[i1 , . . . , ip ]|2

Y

Y

k:αk >0 l:(k−1)J<il 6kJ

−1
πk,i
≥0
l −(k−1)J

(11)

where α1 , . . . , αm are nonnegative integers, Λ(α1 , . . . , αm ) = {(i1 , . . . , ip ) | 1 ≤ i1 <
· · · < ip ≤ mJ; #{l : (k − 1)J < il 6 kJ} = αk , k = 1, . . . , m}, and G[i1 , . . . , ip ] is the
submatrix consisting of the i1 th, . . . , ip th rows of G.
According to Theorem 4.2, the determinant of the Fisher information matrix is an
order-p homogeneous polynomial of w1 , . . . , wm . Another important conclusion is that
the coefficient cα1 ,...,αm must be nonnegative. Actually, the structure of the determinant
can be significantly simplified due to the following results on the coefficient cα1 ,...,αm :
Lemma 4.1. If max1≤i≤m αi ≥ J, then |G[i1 , . . . , ip ]| = 0 for any (i1 , . . . , ip ) ∈ Λ(α1 ,
. . . , αm ). Therefore, cα1 ,...,αm = 0 in this case.
Theorem 4.3. The coefficient cα1 ,...,α
Pm as defined in (11) is nonzero only if the restricted
Fisher information matrix Fres = i:αi >0 Fi is positive definite, where Fi is defined as
in (4).
10

The proofs for Lemma 4.1 and Theorem 4.3 are relegated to the Supplementary
Materials. Combining Theorems 3.4 and 4.3, Theorems 3.6 and 4.3, respectively, we
obtain the following corollaries:
Corollary 4.1. The coefficient cα1 ,...,αm is nonzero only if Hα1 ,...,αm is of full row rank p,
where Hα1 ,...,αm is the submatrix of H after removing all columns associated with xi for
which αi = 0.
Corollary 4.2. The coefficient cα1 ,...,αm = 0 if

max{p1 , . . . , pJ−1 }



pc + 1
kmin =
max{p1 , . . . , pJ−1 , pc + pH }



pc + p1

#{i | αi > 0} ≤ kmin − 1, where
for
for
for
for

npo models;
po models;
ppo models;
ppo with H1 = · · · = HJ−1 .

For typical applications kmin ≥ 2 is the minimal number of experimental settings to
keep the positive definiteness of F, which will be revisited in Section 5.5. As special cases
of ppo models, npo models imply pc = 0 and pH ≤ min{p1 , . . . , pJ−1 }, po models imply
p1 = · · · = pJ−1 = pH = 1, and H1 = · · · = HJ−1 implies p1 = · · · = pJ−1 = pH . That is,
kmin ’s are consistent across different odds models.
In general, the determinant of a ppo model is more complicated than the determinant
of a npo or po model. We provide an example in the Supplementary Materials to illustrate
that cα1 ,...,αm could be nonzero for ppo models with #{i | αi > 0} = pc +pH . Nevertheless,
the following example shows that Lemma 4.1 and Corollary 4.2 may simplify the structure
of |F| significantly.

Example 2.1 (Continued) : In this example, the number of factors is d = 4, and the
experimental settings are xi = (xi2 , xi3 , xi4 , xi5 )T , i = 1, . . . , m. Since p1 = p2 = 3, pc = 2,
and the number of parameters p = p1 + p2 + pc = 8, the minimal number of experimental
settings is m = p1 + pc = 5 according to Theorem 3.6. We consider the simplest case
m = 5. That is,




1 x13 x15
x12 x14

..  , HT =  ..
.. 
HT1 = HT2 =  ... ...
 .
. 
. 
c
1 x53 x55
x52 x54

According to (10) of Theorem 4.2, |F| is an order-8 homogeneous polynomial allocating
to 5 experimental settings which may contain up to (8 + 5 − 1)!/(8!(5 − 1)!) = 465 terms.
However, Lemma 4.1 implies cα1 ,...,α5 6= 0 only if αi ∈ {0, 1, 2}. On the other hand,
Corollary 4.2 says cα1 ,...,α5 6= 0 only if #{i | αi > 0} ≥ p1 + pc = 5, that is, αi > 0 for
each i. Therefore |F| contains only 5!/(3!2!) = 10 nonzero terms with αi ∈ {1, 2}:
|F| =

5
Y
i=1

wi ·

X

ei1 ,i2 ,i3 wi1 wi2 wi3

1≤i1 <i2 <i3 ≤5

for some coefficients ei1 ,i2 ,i3 . Actually, ei1 ,i2 ,i3 can be determined by cα1 ,...,αm in (11) of
Theorem 4.2. For example, e1,2,3 = c2,2,2,1,1 .

11

Remark 4.1. The reformulation of F in Theorem 4.1 enables us to conclude that |F| is an
order-p homogeneous polynomial of allocations in Theorem 4.2. It establishes the foundation for an efficient numerical algorithm for searching D-optimal designs (Section 5).
On the other hand, the simplification of |F| based on Lemma 4.1 and Corollary 4.2 allows us to obtain D-optimal designs analytically for minimally supported designs, which
is critical for investigating their theoretical properties (Section 5.5).
With the aid of Theorem 4.1, Theorem 4.2 is derived in a more straightforward
way. It is broader and with more detailed structures than Theorem 2.2 in Yang et al.
(2017), which focused on po models only. Lemma 4.1 and Corollary 4.2 also extend the
corresponding results in Yang et al. (2017, Lemmas S.4 and S.5 in the Supplementary
Materials).


5

D-optimal Designs

A design here is an allocation (n1 , . . . , nm ) of n experimental units, or proportions
(w1 , . . . , wm ) on experimental settings xi = (xi1 , . . . , xid )T , i = 1, . . . , m. We first
consider the design problem with pre-determined experimental settings. In this case, a
locally D-optimal design is an allocation which maximizes the determinant of the Fisher
information matrix |F|, or equivalently |GT WG| in Theorem 4.2, with specified parameter values. In the last part of this section, we consider the problem of optimizing over
experimental settings as well, using a grid-point search algorithm.

5.1

Design space for multinomial logistic models

A multinomial logistic model connects the categorical probabilities πij ’s and the parameter vector θ via a logit link and the model matrices Xi ’s in (2) consisting of hTj (xi ) and
hTc (xi ) if applicable. Due to the requirements of 0 < πij < 1, j = 1, . . . , J for multinomial
logistic models, from a design point of view, we need to determine the collection of all
feasible experiment settings, known as the design space, X = {x = (x1 , . . . , xd )T | 0 <
πj < 1, j = 1, . . . , J exist uniquely}.
Theorem 5.1. Let aj = hTj (x)β j + hTc (x)ζ, j = 1, . . . , J − 1. The design space X =
{x = (x1 , . . . , xd )T | aj ∈ (−∞, ∞), j = 1, . . . , J − 1} for baseline-category, adjacentcategories, and continuation-ratio logit models; and X = {x = (x1 , . . . , xd )T | −∞ <
a1 < a2 < · · · < aJ−1 < ∞} for cumulative logit models.
The proof of Theorem 5.1 is relegated to the Supplementary Materials. It essentially
puts no restriction on the design space for models other than cumulative logit models,
regardless of odds assumption.
Example 5.1. Trauma clinical trial Chuang-Stein and Agresti (1997) studied a dataset
of trauma patients which have five ordered response categories, death, vegetative
state, major disability, minor disability, and good recovery, describing their
clinical outcomes. These five categories are often called the Glasgow Outcome Scale
(GOS) in the literature of critical care (Jennett and Bond, 1975). In the original clinical
trial, 802 patients were roughly equally assigned to the four treatment groups labeled as
12

placebo, low dose, medium dose, and high dose as reported in Chuang-Stein and Agresti
(1997, Table I).
Although cumulative logit models with proportional odds are commonly used in practice for ordered categorical responses, Agresti (2010) found some strong evidence against
the assumption of proportional odds for the trauma data. We fit the data with different
logit models, as well as different odds assumptions, and confirm that the cumulative link
model with non-proportional odds fits the trauma data the best in terms of AIC and
BIC (see Table 4 in the Supplementary Materials). That is,


πi1 + · · · + πij
log
= βj1 + βj2 xi , i = 1, 2, 3, 4; j = 1, 2, 3, 4
(12)
πi,j+1 + · · · + πi5
There are eight parameters with fitted values β̂ = (β̂11 , β̂21 , β̂31 , β̂41 , β̂12 , β̂22 , β̂32 , β̂42 )T =
(−0.865, −0.094, 0.706, 1.909, −0.113, −0.269, −0.182, −0.119)T . If we treat the fitted
parameter values as the true values, the design space is X = {x ≥ 0 | β11 + β12 x <
β21 + β22 x < β31 + β32 x < β41 + β42 x} = {x ≥ 0 | −9.195 < x < 4.942} = [0, 4.942). It is
not a surprise that the four levels {1, 2, 3, 4} in the original dataset are included in the
design space.


5.2

Reparametrization and D-optimaltiy

In the literature, model parameters other than our θ = (β T1 , · · · , βTJ−1 , ζ)T for linear
predictors ηij = hTj (xi )β j + hTc (xi )ζ were also used, for example, in Perevozskaya et al.
(2003). We use it as an example to show that D-optimal designs won’t be affected by
the choice of parameters. More general discussion is relegated to the Appendix.
Example 2.3 (continued): Perevozskaya et al. (2003) considered the po model:
log

x − αj′
γj (x)
=
1 − γj (x)
β′

j = 2, . . . , J

(13)

where γj (x) = P (Y ≥ j|x). Let us reparametrize this model as
log

γj (x)
= αj + βx
1 − γj (x)

j = 2, . . . , J

(14)

Let θ = (α2 , α3 , β)T be the parameters in (13), and and ϑ = (α2′ , α3′ , β ′ )T be the parameters in (14). Then β = 1/β ′, α2 = −α2′ /β ′, α3 = −α3′ /β ′, and the Jacobian matrix


α′2
0
− β1′
β ′2
′ 

J =  0 − β1′ βα′23 
0
0 − β1′2
Based on Theorem 2.1, the Fisher information Ii (θ) at xi is

 πi1 π2 πi1,2
πi1 πi1,2 πi2,3 πi3
i2,3
−
π
π
π
x
i1
i1,2
i2,3
i
πi2
πi2

 πi1 πi1,2
2 π
πi1,2
i2,3 πi3
πi2,3 πi3

−
π
π
π
x
i3
i1,2
i2,3
i
πi2
πi2
2
2
2
2
πi1 πi1,2 πi2,3 xi πi3 πi1,2 πi2,3 xi (πi1 πi2,3 + πi2 (πi1 − πi3 ) + πi1,2 πi3 )xi
13

where πij,k = πij + πik . It can be verified that Ii (ϑ) = JT Ii (θ)J equals to the corresponding one given by Perevozskaya et al. (2003). For any given design ξ = {(x
P i , wi ), i =
1, . . . , m} with proportions wi ∈ [0, 1], the Fisher information matrix Iξ (ϑ) = m
i=1 wi Ii (ϑ) =
T
2
J Iξ (θ)J. Then |Iξ (ϑ)| = |J| · |Iξ (θ)| and the D-optimal design maximizing |Iξ (θ)| also
maximizes |Iξ (ϑ)|. That is, the D-optimal designs for Models (13) and (14) are the same.


5.3

D-optimal approximate designs

Given distinct experimental settings xi ∈ X , i = 1, . . . , m, we are looking for D-optimal
approximate allocations w = (w1 , . . . , wm )T that maximizes |GT WG| defined in (9).
Note that xi ∈ X implies 0 < πij < 1, j = 1, . . . , J.
The collection
all feasible approximate allocations S = {(w1 , . . . , wm )T ∈ Rm | wi ≥
Pof
m
0, i = 1, . . . , m; i=1 w1 = 1} is a bounded closed convex set. The objective function
f (w) = |GT WG|

is an order-p homogeneous polynomial according to Theorem 4.2. Therefore, a D-optimal
approximate design that maximizes f (w) must exist. For typical applications, we need
designs coming from S+ = {w ∈ S | f (w) > 0}. Due to Theorem 2.1 and log-concavity
of determinant on positive semi-definite matrices, we know f (w) is log-concave (Silvey,
1980; Yang et al., 2017) and S+ is convex. A useful result as a corollary of Theorem 3.4
is as follows:
Corollary 5.1. S+ is nonempty if and only if f (wu ) > 0, where wu = (1/m, . . . , 1/m)T
is the uniform allocation. In this case, f (w) > 0 for any w = (w1 , . . . , wm )T such that
0 < wi < 1, i = 1, . . . , m.
In order to avoid trivial cases, we assume f (wu ) > 0 from now on. Following
Yang et al. (2017, Section 3) (see also Yang and Mandal (2015) and Yang et al. (2016)),
we define fi (z) = f (w1 (1 − z)/(1 − wi ), . . . , wi−1 (1 − z)/(1 − wi ), z, wi+1 (1 − z)/(1 −
wi ), . . . , wm (1 − z)/(1 − wi )) with 0 ≤ z ≤ 1 and w = (w1 , . . . , wm )T ∈ S+ . Parallel to Theorem 3.2.6 in Yang et al. (2017), we obtain the result as follows according to
Theorem 4.2:
Theorem 5.2. Given an approximate allocation w = (w1 , . . . , wm )T ∈ S+ and an i ∈
{1, . . . , m}, for 0 < z < 1,
fi (z) = (1 − z)
fi′ (z) = (1 − z)p−J

J−1
X
j=1

p−J+1

J−1
X

j=0
j−1

bj (j − pz)z

bj z j (1 − z)J−1−j

(15)

(1 − z)J−j−1 − pb0 (1 − z)p−1

(16)

t−1
where b0 = fi (0), (bJ−1 , . . . , b1 )T = B−1
)st is a (J − 1) × (J − 1)
J−1 c, BJ−1 = (s
T
p J−1−p
constant matrix, and c = (c1 , . . . , cJ−1) with cj = (j +1) j
fi (1/(j +1))−j J−1 fi (0),
j = 1, . . . , J − 1.

14

Theorem 5.2 shows that fi (z) is an order-p polynomial of z. Since fi (1) = 0, the
solution to maximization of fi (z), 0 ≤ z ≤ 1 could only occur at z = 0 or 0 < z < 1 such
that fi′ (z) = 0, that is,
J−1
X
j=1

jbj z j−1 (1 − z)J−j−1 = p

J−1
X
j=0

bj z j (1 − z)J−j−1 , 0 < z < 1

(17)

which is an order-(J − 1) polynomial of z. For J ≤ 5, (17) is a polynomial equation of
order 4 or less which can be solved analytically. For J ≥ 6, a quasi-Newton algorithm
could be applied for searching numerical solutions.
In order to find D-optimal designs numerically, we use a lift-one algorithm (see the
Supplementary Materials), which is essentially the same as the one in Yang et al. (2017)
for cumulative link models. The lift-one algorithm is of general-equivalence-theorem
type (Yang et al., 2016, 2017). Its convergence to a global maximum is guaranteed
(Yang and Mandal, 2015).
Example 5.2. Emergence of house flies Zocchi and Atkinson (1999) discussed an
experiment on emergence of house flies originally reported by Itepan (1995). In this
experiment, ni = 500 pupae were exposed to each of seven doses of radiation in Gy,
xi = 80, 100, 120, 140, 160, 180, 200, i = 1, . . . , 7, respectively. After a period of time,
for each set of 500 pupae, the summarized responses include the number yi1 of flies that
died before the opening of the pupae (unopened pupae), the number yi2 of flies out of
opened pupae but died before complete emergence, and the number yi3 = ni − yi1 − yi2
of flies out of opened pupae and completely emerged. Given ni , the summary responses
(yi1 , yi2 , yi3) follow a multinomial distribution and have a clearly nested or hierarchical
structure (see Table 1 of Zocchi and Atkinson (1999) for the experimental data). We
confirm that the continuation-ratio logit model with npo fits the data the best in terms
of AIC and BIC (see Table 5 in the Supplementary Materials), which was adopted by
Zocchi and Atkinson (1999) as follows


 
πi1
πi2
2
log
= β11 + β12 xi + β13 xi , log
= β21 + β22 xi
πi2 + πi3
πi3
The model has 5 parameters with fitted values β̂ = (β̂11 , β̂12 , β̂13 , β̂21 , β̂22 )T = (−1.935,
−0.02642, 0.0003174, −9.159, 0.06386)T .
Table 1: Exact and Approximate Allocations for House Flies Experiment
Dose of radiation (Gy)
Original allocation
D-optimal exact
Original proportion
D-optimal approximate
Bayesian D-optimal
EW D-optimal

80
500
1091
.1429
.3116
.3159
.3120

100
500
0
.1429
0
.0000
0

120
500
1021
.1429
.2917
.2692
.2911

140
500
374
.1429
.1071
.1160
.1087

160
500
1014
.1429
.2896
.2990
.2882

180
500
0
.1429
0
.0000
0

200
500
0
.1429
0
.0000
0

Considering a followup experiment treating the fitted parameter values as the true
values, we obtain the D-optimal approximate design (Table 1) using the lift-one algorithm. The efficiency of the original uniform allocation is (|Foriginal |/|FD−opt |)1/5 =
15

(585317/1480378)1/5 = 83.1%. The D-optimal approximate design only requires 4 rounds
(doses) of experiments. The D-optimal exact design, Bayesian and EW D-optimal designs
listed in Table 1 will be described in Sections 5.4 and 5.6, respectively.


5.4

D-optimal exact designs

In practice, a design problem often targets an integer-valued allocation of n experimental
units on m distinct experimental settings
Pm xi , i = 1, . . . , m. An integer-valued
Pm allocation
T
n = (n1 , . . . , nm ) that maximizes | i=1 ni Fi | as defined in (3) with i=1 ni = n is
known as a D-optimal exact design.
Although different rounding algorithms have been proposed in the design literature
to obtain an exact allocation from an optimal approximate allocation, an algorithm
searching for optimal exact designs directly is still needed (see Yang et al. (2017, Section 4) and reference therein).
P For simplicity, we also denote the objective function
as f (n) = f (n1 , . . . , nm ) = | m
i=1 ni Fi | if there is no ambiguity. Following Yang et al.
(2016) and Yang et al. (2017), we define fij (z) = f (n1 , . . . , ni−1 , z, ni+1 , . . . , nj−1, ni +nj −
z, nj+1 , . . . , nm ) with z = 0, 1, . . . , ni + nj given 1 ≤ i < j ≤ m and n = (n1 , . . . , nm )T .
Similar to Theorem 4.9 in Yang et al. (2017), we obtain the following result from
Theorem 4.2, Lemma 4.1 and Corollary 4.2:
Theorem 5.3. Suppose n = (n1 , . . . , nm )T satisfies f (n) > 0 and ni + nj ≥ q for given
1 ≤ i < j ≤ m, where q = min{2J − 2, p − kmin + 2, p}. Then
fij (z) =

q
X

gs z s ,

z = 0, 1, . . . , ni + nj

(18)

s=0

T
where g0 = fij (0), and g1 , . . . , gq can be obtained using (g1 , . . . , gq )T = B−1
q (d1 , . . . , dq )
t−1
with Bq = (s )st as a q × q constant matrix and ds = (fij (s) − fij (0))/s.

The proof is relegated to the Supplementary Materials. Given n > 0, we assume that
there exists an exact allocation n = (n1 , . . . , nm )T such that f (n) > 0 (otherwise, the
maximization problem is trivial). In this case, if n ≥ m, one may choose any n such that
ni > 0, i = 1, . . . , m. An exchange algorithm (see the Supplementary Materials), which
is essentially the same as the one in Yang et al. (2017, S1.5), is applied for obtaining an
D-optimal exact design.
Example 5.1: Trauma clinical trial (continued) In this experiment, the original
allocation assigned 210, 190, 207, 195 patients to the four treatment groups, Placebo,
Low dose, Medium dose, and High dose, respectively, which is roughly uniform. Using
the exchange algorithm, we find the D-optimal exact design (401, 0, 0, 401), which is
minimally supported and uniform on its supporting points. Compared with the D-optimal
design, the efficiency of the original allocation is only 74.7%.

Example 5.2: Emergence of house flies (continued) In this experiment, the original
allocation of experimental units is 500 for each of the seven doses of radiation. Using
the exchange algorithm, we obtain the D-optimal exact allocation (Table 1), which is
similar to the D-optimal approximate design. Actually, the proportions of the D-optimal
16

exact design (1091, 1021, 374, 1014)/3500 = (0.3117, 0.2917, 0.1069, 0.2897) on the four
support points are roughly the same as the ones in the D-optimal approximate design.
The efficiency of the original allocation is also 83.1% compared with the D-optimal exact
allocation.


5.5

Minimally supported designs

A fundamental question in design theory is, what is the least number of design points
we need to keep F positive definite? In other words, what is the number of distinct
experimental settings for a minimally supported design? According to Theorem 3.4,
two conditions are needed: (1) πij > 0, j = 1, . . . , J for all design point xi , that is,
xi ∈ X , i = 1, . . . , m (see Section 5.1); (2) H is of full row rank p (we assume ni > 0 for
i = 1, . . . , m in this section since we are considering minimal m). Theorem 3.6 provides
a lower bound, denoted as kmin in Corollary 4.2. For typical models proposed in the
literature, kmin is the minimal number of experimental settings since the conditions in
Theorem 3.6 are satisfied. Answers to general cases rely on the specific forms of the
predictor functions hj (and hc if applicable). In general we have kmin < p unless J = 2.
Another question is whether a uniform allocation wu = (1/m, . . . , 1/m)T is D-optimal
given that m is the minimally required number of experimental settings. The answer is
known to be “Yes” for J = 2.
Theorem 5.4. Consider Multinomial logit model (1) with only two response categories
(J = 2). In this case, the minimum number of support points is m = p. The objective
function f (w) ∝ w1 · · · wm and the D-optimal allocation for a minimally supported design
is w = (1/m, . . . , 1/m)T .
It can be verified that with J = 2 all of the four logit models are equivalent to the
usual logistic model for binary response. In this case, po, npo, or ppo are essentially
the same. Theorem 5.4 confirms the corresponding results for binary responses in the
literature (see, for example, Yang and Mandal (2015)). We provide an independent proof
in the Supplementary Materials.
Besides the cases with J = 2, for certain npo models with J ≥ 3, uniform allocations
could still be D-optimal for minimally supported designs. Actually, as a direct conclusion
of Theorems 3.1 and 3.3, we have the following corollary, which confirms the D-optimal
design discovered for the trauma clinical trial (Example 5.1).
Corollary 5.2. Consider multinomial logit models (1) with npo assumption. Suppose
p1 = · · · = pJ−1 and there exist p1 distinct experimental settings such that rank(H1 ) =
· · · rank(HJ−1 ) = p1 . Then the minimal number of experimental settings is m = p1 and
the uniform allocation is D-optimal for a minimally supported design.
According to Corollary 5.2, for “regular” npo models (that is, p1 = · · · = pJ−1 ),
uniform allocations are still D-optimal if restricted on a minimally supported design
even with J ≥ 3. Since Yang et al. (2017) showed that, for cumulative link models with
proportional odds, a D-optimal allocation on a minimally supported design is not uniform
in general. Corollary 5.2 implies that the optimality of uniform allocations depends on
the odds structure.
17

Nevertheless, if the condition p1 = · · · = pJ−1 in Corollary 5.2 is violated, the following
lemma and example represent that even for npo models, the uniform allocations are not
D-optimal in general.
Lemma 5.1. Given 0 < c1 ≤ c2 ≤ c3 , we consider the maximization problem f (w1 , w2 , w3 )
= w1 w2 w3 ( c1 w2 w3 + c2 w1 w3 + c3 w1 w2 ) with respect to 0 ≤ wi ≤ 1 and w1 + w2 + w3 = 1.
Then the solution is w1 = w2 = w3 = 1/3 if and only if c1 = c2 = c3 .
The proof of Lemma 5.1 is relegated to the Supplementary Materials, where analytical
solutions are provided for (w1 , w2 , w3 ) for all possible values of c1 , c2 and c3 .
Example 2.2 (continued): Consider the npo model adopted by Zocchi and Atkinson
(1999) with h1 (xi ) = (1, xi , x2i )T , h2 (xi ) = (1, xi )T , J = 3, p1 = 3, p2 = 2, and p = 5.
According to Corollary 4.2, the minimum number of support points is m = max{p1 , p2 } =
3, which is feasible. The objective function f (w) is an order-5 polynomial with terms
cα1 ,α2 ,α3 w1α1 w2α2 w3α3 . Lemma 4.1 implies that αi ∈ {0, 1, 2}, i = 1, 2, 3 in order to keep
cα1 ,α2 ,α3 6= 0. Combined with Corollary 4.2, we further know αi ∈ {1, 2}, i = 1, 2, 3.
According to Theorem 4.2, the objective function is
f (w1 , w2 , w3 ) = w1 w2 w3 (c122 w2 w3 + c212 w1 w3 + c221 w1 w2 )

(19)

for all the four logit models. Rewriting (c122 , c212 , c221 ) = C · (c1 , c2 , c3 ), it can be verified
that for the continuation-ratio logit model adopted by Zocchi and AtkinsonQ
(1999)
Q3for the
3
2
2
2
house flies experiment (Example 5.2), C = (x1 − x2 ) (x1 − x3 ) (x2 − x3 )
j=1 πij ,
i=1
−1
−1
−1
−1
−1
−1
2
2
2
);
c1 = (x2 − x3 ) (π12 + π13 ), c2 = (x1 − x3 ) (π22 + π23 ), c3 = (x1 − x2 ) (π32 + π33
2
for a cumulative
logit
model
(see,
for
example,
Example
5.1),
C
=
(x
−
x
)
(x
−
1
2
1
Q
−1
−1
x3 )2 (x2 − x3 )2 3i=1 πi1 πi2
πi3 (πi1 + πi2 )2 (πi2 + πi3 )2 , c1 = (x2 − x3 )2 π13
(π11 + π12 )−1 ,
−1
−1
c2 = (x1 − x3 )2 π23
(π21 + π22 )−1 , and c3 = (x1 − x2 )2 π33
(π31 + π32 )−1 . According to
Lemma 5.1, w1 = w2 = w3 = 1/3 is D-optimal if and only if c1 = c2 = c3 , which is
in general not true for both continuation-ratio and cumulative logit models with nonproportional odds.

Yang et al. (2017) showed analytically that uniform allocations are not D-optimal
even for minimally supported designs with po models and cumulative link. We use
Example S.3 in the Supplementary Materials to show that the conclusion is the same
for po models with baseline-category, adjacent-categories, and continuation-ratio logit
links as well. Actually, the objective functions take the same form for the four different
logit models. Since po models are special cases of ppo models, we conclude that uniform
allocations are not D-optimal in general for ppo models as well.

5.6

EW D-optimal designs

Both D-optimal approximate designs and exact designs discussed previously are “locally” D-optimal designs since the values of parameters need to be specified in advance.
Bayesian D-optimality (Chaloner and Verdinelli, 1995), that maximizes E(log |F|) with
a given prior distribution on the unknown parameters, provides an alternative approach.
A drawback of Bayesian approach is its computational intensity since its objective function deals with multiple integrals. An alternative solution is the EW D-optimality
18

(Atkinson et al., 2007; Yang et al., 2016, 2017), which maximizes log |E(F)| or |E(F)|
instead. Among all of the criteria proposed by Atkinson et al. (2007, Table 18.1) as surrogates for Bayesian D-optimality, including − log E(|F|−1 ), − log |E(F−1 )|, and log E|F|,
EW D-optimal design requires the minimum computation. Yang et al. (2016) showed
that an EW D-optimal design could be highly efficient in terms of Bayesian criterion
compared with the Bayesian D-optimal design, while the computational time is essentially the same as a locally D-optimal one. Yang et al. (2017) also used EW-criterion for
cumulative link models with proportional odds and confirmed its high-efficiency.
According
2.1 and Corollary 3.1, the Fisher information matrix in our
Pmto Theorem
T
case is F = i=1 ni Xi Ui Xi , where Ui consists of ust (π i ), the only components involving
model parameters. In order to calculate E(F) with respect to a prior on parameters, we
may calculate E(ust(π i )) first and then the results and algorithms developed for locally
optimal designs can be used for EW D-optimal designs directly.
We provide formulas in both the proof of Theorem 5.1 and Section S.5 (in matrix
form) in the Supplementary Materials for calculating πij ’s given Xi ’s and the parameter
values. The following results provide formulas for calculating ust (π i )’s given πij ’s.
Theorem 5.5. Consider the unified multinomial logit model (1). Recall that ust (π i ) =
cTis diag(π i )−1 cit . For i = 1, . . . , m,
(i) ust (π i ) = uts (π i ), s, t = 1, . . . , J;
(ii) usJ (π i ) = 0 for s = 1, . . . , J − 1 and uJJ (π i ) = 1;
(iii) For s = 1, . . . , J − 1,

πis (1 − πis )


 2
−1
−1
+ πi,s+1
)
γis (1 − γis )2 (πis
uss (π i ) =
γ
(1
−
γ
),

is
is


πis (1 − γis )(1 − γi,s−1)−1
(iv) For 1 ≤ s < t ≤ J − 1,

−πis πit




 −γis γit (1 − γis )(1 − γit )πit−1
0
ust (π i ) =


γis (1 − γit )



0

for
for
for
for

baseline-category,
cumulative,
adjacent-categories,
continuation-ratio;

for
for
for
for
for

baseline-category,
cumulative, t − s = 1,
cumulative, t − s > 1,
adjacent-categories,
continuation-ratio;

where γij = πi1 + · · · + πij , j = 1, . . . , J − 1; γi0 ≡ 0 and γiJ ≡ 1.
According to Theorem 5.1, for baseline-category, adjacent-categories, and continuationratio logit models, there is essentially no restriction on values of parameters for a given
design point x = (x1 , . . . , xd )T . The corresponding E(ust (π i )) can be integrated over
the same range of parameter values for different design points. However, for cumulative logit models, given a design point x, the values of parameters should satisfy
a1 < a2 < · · · < aJ−1 with aj = hTj (x)β j + hTc (x)ζ. The domain of integration in
this case is typically not rectangular.
19

Example 5.1: Trauma clinical trial (continued) In this example, the cumulative
logit model (12) with non-proportional odds was adopted. Given the predetermined
set X = {1, 2, 3, 4} consisting of m = 4 design points, the feasible parameter space is
Θ = {θ = (β11 , β21 , β31 , β41 , β12 , β22 , β32 , β42 )T | β11 + β12 x < β21 + β22 x < β31 + β32 x <
β41 + β42 x, for x ∈ X }, which is not rectangular. Zocchi and Atkinson (1999) used a
multivariate normal prior with its variance-covariance matrix estimated from an initial
dataset. For the computations in this example, we bootstrap the 802 observations from
the initial dataset for 1000 times and denote the corresponding fitted parameters as
θ 1 , . . . , θ 1000 . Then an estimate of the Bayesian criterion φ(w) = E(log |F(w, θ)|) for
design w = (w1 , . . . , wm )T is
1000

1000

m

X
1 X
1 X
φ̂(w) =
log |F(w, θj )| =
log |
wi Fi (θ j )|
1000 j=1
1000 j=1
i=1
The Bayesian D-optimalPdesign listed in Table 2 maximizes φ̂(w), while the EW Dm
optimal
P1000one maximizes | i=1 wi Ê(Fi )| with estimated entry-by-entry expectation Ê(Fi )
= j=1 Fi (θ j )/1000. The two designs are essentially the same in this case.

Table 2: D-optimal Designs for Trauma Clinical Trial
Design point x
Treatment group
Original allocation
D-optimal exact
Original proportion
D-optimal approximate
Bayesian D-optimal
EW D-optimal

1
Placebo
210
401
0.2618
0.5
0.4997
0.5

2
Low
190
0
0.2369
0
1.0 × 10−8
0

3
Medium
207
0
0.2581
0
2.5 × 10−8
0

4
High
195
401
0.2431
0.5
0.5003
0.5

Example 5.2: Emergence of house flies (continued) Similar to the trauma example, we bootstrap the original observations for 1000 times and obtain the corresponding
Bayesian and EW D-optimal designs (see Table 1). In this case, EW D-optimal design is
close to the Bayesian one. In order to check the robustness of these designs towards misspecified parameter values, we consider the same set of bootstrapped parameter vectors
θ 1 , . . . , θ 1000 . For j = 1, . . . , 1000, we apply our lift-one algorithm to find the corresponding D-optimal allocation pj assuming θ j as the true parameter values. The efficiencies
of a target design p with respect to pj is defined as (|F(p, θ j )|/|F(pj , θj )|)1/p with p = 5
parameters in this case. The summary statistics of efficiencies in Table 3 show that
Bayesian and EW D-optimal designs are highly robust in terms of parameter misspecification. EW design is slightly better than Bayesian one. Both of them are much better
than the original uniform design performed in the previous study.


5.7

Grid point search for continuous factors

Both the lift-one and exchange algorithms are for searching optimal design with a predetermined finite set of experimental settings. When a factor is continuous, one common
20

Table 3: Efficiencies of Designs for House Flies Experiment
Design
Bayesian D-optimal
EW D-optimal
Uniform design

Min
0.9912
0.9934
0.7762

1st Quartile
0.9985
0.9991
0.8181

Median
0.9989
0.9996
0.8304

3rd Quartile
0.9992
0.9998
0.8445

Max
0.9998
1.0000
0.8861

practice is to partition the continuous region of the factor levels into finite subintervals
and consider only the grid points as its discrete levels (see, for example, Yang et al.
(2013)).
For example, suppose a continuous factor xi takes values in a continuous but bounded
region, its grid points can be selected equidistantly. If xi should be regarded as unbounded
(which is actually not common in practice), for example, xi ∈ [0, ∞) or (−∞, ∞), a
transformation such as ex /(1 + ex ) may be used to transform the original unbounded
region into a bounded one, like from [0, ∞) to [1/2, 1) or from (−∞, ∞) to (0, 1). Then
the grid points may be selected equidistantly within the transformed region [1/2, 1) or
(0, 1). It works reasonably well with a moderate number of continuous factors (see
Yang et al. (2013) for a discussion on the efficiency of optimal designs based on gridpoint experimental settings).
Once a set of grid points is chosen for each continuous factor, the design problem
becomes an allocation problem on a finite set of design points. All the previous algorithms
and results can be applied accordingly. We use the example of house flies for illustrating
purpose.
Example 5.2: Emergence of house flies (continued) Recall that the sole factor, dose
of radiation, is a continuous factor in nature. In the emergence of house flies experiment,
seven levels ranging from 80 to 200 equidistantly were used, that is, with grid distance
20. Suppose a followup experiment is considered and dose levels with grid distance 5 are
feasible (which technically depends on the sensitivity of the radiation device). Then there
are 25 dose levels available as 80, 85, 90, . . . , 195, 200. The D-optimal approximate design
given the 25 grid-5 design points concentrates on five design points 80, 120, 125, 155, 160
with optimal allocation 0.3163, 0.1429, 0.2003, 0.1683, 0.1723, respectively. If we further
consider the set of 121 grid-1 dose levels, the D-optimal design is supported on five design
points 80, 122, 123, 157, 158 with allocations 0.3163, 0.0786, 0.2636, 0.2206, 0.1209. The
optimal design seems to converge to a three-point design as the grid points are finer and
finer. Actually, if we reallocate the grid-1 D-optimal design, denoted as p1 , into a minimally supported design p1m at dose levels 80, 123, 157 with weights 0.3163, 0.3422, 0.3415,
then the efficiency of p1m compared with p1 is as high as (1503272/1504236)1/5 = 99.99%.


6

Discussion

For EW and Bayesian D-optimal designs, the choice of prior on the parameter space Θ
is critical. Suppose the parameters are θ1 , . . . , θp with individual ranges θi ∈ Ii . As a
common practice for generalized linear models, a uniform prior or independent normal
21

prior may be assumed for a rectangular type of domain I1 × · · · × Ip (Yang et al., 2016).
One issue with multinomial logit models is that the feasible domain of Θ may not be rectangular, at least for cumulative logit models. For the two examples in Section 5.6, we use
bootstrapping strategy to obtain an empirical prior. Zocchi and Atkinson (1999) adopted
a multivariate normal prior with mean vector and covariance structure estimated from
previous experimental data. For the two examples we study, our results are essentially
the same when we use the multivariate normal prior. However, when we use uniform or
independent normal prior, EW D-optimal design is not as good as Bayesian D-optimal
one in terms of robustness. Further investigation would be done towards choices of prior
and criteria other than EW optimality.
When a pilot study or experimental data is not available for locally or EW optimal
designs, sequential design or multistage design may be used to obtain an initial guess or
a reasonable prior for unknown parameters. We refer the readers to a review paper for
generalized linear models (Khuri et al., 2006), which is also applicable for multinomial
logistic models or multivariate generalized linear models.
When confronted with model uncertainty, npo versus ppo for instance, our work may
provide the experimenter the option to choose a design, for instance, that is highly
efficient for both ppo and npo models. Further investigations along this line would be
practically meaningful.

A
A.1

Appendix section
Constant matrix L in equation (1)

1
 1






Lbaseline = 
0 0
0 0

. .
 .. ..

0 0
1 1

1







Lcontinuation = 
0
0

.
 ..

0
1

..

.
1
0
0
..
.

···
···
..
.
···
···

0
1

1
..

.

1 ···
0 1
.. . .
.
.
0 ···
1 ···



0
1 0 ···
0
1 1 0 · · ·
0


 .. .. . .
.. 
. .
.
.


1 1 · · ·
0
1



1
,
L
=
0
1
·
·
·
1
cumulative


0 0 1 · · ·
1


. . .
.. 
.. ...

 .. ..
.

0 0 · · ·
1
0
1
1 1 ···
1


0
1
 1
0



.. 
..

.
.



1 0
1



· · · 1
0
1
,
L
=
adjacent


0
· · · 1
1


.
.. 
..
..
 ..
. .
.


0
0 1
1 1
1 1 ··· 1

22


0
0

.. 
.

0

1

1

.. 
.

1
1

0
0

.. 
.

0








1
1

A.2

J−1
M(HTj ))] in TheFormula for calculating dim[M(HTc ) ∩ (∩j=1
orem 3.5

Theorem A.1. Suppose Hi is of pi × m with rank ri , i = 1, . . . , n. Denote rank((HTi1 ,
. . . , HTik )) = ri1 ,...,ik for any 2 ≤ k ≤ n and 1 ≤ i1 < · · · < ik ≤ n. Then
!
n
n
X
\
X
X
ri1 ,i2 ,i3 − · · · + (−1)n−1 r1,2,...,n (20)
ri1 ,i2 +
dim
M(HTi ) =
ri −
i=1

i=1

i1 <i2

i1 <i2 <i3

The proof of Theorem A.1 is relegated to the Supplementary Materials.

A.3

Reparametrization and D-optimality (continued)

In general, let θ = (θ1 , . . . , θp )T be one set of parameters and ϑ = (ϑ1 , . . . , ϑp )T be
another set of parameters, such that, θl = hl (ϑ), l = 1, · · · , p; the map θ = θ(ϑ) =
(h1 (ϑ), . . . , hp (ϑ))T is one-to-one; hl ’s are differentiable; and the p × p Jacobian matrix
J = (hi (ϑ)/∂ϑj )ij is nonsingular.
Consider a design ξ = {(xi , wi ), i = 1, . . . , m} with the distinct experimental settings
xi ’s and the corresponding proportions wi ∈ [0, 1]. According to Schervish (1995, page
115), the Fisher information matrix Fξ (ϑ) at ϑ and the Fisher information matrix Fξ (θ)
at θ = θ(ϑ) satisfy Fξ (ϑ) = JT Fξ (θ(ϑ))J. Then |Fξ (ϑ)| = |J|2 · |Fξ (θ(ϑ))|, where J
contains no design points but parameters. A locally D-optimal design maximizing |Fξ (ϑ)|
also maximizes |Fξ (θ(ϑ))|. That is, it is mathematically equivalent to find D-optimal
designs for parameters ϑ or θ.
In terms of Bayesian D-optimal criterion, if a prior distribution of ϑ is available, it
induces a prior distribution of θ since θ = θ(ϑ) is one-to-one. Then Eϑ log |Fξ (ϑ)| =
Eϑ log JT Fξ (θ(ϑ))J = Eϑ log |J|2 + Eϑ log |Fξ (θ(ϑ))| = Eϑ log |J|2 + Eθ log |Fξ (θ)|.
Therefore, a Bayesian D-optimal design that maximizes Eθ log |Fξ (θ)| also maximizes
Eϑ log |Fξ (ϑ)|.

References
Agresti, A., 2010: Analysis of Ordinal Categorical Data. Wiley, 2nd ed.
—, 2013: Categorical Data Analysis. Wiley, 3rd ed.
Atkinson, A., A. Donev, and R. Tobias, 2007: Optimum Experimental Designs, with
SAS . Oxford University Press.
Chaloner, K. and I. Verdinelli, 1995: Bayesian experimental design: a review. Statistical
Science, 10, 273–304.
Christensen,
R.,
2015:
Analysis
of
ordinal
data
with
cumulative
link
models
–
estimation
with
the
r-package
ordinal.
http://cran.r-project.org/web/packages/ordinal/vignettes/clm_intro.pdf.

23

Chuang-Stein, C. and A. Agresti, 1997: Tutorial in biostatistics-a review of tests for
detecting a monotone dose-response relationship with ordinal response data. Statistics
in Medicine, 16, 2599–2618.
Dobson, A. and A. Barnett, 2008: An Introduction to Generalized Linear Models. Chapman & Hall/CRC, 3rd ed.
Fedorov, V., 1972: Theory of Optimal Experiments. Academic Press.
Gantmacher, F., 1960: The Theory of Matrices. Chelsea.
Glonek, G. and P. McCullagh, 1995: Multivariate logistic models. Journal of the Royal
Statistical Society, Series B , 57, 533–546.
Itepan, N., 1995: Aumento do periodo de aceitabilidade de pupas de Musca domestica
L., 1758 (Diptera: Muscidae), irradiadas com raios gama, como hospedeiras de parasitoides (Hymenoptera: Pteromalidae). Master’s thesis, Centro de Energia Nuclear na
Agricultura/USP, Piracicaba, SP, Brazil.
Jennett, B. and M. Bond, 1975: Assessment of outcome after severe brain damage.
Lancet, 305, 480484.
Khuri, A., B. Mukherjee, B. Sinha, and M. Ghosh, 2006: Design issues for generalized
linear models: A review. Statistical Science, 21, 376–399.
Kovacs, I., D. Silver, and S. Williams, 1999: Determinants of block matrices and schurs
formula. American Mathematical Monthly, 106, 950–952.
Lall, R., M. Campbell, S. Walters, and K. Morgan, 2002: A review of ordinal regression
models applied on health-related quality of life assessments. Statistical Methods in
Medical Research, 11, 4967.
Liu, I. and A. Agresti, 2005: The analysis of ordered categorical data: An overview and
a survey of recent developments. Test, 14, 1–73.
McCullagh, P., 1980: Regression models for ordinal data. Journal of the Royal Statistical
Society, Series B , 42, 109–142.
McCullagh, P. and J. Nelder, 1989: Generalized Linear Models. Chapman and Hall/CRC,
2nd ed.
Perevozskaya, I., W. Rosenberger, and L. Haines, 2003: Optimal design for the proportional odds model. The Canadian Journal of Statistics, 31, 225–235.
Peterson, B. and F. Harrell, 1990: Partial proportional odds models for ordinal response
variables. Journal of the Royal Statistical Society. Series C (Applied Statistics), 39,
205–217.
Randall, J., 1989: The analysis of sensory data by generalised linear model. Biometrical
Journal , 31, 781–793.
24

Schervish, M., 1995: Theory of Statistics. Springer.
Seber, G., 2008: A Matrix Handbook for Statisticians. Wiley.
Silvey, S., 1980: Optimal Design. Chapman & Hall/CRC.
Stoica, P. and T. Marzetta, 2001: Parameter estimation problems with singular information matrices. IEEE Transactions on Signal Processing, 49, 87–90.
Stokes, M., C. Davis, and G. Koch, 2012: Categorical Data Analysis Using SAS . SAS
Institute, 3rd ed.
Stufken, J. and M. Yang, 2012: Optimal designs for generalized linear models. In Hinkelmann, K., ed., Design and Analysis of Experiments, Volume 3: Special Designs and
Applications, Wiley, chap. 4, pp. 137–165.
Tong, L., H. Volkmer, and J. Yang, 2014: Analytic solutions for d-optimal factorial
designs under generalized linear models. Electronic Journal of Statistics, 8, 1322–1344.
Wu, F.-C., 2008: Simultaneous optimization of robust design with quantitative and ordinal data. International Journal of Industrial Engineering: Theory, Applications and
Practice, 5, 231–238.
Yang, J. and A. Mandal, 2015: D-optimal factorial designs under generalized linear
models. Communications in Statistics - Simulation and Computation, 44, 2264–2277.
Yang, J., A. Mandal, and D. Majumdar, 2016: Optimal designs for 2k factorial experiments with binary response. Statistica Sinica, 26, 385–411.
Yang,
J.,
L.
Tong,
and
A.
Mandal,
2017:
D-optimal
designs with ordered categorical data. Statistica Sinica,
to appear.
http://www3.stat.sinica.edu.tw/preprint/SS-2016-0210_Preprint.pdf.
Yang, M., S. Biedermann, and E. Tang, 2013: On optimal designs for nonlinear models:
a general and efficient algorithm. Journal of the American Statistical Association, 108,
1411–1420.
Yee, T., 2015: Vector Generalized Linear and Additive Models With an Implementation
in R. Springer.
Zocchi, S. and A. Atkinson, 1999: Optimum experimental designs for multinomial logistic
models. Biometrics, 55, 437–444.

25

D-optimal Designs for Multinomial Logistic Models
Xianwei Bu1 , Dibyen Majumdar2 and Jie Yang2
1 AbbVie

Inc. and 2 University of Illinois at Chicago
Supplementary Materials

S.1

List of notations

0k

A vector of k zeros

1D

diag{1T , . . . , 1T }, (J − 1) × m(J − 1)

aj

hTj (x)β j + hTc (x)ζ, j = 1, . . . , J − 1, given x = (x1 , . . . , xd )T

bj

Coefficients in representing fi (z), j = 0, . . . , J − 1

BJ

J × J constant matrix used for deriving the coefficients of fi (z), (st−1 )st

C

J × (2J − 1) constant matrix, same for all the four logit models

c

Vector used for deriving coefficients of fi (z), (c1 , . . . , cJ−1 )T

cij

−1
J × 1 vectors such that (CT D−1
= (ci1 , . . . , ciJ )
i L)

cj

(j + 1)p j J−1−p fi (1/(j + 1)) − j J−1 fi (0), j = 1, . . . , J − 1

cα1 ,...,αm

αm
Coefficient of w1α1 · · · wm
in the determinant of GT WG

d

Total number of design factors

ds

ds = (fij (s) − fij (0))/s, s = 1, . . . , q, for coefficients in fij

Di

diag(Lπ i )

F

Fisher information matrix of the design, F =

f

f (w) = f (w1 , . . . , wm ) =P|GT WG| which is proportional to |F|; or
f (n) = f (n1 , . . . , nm ) = | m
i=1 ni Fi | = |F|

Fi

Pm

i=1

ni Fi

Fisher information matrix at the ith design point

fi

fi (z) = f (w1 (1−z)/(1−wi ), . . . , wi−1 (1−z)/(1−wi ), z, wi+1 (1−z)/(1−
wi ), . . . , wm (1 − z)/(1 − wi )) with 0 ≤ z < 1

fij

fij (z) = f (n1 , . . . , ni−1 , z, ni+1 , . . . , nj−1 , ni + nj − z, nj+1 , . . . , nm ) with
z = 0, 1, . . . , ni + nj

G

Matrix component for Fisher information matrix such that F = nGT WG,
mJ × p
S1

gs

T
g0 = fij (0) and (g1 , . . . , gq )T = B−1
q (d1 , . . . , dq )

H

Matrix component for Fisher information matrix such that F = HUHT ,
consisting of H1 , . . . , HJ−1 and possibly Hc , p × m(J − 1)

Hc

Matrix for the common component of J−1 categories, (hc (x1 ), . . . , hc (xm )),
pc × m

hc (xi )

Vector of common predictors for all of the response categories as known
functions of the ith experimental setting, (h1 (xi ), . . . , hpc (xi ))T

HD

diag{H1 , . . . , HJ−1}, (p1 + · · · + pJ−1 ) × m(J − 1)

hD (xi )

diag{hT1 (xi ), . . . , hTJ−1 (xi )}, (J − 1) × (p1 + · · · + pJ−1 )

Hrow

[Hc , . . . , Hc ], pc × m(J − 1)

hTrow (xi )

[hc (xi ), . . . , hc (xi )]T , (J − 1) × pc

Hj

Matrix for the jth category only, (hj (x1 ), . . . , hj (xm )), pj × m

hj (xi )

Vector of predictors for the jth response category as known functions of
the ith experimental setting, (hj1 (xi ), . . . , hjpj (xi ))T

Ik

The identity matrix of order k

J

Total number of response categories

kmin

Smallest possible #{i | αi > 0} such that cα1 ,...,αm > 0

L

Constant (2J − 1) × J matrix, different for the four logit models

m

Total number of distinct experimental settings or design points

M(H)

Column space of matrix H, that is, the linear subspace spanned by the
columns of H

n

Total number of experimental units, n = n1 + · · · + nm

n

Allocation of experimental units, (n1 , . . . , nm )T , ni ≥ 0,

ni

Number of replicates at the ith experimental setting

p

Total number of parameters

pc

Number of common parameters for J − 1 categories

T
dim ∩J−1
j=1 M(Hj )

pH
pj
q

P

i

Number of parameters for the jth category only

min{2J − 2, p − kmin + 2, p}, upper bound of order of fij (z)
S2

ni = n

S
S+

Collection of all feasible P
approximate allocations, {(w1 , . . . , wm )T ∈
Rm | wi ≥ 0, i = 1, . . . , m; m
i=1 w1 = 1}
Collection of approximate allocations, {w ∈ S | f (w) > 0}

U

Block matrix (Ust )s,t=1,...,J−1 , m(J − 1) × m(J − 1)

Ust

diag{n1 ust(π 1 ), . . . , nm ust (π m )}, m × m

ust (π i )

cTis diag(π i )−1 cit for s, t = 1, . . . , J − 1

w

Real-valued allocation of experimental units, (w1 , . . . , wm )T , wi ≥ 0,
P
i wi = 1

W

diag{w1 diag(π 1 )−1 , . . . , wm diag(π m )−1 }, mJ × mJ

wi

Proportion of experimental units assigned to the ith experimental setting, ni /n

wu

Uniform allocation, (1/m, . . . , 1/m)T

X

Design space, the collection of all design points yielding strictly positive
categorical probabilities of response; or a predetermined set of design
points considered

xi

The ith distinct experimental setting or design point, (xi1 , . . . , xid )T

Xi

Model matrix at the ith design point, J × p, the last row is all 0’s

βj

Vector of parameters for the jth response category only, (βj1 , . . . , βjpj )T

γij

The cumulative probability from the 1th to jth catogory at the ith
experimental setting, γij = πi1 + · · · + πij .

ζ

Vector of common parameters for all of the response categories, (ζ1 , . . . , ζpc )T

ηi

Vector of linear predictors at the ith experimental setting, η i = (ηi1 , . . . , ηiJ )T =
Xi θ with ηiJ ≡ 0

θ

Vector of all parameters, p × 1

Θ

Parameter space, the collection of all feasible parameter vectors

πi

Vector of response category probabilities at the ith experimental setting.
π i = (πi1 , . . . , πiJ )T , πi1 + · · · + πiJ = 1

Λ(α1 , . . . , αm )
{(i1 , . . . , ip ) | 1 ≤ i1 < · · · < ip ≤ mJ; #{l : (k − 1)J < il 6
kJ} = αk , k = 1, . . . , m}
πij

Probability that the response falls into the jth catogory at the ith experimental setting

φ

Bayesian criterion, φ(p) = E(log |F|)
S3

S.2

Formulas of matrix differentiation

According to Seber (2008, Chapter 17)),
∂y
∂xT
∂Ax
∂xT
∂z
∂xT
∂ log y
∂xT

=



∂yi
∂xj



ij

= A
=

∂z
∂y
· T
T
∂y
∂x

= [diag(y)]−1

∂y
∂xT

where x = (xi )i , y = (yi )i , z = (zi )i , and thus log y = (log yi )i are vectors, and A is a
constant matrix.

S.3

Proofs

Proof of Theorem 2.1:
Suppose for distinct xi , i = 1, · · · , m, we have independent multinomial response
Yi = (Yi1 , · · · , YiJ )T ∼ Multinomial(ni ; πi1 , · · · , πiJ )
where ni =

PJ

j=1 Yij .

Then the log-likelihood for the multinomial model is
l(θ) = log L(θ)
m
Y
ni !
YiJ
Yi1
πi1
· · · πiJ
= log
Y ! · · · YiJ !
i=1 i1
= constant +

m
X

YiT log π i

i=1

where log π i = (log πi1 , · · · , log πiJ )T . Then the score vector
∂l
∂θ T

=

m
X

YiT diag(π i )−1

i=1

m

∂π i
∂θ T

X ∂π i
∂l
∂l
= ( T )T =
( T )T diag(π i )−1 Yi
∂θ
∂θ
∂θ
i=1

S4

Using the formulas of matrix differentiation, we get

Lemma S.1.

∂π i
∂π i ∂η i
=
·
T
∂η Ti ∂θ T
∂θ
−1

∂η i
· Xi
=
∂π Ti
−1

∂[CT log(Lπ i )] ∂[log(Lπ i )] ∂[Lπ i ]
· Xi
·
·
=
∂[log(Lπ i )]T
∂[Lπ i ]T
∂π Ti
−1
= CT [diag(Lπ i )]−1 L
Xi
−1
T
π i T diag(π i )−1 (CT D−1
i L) Xi = 0

Proof of Lemma S.1: Recall that 1T π i = πi1 + · · · + πiJ = 1
of Xi is all 0; and



∗ ∗ ··· 0
∗ ∗ ···
∗ ∗ · · · 0
∗ ∗ · · ·



CT =  . .
L =  ..
,



.
.
0 0 ··· 1
1 1 ···

Then

−1
D−1
i = diag(Lπ i )

and


∗ 0 ···
0 ∗ · · ·

=  ..
.

0 0 ···

for each i; the last row

∗
∗



1

 

0
∗ 0 ··· 0


0 
 0 ∗ · · · 0
 =  ..

 

.
1
0 0 ··· 1
1T π i




∗ ··· ∗
∗
∗ · · · ∗
∗



D−1
 and CT D−1
.
i L = 
i L = 
.


. 
T
1


··· ∗
· · · ∗


..
. 
1T

−1
Rewrite (CT D−1
= (ci1 , · · · , ciJ ). Then 1T ci1 = · · · = 1T ci,J−1 = 0 and 1T ciJ = 1
i L)
T
−1
(just check the last row of CT D−1
= (1, · · · , 1), then
i L). Since π i diag(π i )
−1
π i T diag(π i )−1 (CT D−1
= (1, · · · , 1)(ci1 , · · · , ciJ ) = (0, · · · , 0, 1)
i L)
T
−1
Since the last row of Xi is all 0, then π i T diag(π i )−1 (CT D−1
i L) Xi = 0 .

As a direct conclusion of Lemma S.1,
m

E(

X
∂l
−1
T
)
=
ni π i T diag(π i )−1 (CT D−1
i L) Xi = 0
T
∂θ
i=1

S5



Then the Fisher information matrix (see, for example, Schervish (1995, Section 2.3.1))




∂l ∂l
∂l
∂l
F = Cov
=E
,
·
∂θ ∂θ
∂θ ∂θ T
!
m
m
X
X
∂π i T
∂π
j
= E
( T ) diag(π i )−1 Yi ·
YjT diag(π j )−1 T
∂θ
∂θ
i=1
j=1
!
m X
m
X
∂π
∂π i T
j
= E
( T ) diag(π i )−1 Yi YjT diag(π j )−1 T
∂θ
∂θ
i=1 j=1
Since Yi ’s follow independent multinomial distribution, then


2
ni (ni − 1)πi1
+ ni πi1 · · ·
ni (ni − 1)πis πit


..
..
..
E(Yi YiT ) = 

.
.
.
2
ni (ni − 1)πis πit
· · · ni (ni − 1)πiJ + ni πiJ
= ni (ni − 1)π i π Ti + ni diag(π i )

On the other hand, for i 6= j,
E(Yi YjT ) = E(Yi ) · E(YjT ) = ni nj π i π Tj
Then the Fisher information matrix
F=

m
X
∂π i
∂π i
( T )T diag(π i )−1 ni (ni − 1)π i π Ti diag(π i )−1 T
∂θ
∂θ
i=1

m
X
∂π i
∂π i
+
( T )T diag(π i )−1 ni diag(π i )diag(π i )−1 T
∂θ
∂θ
i=1
X ∂π i
∂π j
( T )T diag(π i )−1 ni nj π i π Tj diag(π j )−1 T
+
∂θ
∂θ
i6=j

, (a) + (b) + (c)

where
m
X
∂π i
∂π i
(b) =
( T )T diag(π i )−1 T ni
∂θ
∂θ
i=1
" m
#" m
#T
X ∂π i
X ∂π i
(a) + (c) =
( T )T diag(π i )−1 π i ni
( T )T diag(π i )−1 π i ni
∂θ
∂θ
i=1
i=1
m
X ∂π i
∂π i
( T )T diag(π i )−1 ni π i π i T diag(π i )−1 T
−
∂θ
∂θ
i=1

Actually, let
Ei = π i T diag(π i )−1

∂π i
−1
= π i T diag(π i )−1 (CT D−1
i L) Xi
∂θ T
S6

which is 0T for each i according to Lemma S.1. Then
(a) + (c) =

" m
X

ni ETi

i=1

#"

m
X

ni ETi

i=1

#T

−

m
X

ni ETi Ei = 0J×J

i=1



The arguments above have proved Theorem 2.1.

Proof of Theorem 3.1: Because the last row of Xi consists of all zeros, the entries
in the last row and last column of Ui actually won’t make any difference. In order to
simplify the notations in this proof, we rewrite
hji , hj (xi )
hci , hc (xi )
usti , ust(π i )
J−1
X
us·i ,
usti
u·ti ,
u··i ,

t=1
J−1
X

j = 1, . . . , J − 1; i = 1, . . . , m
i = 1, . . . , m
s, t = 1, . . . , J − 1; i = 1, . . . , m

usti

s=1
J−1
J−1
XX

s = 1, . . . , J − 1;

i = 1, . . . , m

t = 1, . . . , J − 1;

i = 1, . . . , m

usti

i = 1, . . . , m

s=1 t=1

Based on Corollary 3.1, when
information Fi = XTi Ui Xi =

u11i h1i hT1i

..

.

 uJ−1,1,i hJ−1,i hT1i
u·1i hci hT1i

Xi takes partial proportional odds form (2), the Fisher
···
u1,J−1,i h1i hTJ−1,i
u1·i h1i hTci
..
..
..
.
.
.
T
· · · uJ−1,J−1,i hJ−1,i hJ−1,i uJ−1·i hJ−1,i hTci
···
u·J−1,i hci hTJ−1,i
u··i hci hTci

Then the Fisher information matrix F =
m
P



ni u11i h1i hT1i


i=1

..


.
 P
m

ni uJ−1,1,i hJ−1,i hT1i

 i=1

m
P
ni u·1i hci hT1i
i=1

m
P

···
..

i=1

.

···
···

m
P

i=1

or simply









H1
..
Hc

.

···

HJ−1
Hc







U11
..
.
UJ−1,1

Pm

i=1

ni Fi =
m
P

ni u1,J−1,i h1i hTJ−1,i

i=1

..
.

ni uJ−1,J−1,i hJ−1,i hTJ−1,i
m
P

i=1

ni u·J−1,i hci hTJ−1,i

···
..
.
···

U1,J−1
..
.
UJ−1,J−1

S7





m
P

i=1






ni u1·i h1i hTci
..
.

ni uJ−1·i hJ−1,i hTci
m
P

i=1

HT1



..

ni u··i hci hTci

.
HTJ−1


HTc
.. 
. 
HTc













−1
Proof of Theorem 3.2: Recall that (CT D−1
= (ci1 · · · ciJ ) and ust (π i ) = cTis diag(π i )−1 cit ,
i L)
for s, t = 1, . . . , J − 1 and i = 1, . . . , m. Denote
 T

c11


..
.




T

cm1 
 T

 c12



.


..


C̃ = 

T
cm2 

 .
.. 
..
 ..
.
. 


cT

 1,J−1



.
.


.
T
cm,J−1 m(J−1)×mJ


n1 diag(π 1 )−1

..
and W̃ = 
.

nm diag(π m )−1

We claim that U = C̃W̃C̃T . Actually

n1 cT11 diag(π 1 )−1

..
.




..
..
C̃W̃ = 
.
.

n cT
−1
 1 1,J−1diag(π 1 )

..

.





mJ×mJ

nm cTm1 diag(π m )−1
..
.

nm cTm,J−1 diag(π m )−1

and


c11

T
C̃W̃C̃ = C̃W̃ 



=

U11
..
.

UJ−1,1

..

.
cm1

···
..
.
···

···
..
.
···

.
cm,J−1



U1,J−1

..
=U
.
UJ−1,J−1















c1,J−1
..






Note that W̃ is diagonal with positive diagonal entries. Thus W̃ is positive definite. By

S8

adjusting the rows, we can verify that rank(C̃) is the same as rank(C̃′ ), where

 T
c11

 ..

 .

 T

c1,J−1


T


c21


..


.


′
C̃ = 

T
c


2,J−1


..


.



cTm1 



.
.. 


T
cm,J−1



That is, C̃ has full row rank and thus U is positive definite.
Proof of Theorem 3.3:
Q
J−1
|V|, where
Lemma S.2. |U| = ( m
i=1 ni )


V11 · · · V1,J−1


..
..
V =  ...

.
.

VJ−1,1 · · · VJ−1,J−1

u111
· · · u1,J−1,1

..
.



u11m
 .
..

.
= .
.
u
· · · uJ−1,J−1,1
 J−1,1,1

.
..

uJ−1,1,m

..

.
···
..

u1,J−1,m
..
.

.
· · · uJ−1,J−1,m














Kovacs et al. (1999) generalized Schur’s Formula (Gantmacher (1960)) as follows:
Lemma S.3. (Kovacs et al., 1999, Theorem 1)
Assume that M is a k × k block matrix with each block element Aij as an n × n matrix.


A11 · · · A1k

.. 
..
M =  ...
.
. 
Ak1 · · · Akk

If all of Aij ’s commute pairwise , that is, Aij Alm = Alm Aij for all possible pairs of
indices i, j and l, m. Then
|M| =

X

π∈Sk

(sgnπ)A1π(1) A2π(2) · · · Akπ(k)

Here the sum is computed over all permutations π of {1, 2, ..., k}.
S9

(S.1)

In our case, all of Vij ’s are diagonal matrices, so they commute pairwise. Moreover,
the sum of product matrices in Equation (S.1) is a diagonal matrix, in which each element
is the sum of products of the corresponding elements in those matrices. If we apply the
above lemma, we get
|V| =
=

X

π∈SJ −1
m
Y

(sgnπ)V1π(1) V2π(2) · · · VJ−1,π(J−1)

X

i=1 π∈SJ −1

(sgnπ)u1π(1)i u2π(2)i · · · uJ−1,π(J−1),i

Then the following result is obtained:
Q
Lemma S.4. |V| = m
i=1 |Vi |, where

u11 (π i )

..
Vi = 
.


· · · u1,J−1 (π i )

..
..

.
.
uJ−1,1 (π i ) · · · uJ−1,J−1(π i )

Note that Vi defined above is very similar to Ui define in Equation (6).
Q
−1
J
−2
Lemma S.5. |Vi | =
π
· |CT D−1
i L| .
j=1 ij

Proof of Lemma S.5: It can be verified that ciJ = π i .
Since cTij diag(π i )−1 ciJ = cTij 1 = 0 for j = 1, . . . , J − 1 and 1 for j = J, then


h
h
−1 iT
−1 i
Vi 0
T
−1
−1
T
−1
C Di L
diag(π i )
C Di L
=
0T 1



Combining Lemmas S.2, S.4, & S.5, we obtain Theorem 3.3.
Remark S.1. Actually, we provide an explicitQformula for |CT D−1
i L| in (S.18), which
J
can further clarify Lemma S.5 as (1) |Vi | = j=1 πij for baseline-category, adjacent−1 QJ−1 −1 2
2
categories, and continuation-ratio logit models; (2) |Vi | = πiJ
j=1 πij γij (1 − γij ) for
cumulative logit models.
Proof of Theorem 3.5:
The simplest case is npo models whose conclusion is straightforward.
The ppo model is the most general case. In this case, we consider a sequence of linear
subspaces
T
T
{0} ⊂ M(HTc ) ∩ (∩J−1
j=1 M(Hj )) ⊂ M(Hc )
with corresponding dimensions 0 ≤ rc − r0 ≤ rc , rank(Hc ), where r0 = rank(Hc ) −
T
dim[M(HTc ) ∩ (∩J−1
j=1 M(Hj ))]. Then there exist α1 , · · · , αrc −r0 , αrc −r0 +1 , · · · , αrc ∈
T
Rm s.t. {α1 , · · · , αrc −r0 } forms a basis of M(HTc ) ∩ (∩J−1
j=1 M(Hj )) and {α1 , · · · , αrc }
S10

forms a basis of M(HTc ), and by simple operations Hc can be transformed into H∗c =
(α1 , · · · , αrc , 0, · · · , 0)T and Hj can be transformed into
(j)

H∗j = (α1 , · · · , αrc −r0 , αrc −r0 +1 , · · · , αr(j)
, 0, · · · , 0)T
j
where rj = rank(Hj ), j = 1, 2, · · · , J − 1. Then rank(Hppo ) = rank(H∗ppo) with


H∗1


..


.
H∗ppo = 

∗

HJ−1 
H∗c · · · H∗c
p×m(J−1)

Since the first rc − r0 rows of (H∗c , · · · , H∗c ) can be eliminated by applying row operations of H∗j onto it separately, then rank(H∗ppo) = rank(H∗∗
ppo ) where


H∗1


..


.
∗∗
Hppo = 

∗

HJ−1 
H∗∗
· · · H∗∗
c
c
p×m(J−1)

T
∗∗
and H∗∗
c = (0, · · · , 0, αrc −r0 +1 , · · · , αrc , 0, · · · , 0) . Therefore, rank(Hppo ) = rank(Hppo ) 6
r1 + · · · + rJ−1 + r0 .
We claim that the nonzero rows of H∗∗
ppo are linearly independent which will lead
(j)
to the final conclusion. Actually, let’s denote those nonzero rows of H∗∗
ppo as Λi , i =
(j)
1, 2, · · · , rj , j = 1, 2, · · · , J−1 and Λrc −r0 +1 , · · · , Λrc , where Λi is the ith row of (0, · · · , 0, H∗j ,
(j)
∗∗
0, · · · , 0), and Λi is the ith row of (H∗∗
∈ R, i =
c , · · · , Hc ). Suppose there exist ai
1, 2, · · · , rj , j = 1, 2, · · · , J − 1 and ai ∈ R, i = rc − r0 + 1, · · · , rc s.t.

0=

rj
J−1 X
X

(j) (j)
ai Λi

+

j=1 i=1

then for j = 1, . . . , J − 1,
0=

rX
c −r0

(j)
ai αi

+

i=1

rj
X

i=rc −r0 +1

ai αi = −

(j) (j)
ai αi

(j)
ai αi

i=1

−

rj
X

+

rc
X

ai αi

i=rc −r0 +1

i=rc −r0 +1

rX
c −r0

ai Λi

i=rc −r0 +1

which implies for j = 1, . . . , J − 1,
rc
X

rc
X

(j)

i=rc −r0 +1

(j)

ai αi ∈ M(HTc ) ∩ M(HTj )


Prc
T
∩J−1
Thus, i=rc −r0 +1 ai αi ∈ M(HTc )
j=1 M(Hj ) . Then we must have
i=rc −r0 +1 ai αi =
0 since {αrc −r0 +1 , . . . , αrc } and {α1 , . . . , αrc −r0 } are linearly independent. Therefore,
ai = 0 for i = rc − r0 + 1, . . . , rc and thus
Prc

T

0=

rX
c −r0
i=1

(j)
ai αi

+

rj
X

i=rc −r0 +1

S11

(j)

(j)

ai αi

(j)

(j)

It implies ai = 0, i = 1, . . . , rc − r0 , rc − r0 + 1, . . . , rj since {α1 , . . . , αrc −r0 , αrc −r0 +1 ,
(j)
. . . , αrj } are linear independent.
Therefore, the conclusion on ppo models is justified.
Since po models are special cases of ppo models, the corresponding result is a direct
conclusion.

Proof of Theorem A.1:
Recall that dim(M(HTi )) = rank(HTi ) = ri and dim(M(HTi1 ) + · · · + M(HTik )) =
dim(M((HTi1 , · · · , HTik ))) = rank((HTi1 , · · · , HTik )) = ri1 ,...,ik , for i1 < · · · < ik and k =
2, . . . , n, where “+” stands for the sum of two linear subspaces.
First of all, dim(M(HT1 ) ∩ M(HT2 )) = dim(M(HT1 )) + dim(M(HT2 )) − dim(M(HT1 ) +
M(HT2 )) = r1 + r2 − r12 . That is, (20) is true for n = 2.
Suppose (20) is true for n = k. Then for n = k + 1,
T
k
T
T
dim(∩k+1
i=1 M(Hi )) = dim(∩i=1 M(Hi ) ∩ M(Hk+1 ))
= dim(∩ki=1 M(HTi )) + dim(M(HTk+1)) − dim(∩ki=1 M(HTi ) + M(HTk+1))

=

k
X
i=1

ri −

X

16i1 <i2 6k

ri1 i2 + · · · + (−1)k−1 r12···k + rk+1 − △

where
△ = dim(∩ki=1 M(HTi ) + M(HTk+1)) = dim(∩ki=1 M((HTi , HTk+1)))
=

k
X

rank((HTi , HTk+1 )) −

k
X

ri,k+1 −

i=1

+ · · · + (−1)

=

i=1

k−1

X

rank((HTi1 , HTk+1, HTi2 , HTk+1))

16i1 <i2 6k
T
rank((H1 , HTk+1, · · ·

X

16i1 <i2 6k

, HTk , HTk+1 ))

ri1 ,i2 ,k+1 + · · · + (−1)k−1 r1,2,...,k+1

Therefore,
T
dim(∩k+1
i=1 M(Hi ))
k
X
X
=
ri −
ri1 i2 + · · · + (−1)k−1 r12···k + rk+1
16i1 <i2 6k

i=1

−
=

k
X

ri,k+1 +

i=1

k+1
X
i=1

ri −

X

16i1 <i2 6k

X

16i1 <i2 6k+1

ri1 ,i2 ,k+1 + · · · + (−1)k r1,2,...,k+1

ri1 i2 + · · · + (−1)(k+1)−1 r1,2,...,k+1

That is, (20) is true for n = k + 1. By mathematical induction, (20) is true for general
n.

Proof of Corollary 3.3:
S12

Suppose pH > 0. Then there exist m × 1 vectors α1 , · · · , αpH , which form a basis
T
T
of ∩J−1
j=1 M(Hj ). Write Hc = (γ1 , · · · , γpc ) . According to Theorem3.5, if |F| > 0,
T
then r0 = rank(Hc ) = pc , or equivalently, M(HTc ) ∩ ∩J−1
j=1 M(Hj ) = {0}. Then

α1 , · · · , αpH , γ1, · · · , γpc are linearly independent. Thus m ≥ pc + pH .

Proof of Theorem 4.1:
Actually, according to Theorem 3.1, F = HUHT . From the proof of Theorem 3.2,
U = C̃W̃C̃T , where W̃ is a diagonal matrix. Therefore, F = HC̃W̃C̃T HT which leads
to the final result.

Proof of Lemma 4.1: Actually, max1≤i≤m αi ≤ J. Suppose max1≤i≤m αi ≥ J, which
means max1≤i≤m αi = J. Without any loss of generality, we assume α1 = J. Then ij = j
for j = 1, . . . , J.
According to the proof of Lemma S.1, we have 1T cij = 0 for i = 1, . . . , m and
j = 1, . . . , J − 1. Then 1T (c11 + · · · + c1,J−1 ) = 0 and thus 1T G[i1 , . . . , iJ ] = 0. That is,
rank(G[i1 , . . . , iJ ]) ≤ J − 1. Therefore, rank(G[i1 , . . . , ip ]) ≤ p − 1 and |G[i1 , . . . , ip ]| = 0.


Proof of Theorem 4.3: Suppose cα1 ,...,αm 6= 0 for some (α1 , . . . , αm ). Therefore, there
exist (i1 , . . . , ip ) ∈ (α1 , . . . , αm ) such that G[i1 , . . . , ip ] is of full rank p. Without any
loss of generality, we assume α1 ≥ · · · ≥ αk > 0 = αk+1 = · · · = αm , that is, {i | αi >
0} = {1, . . . k}. Consider the submatrix G̃ := G[1, . . . , kJ] which is kJ × p and contains
G[i1 , . . . , ip ] as a submatrix. Then G̃ is of rank p or G̃T is of full row rank p. Write
W̃ = k −1 diag{diag(π 1 )−1 , . . . , diag(π k )−1 }. Then the restricted matrix F := n G̃T W̃G̃
is positive definite. On the other hand, F is the Fisher information matrix nGT WG
as defined in (9) with w1 = · · · = wk = P
1/k and wk+1 = · · · = wm =
According to
P0.
k
k
−1
Theorem 4.1 and Theorem 2.1, F = nk
i=1 Fi is positive
i=1 Fi . Therefore, Fres :=
definite.


Proof of Theorem 5.1
Case 1: Baseline-category logit model for nominal response
The baseline-category logit model for nominal response (Agresti, 2013; Zocchi and Atkinson,
1999) can be extended in general as follows
 
πij
log
= hTj (xi )β j + hTc (xi )ζ , j = 1, . . . , J − 1
(S.2)
πiJ
Lemma S.6. Fixing xi , β j , j = 1, · · · , J − 1 and ζ in Model (S.2), let aj = hTj (xi )β j +
hTc (xi )ζ, j = 1, . . . , J − 1. Then 0 < πij < 1, j = 1, . . . , J exist uniquely if and only if
−∞ < aj < ∞, j = 1, . . . , J − 1. In this case,
πij =



ea j
a
1
e +···+eaJ −1 +1
1
ea1 +···+eaJ −1 +1

1 ≤j ≤ J −1
j=J

(S.3)

Proof of Lemma S.6: Write yj = log πij , j = 1, . . . , J. Then 0 < πij < 1, j = 1, . . . , J
if and only if yj ∈ (−∞, 0), j = 1, . . . , J. In this case, Model (S.2) implies aj = yj − yJ ∈
(−∞, ∞), j = 1, . . . , J − 1.
S13

On the other hand, for any given a1 , . . . , aJ−1 ∈ (−∞, ∞), yj = aj + yJ , j = 1, . . . , J − 1.
Note that
1 =
=
=
=

πi1 + πi2 + · · · + πi,J−1 + πiJ
ey1 + ey2 + · · · + eyJ −1 + eyJ
ea1 +yJ + ea2 +yJ + · · · + eaJ −1 +yJ + eyJ
eyJ (ea1 + ea2 + · · · + eaJ −1 + 1)

Since πij = eyj , we get solutions of πij given in (S.3), and thus πij ∈ (0, 1) exists and is
unique, j = 1, . . . , J.
#
Case 2: Cumulative logit model for ordinal response
The cumulative logit model for ordinal responses (McCullagh, 1980; Christensen,
2015) can be described in general as follows:


πi1 + · · · + πij
log
= hTj (xi )β j + hTc (xi )ζ , j = 1, . . . , J − 1
(S.4)
πi,j+1 + · · · + πiJ
Lemma S.7. Fixing xi , β j , j = 1, · · · , J − 1 and ζ in Model (S.4), let aj = hTj (xi )β j +
hTc (xi )ζ, j = 1, . . . , J − 1. Then 0 < πij < 1, j = 1, . . . , J exist and are unique if and
only if −∞ < a1 < a2 < · · · < aJ−1 < ∞. In this case,
 exp(a1 )
j=1

 1+exp(a1 )
exp(aj−1 )
exp(aj )
(S.5)
πij =
− 1+exp(aj−1 ) 1 < j < J
1+exp(aj )


1
j=J
1+exp(aJ −1 )

Proof of Lemma S.7: Taking j = 1 in Model (S.4), then log (πi1 /(1 − πi1 )) = a1 and
πi1 = exp(a1 )/[1 + exp(a1 )]. Then 0 < πi1 < 1 if and only if −∞ < a1 < ∞. For
j = 2, · · · , J − 1,
exp(aj−1 )
exp(aj )
−
πij =
1 + exp(aj ) 1 + exp(aj−1 )
which implies that πij > 0 if and only if aj > aj−1 . Therefore, πiJ = 1 − (πi1 + · · · +
πi,J−1 ) = 1 − exp(aJ−1 )/[1 + exp(aJ−1 )] = 1/[1 + exp(aJ−1 )], which indicates 0 < πiJ < 1
if and only if −∞ < aJ−1 < ∞. Given πi1 + · · · + πiJ = 1, we have
−∞ < a1 < a2 < · · · < aJ−1 < ∞ ⇔ πij ∈ (0, 1),

j = 1, . . . , J
#

Corollary S.1. For proportional odds model


πi1 + · · · + πij
= βj + hTc (xi )ζ ,
log
πi,j+1 + · · · + πiJ

j = 1, . . . , J − 1

The design space has no restriction since
−∞ < β1 < β2 < · · · < βJ−1 < ∞ ⇐⇒ πij ∈ (0, 1),
S14

j = 1, . . . , J

(S.6)

Case 3: Adjacent-categories logit model for ordinal response
The adjacent-categories logit model for ordinal responses (Liu and Agresti, 2005;
Agresti, 2013) can be extended as follows:


πij
log
= hTj (xi )β j + hTc (xi )ζ , j = 1, . . . , J − 1
(S.7)
πi,j+1
Lemma S.8. Fixing xi , β j , j = 1, · · · , J − 1 and ζ in Model (S.7), let aj = hTj (xi )β j +
hTc (xi )ζ, j = 1, . . . , J − 1. Then 0 < πij < 1, j = 1, . . . , J exist uniquely if and only if
−∞ < aj < ∞, j = 1, . . . , J − 1. In this case,
(
exp(aJ −1 +···+aj )
j = 1, . . . , J − 1
exp(aJ −1 +···+a1 )+exp(aJ −1 +···+a2 )+···+exp(aJ −1 )+1
(S.8)
πij =
1
j=J
exp(aJ −1 +···+a1 )+exp(aJ −1 +···+a2 )+···+exp(aJ −1 )+1
Proof of Lemma S.8: Let yj = log πij . Then 0 < πij < 1, j = 1, . . . , J if and
only if yj ∈ (−∞, 0). In this case, Model (S.7) implies aj = yj − yj+1 ∈ (−∞, ∞),
j = 1, . . . , J − 1.
On the other hand, for any given a1 , . . . , aJ−1 ∈ (−∞, ∞), yj = (aJ−1 + · · · + aj ) + yJ ,
j = 1, . . . , J − 1. Note that
1 = πi1 + πi2 + · · · + πi,J−1 + πiJ
= ey1 + ey2 + · · · + eyJ −1 + eyJ

= eyJ eaJ −1 +···+a1 + eaJ −1 +···+a2 + · · · + eaJ −1 + 1

Since πij = eyj , we get solutions of πij given in (S.8), and thus πij ∈ (0, 1) exists and is
unique, j = 1, . . . , J.
#
Case 4: Continuation-ratio logit model for hierarchical response
The continuation-ratio logit model for hierarchical responses (Agresti, 2013; Zocchi and Atkinson,
1999) can be rewritten in general as follows:


πij
log
= hTj (xi )β j + hTc (xi )ζ , j = 1, . . . , J − 1
(S.9)
πi,j+1 + · · · + πiJ
Lemma S.9. Fixing xi , β j , j = 1, · · · , J − 1 and ζ in Model (S.9), let aj = hTj (xi )β j +
hTc (xi )ζ, j = 1, . . . , J − 1. Then 0 < πij < 1, j = 1, . . . , J exist uniquely if and only if
−∞ < aj < ∞, j = 1, . . . , J − 1. In this case,
 a Qj
(eas + 1)−1 j = 1, . . . , J − 1
ej
(S.10)
πij = QJ−1s=1as
−1
j=J
s=1 (e + 1)

Proof of Lemma S.9: Let yj = log πij . Then 0 < πij < 1, j = 1, . . . , J if and only if
yj ∈ (−∞, 0). In this case, Model (S.9) implies aj = yj − log(eyj+1 + · · · eyJ ) ∈ (−∞, ∞),
j = 1, . . . , J − 1.
On the other hand, for any given a1 , . . . , aJ−1 ∈ (−∞, ∞), it can be verified by induction
that
eyJ −1 = eyJ eaJ −1
eyJ −2 = eyJ eaJ −2 (eaJ −1 + 1)
eyj = eyJ eaj (eaj+1 + 1) · · · (eaJ −1 + 1) , j = J − 3, J − 4, · · · , 1
S15

Therefore, it can be verified that
1 = πi1 + πi2 + · · · + πi,J−1 + πiJ
= ey1 + ey2 + · · · + eyJ −1 + eyJ
= eyJ (ea1 + 1) (ea2 + 1) · · · (eaJ −1 + 1)
Since πij = eyj , we get solutions of πij given in (S.10), and thus πij ∈ (0, 1) exists and is
unique, j = 1, . . . , J.
#
Theorem 5.1 is obtained as a summary of Lemmas S.6, S.7, S.8, & S.9.

Proof of Corollary 5.1: We only need to verity the “only if” part. According to
Theorem 3.4, if f (w) > 0 for some w = (w1 , . . . , wm )T = (n1 , . . . , nm )T /n, then the
corresponding H∗ is of full row rank. Note that H∗ can be obtained from H after
removing the columns of H corresponding to ni = 0. Thus H is of full row rank too,
which corresponds to the uniform allocation. That is, f (wu ) > 0.
In this case, any w = (w1 , . . . , wm )T such that 0 < wi < 1, i = 1, . . . , m leads to f (w) > 0
since it corresponds to the same H matrix.

Proof of Theorem 5.3: According to Theorem 4.2,
X
fij (z) =
coefficient · z αi (ni + nj − z)αj
αi ≥0,αj ≥0,αi +αj ≤p

is a polynomial with nonnegative coefficients, whose order depends on the largest possible
αi + αj . Lemma 4.1 implies that max{αi , αj } ≤ J − 1 for positive coefficients and
Corollary 4.2 further implies that αi + αj ≤ p − (kmin − 2) = p − kmin + 2 for positive
coefficients. Therefore, fij (z) is at most an order-q polynomial of z.

Proof of Theorem 5.4: In this case, the model is essentially a generalized linear
model for binomial response with logit link. Theorem 4.2 says that the objective funcαm
.
tion f (w) = |GT WG| is an order-p polynomial consisting of terms cα1 ,...,αm w1α1 · · · wm
According to Lemma 4.1, cα1 ,...,αm 6= 0 only if αi ∈ {0, 1}, i = 1, . . . , m. Therefore,
in order to keep f (w) > 0, we must have m ≥ p. In other words, a minimally supported design may contain exactly m = p distinct design points or experimental settings.
In this case, the objective function f (w) ∝ w1 · · · wm and the D-optimal allocation is
w = (1/m, . . . , 1/m)T .

Proof of Corollary 5.2: According to Theorem 3.1, F = HUHT . In this case, there
exist m = p1 experimental settings such that rank(H) = p1 (J − 1) = p. On the other
hand, the minimum number of experimental settings is at least max{p1 , . . . , pJ−1 } = p1
based on Corollary 4.2. Therefore, the minimal number is m = p1 . In this case, H is a
square matrix and
!J−1
m
Y
|F| = |H|2 · |U| ∝
wi
i=1

according to Theorem 3.3. Thus, the uniform allocation wu = (1/m, . . . , 1/m)T is Doptimal in this case. Note that m = p1 < p1 (J − 1) = p.


Proof of Lemma 5.1: We actually claim more detailed conclusions as follows:
S16

(i) If c1 = c2 = c3 , then the solution is w1 = w2 = w3 = 1/3.
(ii) If c1 = c2 < c3 , then w1 = w2 >
pw3 > 0. Actually, w1 = w2 = (−2c1 + c3 + ∆1 )/D1
and w3 = c3 /D1 , where ∆1 = 4c21 − c1 c3 + c23 and D1 = −4c1 + 3c3 + 2∆1 .

(iii) If c1 < c2 = c3 , then w1 > w2 = p
w3 > 0. Actually, w1 = (−c1 + 2c3 + ∆2 )/D2 and
w2 = w3 = 3c3 /D2 , where ∆2 = c21 − c1 c3 + 4c23 and D2 = −c1 + 8c3 + ∆2 .

(iv) If c1 < c2 < c3 , then w1 > w2 > w3 > 0. The procedure of obtaining analytic
solutions of w1 , w2, w3 is as follows: (1) obtain y1 from (S.17); (2) obtain y2 from
(S.15); (3) w1 = y1 /(y1 + y2 + 1), w2 = y2 /(y1 + y2 + 1), w3 = 1/(y1 + y2 + 1).

First of all, we only need to consider the cases of 0 < wi < 1, i = 1, 2, 3 (otherwise,
f (w1 , w2 , w3 ) = 0). It can also been verified that 0 < c1 ≤ c2 ≤ c3 implies that w1 ≥
w2 ≥ w3 > 0 (otherwise, for example, if w1 < w2 , one may replace w1 , w2 both with
(w1 + w2 )/2 and strictly increase f ). The same argument implies that if ci = cj , then
wi = wj in the solution.
According to Theorem 5.10 in Yang et al. (2017), (w1 , w2 , w3 )T maximizes f (w1 , w2 , w3 )
if and only if
∂f
∂f
∂f
=
=
∂w1
∂w2
∂w3
which is equivalent to ∂f /∂w1 = ∂f /∂w3 and ∂f /∂w2 = ∂f /∂w3 and thus equivalent to
c3 w1 w2 (w1 − 2w3 ) + 2c2 w1 w3 (w1 − w3 ) = c1 w2 w3 (−2w1 + w3 )

(S.11)

c3 w1 w2 (w2 − 2w3 ) + 2c1 w2 w3 (w2 − w3 ) = c2 w1 w3 (−2w2 + w3 )

(S.12)

c3 y1 y2 (y1 − 2) + 2c2 y1 (y1 − 1) = c1 y2 (−2y1 + 1)
c3 y1 y2 (y2 − 2) + 2c1 y2 (y2 − 1) = c2 y1 (−2y2 + 1)

(S.13)
(S.14)

Following Yang et al. (2016b, Section 5.2), we denote y1 = w1 /w3 > 0 and y2 = w2 /w3 >
0. Actually, w1 ≥ w2 ≥ w3 > 0 implies y1 ≥ y2 ≥ 1. Since w1 + w2 + w3 = 1, it implies
w3 = 1/(y1 + y2 + 1), w1 = y1 /(y1 + y2 + 1), and w2 = y2 /(y1 + y2 + 1). Then (S.11) and
(S.12) are equivalent to

From (S.13) we get y2 [c3 y12 − 2(c3 − c1 )y1 − c1 ] = 2c2 y1 (1 − y1 ). If y1 = 1, then we must
have y2 = 1 and c3 − 2(c3 − c1 ) − c1 = 0, which implies w1 = w2 = w3 = 1/3 and
c1 = c2 = c3 . Actually, we can also verify that c1 = c3 implies y1 = 1.
Now we assume y1 > 1, which implies c1 < c3 . Then
y2 =

c3 y12

2c2 (1 − y1 )y1
− 2(c3 − c1 )y1 − c1

(S.15)

After plugging (S.15) into (S.14), we get
a0 + a1 y1 + a2 y12 + a3 y13 + y14 = 0

(S.16)

where a0 = c21 /c23 > 0, a1 = 4c1 (−2c1 + c2 + 2c3 )/(3c23 ) > 0, a2 = 2(2c21 − 2c1 c2 − 7c1 c3 −
2c2 c3 + 2c23 )/(3c23 ), and a3 = 4(2c1 + c2 − 2c3 )/(3c3 ).
S17

Denote h(y1 ) = a0 + a1 y1 + a2 y12 + a3 y13 + y14 . Note that h(∞) = ∞, h(−c1 /c3 ) =
−c21 (c21 +8c1 c2 −2c1 c3 +8c2 c3 +c23 )/(3c43 ) < 0, h(0) = c21 /c23 > 0, h(1) = −(c1 −c3 )2 /(3c23 ) <
0, and h(∞) = ∞. Then h(y1 ) = 0 yields four real roots in (∞, −c1 /c3 ), (−c1 /c3 , 0),
(0, 1), and (1, ∞), respectively. That is, there is one and only one y1 ∈ (1, ∞).
According to Tong et al. (2014, equation (12)),
√
√
a3
A1
C1
+
,
(S.17)
y1 = − +
4
2
2
where
2a2 a23
G1
A1 = −
+
+
,
3
4
3 × 21/3
4a2 a23
G1
−8a1 + 4a2 a3 − a33
√
C1 = −
+
−
+
,
3
2
3 × 21/3
4 A1
1/3 
1/3

q
q
2
3
2
3
G1 =
F1 − F1 − 4E1
+ F1 + F1 − 4E1
,
E1 = 12a0 + a22 − 3a1 a3 ,
F1 = 27a21 − 72a0 a2 + 2a32 − 9a1 a2 a3 + 27a0 a23 .

The calculation of G1 , A1 , C1 , and y1 are operations among complex numbers, while y1
at the end would be a real number.
The procedure of obtaining analytic solutions of w1 , w2 , w3 would be, (1) obtain y1
from (S.17); (2) obtain y2 from (S.15); (3) w1 = y1 /(y1 + y2 + 1), w2 = y2 /(y1 + y2 + 1),
w3 = 1/(y1 + y2 + 1).
Now we discuss some special cases.
(i) If c1 = c2 < c3 , thenpw1 = w2 and thus y1 = y2 . Both (S.13) and (S.14) yield
4c21 − c1 c3 + c23 ), which implies
y1 = c−1
3 (−2c1 + c3 +
w1 = w2 =

where ∆1 =

−2c1 + c3 + ∆1
,
−4c1 + 3c3 + 2∆1

w3 =

c3
−4c1 + 3c3 + 2∆1

p
4c21 − c1 c3 + c23 . Note that w1 > w3 since ∆1 > 2c1 .

(ii) If c1 < c2 = c3 ,pthen w2 = w3 and thus y2 = 1. From (S.13) we get y1 =
3c−1
c21 − c1 c3 + 4c23 ), which implies
3 (−c1 + 2c3 +
w1 =

where ∆2 =

−c1 + 2c3 + ∆2
,
−c1 + 8c3 + ∆2

w2 = w3 =

3c3
−c1 + 8c3 + ∆2

p
c21 − c1 c3 + 4c23 . Note that w1 > w2 since ∆2 > c1 + c3 .

(iii) If c1 < c2 < c3 , then y1 , y2 and thus w1 , w2 , w3 can be obtained analytically. We
have proven y1 ≥ y2 ≥ 1. Using (S.13) and (S.14), it can be verified that y1 6= y2
unless c1 = c2 ; and y2 6= 1 unless c2 = c3 . That is, y1 > y2 > 1 and w1 > w2 > w3 .



Proof of Theorem 5.5: (i) is straightforward. (ii) follows from the facts in the proof
of Lemma S.1, ciJ ≡ π i ; 1T cij = 0, j = 1, . . . , J − 1; and 1T ciJ = 1. (iii) and (iv) can be
verified using the formulas of cij in Section S.4.

S18

S.4

A key matrix for four different kinds of multinomial logistic
models

There are the four different kinds of multinomial logistic models in the literature: baselinecategory logit model for nominal responses, cumulative logit model for ordinal responses,
adjacent-categories logit model for ordinal responses, and continuation-ratio logit model
−1
for hierarchical responses. According to Theorem 2.1, (CT D−1
is a key matrix we
i L)
must calculate.
Recall that πi1 + · · · + πiJ = 1, i = 1, . . . , m. Then
 1

0 ···
0
− π1iJ
πi1


..
 0 π1 . . .
.
− π1iJ 
i2


.. 
 .. . .
.
(CT D−1
L)
=
.
baseline
i
.
.
 .
0
. 


1
 0 · · · 0 πi,J −1 − π1iJ 
1
1 ···
1
1
J×J

 1
1
1
1
− 1−γi1 − 1−γi1 · · ·
− 1−γ
γi1
i1
1
1
1

 1
− 1−γ
···
− 1−γ
γi2
 γi2
i2
i2 


..
..
..
..
..
(CT D−1

.
.
.
.
i L)cumulative =  .

 1
1
1
1

γ
···
− 1−γ
γ
γ
i,J −1

i,J −1

1

(CT D−1
i L)continuation






=



(CT D−1
i L)adjacent

1

1
πi1

1
− 1−γ
i1
1
πi2

0
..
.
0
1





=




..

.
···
1
1
πi1

0
..
.
0
1

− π1i2
1
πi2
..

.
···
1

i,J −1

···

1
− 1−γ
i1
1
− 1−γ
i2

···
···
..
.

0
···

1
πi,J −1

..

.

0
− π1i3
..
.
0
···

i,J −1

1

1
···
..
.
..
.
1
πi,J −1

1

1

J×J

1
− 1−γ
i1
1
− 1−γ
i2

..
.

1
− 1−γi,J
−1
1

0
.. 
. 


0 

− π1iJ 
1
J×J









J×J

where γij = πi1 + · · · + πij is the cumulative categorical probability, j = 1, . . . , J − 1. The
corresponding inverse matrices are
(CT D−1
L)−1
baseline
 i
2
−πi1 + πi1 −πi1 πi2
 −πi1 πi2 −π 2 + πi2
i2


.
..
.
=
.
.

−πi1 πi,J−1 −πi2 πi,J−1
−πi1 πiJ
−πi2 πiJ

△
= ci1 ci2 · · · ciJ baseline

···
···
..
.

−πi1 πi,J−1
−πi2 πi,J−1
..
.

2
· · · −πi,J−1
+ πi,J−1
···
−πi,J−1 πiJ

S19

πi1
πi2
..
.








πi,J−1 
πiJ
J×J

where (cij )baseline = πij (ej − π i ), j = 1, . . . , J − 1, (ciJ )baseline = π i , and ej is the J × 1
vector with the jth coordinate 1 and all others 0. Recall that π i = (πi1 , . . . , πiJ )T .
−1
(CT D−1
i L)cumulative

γi1 (1 − γi1 )
0

−γi1 (1 − γi1 ) γi2 (1 − γi2 )

=
0
−γi2 (1 − γi2 )


.
..
.

.
.
0
···

△
= ci1 ci2 · · · ciJ cumulative


···
0
πi1
..

..
.
.
πi2 
.. 
..
.
0
. 


..
. γi,J−1 (1 − γi,J−1) πi,J−1 
0 −γi,J−1 (1 − γi,J−1) πiJ
J×J

where (cij )cumulative = γij (1−γij )(ej −ej+1 ) with ej defined as above; and (ciJ )cumulative =
πi .
−1
(CT D−1
i L)continuation

πi1 (1 − γi1 )
0
···
..
 −π π
πi2 (1−γi2 )
.

i1 i2
1−γi1

.
.
..
..
..
=
.


πi2 πi,J −1
···
 −πi1 πi,J−1 − 1−γi1
πi2 πiJ
−πi1 πiJ
− 1−γi1
···

= ci1 ci2 · · · ciJ continuation

0
..
.
0
πi,J −1 (1−γi,J −1 )
1−γi,J −2
πi,J −1 πiJ
− 1−γ
i,J −2

πi1
πi2
..
.









πi,J−1 
πiJ

J×J

where (ci1 )continuation = πi1 (1 − γi1 , −πi2 , . . . , −πiJ )T ,
πij
(0, . . . , 0, 1 − γij , −πi,j+1 , . . . , −πiJ )T with “1 − γij ” being the
(cij )continuation = 1−γi,j−1
jth coordinate, j = 2, . . . , J − 1, and (ciJ )continuation = π i .
−1
(CT D−1
i L)adjacent

(1 − γi1 )πi1 (1 − γi2 )πi1
 −γi1 πi2
(1 − γi2 )πi2


.
..
..
=
.

 −γi1 πi,J−1 −γi2 πi,J−1
−γi1 πiJ
−γi2 πiJ

= ci1 ci2 · · · ciJ adjacent

···
···
..
.

(1 − γi,J−1 )πi1
(1 − γi,J−1 )πi2
..
.

· · · (1 − γi,J−1)πi,J−1
···
−γi,J−1 πiJ

πi1
πi2
..
.








πi,J−1 
πiJ
J×J

where (cij )adjacent = ((1−γij )πi1 , . . . , (1−γij )πij , −γij πi,j+1 , . . . , −γij πiJ )T , j = 1, . . . , J −
1, and (ciJ )adjacent = π i .
For certain applications, we need to know |CT D−1
i L| (see, for example, Lemma S.5).
Since adding to one row (column) any multiple of another row (column) does not change
the determinant (see, for example, 4.28(f) in Seber (2008, page 58)), we may (1) do row
operations on (CT D−1
i L)baseline and change it into an upper triangular matrix with diag−1
−1
−1
onal entries πi1 , . . . , πiJ
; (2) do row operations on (CT D−1
i L)cumulative and change it into
S20

an upper triangular matrix with diagonal entries γi1 (1 − γi1 ), . . . , γi,J−1(1 − γi,J−1 ), 1; (3)
do column operations on (CT D−1
i L)adjacent and change it into a lower triangular matrix
−1
−1
with diagonal entries πi1
, . . . , πiJ
; and (4) do column operations on (CT D−1
i L)continuation
−1
−1
and change it into a lower triangular matrix with diagonal entries πi1 , . . . , πiJ
. Therefore,

QJ
−1
for baseline-category,

j=1 πij


adjacent-categories,
(S.18)
|CT D−1
i L| =
and continuation-ratio logit models


 QJ−1 −1
−1
for cumulative logit models
j=1 γij (1 − γij )
Apparently, |CT D−1
i L| > 0 as long as πij > 0 for all j = 1, . . . , J.

S.5

Formulas for calculating πij ’s from Xi’s

Following the notations in model (1), η i = Xi θ = CT log(Lπ i ). The formulas towards
calculating πij ’s as follows:
(1) Baseline-category logit model





1
0
1
−1

 1
 1
0
−1








.. 
.. 
.
..
.
·
exp(η
)
·
log
log(π i ) = 

.
.
.
. 
i 








1 0
1 −1
1 1 · · · 1 1 J×J
−1 J×J

(2) Adjacent-categories logit model

log(π i ) =




−1
 1
−1



.. 
..
·

.
. 




1 −1
−1 J×J



1
1
0


 1
0





.. 
..
·
exp
log 


.
.





1 0
1 1 · · · 1 1 J×J
1

S21

1 ···
1 ···
..
.



1 0

1 0




.. ..
·
η


i
. .



1 0
1 J×J

(3) Continuation-ratio logit model

1
1


log(π i ) = η i −  ...

1
1

1
0


log 

0
0

1
.
1 ··· 1
1 ··· 1 1

(4) Cumulative logit model


1
1


log  ...

1
0



1
..

.
1 ···
0 ···

1
0 1








J×J

..









·

J×J


0 ··· 0 1

1 · · · 0 1




..
·
exp(η
)

.
i 



0 · · · 1 1
0 · · · 0 1 J×J




1

 
 
· πi = 
 

0

..

−1
.

···


1





· log 
1





1
0

..

.

..

.

0

..

.

···






−1
−1 J×2(J−1)



0

.. 

.



1 0


·
exp(η
)

i
1



.. 

.

1 1 2(J−1)×J

Note that Xi θ in the above models could be po, npo, or ppo.

S.6

Algorithms

Lift-one algorithm for D-optimal allocation w = (w1 , . . . , wm )T :
1◦ Start with an arbitrary allocation w0 = (w1 , . . . , wm )T satisfying 0 < wi < 1,
i = 1, . . . , m and compute f (w0 ).
2◦ Set up a random order of i going through {1, 2, . . . , m}.
3◦ For each i, determine fi (z) according to Theorem 5.2. In this step, J determinants
fi (0), fi (1/2), fi (1/3), . . . , fi (1/J) are calculated.
4◦ Use quasi-Newton algorithm to find z∗ maximizing fi (z) with 0 ≤ z ≤ 1. If
(i)
fi (z∗ ) ≤ fi (0), let z∗ = 0. Define w∗ = (w1 (1 − z∗ )/(1 − wi ), . . . , wi−1 (1 − z∗ )/(1 −
(i)
wi ), z∗ , wi+1 (1 − z∗ )/(1 − wi ), . . . , wm (1 − z∗ )/(1 − wi ))T . Note that f (w∗ ) = fi (z∗ ).
(i)

(i)

5◦ Replace w0 with w∗ , and f (w0 ) with f (w∗ ).
S22

(i)

6◦ Repeat 2◦ ∼ 5◦ until convergence, that is, f (w0 ) = f (w∗ ) for each i.
Exchange algorithm for D-optimal allocation (n1 , . . . , nm )T given n > 0:
1◦ Start with an initial allocation n = (n1 , . . . , nm )T such that f (n) > 0.
2◦ Set up a random order of (i, j) going through all pairs {(1, 2), (1, 3), . . . , (1, m),
(2, 3), . . . , (m − 1, m)}.
3◦ For each (i, j), let c = ni + nj . If c = 0, let n∗ij = n. Otherwise, there are
two cases. Case one: 0 < c ≤ q, we calculate fij (z) for z = 0, 1, . . . , c directly
and find z ∗ which maximizes fij (z). Case two: c > q, we first calculate fij (z) for
z = 0, 1, . . . , q; secondly determine g0 , g1 , . . . , gq in (18) according to Theorem 5.3;
thirdly calculate fij (z) for z = q+1, . . . , c based on (18); fourthly find z ∗ maximizing
fij (z) for z = 0, . . . , c. For both cases, we define
n∗ij = (n1 , . . . , ni−1 , z ∗ , ni+1 , . . . , nj−1 , c − z ∗ , nj+1 , . . . , nm )T
Note that f (n∗ij ) = fij (z ∗ ) ≥ f (n) > 0. If f (n∗ij ) > f (n), replace n with n∗ij , and
f (n) with f (n∗ij ).
4◦ Repeat 2◦ ∼ 3◦ until convergence, that is, f (n∗ij ) = f (n) in step 3◦ for any (i, j).

S.7

Model selection

See Tables 4 and 5.
Table 4: Model Comparison for the Trauma Data

AIC
BIC

Cumulative
po

Cumulative
npo

Continuation
po

Continuation
npo

Adjacent
po

Adjacent
npo

107.75
104.68

99.41
94.51

108.98
105.91

101.36
96.45

107.67
104.60

101.54
96.63

Table 5: Model Comparison for House Flies Data (po and npo)

AIC
BIC

Cumulative
po

Cumulative
npo

Continuation
po

Continuation
npo

Adjacent
po

Adjacent
npo

195.87
195.71

121.17
120.96

116.40
116.24

114.42
114.20

209.64
209.47

194.47
194.25

S23

S.8

More examples

Example S.1. (For Section 3.3) Consider an experiment with three factors (d = 3), three
response categories (J = 3), and three distinct experimental settings (m = 3). Denote
the experimental settings as xi = (xi1 , xi2 , xi3 )T , i = 1, 2, 3. Consider a multinomial
logistic model with ppo such that


 


1 x11
1
x12 x13
HT1 =  1 x21  , HT2 =  1  , HTc =  x22 x23 
1 x31
1
x32 x33
That is, p1 = 2, p2 = 1, pc = 2, pH = 1,
p1 + p2 + pc = 5 parameters. In this case,

1
1
 x11 x21

0
H=
 0
 x12 x22
x13 x23

max{p1 , p2 , pc + pH } = 3, and there are p =
1
0
0
0
x31 0
0
0
0
1
1
1
x32 x12 x22 x32
x33 x13 x23 x33








is 5×6. It can be verified that rank(H) = 5 for general x1 , x2 , x3 using Theorem 3.5. That
is, the minimal number of experimental settings in this case is m = max{p1 , . . . , pJ−1 , pc +
pH } = 3.

Example S.2. (For Section 4) Consider an example with responses in J = 4 categories,
d = 5 factors, and m = 5 distinct experimental settings xi = (xi,1 , . . . , xi,5 )T , i = 1, . . . , 5.
Suppose a multinomial logistic model with




HT1 = 


1 x11
1 x21
..
..
.
.
1 x51

x12
x22
..
.
x52









 , HT2 = 



1
1
..
.

x11
x21
..
.

1

x51









 , HT3 = 



1
1
..
.
1









 , HTc = 



x13
x23
..
.

x14
x24
..
.

x15
x25
..
.

x53

x54

x55







is used. That is, p1 = 3, p2 = 2, p3 = 1, pH = 1, pc = 3, and p = 9. In this
case, G defined in Theorem 4.1 is 20 × 9 and pc + pH = 4 is the minimum number
of #{i | αi > 0} to keep |G[i1 , . . . , ip ]| =
6 0 if (i1 , . . . , ip ) ∈ Λ(α1 , . . . , αm ). Actually,
(i1 , . . . , i9 ) = (1, 2, 3, 6, 7, 8, 10, 11, 12) ∈ Λ(3, 3, 3, 0, 0) leads to rank(G[i1 , . . . , i9 ]) = 8,
while (1, 2, 5, 6, 9, 10, 13, 14, 15) ∈ Λ(2, 2, 2, 3, 0) leads to rank(G[i1 , . . . , i9 ]) = 9. Therefore, |G[i1 , . . . , i9 ]| =
6 0 in general if (i1 , . . . , i9 ) ∈ Λ(2, 2, 2, 3, 0) for such a ppo model.

Example S.3. (For Section 5.5) Consider a multinomial logistic model with proportional odds for responses with J = 3 categories, d = 1 factors, and m = 2 distinct
experimental settings x1 , x2 . Same as Example 2.3, the parameters are β1 , β2 , ζ1 and the
linear predictors
ηi1 = β1 + ζ1 xi , ηi2 = β2 + ζ1 xi , i = 1, 2.
According to Theorem 4.2, the objective function of allocation (w1 , w2 ) is an order-3
homogeneous polynomial of w1 , w2 consisting of monomials cα1 ,α2 w1α1 w2α2 with coefficients
S24

cα1 ,α2 ≥ 0. Based on Lemma 4.1 and Corollary 4.2, cα1 ,α2 > 0 only if max{α1 , α2 } ≤ 2
and #{i | αi > 0} = 2, which implies (α1 , α2 ) is either (2, 1) or (1, 2). That is, the
objective function is
f (w1 , w2 ) = w1 w2 (c21 w1 + c12 w2 ),
which takes the same form as in Corollary 5.2 in Yang et al. (2017). If we rewrite c21 =
C · c2 and c12 = C · c1 , that is, f (w1 , w2) = C · w1 w2 (c2 w1 + c1 w2 ), then for a baselinecategory logit model, C = π13 π23 (x1 − x2 )2 , c2 = π11 π12 (1 − π23 ), c1 = π21 π22 (1 − π13 );
−1
−1
for a cumulative logit model, C = π12
(1 − π13 )(1 − π11 )π22
(1 − π23 )(1 − π21 )(x1 − x2 )2 ,
c2 = π11 (1 − π11 )π13 (1 − π13 )π22 (1 − π22 ), c1 = π12 (1 − π12 )π21 (1 − π21 )π23 (1 − π23 ); for
an adjacent-categories logit model, C = (x1 − x2 )2 , c2 = π11 π12 π13 (π21 π22 + π22 π23 +
4π21 π23 ), c1 = π21 π22 π23 (π11 π12 + π12 π13 + 4π11 π13 ); for a continuation-ratio logit model,
C = (1 − π11 )−1 (1 − π21 )−1 (x1 − x2 )2 , c2 = π11 π12 π13 (1 − π11 )[π22 π23 + π21 (1 − π21 )2 ],
c1 = π21 π22 π23 (1 − π21 )[π12 π13 + π11 (1 − π11 )2 ]. According to Corollary 5.2 in Yang et al.
(2017), the uniform allocation w1∗ = w2∗ = 1/2 is D-optimal if and only if c1 = c2 , which
is not true in general for all the four logit models.


S25

