Multivariate statistics for Data science Notes

Multivariate Methods

$y_1,...,y_p$ are possibly correlated random variables with means ${\mu }_1,...,{\mu }_p$

$$y=\left(\begin{array}{c}{y}_1\\ .\\ y_p\end{array}\right)$$

$$E\left(y\right)=\left(\begin{array}{c}{\mu }_1\\ .\\ {\mu }_p\end{array}\right)$$

Let ${\sigma }_{ij}=cov(y_i,y_j)$ for $i,j=1,\dots ,p$

$$\Sigma =\left({\sigma }_{ij}\right)=\left(\begin{array}{cccc}{\sigma }_{11}& {\sigma }_{22}& ...& {\sigma }_{1p}\\ {\sigma }_{21}& {\sigma }_{22}& ...& {\sigma }_{2p}\\ .& .& .& .\\ {\sigma }_{p1}& {\sigma }_{p2}& ...& {\sigma }_{pp}\end{array}\right)$$

where $\Sigma$ (symmetric) is the variance-covariance or dispersion matrix

Let $u_{p\times 1}$ and $v_{q\times 1}$ be random vectors with means ${\mu }_u$ and ${\mu }_v$ . Then

$${\Sigma }_{uv}=cov\left(u,v\right)=E\left[\right(u-{\mu }_u\left)\right(v-{\mu }_v{)}^{\prime }]$$

in which ${\Sigma }_{uv}\ne {\Sigma }_{vu}$ and ${\Sigma }_{uv}={\Sigma }_{vu}^{\prime }$

Properties of Covariance Matrices

1. Symmetric ${\Sigma }^{\prime }=\Sigma$
2. Non-negative definite $a^{\prime }\Sigma a\ge 0$ for any $a\in R^p$, which is equivalent to eigenvalues of $\Sigma$, ${\lambda }_1\ge {\lambda }_2\ge ...\ge {\lambda }_p\ge 0$
3. $\left|\Sigma \right|={\lambda }_1{\lambda }_2...{\lambda }_p\ge 0$ (generalized variance) (the bigger this number is, the more variation there is
4. $trace\left(\Sigma \right)=tr\left(\Sigma \right)={\lambda }_1+...+{\lambda }_p={\sigma }_{11}+...+{\sigma }_{pp}=$ sum of variance (total variance)

Note:

• $\Sigma$ is typically required to be positive definite, which means all eigenvalues are positive, and $\Sigma$ has an inverse ${\Sigma }^{-1}$ such that ${\Sigma }^{-1}\Sigma =I_{p\times p}={\Sigma \Sigma }^{-1}$

Correlation Matrices

$${\rho }_{ij}=\frac{{\sigma }_{ij}}{\sqrt{{\sigma }_{ii}{\sigma }_{jj}}}$$

$$R=\left(\begin{array}{cccc}{\rho }_{11}& {\rho }_{12}& ...& {\rho }_{1p}\\ {\rho }_{21}& {\rho }_{22}& ...& {\rho }_{2p}\\ .& .& .& .\\ {\rho }_{p1}& {\rho }_{p2}& ...& {\rho }_{pp}\end{array}\right)$$

where ${\rho }_{ij}$ is the correlation, and ${\rho }_{ii}=1$ for all i

Alternatively,

$$R=\left[diag\right(\Sigma \left){]}^{-1/2}\Sigma \right[diag\left(\Sigma \right){]}^{-1/2}$$

where $diag\left(\Sigma \right)$ is the matrix which has the ${\sigma }_{ii}$’s on the diagonal and 0’s elsewhere

and $A^{1/2}$ (the square root of a symmetric matrix) is a symmetric matrix such as $A=A^{1/2}{A}^{1/2}$

Equalities

Let

• $x$ and $y$ be random vectors with means ${\mu }_x$ and ${\mu }_y$ and variance -variance matrices ${\Sigma }_x$ and ${\Sigma }_y$.

• $A$ and $B$ be matrices of constants and $c$ and $d$ be vectors of constants

Then

• $E\left(\mathrm{Ay}+c\right)=A{\mu }_y+c$

• $var\left(\mathrm{Ay}+c\right)=Avar\left(y\right)A^{\prime }={A{\Sigma }_{y}A}^{\prime }$

• $cov\left(\mathrm{Ay}+c,\mathrm{By}+d\right)={A{\Sigma }_{y}B}^{\prime }$

• $E\left(\mathrm{Ay}+\mathrm{Bx}+c\right)=A{\mu }_y+B{\mu }_x+c$

• $var\left(\mathrm{Ay}+\mathrm{Bx}+c\right)=A{\Sigma }_y{A}^{\prime }+B{\Sigma }_x{B}^{\prime }+A{\Sigma }_{\mathrm{yx}}{B}^{\prime }+B{\Sigma }_{\mathrm{yx}}^{\prime }{A}^{\prime }$

Multivariate Normal Distribution

Let $y$ be a multivariate normal (MVN) random variable with mean $\mu$ and variance $\Sigma$. Then the density of $y$ is

$$f\left(y\right)=\frac{1}{\left(2\pi {)}^{p/2}\right|\Sigma {|}^{1/2}}\exp (-\frac{1}{2}(y-\mu {)}^{\prime }{\Sigma }^{-1}(y-\mu ))$$

$y\sim N_p(\mu ,\Sigma )$

Properties of MVN

• Let $A_{r\times p}$ be a fixed matrix. Then $\mathrm{Ay}\sim N_r\left(A\mu ,A\Sigma A^{\prime }\right)$ . $r\le p$ and all rows of $A$ must be linearly independent to guarantee that ${A\Sigma A}^{\prime }$ is non-singular.

• Let $G$ be a matrix such that ${\Sigma }^{-1}={\mathrm{GG}}^{\prime }$. Then $G^{\prime }y\sim N_p\left(G^{\prime }\mu ,I\right)$ and $G^{\prime }(y-\mu )\sim N_p(0,I)$

• Any fixed linear combination of $y_1,...,y_p$ (say $c^{\prime }y$) follows $c^{\prime }y\sim N_1\left(c^{\prime }\mu ,c^{\prime }\Sigma c\right)$

• Define a partition, $[y_1^{\prime },y_2^{\prime }{]}^{\prime }$ where

  • $y_1$ is $p_1\times 1$

  • $y_2$ is $p_2\times 1$,

  • $p_1+p_2=p$

  • $p_1,p_2\ge 1$ Then

$$\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)\sim N(\left(\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right),\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right))$$

