# Statistics: Ellipses

### From MathWiki

Table of contents |

## Ellipses

For a positive definite matrix we use to denote the ellipse .

Note that every non-negative definite matrix can be factored as . The matrix can be selected so that it is square. It will be non-singular if and only if is non-singular. A definition of an ellipse that can be used for all non-negative definite matrices and that corresponds to the previous definition in the case of positive-definite matrices is

where is a unit sphere of conformable dimension.

A convenient notation in terms of *W* is

Note that under a linear transformation given by a matrix
*L*
the image of the ellipse is:

### Conjugate axes

The columns of the matrix
form a set of *conjugate axes* of the ellipse. Each vector
lies on the ellipse and the tangent space at that point is parallel to the span of all the other column vectors of

For this is illustrated in the following figure

in which

and

If is any other factorization of
then the columns of
have the same properties as the columns of
.

Consider the inner product space with inner product matrix and inner product

- .

Since we see that and are orthogonal unit vectors (in fact, an orthonormal basis) in this inner product:

Particular factorizations yield interesting sets of conjugate axes. In the figure on the right

the axes for the Choleski factorization where
is lower triangular are shown in green and the axes for the *principal component factorization* (where
and
is the spectral decomposition of
), are shown in magenta.
The figures were produced with a script in R.

## Hypothesis and error ellipses

Let be a positive definite matrix and be a non-negative definite matrix.

There exists a matrix and a non-negative diagonal matrix such that:

and

- .

The diagonal elements of , , are the eigenvalues of relative to . An R function to compute is given below.

In this discussion and can represent the and an matrices in multivariate linear regression both divided by a common factor, perhaps the degrees of freedom for so that is an unbiased estimate of the residual variance matrix .

## Expected value of SSH

We can find expressions for in terms of the data and find expected values for under fixed effects or under random effects models.

### Univariate response

Consider first the expected value of with a univariate response and a one-way Anova model with random group means.

The following formula is used repeatedly to find the expected value of a quadratic form. If is a random vector with and , and is a conformable matrix then

With a univariate response and groups we can use an indicator design matrix and write:

where

and

and is independent of .

Thus

- and

Consequently

= = = = = =

where is the group-size weighted mean of group sizes. With equal groups, and

Thus

= = =

### Multivariate response

Consider one-way MANOVA with random group means:

and

- .

With we have

The analogous results are:

and

Note that

and that the group-size weighted average of these variances is:

= = = with equal groups

The following table shows expectations of combinations of and of the form :

1 0 0 1 0

0 with equal groups

with equal groups 0

## Definition of SSH with one-way MANOVA

Consider one-way MANOVA with groups, each with observations. Let and let

be the design matrix.

Let be the residual projection matrix defined by

In the case of one-way MANOVA there is an easy way to derive an expression for SSH:

### Definition of SSH with a general linear hypothesis

Given a univariate linear model

with *X*_{i} of dimension , the test for is based on

where

and

We can show that

where is the partial SSCP matrix for regressing on and is the matrix of residuals of regressing on .

Thus is the sum of squared deviations (assuming that spans the intercept) of 'y hat' in the added-variable plot for . It is the variability in 'explained' by beyond that 'explained' by . is the sum of squared deviations of residuals in the same plot.

Note that

The matrix of residuals, , is obtained by regression on and, thus, has degrees of freedom. We could therefore think of

with

as an estimate of the variance of contribution of to the variability of of .

### Computing eigenvalues

To understand how eigenvalues are computed in car:::print.Anova.mlm, let

Since

contains the eigenvalues of .

Thus the eigenvalues and eigenvectors of relative can be obtained with

> rel.eval ( H, E ) > rel.evac ( H, E )

where

