One sample tests on mean vectors

Given data ${\bf X}_1,\ldots,{\bf X}_n$ iid $MVN(\boldsymbol\mu,\boldsymbol\Sigma)$ test

$$H_o: \boldsymbol\mu=\boldsymbol\mu_o$$

by computing

$$T^2 = n(\bar{{\bf X}} -\boldsymbol\mu_o)^T {{\bf S}}^{-1}(\bar{{\bf X}} -\boldsymbol\mu_o)$$

and getting $P$-values from $F$ distribution using theorem.

Example: no realistic ones. This hypothesis is not intrinsically useful. However: other tests can sometimes be reduced to it.

Example: Ten water samples split in half. One half of each to each of two labs. Measure biological oxygen demand (BOD) and suspended solids (SS). For sample $i$ let $X_{i1}$ be BOD for lab A, $X_{i2}$ be SS for lab A, $X_{i3}$ be BOD for lab B and $X_{i4}$ be SS for lab B. Question: are labs measuring the same thing? Is there bias in one or the other?

Notation ${\bf X}_i$ is vector of 4 measurements on sample $i$.

Data:

	Lab A		Lab B
Sample	BOD	SS	BOD	SS
1	6	27	25	15
2	6	23	28	13
3	18	64	36	22
4	8	44	35	29
5	11	30	15	31
6	34	75	44	64
7	28	26	42	30
8	71	124	54	64
9	43	54	34	56
10	33	30	29	20
11	20	14	39	21

Model: ${\bf X}_1,\ldots,{\bf X}_{11}$ are iid $MVN_4(\boldsymbol\mu,\boldsymbol\Sigma)$.

Multivariate problem because: not able to assume independence between any two measurements on same sample.

Potential sub-model: each measurement is

true mmnt + lab bias + mmnt error.

Model for measurement error vector ${\bf U}_i$ is multivariate normal mean 0 and diagonal covariance matrix $\boldsymbol\Sigma_{\bf U}$.

Lab bias is unknown vector $\boldsymbol\beta$.

True measurement should be same for both labs so has form

$$[Y_{i1},Y_{i2},Y_{i1},Y_{i2}]$$

where $Y_{i1},Y_{i2}$ are iid bivariate normal with unknown means $\theta_1,\theta_2$ and unknown $2 \times 2$ variance covariance $\boldsymbol\Sigma_{\bf Y}$.

This would give structured model

$${\bf X}_i = {\bf C}{\bf Y}+ \boldsymbol\beta + {\bf U}$$

where

$${\bf C}= \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \\ 1 & 0 \\ 0 & 1\end{array}\right]$$

This model has variance covariance matrix

$$\boldsymbol\Sigma_{\bf X}= {\bf C}\boldsymbol\Sigma_{\bf Y}{\bf C}^T + \boldsymbol\Sigma_{\bf U}$$

Notice that this matrix has only 7 parameters: four for the diagonal entries in $\boldsymbol\Sigma_{\bf U}$ and 3 for the entries in $\boldsymbol\Sigma_{\bf Y}$.

We skip this model and let $\boldsymbol\Sigma_{\bf X}$ be unrestricted.

Question of interest:

$$H_o: \mu_1=\mu_3$$    and $$\mu_2 = \mu_4$$

Reduction: partition ${\bf X}_i$ as

$$\left[\begin{array}{c} {\bf U}_i \\ {\bf V}_i \end{array}\right]$$

where ${\bf U}_i$ and ${\bf V}_i$ each have two components.

Define ${\bf W}_i={\bf U}_i$. Then our model makes ${\bf W}_i$ iid $MVN_2(\boldsymbol\mu_{\bf W}, \boldsymbol\Sigma_{\bf W})$. Our hypothesis is

$$H_o: \boldsymbol\mu_{\bf W}= 0$$

Carrying out our test in SPlus:

Working on CSS unix workstation:

Start SPlus then read in, print out data:

[61]ehlehl% mkdir .Data
[62]ehlehl% Splus
S-PLUS : Copyright (c) 1988, 1996 MathSoft, Inc.
S : Copyright AT&T.
Version 3.4 Release 1 for Sun SPARC, SunOS 5.3 : 1996 
Working data will be in .Data 
> # Read in and print out data
> eff <- read.table("effluent.dat",header=T)
> eff
   BODLabA SSLabA BODLabB SSLabB 
 1       6     27      25     15
 2       6     23      28     13
 3      18     64      36     22
 4       8     44      35     29
 5      11     30      15     31
 6      34     75      44     64
 7      28     26      42     30
 8      71    124      54     64
 9      43     54      34     56
10      33     30      29     20
11      20     14      39     21

Do some graphical preliminary analysis.

Look for non-normality, non-linearity, outliers.

Make plots on screen or saved in file.

> # Make pairwise scatterplots on screen using
> # motif graphics device and then in a postscript
> # file.
> motif()
> pairs(eff)
> postscript("pairs.ps",horizontal=F,
+   height=6,width=6)
> pairs(eff)
> dev.off()
Generated postscript file "pairs.ps".
 motif 
     2

Examine correlations

> cor(eff)
          BODLabA    SSLabA   BODLabB    SSLabB 
BODLabA 0.9999999 0.7807413 0.7228161 0.7886035
 SSLabA 0.7807413 1.0000000 0.6771183 0.7896656
BODLabB 0.7228161 0.6771183 1.0000001 0.6038079
 SSLabB 0.7886035 0.7896656 0.6038079 1.0000001

Notice high correlations.

Mostly caused by variation in true levels from sample to sample.

Get partial correlations.

Adjust for overall BOD and SS content of sample.

> aug <- cbind(eff,(eff[,1]+eff[,3])/2,
+                  (eff[,2]+eff[,4])/2)
> aug
   BODLabA SSLabA BODLabB SSLabB   X2   X3 
 1       6     27      25     15 15.5 21.0
 2       6     23      28     13 17.0 18.0
 3      18     64      36     22 27.0 43.0
 4       8     44      35     29 21.5 36.5
 5      11     30      15     31 13.0 30.5
 6      34     75      44     64 39.0 69.5
 7      28     26      42     30 35.0 28.0
 8      71    124      54     64 62.5 94.0
 9      43     54      34     56 38.5 55.0
10      33     30      29     20 31.0 25.0
11      20     14      39     21 29.5 17.5
> bigS <- var(aug)

Now compute partial correlations for first four variables given means of BOD and SS:

> S11 <- bigS[1:4,1:4]
> S12 <- bigS[1:4,5:6]
> S21 <- bigS[5:6,1:4]
> S22 <- bigS[5:6,5:6]
> S11dot2 <- S11 - S12 %*% solve(S22,S21)
> S11dot2
           BODLabA     SSLabA    BODLabB     SSLabB 
