The Growth Curve Model and its Extension ^{1}

1. Models Consider the general form of the multivariate normal linear model which is summarized in the following:

Y : n × p = ∑_{i=1}^{m} A_{i}β_{i}C_{i} + ϵ, with ϵ ∼ N(0, I_{n} ⊗ Σ)  (1)

where Y is an observation matrix, and A_{i} and C_{i} are known n × a_{i} and c_{i} × p matrices of the full rank, and β_{i} are unknown a_{i} × c_{i} matrices, i = 1, 2, …, m. We also assume that the rows of ϵ are independently p-variate normally distributed with a covariance matrix Σ. Then following are examples of (??).

1. Multivariate regression model or MANOVA: Let m = 1 and C_{1} = I. Then

𝔼[Y] = Aβ.

2. The growth curve model or GMANOVA ( Potthoff and Roy (1964)): ^{2} Let m = 1.

Then

𝔼[Y] = AβC.

3. Mixture of MANOVA and GMANOVA (Chinchilli and Elswick (1985)): ^{3} Let m = 2

and C_{1} = I. Then

𝔼[Y] = A_{1}β_{1} + A_{2}β_{2}C_{2},

where [A_{1}, A_{2}] is of full rank.

4. Extended GMANOVA model with the nested subspace condition (von Rosen (1989)): 4 Assume that the nested subspace condition, i.e., span(A_{1}) ⊇ span(A_{2}) ⊇ ⋯ ⊇ span(A_{m}), is hold, where span(·) denotes the column vector space.

1Statistics 542 (March 11, 1998)

2Consider a groups of animals, with n_{j} animals in the jth group, and with each group being subjected to a different treatment. Animals in all groups are measured at the same p points in time t_{1}, t_{2}, …, t_{p} and are assumed to have the same covariance matrix Σ. The growth curve associated with the jth group is β_{j0} + β_{j1}t + β_{j2}t^{2} + ⋯ + β_{jc}t^{c}. The matrix A will contain c columns and will consist of n_{1} rows (1, 0, …, 0), n_{2} rows (0, 1, 0, …, 0), …, n_{a} rows (0, …, 0, 1) and C_{ij} = t_{i}^{j}, 1 ≤ i ≤ c, 1 ≤ j ≤ p).

3The response of interest is a score on a psycho-motor testing device.There were four different groups: control ( 6 patients), low dose (14 patients), medium dose (15 patients), and high dose (10 patients). A baseline measurement was taken at day 0 (pre-treatment), and then measurements were taken on a daily basis for ten consecutive days (post-treatment). Then dimensions of Y, A_{1}, β_{1}, A_{2}, β_{2}, C_{2} are 45 × 10, 45 × 1, 1 × 10, 45 × 4, 4 × 2, 2 × 10.

4Wesolowska-Janczarek (1996) considered the model of the form Y = A_{1}β_{1}C_{2}+1_{n}β_{2}C_{2}+ϵ where Y : n×p is an observation matrix, obtained on n experimental units in p time points, A_{1} is a design matrix which divides the observations into a_{1} groups, and C_{2} is a matrix of concomitant variables.

5. Sum of profiles with changing covariates (Patel(1986)) ^{5}: Let C_{i} = I and β_{i}(i ≥ 2) are diagonal matrices of unknown parameters. Then

𝔼[Y] = A_{1}β_{1} + ∑_{i=2}^{m} A_{i}β_{i}.

2. A canonical reduction of the model \underline{MANOVA}. Let H be an (n − a) × n matrix such that

HA = 0 and HH′ = I_{n−a},

so that the rows of H form an orthogonal basis for {span(A)}^{⊥}, the orthogonal complement of the range of A. Hence H′H = I − A(A′A)^{−1}A′. Now, the distribution of the matrix

[ \begin{matrix}X \\ Z \end{matrix} ] = [ \begin{matrix}(A′A)^{−1/2}A′ \\ H \end{matrix} ] Y

is normal with the mean

𝔼 [ \begin{matrix}X \\ Z \end{matrix} ] = [ \begin{matrix}(A′A)^{−1/2}A′Aβ \\ HAβ \end{matrix} ] = [ 0^{(A′A)^{1/2}β} ] ,

and the covariance matrix

([ \begin{matrix}(A′A)^{−1}A′ \\ H \end{matrix} ] ⊗ I_{n}) (I ⊗ Σ) ([ A(A′A)^{−1/2} ⋮ H′ ] ⊗ I_{n}) = I ⊗ Σ.

The maximum likelihood estimators of β and Σ are given by

β̂ = (A′A)^{−1}A′Y

and

nΣ̂ = Z′Z = Y′H′HY = Y′(I − A(A′A)^{−1}A′)Y = (Y − Aβ̂)′(Y − Aβ̂).

Note that (β̂, Σ̂) is sufficient for (β, Σ) (also complete ) and Z′Z ∼ W_{p}(n − a, Σ).

\underline{GMANOVA}. Let D : p × (p − c) be a matrix such that

CD = 0 and D′D = I_{p−c}

5In Patel’s model, A_{1} is a design matrix and A_{i}(i ≥ 2) are matrices of values of covariates. Patel proposed an iterative scheme for the estimation of β_{1} both in the unconstrained case and subject to the general linear hypothesis D_{1}β_{1}E_{1} = 0. Later Verblya(1986) pointed out that this model can be written as seemingly unrelated regressions