• The marginal distributions of $y_1$ and $y_2$ are $y_1\sim N_{p1}\left({\mu }_1,{\Sigma }_{11}\right)$ and $y_2\sim N_{p2}\left({\mu }_2,{\Sigma }_{22}\right)$

• Individual components $y_1,...,y_p$ are all normally distributed $y_i\sim N_1({\mu }_i,{\sigma }_{ii})$

• The conditional distribution of $y_1$ and $y_2$ is normal

  • $y_1|y_2\sim N_{p1}({\mu }_1+{\Sigma }_{12}{\Sigma }_{22}^{-1}(y_2-{\mu }_2),{\Sigma }_{11}-{\Sigma }_{12}{\Sigma }_{22}^{-1}{\sigma }_{21})$

    • In this formula, we see if we know (have info about) $y_2$, we can re-weight $y_1$ ’s mean, and the variance is reduced because we know more about $y_1$ because we know $y_2$

  • which is analogous to $y_2|y_1$. And $y_1$ and $y_2$ are independently distrusted only if ${\Sigma }_{12}=0$

• If $y\sim N\left(\mu ,\Sigma \right)$ and $\Sigma$ is positive definite, then $(y-\mu {)}^{\prime }{\Sigma }^{-1}(y-\mu )\sim {\chi }_{\left(p\right)}^2$

• If $y_i$ are independent $N_p({\mu }_i,{\Sigma }_i)$ random variables, then for fixed matrices $A_{i(m\times p)}$, ${\sum }_{i=1}^k{A}_i{y}_i\sim N_m({\sum }_{i=1}^k{A}_i{\mu }_i,{\sum }_{i=1}^k{A}_i{\Sigma }_i{A}_i)$

Multiple Regression

