Geometry, Statistics and Decision Making in Gene Therapy

Victor Pătrângenaru

Abstract

One emphasizes the role of geometry in normal parametric statistics and introduces the pi-beta distributions. Using a standard geometric interpretation of a random effect model, one explains the origin of tests for means useful in a statistical analysis of the effect of adenovirus treatments of tumors in vitro and in vivo.

Introduction

This paper is addressed to a large audience, some of which at times has to teach an introductory course in analytic geometry. It the past years some departments of mathematics, have eliminated such courses from their undergraduate curriculum.

The present paper comes with a very concrete example as of why a multidimensional geometry course is necessary both in a mathematics or statistics department.

One of the most important results in probability is the central limit theorem (CLT). The CLT states that variates which are sums of many independent and identically distributed effects tend to be normally distributed as the number of effects becomes large [10].

A normally distributed random variable has the probability density function

(0.1)  dp = \frac{1}{\sqrt{2πσ}} exp (−\frac{1}{2σ^{2}} (x − μ^{2})) dx.

When possibly related measurements such as heights and weights of individuals of same group are made, the quadratic term under the exponential in (0.1) is substituted with a quadratic form in more variables. It does not take an eye of an expert to see that the hidden correlation parameters (coefficients of the quadratic forms) of these measurements can be expressed in terms of some Euclidean objects such as various projections, volumes, etc. It is less obvious , still straightforward [2], [10] to find the distribution function of these measurements, under certain linear hypothesis. This will be taken care of in the first two sections. The Wishart distributions are derived in section 1. In section 2, we will make standard geometric assumptions on the covariance matrices, to write computable tests in terms of ANOVA (analysis of variance) tables. These tests are used to decide to what extent the mean value of a multivariate normal distribution of an orthogonal model, lies in a vector subspace L versus being in a larger subspace M. In very few cases such a test may be based on the popular F- distribution; in general they are unknown, and under additional splitting constraints on L and M, one is lead to pi-beta tests, that is products tests based on powers of independently beta distributed r.v.’s, which are rarely tabulated, leading to interesting numerical problems.

Balkan Journal of Geometry and Its Applications, Vol.2, No.2, 1997, pp. 83-100 \begin{matrix}c \\ ◯ \end{matrix}Balkan Society of Geometers, Geometry Balkan Press

Which brings us to the second goal of this paper, a question regarding decision making in recent gene therapy models experimented at Dr. Milton W. Taylor’s biology lab [11]. We will put it in two ways:

One is do decide if there is a significant difference at the end of the application of various treatments in vitro of randomly the MDA-MB-435 breast cancer cells. This leads to an F test.

The other concerns the difference in time between treatments with recombinant adenovirus on K562 melanoma cells in vivo, in animal models. We will explain in details the design we use, which is called by some a 3 way ANOVA. This leads to a pi-beta test.

It is interesting to note that other modern approved immunotherapeutical protocols [9], using T IL ( tumor infiltrating leukocytes ) who have an established 40% response in patients with melanoma or kidney cancer, were not successful in breast carcinoma [6], which makes the successful preclinical trials in [11] even more important.

The relationship between geometry and statistics, goes beyond this introductory paper [5]. Parametric statistics is viewed by some researchers as the study of statistical manifolds, which are Riemannian manifolds of distributions, where the metric is the Fisher information on the space of parameters [1] .

1 The geometry of multivariate normal distribution

Multivariate analysis deals with ordered lists of data, representing a number of aspects of the same phenomenon. These numbers or measurements are indexed by the ordered set I, of size p, so that one such list can be regarded as a matrixvalued random variable (r.v.) X = (X_{n}), where X_{n} are independent multivariate normal distributed r.v.’s, with same covariance matrix Σ = (σ^{ij})_{(i,j)∈I×I} , that is has the probability density function (p.d.f.) :

(1.1)  (2π)^{− \frac{p}{2}} det(Σ)^{− \frac{1}{2}} exp (−\frac{1}{2} < x_{n} − μ_{n}, Σ^{−1}(x_{n} − μ_{n}) >) ,

where <, > is the usual scalar product on R^{p}.

Assume X_{n} are identically distributed. We would like to make inferences about the ”hidden parameters” (μ, Σ) in R^{p} × P (I) from the given numerical data (x_{n})_{n=\overline{1,N}}.

In this respect we will follow the maximum likelihood estimator (MLE) [2] approach .

The joint distribution function of X = (X_{n})_{n=\overline{1,N}}, is a multiple of the so called the likelihood function

(1.2) L(μ, Σ; x) = det(Σ)^{− \frac{N}{2}} exp (−\frac{1}{2} ∑_{n=\overline{1,N}} < x_{n} − μ, Σ^{−1}(x_{n} − μ) >) .

The MLE of (μ̂, Σ̂) = (μ̂, Σ̂)(x) is the point of maximum of L(·, ·; x) and may be easily determined with a geometric argument. Indeed, for Σ fixed, define on R^{p} the scalar product < ·, · >_{Σ}:=< ·, Σ^{−1}(·) >, which obviously extends to a scalar product on (R^{p})^{N}, labeled << ·, · >>_{Σ} .

Let Δ = Δ((R^{p})^{N}) be the diagonal of elements of the form (x, x, ..., x), x ∈ R^{p}, and let P be the orthoprojection of ((R^{p})^{N}, << ·, · >>_{Σ}) on Δ, which is actually independent of Σ.

Let \overline{x} be the baricenter with equal weights of the R^{p} - components of x = (x_{n})_{n=\overline{1,N}}. Then

(1.3)  P ((x_{n})_{n=\overline{1,N}}) = (\overline{x})_{n=\overline{1,N}}.

Q = Id − P is the orthoprojection on the orthocomplement Δ^{⊥} (again w.r.t. any of the scalar products << ·, · >>_{Σ}). It is obvious that

(1.4)  ∥x − (μ, μ, ..., μ)∥_{Σ}^{2} ≥ ∥Q(x)∥_{Σ}^{2}

and the equality holds iff μ = \overline{x}. The value of the likelihood function at such points is