1Since [ \begin{matrix}C \\ D′ \end{matrix} ] [ C′(CC′)^{−1} ⋮ D ] = I_{p}, the transpose of this equation gives C′(CC′)^{−1}C + DD′ = I.

Then the distribution of the matrix

[ \begin{matrix}X_{1} ⋮ & X_{2} \\ . & . & . & . & . & . & . & . \\ Z \end{matrix} ] = [ \begin{matrix}(A′A)^{−1/2}A′ \\ H \end{matrix} ] Y [ C′(CC′)^{−1} ⋮ D ]

is normal with the mean

𝔼 [ \begin{matrix}X_{1} ⋮ & X_{2} \\ . & . & . & . & . & . & . & . \\ Z \end{matrix} ] = [ \begin{matrix}(A′A)^{1/2}β ⋮ & & & 0 \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 0 \end{matrix} ]

and the covariance matrix

([ \begin{matrix}(A′A)^{−1/2}A′ \\ H \end{matrix} ] ⊗ [ \begin{matrix}(CC′)^{−1}C \\ D′ \end{matrix} ]) (I_{n} ⊗ Σ) ([ A(A′A)^{−1/2} ⋮H′ ] ⊗ [ C′(CC′)^{−1} ⋮D ])

= ([ \begin{matrix}(A′A)^{−1/2}A′ \\ H \end{matrix} ] [ A(A′A)^{−1/2} ⋮H′ ]) ⊗ ([ \begin{matrix}(CC′)^{−1}C \\ D′ \end{matrix} ] Σ [ C′(CC′)^{−1} ⋮D ])

= I_{n} ⊗ [ \begin{matrix}(CC′)^{−1}CΣC′(CC′)^{−1} & (CC′)^{−1}CΣD \\ D′ΣC′(CC′)^{−1} & D′ΣD \end{matrix} ] = I_{n} ⊗ Σ̃ ( say ).

Partition Σ̃ as

Σ̃ = [ \begin{matrix}Σ̃_{11} & Σ̃_{12} \\ Σ̃_{21} & Σ̃_{22} \end{matrix} ]

where Σ̃_{11} is c × c matrix and put μ = (A′A)^{1/2}β. Recall the relation ^{6}

Σ̃ = [ 0 I_{p−c}^{I_{c} Σ̃_{12}Σ̃_{22}^{−1}} ] [ \begin{matrix}Σ̃_{11.2} & 0 \\ 0 & Σ̃_{22} \end{matrix} ] [ I_{c} 0_{Σ̃_{22}^{−1}Σ̃_{21} I_{p−c}} ] ,

where Σ̃_{11.2} = Σ̃_{11} − Σ̃_{12}Σ̃_{22}^{−1}Σ̃_{21}. Then

tr ([ \begin{matrix}X_{1} − μ ⋮ & & X_{2} \\ . & . & . & . & . & . & . & . & . & . & . & . & . \\ Z \end{matrix} ] Σ̃^{−1} [⋯]′)

= tr ([ X_{1} − μ ⋮ X_{2} ] Σ̃^{−1} [ \begin{matrix}(X_{1} − μ)′ \\ X′_{2} \end{matrix} ]) + tr(Z′ZΣ̃^{−1})

= tr ({X_{1} − X_{2}Σ̃_{22}^{−1}Σ̃_{21} − μ}Σ̃_{11.2}^{−1}{⋯}′) + tr(S̃_{11.2}Σ̃_{11.2}^{−1})

+tr((S̃_{22} + X′_{2}X_{2})Σ̃_{22}^{−1}) + tr(S̃_{22}{S̃_{22}^{−1}S̃_{21} − Σ̃_{22}^{−1}Σ̃_{21}}Σ̃_{11.2}^{−1}{⋯}′),  (2)
Σ̃^{−1} = [ I_{c}  0_{−Σ̃_{22}^{−1}Σ̃_{21} I_{p−c}} ] [ \begin{matrix}Σ̃_{11.2}^{−1} & 0 \\ 0 & Σ̃_{22}^{−1} \end{matrix} ] [ 0 I_{p−c}^{I_{c} Σ̃_{12}Σ̃_{22}^{−1}} ] .

6Inverse formula is given by

where S̃ = Z′Z and the partition of S̃ is defined like Σ̃. Hence the log-likelihood function is

−\frac{1}{2}\{ RHS of (??) \} − \frac{n}{2} log det(Σ̃_{11.2}) − \frac{n}{2} log det(Σ̃_{22}),

where parameters are now (μ, Σ̃_{11.2}, Σ̃_{22}^{−1}Σ̃_{21}, Σ̃_{22}) = (μ, Λ_{1}, R, Λ_{2}) ^{7} ( say ). The maximum likelihood estimators of the parameters are given by

μ̂ = X_{1} − X_{2}S̃_{22}^{−1}S̃_{21}
Λ̂_{1} = \frac{1}{n}S̃_{11.2} R̂ = S̃_{22}^{−1}S̃_{21} Λ̂_{2} = \frac{1}{n}(S̃_{22} + X′_{2}X_{2}),  (3)

where