BODLabA  24.804665  -7.418491 -24.804665   7.418491
 SSLabA  -7.418491  59.142084   7.418491 -59.142084
BODLabB -24.804665   7.418491  24.804665  -7.418491
 SSLabB   7.418491 -59.142084  -7.418491  59.142084
> S11dot2SD <- diag(sqrt(diag(S11dot2)))
> S11dot2SD
         [,1]     [,2]     [,3]     [,4] 
[1,] 4.980428 0.000000 0.000000 0.000000
[2,] 0.000000 7.690389 0.000000 0.000000
[3,] 0.000000 0.000000 4.980428 0.000000
[4,] 0.000000 0.000000 0.000000 7.690389
> R11dot2 <- solve(S11dot2SD)%*% 
+           S11dot2%*%solve(S11dot2SD)
> R11dot2
          [,1]      [,2]      [,3]      [,4] 
[1,]  1.000000 -0.193687 -1.000000  0.193687
[2,] -0.193687  1.000000  0.193687 -1.000000
[3,] -1.000000  0.193687  1.000000 -0.193687
[4,]  0.193687 -1.000000 -0.193687  1.000000

Notice little residual correlation. Carry out Hotelling's $T^2$ test of $H_0: \boldsymbol\mu_{\bf W}=0$.

> w <- eff[,1:2]-eff[3:4]
> dimnames(w)<-list(NULL,c("BODdiff","SSdiff"))
> w
      BODdiff SSdiff 
 [1,]     -19     12
 [2,]     -22     10
  etc
 [8,]      17     60
  etc
> Sw <- var(w)
> cor(w)
          BODdiff    SSdiff 
BODdiff 1.0000001 0.3057682
 SSdiff 0.3057682 1.0000000
> mw <- apply(w,2,mean)
> mw
   BODdiff   SSdiff 
 -9.363636 13.27273
> Tsq <- 11*mw%*%solve(Sw,mw)
> Tsq
         [,1] 
[1,] 13.63931
> FfromTsq <- (11-2)*Tsq/(2*(11-1))
> FfromTsq
        [,1] 
[1,] 6.13769
> 1-pf(FfromTsq,2,9)
[1] 0.02082779

Conclusion: Pretty clear evidence of difference in mean level between labs. Which measurement causes the difference?

> TBOD <- sqrt(11)*mw[1]/sqrt(Sw[1,1])
> TBOD
   BODdiff 
 -2.200071
> 2*pt(TBOD,1)
   BODdiff 
 0.2715917
> 2*pt(TBOD,10)
    BODdiff 
 0.05243474
> TSS <- sqrt(11)*mw[2]/sqrt(Sw[2,2])
> TSS
  SSdiff 
 2.15153
> 2*pt(-TSS,10)
     SSdiff 
 0.05691733
> postscript("differences.ps",
+      horizontal=F,height=6,width=6)
> plot(w)
> abline(h=0)
> abline(v=0)
> dev.off()

Conclusion? Neither? Not a problem - summarizes evidence!

Problem: several tests at level 0.05 on same data. Simultaneous or Multiple comparisons.

In general can test hypothesis $H_o: {\bf C}\boldsymbol\mu=0$ by computing ${\bf Y}_i = {\bf C}{\bf X}_i$ and then testing $H_o: \boldsymbol\mu_{\bf Y}= 0$ using Hotelling's $T^2$.

Simultaneous confidence intervals

Confidence interval for ${\bf a}^T \boldsymbol\mu$:

$${\bf a}^T\bar{{\bf X}} \pm t_{n-1,\alpha/2}\sqrt{{\bf a}^T {\bf S}{\bf a}}$$

Based on $t$ distribution.

Give coverage intervals for 6 parameters of interest: 4 entries in $\boldsymbol\mu$ and $\boldsymbol\mu_1-\boldsymbol\mu_3$ and $\boldsymbol\mu_2-\boldsymbol\mu_4$

$\boldsymbol\mu_1$	$25.27 \pm 2.23\times 19.68/\sqrt{11}$
$\boldsymbol\mu_2$	$46.45 \pm 2.23\times 31.84/\sqrt{11}$
$\boldsymbol\mu_3$	$34.64 \pm 2.23\times 10.45/\sqrt{11}$
$\boldsymbol\mu_4$	$33.18 \pm 2.23\times 19.07/\sqrt{11}$
$\boldsymbol\mu_1-\boldsymbol\mu_3$	$-9.36 \pm 2.23\times 14.12/\sqrt{11}$
$\boldsymbol\mu_2-\boldsymbol\mu_4$	$13.27 \pm 2.23\times 20.46/\sqrt{11}$

Problem: each confidence interval has 5% error rate. Pick out last interval (on basis of looking most interesting) and ask about error rate?

Solution: adjust 2.23, $t$ multiplier to get

$$P($$all intervals cover truth$$) \ge 0.95$$

Rao or Scheffé type intervals

Based on inequality:

$$\left\vert{\bf a}^T {\bf b} \right\vert^2 \le {\bf a}^T {\bf M} {\bf a}{\bf b} ^T{\bf M}^{-1}{\bf b}$$

for any symmetric non-singular matrix ${\bf M}$.

Proof by Cauchy Schwarz: inner product of vectors ${\bf M}^{1/2} {\bf a}$ and ${\bf M}^{-1/2} {\bf B}$.

Put ${\bf b} = n^{-1/2}(\bar{{\bf X}} - \boldsymbol\mu)$ and ${\bf M}= {\bf S}$ to get

$$\left\vert n^{1/2}({\bf a}^T\bar{{\bf X}} - {\bf a}^T\boldsymbol\mu)\right\vert^2 \le {\bf a}^T{\bf S}{\bf a} T^2$$

This inequality is true for all $\bf a$. Thus the event that there is any $\bf a$ such that

$$\frac{({\bf a}^T\bar{{\bf X}} - {\bf a}^T\boldsymbol\mu)^2}{{\bf a}^T{\bf S}{\bf a}} > c$$

is a subset of the event

$$T^2 > c$$

Adjust $c$ to make the latter event have probability $\alpha$ by taking

$$c=\frac{p(n-1)}{n-p}F_{p,n-p,\alpha}.$$

Then the probability that every one of the uncountably many confidence intervals

$${\bf a}^T \bar{{\bf X}} \pm \sqrt{c} \sqrt{{\bf a}^T{\bf S}{\bf a}}$$

covers the corresponding true parameter value is at least $1-\alpha$.

In fact the probability of this happening is exactly equal to $1-\alpha$ because for each data set the supremum of

$$\frac{({\bf a}^T\bar{{\bf X}} - {\bf a}^T\boldsymbol\mu)^2}{{\bf a}^T{\bf S}{\bf a} }$$

over all $\bf a$ is $T^2$.

Our case

$$\sqrt{c} = \frac{4(10)}{7} F_{4,7,0.05} = 4.85$$