$$\left(\begin{array}{c}Y\\ x\end{array}\right)\sim N_{p+1}(\left[\begin{array}{c}{\mu }_y\\ {\mu }_x\end{array}\right],\left[\begin{array}{cc}{\sigma }_Y^2& {\Sigma }_{yx}\\ {\Sigma }_{yx}& {\Sigma }_{xx}\end{array}\right])$$

The conditional distribution of Y given x follows a univariate normal distribution with

$$\begin{array}{cc}E\left(Y\right|x)& ={\mu }_y+{\Sigma }_{yx}{\Sigma }_{xx}^{-1}(x-{\mu }_x)\\ & ={\mu }_y-{\Sigma }_{yx}{\Sigma }_{xx}^{-1}{\mu }_x+{\Sigma }_{yx}{\Sigma }_{xx}^{-1}x\\ & ={\beta }_0+{\beta }^{\prime }x\end{array}$$

where $\beta =({\beta }_1,...,{\beta }_p{)}^{\prime }={\Sigma }_{xx}^{-1}{\Sigma }_{yx}^{\prime }$ (e.g., analogous to $(x^{\prime }x{)}^{-1}{x}^{\prime }y$ but not the same if we consider $Y_i$ and $x_i$, $i=1,..,n$ and use the empirical covariance formula: $var\left(Y\right|x)={\sigma }_Y^2-{\Sigma }_{\mathrm{yx}}{\Sigma }_{\mathrm{xx}}^{-1}{\Sigma }_{\mathrm{yx}}^{\prime }$)

Samples from Multivariate Normal Populations

A random sample of size n, $y_1,..,y_n$ from $N_p(\mu ,\Sigma )$. Then

• Since $y_1,...,y_n$ are iid, their sample mean, $\overline{y}={\sum }_{i=1}^n{y}_i/n\sim N_p(\mu ,\Sigma /n)$. that is, $\overline{y}$ is an unbiased estimator of $\mu$

• The $p\times p$ sample variance-covariance matrix, $S$ is $S=\frac{1}{n-1}{\sum }_{i=1}^n(y_i-\overline{y})(y_i-\overline{y}{)}^{\prime }=\frac{1}{n-1}({\sum }_{i=1}^n{y}_i{y}_i^{\prime }-n\overline{y}{\overline{y}}^{\prime })$

  • where $S$ is symmetric, unbiased estimator of $\Sigma$ and has $p(p+1)/2$ random variables.

• $(n-1)S\sim W_p(n-1,\Sigma )$ is a Wishart distribution with n-1 degrees of freedom and expectation $(n-1)\Sigma$. The Wishart distribution is a multivariate extension of the Chi-squared distribution.

• $\overline{y}$ and $S$ are independent

• $\overline{y}$ and $S$ are sufficient statistics. (All of the info in the data about $\mu$ and $\Sigma$ is contained in $\overline{y}$ and $S$ , regardless of sample size).

Large Sample Properties

$y_1,...,y_n$ are a random sample from some population with mean $\mu$ and variance-covariance matrix $\Sigma$

• $\overline{y}$ is a consistent estimator for $\mu$

• $S$ is a consistent estimator for $\Sigma$

• Multivariate Central Limit Theorem: Similar to the univariate case, $\sqrt{n}(\overline{y}-\mu )\dot{\sim }{N}_p\left(0,\Sigma \right)$ where n is large relative to p ($n\ge 25p$), which is equivalent to $\overline{y}\dot{\sim }{N}_p(\mu ,\Sigma /n)$

• Wald’s Theorem: $n(\overline{y}-\mu {)}^{\prime }{S}^{-1}(\overline{y}-\mu )$ when n is large relative to p.

Maximum Likelihood Estimation for MVN

Suppose iid $y_1,...y_n\sim N_p(\mu ,\Sigma )$, the likelihood function for the data is

$$\begin{array}{cc}L(\mu ,\Sigma )& =\prod _{j=1}^n\left(\frac{1}{\left(2\pi {)}^{p/2}\right|\Sigma {|}^{1/2}}\exp \right(-\frac{1}{2}(y_j-\mu {)}^{\prime }{\Sigma }^{-1})(y_j-\mu ))\\ & =\frac{1}{\left(2\pi {)}^{np/2}\right|\Sigma {|}^{n/2}}\exp (-\frac{1}{2}\sum _{j=1}^n(y_j-\mu {)}^{\prime }{\Sigma }^{-1}\left)\right(y_j-\mu )\end{array}$$