(1.5)  L(\overline{x}, Σ; x) = det(Σ)^{− \frac{N}{2}} exp(−\frac{1}{2} ∥Q(x)∥_{Σ}^{2}).

Assume \{x_{n} − \overline{x}\}_{n=\overline{1,N}} spans R^{p}, which happens iff the rows of Q(x) are l.i., or if the Gramm determinant det(Q(x) Q(x)^{t}) is positive. Then if Ω is a positive definite square root of Σ^{−1}, then

(1.6)  L(\overline{x}, Σ; x) =

= det(Q(x) Q(x)^{t})^{− \frac{N}{2}} det(ΩQ(x) Q(x)^{t}Ω)^{\frac{N}{2}} exp(−\frac{1}{2}T r(ΩQ(x) Q(x)^{t}Ω) ).

If we set Λ = ΩQ(x) Q(x)^{t}Ω, then L(\overline{x}, Σ; x) is proportional with

G(Λ) = det(Λ)^{\frac{N}{2}} exp(−\frac{1}{2}T rΛ)

and since Λ is symmetric, there is an orthogonal matrix T, such that T ΛT^{t} = Δ is diagonal. Then G(Λ) = G(Δ) and w.l.o.g., one may suppose that Λ = diag Λ_{i}. Therefore

g(Λ) = ln G(Λ) = ∑_{i∈I} (\frac{N}{2} ln Λ_{i} − \frac{1}{2}Λ_{i}) .

The max of u(t) = \frac{1}{2}(N ln t − t) occurs for t = N, and therefore the max of

g(Λ) occurs when Λ_{i} = N, that is Λ = N Id, or Σ^{−1}Q(x) Q(x)^{t} = N Id.

Let Φ : R^{p} ⊗ R^{N} → (R^{p})^{N} be the isomorphism Φ(x ⊗ y) = (y^{1}x, ..., y^{N}x), and for a vector subspace V of R^{p} consider S = Φ(R^{p} ⊗ V). Since ΣS = S, for any Σ in P (I), S^{⊥} is independent of Σ.

Assume the mean \underline{μ} = (μ_{n}) of the distribution in (1.1) is in S. Along the same lines, if P_{S} is the orthoprojector on S, and Q_{S} = 1 − P_{S}, we get at no additional effort the following classical result: Theorem 1.1 If N ≥ dim M + p. Then the MLE of L(\underline{μ}, Σ; x) exists in general and \underline{μ̂} = P_{S}, Σ̂ = \frac{1}{N} Q_{S}(x) Q_{S}(x)^{t}.

Let D_{S} = {x ∈ R^{p} | Q_{S}(x) Q_{S}(x)^{t}is positive definite} and assume F_{S} associates the MLE with data, that is

F_{S} : D_{S} → S × P (I)

F_{S}(x) = (P_{S}(x), \frac{1}{N} Q_{S}(x) Q_{S}(x)^{t}).

In order to answer specific inference questions, one has to find the distribution function of F_{S}(x). Behind the answer to this standard question, there is again a nice geometric idea. If Ω is positive definite on W ⊕U and W, U are Ω -invariant then for any W ⊕U-valued multivariate normally distributed r.v. Y of covariance matrix Ω, the U and W parts are normally distributed and independent as r.v.’s. It turns out that P_{S}(x) is normally distributed and furthermore \frac{1}{N} Q_{S}(x) Q_{S}(x)^{t} and P_{S}(x) are independent. Moreover if we select a basis in R^{N} , the first vectors of which are in M, then S = (R^{p})^{d} × {0} and S^{⊥} = {0} × (R^{p})^{N−d}. As such, w.l.o.g. one may assume that z = Q_{S}(x) = (x_{1}, x_{2}, ..., x_{N−d}) is a sample of N − d independent N(0, Σ)-distributed r.v’s. As such we are interested only in a distribution of θ = t(z) = z z^{t}, given that the joint distribution of z is

d N(z) = det(Σ)^{− \frac{m}{2}} exp (−\frac{1}{2}Tr(Σ^{−1}z z^{t})) dλ(z), m = N − d.

By the change of variable formula, θ has the density

(1.7)  det(Σ)^{− \frac{m}{2}} exp (−\frac{1}{2}Tr(Σ^{−1}Θ)) dt(λ)(Θ).

In order to determine dt(λ), we make a digression on invariant measures we learned from the senior author of [4].We are looking for GL(p) -invariant measure on P (I), given that GL(p) acts transitively on P (I) by

(1.8)  (A, Σ) → AΣA^{t}.

Let G be a group, acting on Y, on the left. If μ is a measure on Y, and g is a transformation of Y, g^{−1}μ is the measure defined by g^{−1}μ (M) =: μ(gM).

Definition a. A measure μ on Y is relatively invariant, with multiplicator χ, if for any g in G, g^{−1}μ = χ(g)μ. Since χ has to be positive, and it follows that χ(1_{G}) = 1. b. A measure μ on Y is invariant, if μ is relatively invariant with multiplicator 1. Example 1.1 The Lebesgue measure on R^{p} is relatively invariant w.r.t. the natural action of the affine group GL(p) × R^{p} on R^{p} with multiplicator

χ((g, u)) = det(g).

Remark 1.1. Assume μ is a relatively invariant measure on Y with multiplicator χ, and n : Y → R is a measurable function, such that n(gx) = χ(g)n(x). Then n^{−1}μ is an invariant measure. Lemma 1.1. The Lebesgue measure λ on P (I), is relatively invariant w.r.t. the natural action (1.6) of GL(p) on P (p) with multiplicator χ(g) = (det g)^{p+1}. Proof (sketch). Assume A = cId, then A · Θ = AΘA^{t} = c^{2}Θ . In this case

A^{−1}λ(S) = ∫_{c^{2}S} dτ_{1}^{1}...dτ_{p}^{p} = c^{2(p+1)p/2} ∫_{S} dτ_{1}^{1}...dτ_{p}^{p} = (det A)^{p+1}λ(S).

It is obvious then that n(Θ) = (det Θ)^{\frac{p+1}{2}} satisfies the conditions in Remark 1.1. and therefore Corollary 1.1. Let (A, Θ) → AΘA^{t} be the transitive action of GL(p) on P (I) and λ be the Lebesgue measure on P (I) = GL(p)/O(p). Then density function relatively to the Lebesgue measure given by (det Θ)^{− \frac{(p+1)}{2}} defines GL(p) invariant measure on P (I).