S̃_{11.2} ∼ W_{c}(n−a−p+c, Σ̃_{11.2}), S̃_{22}+X′_{2}X_{2} ∼ W_{p−c}(n, Σ̃), S̃_{22}^{−1}S̃_{21}|S̃_{22} ∼ N(Σ̃_{22}^{−1}Σ̃_{21}, S̃_{22}^{−1}⊗Σ̃_{11.2}).

Furthermore, simple reconstruction gives

nˆΣ̃ = [ 0 I_{p−c}^{I_{c} S̃_{12}S̃_{22}^{−1}} ] [ \begin{matrix}S̃_{11.2} & 0 \\ 0 & S̃_{22} + X′_{2}X_{2} \end{matrix} ] [ I_{c} 0_{S̃_{22}^{−1}S̃_{21} I_{p−c}} ] .

By the virtue of the relation

Σ^{−1} − Σ^{−1}C′(CΣ^{−1}C′)^{−1}CΣ^{−1} = D(D′ΣD)^{−1}D′ for Σ > 0,  (4)

and that the estimator for Σ is

Σ̂ = [ C′ D ] ˆ Σ̃ [ \begin{matrix}C \\ D′ \end{matrix} ] = S + SD(D′SD)^{−1}D′Y′A(A′A)^{−1}A′Y D(D′SD)^{−1}D′S,

[ C′ ⋮D ] [ 0 I_{p−c}^{I_{c} S̃_{12}S̃_{22}^{−1}} ] = [ C′ ⋮SD(D′SD)^{−1} ] .

7Note that (β, Σ) is uniquely determined by these parameters.

10Since Σ is positive definite, there exists a nonsingular matrix Q : p × p such that Σ = QQ′. Similarly, (C′Σ^{−1}C)^{−1} = EE′ and (D′ΣD)^{−1} = F F′, where E : c × c and F : (p − c) × (p − c). Now put Δ : n × n = [ Q^{−1}CE ⋮ QDF ]. Then Δ′Δ = I. Hence,

I = ΔΔ′ = Q^{−1}C(EE′)C′Q^{′−1} + QD(F F′)D′Q′.

11By (??) and the definition of S and S̃,

S̃_{11.2} = (CS^{−1}C′)^{−1}, S̃_{22} = D′SD, X′_{2}X_{2} = D′Y′A(A′A)^{−1}A′Y D.

From the identity I_{p} − C′(CC′)^{−1}C = DD′, we may observe that C′S̃_{12}S̃_{22}^{−1} + D = SD(D′SD)^{−1}, from which it follows that

where S = Y′[I − A(A′A)^{−1}A′]Y. Again apply (??) to LHS of equation above to see that

Σ̂ = \frac{1}{n}(Y − Aβ̂C)′(Y − Aβ̂C).

Note that (β̂, Σ̂) is sufficient for (β, Σ) but it is not complete. So no Uniformally Minimum Variance Unbiased estimator exists. See von Rosen (1991) for an example of an unbiased estimator. ^{8}

\underline{Mixture of MANOVA and GMANOVA}. Let P : n × n = I − A_{1}(A′_{1}A_{1})^{−1}A′_{1} and let H be an (n − a_{1} − a_{2}) × n matrix such that H[A_{1}, A_{2}] = 0 and HH′ = I_{n−a_{1}−a_{2}}. Hence

H′H = I_{n} − [A_{1}, A_{2}] \{[ \begin{matrix}A′_{1} \\ A′_{2} \end{matrix} ] [A_{1}, A_{2}]\}^{−1} [ \begin{matrix}A′_{1} \\ A′_{2} \end{matrix} ] .

Let D_{2} : p × (p − c_{2}) be a matrix such that C_{2}D_{2} = 0 and D′_{2}D_{2} = I_{p−c_{2}}. The distribution of the matrix

[ \begin{matrix}X_{11} & X_{12} \\ X_{21} & X_{22} \\ . & . & . & . & . & . & . & . & . \\ Z \end{matrix} ] = [ \begin{matrix}(A′_{1}A_{1})^{−1/2}A′_{1} \\ (A′_{2}P_{2}A_{2})^{−1/2}A′_{2}P \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ H \end{matrix} ] Y [C′_{2}(C_{2}C′_{2})^{−1} ⋮D_{2}].

Use these equations to see that

nΣ̂ = C′(CS^{−1}C′)^{−1}C + SD(D′SD)^{−1}{D′SD + DY A(A′A)^{−1}A′Y D′}(D′SD)^{−1}D′S

and that

C′(CS^{−1}C′)^{−1}C = S − SD(D′SD)^{−1}D′S.

Combine these to obtain the desired result. Note that this estimate does not depend on choice of D. 12By (??) and S = Y′(I − A(A′A)^{−1}A′)Y,

S + SD(D′SD)^{−1}D′Y′A(A′A)^{−1}A′Y D(D′SD)^{−1}D′S = Y′(I − A(A′A)^{−1}A′)Y + [I − C′(CS^{−1}C′)^{−1}CS^{−1}]Y′A(A′A)^{−1}A′Y [I − C′(CS^{−1}C′)^{−1}CS^{−1}]′

= {Y − A(A′A)^{−1}A′Y S^{−1}C′(CS^{−1}C′)^{−1}C}′{Y − A(A′A)^{−1}A′Y S^{−1}C′(CS^{−1}C′)^{−1}C}.

