Let us revisit the original hypothesis of interest, as below

\(H_0\colon \boldsymbol{\mu} = \boldsymbol{\mu_0}\) against \(H_a\colon \boldsymbol{\mu} \ne \boldsymbol{\mu_a}\)

Instead of testing the null hypothesis for the ratios of the means over their hypothesized means are all equal to one, profile analysis involves testing the null hypothesis that all of these ratios are equal to one another, but not necessarily equal to 1.

After rejecting

\(H_0\colon \dfrac{\mu_1}{\mu^0_1} = \dfrac{\mu_2}{\mu^0_2} = \dots = \dfrac{\mu_p}{\mu^0_p} = 1\)

we may wish to test

\(H_0\colon \dfrac{\mu_1}{\mu^0_1} = \dfrac{\mu_2}{\mu^0_2} = \dots = \dfrac{\mu_p}{\mu^0_p}\)

Profile Analysis can be carried out using the following procedure.

Step 1: Compute the differences between the successive ratios. That is we take the ratio of the j + 1-th variable over its hypothesized mean and subtract this from the ratio of jth variable over its hypothesized mean as shown below:

\(D_{ij} = \dfrac{X_{ij+1}}{\mu^0_{j+1}}-\dfrac{X_{ij}}{\mu^0_j}\)

We call this ratio \(D_{ij}\) for observation i.

\(H_0\colon \boldsymbol{\mu}_D = \mathbf{0}\)

Step 2: Apply Hotelling’s \(T^{2}\) test to the data \(D_{ij}\) to test the null hypothesis that the mean of these differences is equal to 0.

\(H_0\colon \boldsymbol{\mu}_D = \mathbf{0}\)

options ls=78;
title "Profile Analysis - Women's Nutrition Data";

 /* After reading in the data, each of the original
  * variables is divided by its null value to convert 
  * to a common scale without specific units.
  * The differences are then defined for each successive
  * pair of variables.
  */

data nutrient;
  infile "D:\Statistics\STAT 505\data\nutrient.csv" firstobs=2 delimiter=','
  input id calcium iron protein a c;
  calcium=calcium/1000;
  iron=iron/15;
  protein=protein/60;
  a=a/800;
  c=c/75;
  diff1=iron-calcium;
  diff2=protein-iron;
  diff3=a-protein;
  diff4=c-a;
  run;

 /* The iml code to compute the hotelling t2 statistic
  * is similar to that for the one-sample t2 statistic
  * except that by working with differences of variable,
  * the null values are all 0s, and the corresponding
  * null hypothesis is that all variable means are equal
  * to each other.
  */