We recall that the Lebesgue measure λ is GL(p) invariant, that is if A is in GL(p) transformation will justify the formula:

A^{−1}λ = | det A|^{p}λ.

Notice that GL(p) acts on the left on (R^{p})^{m}. Since t is equivariant, t(λ) is relatively invariant with multiplicator

χ(A) = | det A|^{m}

since

A^{−1}t(λ) (G) = tλ(A G) = ∫_{{z,zz^{t}∈AGA^{t}}} dλ_{z} = ∫_{{z,A^{−1}z(A^{−1}z)^{t}∈G}} dλ_{z} =

= | det A|^{m} ∫_{{u,uu^{t}∈G}} dλ_{u} = | det A|^{m}t(λ) (G).

Let then n : P (I) → R, n(Θ) := (det Θ)^{\frac{m}{2}} , satisfies to the condition in Remark 1.1, and thus (det Θ)^{− \frac{m}{2}} t(λ) is an invariant measure on P (I). Two invariant measures are equal up to a constant on each orbit of the group action. In our case since P (I) is homogeneous for the action (1.6) is transitive, it follows that

dt(λ)_{Σ} = (det Θ)^{\frac{m}{2}} (det Θ)^{− \frac{p+1}{2}} dλ_{P (I)}.

Then, due to (1.8), we have shown that Proposition 1.2. If z is a sample of N(0, Σ), then the p.d.f. of z z^{t} := Θ is proportional to

(1.9)  det(Σ)^{− \frac{m}{2}} det(Θ)^{\frac{m−p−1}{2}} exp(−\frac{1}{2}T r(Σ^{−1}Θ)).

The probability distribution proportional to (1.9), is named Wishart distribution [2] with n degrees of freedom, and form parameter Σ, and is labeled as W (Σ, m). Corollary 1.2. Assume X = (X_{n})_{n=\overline{1,N}} is a sample of N(\underline{μ}, Σ), μ ∈ S. Then the sample variance \frac{N}{N−d} Σ̂ is distributed as W (\frac{1}{N−d}Σ, N − d).

2 Orthogonal designs

Many problems in multivariate statistics are concerned with test for means. One would like to estimate from experimental data the subspace S, where the parameter \underline{μ} of the N(\underline{μ}, Σ) distributed r.v. X = (X_{i})_{i=\overline{1,N}}. Such tests can be expressed in terms of some Wishart distributions. Still given the peculiar form of the subspace S, we would like to relax the conditions on Σ. We will assume that X = (X_{i})_{i=\overline{1,N}} is a sample from R^{p}-valued r.v. that is N(μ, Σ) - distributed.

a. To start with let us first assume that Σ has the simplest possible form, that is Σ = σ^{2}Id. In this case the likelihood function is

(2.1)  L(μ, σ; x) = σ^{−Np} exp (−\frac{σ^{−2}}{2} ∑ ∥x_{n} − μ∥^{2}) .

Let L ⊆ M be two vector subspace of R^{p}. We would like to test to what extent the hypothesis H : μ is in L, versus μ is in M, holds true ([10], p.33 ). Let pr_{L} denote the orthoprojection of R^{p} on L and \overline{x} as before denote the baricenter with equal weights of (x_{n})_{n=\overline{1,N}}. Remark 2.1. Since there is only one variance component, one may substitute list the data in a single observation of size N p, and substitute L by the diagonal of L^{N}, as subspace of (R^{p})^{N}. As such w.l.o.g., one may assume that N = 1.

A straightforward computation shows that under the hypothesis that μ is in L, the MLE of (μ , σ^{2}) is

(2.2)  (μ̂_{L}, σ̂_{L}^{2}) = (pr_{L}\overline{x}, \frac{1}{p}∥Q_{L}(x)∥^{2})

and the corresponding maximum of the likelihood function is

(2.3)  L̂_{L} = (σ̂_{L}^{2})^{− \frac{p}{2}} exp(−\frac{p}{2}).

The components of x are uncorrelated, and being normally distributed, they are independent. If d = dim L, since σ^{−1}Q_{L}(x) is a vector valued r.v. all of that splits into a sum of independent unit normally distributed r.v.’s \frac{1}{p−d}σ^{−2}∥Q_{L}(x)∥^{2} follows a chi-square distribution with p − d degrees of freedom, χ_{p−d}^{2}. The likelihood ratio associated whit the hypothesis H, is

L̂_{L}/L̂_{M} = (σ̂_{M}^{2}/σ̂_{L}^{2})^{p/2}.

Before we proceed further we recall some facts on beta distributions. The beta integrals with form parameters (a, b) are B(a, b) = ∫_{0}^{1} r^{a−1}(1 − r)^{b−1}dr. The beta distribution with parameters (a, b) has the density function on (0, 1) :

B(a, b)^{−1}r^{a−1}(1 − r)^{b−1}.

Recall also that Γ(a) = ∫_{0}^{∞} x^{a−1} exp(−x)dx and the gamma distribution with form parameter a has the density function on (0, ∞) :

Γ(a)^{−1}x^{a−1} exp(−x).

A straightforward calculation shows that if x, y are independent r.v’s with gamma distributions of form parameters a, b, then x + y, r = x/(x + y) are also independent and r has a beta distribution with form parameters (a, b). In particular, since a χ_{m}^{2} distributed r.v, follows a gamma distribution with form parameter \frac{m}{2}, from the above argument, it turns out that

Proposition 2.1. σ̂_{M}^{2}/σ̂_{L}^{2} has a distribution function g on the interval (0, codimM/codimL)given by

(2.4)  g(u) = codimL/codimMf(u codimL/codimM),

where f is a beta distribution with form parameters (codimM, codim_{M}L). Proof. We have to add only the remark that ∥Q_{M}(x)∥^{2} and ∥Q_{L}(x) − Q_{M}(x)∥^{2} are independent and orthogonal.