8The MLE Σ̂ is biased and a trivial example of unbiased estimator is {1/(n−a)}S. But von Rosen (1991) gave another unbiased estimator based on the MLE, i.e.,

Σ̂_{U} = Σ̂ + uC′(CΣ̂^{−1}C′)^{−1}C,

where

u = \frac{a(n − a − 2(p − c) − 1}{(n − a − p + c − 1)(n − a − p + c)}.

is normal with mean

𝔼 [ \begin{matrix}X_{11} & X_{12} \\ X_{21} & X_{22} \\ . & . & . & . & . & . & . & . & . \\ Z \end{matrix} ] = [ \begin{matrix}(A′_{1}A_{1})^{1/2}β_{1}C′_{2}(C_{2}C′_{2})^{−1} + (A′_{1}A_{1})^{−1/2}A′_{1}A_{2}β_{2} & & & (A′_{1}A_{1})^{1/2}β_{1}D_{2} \\ (A′_{2}P A_{2})^{1/2}β_{2} & & & 0 \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 0 \end{matrix} ] .

Similar argument as that for GMANOVA gives that the maximum likelihood estimator of (β_{1}, β_{2}) is given by

β̂_{1} = (A′_{1}A_{1})^{−1}A′_{1}{Y − A_{2}β̂_{2}C_{2}}
β̂_{2} = (A′_{2}P A_{2})^{−1}A′_{2}P Y S^{−1}C′_{2}(C_{2}S^{−1}C′_{2})^{−1},

where S = Y′ [I_{n} − [A_{1}, A_{2}] \{[ \begin{matrix}A′_{1} \\ A′_{2} \end{matrix} ] [A_{1}, A_{2}]\}^{−1} [ \begin{matrix}A′_{1} \\ A′_{2} \end{matrix} ]] Y.

\underline{Extended GMANOVA due to von Rosen}. For i = 1, 2, …, m, let

Q_{i} : a_{i} × n = (A′_{i}A_{i})^{−1/2}A′_{i};
H_{i} : (n − a_{i}) × n  is a matrix such that H_{i}H′_{i} = I_{n−a_{i}} and H′_{i}H_{i} = I − Q′_{i}Q_{i};
H_{i+1}^{∗} : (a_{i} − a_{i+1}) × n  is a matrix such that H_{i+1}^{∗}(H_{i+1}^{∗})′ = I_{a_{i}−a_{i+1}} and
(H_{i+1}^{∗})′H_{i+1}^{∗} = Q_{i}Q_{i} − Q′_{i+1}Q_{i+1} with H_{m}^{∗} = 0;  (5)
C̃_{i} : c_{i} × (p − c_{i−1}^{+}) = C_{i}D_{1}D_{2} × ⋯ × D_{i−1};
D_{i} : (p − c_{i−1}^{+}) × (p − c_{i}^{+}) is a matrix such that C_{i}D_{i} = 0, D′_{i}D_{i} = I_{p−c_{i}^{+}} for i ≥ 1
with D_{0} = I_{p} and c_{i}^{+} = ∑_{j=1}^{i} c_{i} for i ≥ 1 and c_{0}^{+} = 0. (6)

Furthermore put, for i = 1, 2, …, m,

X_{[i]} : a_{i} × c_{i} = Q_{i}Y D_{0}D_{1} × ⋯ × D_{i−1}C̃′_{i}(C̃_{i}C̃′_{i})^{−1};  (7)
X_{<i>} : a_{i} × (p − c_{i}^{+}) = Q_{i}Y D_{0}D_{1} × ⋯ × D_{i};  (8)
Z_{[i]} : (n − a_{i}) × c_{i} = H_{i}Y D_{0}D_{1} × ⋯ × D_{i−1}C̃′_{i}(C̃_{i}C̃′_{i})^{−1};  (9)
Z_{<i>} : (n − a_{i}) × (p − c_{i}^{+}) = H_{i}Y D_{0}D_{1} × ⋯ × D_{i};  (10)
Z_{<i>}^{∗} : (a_{i} − a_{i}) × (p − c_{i}^{+}) = H_{i}^{∗}Y D_{0}D_{1} × ⋯ × D_{i};  (11)
μ_{[i]} : a_{i} × c_{i} = 𝔼[X_{[i]}]

= (A′_{i}A_{i})^{1/2}β_{i} + Q_{i}( ∑_{j=i+1}^{m} A_{j}β_{j}C_{j})D_{0}D_{1} × ⋯ × D_{i−1}C̃′_{j}(C̃_{j}C̃′_{j})^{−1};

μ_{<i>} : a_{i} × (p − c_{i}^{+}) = 𝔼[X_{<i>}]

= Q_{i} ∑_{j=i+1}^{m} A_{j}β_{j}C_{j}D_{0}D_{1} × ⋯ D_{i}

= Q_{i}A_{i+1}β_{i+1}C̃_{i+1} + Q_{i}( ∑_{j=i+2}^{m} A_{j}β_{j}C_{j})D_{0}D_{1} × ⋯ × D_{i};

Σ_{[i]} : c_{i} × c_{i} = cov[X_{[i]}, X_{[i]}];
Σ_{[i>} : c_{i} × (p − c_{i}^{+}) = Σ_{<i]} = cov[X_{[i]}, X_{<i>}];
Σ_{<i>} : (p − c_{i}^{+}) × (p − c_{i}^{+}) = cov[X_{<i>}, X_{<i>}].

If we put Σ_{<0>} = Σ and

Σ_{<i−1>} : (p − c_{i−1}^{+}) × (p − c_{i−1}^{+}) = [ \begin{matrix}Σ_{[i]} & Σ_{<i]} \\ Σ_{<i]} & Σ_{<i>} \end{matrix} ] for i = 1, 2, …, m,