Coverage probability of single interval using $\sqrt{c}4.85$? From $t$ distribution:

$$99.93\%$$

Probability all 6 intervals would cover using $\sqrt{c}4.85$?

Use Bonferroni inequality:

$$P(\cup A_i) \le \sum P(A_i)$$

Simultaneous coverage probability of 6 intervals using $\sqrt{c}4.85$

$$\ge 1-6*(1-0.9993) = 99.59\%$$

Usually just use

$$\sqrt{c} =t_{n-1,\alpha/12} =3.28$$

General Bonferroni strategy. If we want intervals for $\theta_1, \ldots,\theta_k$ get interval for $\theta_i$ at level $1-\alpha/k$. Simultaneous coverage probability is at least $1-\alpha$. Notice that Bonferroni narrower in our example unless $0.0007k=0.5$ giving $k> 71$.

Motivations for $T^2$:

1: Hypothesis $H_o: \boldsymbol\mu=0$ is true iff all hypotheses $H_{o{\bf a}}: {\bf a}^T\boldsymbol\mu=0$ are true. Natural test for $H_{o{\bf a}}$ rejects if

$$t({\bf a}) = \frac{n^{1/2}{\bf a}^T(\bar{{\bf X}}-\boldsymbol\mu)}{ {\bf a}^T {\bf S}{\bf a}}$$

large. So take largest test statistic.

Fact:

$$\sup_{\bf a} t^2({\bf a}) = T^2$$

Proof: like calculation of maximal correlation.

2: likelihood ratio method.

Compute

$$\ell(\hat{\boldsymbol\mu},\hat{\boldsymbol\Sigma}) - \ell( \hat{\boldsymbol\mu}_o,\hat{\boldsymbol\Sigma}_o)$$

where the subscript $o$ indicates estimation assuming $H_o$.

In our case to test $H_o: \boldsymbol\mu=0$ find

$$\hat{\boldsymbol\mu}_o = 0 \quad \hat{\boldsymbol\Sigma}_o = \sum {\bf X}_i{\bf X}_i^T/n$$

and

$$\ell(\hat{\boldsymbol\mu},\hat{\boldsymbol\Sigma}) - \ell( \hat{\... ...o) = n \log\{\det(\hat{\boldsymbol\Sigma})/\det(\hat{\boldsymbol\Sigma}_0)\}/2$$

Now write

$$\sum {\bf X}_i{\bf X}_i^T = \sum({\bf X}_i-\bar{{\bf X}})({\bf X}_i-\bar{{\bf X}})^T + n\bar{{\bf X}}\bar{{\bf X}}^T$$

Use formula:

$$\det({\bf A}+{\bf v}{\bf v}^T) = \det({\bf A})(1+{\bf v}^T{\bf A}^{-1}{\bf v})$$

to get

$$\det(n\hat{\boldsymbol\Sigma}_0) =\det(n\hat{\boldsymbol\Sigma})( 1+n\bar{{\bf X}}^T\left(n\hat{\boldsymbol\Sigma}\right)^{-1}\bar{{\bf X}})$$

so that the ratio of determinants is a monotone increasing function of $T^2$.

Again conclude: likelihood ratio test rejects for $T^2>c$ where $c$ chosen to make level $\alpha$.

3: compare estimates of $\boldsymbol\Sigma$.

In univariate regression $F$ tests to compare a restricted model with a full model have form

$$\frac{ESS_\text{Restricted}-ESS_\text{Full}}{ESS_\text{Full}}\frac{df_\text{Error}}{df_\text{difference}}$$

This is a monotone function of

$$\frac{\hat\sigma^2_\text{Restricted}}{\hat\sigma^2_\text{Full}}$$

where $ESS$ denotes an error sum of squares and $\hat\sigma^2$ and estimate of the residual variance - $ESS/df$.

Here: substitute matrices.

Analogue of ESS for full model:

$${\bf E}$$

Analogue of ESS for reduced model:

$${\bf E}+{\bf H}$$

(This defined ${\bf H}$ to be the change in the Sum of Squares and Cross Products matrix.)

In 1 sample example:

$${\bf E}= \sum ({\bf X}_i-\bar{{\bf X}})({\bf X}_i-\bar{{\bf X}})^T$$

and

$${\bf E}+{\bf H}= \sum {\bf X}_i{\bf X}_i^T$$

Test of $\boldsymbol\mu=0$ based on comparing

$${\bf H}=\sum {\bf X}_i{\bf X}_i^T - \sum ({\bf X}_i-\bar{{\bf X}})({\bf X}_i-\bar{{\bf X}})^T = n\bar{{\bf X}}\bar{{\bf X}}^T$$

to

$${\bf E}= \sum ({\bf X}_i-\bar{{\bf X}})({\bf X}_i-\bar{{\bf X}})^T= (n-1){\bf S}$$

To make comparison. If null true

$${\rm E}(n\bar{{\bf X}}\bar{{\bf X}}^T) = \boldsymbol\Sigma$$

and

$${\rm E}({\bf S}) = \boldsymbol\Sigma$$

so try tests based on closeness of

$${\bf S}^{-1}n\bar{{\bf X}}\bar{{\bf X}}^T$$

to identity matrix.

Measures of size based on eigenvalues of

$${\bf S}^{-1}n\bar{{\bf X}}\bar{{\bf X}}^T$$

which are same as eigenvalues of

$${\bf S}^{-1/2}n\bar{{\bf X}}\bar{{\bf X}}^T{\bf S}^{-1/2}$$

Suggested size measures for ${\bf A}-{\bf I}$:

• ${\rm trace}({\bf A}-{\bf I})$ (= sum of eigenvalues).

• $\det({\bf A}-{\bf I})$ (= product of eigenvalues).

• maximum eigenvalue of ${\bf A}-{\bf I}$.

For our matrix: eigenvalues all 0 except for one. (So really-matrix not close to ${\bf I}$.)

Largest eigenvalue is

$$T^2 = n\bar{{\bf X}}{\bf S}^{-1}\bar{{\bf X}}$$

But: see two sample problem for precise tests based on suggestions.

Two sample problem

Given data ${\bf X}_{ij}; j=1,\ldots,n_i;i=1,2$. Model ${\bf X}_{ij} \sim MVN_p(\boldsymbol\mu_i,\boldsymbol\Sigma_i)$, independent.

Test $H_o: \boldsymbol\mu_1 = \boldsymbol\mu_2$.

Case 1: for motivation only. $\boldsymbol\Sigma_i$ known $i=1,2$.

Natural test statistic: based on

$${\bf D}= \bar{{\bf X}}_1 -\bar{{\bf X}}_2$$

which has $MVN(\boldsymbol\mu_{\bf D},\boldsymbol\Sigma_{\bf D})$ where

$$\boldsymbol\mu_{\bf D}= \boldsymbol\mu_1-\boldsymbol\mu_2$$

and

$$\boldsymbol\Sigma_{\bf D}= n_1^{-1} \boldsymbol\Sigma_1+ n_2^{-1}\boldsymbol\Sigma_2$$

So

$${\bf D}^T \left({\boldsymbol\Sigma}_1/n_1+{\boldsymbol\Sigma}_2/n_2\right)^{-1} {\bf D}$$

has a $\chi_p^2$ distribution if null true.

If $\boldsymbol\Sigma_i$ not known must estimate. No universally agreed best procedure (even for $p=1$ -- called Behrens-Fisher problem).

Usually: assume $\boldsymbol\Sigma_1 = \boldsymbol\Sigma_2$.

If so: MLE of $\boldsymbol\mu_i$ is $\bar{{\bf X}}_i$ and of $\boldsymbol\Sigma$ is

$$\frac{\sum_{ij} ({\bf X}_{ij}-\bar{{\bf X}}_i)({\bf X}_{ij}-\bar{{\bf X}}_i)^T}{n_1+n_2}$$

Usual estimate of $\boldsymbol\Sigma$ is

$${\bf S}_$$Pooled$$= \frac{\sum_{ij} ({\bf X}_{ij}-\bar{{\bf X}}_i)({\bf X}_{ij}-\bar{{\bf X}}_i)^T}{n_1+n_2-2}$$

which is unbiased.

Possible test developments:

1) By analogy with 1 sample:

$$T^2 = \frac{n_1n_2}{n_1+n_2} {\bf D}^T {\bf S}^{-1}_$$Pooled$${\bf D}$$

which has the distribution

$$\frac{n_1+n_2-1-p}{p(n_1+n_2-2)} T^2 \sim F_{p,n_1+n_2-1-p}$$

2) Union-intersection: test of ${\bf a}^t({\boldsymbol\mu}_1-{\boldsymbol\mu}_2) = 0$ based on

$$t_a = \sqrt{\frac{n_1n_2}{n_1+n_2}}\frac{{\bf a}^T {\bf D}}{\sqrt{{\bf a}^T {\bf S}{\bf a}}}$$

Maximize $t^2$ over all $\bf a$.

Get

$$T^2 = \frac{n_1n_2}{n_1+n_2} {\bf D}^T {\bf S}^{-1}{\bf D}$$

3) Likelihood ratio: the MLE of ${\boldsymbol\Sigma}$ for the unrestricted model is

$$\frac{n_1+n_2-2}{n_1+n_2} {\bf S}$$

Under the hypothesis ${\boldsymbol\mu}_1={\boldsymbol\mu}_2$ the mle of ${\boldsymbol\Sigma}$ is

$$\frac{\sum_{ij} ({\bf X}_{ij} -\bar{\bar{{\bf X}}})({\bf X}_{ij} -\bar{\bar{{\bf X}}})^T}{n_1+n_2}$$

where

$$\bar{\bar{{\bf X}}} = \frac{ n_1\bar{{\bf X}}_1 + n_2\bar{{\bf X}}_2}{n_1+n_2}$$

This simplifies to

$$\frac{{\bf E}+{\bf H}}{n_1+n_2}$$

The log-likelihood ratio is a multiple of

$$\log\det\hat{\boldsymbol\Sigma}_$$Full$$- \log\det\hat{\boldsymbol\Sigma}_$$Restricted$$$$

which is a function of

$$\log\left\{\det{\bf E}/\det({\bf E}+{\bf H})\right\}$$

or equivalently a function of Wilk's $\Lambda$:

$$\Lambda = \frac{\det{\bf E}}{\det({\bf E}+{\bf H})} = \frac{1}{\det({\bf H}{\bf E}^{-1}+{\bf I})}$$

Compute det: multiply together eigenvalues.

If $\lambda_i$ are the eigenvalues of ${\bf H}{\bf E}^{-1}$ then

$$\Lambda = \frac{1}{\prod(1+\lambda_i)}$$

Two sample analysis in SAS on css network

• Type sas to start system.

• Several windows open. Go to Program Editor.

• Under file menu open file with sas code. Contents of sas2sample.sas

    data long;
    infile 'tab57sh';
    input group a b c;
run;
proc print;
run;
proc glm;
    class group;
    model a b c = group;
    manova h=group / printh printe;
run;

Notes:

1) First 4 lines form DATA step:

a) creates data set named long by reading in 4 columns of data from file named tab57sh stored in same directory as I was in when I typed sas.

b) Calls variables group (=1 or 2 as label for the two groups) and a, b, c which are names for the 3 test scores for each subject.

2) Next two lines: print out data: result is (slightly edited)

    Obs    group     a     b     c
      1      1      19    20    18
      2      1      20    21    19
      3      1      19    22    22
     etc till
     11      2      15    17    15
     12      2      13    14    14
     13      2      14    16    13

3) Then use proc glm to do analysis:

a) class group declares that the variable group defines levels of a categorical variable.

b) model statement says to regress the variables a, b, c on variable group.

c) manova statement says to do both 3 univariate regressions and a mulivariate regression and to print out the ${\bf H}$ and ${\bf E}$ matrices where ${\bf H}$ is the matrix corresponding to the presence of the factor group in the model.

Output of MANOVA: First univariate results

    The GLM Procedure
 Class Level Information
        Class         Levels    Values
        group              2    1 2
 Number of observations    13
Dependent Variable: a
                Sum of
Source    DF   Squares  Mean Square F Value Pr > F
Model      1  54.276923  54.276923   19.38 0.0011
Error     11  30.800000   2.800000
Corrd Tot 12  85.076923
R-Square   Coeff Var      Root MSE        a Mean
0.637975    10.21275      1.673320      16.38462
Source    DF  Type ISS  Mean Square F Value Pr > F
group      1  54.276923  54.276923   19.38 0.0011
Source    DF  TypeIIISS Mean Square F Value Pr > F
group      1  54.276923  54.276923   19.38 0.0011

Dependent Variable: b
                Sum of
Source    DF   Squares  Mean Square F Value Pr > F
Model      1  70.892308  70.892308  34.20 0.0001
Error     11  22.800000   2.072727
Corrd Tot 12  93.692308

Dependent Variable: c
                Sum of
Source    DF   Squares Mean Square F Value Pr > F
Model      1   94.77692  94.77692   39.64 <.0001
Error     11   26.30000   2.39090
Corrd Tot 12  121.07692

The matrices ${\bf E}$ and ${\bf H}$.

  E = Error SSCP Matrix
   a     b     c
a 30.8  12.2  10.2
b 12.2  22.8   3.8
c 10.2   3.8  26.3
Partial Correlation Coefficients from 
    the Error SSCP Matrix / Prob > |r|
