Hotelling test from ICSNP

HotellingsT2(cbind(xs, ys, zs) ~ group)

    Hotelling's two sample T2-test

data:  cbind(xs, ys, zs) by group
T.2 = 2.7408, df1 = 3, df2 = 696, p-value = 0.04244
alternative hypothesis: true location difference is not equal to c(0,0,0)

hotelling.test from Hotelling package

Previously, the Hotelling package had a function called hotel.test that worked differently from HotellingsT2. In recent versions this has been replaced by hotelling.test, that permits specifying the dependent measures either as a matrix or with \(+\) syntax on the left of the ~ sign:

## From Hotelling package
print(hotelling.test(cbind(xs, ys, zs) ~ group))  ##This won't work in earlier versions of hotel.test.

Test stat:  8.2459 
Numerator df:  3 
Denominator df:  696 
P-value:  0.04244 

print(hotelling.test(xs + ys + zs ~ group))  ##correct; value is the same as ICSNP

Test stat:  8.2459 
Numerator df:  3 
Denominator df:  696 
P-value:  0.04244 

Notice that the sign and scale of the components don’t matter– z is negatively related to the others, but if we reverse code it, the results are the same: ##Negatively code z

z2 <- 5 - zs

print(hotelling.test(xs + ys + z2 ~ group))

Test stat:  8.2459 
Numerator df:  3 
Denominator df:  696 
P-value:  0.04244 

print(HotellingsT2(cbind(xs, ys, z2) ~ group))

    Hotelling's two sample T2-test

data:  cbind(xs, ys, z2) by group
T.2 = 2.7408, df1 = 3, df2 = 696, p-value = 0.04244
alternative hypothesis: true location difference is not equal to c(0,0,0)

Similarly, if we reduce scale of one of the components (y), everything stays the same:

y2 <- ys/10

print(hotelling.test(xs + y2 + z2 ~ group))

Test stat:  8.2459 
Numerator df:  3 
Denominator df:  696 
P-value:  0.04244 

print(HotellingsT2(cbind(xs, y2, z2) ~ group))

    Hotelling's two sample T2-test

data:  cbind(xs, y2, z2) by group
T.2 = 2.7408, df1 = 3, df2 = 696, p-value = 0.04244
alternative hypothesis: true location difference is not equal to c(0,0,0)

So, Hotelling’s \(T^2\) provides a way of creating sort of a composite score, without worrying about the scale or coding direction of the components. Notice that each individual component did not differ significantly across groups, but the set of three did. This also protects against multiple testing–if we’d measured 10 things and expected them all to differ a little but none did, we’d still have a high chance of a type-II false alarm error. Here, we do just one test. However, notice that the \(T^2\) was not as strong as our “theory-based” composite score above, so we could do better if we have a principled method for combining values (and we were right).

Adding garbage does not help!

xx <- rnorm(700)
yy <- rnorm(700)
print(hotelling.test(xs + y2 + z2 + xx + yy ~ group))

Test stat:  9.3857 
Numerator df:  5 
Denominator df:  694 
P-value:  0.09805 

Roy’s greatest characteristic root

• Tests for differences on only the first discriminant function
• Most appropriate when DVs are strongly interrelated on a single dimension
• Highly sensitive to violation of assumptions - most powerful when all assumptions are met

Wilks’ lambda (\(\Lambda\))

• Most commonly used statistic for overall significance
• Considers differences over all the characteristic roots
• The smaller the value of Wilks’ lambda, the larger the between-groups dispersion

Hotelling’s trace

• Considers differences over all the characteristic roots

Pillai’s criterion

• Considers differences over all the characteristic roots
• More robust than Wilks’; should be used when sample size decreases, unequal cell sizes or homogeneity of covariances is violated

Typically, packages will report Pillai’s criterion or Pillai’s trace. It is the sum of \((\lambda)/(1-\lambda)\), which should remind you of the effect size \(f^2\) used in calculating power, since the eigenvalue is a measure of the proportion of variance. Hotelling’s Trace (or Hotelling-Lawley’s Trace) is the sum of the eigenvalues; Wilks \(\lambda\) is the sum of \(1/(1+\lambda_i)\). Roy’s greatest root is the largest \(\lambda\); essentially considering just the first principle component.

Each of these tests are available for specifying in Anova, but “Pillai” is the default:

Manova(lm1, test.statistic = "Wilks")

Type II MANOVA Tests: Wilks test statistic
      Df test stat approx F num Df den Df   Pr(>F)   