proc iml;
  start hotel;
    mu0={0,0,0,0};
    one=j(nrow(x),1,1);
    ident=i(nrow(x));
    ybar=x`*one/nrow(x);
    s=x`*(ident-one*one`/nrow(x))*x/(nrow(x)-1.0);
    print mu0 ybar;
    print s;
    t2=nrow(x)*(ybar-mu0)`*inv(s)*(ybar-mu0);
    f=(nrow(x)-ncol(x))*t2/ncol(x)/(nrow(x)-1);
    df1=ncol(x);
    df2=nrow(x)-ncol(x);
    p=1-probf(f,df1,df2);
    print t2 f df1 df2 p;
  finish;
  use nutrient;
  read all var{diff1 diff2 diff3 diff4} into x;
  run hotel;

When the hypothesis of equality of two independent population means is rejected, one would like to know on account of which variable the hypothesis is rejected. To assess that, we go back to the example of Swiss Banknotes.

To assess which variable these notes differ on we will consider the \(( 1 - \alpha ) \times 100 \%\) Confidence Ellipse for the difference in the population mean vectors for the two populations of banknotes, \(\boldsymbol{\mu_{1}}\) - \(\boldsymbol{\mu_{2}}\). This ellipse is given by the set of \(\boldsymbol{\mu_{1}}\) satisfying the expression:

\(\mathbf{(\bar{x}_1-\bar{x}_2-(\pmb{\mu}_1-\pmb{\mu}_2))}'\left\{\mathbf{S}_p \left(\dfrac{1}{n_1}+\dfrac{1}{n_2}\right) \right\}^{-1}\mathbf{(\bar{x}_1-\bar{x}_2-(\pmb{\mu}_1-\pmb{\mu}_2))} \le \dfrac{p(n_1+n_2-2)}{n_1+n_2-p-1}F_{p,n_1+n_2-p-1,\alpha}\)

To understand the geometry of this ellipse, let

\(\lambda_1, \lambda_2, \dots, \lambda_p\)

denote the eigenvalues of the pooled variance-covariance matrix \(S_{p}\), and let

\(\mathbf{e}_{1}\), \(\mathbf{e}_{2}\), ..., \(\mathbf{e}_{p}\)

denote the corresponding eigenvectors. Then the \(k^{th}\) axis of this p dimensional ellipse points in the direction specified by the \(k^{th}\) eigenvector, \(\mathbf{e}_{k}\) And, it has a half-length given by the expression below:

\[l_k = \sqrt{\lambda_k\frac{p(n_1+n_2-2)}{n_1+n_2-p-1}\left(\frac{1}{n_1}+\frac{1}{n_2}\right)F_{p,n_1+n_2-p-1,\alpha}}\]

Note, again, that this is a function of the number of variables, p, the sample sizes \(n_{1}\) and \(n_{2}\), and the critical value from the F-table.

The \(( 1 - \alpha ) \times 100 \%\) confidence ellipse yields simultaneous \(( 1 - \alpha ) \times 100 \%\) confidence intervals for all linear combinations of the form given in the expression below:

\(c_1(\mu_{11}-\mu_{21})+c_2(\mu_{12}-\mu_{22})+\dots+c_p(\mu_{1p}-\mu_{2p}) = \sum_{k=1}^{p}c_k(\mu_{1k}-\mu_{2k}) = \mathbf{c'(\mu_1-\mu_k)}\)

So, these are all linear combinations of the differences in the sample means between the two populations where we are taking linear combinations across variables. These simultaneous confidence intervals are given by the expression below:

\(\sum_{k=1}^{p}c_k(\bar{x}_{1k}-\bar{x}_{2k}) \pm \sqrt{\dfrac{p(n_1+n_2-2)}{n_1+n_2-p-1}F_{p,n_1+n_2-p-1,\alpha}}\sqrt{\left(\dfrac{1}{n_1}+\dfrac{1}{n_2}\right) \sum_{k=1}^{p}\sum_{l=1}^{p}c_kc_ls^{(p)}_{kl}}\)

Here, the terms \(s^{(p)}_{kl}\) denote the pooled covariances between variables k and l.

Interpretation

The interpretation of these confidence intervals is the same as that for the one sample of Hotelling's T-square. Here, we are \(( 1 - \alpha ) \times 100 \%\) confident that all of these intervals cover their respective linear combinations of the differences between the means of the two populations. In particular, we are also \(( 1 - \alpha ) \times 100 \%\) confident that all of the intervals of the individual variables also cover their respective differences between the population means. For the individual variables, if we are looking at, say, the \(k^{th}\) individual variable, then we have the difference between the sample means for that variable, k, plus or minus the same radical term that we had in the expression previously, times the standard error of that difference between the sample means for the \(k^{th}\) variable. The latter involves the inverses of the sample sizes and the pooled variance for variable k.

\(\bar{x}_{1k}-\bar{x}_{2k} \pm \sqrt{\dfrac{p(n_1+n_2-2)}{n_1+n_2-p-1}F_{p,n_1+n_2-p-1,\alpha}}\sqrt{\left(\dfrac{1}{n_1}+\dfrac{1}{n_2}\right) s^2_k}\)

So, here, \(s^2_k\) is the pooled variance for variable k. These intervals are called simultaneous confidence intervals.

Let's work through an example of their calculation using the Swiss Banknotes data.

As in the one-sample case, the simultaneous confidence intervals should be computed only when we are interested in linear combinations of the variables. If the only intervals of interest, however, are the confidence intervals for the individual variables with no linear combinations, then a better approach is to calculate the Bonferroni corrected confidence intervals as given in the expression below:

\(\bar{x}_{1k} - \bar{x}_{2k} \pm t_{n_1+n_2-2, \frac{\alpha}{2p}}\sqrt{s^2_k\left(\dfrac{1}{n_1}+\frac{1}{n_2}\right)}\)

 Carrying out the math we end up with an interval that runs from 0.006 to 0.286 as shown below:

\(\bar{x}_{1k} - \bar{x}_{2k} \pm t_{2n-2, \frac{\alpha}{2p}}\sqrt{\frac{2s^2_k}{n}}\)

\(214.959 - 214.813 \pm \underset{2.665}{\underbrace{t_{2\times 100-2, \frac{0.05}{2 \times 6}}}}\sqrt{\dfrac{2 \times 0.137}{100}}\)

\((0.006, 0.286)\)

Profile plots provide another useful graphical summary of the data. These are only meaningful if all variables have the same units of measurement. They are not meaningful if the variables have different units of measurement. For example, some variables may be measured in grams while other variables are measured in centimeters. In this case, profile plots should not be constructed.

• In the traditional profile plot, the samples mean for each group is plotted against the variables.
• For the bank notes, it is preferable to subtract the government specifications before carrying out the analyses.

options ls=78;
title "Profile Plot - Swiss Bank Notes";

 /* After reading in the swiss data, where each variable is
  * originally in its own column, the next statements stack the data
  * so that all variable names are in one column called 'variable',
  * and all response values minus their null values
  * are in another column called 'x'.
  * The subtraction is useful for plotting purposes because
  * it allows for a common reference value of 0 for all variables.
  */

data swiss;
  infile "D:\Statistics\STAT 505\data\swiss3.csv" firstobs=2 delimeter=',';
  input type $ length left right bottom top diag;
  variable="length"; x=length-215.0; output;
  variable="left  "; x=left-130.0;   output;
  variable="right "; x=right-130.0;  output;
  variable="bottom"; x=bottom-9.0;   output;
  variable="top   "; x=top-10.0;     output;
  variable="diag  "; x=diag-141.354; output;
  run;

proc sort;
  by type variable;
  run;

 /* The means procedure calculates and saves mean for
  * each variable and saves the results in a new data set 'a'
  * for use in the steps below.
  * /

proc means data=swiss;
  by type variable;
  var x;
  output out=a mean=xbar;
  run;

 /* The axis commands define the size of the plotting window.
  * The horizontal axis is of the variables, and the vertical
  * axis is used for the mean values.
  * /

proc gplot;
  axis1 length=4 in label=("Mean");
  axis2 length=6 in;
  plot xbar*variable=type / vaxis=axis1 haxis=axis2;
  symbol1 v=J f=special h=2 i=join color=black;
  symbol2 v=K f=special h=2 i=join color=black;
  run;

In carrying out any statistical analysis it is always important to consider the assumptions for the analysis and confirm that all assumptions are satisfied.

Let's recall the four assumptions underlying Hotelling's T-square test.

1. The data from population i is sampled from a population with mean vector \(\boldsymbol{\mu}_{i}\).
2. The data from both populations have a common variance-covariance matrix \(Σ\)
3. Independence. The subjects from both populations are independently sampled.
   Note! This does not mean that the variables are independent of one another
4. Normality. Both populations are multivariate normally distributed.

The following will consider each of these assumptions separately, and methods for diagnosing their validity.

1. Assumption 1: The data from population i is sampled from a population mean vector \(\boldsymbol{\mu}_{i}\).

   • This assumption essentially means that there are no subpopulations with a different population mean vectors.
   • In our current example, this might be violated if the counterfeit notes were produced by more than one counterfeiter.
   • Generally, if you have randomized experiments, this assumption is not of any concern. However, in the current application, we would have to ask the police investigators whether more than one counterfeiter might be present.

2. Assumption 2: For now we will skip Assumption 2 and return to it at a later time.

3. Assumption 3: Independence

   • Says the subjects for each population were independently sampled. This does not mean that the variables are independent of one another.
   • This assumption may be violated for three different reasons:
     • Clustered data: If bank notes are produced in batches, then the data may be clustered. In this case, the notes sampled within a batch may be correlated with one another.
     • Time-series data: If the notes are produced in some order over time, there might be some temporal correlation between notes over time. The notes produced at times close to one another may be more similar. This could result in temporal correlation violating the assumptions of the analysis.
     • Spatial data: If the data were collected over space, we may encounter some spatial correlation.
     Note! the results of Hotelling's T-square are not generally robust to violations of independence.

4. Assumption 4: Multivariate Normality

   To assess this assumption we can produce the following diagnostic procedures:

   • Produce histograms for each variable. We should look for a symmetric distribution.
   • Produce scatter plots for each pair of variables. Under multivariate normality, we should see an elliptical cloud of points.
   • Produce a three-dimensional rotating scatter plot. Again, we should see an elliptical cloud of points.

So, in general, Hotelling's T-square is not going to be sensitive to violations of this assumption.

Now let us return to assumption 2.

Assumption 2. The data from both populations have a common variance-covariance matrix \(Σ\).

This assumption may be assessed using Box's Test.

Suppose that the data from population i have variance-covariance matrix \(\Sigma_i\); for population i = 1, 2. Need to test the null hypothesis that \(\Sigma_1\) is equal to \(\Sigma_2\) against the general alternative that they are not equal as shown below:

\(H_0\colon \Sigma_1 = \Sigma_2\) against \(H_a\colon \Sigma_1 \ne \Sigma_2\)

Here, the alternative is that the variance-covariance matrices differ in at least one of their elements.

The test statistic for Box's Test is given by L-prime as shown below:

\(L' = c\{(n_1+n_2-2)\log{|\mathbf{S}_p|}- (n_1-1)\log{|\mathbf{S}_1|} - (n_2-1)\log{|\mathbf{S}_2|}\}\)

This involves a finite population correction factor c, which is given below.

The finite population correction factor, c, is given below:

\(c = 1-\dfrac{2p^2+3p-1}{6(p+1)}\left\{\dfrac{1}{n_1-1}+\dfrac{1}{n_2-1} - \dfrac{1}{n_1+n_2-2}\right\}\)

It is a function of the number of variables p, and the sample sizes \(n_{1}\) and \(n_{2}\).

Under the null hypothesis, \(H_{0}\colon \Sigma_{1} = \Sigma_{2} \), Box's test statistic is approximately chi-square distributed with p(p + 1)/2 degrees of freedom. That is,

\(L' \overset{\cdot}{\sim} \chi^2_{\dfrac{p(p+1)}{2}}\)

The degrees of freedom is equal to the number of unique elements in the variance-covariance matrix (taking into account that this matrix is symmetric). We will reject \(H_o\) at level \(\alpha\) if the test statistic exceeds the critical value from the chi-square table evaluated at level \(\alpha\).

\(L' > \chi^2_{\dfrac{p(p+1)}{2}, \alpha}\)

options ls=78;
title "Bartlett's Test - Swiss Bank Notes";

data swiss;
  infile "D:\Statistics\STAT 505\data\swiss3.csv" firstobs=2 delimeter=',';
  input type $ length left right bottom top diag;
  run;

 /* The discrim procedure is called with the pool=test
  * option to produce Bartlett's test for equal 
  * covariance matrices. The remaining parts of the
  * output are not used.
  * The class statement defines the grouping variable,
  * which is the type of note, and all response
  * variables are specified  in the var statement.
  */

proc discrim pool=test;
  class type;
  var length left right bottom top diag;
  run;

The following considers a test for equality of the population mean vectors when the variance-covariance matrices are not equal.

Here we will consider the modified Hotelling's T-square test statistic given in the expression below:

\(T^2 = \mathbf{(\bar{x}_1-\bar{x}_2)}'\left\{\dfrac{1}{n_1}\mathbf{S}_1+\dfrac{1}{n_2}\mathbf{S}_2\right\}^{-1}\mathbf{(\bar{x}_1-\bar{x}_2)}\)

Again, this is a function of the differences between the sample means for the two populations. Instead of being a function of the pooled variance-covariance matrix, we can see that the modified test statistic is written as a function of the sample variance-covariance matrix, \(\mathbf{S}_{1}\), for the first population and the sample variance-covariance matrix, \(\mathbf{S}_{2}\), for the second population. It is also a function of the sample sizes \(n_{1}\) and \(n_{2}\).

For large samples, that is if both samples are large, \(T^{2}\) is approximately chi-square distributed with p d.f. We will reject \(H_{0}\): \(\boldsymbol{\mu}_{1}\) = \(\boldsymbol{\mu}_{2}\) at level \(α\) if \(T^{2}\) exceeds the critical value from the chi-square table with p d.f. evaluated at level \(α\).

\(T^2 > \chi^2_{p, \alpha}\)

For small samples, we can calculate an F transformation as before using the formula below.

\(F = \dfrac{n_1+n_2-p-1}{p(n_1+n_2-2)}T^2\textbf{ } \overset{\cdot}{\sim}\textbf{ } F_{p,\nu}\)

This formula is a function of sample sizes \(n_{1}\) and \(n_{2}\), and the number of variables p. Under the null hypothesis, this will be F-distributed with p and approximately ν degrees of freedom, where 1 divided by ν is given by the formula below:

\( \dfrac{1}{\nu} = \sum_{i=1}^{2}\frac{1}{n_i-1} \left\{ \dfrac{\mathbf{(\bar{x}_1-\bar{x}_2)}'\mathbf{S}_T^{-1}(\dfrac{1}{n_i}\mathbf{S}_i)\mathbf{S}_T^{-1}\mathbf{(\bar{x}_1-\bar{x}_2)}}{T^2} \right\} ^2 \)

This involves summing over the two samples of banknotes, a function of the number of observations of each sample, the difference in the sample mean vectors, the sample variance-covariance matrix for each of the individual samples, as well as a new matrix \(\mathbf{S}_{T}\) which is given by the expression below:

\(\mathbf{S_T} = \dfrac{1}{n_1}\mathbf{S_1} + \dfrac{1}{n_2}\mathbf{S}_2\)

We will reject \(H_{0}\) \colon \(\mu_{1}\) = \(\mu_{2}\) at level \(α\) if the F-value exceeds the critical value from the F-table with p and ν degrees of freedom evaluated at level \(α\).

\(F > F_{p,\nu, \alpha}\)

A reference for this particular test is given in Seber, G.A.F. 1984. Multivariate Observations. Wiley, New York.

options ls=78;
title "2-Sample Hotellings T2 - Swiss Bank Notes (unequal variances)";

data swiss;
  infile "D:\Statistics\STAT 505\data\swiss3.csv" firstobs=2 delimeter=',';
  input type $ length left right bottom top diag;
  run;

 /* The iml code below defines and executes the 'hotel2m' module 
  * for calculating the two-sample Hotelling T2 test statistic,
  * where here the sample covariances are not pooled.
  * The commands between 'start' and 'finish' define the 
  * calculations of the module for two input vectors 'x1' and 'x2',
  * which have the same variables but correspond to two separate groups.
  * Note that s1 and s2 are not pooled in these calculations, and the
  * resulting degrees of freedom are considerably more involved.
  * The 'use' statement makes the 'swiss' data set available, from 
  * which all the variables are taken. The variables are then read 
  * separately into the vectors 'x1' and 'x2' for each group, and
  * finally the 'hotel2' module is called.
  */

proc iml;
  start hotel2m;
    n1=nrow(x1);
    n2=nrow(x2);
    k=ncol(x1);
    one1=j(n1,1,1);
    one2=j(n2,1,1);
    ident1=i(n1);
    ident2=i(n2);
    ybar1=x1`*one1/n1;
    s1=x1`*(ident1-one1*one1`/n1)*x1/(n1-1.0);
    print n1 ybar1;
    print s1;
    ybar2=x2`*one2/n2;
    s2=x2`*(ident2-one2*one2`/n2)*x2/(n2-1.0);
    st=s1/n1+s2/n2;
    print n2 ybar2;
    print s2;
    t2=(ybar1-ybar2)`*inv(st)*(ybar1-ybar2);
    df1=k;
    p=1-probchi(t2,df1);
    print t2 df1 p;
    f=(n1+n2-k-1)*t2/k/(n1+n2-2);
    temp=((ybar1-ybar2)`*inv(st)*(s1/n1)*inv(st)*(ybar1-ybar2)/t2)**2/(n1-1);
    temp=temp+((ybar1-ybar2)`*inv(st)*(s2/n2)*inv(st)*(ybar1-ybar2)/t2)**2/(n2-1);
    df2=1/temp;
    p=1-probf(f,df1,df2);
    print f df1 df2 p;
  finish;
  use swiss;
    read all var{length left right bottom top diag} where (type="real") into x1;
    read all var{length left right bottom top diag} where (type="fake") into x2;
  run hotel2m;
  

As before, the next step is to determine how these notes differ. This may be carried out using the simultaneous \((1 - α) × 100\%\) confidence intervals.

For Large Samples: simultaneous \((1 - α) × 100\%\) confidence intervals may be calculated using the expression below:

\(\bar{x}_{1k} - \bar{x}_{2k} \pm \sqrt{\chi^2_{p,\alpha}\left(\dfrac{s^2_{1k}}{n_1} + \dfrac{s^2_{2k}}{n_2}\right)}\)

This involves the differences in the sample means for the kth variable, plus or minus the square root of the critical value from the chi-square table times the sum of the sample variances divided by their respective sample sizes.

For Small Samples: it is better to use the expression below:

\(\bar{x}_{1k} - \bar{x}_{2k} \pm \sqrt{\dfrac{p(n_1+n_2-2)}{n_1+n_2-p-1}F_{p,\nu,\alpha}}\sqrt{\left(\dfrac{s^2_{1k}}{n_1} + \dfrac{s^2_{2k}}{n_2}\right)}\)

Basically, the chi-square value and the square root are replaced by the critical value from the F-table, times a function of the number of variables, p, and the sample sizes n₁ and n₂.

In this lesson we learned about:

• Profile analysis for one-sample mean and how to carry it out in SAS
• Computing simultaneous and Bonferroni confidence intervals for the differences between sample means for each variable;
• Drawing conclusions from those confidence intervals (make sure that you indicate which population has the larger mean for each significant variable); 
• Profile plots for two-sample means
• Methods for diagnosing the assumptions of the two-sample Hotelling's T-square test.

In practice, the data analyses should proceed as follows:

1. Step 1. For small samples, use histograms, scatterplots, and rotating scatterplots to assess the multivariate normality of the data. If the data do not appear to be normally distributed, apply appropriate normalizing transformations.

2. Step 2. Use Bartlett's test to assess the assumption that the population variance-covariance matrices are homogeneous.

3. Step 3. Carry out the two-sample Hotelling's T-square test for equality of the population mean vectors. If Bartlett's test in Step 2 is significant, use the modified two-sample Hotelling's T-square test. If the two-sample Hotelling's T-square test is not significant, conclude that there is no statistically significant evidence that the two populations have different mean vectors and stop. Otherwise, go to Step 4.

4. Step 4. Compute either simultaneous or Bonferroni confidence intervals. For the significant variables, draw conclusions regarding which population has the larger mean.