DF = 11     a              b              c
a    1.000000       0.460381       0.358383
                      0.1320         0.2527
b    0.460381       1.000000       0.155181
       0.1320                        0.6301
c    0.358383       0.155181       1.000000
       0.2527         0.6301
H = Type III SSCP Matrix for group
     a                 b                 c
a   54.276923077   62.030769231   71.723076923
b   62.030769231   70.892307692   81.969230769
c   71.723076923   81.969230769   94.776923077

The eigenvalues of ${\bf E}^{-1}{\bf H}$.

Characteristic Roots and Vectors of: E Inverse * H
H = Type III SSCP Matrix for group
E = Error SSCP Matrix
Characteristic       Characteristic Vector V'EV=1
Root     Percent        a           b           c
5.816159 100.00  0.00403434  0.12874606 0.13332232
0.000000   0.00 -0.09464169 -0.10311602 0.16080216
0.000000   0.00 -0.19278508  0.16868694 0.00000000
MANOVA Test Criteria and Exact F Statistics 
 for the Hypothesis of No Overall group Effect
H = Type III SSCP Matrix for group
E = Error SSCP Matrix
S=1    M=0.5    N=3.5
Statistic            Value   F   NumDF DenDF Pr > F
Wilks' Lambda       0.1467 17.45    3     9  0.0004
Pillai's Trace      0.8533 17.45    3     9  0.0004
Hotelling-Lawley Tr 5.8162 17.45    3     9  0.0004
Roy's Greatest Root 5.8162 17.45    3     9  0.0004

Things to notice:

1. The conclusion is clear. The mean vectors for the two groups are not the same.

2. The four statistics have the following definitions in terms of eigenvalues of ${\bf E}^{-1}{\bf H}$:

   Wilk's Lambda:

   $$\frac{1}{\prod(1+\lambda_i)}=\frac{1}{6.816}$$
   Pillai's trace:
   $${\rm trace}({\bf H}({\bf H}+{\bf E})^{-1}) = \sum\frac{\lambda_i}{1+\lambda_i} = \frac{5.816}{6.816}$$
   Hotelling-Lawley trace:
   $${\rm trace}({\bf H}{\bf E}^{-1}) = \sum\lambda_i = 5.816$$
   Roy's greatest Root:
   $$\max\{\lambda_i\} = 5.816$$

1 way layout

Also called $m$ sample problem.

Data ${\bf X}_{ij},j=1,\ldots,n_i;i=1,\ldots, m$.

Model ${\bf X}_{ij}$ independent $MVN_p({\boldsymbol\mu}_i,{\boldsymbol\Sigma})$.

First problem of interest: test

$$H_o: {\boldsymbol\mu}_1= \cdots = {\boldsymbol\mu}_m$$

Based on ${\bf E}$ and ${\bf H}$. MLE of ${\boldsymbol\mu}_i$ is $\bar{{\bf X}}_i$.

$${\bf E}= \sum_{ij} ({\bf X}_{ij} - \bar{{\bf X}}_i)({\bf X}_{ij} - \bar{{\bf X}}_i)^T$$

Under $H_o$ MLE of ${\boldsymbol\mu}$, the common value of the ${\boldsymbol\mu}_i$ is

$$\bar{\bar{{\bf X}}} = \frac{\sum_{ij} {\bf X}_{ij}}{\sum_i n_i}$$

So

$${\bf E}+{\bf H}= \sum_{ij} ({\bf X}_{ij} - \bar{\bar{{\bf X}}})({\bf X}_{ij} - \bar{\bar{{\bf X}}})^T$$

This makes

$${\bf H}= \sum_{ij}(\bar{{\bf X}}_i - \bar{\bar{{\bf X}}})(\bar{{\bf X}}_i - \bar{\bar{{\bf X}}})^T$$

Notice can do sum over $j$ to get factor of $n_i$:

$${\bf H}= \sum_i n_i(\bar{{\bf X}}_i - \bar{\bar{{\bf X}}})(\bar{{\bf X}}_i - \bar{\bar{{\bf X}}})^T$$

Note: rank of ${\bf H}$ is minimum of $p$ and $m-1$. The data

1 19 20 18
1 20 21 19
1 19 22 22
1 18 19 21
1 16 18 20
1 17 22 19
1 20 19 20
1 15 19 19
2 12 14 12
2 15 15 17
2 15 17 15
2 13 14 14
2 14 16 13
3 15 14 17
3 13 14 15
3 12 15 15
3 12 13 13
4  8  9 10
4 10 10 12
4 11 10 10
4 11  7 12

Code

    data three;
    infile 'tab57for3sams';
    input group a b c;
run;
proc print;
run;
proc glm;
    class group;
    model a b c = group;
    manova h=group / printh printe;
run;
    data four;
    infile 'table5.7';
    input group a b c;
run;
proc print;
run;
proc glm;
    class group;
    model a b c = group;
    manova h=group / printh printe;
run;

Pieces of output: first set of code does first 3 groups.

So: ${\bf H}$ has rank 2.

Characteristic Roots & Vectors of: E Inverse * H
 Characteristic     Characteristic Vector V'EV=1
      Root Percent     a        b        c
6.90568180  96.94   0.01115  0.14375 0.08795
0.21795125   3.06  -0.07763 -0.09587 0.16926
0.00000000   0.00  -0.18231  0.13542 0.02083
S=2    M=0    N=5
Statistic         Value   F   NumDF  Den DF Pr > F
Wilks'           0.1039  8.41   6      24   <.0001
Pillai's         1.0525  4.81   6      26   0.0020
Hotelling-Lawley 7.1236 13.79   6   14.353  <.0001
Roy's            6.9057 29.92   3      13   <.0001
NOTE: F Statistic for Roy's is an upper bound.
NOTE: F Statistic for Wilks'is exact.

Notice two eigenvalues not 0. Note that exact distribution for Wilk's Lambda is available. Now 4 groups

      Root Percent      a       b       c
15.3752900  98.30   0.01128  0.13817 0.08126
0.2307260    1.48  -0.04456 -0.09323 0.15451
0.0356937    0.23  -0.17289  0.09020 0.04777
S=3    M=-0.5    N=6.5
Statistic          Value    F    NumDF Den DF Pr > F
Wilks'         0.04790913 10.12    9   36.657 <.0001
Pillai's       1.16086747  3.58    9       51 0.0016
Hot'ng-Lawley 15.64170973 25.02    9   20.608 <.0001
Roy's         15.37528995 87.13    3       17 <.0001
NOTE: F Statistic for Roy's is an upper bound.

Other Hypotheses

How do mean vectors differ? One possibility:

$${\boldsymbol\mu}_{ik} - {\boldsymbol\mu}_{jk} = c_{i} -c_{j}$$

for constants $c_{i}$ and $c_j$ which do not depend on $k$. This is an additive model for the means.

Test ?