> rel.eval <- function ( H, E ) { # This code is adapted from car:::print.Anova.mlm Re( eigen( qr.coef( qr(E), H ), only.values = TRUE )$values ) } > rel.evec <- function ( H, E ) { Re(eigen( t(qr.coef( qr(E), H )))$vectors) }

The eigenvectors are not, in general, a factor of . To obtain a factor of one can use:

> rel.fac <- function( H, E ) { Bp <- fac(E) B.qr <- qr(t(Bp)) H.star <- qr.coef( B.qr, t( qr.coef( B.qr, H ) )) # B^{-1} H B'^{-1} ev <- eigen( H.star , symmetric = T) ret <- t(ev$vectors) %*% Bp attr(ret,"lambda") <- ev$values ret }

Note that the numeric properties of this function could be improved.

## Visualizing significance

If we plot the ellipses and , the ellipse will protrude beyond along the direction , where iff .

### Roy's largest root test

The F-test based on Roy's largest root uses the approximation ( e.g., see SAS documentation for multivariate tests (*http://post.queensu.ca:8080/SASDoc/getDoc/en/statug.hlp/introreg_sect21.htm*) )

with degrees of freedom and where

where is the number of degrees of freedom for the hypothesis and

where is the number of degrees of freedom for the error matrix .

The critical value for for an level test, can be computed as

> lambda.crit <- function( alpha, p , q, nu.E) { d <- max( p, q) nu <- nu.E - d + q (d/nu) * qf( alpha, d , nu, lower.tail = F) }

Thus, the ellipse will have the property that it protrudes beyond if and only if the Roy test for the hypothesis whose matrix is is significant at level .

### Hotelling-Lawley trace

The Hotelling-Lawley trace statistic is based on the sum of the eigenvalues of . This suggests another criterion for judging visual significance based on the average eigenvalue, .

One approximate F-test for the HLT is

where , are the degrees of freedom for the numerator and denominator, and , .

Thus, the critical value for can be computed as

> HLT.crit <- function ( alpha, p, dfh, dfe) { s <- min(p, dfh) m <- (abs(p-dfh)-1)/2 n <- (dfe-p-1)/2 df1 <- 2*m + s + 1 df2 <- 2*(s*n +1) s * (df1/df2) * qf(alpha, df1, df2, lower.tail=FALSE) }

## When p > 3

When the number of response variables involved in a hypothesis is greater than the number of dimensions plotted, one can use the vectors of the matrix to construct linear combinations of the response variables that will show maximal significance of the hypothesis. These linear combinations depend on both and and thus could not be used to display results for many hypotheses simultaneously.

By projecting into the space corresponding to the three largest relative eigenvalues, we scale the projected matrix, , so its ellipse will protruded outside the ellipse of the projected matrix, , if and only if the overall Roy test is significant.

A convenient choice of is

where

Then

- and

**UNFINISHED**

### Note on computing matrix factors

A Choleski root, above, can be obtained in R with

> chol(W)

when is positive-definite.

The singular value decomposition also works for singular and ill-conditioned matrices and can be obtained with

> fac(W)

where

> fac <- function(W) { # assumes W is non-negative definite # returns the 'right' of W, D^(1/2) Gamma', where W = Gamma D Gamma' xx <- svd(W, nu=0) t(xx$v) * sqrt( pmax( xx$d,0) ) }

Thus, an ellipse can be generated in R with:

> ell( mu, W )

where

> ell <- function( mu, W ) { p <- nrow(W) t( mu + t( sphere(p) %*% fac(W) ) ) }where

sphere(p)is unit sphere of the right dimension. Of course, an implementation of this function would need to attend to many additional details.

To the extent that the sphere is isotropic, the generated ellipse is independent of the choice of factor for *W*.

If the sphere object is augmented with diameters along usual axes, then the image of the diameters will be columns of which will form a set of conjugate vectors of the ellipse. A set of conjugate vectors of an ellipse has the property that each vector is on the surface of the ellipse and any subset of vectors spans a space that is parallel to the tangent spaces at any of the other vectors. The set of linear combinations of the vectors with coefficients in is a parallelogram that contains the ellipse and whose faces are tangent spaces to the ellipse.

With two symmetric matrices, and , with positive definite, we can find a factor of whose columns will also be conjugate *directions* for . The following function will produce an 'eigenvector' based right factor of if called with a single argument. With two matrix arguments, it returns a *right* factor of whose *rows* are also conjugate directions for . With two matrix arguments, the first argument must be positive definite.

fac <- function( E, # non-neg. def. matrix to factor -- should be pos. def. if H is not missing H, # if not missing, factor of E will consist of axes conjugate relative to H also which = 1 # which = 2 returns a factor of H ) { xx <- eigen(E,symmetric = T) Efac <- t(xx$vector) * sqrt( pmax( Re( xx$values ), 0 )) # right factor of E if ( missing(H) ) return(Efac) Efac.qr <- qr( t( Efac ) ) H.star <- qr.coef( Efac.qr, t( qr.coef( Efac.qr, H))) # Efac^{-1} H Efac'^{-1} ev <- eigen( H.star, symmetric = T ) ret <- t(ev$vectors) %*% Efac if ( which == 2 ) ret <- ret * ev$values attr( ret, "lambda" ) <- ev$values ret }

## Canonical transformations

See Statistics: Canonical transformations for HEplots