Still the beta distributions are seldom used in this case. Historically another distribution related to the beta distribution appeared first. The larger the ratio σ̂_{L}^{2} / σ̂_{M}^{2}, the less likely the hypothesis that μ is in L. This ratio is 1+ ∥ Q_{L}(x)− Q_{M}(x) ∥^{2} / ∥ Q_{M}(x) ∥^{2} .

Assume dim M = m. The hypothesis H is rejected if the ratio of orthoprojections R(x) =∥ Q_{L}(x) − Q_{M}(x) ∥^{2} / ∥ Q_{M}(x) ∥^{2} is large enough. Since Q_{L}(x) − Q_{M}(x), Q_{M}(x) are independent as random variables, and \frac{1}{m−d}σ^{−2} ∥ Q_{L}(x)−Q_{M}(x) ∥^{2} follows a χ^{2}(m−d), and since the denominator and numerator

of R(x) are independent r.v.’s, the ratio F (x) = \frac{| I | −d}{m − d}R(x), as a quotient of

two independent chi-square distributions with prescribed degrees of freedom, will follow an F distribution with bidegree (| I | −m, m − d), F_{(codimM,codim_{M} L)}.

An F_{a,b} distribution function on (0, ∞) is given by

f(x) = B(a, b)^{−1}a^{a}b^{b}x^{a−1}/(ax + b)^{a+b} , x > 0

for small bidegrees, their cumulative distributions are tabulated in textbooks. Proposition 2.2. Let α ∈ (0, 1). H is rejected at level the confidence level 1−α, if F (x) > F_{α}, where F_{α} is such that ∫_{F_{α}}^{∞} f(x)dx = α and f is the distribution function of F_{codim M,codim_{M} L}.

b. The symbol ⊥ stands for a direct sum of orthogonal subspaces. Assume that the sample space has an orthogonal decomposition R^{I} = S_{1} ⊥ S_{2} ⊥ … ⊥ S_{c}, and w.r.t. this decomposition, the components of the vector valued random variable are uncorrelated. Also assume for each j, the S_{j} - component has the covariance matrix Σ_{j} = σ_{j}^{2} Id. Such a decomposition is usually called in statistics orthogonal design. Note that for an orthogonal design, the probability distribution function obviously factors

(2.5)  L(μ, Σ; x) = Π_{j}L_{j}(μ_{j}, Σ_{j}; x_{j})

and each of the factors L_{j} is given by a (2.1) type of expression. It is then obvious that if L has a decomposition L =⊥ L_{j}, with L_{j} ⊆ S_{j}, then by a similar argument as in Rem. 2.1, one may assume that N = 1. For each j ,

(2.6)  (μ̂_{L_{j}} , σ̂_{jL_{j}}^{2}) = (pr_{L_{j}} \overline{x} , \frac{1}{|I_{j}|}∥_{j}Q_{L_{j}}(x)∥^{2}) .

In (2.6), _{j}Q_{L_{j}} represents the orthoprojection on S_{j} ∩L_{j}^{⊥}. Note that _{j}Q_{L_{j}}(x) is a vector valued ”unitary” normally distributed random variable of rank | I_{j} | −d_{j}, where d_{j} = dim L_{j} and | I_{j} |= dim S_{j} and \frac{1}{|I_{j}|−d_{j}} σ_{j}^{−2} ∥_{j}Q_{L_{j}}(x)∥^{2} follows a X_{|I_{j}|−d_{j}}^{2} distribution. The likelihood ratio associated with the hypothesis H in a., where M is also assumed to have an orthogonal splitting M =⊥ M_{j},

(2.7)  p(x) = L̂_{L}/L̂_{M} = Π_{j} (σ̂_{jM_{j}}^{2}/σ̂_{jL_{j}}^{2})^{|I_{j}|/2} .

According to Proposition 2.1, each of the factors σ̂_{jM_{j}}^{2}/σ̂_{jL_{j}}^{2} involved in the likelihood ratio, has up to a constant a beta distribution of form parameters (| I_{j} | −m_{j}, m_{j} − d_{j}). It is also obvious that these factors are independent. We set

F_{j}(x) = \frac{| I_{j} | −d_{j}}{m_{j} − d_{j}} (∥Q_{L_{j}}(x) − Q_{M_{j}}(x)∥^{2} / ∥Q_{M_{j}}(x)∥^{2}) .

There is only one situation in which the question whether H is true has a straightforward answer, that is when L and M differ only on one component of the decomposition of S. Then we may use an F -test. Proposition 2.3. Let α ∈ (0, 1). H_{j_{0}} be the hypothesis H in which L_{j} = M_{j}, if j ≠ j_{0}. H_{j_{0}} is rejected at level the confidence level 1 − α, if F_{j_{0}}(x) > F_{α}, where

F_{α} is such that ∫_{F}^{∞}f(x)dx = α, where f is the distribution of a

F_{|I_{j_{0}} |} − m_{j_{0}}, m_{j_{0}} − d_{j_{0}}.

The situation encountered in Proposition 2.3 occurs seldom in practical situations. For the general case, we need a definition. Let s be a notation for the multi index (s_{j})_{j=\overline{1,r}}. Definition 2.1. A pi-beta distribution with parameters the multi indices (a,b) and powers the multi index n is a distribution of a product r independent powers of beta distributed random variables, the j- factor, being the n_{j} power of a r.v. following a beta distribution with form parameters (a_{j}, b_{j}). Such a distribution is denoted by πβ(a, b; n).

Let b be the p.d.f. of πβ(a, b; n) distributed r.v. and let p_{a,b;n}(α) = t be a

value such that ∫_{0}^{t} b(s)ds = α. We proved the following

Theorem 2.1.. Assume L =⊥ L_{j} ⊆ M =⊥ M_{j} and α ∈ (0, 1). Then H is rejected at level the confidence level 1 − α, if p(x) < p_{a,b;n}(α).

In concrete situations, since tables of πβ_{(}a, b; n) distributions are scarce one runs into numerical problems, solvable on a computer algebra system.

3 Two and three way layouts with applications to genetherapy