Define $\alpha = \sum_{ij} {\boldsymbol\mu}_{ij} / (pk)$ Then put

$$\beta_i = \sum_j{\boldsymbol\mu}_{ij} /p -\alpha,$$

$$\gamma_j = \sum_i{\boldsymbol\mu}_{ij}/k -\alpha$$

$$\tau_{ij} = {\boldsymbol\mu}_{ij} -\beta_i-\gamma_j - \alpha$$

If the $\tau_{ij}$ are all 0 then

$${\boldsymbol\mu}_{ik} - {\boldsymbol\mu}_{jk} = \beta_i-\beta_j$$

so we test the hypothesis that all $\tau_{ij}$ are 0.

Univariate Two Way Anova

Data $Y_{ijk}$

$k=1,\ldots,n_{ij};j=1,\ldots,p;i=1,\ldots,m$.

Model: independent, $Y_{ijk}\sim N(\mu_{ij},\sigma^2)$.

Note: this is the fixed effects model.

Usual approach: define grand mean, main effects, interactions:

$$\mu = \sum_{ijk} \mu_{ij} /\sum_{ij}n_{ij}$$

$$\alpha_i = \sum_{jk}\mu_{ij} /\sum_j n_{ij} - \mu$$

$$\beta_j = \sum_{ik}\mu_{ij} /\sum_i n_{ij} - \mu$$

$$\gamma_{ij} = \mu_{ij} -(\mu+\alpha_i +\beta_j)$$

Test additive effects: $\gamma_{ij} = 0$ for all $i,j$.

Usual test based on ANOVA:

Stack observations $Y_{ijk}$ into vector ${\bf Y}$, say.

Estimate $\mu$, $\alpha_i$, etc by least squares.

Form vectors with entries $\hat\mu$, $\hat\alpha_i$ etc.

Write

$${\bf Y} = \hat{\boldsymbol\mu}+ \hat{\boldsymbol\alpha} + \hat{\boldsymbol\beta} + \hat{\boldsymbol\gamma} + \hat{\boldsymbol\epsilon}$$

This defines the vector of fitted residuals $\hat{\boldsymbol\epsilon}$.

Fact: all vectors on RHS are independent and orthogonal. So:

$$\vert\vert{\bf Y}\vert\vert^2 = \vert\vert \hat{\boldsymbol\mu}\v... ...oldsymbol\gamma}\vert\vert^2 + \vert\vert\hat{\boldsymbol\epsilon}\vert\vert^2$$

This is the ANOVA table. Usually we defined the corrected total sum of squares to be

$$\vert\vert{\bf Y}\vert\vert^2 - \vert\vert \hat{\boldsymbol\mu}\vert\vert^2$$

Our problem is like this one BUT the errors are not modeled as independent.

In the analogy:

$i$ labels group.

$j$ labels the columns: ie $j$ is a, b, c.

$k$ runs from 1 to $n_{ij}=n_i$.

But

