ARTICLE IN PRESS
Statistics & Probability Letters 77 (2007) 426–430 www.elsevier.com/locate/stapro
Optimal designs in multivariate linear models A. Markiewicz, A. Szczepan´ska Department of Mathematical and Statistical Methods, Agricultural University of Poznan´, Wojska Polskiego 28, PL-60637Poznan´, Poland Received 10 June 2005; received in revised form 27 June 2006; accepted 7 August 2006 Available online 5 September 2006
Abstract The purpose of this paper is to study optimality of an experimental design under the multivariate models with a known or unknown dispersion matrix. In the case of unknown dispersion matrix optimality is considered with respect to the precision in maximum likelihood estimation. We show relations between optimality of designs in univariate models and in their multivariate extensions. r 2006 Elsevier B.V. All rights reserved. MSC: primary 62K05; secondary 62K10 Keywords: Multivariate model; Growth curve model; Estimable function; Optimal design; Universal optimality
1. Introduction Consider a linear model associated with a design d 2 D y ¼ A1;d b1 þ A2 b2 þ e;
DðeÞ ¼ V,
(1)
where A1;d 2 Rnv , depending on d, and A2 2 Rnb , the same for every design, are known, while e 2 Rn is an n-dimensional vector of random errors with EðeÞ ¼ 0 and the matrix V 2 R4 n is a known n n positive definite dispersion matrix. The vector y 2 Rn is an observable random vector, which depends linearly on several parameters. The vector y represents measurements on a single response variable y. Here, b1 2 Rv is the vector of parameters of interest and b2 2 Rb is the vector of nuisance parameters. If we are interested in measuring q response variables instead of one variable on each of sampling units we can consider a multivariate linear model associated with a design d 2 D: Y ¼ A1;d B1 C1 þ A2 B2 C2 þ E.
(2)
Matrices A1;d and A2 are the same as in model (1) and matrices C1 2 Rpq and C2 2 Rsq are known. The matrix Y 2 Rnq is an observations matrix, while B1 2 Rvp is the matrix of parameters of interest and B2 2 Rbs is the matrix of nuisance parameters. The matrix E 2 Rnq is a matrix of random errors with expectation EðEÞ ¼ 0 and with dispersion matrix DðvecðEÞÞ ¼ R V, where vecðAÞ denotes the vector from Corresponding author.
E-mail address:
[email protected] (A. Markiewicz). 0167-7152/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.spl.2006.08.010
ARTICLE IN PRESS A. Markiewicz, A. Szczepan´ska / Statistics & Probability Letters 77 (2007) 426–430
427
Rmn formed by writing the columns of A 2 Rnm one under the other in the sequence and denotes the Kronecker product. The matrix V 2 R4 n is a known n n positive definite dispersion matrix of columns of E while the matrix R 2 R4 is a known or unknown dispersion matrix of rows of E. q Model (2) with restrictions or when no restrictions exist was studied by Kollo and von Rosen (2005). In the literature model (2) is usually considered with an additional condition that RðC01 Þ RðC02 Þ, i.e. the model: Y ¼ A1;d B1 C1 þ A2 B2 C2 þ E;
RðC01 Þ RðC02 Þ.
(3)
Chinchilli and Elswick (1985) studied model (3) reduced to the MANOVA-GMANOVA, i.e. the model in which C2 ¼ Iq , where the matrix Iq is the q q identity matrix. Fujikoshi and Satoh (1996) considered model (3) as an extended growth curve model in which a polynomial curve of degree p 1 is to be fit to q responses over time. Moreover, this type of the growth curve model was studied by Yokoyama and Fujikoshi (1992) with C2 ¼ Iq and by Srivastava and Khatri (1979) with C2 ¼ C1 . In this paper we are interested in the optimality of designs w.r.t. the estimation of B1 in model (3). Our purpose is to show that designs which are universally optimal in the Kiefer sense (see Kiefer, 1975) in univariate model (1) are also optimal in a certain sense in multivariate model (3) when the dispersion matrix R is known as well as when it is unknown. In Section 2 we recall the criteria for estimability in models (1)–(3). In further part of the paper we study designs optimality assuming that the matrix E of random errors is normally distributed (vec ENð0; R VÞ). In Section 3 we consider the Kiefer optimality w.r.t. the estimation of B1 in model (3) with known dispersion matrix R. The case of R unknown is considered in Section 4 where we define and study the Kiefer optimality w.r.t. the precision in the maximum likelihood estimation. 2. Estimability We are interested in optimality of designs for the estimation of B1 in model (3), i.e. the estimation of all estimable functions of the form: KB1 L, where K 2 Rgv and L 2 Rph . First, we derive a general condition for estimability of parametric functions in model (2). Observe that using the vec operator and its property: vecðABCÞ ¼ ðC0 AÞ vec B, model (2) can be represented in an univariate manner as vec Y ¼ ðC01 A1;d Þ vec B1 þ ðC02 A2 Þ vec B2 þ vec E. 0
(4)
0
Let PL ¼ LðL LÞ L and QL ¼ Im PL denote the orthogonal projectors on RðLÞ and the orthocomplement of RðLÞ, respectively, where RðLÞ stands for the range of a given matrix L 2 Rmn . Following Baksalary (1984), a vector of parametric functions ðL0 KÞvecB1 in model (4) is estimable iff RðL K0 Þ RððC1 A01;d ÞQC 02 A2 Þ. Using the above estimability criterion we can formulate the following theorem. Theorem 1. Let K 2 Rgv and L 2 Rph . A matrix of linear parametric functions KB1 L is estimable in model (2) if and only if RðL K0 Þ RðC1 QC 02 A01;d PA2 þ C1 A01;d QA2 Þ.
(5)
Note that estimability criterion presented in Theorem 1 can be simplified in model (3). Corollary. Let K 2 Rgv and L 2 Rph . A matrix of parametric functions KB1 L is estimable in model (3) if and only if RðK0 Þ RðA01 QA2 Þ
and
RðLÞ RðC1 Þ.
(6)
3. Optimality in models with the known dispersion matrix In this section we consider optimality of designs for the estimation of B1 in model (3), where the matrix of random errors is normally distributed with mean 0 and the known dispersion matrix R V (vec ENð0; R VÞ). To derive optimal designs in a multivariate model with the known dispersion matrix we apply an optimality criterion from a univariate model on univariate formulation (4) of model (3) and determine the information matrix for the estimation of B1 in model (3) using the following lemma.
ARTICLE IN PRESS A. Markiewicz, A. Szczepan´ska / Statistics & Probability Letters 77 (2007) 426–430
428
Lemma 1. Let R 2 Rnv , S 2 Rnb , T 2 Rpq , U 2 Rsq and RðT0 Þ RðU0 Þ, then QU 0 S ðT0 RÞ ¼ T0 QS R.
(7)
Using the vec operator we can present model (2) in a univariate manner as model (4) with a known dispersion matrix. Following Markiewicz (2001), the information matrix in model (4) can be written: Nd ¼ ðC1 R1=2 A01 V1=2 ÞQðS1=2 C 0 ÞðV 1=2 A2 Þ ðR1=2 C01 V1=2 A1 Þ. 2
(8)
Using Lemma 1 the above information matrix (8) assumes simpler form in model (3) as Nd ¼ ðC1 R1 C01 Þ ðA01;d QA2 ðQA2 VQA2 Þ QA2 A1;d Þ ¼ F Cd ,
(9)
where Cd is the information matrix for d design in model (1) (see Markiewicz, 2001 and references given there). A g-inverse of the information matrix is a dispersion matrix of a generalized least-squares estimator of B1 and can be written as b 1 Þ ¼ ðC1 R1 C0 Þ ðA0 QA ðQA VQA Þ QA A1;d Þ ¼ F C , DðB 1 1;d d 2 2 2 2
(10)
b 1 is a generalized least-squares estimator of B1 taking the form where B b 1 ¼ ðA0 QA ðQA VQA Þ QA A1;d Þ A0 QA ðQA VQA Þ QA YR1 C0 ðC1 R1 C0 Þ . B 1;d 1;d 1 1 2 2 2 2 2 2 2 2
(11)
Note that for a matrix of estimable functions KB1 L conditions (6) are fulfilled and its generalized least-squares c1 L and its dispersion matrix are invariant w.r.t. a choice of g-inverses involved in (10) and (11). estimator KB In the paper we consider the Kiefer optimality relative to the group H which is based on the Kiefer ordering defined as follows. Given two symmetric matrices A and B 2 Sym(k) and a group H OrtðkÞ, we say that A is below B w.r.t. the Kiefer ordering relative to H (B is more informative than A) and we write A5H B when B is better in the Loewner ordering than some matrix D which is H-majorized by A, i.e. A5H B () DL B
for some D 2 SymðkÞ such that D H A,
where D H A means that D 2 conv½HAH0 : H 2 H while DL B means that B D 2 NNDðkÞ. Definition 1. A design d is called Kiefer optimal relative to the group H when its information matrix Zd is H- invariant and Zd 5H Zd
for all d 2 D.
Note that the Kiefer optimality implies the universal optimality; see Pukelsheim (1993, Chapter 14). Let design d be Kiefer optimal relative to the group Iq H in model (3). Then from (9) we have the following equivalences: Nd ¼ F Cd is Iq H-invariant()Cd is H-invariant and Nd 5I q H Nd ()Cd 5H Cd . From the above conditions we get the following relation between optimality of a design in univariate model (1) and its multivariate extension (3). Theorem 2. The following statements are equivalent: (a) a design d is Kiefer optimal relative to the group H for the estimation of b1 in model (1); (b) a design d is Kiefer optimal relative to the group Iq H for the estimation of B1 in model (3). To illustrate applicability of Theorem 2, consider the following example. Let D be the class of block designs with this same number k of plots in each of b blocks on which vXk treatments are allocated, and let
ARTICLE IN PRESS A. Markiewicz, A. Szczepan´ska / Statistics & Probability Letters 77 (2007) 426–430
429
y ¼ A1;d b1 þ ðIb 1k Þb2 þ e, CovðeÞ ¼ I be the associated model. Suppose there exists a design d 2 D which information matrix for the estimation of b1 , Cd , is Pv -invariant (completely symmetric), where Pv is the group of v v permutation matrices. Moreover, let tr Cd ptr Cd for all d 2 D, which implies that Cd 5Pv Cd for all d 2 D, i.e. that d is Kiefer optimal relative to the group Pv for the estimation of b1 . This property is fulfilled by BIB designs. According to Theorem 2, a BIB design is also Kiefer optimal relative to the group Iq Pv for the estimation of B1 in model (3) with A2 ¼ ðIb 1k Þ, V ¼ In , and R known. 4. Optimality in models with the unknown dispersion matrix Usually, the dispersion matrix R in model (3) is unknown. In such case the best linear unbiased estimator of estimable functions of B1 does not exist and the maximum likelihood estimation is derived. The maximum likelihood estimator of B1 in model (3) under normality takes the form e 1 ¼ ðA0 QA ðQA VQA Þ QA A1;d Þ A0 QA ðQA VQA Þ QA YS1 C0 ðC1 S1 C0 Þ , B 1;d 1;d 1 1 2 2 2 2 2 2 2 2
(12)
where S ¼ Y0 QA2 QQA A1;d QA2 Y (see Kollo and von Rosen, 2005; for more details). 2 e 1 is The unconditional covariance matrix of B e 1 Þ ¼ ðC0 R1 C1 Þþ ðA0 QA ðQA VQA Þ QA A1;d Þþ ½md 1 =½md q þ rðC1 Þ 1
DðB 1 1;d 2 2 2 2 ¼ Pþ Cþ d ½md 1 =½md q þ rðC1 Þ 1
ð13Þ
(see Grizzle and Allen, 1980), where Cd is the information matrix for d design in model (1), . md ¼ rðA1;d Þ þ rðA2 Þ rðA1;d ..A2 Þ, and rð:Þ denotes the rank of a matrix argument. We consider optimality w.r.t. the estimation of all estimable functions of the form: KB1 L, where K 2 Rgv and L 2 Rph . Similarly as in model (3) with R known, the maximum likelihood estimator of a matrix of f1 L and its dispersion matrix are invariant w.r.t. the choice of g-inverses involved in (12) estimable functions KB and (13), respectively. The Moore–Penrose inverse of the dispersion matrix (13) is the precision matrix in maximum likelihood estimation, which depends on the unknown matrix R, and is equal e 1 Þ þ ¼ ðC0 R1 C1 Þ ðA0 QA ðQA VQA Þ QA A1;d Þ½md q þ rðC1 Þ 1 =½md 1
Wd ¼ ½DðB 1 1;d 2 2 2 2 ¼ P Cd ½md q þ rðC1 Þ 1 =½md 1 .
ð14Þ
When R is unknown we may consider the Kiefer optimality relative to the group H w.r.t. the precision in the maximum likelihood estimation defined as follows. Definition 2. A design d is called the Kiefer optimal relative to the group H w.r.t. the precision in the maximum likelihood estimation when the precision matrix Wd is H-invariant and Wd 5H Wd for all d 2 D. Noting that precision matrix (14) has a similar form as information matrix (9) in model (3) with R known, we get a counterpart of Theorem 2. Theorem 3. The following statements are equivalent: (a) a design d is the Kiefer optimal relative to the group H for the estimation of b1 in model (1); (b) a design d is the Kiefer optimal relative to the group Iq H w.r.t. the precision in the maximum likelihood estimation of B1 in model (3). Observe that similarly as in the case of R known (Section 3), the following relation holds, Wd 5I q H Wd ()Cd 5H Cd .
ARTICLE IN PRESS A. Markiewicz, A. Szczepan´ska / Statistics & Probability Letters 77 (2007) 426–430
430
In consequence, a BIB design is also Kiefer optimal relative to the group Iq Pt w.r.t. the precision in the maximum likelihood estimation of B1 in the model Y ¼ A1;d B1 C1 þ ðIb 1k ÞB2 C2 þ E;
DðvecðEÞÞ ¼ R I
with R unknown.
References Baksalary, J.K., 1984. A study of the equivalence between a Gauss–Markoff model and its argumentation by nuisance parameters. Math. Operationsforsch. Statist. 15, 3–35. Chinchilli, V.M., Elswick, R.K., 1985. A mixture of the MANOVA and GMANOVA models. Comm. Statist. Theory Methods 14 (12), 3075–3089. Grizzle, J.E., Allen, D.M., 1980. Analysis of growth and dose response curves. Biometrics 25, 357–381. Fujikoshi, Y., Satoh, K., 1996. Estimation and model selection in the extended growth curve model. Hiroshima Math. J. 26, 635–647. Kiefer, J., 1975. Construction and optimality of generalized Youden designs. In: Srivastava, J.N. (Ed.), A Survey of Statistical Design and Linear Models. North-Holland, Amsterdam, pp. 333–353. Kollo, T., von Rosen, D., 2005. Advanced Multivariate Statistics with Matrices. Springer, Dordrecht. Markiewicz, A., 2001. On dependence structures preserving optimality. Statist. Probab. Lett. 53, 415–419. Pukelsheim, F., 1993. Optimal Design of Experiments. Wiley, New York. Srivastava, M.S., Khatri, C.G., 1979. An Introduction to Multivariate Statistics. North-Holland, New York. Yokoyama, T., Fujikoshi, Y., 1992. Tests for random-effects covariance structures in the growth curve model with covariates. Hiroshima Math. J. 22, 195–202.