In this last section, we give two examples of orthogonal designs, as an application of the theory considered in the previous section. a. The two way layout (with one observation per cell) is a classic example of orthogonal model [7]. I = R × C, where R stands for rows and C for columns; identify R^{I} with R × C real matrices, with the convention (r, c) = rc. The natural projection R × C → R, rc → r, induces the inclusion of R^{R} → R^{I}, given by (x_{r}) → (x_{rc}), where for each c, X_{rc} = x_{r}. The image of R^{R} under this monomorphism is L_{R}, the set of all R × C matrices with equal components within each row. Similarly, L_{C} is set of all R×C matrices with equal components within each column, and L_{0} = the set of R ×C matrices with equal entries. Also let L_{R+C} = L_{R} + L_{C}. Thus the main subspaces for this design are the vertices of the lattice

(3.1)

The scalar product in L = R^{I} is in matrix notation x · y = Tr(xy^{t}) which induces the following orthogonal decomposition w.r.t.:

(3.2)  R^{I} = L_{0} ⊥ (L_{R} ∩ L_{0}^{⊥}) ⊥ (L_{C} ∩ L_{0}^{⊥}) ⊥ L_{R+C}^{⊥}.

Note the following useful decompositions :

L_{R} = L_{0} ⊥ (L_{R} ∩ L_{0}^{⊥})
L_{C} = L_{0} ⊥ (L_{C} ∩ L_{0}^{⊥})
L_{R+C} = L_{0} ⊥ (L_{R} ∩ L_{0}^{⊥}) ⊥ (L_{C} ∩ L_{0}^{⊥}).

If we set I to be a matrix with all entries = 1, the projection R_{0} := P_{0}, is given by

(3.3)  R_{0}(x) = \frac{Tr(xI^{t})}{| R | | C |}I = x̄_{..}I = (x̄_{..}).

One dot as index means average of the data, w.r.t. that index.

The orhtoprojections onto the other summands in (3.2) are

R_{R} := P_{R} − P_{0}, R_{C} := P_{C} − P_{0}, R_{1} := IR^{I} − P_{R+C} = IR^{I} − (P_{R} + P_{C} − P_{0}), R_{0}x := (x̄_{..}), is the contrast vector for the level R_{R}x = (x̄_{r.} − x̄_{..}), is the contrast vector for the rows R_{C}x = (x̄_{.c} − x̄_{..}), is the contrast vector for the columns R_{1}x = (x_{rc} − x̄_{r.} − x̄_{.c} + x̄_{..}), is the contrast vector for the interaction.

Their lengths are called effects, and are key elements for the test statistics, are usually found in the so called ANOVA tables under various abbreviations, coming from their analytical formulae.

degrees of freedom  SSD
| I | − | R | − | C | +1 ∥R_{1}x∥^{2} = SSD_{1}
| R | −1  ∥R_{R}x∥^{2} = SSD_{R}
| C | −1  ∥R_{C}x∥^{2} = SSD_{C}
1  ∥R_{0}x∥^{2} = SSD_{0}

The covariance matrix has the decomposition

Σ(σ_{1}^{2}, σ_{R}^{2}, σ_{C}^{2}, σ_{0}^{2}) = Σ = σ_{1}^{2}R_{1} + σ_{R}^{2}R_{R} + σ_{C}^{2}R_{C} + σ_{0}^{2}R_{0},

where σ_{1}^{2}, σ_{R}^{2}, σ_{C}^{2}, σ_{0}^{2} are the covariance components, and various orthogonals submodels , can be obtained by setting them equal in accordance with the partial order of the lattice (3.1). They are

1. (σ_{1}^{2} = σ_{R}^{2} = σ_{C}^{2} = σ_{0}^{2})

2. (σ_{1}^{2} = σ_{R}^{2} = σ_{C}^{2}, σ_{0}^{2} > 0)

3. (σ_{1}^{2} > 0, σ_{R}^{2} = σ_{C}^{2} = σ_{0}^{2})

4. (σ_{1}^{2}, σ_{R}^{2} > 0, σ_{C}^{2} = σ_{0}^{2})

5. (σ_{1}^{2}, σ_{C}^{2} > 0, σ_{0}^{2} = σ_{R}^{2})

6. (σ_{1}^{2} = σ_{R}^{2}, σ_{C}^{2} = σ_{0}^{2} > 0)

7. (σ_{1}^{2} = σ_{C}^{2}, σ_{R}^{2} = σ_{0}^{2})

8. (σ_{1}^{2} = σ_{R}^{2}, σ_{C}^{2}, σ_{0}^{2} > 0)

9. (σ_{1}^{2} = σ_{C}^{2}, σ_{R}^{2}, σ_{0}^{2} > 0)

10. (σ_{1}^{2}, σ_{R}^{2}, σ_{C}^{2}, σ_{0}^{2} > 0)

Note the covariances are:

Σ_{rc, r′c′} = \{ \begin{matrix}| I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) & r ≠ r′, c ≠ c′ \\ | I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) + | c |^{−1} (σ_{R}^{2} − σ_{1}^{2}) & r = r′, c ≠ c′ \\ | I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) + | R |^{−1} (σ_{C}^{2} − σ_{1}^{2}) & r ≠ r′, c = c′ \\ | I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) + | c |^{−1} (σ_{R}^{2} − σ_{1}^{2}) \\ + | R |^{−1} (σ_{C}^{2} − σ_{1}^{2}) + σ_{1}^{2} & r = r′, c = c′ . \end{matrix}

The covariance structures together with the five subspaces define all the orthogonal models in this variance component design. There are 50 of them. There is a total of 13 models with unique solution for the MLE. As far as the likelihood ratio test is concerned there are even less cases.

Out of these models those that are not boxed, are uninteresting ( general linear model ), those boxed are ill posed from mathematical or random effect standpoint, those double boxed will be studied. The names of these diagrams show which of the variance component models are equal, and tell the subspace of means. For example 0R, C1, C is a distribution, with σ_{1}^{2} = σ_{C}^{2}, σ_{R}^{2} = σ_{0}^{2}, μ ∈ L_{C}.

An example of random effect problem in two indices is the following. Assume the r.v.’s X_{rc} are given by

X_{rc} = Y + Y_{r} + Y_{c} + Y_{rc},