$${\rm Cov}(Y_{ijk},Y_{i^\prime j^\prime k^\prime}) = \begin{... ... i = i^\prime, k=k^\prime \\ 0 & \text{otherwise} \end{cases}$$

Now do analysis in SAS.

Tell SAS that the variables A, B and C are repeated measurements of the same quantity.

proc glm;
    class group;
    model a b c = group;
    repeated scale;
run;

The results are as follows:

General Linear Models Procedure
Repeated Measures Analysis of Variance
Repeated Measures Level Information
Dependent Variable  A   B    C
Level of SCALE      1   2    3

Manova Test Criteria and Exact F 
   Statistics for the Hypothesis of no
   SCALE Effect
H = Type III SS&CP Matrix for SCALE   
               E = Error SS&CP Matrix
S=1    M=0    N=7
Statistic
                  Value     F   NumDF DenDF Pr > F
Wilks' Lambda    0.56373  6.1912   2    16  0.0102
Pillai's Trace   0.43627  6.1912   2    16  0.0102
Hotelling-Lawley 0.77390  6.1912   2    16  0.0102
Roy's            0.77390  6.1912   2    16  0.0102

Note: should look at interactions first.

Manova Test Criteria and F Approximations 
 for the Hypothesis of no SCALE*GROUP Effect
S=2    M=0    N=7
Statistic          Value   F   NumDF DenDF Pr > F
Wilks' Lambda    0.56333 1.7725  6    32  0.1364
Pillai's Trace   0.48726 1.8253  6    34  0.1234
Hotelling-Lawley 0.68534 1.7134  6    30  0.1522
Roy's            0.50885 2.8835  3    17  0.0662
NOTE: F Statistic for Roy's Greatest 
        Root is an upper bound.
NOTE: F Statistic for Wilks' Lambda is exact.

Only weak evidence of interaction. Repeated statement: univariate anova. Results:

Repeated Measures Analysis of Variance
Tests of Hypotheses for Between Subjects Effects
Source DF Type III SS Mean Square  F    Pr > F
GROUP   3  743.900000  247.966667 70.93 0.0001
Error  17   59.433333    3.496078
Repeated Measures Analysis of Variance
Univariate Tests of Hypotheses for 
          Within Subject Effects
Source: SCALE                    Adj Pr > F
DF TypeIIISS MS      F   Pr > F G - G  H - F
 2  16.624   8.312  5.39 0.0093 0.0101 0.0093
Source: SCALE*GROUP
DF TypeIII   MS    F   Pr > F G - G  H - F
6   18.9619 3.160 2.05 0.0860 0.0889 0.0860
Source: Error(SCALE)
DF TypeIII SS     Mean Square
34  52.4667     1.54313725
Greenhouse-Geisser Epsilon = 0.9664
       Huynh-Feldt Epsilon = 1.2806

Greenhouse-Geisser, Huynh-Feldt test to see if ${\boldsymbol\Sigma}$ has certain structure.

Return to 2 way anova model. Express as:

$$Y_{ijk} = \mu+\alpha_i+\beta_j+ \gamma_{ij}+\epsilon_{ijk}$$

For fixed effects model is $\epsilon_{ijk}$ iid $N(0,\sigma^2)$.

For MANOVA model vector of $\epsilon_{ijk}$ is MVN but with covariance as for $Y$.

Intermediate model. Put in subject effect.

Assume

$$\epsilon_{ijk} = \delta_{ik} + u_{ijk}$$

where $u_{ijk}$ iid $N(0,\sigma^2)$ and $\delta_{ik}$ are iid $N(0,\tau^2)$. Then

$$\begin{multline*} {\rm Cov}(Y_{ijk},Y_{i^\prime j^\prime k^\prime}) = \\ \begin... ..., k=k^\prime, j\neq j^\prime \\ 0 & \text{otherwise} \end{cases}\end{multline*}$$

This model is usually not fitted by maximum likelihood but by analyzing the behaviour of the ANOVA table under this model.

Essentially model says

$${\boldsymbol\Sigma}= \tau^2 {\bf 1}{\bf 1}^T + \sigma^2 {\bf I}$$

GG, HF test for slightly more general pattern for ${\boldsymbol\Sigma}$.

Do univariate anova: The data reordered:

1 1 1 19
1 1 2 20
1 1 3 18
2 1 1 20
2 1 2 21
2 1 3 19
 et cetera
2 4 2 10
2 4 3 12
3 4 1 11
3 4 2 10
3 4 3 10
4 4 1 11
4 4 2 7
4 4 3 12

The four columns are now labels for subject number, group, scale (a, b or c) and the response. The sas commands:

    data long;
    infile 'table5.7uni';
    input subject group scale score;
run;
proc print;
run;
proc glm;
    class group;
    class scale;
    class subject;
    model score =group subject(group) 
                 scale group*scale;
    random subject(group) ;
run;

Some of the output:

Dependent Variable: SCORE
               Sum of    Mean
Source   DF    Squares  Square F    Pr > F
Model    28   843.5333 30.126 19.52 0.0001
Error    34    52.4667  1.543
Total  62   896.0000
Root MSE      SCORE Mean
1.242231        15.33333


Source         DF  TypeISS    MS      F     Pr > F

GROUP           3 743.9000 247.9667 160.69  0.0001
SUBJECT(GROUP) 17  59.4333   3.4961   2.27  0.0208
SCALE           2  21.2381  10.6190   6.88  0.0031
GROUP*SCALE     6  18.9620   3.1603   2.05  0.0860

Source         DF TypeIIISS     MS    F     Pr > F
GROUP           3 743.9000  247.9667 160.69 0.0001
SUBJECT(GROUP) 17  59.4333    3.4961   2.27 0.0208
SCALE           2  16.6242    8.3121   5.39 0.0093
GROUP*SCALE     6  18.9619    3.1603   2.05 0.0860

Source         Type III Expected Mean Square

GROUP          Var(Error) + 3 Var(SUBJECT(GROUP)) 
                    + Q(GROUP,GROUP*SCALE)
SUBJECT(GROUP) Var(Error) + 3 Var(SUBJECT(GROUP))
SCALE          Var(Error) + Q(SCALE,GROUP*SCALE)
GROUP*SCALE    Var(Error) + Q(GROUP*SCALE)

Type I Sums of Squares:

• Sequential sums of squares.

• Each line is a sum of squares comparing the model with effects listed above to one with one extra effect.

• Depend on order terms listed in model.

Type III Sums of Squares:

• Roughly: each line compares model with all other effects in model.

• In unbalanced designs be careful about the differences between Types II, III and IV.

Notice hypothesis of no group by scale interaction is acceptable.

Under the assumption of no such group by scale interaction the hypothesis of no group effect is tested by dividing group ms by subject(group) ms.

Value is 70.9 on 3,17 degrees of freedom.

This is NOT the F value in the table above since the table above is for FIXED effects.

Notice that the sums of squares in this table match those produced in the repeated measures ANOVA. This is not accidental.

Two Way MANOVA

Data ${\bf Y}_{ijk}$

$k=1,\ldots,n_{ij};j=1,\ldots,p;i=1,\ldots,m$.

Model: independent, ${\bf Y}_{ijk}\sim MVN_p({\boldsymbol\mu}_{ij},{\boldsymbol\Sigma}^2)$.

Note: fixed effects model.

Usual approach: define grand mean, main effects, interactions:

$$\bar{\bar{{\boldsymbol\mu}}} = \sum_{ijk} {\boldsymbol\mu}_{ij} /\sum_{ij}n_{ij}$$

$$\boldsymbol\alpha_i = \sum_{jk}{\boldsymbol\mu}_{ij} /\sum_j n_{ij} - \bar{\bar{{\boldsymbol\mu}}}$$

$$\boldsymbol\beta_i = \sum_{ik}{\boldsymbol\mu}_{ij} /\sum_i n_{ij} - \bar{\bar{{\boldsymbol\mu}}}$$

$$\boldsymbol\gamma_{ij} = {\boldsymbol\mu}_{ij} - (\bar{\bar{{\boldsymbol\mu}}}+ \boldsymbol\alpha_i +\boldsymbol\beta_j)$$

Test additive effects: $\gamma_{ij} = 0$ for all $i,j$.

Exactly parallel to univariate 2 way ANOVA:

Estimate $\mu$, $\alpha_i$, etc by least squares.

Formulas exactly like univariate.

Form vectors with entries $\hat\mu$, $\hat\alpha_i$ etc.

Write

$${\bf Y} = \hat{\boldsymbol\mu}+ \hat{\boldsymbol\alpha} + \hat{\boldsymbol\beta} + \hat{\boldsymbol\gamma} + \hat{\boldsymbol\epsilon}$$

but now in matrix form; one column for each response variable, one row for each case.

This defines the matrix of fitted residuals $\hat{\boldsymbol\epsilon}$.

Fact:

$${\bf Y}^T{\bf Y}= \hat{\boldsymbol\mu}^T \hat{\boldsymbol\mu} + ... ...\hat{\boldsymbol\gamma} + \hat{\boldsymbol\epsilon}^T\hat{\boldsymbol\epsilon}$$

This is like the ANOVA table. Read off formulas for various ${\bf H}$ matrices.

Also

$${\bf E}= \hat{\boldsymbol\epsilon}^T\hat{\boldsymbol\epsilon}$$

The SAS commands for a two way analysis of variance with 3 response variables.

data mas;
    infile 'mas';
    input row column y1 y2 y3;
proc print;
proc glm;
     class row column;
     model y1-y3 = row | column;
           manova h=_all_ / printh printe;
run;

Features of the code:

• The model statement requests that effects for row, column and the interaction of row and column be included in the model.

• h=_all_ requests that all the effects in the model be tested for.

The data:

OBS   ROW    COLUMN     Y1      Y2      Y3
 1     1        1      18.2    16.5    0.2
 2     1        1      18.7    19.5    0.3
 3     1        1      19.5    19.8    0.2
 4     1        2      19.2    19.5    0.2
 5     1        2      18.4    19.8    0.2
 6     1        2      20.7    19.4    0.2
 7     2        1      21.3    23.3    0.3
 8     2        1      19.6    22.3    0.5
 9     2        1      20.2    19.0    0.4
10     2        2      18.9    22.0    0.3
11     2        2      20.7    21.1    0.2
12     2        2      21.6    20.3    0.2
13     3        1      20.7    16.7    0.3
14     3        1      21.0    19.3    0.4
15     3        1      17.2    15.9    0.3
16     3        2      20.2    19.0    0.2
17     3        2      18.4    17.9    0.3
18     3        2      20.9    19.9    0.2

Here is the matrix ${\bf E}$:

Error SS&CP Matrix
     Y1           Y2         Y3
Y1 21.106667   8.783333  -0.336667
Y2  8.783333  26.646667   0.173333
Y3 -0.336667   0.173333   0.046667

and the matrix used to test the hypothesis of no interactions, that is, $\boldsymbol\gamma_{ij}=0$ for all $i,j$:

H = Type III SS&CP Matrix for ROW*COLUMN
        Y1      Y2         Y3
Y1   0.28778   0.435     0.06
Y2   0.435     3.2233    0.13667
Y3   0.06      0.13667   0.01333

This leads to the test statistics:

S=2    M=0    N=4
Statistic  Value    F   NumDF DenDF Pr > F
Wilks'    0.6576  0.7771  6     20  0.5973
Pillai's  0.3673  0.8249  6     22  0.5629
Hot'g-L'y 0.4827  0.7241  6     18  0.6359
Roy's     0.3839  1.4077  3     11  0.2926

Conclusion: there is no evidence that interaction terms are needed. In other words the effect of changing the level of the row variable does not depend on the level of the column variable, and vice versa.

No interaction: investigate main effects.

H = Type III SS&CP Matrix for COLUMN
      Y1      Y2        Y3
Y1  0.37556  0.9533   -0.13
Y2  0.95333  2.42     -0.33
Y3 -0.13    -0.33      0.045
S=1    M=0.5    N=4
Statistic   Value     F  NumDF DenDF Pr > F
Wilks'    0.41672  4.6656  3    10   0.0275
Pillai's  0.58328  4.6656  3    10   0.0275
Hot'g-L'y 1.39968  4.6656  3    10   0.0275
Roy's     1.39968  4.6656  3    10   0.0275
H = Type III SS&CP Matrix for ROW
    Y1       Y2        Y3
Y1 4.8144  8.68944  0.37889
Y2 8.6894 32.68778  0.53556
Y3 0.3789  0.53556  0.03111
S=2    M=0    N=4
Statistic  Value    F  NumDF DenDF Pr>F
Wilks'    0.2144 3.8652   6   20 0.0101
Pillai's  1.0680 4.2019   6   22 0.0058
Hot'g-L'y 2.3465 3.5198   6   18 0.0175
Roy's     1.4169 5.1953   3   11 0.0177

Conclusions: both row and column effects appear to exist.

Rerun model without interactions.

Common tactic to increase degrees of freedom for estimation of error variance.

data mas;
    infile 'mas';
    input row column y1 y2 y3;
proc print;
proc glm;
     class row column;
     model y1-y3 = row column;
           manova h=_all_ / printh printe;
run;

Notice absence of | which means no interaction included in model.

Effect: only ${\bf E}$ changed. Pool ${\bf E}$ above with ${\bf H}$ for interaction.

E = Error SSCP Matrix
       y1       y2       y3
y1  21.39444  9.21833 -0.27667
y2   9.21833 29.87     0.31
y3  -0.27667  0.31     0.06
H = Type III SSCP Matrix for row
E = Error SSCP Matrix
S=2    M=0    N=5
Statistic    Value   F  NumDF DenDF Pr > F
Wilks'       0.2630 3.80  6    24   0.0084
Pillai's     0.9700 4.08  6    26   0.0052
Hot'g-Lawley 1.9162 3.71  6 14.353  0.0195
Roy's        1.1377 4.93  3    13   0.0168

Compare:

$$\left[\begin{array}{rrr} 21.39 & 9.22 & -0.28 \\ 9.22 & 29.87 & 0.31 \\ -0.28 & 0.31 & 0.06 \end{array}\right] =$$

$$\left[\begin{array}{rrr} 0.29 & 0.44 & 0.06 \\ 0.44 & 3.22 & 0.... ...8& -0.34 \\ 8.78 & 26.65 & 0.17 \\ -0.34 & 0.17 & 0.05 \end{array}\right]$$

Tests for new model use pooled ${\bf E}$. Further topics in MANOVA:

1) In 1 way manova there is a likelihood ratio test for the hypothesis of homogeneous variances against the alternative of a different ${\boldsymbol\Sigma}$ in each group.