Then, the MLEs are

$$\hat{\mu }=\overline{y}$$

$$\hat{\Sigma }=\frac{n-1}{n}S$$

using derivatives of the log of the likelihood function with respect to $\mu$ and $\Sigma$

Properties of MLEs

• Invariance: If $\hat{\theta }$ is the MLE of $\theta$, then the MLE of $h\left(\theta \right)$ is $h\left(\hat{\theta }\right)$ for any function h(.)

• Consistency: MLEs are consistent estimators, but they are usually biased

• Efficiency: MLEs are efficient estimators (no other estimator has a smaller variance for large samples)

• Asymptotic normality: Suppose that ${\hat{\theta }}_n$ is the MLE for $\theta$ based upon n independent observations. Then ${\hat{\theta }}_n\dot{\sim }N(\theta ,H^{-1})$

  • $H$ is the Fisher Information Matrix, which contains the expected values of the second partial derivatives fo the log-likelihood function. the (i,j)th element of $H$ is $-E\left(\frac{{\partial }^{2}l\left(\theta \right)}{\partial {\theta }_i\partial {\theta }_j}\right)$

  • we can estimate $H$ by finding the form determined above, and evaluate it at $\theta ={\hat{\theta }}_n$

• Likelihood ratio testing: for some null hypothesis, $H_0$ we can form a likelihood ratio test

  • The statistic is: $\Lambda =\frac{{max}_{H_0}l(\mu ,\Sigma |Y)}{maxl(\mu ,\Sigma |Y)}$

  • For large n, $-2\log \Lambda \sim {\chi }_{\left(v\right)}^2$ where v is the number of parameters in the unrestricted space minus the number of parameters under $H_0$

Test of Multivariate Normality

• Check univariate normality for each trait (X) separately

  • Can check [Normality Assessment]

  • The good thing is that if any of the univariate trait is not normal, then the joint distribution is not normal (see again [m]). If a joint multivariate distribution is normal, then the marginal distribution has to be normal.

  • However, marginal normality of all traits does not imply joint MVN

  • Easily rule out multivariate normality, but not easy to prove it

• Mardia’s tests for multivariate normality

  • Multivariate skewness is$${\beta }_{1,p}=E\left[\right(y-\mu {)}^{\prime }{\Sigma }^{-1}(x-\mu ){]}^3$$

  • where $x$ and $y$ are independent, but have the same distribution (note: $\beta$ here is not regression coefficient)

  • Multivariate kurtosis is defined as

  • $${\beta }_{2,p}-E\left[\right(y-\mu {)}^{\prime }{\Sigma }^{-1}(x-\mu ){]}^2$$

  • For the MVN distribution, we have ${\beta }_{1,p}=0$ and ${\beta }_{2,p}=p(p+2)$

  • For a sample of size n, we can estimate

    $${\hat{\beta }}_{1,p}=\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n{g}_{ij}^2$$

    $${\hat{\beta }}_{2,p}=\frac{1}{n}\sum _{i=1}^n{g}_{ii}^2$$

    • where $g_{ij}=(y_i-\overline{y}{)}^{\prime }{S}^{-1}(y_j-\overline{y})$. Note: $g_{ii}=d_i^2$ where $d_i^2$ is the Mahalanobis distance

  • [@MARDIA_1970] shows for large n

    $${\kappa }_1=\frac{n{\hat{\beta }}_{1,p}}{6}\dot{\sim }{\chi }_{p(p+1)(p+2)/6}^2$$

    $${\kappa }_2=\frac{{\hat{\beta }}_{2,p}-p(p+2)}{\sqrt{8p(p+2)/n}}\sim N(0,1)$$

    • Hence, we can use ${\kappa }_1$ and ${\kappa }_2$ to test the null hypothesis of MVN.

    • When the data are non-normal, normal theory tests on the mean are sensitive to ${\beta }_{1,p}$ , while tests on the covariance are sensitive to ${\beta }_{2,p}$