where Y, Y_{r}, Y_{c}, Y_{rc} are independent normally distributed, Y_{r} from N(α_{r}, ω_{r}), Y_{c} from N(α_{c}, ω_{c}), Y from N(α, ω_{0}), and Y_{rc} from N(α_{rc}, ω_{1}), ω_{1} and the covariance matrix for X’s is

Σ(X_{rc}, X_{r′c′}) = \{ \begin{matrix}ω_{0} & r ≠ r′, c ≠ c′ \\ ω_{0} + ω_{r} & r = r′, c ≠ c′ \\ ω_{0} + ω_{c} & r ≠ r′, c = c′ \\ ω + ω + ω + ω & r = r′, c = c′. \end{matrix}

We identify the two way layout with the random effect model,

| I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) = ω_{0}

| I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) + | c |^{−1} (σ_{R}^{2} − σ_{1}^{2}) = ω_{0} + ω_{r}

| I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) + | R |^{−1} (σ_{C}^{2} − σ_{1}^{2}) = ω_{0} + ω_{c}

| I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) + | c |^{−1} (σ_{R}^{2} − σ_{1}^{2}) + | R |^{−1} (σ_{C}^{2} − σ_{1}^{2}) + σ_{1}^{2} = = ω_{0} + ω_{r} + ω_{c} + ω_{1}.

That is

| I |^{−1} (σ_{1}^{2} − σ_{R}^{2} − σ_{C}^{2} + σ_{0}^{2}) = ω_{0}
| c |^{−1} (σ_{R}^{2} − σ_{1}^{2}) = ω_{r}
| R |^{−1} (σ_{C}^{2} − σ_{1}^{2}) = ω_{c}
σ_{1}^{2} = ω_{1}.

Reading off this correspondence, we see that

R1  is  ω_{r} = 0
C1  is  ω_{c} = 0
0RC1 is ω_{r} = ω_{c} = ω_{0}
0RC is
| I |^{−1} (σ_{1}^{2} − σ_{RC0}^{2}) = ω_{0}
− | c |^{−1} (σ_{1}^{2} − σ_{RC0}^{2}) = ω_{r}
− | R |^{−1} (σ_{1}^{2} − σ_{RC0}^{2}) = ω_{c},

that is | I |^{−1} ω_{0} = − | c |^{−1} ω_{r} = − | R |^{−1} ω_{c}, which is quite strange;

\begin{matrix}oC, R1 & is & ω_{r} = ω_{0} = 0 \\ oR, C1 & is & ω_{c} = ω_{0} = 0. \end{matrix}

As an example one may consider the following problem, in connection with the gene therapy data [11]: assume | c | gene therapy treatments from the same category ( modified adenovirus vectors Ad5-lfu, Ad5-ifn, etc ) are injected in | R | cultures of MDA-MB-435 breast carcinoma . The data, measures the area of the tumor per culture treated after a fixed period of time. We would like to test the following simple hypothesis H: all the cultures grow in the same way independently on adenovirus type, versus K: cultures react the same way to each of the adenoviruses, but different adenoviruses yield different reactions. Their reactions to these adenoviruses are the | c | | R | observations which are assumed to be modeled in the random effect form:

X_{rc} = α + Y_{r} + α_{c} + Y_{rc},

or

X_{rc} = Y + Y_{r} + Y_{c} + Y_{rc}, with ω_{c} = ω_{0} = 0.

We have seen before that this is the mathematical model oR, C1 .

Let α̂ (respectively α̃) be value of α associated with the MLE’s of the regular model oR, C1, C ( respectively oR, C1, 0) . ˆ The evaluation spaces are

- for ξ̂: L_{0} ⊥ (L_{C} ∩ L_{0}^{⊥}) and ξ̂ = P_{C}(x) = (x̄_{.c}),

- for σ̂_{Ro}^{2} : L_{0} ⊥ (L_{R} ∩ L_{0}^{⊥}) and σ̂_{Ro}^{2} = \frac{1}{|R|}∥R_{R}x∥^{2} = \frac{1}{|R|} ∑_{rc}(x̄_{r.} − x̄_{..})^{2},

- for σ̂_{1C}^{2} and σ̂_{1C}^{2} = \frac{1}{|I|−|R|} : L_{R+C}^{⊥} and

∥R_{1}x∥^{2} = \frac{1}{| R | (| C | −1)} ∑_{rc}(x_{rc} − x̄_{r.} − x̄_{.c} + x̄_{..})^{2}.

˜ The evaluation spaces are

- for ξ̃: L_{0} and ξ̃ = P_{0}(x) = (x̄_{..}),

- for σ̃_{Ro}^{2} : L_{R} ∩ L_{0}^{⊥} and σ̃_{Ro}^{2} = \frac{1}{|R|}∥R_{R}x∥^{2} = \frac{1}{|R|} ∑_{rc}(x̄_{r.} − x̄_{..})^{2},

- for L_{R}^{⊥} : σ̃_{1C}^{2} and

σ̃_{1C}^{2} = \frac{1}{| I | − | R |}(∥R_{1}x∥^{2} + ∥R_{C}x∥^{2}) = \frac{1}{| I | − | R |} ∑_{rc}(x_{rc} − x̄_{r.} − x̄_{.c} + x̄_{..})^{2}.

(3.4)  R^{I} = L_{0} ⊥ (L_{R} ∩ L_{0}^{⊥}) ⊥ (L_{C} ∩ L_{0}^{⊥}) ⊥ L_{R+C}^{⊥}.

In the notation of section 2, R^{I} = S_{1} ⊥ S_{2}, S_{1} = L_{0} ⊥ (L_{R} ∩ L_{0}^{⊥}) = L_{R},

S_{2} = (L_{C} ∩ L_{0}^{⊥}) ⊥ L_{R+C}^{⊥}, L_{1} = L_{0}, L_{2} = 0, M_{1} = L_{0} ⊥ (L_{R} ∩ L_{0}^{⊥}), M_{2} = 0.

Since the evaluation spaces differ only on one component, we may use in this case an F -test for means.

Let us derive this test. The likelihood ratio is

(3.5)  L̃/L̂ = (∥R_{1}x∥^{2}/∥R_{1}x∥^{2} + ∥R_{C}x∥^{2})^{|I|−|R|/2} .