we can get the recursive formula for the covariance matrix

Σ_{<i−1>} = ( \begin{matrix}(C̃_{i}C̃′_{i})^{−1}C̃_{i}Σ_{<i−1>}C̃′_{i}(C̃_{i}C̃′_{i})^{−1} & (C̃_{i}C̃′_{i})^{−1}C̃_{i}Σ_{<i−1>}D_{i} \\ D′_{i}Σ_{<i−1>}C̃′_{i}(C̃_{i}C̃′_{i})^{−1} & D′_{i}Σ_{<i−1>}D_{i} \end{matrix} ) = ( \begin{matrix}Σ_{[i]} & Σ_{[i>} \\ Σ_{<i]} & Σ_{<i>} \end{matrix} ) .

Furthermore, we put

μ′_{[i]·} = μ′_{[i]} − Σ_{[i>}Σ_{<i>}^{−1}μ′_{<i>},
P_{i} : (p − c_{i−1}^{+}) × (p − c_{i−1}^{+}) = Σ_{<i−1>}^{−1}C̃′_{i}(C̃_{i}Σ_{<i−1>}^{−1}C̃′_{i})^{−1}C̃_{i}.

Then we have, for i = 1, 2, …, m,

X_{i} : a_{i} × (p − c_{i−1}) = (X_{[i]}, X_{<i>}) ∼ N([μ_{[i]}, μ_{<i>}].I_{a_{i}} ⊗ Σ_{<i−1>}),
Z_{i} : (n − a_{i}) × (p − c_{i−1}) = (Z_{[i]}, Z_{<i>}) ∼ N([0, I_{n−a_{i}} ⊗ Σ_{<i−1>}),  (12)

and

(X′_{[i]}, Z′_{[i]})′ : n × c_{i} ∼ N(μ′_{[i]}, 0′)′, I_{n} ⊗ Σ_{[i]}),
(X′_{<m>}, Z′_{<m>})′ : n × (p − c_{m}^{+})_{i} ∼ N(0, I_{n} ⊗ Σ_{<m>}).

For i = 1, 2, …, m + 1, put

S_{<i−1>} + U_{<i−1>} : (p − c_{i−1}^{+}) × (p − c_{i−1}^{+}) = [ \begin{matrix}S_{[i]} & S_{[i>} \\ S_{<i]} & S_{<i>} \end{matrix} ] = [ \begin{matrix}Z′_{[i]}Z_{[i]} & Z′_{[i]}Z_{[i>} \\ Z_{<i]}Z_{[i]} & Z′_{<i>}Z_{<i>} \end{matrix} ] (13)

and

U_{<i>} : (p − c_{i}^{+}) × (p − c_{i}^{+}) = (Z_{<i>}^{∗})′Z_{<i>}^{∗},

where S_{[i]} are c_{i} × c_{i} matrices and U_{<0>} = 0. Then we have, for i = 1, 2, …, m + 1,

S_{<i−1>} + U_{<i−1>} ∼ W_{p−c_{i−1}^{+}}(n − a_{i}, Σ_{<i−1>}) and S_{<i−1>} ∼ W_{p−c_{i−1}^{+}}(n − a_{i−1}, Σ_{<i−1>})

where a_{m+1} = 0.

Now we decompose the density function of Y into canonical form whose parameter space is expressed as a product of simpler subspaces. This parameterization is regarded as Mparameter space in [AMP]. Theorem. The density of Y can be factored as

∏_{i=1}^{m+1} \{ | Σ_{[i]·} |^{−a_{i}/2} exp (−\frac{1}{2}tr [\{X_{[i]} + X_{<i>}Σ_{<i>}^{−1}Σ_{<i]} − μ_{[i]·}\}Σ_{[i]·}^{−1}\{⋯\}′]) dX_{[i]}

× | Σ_{[i]·} |^{−(p−c_{i}^{+})/2}| S_{<i>} |^{a_{i}/2} exp (−\frac{1}{2}tr [S_{<i>}\{S_{<i>}^{−1}S_{<i]} − Σ_{<i>}^{−1}Σ_{<i]}\}Σ_{[i]·}^{−1}\{⋯\}′]) d(S_{<i>}^{−1}S_{<i]})

× | Σ_{[i]·} |^{−\{n−a_{i}−(p−c_{i}^{+})\}/2}| S_{[i]·} |^{\{n−a_{i}−(p−c_{i}i^{+})−c_{i}−1\}/2} exp (−\frac{1}{2}tr [S_{[i]·}Σ_{[i]·}^{−1}]) d(S_{[i]·})\},  (14)

where S_{[m+1]·} = S_{<m>}, X_{[m+1]} = 0, S_{<m+1]} = 0 and a_{m+1} = 0. Proof. For i = 1, 2, …, m, let