• Alternatively, Doornik-Hansen test for multivariate normality [@doornik2008]

• Chi-square Q-Q plot

  • Let $y_i,i=1,...,n$ be a random sample sample from $N_p(\mu ,\Sigma )$

  • Then $z_i={\Sigma }^{-1/2}(y_i-\mu ),i=1,...,n$ are iid $N_p(0,I)$. Thus, $d_i^2=z_i^{\prime }{z}_i\sim {\chi }_p^2,i=1,...,n$

  • plot the ordered $d_i^2$ values against the qualities of the ${\chi }_p^2$ distribution. When normality holds, the plot should approximately resemble a straight lien passing through the origin at a 45 degree

  • it requires large sample size (i.e., sensitive to sample size). Even if we generate data from a MVN, the tail of the Chi-square Q-Q plot can still be out of line.

• If the data are not normal, we can

  • ignore it

  • use nonparametric methods

  • use models based upon an approximate distribution (e.g., GLMM)

  • try performing a transformation

Mean Vector Inference

In the univariate normal distribution, we test $H_0:\mu ={\mu }_0$ by using

$$T=\frac{\overline{y}-{\mu }_0}{s/\sqrt{n}}\sim t_{n-1}$$

under the null hypothesis. And reject the null if $\left|T\right|$ is large relative to $t_{(1-\alpha /2,n-1)}$ because it means that seeing a value as large as what we observed is rare if the null is true

Equivalently,

$$T^2=\frac{(\overline{y}-{\mu }_0{)}^2}{s^2/n}=n(\overline{y}-{\mu }_0)\left(s^2{)}^{-1}\right(\overline{y}-{\mu }_0)\sim f_{(1,n-1)}$$

Natural Multivariate Generalization

$$H_0:\mu ={\mu }_0{H}_a:\mu \ne {\mu }_0$$

Define Hotelling’s $T^2$ by

$$T^2=n(\overline{y}-{\mu }_0{)}^{\prime }{S}^{-1}(\overline{y}-{\mu }_0)$$

which can be viewed as a generalized distance between $\overline{y}$ and ${\mu }_0$

Under the assumption of normality,

$$F=\frac{n-p}{(n-1)p}{T}^2\sim f_{(p,n-p)}$$

and reject the null hypothesis when $F>f_{(1-\alpha ,p,n-p)}$

• The $T^2$ test is invariant to changes in measurement units.

  • If $z=\mathrm{Cy}+d$ where $C$ and $d$ do not depend on $y$, then $T^2\left(z\right)-T^2\left(y\right)$

• The $T^2$ test can be derived as a likelihood ratio test of $H_0:\mu ={\mu }_0$

Confidence Intervals

Confidence Region

An “exact” $100(1-\alpha )\%$ confidence region for $\mu$ is the set of all vectors, $v$, which are “close enough” to the observed mean vector, $\overline{y}$ to satisfy

$$n(\overline{y}-{\mu }_0{)}^{\prime }{S}^{-1}(\overline{y}-{\mu }_0)\le \frac{(n-1)p}{n-p}{f}_{(1-\alpha ,p,n-p)}$$

• $v$ are just the mean vectors that are not rejected by the $T^2$ test when $\overline{y}$ is observed.

In case that you have 2 parameters, the confidence region is a “hyper-ellipsoid”.

In this region, it consists of all ${\mu }_0$ vectors for which the $T^2$ test would not reject $H_0$ at significance level $\alpha$

Even though the confidence region better assesses the joint knowledge concerning plausible values of $\mu$ , people typically include confidence statement about the individual component means. We’d like all of the separate confidence statements to hold simultaneously with a specified high probability. Simultaneous confidence intervals: intervals against any statement being incorrect

Simultaneous Confidence Statements

• Intervals based on a rectangular confidence region by projecting the previous region onto the coordinate axes:

$${\overline{y}}_i\pm \sqrt{\frac{(n-1)p}{n-p}{f}_{(1-\alpha ,p,n-p)}\frac{s_{ii}}{n}}$$

for all $i=1,..,p$

which implied confidence region is conservative; it has at least $100(1-\alpha )\%$