val    1   0.97832   5.4857      2    495 0.004402 **
group  2   0.96413   4.5619      4    990 0.001180 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Manova(lm1, test.statistic = "Hotelling-Lawley")

Type II MANOVA Tests: Hotelling-Lawley test statistic
      Df test stat approx F num Df den Df   Pr(>F)   
val    1  0.022164   5.4857      2    495 0.004402 **
group  2  0.037162   4.5895      4    988 0.001124 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Manova(lm1, test.statistic = "Roy")

Type II MANOVA Tests: Roy test statistic
      Df test stat approx F num Df den Df    Pr(>F)    
val    1  0.022164   5.4857      2    495 0.0044017 ** 
group  2  0.036008   8.9300      2    496 0.0001549 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The different statistics give different interpretations. Roy’s greatest root essentially asks whether there is a difference on the primary dimension of shared variance of the dependent variables. The others consider all the dimensions of shared variance, but combine these differently.

Ultimately, the MANOVA often hides more than it reveals. Many statisticians have recommended you avoid using it unless you really need to, but it can come in handy, and understanding is important for evaluating other research.

Conceptual relationship to PCA

Conceptually, MANOVA is sort of like doing something like the following, which performs a regression on dimensions of a principal components analysis. If you are trying to do something like this, it would usually be better to use the MANOVA directly rather than trying to figure out something ad hoc like predicting eigen vectors ‘by hand’:

E <- eigen(cor((dat1[, 3:4])))
print(E)

eigen() decomposition
$values
[1] 1.0331709 0.9668291

$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

p1 <- as.matrix(dat1[, 3:4]) %*% as.matrix(E$vectors[, 1])
p2 <- as.matrix(dat1[, 3:4]) %*% as.matrix(E$vectors[, 2])

m1 <- (lm(p1 ~ val + group))
summary(m1)

Call:
lm(formula = p1 ~ val + group)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.5333  -2.4216  -0.1097   2.7163  11.3580 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 49.84928    0.31788 156.818   <2e-16 ***
val         -0.09538    0.18372  -0.519    0.604    
groupB      -0.33967    0.44495  -0.763    0.446    
groupC      -0.04186    0.44686  -0.094    0.925    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.057 on 496 degrees of freedom
Multiple R-squared:  0.001911,  Adjusted R-squared:  -0.004126 
F-statistic: 0.3165 on 3 and 496 DF,  p-value: 0.8134

Anova(m1)

Anova Table (Type II tests)

Response: p1
          Sum Sq  Df F value Pr(>F)
val          4.4   1  0.2695 0.6039
group       11.5   2  0.3500 0.7049
Residuals 8162.5 496               

m2 <- lm(p2 ~ val + group)
summary(m2)

Call:
lm(formula = p2 ~ val + group)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.4104  -2.4849  -0.2255   2.5036  13.4315 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.0150     0.2990  23.458  < 2e-16 ***
val          -0.5658     0.1728  -3.274  0.00113 ** 
groupB        1.3304     0.4186   3.178  0.00157 ** 
groupC        1.6798     0.4204   3.996 7.42e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.816 on 496 degrees of freedom
Multiple R-squared:  0.05608,   Adjusted R-squared:  0.05037 
F-statistic: 9.822 on 3 and 496 DF,  p-value: 2.649e-06

Anova(m2)

Anova Table (Type II tests)

Response: p2
          Sum Sq  Df F value    Pr(>F)    
val        156.1   1 10.7188 0.0011346 ** 
group      258.3   2  8.8685 0.0001643 ***
Residuals 7223.8 496                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notice that for two groups, the eigenvectors are essentially the mean and the difference between the observations. As the number of outcomes increases, these would be less constrained. If you look at the F tests for each predictor, you will see that the values are very similar to those produced in the PCA version (either on average, or max values, etc.). But this means that, depending on the test statistic, a MANOVA might also be considering strange and non-intuitive combinations of variables. Importantly, the group of dependent measures do not need to be correlated. Consider the following dependent variables. Here, out1 and out2 are independent, but group is related to the difference between them. Group is thus also indirectly related to each one, because the difference is larger when one is large or when the other is small (or both). Group2 and out3 are unrelated to the others.