Ỹ_{i} : a_{i−1} × (p − c_{i−1}^{+}) = Q_{i−1}(Y − ∑_{j=i}^{m} A_{j}β_{j}C_{j})D_{0} × ⋯ × D_{i−1};

Ỹ_{[i]} : a_{i} × c_{i} = Q_{i}(Y − ∑_{j=i}^{m} A_{j}β_{j}C_{j})D_{0} × ⋯ × D_{i−1}C̃′_{i}(C̃_{i}C̃′_{i})^{−1};

Ỹ_{<i>} : a_{i−1} × (p − c_{i}^{+}) = Q_{i}(Y − ∑_{j=i}^{m} A_{j}β_{j}C_{j})D_{0} × ⋯ × D_{i}.

From the orthogonality of H_{i} to A_{i}, …, A_{m} and of D_{i} to C_{i} and the definition of (??)-(??),

Z_{[i]} = H_{i}Ỹ_{1}D_{0}D_{1} × ⋯ × D_{i−1}C̃′_{i}(C̃_{i}C̃′_{i})^{−1};
Z_{<i>} = H_{i}Ỹ_{1}D_{0}D_{1} × ⋯ × D_{i};
Z_{<i>}^{∗} = H_{i}^{∗}Ỹ_{1}D_{0}D_{1} × ⋯ × D_{i}.

From

Q′_{1}Q_{1} + H′_{1}H_{1} = I_{n} and C′_{1}(C_{1}C′_{1})^{−1}C_{1} + D_{1}D′_{1} = I_{p},

we have

tr {Ỹ_{1}Σ_{<0>}^{−1}Ỹ′_{1}} = tr {Ỹ Σ^{−1}Ỹ′}
= tr \{ [ Ỹ_{[1]} Ỹ_{<1>} ] [ \begin{matrix}C_{1}Σ^{−1}C′_{1} & C_{1}Σ^{−1}D_{1} \\ D′_{1}Σ^{−1}C′_{1} & D′_{1}Σ^{−1}D_{1} \end{matrix} ] [ Ỹ′_{[1]} ; Ỹ′_{<1>} ] \}
+tr \{ [ Z_{[1]} Z_{<1>} ] [ \begin{matrix}C_{1}Σ^{−1}C′_{1} & C_{1}Σ^{−1}D_{1} \\ D′_{1}Σ^{−1}C′_{1} & D′_{1}Σ^{−1}D_{1} \end{matrix} ] [ \begin{matrix}Z′_{[1]} \\ Z′_{<1>} \end{matrix} ] \}. (15)

Noting that

[ \begin{matrix}C_{1}Σ^{−1}C′_{1} & C_{1}Σ^{−1}D_{1} \\ D′_{1}Σ^{−1}C′_{1} & D′_{1}Σ^{−1}D_{1} \end{matrix} ] = [ \begin{matrix}I & 0 \\ D′_{1}Σ^{−1}C′_{1}(C_{1}Σ^{−1}C′_{1})^{−1} & I \end{matrix} ] [ \begin{matrix}C_{1}Σ^{−1}C′_{1} & 0 \\ 0 & (D′_{1}ΣD_{1})^{−1} \end{matrix} ]
× [ \begin{matrix}I & 0 \\ D′_{1}Σ^{−1}C′_{1}(C_{1}Σ^{−1}C′_{1})^{−1} & I \end{matrix} ]′ ,  (16)
[ \begin{matrix}Σ_{[1]} & Σ_{[1>} \\ Σ_{<1]} & Σ_{<1>} \end{matrix} ] = [ \begin{matrix}C_{1}Σ^{−1}C′_{1} & C_{1}Σ^{−1}D_{1} \\ D′_{1}Σ^{−1}C′_{1} & D′_{1}Σ^{−1}D_{1} \end{matrix} ]^{−1} ,  (17)
Σ_{[1]·}^{−1} = C_{1}Σ^{−1}C′_{1}.

we can see that the first term of the right hand side in (??) equals

tr [{X̃_{[1]} + X̃_{<1>}D′_{1}Σ^{−1}C′_{1}(C_{1}Σ^{−1}C_{1})^{−1}}C_{1}Σ^{−1}C′_{1}{⋯}′] + tr [Ỹ_{<1>}(D′_{1}ΣD_{1})^{−1}Ỹ′_{<1>}]

= tr [{X_{[1]} + X_{<1>}Σ_{<1>}^{−1}Σ_{<1]} − μ_{[1]·}}Σ_{[1]·}^{−1}{⋯}′] + tr [Ỹ_{<1>}Σ_{<1>}^{−1}(Ỹ_{<1>})′]
= tr [{X_{[1]} + X_{<1>}Σ_{<1>}^{−1}Σ_{<1]} − μ_{[1]·}}Σ_{[1]·}^{−1}{⋯}′] + tr [Ỹ_{2}Σ_{<1>}^{−1}Ỹ′_{2}]
+tr [Z_{<1>}^{∗}Σ_{<1>}^{−1}(Z_{<1>}^{∗})′]  (18)

The last equality in (??) follows from Q′_{1}Q_{1} = (H_{1}^{∗})′H_{1}^{∗}+Q′_{2}Q_{2}. From (??)-(??), we conclude that