Generally, simultaneous $100(1-\alpha )\%$ confidence intervals for all linear combinations , $a$ of the elements of the mean vector are given by

$$a^{\prime }\overline{y}\pm \sqrt{\frac{(n-1)p}{n-p}{f}_{(1-\alpha ,p,n-p)}\frac{a^{\prime }\mathrm{Sa}}{n}}$$

• works for any arbitrary linear combination $a^{\prime }\mu =a_1{\mu }_1+...+a_p{\mu }_p$, which is a projection onto the axis in the direction of $a$

• These intervals have the property that the probability that at least one such interval does not contain the appropriate $a^{\prime }\mu$ is no more than $\alpha$

• These types of intervals can be used for “data snooping” (like [Scheffe])

One $\mu$ at a time

• One at a time confidence intervals:

$${\overline{y}}_i\pm t_{(1-\alpha /2,n-1}\sqrt{\frac{s_{ii}}{n}}$$

• Each of these intervals has a probability of $1-\alpha$ of covering the appropriate ${\mu }_i$

• But they ignore the covariance structure of the $p$ variables

• If we only care about $k$ simultaneous intervals, we can use “one at a time” method with the [Bonferroni] correction.

• This method gets more conservative as the number of intervals $k$ increases.

General Hypothesis Testing

One-sample Tests

$$H_0:C\mu =0$$

where

• $C$ is a $c\times p$ matrix of rank c where $c\le p$

We can test this hypothesis using the following statistic

$$F=\frac{n-c}{(n-1)c}{T}^2$$

where $T^2=n\left(C\overline{y}{)}^{\prime }\right({\mathrm{CSC}}^{\prime }{)}^{-1}\left(C\overline{y}\right)$

Example:

$$H_0:{\mu }_1={\mu }_2=...={\mu }_p$$

Equivalently,

$${\mu }_1-{\mu }_2=0⋮{\mu }_{p-1}-{\mu }_p=0$$

a total of $p-1$ tests. Hence, we have $C$ as the $p-1\times p$ matrix

$$C=\left(\begin{array}{ccccc}1& -1& 0& \dots & 0\\ 0& 1& -1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 1& -1\end{array}\right)$$

number of rows = $c=p-1$

Equivalently, we can also compare all of the other means to the first mean. Then, we test ${\mu }_1-{\mu }_2=0,{\mu }_1-{\mu }_3=0,...,{\mu }_1-{\mu }_p=0$, the $(p-1)\times p$ matrix $C$ is

$$C=\left(\begin{array}{ccccc}-1& 1& 0& \dots & 0\\ -1& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -1& 0& \dots & 0& 1\end{array}\right)$$

The value of $T^2$ is invariant to these equivalent choices of $C$

This is often used for repeated measures designs, where each subject receives each treatment once over successive periods of time (all treatments are administered to each unit).

Example:

Let $y_{ij}$ be the response from subject i at time j for $i=1,..,n,j=1,...,T$. In this case, $y_i=(y_{i1},...,y_{iT}{)}^{\prime },i=1,...,n$ are a random sample from $N_T(\mu ,\Sigma )$

Let $n=8$ subjects, $T=6$. We are interested in ${\mu }_1,..,{\mu }_6$

$$H_0:{\mu }_1={\mu }_2=...={\mu }_6$$

Equivalently,

$${\mu }_1-{\mu }_2=0{\mu }_2-{\mu }_3=0...{\mu }_5-{\mu }_6=0$$

We can test orthogonal polynomials for 4 equally spaced time points. To test for example the null hypothesis that quadratic and cubic effects are jointly equal to 0, we would define $C$

$$C=\left(\begin{array}{cccc}1& -1& -1& 1\\ -1& 3& -3& 1\end{array}\right)$$

Two-Sample Tests

Consider the analogous two sample multivariate tests.

Example: we have data on two independent random samples, one sample from each of two populations

$$y_{1i}\sim N_p\left({\mu }_1,\Sigma \right)y_{2j}\sim N_p\left({\mu }_2,\Sigma \right)$$

We assume

• normality

• equal variance-covariance matrices

• independent random samples

We can summarize our data using the sufficient statistics ${\overline{y}}_1,S_1,{\overline{y}}_2,S_2$ with respective sample sizes, $n_1,n_2$

Since we assume that ${\Sigma }_1={\Sigma }_2=\Sigma$, compute a pooled estimate of the variance-covariance matrix on $n_1+n_2-2$ df

$$S=\frac{(n_1-1)S_1+(n_2-1)S_2}{(n_1-1)+(n_2-1)}$$

$$H_0:{\mu }_1={\mu }_2{H}_a:{\mu }_1\ne {\mu }_2$$

At least one element of the mean vectors is different

We use

• ${\overline{y}}_1-{\overline{y}}_2$ to estimate ${\mu }_1-{\mu }_2$

• $S$ to estimate $\Sigma$

  Note: because we assume the two populations are independent, there is no covariance

  $cov({\overline{y}}_1-{\overline{y}}_2)=var\left({\overline{y}}_1\right)+var\left({\overline{y}}_2\right)=\frac{{\Sigma }_1}{n_1}+\frac{{\Sigma }_2}{n_2}=\Sigma (\frac{1}{n_1}+\frac{1}{n_2})$

Reject $H_0$ if

$$\begin{array}{cc}{T}^2& =({\overline{y}}_1-{\overline{y}}_2{)}^{\prime }\{S(\frac{1}{n_1}+\frac{1}{n_2}){\}}^{-1}({\overline{y}}_1-{\overline{y}}_2)\\ & =\frac{n_1{n}_2}{n_1+n_2}({\overline{y}}_1-{\overline{y}}_2{)}^{\prime }\{S{\}}^{-1}({\overline{y}}_1-{\overline{y}}_2)\\ & \ge \frac{(n_1+n_2-2)p}{n_1+n_2-p-1}{f}_{(1-\alpha ,n_1+n_2-p-1)}\end{array}$$

or equivalently, if

$$F=\frac{n_1+n_2-p-1}{(n_1+n_2-2)p}{T}^2\ge f_{(1-\alpha ,p,n_1+n_2-p-1)}$$

A $100(1-\alpha )\%$ confidence region for ${\mu }_1-{\mu }_2$ consists of all vector $\delta$ which satisfy

$$\frac{n_1{n}_2}{n_1+n_2}({\overline{y}}_1-{\overline{y}}_2-\delta {)}^{\prime }{S}^{-1}({\overline{y}}_1-{\overline{y}}_2-\delta )\le \frac{(n_1+n_2-2)p}{n_1+n_2-p-1}{f}_{(1-\alpha ,p,n_1+n_2-p-1)}$$

The simultaneous confidence intervals for all linear combinations of ${\mu }_1-{\mu }_2$ have the form

$$a^{\prime }({\overline{y}}_1-{\overline{y}}_2)\pm \sqrt{\frac{(n_1+n_2-2)p}{n_1+n_2-p-1}}{f}_{(1-\alpha ,p,n_1+n_2-p-1)}\times \sqrt{a^{\prime }\mathrm{Sa}(\frac{1}{n_1}+\frac{1}{n_2})}$$

Bonferroni intervals, for k combinations

$$({\overline{y}}_{1i}-{\overline{y}}_{2i})\pm t_{(1-\alpha /2k,n_1+n_2-2)}\sqrt{(\frac{1}{n_1}+\frac{1}{n_2})s_{ii}}$$

Model Assumptions

If model assumption are not met

• Unequal Covariance Matrices

  • If $n_1=n_2$ (large samples) there is little effect on the Type I error rate and power fo the two sample test

  • If $n_1>n_2$ and the eigenvalues of ${\Sigma }_1{\Sigma }_2^{-1}$ are less than 1, the Type I error level is inflated

  • If $n_1>n_2$ and some eigenvalues of ${\Sigma }_1{\Sigma }_2^{-1}$ are greater than 1, the Type I error rate is too small, leading to a reduction in power

• Sample Not Normal

  • Type I error level of the two sample $T^2$ test isn’t much affect by moderate departures from normality if the two populations being sampled have similar distributions

  • One sample $T^2$ test is much more sensitive to lack of normality, especially when the distribution is skewed.

  • Intuitively, you can think that in one sample your distribution will be sensitive, but the distribution of the difference between two similar distributions will not be as sensitive.

  • Solutions:

    • Transform to make the data more normal

    • Large large samples, use the ${\chi }^2$ (Wald) test, in which populations don’t need to be normal, or equal sample sizes, or equal variance-covariance matrices

      • $H_0:{\mu }_1-{\mu }_2=0$ use $({\overline{y}}_1-{\overline{y}}_2{)}^{\prime }(\frac{1}{n_1}{S}_1+\frac{1}{n_2}{S}_2{)}^{-1}({\overline{y}}_1-{\overline{y}}_2)\dot{\sim }{\chi }_{\left(p\right)}^2$

Equal Covariance Matrices Tests

With independent random samples from k populations of p-dimensional vectors. We compute the sample covariance matrix for each, $S_i$, where $i=1,...,k$

$$H_0:{\Sigma }_1={\Sigma }_2=\dots ={\Sigma }_k=\Sigma H_a:\text{at least 2 are different}$$

Assume $H_0$ is true, we would use a pooled estimate of the common covariance matrix, $\Sigma$

$$S=\frac{{\sum }_{i=1}^k(n_i-1)S_i}{{\sum }_{i=1}^k(n_i-1)}$$

with ${\sum }_{i=1}^k(n_i-1)$

Bartlett’s Test

(a modification of the likelihood ratio test). Define

$$N=\sum _{i=1}^k{n}_i$$

and (note: $\left|\right|$ are determinants here, not absolute value)

$$M=(N-k)\log \left|S\right|-\sum _{i=1}^k(n_i-1)\log \left|S_i\right|$$

$$C^{-1}=1-\frac{2p^2+3p-1}{6(p+1)(k-1)}\left\{\sum _{i=1}^k\right(\frac{1}{n_i-1})-\frac{1}{N-k}\}$$

• Reject $H_0$ when $M{C}^{-1}>{\chi }_{1-\alpha ,(k-1)p(p+1)/2}^2$

• If not all samples are from normal populations, $M{C}^{-1}$ has a distribution which is often shifted to the right of the nominal ${\chi }^2$ distribution, which means $H_0$ is often rejected even when it is true (the Type I error level is inflated). Hence, it is better to test individual normality first, or then multivariate normality before you do Bartlett’s test.

Two-Sample Repeated Measurements

• Define $y_{hi}=(y_{hi1},...,y_{hit}{)}^{\prime }$ to be the observations from the i-th subject in the h-th group for times 1 through T

• Assume that $y_{11},...,y_{1n_1}$ are iid $N_t({\mu }_1,\Sigma )$ and that $y_{21},...,y_{2n_2}$ are iid $N_t({\mu }_2,\Sigma )$

• $H_0:C({\mu }_1-{\mu }_2)=0_c$ where $C$ is a $c\times t$ matrix of rank $c$ where $c\le t$

• The test statistic has the form

$$T^2=\frac{n_1{n}_2}{n_1+n_2}({\overline{y}}_1-{\overline{y}}_2{)}^{\prime }{C}^{\prime }({\mathrm{CSC}}^{\prime }{)}^{-1}C({\overline{y}}_1-{\overline{y}}_2)$$

where $S$ is the pooled covariance estimate. Then,

$$F=\frac{n_1+n_2-c-1}{(n_1+n_2-2)c}{T}^2\sim f_{(c,n_1+n_2-c-1)}$$

when $H_0$ is true

If the null hypothesis $H_0:{\mu }_1={\mu }_2$ is rejected. A weaker hypothesis is that the profiles for the two groups are parallel.

$${\mu }_{11}-{\mu }_{21}={\mu }_{12}-{\mu }_{22}⋮{\mu }_{1t-1}-{\mu }_{2t-1}={\mu }_{1t}-{\mu }_{2t}$$

The null hypothesis matrix term is then

$H_0:C({\mu }_1-{\mu }_2)=0_c$ , where $c=t-1$ and

$$C={\left(\begin{array}{ccccc}1& -1& 0& \dots & 0\\ 0& 1& -1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & -1\end{array}\right)}_{(t-1)\times t}$$