The data are stored in the file that can be downloaded: nutrient.csv

options ls=78;
title "Hotellings T2 - Women's Nutrition Data";

data nutrient;
  infile "D:\Statistics\STAT 505\data\nutrient.csv" firstobs=2 delimiter=',';
  input id calcium iron protein a c;
  run;
  
  /* iml code to compute the hotelling t2 statistic
  * hotel is the name of the module we define here
  * mu0 is the null vector
  * one is a vector of 1s
  * ident is the identity matrix
  * ybar is the vector of sample means
  * s is the sample covariance matrix
  * t2 is the squared statistical distance between ybar and mu0
  * f is the final form of the t2 statistic after scaling
  * to have an f-distribution
  * the module definition is ended with the 'finish' statement
  * use nutrient makes the data set 'nutrient' available
  * the variables from nutrient are input to x and hotel module is called
  */

proc iml;
  start hotel;
    mu0={1000, 15, 60, 800, 75};
    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{calcium iron protein a c} into x;
  run hotel;

A one-at-a-time 95% confidence interval for calcium is given by the following where values are substituted into the formula and the calculations are worked out as shown below:

\(624.04925 \pm t_{736, 0.025}\sqrt{157829.4/737}\)

\(624.04925 \pm 1.96 \times 14.63390424\)

\(624.04925 \pm 28.68245\)

\((595.3668, 652.7317)\)

The one-at-a-time confidence intervals are summarized in the table below:

Looking at this table, it appears that the average daily intakes of calcium and iron are too low (because the intervals fall below the recommended intakes of these variables), and the average daily intake of protein is too high (because the interval falls above the recommended intake of protein).

Problem: The problem with these one-at-a-time intervals is that they do not control for a family-wide error rate.

Consequence: We are less than 95% confident that all of the intervals simultaneously cover their respective means.

To fix this problem we can calculate a \((1 - \alpha) × 100\%\) Confidence Ellipse for the population mean vector \(\boldsymbol{\mu}\). To calculate this confidence ellipse we must recall that for independent random observations from a multivariate normal distribution with mean vector \(\mu\) and variance-covariance matrix \(Σ\), the F-statistic, (shown below), is F-distributed with p and n-p degrees of freedom:

\(F = n\mathbf{(\bar{x} - \boldsymbol{\mu})'S^{-1}(\bar{x} - \boldsymbol{\mu})}\dfrac{(n-p)}{p(n-1)} \sim F_{p,n-p}\)

This next expression says that the probability that n times the squared Mahalanobis distance between the sample mean vector, \(\boldsymbol{\bar{x}}\), and the population mean vector \(\boldsymbol{\mu}\) is less than or equal to p times n-1 times the critical value from the F-table divided by n-p is equal to \(1 - α\).

\(\text{Pr}\{n\mathbf{(\bar{x} - \boldsymbol{\mu})'S^{-1}(\bar{x} - \boldsymbol{\mu})} \le \dfrac{p(n-1)}{(n-p)}F_{p,n-p,\alpha}\} = 1-\alpha\)

Here the squared Mahalanobis distance between \(\bar{\mathbf{x}}\) and \(\boldsymbol{\mu}\) is being used.

In particular, this is the \((1 - \alpha) × 100%\) confidence ellipse for the population mean, \(\boldsymbol{\mu}\).

Interpretation

The \((1 - \alpha) × 100\%\) confidence ellipse is very similar to the prediction ellipse that we discussed earlier in our discussion of the multivariate normal distribution. A \((1 - \alpha) × 100\%\) confidence ellipse yields simultaneous \((1 - \alpha) × 100\%\) confidence intervals for all linear combinations of the variable means. Consider linear combinations of population means as below:

\(c_1\mu_1 + c_2\mu_2 + \dots c_p \mu_p = \sum_{j=1}^{p}c_j\mu_j = \mathbf{c'\boldsymbol{\mu}}\)

The simultaneous \((1 - \alpha) × 100\%\) confidence intervals for the above are given by the expression below

\(\sum_{j=1}^{p}c_j\bar{x}_j \pm \sqrt{\frac{p(n-1)}{(n-p)}F_{p, n-p, \alpha}}\sqrt{\frac{1}{n}\sum_{j=1}^{p}\sum_{k=1}^{p}c_jc_ks_{jk}}\)

In terms of interpreting the \((1 - \alpha) × 100\%\) confidence ellipse, we can say that we are \((1 - \alpha) × 100%\) confident that all such confidence intervals cover their respective linear combinations of the treatment means, regardless of what linear combinations we may wish to consider. In particular, we can consider the trivial linear combinations which correspond to the individual variables. So this says that we going to be also \((1 - \alpha) × 100\%\) confident that all of the intervals given in the expression below:

\(\bar{x}_j \pm \sqrt{\frac{p(n-1)}{(n-p)}F_{p, n-p, \alpha}}\sqrt{\frac{s^2_j}{n}}\)

cover their respective treatment population means. These intervals are called simultaneous confidence intervals.

Simultaneous confidence intervals are computed using hand calculations. For calcium, we substituted the following values: The sample mean was 624.04925. We have p = 5 variables, a total sample size of n = 737, and if we look up in the F-table for 5 and 732 degrees of freedom for alpha = 0.05, the critical value is 2.21. The standard error of the sample mean for calcium is equal to the square root of the sample variance for calcium, 157,829.44, divided by the sample size, 737. The math is carried out to obtain an interval running from 575.27 to approximately 672.83 as shown below:

\(624.04925 \pm \sqrt{\frac{5(737-1)}{737-5}\underset{2.21}{\underbrace{F_{5,737-5,0.05}}}}\sqrt{157829.4/737}\)

\(624.04925 \pm 3.333224 \times 14.63390424\)

\(624.04925 \pm 48.77808\)

\((575.27117, 672.82733)\)

options ls=78;
title "Confidence Intervals - Women's Nutrition Data";

/* %let allows the p variable to be used throughout the code below
  * After reading in the nutrient 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 are in another column called 'x'.
  * This format is used for the calculations that follow.
  */

%let p=5;
data nutrient;
  infile "D:\Statistics\STAT 505\data\nutrient.csv" firstobs=2 delimiter=',';
  input id calcium iron protein a c;
  variable="calcium"; x=calcium; output;
  variable="iron";    x=iron;    output;
  variable="protein"; x=protein; output;
  variable="vit a";   x=a;       output;
  variable="vit c";   x=c;       output;
  keep variable x;
  run;

proc sort;
  by variable;
  run;

 /* The means procedure calculates and saves the sample size,
  * mean, and variance for each variable. It then saves these results 
  * in a new data set 'a' for use in the final step below.
  * /

proc means noprint;
  by variable;
  var x;
  output out=a n=n mean=xbar var=s2;
  run;

 /* The data step here is used to calculate the confidence interval
  * limits from the statistics calculated in the data set 'a'.
  * The values 't1', 'tb', and 'f' are the critical values used in the 
  * one-at-a-time, Bonferroni, and F intervals, respectively.
  * /

data b;
  set a;
  t1=tinv(1-0.025,n-1);
  tb=tinv(1-0.025/&p,n-1);
  f=finv(0.95,&p,n-&p);
  loone=xbar-t1*sqrt(s2/n);
  upone=xbar+t1*sqrt(s2/n);
  losim=xbar-sqrt(&p*(n-1)*f*s2/(n-&p)/n);
  upsim=xbar+sqrt(&p*(n-1)*f*s2/(n-&p)/n);
  lobon=xbar-tb*sqrt(s2/n);
  upbon=xbar+tb*sqrt(s2/n);
  run;

proc print data=b;
  run;

The following table gives the confidence intervals:

Looking at these simultaneous confidence intervals we can see that the upper bound of the interval for calcium falls below the recommended daily intake of 1000 mg. Similarly, the upper bound for iron also falls below the recommended daily intake of 15 mg. Conversely, the lower bound for protein falls above the recommended daily intake of 60 g of protein. The intervals for both vitamins A and C both contain the recommended daily intake for these two vitamins.

Therefore, we can conclude that the average daily intake of calcium and iron fall below the recommended levels, and the average daily intake of protein exceeds the recommended level.

Here, the 95% confidence interval for calcium under the Bonferroni correction is calculated:

\(624.04925 \pm t_{737-1, \frac{0.05}{2 \times 5}}\sqrt{157829.4/737}\)

\(624.04925 \pm 2.576 \times 14.63390424\)

\(624.04925 \pm 37.69694\)

\((586.35231, 661.74619)\)

This calculation uses the values for the sample mean for calcium, 624.04925, the critical value from the t-table, with n-1 degrees of freedom, evaluated at alpha over 2 times p, (0.05 divided 2 times 5, or .005). This critical value turns out to be 2.576 from the t-table. The standard error is calculated by taking the square root of the sample variance, (157,829.44), divided by the sample size, 737.

Carrying out the math, the interval goes from 586.35 to 661.75.

lobon=xbar-tb*sqrt(s2/n);
upbon=xbar+tb*sqrt(s2/n); 

A sample of husband and wife pairs are asked to respond to each of the following questions:

1. What is the level of passionate love you feel for your partner?
2. What is the level of passionate love your partner feels for you?
3. What is the level of companionate love you feel for your partner?
4. What is the level of companionate love your partner feels for you?

A total of 30 married couples were questioned. Responses were recorded on a five-point scale. Responses included the following values:

1. None at all
2. Very little
3. Some
4. A great deal
5. Tremendous amount

We will try to address two separate questions from these data.

1. Do the husbands respond to the questions in the same way as their wives?
2. Do the husbands and wives accurately perceive the responses of their spouses?

We shall address Question 1 first...

Question 1: Do the husbands respond to the questions in the same way as their wives?

Before considering the multivariate case let's review the univariate approach to answering this question. In this case, we will compare the responses to a single question.

Univariate Paired t-test Case: Consider comparing the responses to a single question.

A notation will be as follows:

• \(X _ { 1 i }\) = response of husband i - the first member of pair i
• \(X _ { 2 i }\) = response of wife i - the second member of pair i
• \(\mu _ { 1 }\) = population mean for the husbands - the first population
• \(\mu _ { 2 }\) = population mean for the wives - the second population

Our objective here is to test the null hypothesis that the population means are equal against the alternative hypothesis that means are not equal, as described in the expression below:

\(H_0\colon \mu_1 =\mu_2 \) against \(H_a\colon \mu_1 \ne \mu_2\)

In the univariate course, you learned that the null hypothesis is tested as follows. First, we define \(Y _ { i }\) to be the difference in responses for the \(i^{th}\) pair of observations. In this case, this will be the difference between husband i and wife i. Likewise, we can also define \(\mu _ { Y }\) to be the population mean of these differences, which is the same as the difference between the population means \(\mu _ { 1 }\) and \(\mu _ { 2 }\), both as noted below:

\(Y_i = X_{1i}-X_{2i}\) and \(\mu_Y = \mu_1-\mu_2\)

Testing the null hypothesis for the equality of the population means is going to be equivalent to testing the null hypothesis that \(\mu _ { Y }\) is equal to 0 against the general alternative that \(\mu _ { Y }\) is not equal to 0.

\(H_0\colon \mu_Y =0 \) against \(H_a\colon \mu_Y \ne 0\)

This hypothesis is tested using the paired t-test.

We will define \(\bar{y}\) to be the sample mean of the \(Y _ { i }\)'s:

\(\bar{y} = \dfrac{1}{n}\sum_{i=1}^{n}Y_i\)

We will also define \(s_{2}Y\) to be the sample variance of the \(Y _ { i }\)'s:

\(s^2_Y = \dfrac{\sum_{i=1}^{n}Y^2_i - (\sum_{i=1}^{n}Y_i)^2/n}{n-1}\)

We will make the usual four assumptions in doing this:

1. The \(Y _ { i }\)'s have common mean \(\mu _ { Y }\)
2. Homoskedasticity. The \(Y _ { i }\)'s have common variance \(\sigma^2_Y\).
3. Independence. The \(Y _ { i }\)'s are independently sampled.
4. Normality. The \(Y _ { i }\)'s are normally distributed.

The test statistic is a t-statistic which is, in this case, equal to the sample mean divided by the standard error as shown below:

\[t = \frac{\bar{y}}{\sqrt{s^2/n}} \sim t_{n-1}\]

Under the null hypothesis, \(H _ { o }\) this test statistic is going to be t-distributed with n-1 degrees of freedom and we will reject \(H _ { o }\) at level \(α\) if the absolute value of the t-value exceeds the critical value from the t-distribution with n-1 degrees of freedom evaluated at \(α\) over 2.

\(|t| > t_{n-1, \alpha/2}\)

Now let's consider the multivariate case. All scalar observations will be replaced by vectors of observations. As a result, we will use the notation that follows:

The scalar population means are replaced by the population mean vectors so that \(\mu _{ 1 }\) = population mean vector for husbands and \(\mu _{ 2 }\) = population mean vector for the wives.

Here we are interested in testing the null hypothesis that the population mean vectors are equal against the general alternative that these mean vectors are not equal.

\(H_0\colon \boldsymbol{\mu_1} = \boldsymbol{\mu_2}\) against \(H_a\colon \boldsymbol{\mu_1} \ne \boldsymbol{\mu_2}\)

Under the null hypothesis, the two mean vectors are equal element by element. As with the one-sample univariate case, we are going to look at the differences between the observations. We will define the vector \(Y _{ 1 }\) for the \(i ^{ th }\) couple to be equal to the vector \(X _{ 1i }\) for the \(i ^{ th }\) husband minus the vector \(X _{ 2i }\) for the \(i ^{ th }\) wife. Then we will also, likewise, define the vector \(\mu _{ Y }\) to be the difference between the vector \(\mu _{ 1 }\) and the vector \(\mu _{ 2 }\).

\(\mathbf{Y}_i = \mathbf{X}_{1i}-\mathbf{X}_{2i}\) and \(\boldsymbol{\mu}_Y = \boldsymbol{\mu}_1-\boldsymbol{\mu}_2\)

Testing the above null hypothesis is going to be equivalent to testing the null hypothesis that the population mean vector \(\mu _{ Y }\) is equal to 0. That is, all of its elements are equal to 0. This is tested against the alternative that the vector \(\mu _{ Y }\) is not equal to 0, that is at least one of the elements is not equal to 0.

\(H_0\colon \boldsymbol{\mu}_Y = \mathbf{0}\) against \(H_a\colon \boldsymbol{\mu}_Y \ne \mathbf{0}\)

This hypothesis is tested using the paired Hotelling's \(T^{ 2 }\) test.

As before, we will define \(\mathbf{\bar{Y}}\) to denote the sample mean vector of the vectors \(Y_{ i }\).

\(\mathbf{\bar{Y}} = \dfrac{1}{n}\sum_{i=1}^{n}\mathbf{Y}_i\)

And, we will define \(S _{ Y }\) to denote the sample variance-covariance matrix of the vectors \(Y_{ i }\).

\(\mathbf{S}_Y = \dfrac{1}{n-1}\sum_{i=1}^{n}\mathbf{(Y_i-\bar{Y})(Y_i-\bar{Y})'}\)

The assumptions are similar to the assumptions made for the one-sample Hotelling's T-square test:

1. The vectors \(Y _{ 1 }\)'s have a common population mean vector \(\mu _{ Y }\), which essentially means that there are no sub-populations with mean vectors.
2. The vectors \(Y _{ 1 }\)'s have common variance-covariance matrix \(\Sigma_Y\).
3. Independence. The \(Y _{ 1 }\)'s are independently sampled. In this case, independence among the couples in this study.
4. Normality. The \(Y _{ 1 }\)'s are multivariate normally distributed.

Paired Hotelling's T-Square test statistic is given by the expression below:

\(T^2 = n\bar{\mathbf{Y}}'\mathbf{S}_Y^{-1}\mathbf{\bar{Y}}\)

It is a function of sample size n, the sample mean vectors, \(\mathbf{\bar{Y}}\), and the inverse of the variance-covariance matrix \(\mathbf{S} _{ Y }\).

Then we will define an F-statistic as given in the expression below:

\(F = \dfrac{n-p}{p(n-1)}T^2 \sim F_{p, n-p}\)

Under the null hypothesis, \(H _{ o } \colon \mu _{ Y } = 0\), this will have an F-distribution with p and n-p degrees of freedom. We will reject \(H _{ o }\) at level \(\alpha\) if the F-value exceeds the value from F-value with p and n-p degrees of freedom, evaluated at level \(\alpha\).

\(F > F_{p,n-p,\alpha}\)

Let us return to our example...

Download the text file containing the data here: spouse.csv

options ls=78;
title "Paired Hotelling's T-Square";

 /* The differences 'd1' through 'd4' are defined
  * and added to the data set.
  */

data spouse;
  infile "D:\Statistics\STAT 505\data\spouse.csv" firstobs=2 delimiter=',';
  input h1 h2 h3 h4 w1 w2 w3 w4;
  d1=h1-w1;
  d2=h2-w2;
  d3=h3-w3;
  d4=h4-w4;
  run;

proc print data=spouse;
  run;

 /* The iml code below defines and executes the 'hotel' module 
  * for calculating the one-sample Hotelling T2 test statistic. 
  * The calculations are based on the differences, which is why 
  * there is a single null vector consisting of only 0s.
  * The commands between 'start' and 'finish' define the 
  * calculations of the module for an input vector 'x'.
  * The 'use' statement makes the 'spouse' data set available, from 
  * which the difference variables are taken. The difference variables 
  * are then read into the vector 'x' before the 'hotel' module is called.
  */

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 spouse;
  read all var{d1 d2 d3 d4} into x;
  run hotel;

options ls=78;
title "Confidence Intervals - Spouse Data";


 /* %let allows the p variable to be used throughout the code below
  * After reading in the spouse data, where each variable is
  * originally in its own column, the next statements define difference
  * variables between husbands and wives, and they stack the data
  * so that all group labels (1 through 4) are in one column called 'variable',
  * and all differences are in another column called 'diff'.
  * This format is used for the calculations that follow.
  */

%let p=4;
data spouse;
  infile "D:\Statistics\STAT 505\data\spouse.csv" firstobs=2 delimiter=',';
  input h1 h2 h3 h4 w1 w2 w3 w4;
  variable=1; diff=h1-w1; output;
  variable=2; diff=h2-w2; output;
  variable=3; diff=h3-w3; output;
  variable=4; diff=h4-w4; output;
  drop h1 h2 h3 h4 w1 w2 w3 w4;
  run;

proc sort;
  by variable;
  run;

 /* The means procedure calculates and saves the sample size,
  * mean, and variance for each variable. It then saves these results 
  * in a new data set 'a' for use in the final step below.
  */

proc means data=spouse noprint;
  by variable;
  var diff;
  output out=a n=n mean=xbar var=s2;
  run;

 /* The data step here is used to calculate the confidence interval
  * limits from the statistics calculated in the data set 'a'.
  * The values 't' and'f' are the critical values used in the 
  * Bonferroni and F intervals, respectively.
  */

data b;
  set a;
  f=finv(0.95,&p,n-&p);
  t=tinv(1-0.025/&p,n-1);
  losim=xbar-sqrt(&p*(n-1)*f*s2/(n-&p)/n); 
  upsim=xbar+sqrt(&p*(n-1)*f*s2/(n-&p)/n); 
  lobon=xbar-t*sqrt(s2/n);
  upbon=xbar+t*sqrt(s2/n);
  run;

proc print data=b;
  run;

Analysis

The simultaneous confidence intervals are plotted using Profile Plots.

options ls=78;
title "Profile Plot - Spouse Data";

 /* %let allows the p variable to be used throughout the code below
  * After reading in the spouse data, where each variable is
  * originally in its own column, the next statements define difference
  * variables between husbands and wives, and they stack the data
  * so that all group labels (1 through 4) are in one column called 'variable',
  * and all differences are in another column called 'diff'.
  * This format is used for the calculations that follow.
  */

%let p=4;
data spouse;
  infile "D:\Statistics\STAT 505\data\spouse.csv" firstobs=2 delimiter=','
  input h1 h2 h3 h4 w1 w2 w3 w4;
  variable=1; diff=h1-w1; output;
  variable=2; diff=h2-w2; output;
  variable=3; diff=h3-w3; output;
  variable=4; diff=h4-w4; output;
  drop h1 h2 h3 h4 w1 w2 w3 w4;
  run;

proc sort;
  by variable;
  run;

 /* The means procedure calculates and saves the sample size,
  * mean, and variance for each variable. It then saves these results 
  * in a new data set 'a' for use in the final step below.
  * /

proc means data=spouse;
  by variable;
  var diff;
  output out=a n=n mean=xbar var=s2;
  run;

 /* The data step here is used to calculate the simultaneous
  * confidence intervals based on the F-multiplier
  * from the statistics calculated in the data set 'a'.
  * /

data b;
  set a;
  f=finv(0.95,&p,n-&p);
  diff=xbar; output;
  diff=xbar-sqrt(&p*(n-1)*f*s2/(n-&p)/n); output;
  diff=xbar+sqrt(&p*(n-1)*f*s2/(n-&p)/n); output;
  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 confidence limits.
  * The reference line of 0 corresponds to the null value of the 
  * difference for each variable.
  * /

proc gplot data=b;
  axis1 length=4 in;
  axis2 length=6 in;
  plot diff*variable / vaxis=axis1 haxis=axis2 vref=0 lvref=21;
  symbol v=none i=hilot color=black;
  run;

Next, we will return to the second question posed earlier in this lesson.

Question 2: Do the husbands and wives accurately perceive the responses of their spouses?

To understand this question, let us return to the four questions asked of each husband and wife pair. The questions were:

1. What is the level of passionate love you feel for your partner?
2. What is the level of passionate love your partner feels for you?
3. What is the level of companionate love you feel for your partner?
4. What is the level of companionate love your partner feels for you?

Notice that these questions are all paired. The odd-numbered questions ask about how each person feels about their spouse, while the even-numbered questions ask how each person thinks their spouse feels towards them. The question that we are investigating now asks about perception, so here we are trying to see if the husbands accurately perceive the responses of their wives and conversely to the wives accurately perceive the responses of their husbands.

In more detail, we may ask:

In this case, we are asking if the wife accurately perceives the amount of passionate love her husband feels toward her.

Here, we are asking if the husband accurately perceives the amount of passionate love his wife feels toward him.

• Does the husband's answer to question 1 match the wife's answer to question 2?
• Secondly, does the wife's answer to question 1 match the husband's answer to question 2?
• Similarly, does the husband's answer to question 3 match the wife's answer to question 4?
• Does the wife's answer to question 3 match the husband's answer to question 4?

To address the research question we need to define four new variables as follows:

• \(Z_{i1} = X_{1i1} - X_{2i2}\) - (for the \(i^{th}\) couple, the husband's response to question 1 minus the wife's response to question 2.)
• \(Z_{i2} = X_{1i2} - X_{2i1}\) - (for the \(i^{th}\) couple, the husband's response to question 2 minus the wife's response to question 1.)
• \(Z_{i3} = X_{1i3} - X_{2i4}\) - (for the \(i^{th}\) couple, the husband's response to question 3 minus the wife's response to question 4.)
• \(Z_{i4} = X_{1i4} - X_{2i3}\) - (for the \(i^{th}\) couple, the husband's response to question 4 minus the wife's response to question 3.)

These Z's can then be collected into a vector. We can then calculate the sample mean for that vector...

\(\mathbf{\bar{z}} = \dfrac{1}{n}\sum_{i=1}^{n}\mathbf{Z}_i\)

as well as the sample variance-covariance matrix...

\(\mathbf{S}_Z = \dfrac{1}{n-1}\sum_{i=1}^{n}\mathbf{(Z_i-\bar{z})(Z_i-\bar{z})'}\)

Here we will make the usual assumptions about the vector \(Z_{i}\) containing the responses for the \(i^{th}\) couple:

1. The \(Z_{i}\)'s have common mean vector \(\mu_{Z}\)
2. The \(Z_{i}\)'s have common variance-covariance matrix \(\Sigma_{Z}\)
3. Independence. The \(Z_{i}\)'s are independently sampled.
4. Multivariate Normality. The \(Z_{i}\)'s are multivariate normally distributed.

Question 2 is equivalent to testing the null hypothesis that the mean \(\mu_{Z}\) is equal to 0, against the alternative that \(\mu_{Z}\) is not equal to 0 as expressed below:

\(H_0\colon \boldsymbol{\mu}_Z = \mathbf{0}\) against \(H_a\colon \boldsymbol{\mu}_Z \ne \mathbf{0}\)

We may then carry out the Paired Hotelling's T-Square test using the usual formula with the sample mean vector z-bar replacing the mean vector y-bar from our previous example, and the sample variance-covariance matrix \(S_{Z}\) replacing the sample variance-covariance matrix \(S_{Y}\) also from our previous example:

\(n\mathbf{\bar{z}'S^{-1}_Z\bar{z}}\)

We can then form the F-statistic as before:

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

And, under \(H_{o} \colon \mu_{Z} = 0\) we will reject the null hypothesis \(H_{o}\) at level \(\alpha\) if this F-statistic exceeds the critical value from the F-table with p and n-p degrees of freedom evaluated at \(α\).

\(F > F_{p, n-p, \alpha}\)

options ls=78;
title "Paired Hotelling's T-Square: Matching Perceptions";

 /* The differences 'd1' through 'd4' are defined
  * and added to the data set.
  */

data spouse;
  infile "D:\Statistics\STAT 505\data\spouse.csv" firstobs=2 delimiter=','
  input h1 h2 h3 h4 w1 w2 w3 w4;
  d1=h1-w2;
  d2=h2-w1;
  d3=h3-w4;
  d4=h4-w3;
  run;

proc print data=spouse;
  run;

 /* The iml code below defines and executes the 'hotel' module 
  * for calculating the one-sample Hotelling T2 test statistic. 
  * The calculations are based on the differences, which is why 
  * there is a single null vector consisting of only 0s.
  * The commands between 'start' and 'finish' define the 
  * calculations of the module for an input vector 'x'.
  * The 'use' statement makes the 'spouse' data set available, from 
  * which the difference variables are taken. The difference variables 
  * are then read into the vector 'x' before the 'hotel' module is called.
  */

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 spouse;
  read all var{d1 d2 d3 d4} into x;
  run hotel;

An example where we are sampling from two distinct populations occurs with 1000 franc Swiss Bank Notes.

1. The first population is the population of Genuine Bank Notes
2. The second population is the population of Counterfeit Bank Notes

While the diagram shows a more recent version issue of the 1000 franc note, it depicts the different measurement locations taken in this study. For both populations of banknotes the following measurements were taken:

• Length of the note
• Width of the Left-Hand side of the note
• Width of the Right-Hand side of the note
• Width of the Bottom Margin
• Width of the Top Margin
• Diagonal Length of Printed Area

Objective: To determine if counterfeit notes can be distinguished from genuine Swiss banknotes.

This is essential for the police if they wish to use these kinds of measurements to determine if a banknote is genuine or not and help solve counterfeiting crimes.

Before considering the multivariate case, we will first consider the univariate case...