tr [Ỹ_{1}Σ_{<0>}^{−1}Ỹ′_{1}] = tr [{X_{[1]} + X_{<1>}Σ_{<1>}^{−1}Σ_{<1]} − μ_{[1]·}}Σ_{[1]·}^{−1}{⋯}′]
+tr [S_{[1]·}Σ_{[i]·}^{−1}] + tr [S_{<1>}{S_{<1>}^{−1}S_{<1]} − Σ_{<1>}^{−1}Σ_{<1]}}Σ_{[1]·}^{−1}{⋯}′]
+tr [(S_{<1>} + U_{<1>})Σ_{<1>}^{−1}] + tr [Ỹ_{2}Σ_{<1>}^{−1}Ỹ′_{2}].  (19)

Finally using the standard decomposition

| Σ |= ∏_{1}^{m+1} | Σ_{[i]·} |

and applying similar argument through (??) to (??), the formula (??) can be factored into (??).

From Theorem, we can see that a decomposition of the model (??) can be made as

X_{[i]} : a_{i} × c_{i} given X_{<i>} ∼ N(μ_{[i]·} + X_{<i>}γ_{<i]}, I_{a_{i}} ⊗ Σ_{[i]·})
X_{<i>} : a_{i} × (p − c_{i}^{+}) ∼ N(μ_{<i>}, I_{a_{i}} ⊗ Σ_{<i>})
F_{<i]}; (p − c_{i}^{+}) × c_{i} = S_{<i]}S_{[i]} given S_{<i>} ∼ N(γ_{<i]}, S_{<i>}^{−1} ⊗ Σ_{[i]·})
S_{[i]·} = S_{[i]} − S_{[i>}S_{<i>}^{−1}S_{<i]} ∼ W_{c_{i}}(n − a_{i} − (p − c_{i}^{+}), Σ_{[i]·}),  (20)

where (X_{[i]}, X_{<i>}), (F_{<i]}, S_{<i>}), S_{[i]·} are mutually independent.

In [AMP], the family of parameters (μ_{[i]·}, Σ_{<i>}Σ_{<i]}, Σ_{[i]·}) is termed as M-parameters. As noted in [AMP], there is a one-to-one correspondence between the original parameters and M-parameters and reconstruction algorithm from M-parameter to the original parameter is given in [AMP].

3. Further discussion on GMANOVA

\underline{Canonical reduction of testing problem}^{9} Consider the problem of testing the null hypothesis

H : EβF = 0, v.s. K : EβF ≠ 0,

where E : e × a and F : c × f. Let M : (a − e) × a and N : c × (c − f) such that

E^{∗}M′ = 0, MM′ = I_{a−e}, F′N = 0, N′N = I_{c−f}

where E^{∗} = E(A′A)^{−1/2}. Put

W = [ \begin{matrix}W_{11} & W_{12} & W_{13} \\ W_{21} & W_{22} & W_{23} \\ W_{31} & W_{32} & W_{33} \end{matrix} ] = [ \begin{matrix}(E^{∗}E^{∗}′)^{−1/2}E^{∗} & 0 \\ M & 0 \\ 0 & 0 \end{matrix} ] [ \begin{matrix}X_{1} ⋮ & X_{2} \\ . & . & . & . & . & . & . & . \\ Z \end{matrix} ] [ \begin{matrix}N & F (F′F)^{−1/2} & 0 \\ 0 & 0 & 0 \end{matrix} ]

and

[ \begin{matrix}β_{11}^{∗} & β_{12}^{∗} & 0 \\ β_{21}^{∗} & β_{22}^{∗} & 0 \\ 0 & 0 & 0 \end{matrix} ] = [ \begin{matrix}(E^{∗}E^{∗}′)^{−1/2}E^{∗} & 0 \\ M & 0 \\ 0 & 0 \end{matrix} ] [ \begin{matrix}(A′A)^{1/2}β ⋮ & & & 0 \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 0 \end{matrix} ] [ \begin{matrix}N & F (F′F)^{−1/2} & 0 \\ 0 & 0 & 0 \end{matrix} ] .

Then null hypothesis H : EβF = 0 is equivalent to H : β_{12}^{∗} = 0.

Now consider the the groups which leave invariant the testing problem, i.e.,

U = \{U : p × p = [ \begin{matrix}U_{11} & 0 & 0 \\ U_{21} & U_{22} & 0 \\ U_{31} & U_{32} & U_{33} \end{matrix} ] |U_{11} ∈ GL(c − f), U_{22} ∈ GL(f), U_{33} ∈ GL(p − c)\}

B = \{B : n × p = [ \begin{matrix}B_{11} & 0 & 0 \\ B_{21} & B_{22} & 0 \\ 0 & 0 & 0 \end{matrix} ] |B_{11} ∈ M(e × (c − f)), B_{22} ∈ M((e − a) × f)\}

O = \{O : n × n = [ \begin{matrix}O_{11} & 0 & 0 \\ 0 & O_{22} & 0 \\ 0 & 0 & O_{33} \end{matrix} ] |O_{11} ∈ O(e), O_{22} ∈ O(e − a), O_{33} ∈ O(n − a)\} .

9See Gleser and Olkin (1970), Muirhead (1982), p. 520-525 and Kariya (1985), Chapter 3.

Consider the action

Z ⟶ OW U + B.

To obtain analytically tractable maximal invariant statistics, consider a subgroup, O×U ×B with O_{11} = I, O_{22} = I. Then a maximal invariant statistics under the subgroup is