2) Univariate analyses have advantage of increased power; price is increased assumptions.

3) In SAS proc glm the statements means may be used to generate multiple comparisons procedures. SAS implements many such (Bonferroni, Scheffé, Tukey, ...

4) The random statement causes production of a table of expected mean squares.

Profile Analysis

Return to repeated measures analysis of 1 way layout.

Plot variable mean versus variable level for each group.

Put group mean vectors in rows of 4$\times$3 matrix tab5.7means then:

plot(c(1,1,1,1,2,2,2,2,3,3,3,3),tab5.7means,
     type='n',xaxt='n',xlab="",ylab="Mean")
 lines(1:3,tab5.7means[1,])
 lines(1:3,tab5.7means[2,])
 lines(1:3,tab5.7means[3,])
 lines(1:3,tab5.7means[4,])
axis(side=1,at=1:3,labels=c("A","B","C"))

Tested hypothesis of parallel profiles by testing for group by scale interaction using repeated statement.

If accepted use same output to test no group effect.

Recall output from SAS

Univariate Tests of Hypotheses for
          Within Subject Effects
Source: SCALE*GROUP
DF TypeIII   MS    F   Pr > F G - G  H - F
6   18.9619 3.160 2.05 0.0860 0.0889 0.0860

As before: weak evidence of non-parallel profiles. Notice adjustments in univariate test designed to improve $F$ approximations.

Test for no group effect assuming no interaction: model is

$${\boldsymbol\mu}_i-{\boldsymbol\mu}_j = c_j {\bf 1}^T$$

Do one way anova on sum of the three scores to test hypothesis all $c_j=0$.

Repeated Measures Analysis of Variance
Tests of Hypotheses: Between Subjects Effects
Source DF Type III SS Mean Square  F    Pr > F
GROUP   3  743.900000  247.966667 70.93 0.0001
Error  17   59.433333    3.496078

New data set (to match sas total variables, divide by $\sqrt{3}$:

1 32.908965343808667
1 34.641016151377542
1 36.373066958946424
1 33.48631561299829
   et cetera
4 18.475208614068023
4 17.897858344878397
4 17.320508075688771

Sas code:

    data profile;
    infile 'prof.dat';
    input group  average;
run;
proc glm;
    class group;
    model average = group;
run;

This gives the output

Dependent Variable: average
              Sum of
Source    DF Squares Mean Sq   F    Pr > F
Model      3 743.90  247.9667 70.93 <.0001
Error     17  59.43    3.4961
Corrd Tot 20 803.33

Notice agreement with SAS output. Conclusion: profiles are credibly parallel in the 4 groups but not co-incident. Notice small sample size. Had they been co-incident we might test hypothesis of constant profiles: no scale effect.