Some statisticians have explored procedures like this for doing ANOVA on PCA (e.g., https://iopscience.iop.org/article/10.1088/1742-6596/1899/1/012103/pdf). According to Abapihi et al. (2021), practitioners sometimes have found the MANOVA difficult to apply because of its strict assumptions, and propose ANOVA on the first principal component as an alternative. It seems like this would be successful if your first dimension is highly interpretable, but I wouldn’t recommend it unless you take the PCA seriously first. For example, if your first component is highly interpretable and accounts for a substantial amount of variance, this might be easily justified, but it is mostly a similar procedure as Roy’s GCR statistic. Other MANOVA statistics will incorporate all dimensions.

set.seed(150)
yout1 <- rnorm(150)
yout2 <- rnorm(150)
yout3 <- rnorm(150)
delta <- ((yout1 - mean(yout1)) - (yout2 - mean(yout2))) + rnorm(150) * 5
group1 <- floor(rank(delta)/50) + 1
group2 <- sample(1:5, replace = T, size = 150)
out <- cbind(yout1, yout2, yout3)
summary(lm(out ~ group1))

Response yout1 :

Call:
lm(formula = yout1 ~ group1)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.41516 -0.71119  0.00642  0.75209  2.45372 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.21468    0.21687   -0.99    0.324
group1       0.14901    0.09933    1.50    0.136

Residual standard error: 1.008 on 148 degrees of freedom
Multiple R-squared:  0.01498,   Adjusted R-squared:  0.008322 
F-statistic:  2.25 on 1 and 148 DF,  p-value: 0.1357


Response yout2 :

Call:
lm(formula = yout2 ~ group1)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.80585 -0.52989 -0.04297  0.60846  2.24696 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.40930    0.20202   2.026   0.0446 *
group1      -0.23560    0.09253  -2.546   0.0119 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9388 on 148 degrees of freedom
Multiple R-squared:  0.04197,   Adjusted R-squared:  0.0355 
F-statistic: 6.484 on 1 and 148 DF,  p-value: 0.01191


Response yout3 :

Call:
lm(formula = yout3 ~ group1)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.99510 -0.70386 -0.03446  0.70520  2.32555 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.05673    0.21535   0.263    0.793
group1      -0.03119    0.09864  -0.316    0.752

Residual standard error: 1.001 on 148 degrees of freedom
Multiple R-squared:  0.0006752, Adjusted R-squared:  -0.006077 
F-statistic: 0.09999 on 1 and 148 DF,  p-value: 0.7523

Manova(lm(out ~ group1 + group2), test.statistic = "Roy")

Type II MANOVA Tests: Roy test statistic
       Df test stat approx F num Df den Df  Pr(>F)  
group1  1  0.053897  2.60501      3    145 0.05412 .
group2  1  0.016187  0.78238      3    145 0.50559  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Manova(lm(out ~ group1 + group2), test.statistic = "Pillai")

Type II MANOVA Tests: Pillai test statistic
       Df test stat approx F num Df den Df  Pr(>F)  
group1  1  0.051140  2.60501      3    145 0.05412 .
group2  1  0.015929  0.78238      3    145 0.50559  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Manova(lm(out ~ group1 + group2), test.statistic = "Wilks")

Type II MANOVA Tests: Wilks test statistic
       Df test stat approx F num Df den Df  Pr(>F)  
group1  1   0.94886  2.60501      3    145 0.05412 .
group2  1   0.98407  0.78238      3    145 0.50559  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Relation to Repeated Measures ANOVA

This all seems like it is a bit like a repeated-measures ANOVA. By putting the data in ‘long’ format, we can just combine out1 and out2 into a single DV, then predict the mean difference between the two with a categorical value:

dat2 <- data.frame(val, group, out = c(out1, out2), measure = rep(c("out1", "out2"),
    500))


model.rm <- lm(out ~ val + group + measure, data = dat2)
summary(model.rm)

Call:
lm(formula = out ~ val + group + measure, data = dat2)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.3432  -5.8142   0.1021   5.6072  18.0186 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 35.38334    0.44659  79.231   <2e-16 ***
val         -0.06344    0.22260  -0.285    0.776    
groupB      -0.25003    0.53912  -0.464    0.643    
groupC      -0.03096    0.54119  -0.057    0.954    
measureout2 -0.26154    0.43988  -0.595    0.552    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.948 on 995 degrees of freedom
Multiple R-squared:  0.0006802, Adjusted R-squared:  -0.003337 
F-statistic: 0.1693 on 4 and 995 DF,  p-value: 0.9541

This might be OK is some situations. But this is assuming that the two measures differ just by a constant. What if they have different scales, different variances, different directional relationships? Just like the t-test above, if we had a theory that could be used to combine these, it might be a better approach, but in this case, we are restricted in the inferences we can make. Notice that for each sub-model, the effects are in different directions, and so cancel out one another if we do it this way. We could sort of do this by predicting all the interactions too, but this begins to get very complicated such that a MANOVA might be preferred.

Relationship to Logistic Regression, PCA, LDA, and Classification

As we can see, although MANOVA seems like it is just a simple extension of ANOVA, it relates to many multi-variate concepts. For Hotelling’s \(T^2\) version, we are attempting to find whether there is a difference between two groups on multiple measures. We could easily turn this around; predicting group membership based on those multiple `predictors’. This is what methods such as logistic regression and linear discriminant analysis do. As we extend to multiple categories in MANOVA, this becomes similar to what some machine classification schemes permit–classifying categories based on multiple observed ‘features’. Finally, we can understand MANOVA as sort of doing regressions on each principle component, and determining a reasonable means of combining those regressions at the end of the day.

library(psych)
library(car)
data(bfi)

dplyr::select(bfi, A1:O5)

      A1 A2 A3 A4 A5 C1 C2 C3 C4 C5 E1 E2 E3 E4 E5 N1 N2 N3 N4 N5 O1 O2 O3 O4
61617  2  4  3  4  4  2  3  3  4  4  3  3  3  4  4  3  4  2  2  3  3  6  3  4
61618  2  4  5  2  5  5  4  4  3  4  1  1  6  4  3  3  3  3  5  5  4  2  4  3
61620  5  4  5  4  4  4  5  4  2  5  2  4  4  4  5  4  5  4  2  3  4  2  5  5
      O5
61617  3
61618  3
61620  2
 [ reached 'max' / getOption("max.print") -- omitted 2797 rows ]

lm1 <- lm(as.matrix(dplyr::select(bfi, A1:O5)) ~ bfi$gender)
Manova(lm1)

Type II MANOVA Tests: Pillai test statistic
           Df test stat approx F num Df den Df    Pr(>F)    
bfi$gender  1   0.14813   16.763     25   2410 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Packages and resources

• https://stats.idre.ucla.edu/r/dae/canonical-correlation-analysis/
• cancor function in stats package
• CCA package
• candisc package

library(CCA)
matcor(out, cbind(group1, group2))

$Xcor
            yout1       yout2       yout3
yout1  1.00000000 -0.06856134 -0.02984807
yout2 -0.06856134  1.00000000 -0.03009671
yout3 -0.02984807 -0.03009671  1.00000000

$Ycor
          group1    group2
group1 1.0000000 0.1070147
group2 0.1070147 1.0000000

$XYcor
             yout1       yout2       yout3      group1      group2
yout1   1.00000000 -0.06856134 -0.02984807  0.12238444 -0.03493933
yout2  -0.06856134  1.00000000 -0.03009671 -0.20486423 -0.12148177
yout3  -0.02984807 -0.03009671  1.00000000 -0.02598397 -0.04842062
group1  0.12238444 -0.20486423 -0.02598397  1.00000000  0.10701468
group2 -0.03493933 -0.12148177 -0.04842062  0.10701468  1.00000000

img.matcor(matcor(out, cbind(group1, group2)))

img.matcor(matcor(out, cbind(group1, group2)), type = 2)

x <- cc(out, cbind(group1, group2))
print(x)

$cor
[1] 0.24555116 0.09975194

$names
$names$Xnames
[1] "yout1" "yout2" "yout3"

$names$Ynames
[1] "group1" "group2"

$names$ind.names
NULL


$xcoef
            [,1]       [,2]
yout1 -0.3316122 -0.8923137
yout2  0.9470580 -0.3420840
yout3  0.1802280 -0.3832381

$ycoef
             [,1]       [,2]
group1 -1.0892190 -0.5270819
group2 -0.2316294  0.6476470

$scores
$scores$xscores
              [,1]         [,2]
  [1,]  1.36305021  0.995883466
  [2,] -0.56565672  0.507227219
  [3,] -0.91210575  1.338943099
  [4,]  0.04287964  0.109404897
  [5,] -0.40624108  0.911386088
  [6,] -0.94537667  0.595946887
  [7,]  0.05491903 -0.853656228
  [8,]  1.80673300 -0.959494532
  [9,]  0.01749876 -0.286410996
 [10,] -0.35742517  0.791743077
 [11,]  0.50436696 -0.109395258
 [12,]  0.45791621  1.929404728
 [13,]  1.84623811 -0.978375490
 [14,] -0.13690027  1.415658235
 [15,] -0.48291423  0.965121629
 [16,] -1.56541826 -0.654126001
 [17,]  0.12005283  0.220043991
 [18,] -0.81618654 -1.359155663
 [19,] -1.60095541  0.393913477
 [20,]  0.83146701  1.820576387
 [21,]  0.92091739 -1.155230595
 [22,] -0.68469957 -1.096972683
 [23,] -0.25309080 -0.714380944
 [24,]  0.57425159 -0.109118293
 [25,] -1.05845240 -2.193440426
 [26,]  0.19111584 -0.911822474
 [27,]  0.03437887 -1.641895806
 [28,]  0.20492095  1.615580213
 [29,] -0.30816833 -1.963921730
 [30,] -0.45679786  0.930771364
 [31,] -1.73536337 -0.351629856
 [32,]  0.55155071  2.275276390
 [33,] -0.41845943 -0.373586254
 [34,] -0.10720517 -1.095527309
 [35,] -0.55227934  0.573454033
 [36,]  1.25759802 -1.056976195
 [37,]  0.62647125 -0.771836960
 [ reached getOption("max.print") -- omitted 113 rows ]

$scores$yscores
               [,1]        [,2]
  [1,]  1.313293056 -0.02798816
  [2,]  1.313293056 -0.02798816
  [3,] -0.633515672 -1.72979909
  [4,] -1.096774378 -0.43450501
  [5,] -1.328403731  0.21314202
  [6,] -0.470814043  1.38787100
  [7,] -0.633515672 -1.72979909
  [8,] -0.239184690  0.74022397
  [9,] -0.633515672 -1.72979909
 [10,]  1.081663703  0.61965887
 [11,]  0.455703369 -1.20271714
 [12,]  0.224074016 -0.55507011
 [13,]  1.081663703  0.61965887
 [14,]  1.544922409 -0.67563520
 [15,]  1.081663703  0.61965887
 [16,] -0.633515672 -1.72979909
 [17,] -1.560033084  0.86078906
 [18,] -1.560033084  0.86078906
 [19,] -1.328403731  0.21314202
 [20,] -0.239184690  0.74022397
 [21,] -0.007555337  0.09257693
 [22,] -1.096774378 -0.43450501
 [23,]  0.850034350  1.26730591
 [24,]  0.618404997  1.91495295
 [25,] -2.417622771 -0.31393992
 [26,]  0.455703369 -1.20271714
 [27,]  1.544922409 -0.67563520
 [28,] -1.096774378 -0.43450501
 [29,]  0.224074016 -0.55507011
 [30,] -0.865145025 -1.08215205
 [31,]  1.313293056 -0.02798816
 [32,]  1.081663703  0.61965887
 [33,] -0.470814043  1.38787100
 [34,] -1.328403731  0.21314202
 [35,]  1.544922409 -0.67563520
 [36,] -0.007555337  0.09257693
 [37,] -0.007555337  0.09257693
 [ reached getOption("max.print") -- omitted 113 rows ]

$scores$corr.X.xscores
            [,1]       [,2]
yout1 -0.4030366 -0.8692130
yout2  0.9228995 -0.2535793
yout3  0.1625920 -0.3455781

$scores$corr.Y.xscores
             [,1]        [,2]
group1 -0.2312088 -0.03359224
group2 -0.1069591  0.08979132

$scores$corr.X.yscores
             [,1]        [,2]
yout1 -0.09896611 -0.08670568
yout2  0.22661905 -0.02529503
yout3  0.03992466 -0.03447209

$scores$corr.Y.yscores
             [,1]       [,2]
group1 -0.9415913 -0.3367578
group2 -0.4355880  0.9001461

plot(x$scores$xscores, x$scores$yscores)

Here, we see the canonical correlation uses two dimensions, because group only has two variables. The first dimension finds a correlation between groups of .24, and the second .09.
The first dimension correlates negatively with out1, positively with out2, and not with out3. The second dimension correlates with out3 only. On the other side, group 1 and group 2 correlate (negatively) with dimension 1, and group1 is positively related to dimension 2 but group2 is negatively related to dimension 2.

We can sort of look at what is going on by plotting the first dimension xscores and yscores for each person—this is how each person maps into each of the first canonical correlation dimensions. We can see the correlation is .24, which is a small but probably significant relationship.

cor(x$scores$xscores[, 1], x$scores$yscores[, 1])

[1] 0.2455512

plot(x$scores$xscores[, 1], x$scores$yscores[, 1])