(W_{12}, W_{13}) [ \begin{matrix}V_{22} & V_{23} \\ V_{32} & V_{33} \end{matrix} ]^{−1} (W_{12}, W_{13})′, [ \begin{matrix}W_{13} \\ W_{23} \end{matrix} ] V_{33}^{−1} [ W′_{13} W′_{23} ] ,

where V_{ij} = W′_{3i}W_{3j}. We may show that likelihood ratio statistics for testing H : β_{12}^{∗} = 0 against K : H : β_{12}^{∗} ≠ 0 is

λ^{2/n} = \frac{det(I + W_{13}V_{33}^{−1}W′_{13})}{det (I + (W_{12}, W_{13}) [ \begin{matrix}V_{22} & V_{23} \\ V_{32} & V_{33} \end{matrix} ]^{−1} (W_{12}, W_{13})′) .}

The LRS λ is a function of (W, V) only through the maximal invariant statistics, but it does not involve W_{23} whereas the maximal invariant does.

\underline{Inadmissibility of the MLE’s} . First let us consider the estimation of β. Then

X_{1} : a × c|X_{2} ∼ N(β + X_{2}R, I_{a} ⊗ Σ_{11.2}),  X_{2} : a × (p − c) ∼ N(0, I_{a} ⊗ Σ_{22}),
S_{11.2} ∼ W_{c}(n − a − p + c, Σ_{11.2}),  S_{22}^{−1}S_{21}|S_{22} ∼ N(R, S_{22}^{−1} ⊗ Σ_{11.2}).

Consider a loss function L(β̂, β) = tr{(β̂ − β)Σ_{11.2}^{−1}(β̂ − β)′}. Estimators we look at are of the form

β̂_{G} = X_{1} − X_{2}{S_{22}^{−1}S_{21} + G(S_{22}^{−1}S_{21}, S_{22}, S_{11.2}, X_{2})}.

Then the risk function is written as

𝔼[tr{(β̂_{G} − β)Σ_{11.2}^{−1}(β̂_{G} − β)′}]
= 𝔼^{X_{2}}𝔼^{X_{1}|X_{2}}[{X_{1} − (β + X_{2}R) − X_{2}(S_{22}^{−1}S_{21} + G − R)}Σ_{11.2}^{−1}{⋯}′]
= 𝔼[{X_{1} − β − X_{2}R}Σ_{11.2}^{−1}{⋯}′]
+𝔼[X_{2}X′_{2}{S_{22}^{−1}S_{21} + G − R}Σ_{11.2}^{−1}{S_{22}^{−1}S_{21} + G − R}′].  (21)

The second term of RHS of (??) is the expectation of weighted sum of squared error of estimator S_{22}^{−1}S_{21} + G based on S_{22}^{−1}S_{21}|S_{22} ∼ N_{p−c}(R, S_{22}^{−1} ⊗ Σ) and S_{11.2} ∼ W_{c}(n − a − p + c, Σ_{11.2}), given X_{2}, S_{22}. Hence we may show that β̂ is inadmissible, using the argument to prove the inadmissibility of mean vector of normal distribution. Careful application of argument by Brown (1990) ( general discussion on inadmissibility of estimators which depend on the distribution of the ancillary statistics ) may leads to the inadmissibility of β̂ when β is a scalar and p − c ≥ 3. Also see Kariya et al (1996) for extensive discussion on shrinkage estimators of regression coefficient in GMANOVA model.

Estimation of covariance matrix in GMANOVA can be reduced to estimation of covariance matrix under monotone missing data (simplest one!).

References

Andersson, S.T., Marden, J.I. and Perlman, M.D. (1993). Totally ordered multivariate linear models. Sankhyā 55 Series A 370-394.

Brown, L.D. (1990). An ancillarity paradox which appears in multiple linear regression Ann. Statist. 18, 471-538.

Chinchilli, V.M. and Elswick, R.K. (1985). A mixture of the manova and gmanova models. Comm. Statist. -Theor. Meth. 14, 3075-3089.

Gleser, L.J. and Olkin, I. (1970). Linear models in multivariate analysis. Essays in Probability and Statistics. 267-292, Wiley, New York.

Kariya, T. (1985). Testing in the multivariate general linear model. Kinokuniya company LTD. Tokyo, Japan.

Kariya, T., Konno, Y., Strawderman, W.E. (1996). Double shrinkage estimators in the GMANOVA model. J. multivariate analysis 56, 245-258.

Muirhead, R.J. (1982). Aspects of multivariate statistical theory. John Wiley & Sons, INC.

Patel, H.I. (1986). Analysis of the repeated measures designs with changing covariates in clinical trials. Biometrika 73, 707-715.

Potthoff, R.F. and Roy, S.N. (1964). A generalized multivariate analysis of variance model useful especially for growth curve problems Biometrika 51, 313-326.

von Rosen, D. (1989). Maximum likelihood estimators in multivariate linear normal models. J. multivariate analysis 31, 187-200.

von Rosen, D. (1991). Moments of maximum likelihood estimators in the growth curve model. statistics 22, 111-131.

Verblya, A.P. (1988). Analysis of repeated designs with changing covariates. Biometrika 75, 172-174.

Verblya, A.P. and Venables, W.N. (1988). An extension of the growth curve model. Biometrika 75, 129-138.