The hypothesis H : x ∈ L_{0} is rejected if

∥R_{C}x∥^{2}∥R_{1}x∥^{−2} > (\frac{| C | −1}{| I | − | R | − | C | +1})^{−1} F_{|C|−1, |I|−|R|−|C|+1}(1 − α).

b. Related to the gene therapy data [11], we consider now question of qualitative responses of the adenoviruses in vivo when equal amounts of identical breast cancer tumor cells are subcutaneously injected in breasts of | M | | T | female mice and then they are locally treated with | T | recombinant adenoviruses. One treatment is administered to an equal number of | M | female mice. The superficial area of the tumor is measured daily. One would like to compare the reactions to different treatments. Like in the previous example the numerical computations are deferred to a future paper. The index set is therefore a subset of I = D × T × M, so that dtm is the index of the measurement on day d of the m-th mouse tested with the treatment t.

| I |=| D | | T | | M | , D = days, T = treatments, M = mice.

The area measured is modeled by a random variable X_{dtm}. The index ordered set is (dtm) and we under the experiment the mouse is assumed to be random, as such we use a random effect model

(3.6)  X_{dtm} = α + α_{d} + α_{t} + α_{dt} + Y_{tm} + Y_{dtm},

where

Y_{tm} ∼ N(α_{tm}, ω_{T M}), ω_{T M} ≥ 0

and

Y_{dtm} ∼ N(α_{dtm}, ω_{DT M}), ω_{DT M} > 0,

are i.d.r.v. ’s (usually α_{tm} = α_{dtm} = 0). Since we have Y_{tm} and Y_{dtm} random, this is part of a three way layout

X_{dtm} = Y + Y_{d} + Y_{m} + Y_{t} + Y_{dt} + Y_{tm} + Y_{dm} + Y_{dtm}.

The correlation of the rv’s in (3.6) is

(3.7)  Σ_{dtm, d′t′m′} = \{ \begin{matrix}0 & tm ≠ t′m′ \\ ω_{T M} & tm = t′m′, d ≠ d′ \\ ω_{T M} + ω_{DT M} & dtm = d′t′m′. \end{matrix}

Like in the previous example, we would like to test the hypothesis H_{0}: the mice react the same way to all treatments, in time, that is E[X_{dtm}] = α, versus H_{1}: the mice react the same way, but in time there is a treatment effect, that is E[X_{dtm}] = α + α_{dt}. The random effect model:

X_{dtm} = α + α_{d} + α_{t} + α_{dt} + Y_{tm} + Y_{dtm}

is obviously a submodel of

X_{dtm} = Y + Y_{d} + Y_{m} + Y_{t} + Y_{dt} + Y_{tm} + Y_{dm} + Y_{dtm}.

The idea is again to identify an orthogonal variance component submodel of a three way layout with the random effect model. We use the lattice formulation in [3], but our presentation is self contained. Consider the orthogonal design

This is a sub lattice of a lattice associated with a three way layout. In the subspace language, testing problem can be translated as follows: H_{0}: The mean μ is in L_{0}; versus H_{1}: The mean μ in is L_{DT}.

The decomposition of R^{I} in terms of contrast (orthogonal) subspaces is

(3.9)  \begin{matrix}R^{I} = C_{0} ⊥ C_{D} ⊥ C_{T} ⊥ C_{DT} ⊥ C_{T M} ⊥ C_{DT M} = \\ = L_{0} ⊥ (L_{D} ∩ L_{0}^{⊥}) ⊥ (L_{T} ∩ L_{0}^{⊥}) ⊥ ((L_{D} + L_{T})^{⊥} ∩ L_{DT}) ⊥ \\ ⊥ (L_{DT}^{⊥} ∩ L_{T M}) ⊥ (L_{DT M} ∩ (L_{DT} + L_{T M})^{⊥}). \end{matrix}

Let P_{0}, P_{D}, P_{T} , P_{DT} , P_{T M}, P_{DT M} = Id be the orthoprojections on the various subspaces. They are given by

(3.10)

P_{0}(x) = (x̄_{...}) =| I |^{−1} (s_{0}), where s_{0} = ∑ x_{dtm}
P_{D}(x) = (x̄_{d..}) = (| T | | M |)^{−1}(s_{d..}), where s_{d..} = ∑_{tm} x_{dtm}
P_{T} (x) = (x̄_{.t.}) = (| D | | M |)^{−1}(s_{.t.}), where s_{.t.} = ∑_{dm} x_{dtm}
P_{DT} (x) = (x̄_{dt.}) = (| M |)^{−1}(s_{dt.}), where s_{dt.} = ∑_{m} x_{dtm}
P_{T M}(x) = (x̄_{.tm}) = (| T |)^{−1}(s_{.tm}), where s_{.tm} = ∑_{d} x_{dtm}
P_{DT M}(x) = x = (x_{dtm}).

The contrast vectors (3.11) R_{0}(x) = (x̄_{...}) R_{D}(x) = P_{D}(x) − P_{0}(x) = (x̄_{d..} − x̄_{...}) R_{T} (x) = P_{T} (x) − P_{0}(x) = (x̄_{.t.} − x̄_{...}) R_{DT} (x) = P_{DT} (x) − P_{D}(x) − P_{T} (x) + P_{0}(x) = (x̄_{dt.} − x̄_{d..} − x̄_{.t.} + x̄_{...}) L_{T M} = (L_{DT} ∩ L_{T M}) ⊥ (L_{DT}^{⊥} ∩ L_{T M}) = L_{T} ⊥ (L_{DT}^{⊥} ∩ L_{T M}) R_{T M}(x) = P_{T M}(x) − P_{T} (x) = (x̄_{.tm} − x̄_{.t.}) R_{DT M}(x) = (Id(x) + P_{T} (x) − P_{T M}(x) − P_{DT} (x)) = (x_{dtm} − x̄_{.tm} − x̄_{dt.} + x̄_{.t.}). In the three way layout the covariance matrix has 6 variance components (3.12)

∑ = σ_{0}^{2}R_{0} + σ_{D}^{2}R_{D} + σ_{T}^{2} R_{T} + σ_{DT}^{2} R_{DT} + σ_{T M}^{2}R_{T M} + σ_{DT M}^{2}R_{DT M} =
= (σ_{0}^{2} − σ_{D}^{2} − σ_{T}^{2} + σ_{DT}^{2} )P_{0} + (σ_{D}^{2} − σ_{DT}^{2} )P_{D}+
+(σ_{T}^{2} − σ_{DT}^{2} − σ_{T M}^{2} + σ_{DT M}^{2})P_{T} + (σ_{DT}^{2} − σ_{DT M}^{2})P_{DT} +
+(σ_{T M}^{2} − σ_{DT M}^{2})P_{T M} + σ_{DT M}^{2}Id.

If one identifies the matrix ∑ in (3.7) with a linear endomorphism of R^{I},

(3.13)  ∑ =| D | ω_{T M}P_{T M} + ω_{DT M}Id.

The models are the same if

(3.14)  \begin{matrix}σ_{0}^{2} − σ_{D}^{2} − σ_{T}^{2} + σ_{DT}^{2} = 0, σ_{D}^{2} − σ_{DT}^{2} = 0 \\ σ_{T}^{2} − σ_{DT}^{2} − σ_{T M}^{2} + σ_{DT M}^{2} = 0, σ_{DT}^{2} − σ_{DT M}^{2} = 0 \\ σ_{T M}^{2} − σ_{DT M}^{2} =| D | ω_{T M}, σ_{DT M}^{2} = ω_{DT M}. \end{matrix}

or

| D |^{−1} (σ_{T M}^{2} − σ_{DT M}^{2}) = ω_{T M}, σ_{DT M}^{2} = ω_{DT M}.

One can treat the random effect model as a 2 variance components model,

\begin{matrix}σ_{0}^{2} = σ_{T}^{2} = σ_{T M}^{2} = ω_{DT M}+ | D | ω_{T M} \\ σ_{D}^{2} = σ_{DT}^{2} = σ_{DT M}^{2} = ω_{DT M}. \end{matrix}

In the notation of section 2, R^{I} = S_{1} ⊥ S_{2}, S_{1} = C_{0} ⊥ C_{T} ⊥ C_{T M},

S_{2} = C_{D} ⊥ C_{DT} ⊥ C_{DT M}
L_{1} = C_{0}, L_{2} = 0,
M_{1} = C_{0} ⊥ C_{T} , M_{2} = C_{D} ⊥ C_{DT} .

As such although the question in the examples in 3.a. and 3.b. are the same, using the results in the previous section, we see that the test statistics has to be based on a product of powers two independent beta’s. Let h be the distribution function of a product of two independent distributions, a beta with form parameters (| T | | M | − | T |, | T | −1) by the | D | −1 power of a beta with form parameters

(| D | | T | | M | − | T | | M | − | D | | T | + | T |, | D || T | − | T |).

Theorem 3.1. H_{0} is rejected at level of confidence 1 − α, if

(∥R_{DT M}x∥^{2} / (∥R_{D}x∥^{2} + ∥R_{DT} x∥^{2} + ∥R_{DT M}x∥^{2}))^{(|D|−1)} ·
· (∥R_{T M}x∥^{2} / (∥R_{T} x∥^{2} + ∥R_{T M}x∥^{2})) < h_{α},

where

∫_{0}^{h_{α}} h(x)dx = α.

Remark 3.1. Under this simple hypothesis testing problem, the main obstacles are numerical rather then conceptual. The key steps in solving such a problem are not only hidden in tables for the cumulative distributions of products of powers of beta distributed r.v.’ s, but equally for concrete questions the numerical work stays behind the computation of various orthogonal components of the projections of the data x onto the corresponding summands L_{j}.

Acknowledgement. I am very grateful to Steen Andersson, for sharing his insight on multivariate normal statistics and orthogonal models and to Victor Goodman for useful conversations.

References

[1] S. Amari, Differential geometry in statistical inference, Bull. Inst. Internat. Stat., 52 (1987), no 2, 321-338.

[2] T.W. Anderson, An introduction to multivariate statistical analysis, 2-nd edition, Wiley Series in Mathematical Probability and Statistics, 1984.

[3] S.A. Andersson, The lattice structure of orthogonal linear models and orthogonal variance component models, Scand. J. Statist. 17 (1990), 287-319.

[4] S.A. Andersson, J. Madsen, Symmetry and lattice conditional independence in a multivariate normal distribution, to appear in Ann. of Statistics.

[5] O.E. Barndorff - Nielsen, P.E. Jupp, Statistical inference and differential geometry - some recent developments, Multivariate analysis: future direction (Barcelona, 1992), 385-396, North Holland Ser. in Statist. Prob., Amsterdam, 1993.

[6] B.J. Coventry, S.C.Weeks, S.E. Heckford, P.J. Sykes, J.M. Skinner, Lack of IL-2 cytokine expression despite IL-2 messanger RNA transcription in TIL in primary human breast cancer, J-Immunol. 1996, May 1, 156 (9): 3486-3492.

[7] I. Guttman, Linear models: an introduction, Wiley Ser. in Prob. and Math. Stat., 1982.

[8] D. Morrison, Multivariate Statistical Methods, McGraw Hill Ser. in Prob. and Stat., Third Ed., 1990.

[9] P.F. Robbins, M.Li.Y.F. El-Gamil, Y. Kawakami, D. Loftus, E. Appella, S.A. Rosenberg, Breast cancers recognized by TIL, J. Exp. Med., 1996 Mar. 1; 183 (3); 1185-1192.

[10] H. Scheffé, The analysis of variance, J. Wiley and Sons, Inc., 1959.

[11] M.W. Taylor, J.F. Zhang, C. Hu, Y. Geng, J. Selm, S.B. Klein, A. Orazi, Treatment of a human breast cancer xenograft with an adenovirus vector containing an interferon gene results in rapid regression due to viral oncolysis and gene therapy, Proc. Natl. Acad. Sci. USA, vol. 93, 4513-4518, April 1996.

Indiana University Department of Mathematics Bloomington, IN 47405-5701, USA e-mail: vpatrang@indiana.edu