title | author | date |
---|---|---|
MultivariaTeach |
Ben Goren |
2023 |
MultivariaTeach is designed to help students of multivariate statistics understand and perform a selection of fundamental tests. Where many statistical packages provide "helpful" features that automatically perform various functions or hide "useless" information, the goal of this package is to provide a selection of simple tools which can be easily put together as desired. For example, rather than use a command to indicate you wish to perform a repeated measures analysis analysis of variance (ANOVA), there is a tool to create a matrix of polynomial contrasts which you can then use in the rest of the analysis.
In addition to Python, you will need the standard numpy
,
pandas
, and scipy
packages installed (which is almost
certainly already the case). Almost all of the mathematical
calculations are done with numpy
; the pandas
library is used
in the initial transformation of data into the necessary format;
and scipy
is used to calculate sklearn
package for its datasets
.
It is assumed that you are a graduate-level student in a multivariate statistics course. The math won't make much sense to those without the necessary academic prerequisites; and, to those who have completed such a course, the more common statistical analysis software packages (such as SAS or R) are almost certainly more practical.
To install the multivariateach
package, simply run the following
command:
pip install multivariateach
Before detailing how to use MultivariaTeach, it should be noted that the math it implements (described below) is heavily influenced by and modeled on the documentation provided online by SAS.
However, that documentation is not entirely consistent. Within a single section, there are no conflicts, but the same operations are not uncommonly described with different notation or logic across sections.
Below, you will find references to the SAS documentation, including where conflicts can be found and the particular implementation used by MultivariaTeach.
As much as possible, notation in this document and naming conventions in the code follow those provided by SAS in its documentation.
For those familiar with SAS, MultivariTeach provides much (but
certainly not all!) of the basic MANOVA (and related)
functionality of PROC GLM
; that is, of testing the null
hypothesis
In most software packages, performing a MANOVA test is done by
loading all the data into an object and then writing out a textual
command mimicking mathematical syntax to indicate the test you
wish to perform. For example, consider Ronald Fisher's famous 1936
iris flower data set, with the species in the first column named
group
and the four measurements in the remaining columns,
labeled x1
, x2
, x3
, and x4
. In Python, you can easily load
the data as follows:
import pandas as pd
from sklearn import datasets
iris = datasets.load_iris()
data = pd.DataFrame(np.hstack([iris.target.reshape(-1, 1), iris.data]))
data.columns = ['group', 'x1', 'x2', 'x3', 'x4']
Elsewhere, you might use something in the spirit of manova data=data, test='group ~ x1 + x2 + x3 + x4'
. The software then
"helpfully" does various transformations to the data and reports
the various test statistics.
Here, instead, you must yourself prepare the
In our model,
While you are welcome to construct all four matrices by hand, MultivariaTeach provides tools (described below) to help construct them.
When you perform the MANOVA, the function first finds
$$
\mathbf{b} = (\mathbf{X}^\intercal \mathbf{X})^- \mathbf{X}^\intercal \mathbf{Y},
$$
where
Then, the function calculates:
$$
\mathbf{H} = \mathbf{M}^\intercal (\mathbf{L}\mathbf{b})^\intercal(\mathbf{L}(\mathbf{X}^\intercal\mathbf{X})^-\mathbf{L}^\intercal)^{-1}(\mathbf{L}\mathbf{b})\mathbf{M},
$$
where
Next, the error SSCP matrix
Note that these calculations match those in the SAS documentation of the Multivariate Analysis of Variance for the GLM Procedure. However, in the SAS documentation of Multivariate Tests in its Introduction to Regression Procedures, these matrices are given as a general case using weighted regression capable of testing a null hypothesis where the right-hand-side is non-zero. Of course, when the weights are equal and the right-hand-size is zero, both calculations are equivalent.
First, as usual, make sure you load the package at the top of your Python file:
import multivariateach as mt
This import
statement would typically be included with the lines
above importing pandas
and sklearn
. We can now create our
X = mt.create_design_matrix(data, 'group')
Y = mt.create_response_matrix(data, ['x1', 'x2', 'x3', 'x4'])
You can inspect these two matrices, if you wish, before performing the analysis.
For this example, we conduct a "Type III" hypothesis test as
described in the Type III and IV SS and Estimable
Functions,
section of the SAS manual on The Four Types of Estimable
Functions.
This choice of hypothesis dictates our construction of
Creating such a matrix with MultivariaTeach is trivial:
L = mt.create_type_iii_hypothesis_contrast(X)
In many MANOVA tests, there is no need to examine transformations
of variables, in which case one uses numpy
, another library which you
should import at the top of your file:
import numpy as np
M = np.eye(Y.shape[1])
With all the building blocks assembled, you can now perform the MANOVA:
results = mt.run_manova(X, Y, L, M)
Unlike with most other software, you are not presented with a
pretty-formatted table of results; instead, the results
object
can now be itself queried for \ldots\ well, for anything and
everything. So, for example, a simple-but-complete test might look
like this:
import numpy as np
import pandas as pd
from sklearn import datasets
import multivariateach as mt
iris = datasets.load_iris() # Fisher's 1936 flower data
# Extract the parts we want in the structure we need
data = pd.DataFrame(np.hstack([iris.target.reshape(-1, 1), iris.data]))
# Give names to the columns
data.columns = ['group', 'x1', 'x2', 'x3', 'x4']
X = mt.create_design_matrix(data, 'group')
Y = mt.create_response_matrix(data, ['x1', 'x2', 'x3', 'x4'])
L = mt.create_type_iii_hypothesis_contrast(X)
M = np.eye(Y.shape[1])
results = mt.run_manova(X, Y, L, M)
print("L:\n", L)
print("M:\n", M)
print("b (aka beta hat):\n", results.b)
print("E (error SSCP):\n", results.E)
print("H (hypothesis SSCP):\n", results.H)
print(
f"Wilks's :\t{results.wilks_lambda.statistic}\n"
f"\tF stat:\t{results.wilks_lambda.F}\n"
f"\tNum DF:\t{results.wilks_lambda.df_n}\n"
f"\tDen DF:\t{results.wilks_lambda.df_d}\n"
f"\tPr > F:\t{results.wilks_lambda.p_value}"
)
print(
f"Pillai's :\t{results.pillais_trace.statistic}\n"
f"\tF stat:\t{results.pillais_trace.F}\n"
f"\tNum DF:\t{results.pillais_trace.df_n}\n"
f"\tDen DF:\t{results.pillais_trace.df_d}\n"
f"\tPr > F:\t{results.pillais_trace.p_value}"
)
print(
f"H-L Trace:\t{results.hotelling_lawley_trace.statistic}\n"
f"\tF stat:\t{results.hotelling_lawley_trace.F}\n"
f"\tNum DF:\t{results.hotelling_lawley_trace.df_n}\n"
f"\tDen DF:\t{results.hotelling_lawley_trace.df_d}\n"
f"\tPr > F:\t{results.hotelling_lawley_trace.p_value}"
)
print(
f"Roy's Root:\t{results.roys_largest_root.statistic}\n"
f"\tF stat:\t{results.roys_largest_root.F}\n"
f"\tNum DF:\t{results.roys_largest_root.df_n}\n"
f"\tDen DF:\t{results.roys_largest_root.df_d}\n"
f"\tPr > F:\t{results.roys_largest_root.p_value}"
)
This will create the following output:
L:
[[ 0. -1. 1. 0.]
[ 0. 0. -1. 1.]]
M:
[[1. 0. 0. 0.]
[0. 1. 0. 0.]
[0. 0. 1. 0.]
[0. 0. 0. 1.]]
b (aka beta hat):
[[ 4.3825 2.293 2.8185 0.8995]
[ 0.6235 1.135 -1.3565 -0.6535]
[ 1.5535 0.477 1.4415 0.4265]
[ 2.2055 0.681 2.7335 1.1265]]
E (error SSCP):
[[38.9562 13.63 24.6246 5.645 ]
[13.63 16.962 8.1208 4.8084]
[24.6246 8.1208 27.2226 6.2718]
[ 5.645 4.8084 6.2718 6.1566]]
H (hypothesis SSCP):
[[ 63.21213333 -19.95266667 165.2484 71.27933333]
[-19.95266667 11.34493333 -57.2396 -22.93266667]
[165.2484 -57.2396 437.1028 186.774 ]
[ 71.27933333 -22.93266667 186.774 80.41333333]]
Wilks's : 0.023438630650878298
F stat: 199.14534354008438
Num DF: 8
Den DF: 288.0
Pr > F: 1.3650058325896826e-112
Pillai's : 1.1918988250414653
F stat: 53.46648878461304
Num DF: 8.0
Den DF: 290.0
Pr > F: 9.74216271943989e-53
H-L Trace: 32.47732024090117
F stat: 582.1970181105928
Num DF: 8
Den DF: 203.40239043824732
Pr > F: 1.0774199831263415e-135
Roy's Root: 32.19192919827811
F stat: 1166.9574334375816
Num DF: 4
Den DF: 145
Pr > F: 3.787297649634471e-109
With luck, by now the code above and its results are obvious and self-explanatory.
Also note that the
To reproduce these results in SAS, one of course needs the data. The following Python code will output the data in a format suitable to copy and paste into SAS:
from sklearn import datasets
iris = datasets.load_iris()
data = iris.data
target = iris.target
target_names = iris.target_names
# Convert numerical target values to their corresponding string labels
species = [target_names[i] for i in target]
# Combine data and species into a single list of tuples
combined_data = list(zip(data, species))
# Print the combined data in a SAS-compatible format
for row in combined_data:
print(
"{:.1f} {:.1f} {:.1f} {:.1f} {}".format(
row[0][0], row[0][1], row[0][2], row[0][3], row[1]
)
)
Then, the analysis may be performed in SAS by replacing the placeholder data below with the full set:
data iris;
input Sepal_Length Sepal_Width Petal_Length Petal_Width Species $;
datalines;
5.1 3.5 1.4 0.2 Setosa
4.9 3.0 1.4 0.2 Setosa
4.7 3.2 1.3 0.2 Setosa
...
7.4 2.8 6.1 1.9 Virginica
7.9 3.8 6.4 2.0 Virginica
6.4 2.8 5.6 2.2 Virginica
;
run;
proc glm data=iris;
class Species;
model Sepal_Length Sepal_Width Petal_Length Petal_Width = Species / nouni;
manova h=Species;
run;
quit;
The test statistics above are calculated as described in the SAS manual section on Multivariate Tests mentioned above. Note again that this section documents the general case where one may test hypotheses where the right-hand-side is non-zero with weighted regression; this is not supported by MultivariaTeach.
The following variables are used by some or all of the calculations of the various statistics:
-
$p$ : the rank of$\mathbf{H}+\mathbf{E}$ ; -
$q$ : the rank of$\mathbf{L}(\mathbf{X^\intercal \mathbf{X})^- \mathbf{L}^\intercal}$ ; -
$v$ : the error degrees of freedom, calculated as the total observations less the rank of$\mathbf{X}$ ; -
$s$ : the lesser of$p$ and$q$ ; -
$m = \frac{\lvert p-q \rvert - 1}{2}$ ; and -
$n = \frac{v-p-1}{2}$ .
The statistic is calculated as:
$$
\lambda = \frac{ \lvert \mathbf{E} \rvert }{ \lvert \mathbf{E} + \mathbf{H} \rvert }.
$$
Then, let
If
$$
V = \operatorname{tr}\left(\mathbf{H}{(\mathbf{H}+\mathbf{E})^{-1}}\right),
$$
then
$$
F = \frac{2n+s+1}{2m+s+1} \cdot \frac{V}{s-V}
$$
has
Let $$ U = \operatorname{tr}(\mathbf{E}^{-1}\mathbf{H}). $$
Now, for
Otherwise, when
Last, let
One can use MANOVA with an
The SAS manual section on Repeated Measures Analysis of
Variance
provides as an example five levels of a drug corresponding to 1,
2, 5, 10, and 20 milligrams administered to different groups; that
same
levels = np.array([1, 2, 5, 10, 20])
degree = 5
M = mt.create_orthopolynomial_contrasts(levels, degree)
Note that the SAS documentation gives the transpose of
print(M.T)
which outputs:
[[-0.425 -0.3606 -0.1674 0.1545 0.7984]
[ 0.4349 0.2073 -0.3252 -0.7116 0.3946]
[-0.4331 0.1366 0.7253 -0.5108 0.0821]
[-0.4926 0.78 -0.3744 0.0936 -0.0066]]
(Also note that there are differences in some of the signs of the elements of the matrix; these differences are inconsequential and represent different orientations of equivalent orthonormal bases.)
The iris data, of course, is not a repeated measures experiment,
so performing such a test with the data would not make sense;
however, the following repeated measures example is from Littell,
Pendergast, and Natarajan (2000)
of an examination of
the effects of two drugs (a
and c
) and a placebo (p
), with
"FEV1" measurements of respiratory ability taken over eight
one-hour intervals after treatment. The first column is the
patient ID, which is not of interest in this analysis. The second
is of a baseline FEV1 measurement made before administering the
drug; one could include it in an analysis, but we do not do so
here. The next eight columns are the after-administration
measurements; and the last is the drug.
We load the data and perform the analysis as above with the Iris
data, but with M
as polynomial contrasts rather than the
identity matrix:
import numpy as np
import pandas as pd
import multivariateach as mt
raw_data = [
(201, 2.46, 2.68, 2.76, 2.50, 2.30, 2.14, 2.40, 2.33, 2.20, 'a'),
(202, 3.50, 3.95, 3.65, 2.93, 2.53, 3.04, 3.37, 3.14, 2.62, 'a'),
(203, 1.96, 2.28, 2.34, 2.29, 2.43, 2.06, 2.18, 2.28, 2.29, 'a'),
(204, 3.44, 4.08, 3.87, 3.79, 3.30, 3.80, 3.24, 2.98, 2.91, 'a'),
(205, 2.80, 4.09, 3.90, 3.54, 3.35, 3.15, 3.23, 3.46, 3.27, 'a'),
(206, 2.36, 3.79, 3.97, 3.78, 3.69, 3.31, 2.83, 2.72, 3.00, 'a'),
(207, 1.77, 3.82, 3.44, 3.46, 3.02, 2.98, 3.10, 2.79, 2.88, 'a'),
(208, 2.64, 3.67, 3.47, 3.19, 2.19, 2.85, 2.68, 2.60, 2.73, 'a'),
(209, 2.30, 4.12, 3.71, 3.57, 3.49, 3.64, 3.38, 2.28, 3.72, 'a'),
(210, 2.27, 2.77, 2.77, 2.75, 2.75, 2.71, 2.75, 2.52, 2.60, 'a'),
(211, 2.44, 3.77, 3.73, 3.67, 3.56, 3.59, 3.35, 3.32, 3.18, 'a'),
(212, 2.04, 2.00, 1.91, 1.88, 2.09, 2.08, 1.98, 1.70, 1.40, 'a'),
(214, 2.77, 3.36, 3.42, 3.28, 3.30, 3.31, 2.99, 3.01, 3.08, 'a'),
(215, 2.96, 4.31, 4.02, 3.38, 3.31, 3.46, 3.49, 3.38, 3.35, 'a'),
(216, 3.11, 3.88, 3.92, 3.71, 3.59, 3.57, 3.48, 3.42, 3.63, 'a'),
(217, 1.47, 1.97, 1.90, 1.45, 1.45, 1.24, 1.24, 1.17, 1.27, 'a'),
(218, 2.73, 2.91, 2.99, 2.87, 2.88, 2.84, 2.67, 2.69, 2.77, 'a'),
(219, 3.25, 3.59, 3.54, 3.17, 2.92, 3.48, 3.05, 3.27, 2.96, 'a'),
(220, 2.73, 2.88, 3.06, 2.75, 2.71, 2.83, 2.58, 2.68, 2.42, 'a'),
(221, 3.30, 4.04, 3.94, 3.84, 3.99, 3.90, 3.89, 3.89, 2.98, 'a'),
(222, 2.85, 3.38, 3.42, 3.28, 2.94, 2.96, 3.12, 2.98, 2.99, 'a'),
(223, 2.72, 4.49, 4.35, 4.38, 4.36, 3.77, 4.23, 3.83, 3.89, 'a'),
(224, 3.68, 4.17, 4.30, 4.16, 4.07, 3.87, 3.87, 3.85, 3.82, 'a'),
(232, 2.49, 3.73, 3.51, 3.16, 3.26, 3.07, 2.77, 2.92, 3.00, 'a'),
(201, 2.30, 3.41, 3.48, 3.41, 3.49, 3.33, 3.20, 3.07, 3.15, 'c'),
(202, 2.91, 3.92, 4.02, 4.04, 3.64, 3.29, 3.10, 2.70, 2.69, 'c'),
(203, 2.08, 2.52, 2.44, 2.27, 2.23, 2.01, 2.26, 2.34, 2.44, 'c'),
(204, 3.02, 4.43, 4.30, 4.08, 4.01, 3.62, 3.23, 2.46, 2.97, 'c'),
(205, 3.26, 4.55, 4.58, 4.44, 4.04, 4.33, 3.87, 3.75, 3.81, 'c'),
(206, 2.29, 4.25, 4.37, 4.10, 4.20, 3.84, 3.43, 3.79, 3.74, 'c'),
(207, 1.96, 3.00, 2.80, 2.59, 2.42, 1.61, 1.83, 1.21, 1.50, 'c'),
(208, 2.70, 4.06, 3.98, 4.06, 3.93, 3.61, 2.91, 2.07, 2.67, 'c'),
(209, 2.50, 4.37, 4.06, 3.68, 3.64, 3.17, 3.37, 3.20, 3.25, 'c'),
(210, 2.35, 2.83, 2.79, 2.82, 2.79, 2.80, 2.76, 2.64, 2.69, 'c'),
(211, 2.34, 4.06, 3.68, 3.59, 3.27, 2.60, 2.72, 2.22, 2.68, 'c'),
(212, 2.20, 2.82, 1.90, 2.57, 2.30, 1.67, 1.90, 2.07, 1.76, 'c'),
(214, 2.78, 3.18, 3.13, 3.11, 2.97, 3.06, 3.27, 3.24, 3.33, 'c'),
(215, 3.43, 4.39, 4.63, 4.19, 4.00, 4.01, 3.66, 3.47, 3.22, 'c'),
(216, 3.07, 3.90, 3.98, 4.09, 4.03, 4.07, 3.56, 3.83, 3.75, 'c'),
(217, 1.21, 2.31, 2.19, 2.21, 2.09, 1.75, 1.72, 1.80, 1.36, 'c'),
(218, 2.60, 3.19, 3.18, 3.15, 3.14, 3.08, 2.96, 2.97, 2.85, 'c'),
(219, 2.61, 3.54, 3.45, 3.25, 3.01, 3.07, 2.65, 2.47, 2.55, 'c'),
(220, 2.48, 2.99, 3.02, 3.02, 2.94, 2.69, 2.66, 2.68, 2.70, 'c'),
(221, 3.73, 4.37, 4.20, 4.17, 4.19, 4.07, 3.86, 3.89, 3.89, 'c'),
(222, 2.54, 3.26, 3.39, 3.27, 3.20, 3.32, 3.09, 3.25, 3.15, 'c'),
(223, 2.83, 4.72, 4.97, 4.99, 4.96, 4.95, 4.82, 4.56, 4.49, 'c'),
(224, 3.47, 4.27, 4.50, 4.34, 4.00, 4.11, 3.93, 3.68, 3.77, 'c'),
(232, 2.79, 4.10, 3.85, 4.27, 4.01, 3.78, 3.14, 3.94, 3.69, 'c'),
(201, 2.14, 2.36, 2.36, 2.28, 2.35, 2.31, 2.62, 2.12, 2.42, 'p'),
(202, 3.37, 3.03, 3.02, 3.19, 2.98, 3.01, 2.75, 2.70, 2.84, 'p'),
(203, 1.88, 1.99, 1.62, 1.65, 1.68, 1.65, 1.85, 1.96, 1.30, 'p'),
(204, 3.10, 3.24, 3.37, 3.54, 3.31, 2.81, 3.58, 3.76, 3.05, 'p'),
(205, 2.91, 3.35, 3.92, 3.69, 3.97, 3.94, 3.63, 2.92, 3.31, 'p'),
(206, 2.29, 3.04, 3.28, 3.17, 2.99, 3.31, 3.21, 2.98, 2.82, 'p'),
(207, 2.20, 2.46, 3.22, 2.65, 3.02, 2.25, 1.50, 2.37, 1.94, 'p'),
(208, 2.70, 2.85, 2.81, 2.96, 2.69, 2.18, 1.91, 2.21, 1.71, 'p'),
(209, 2.25, 3.45, 3.48, 3.80, 3.60, 2.83, 3.17, 3.22, 3.13, 'p'),
(210, 2.48, 2.56, 2.52, 2.67, 2.60, 2.68, 2.64, 2.65, 2.61, 'p'),
(211, 2.12, 2.19, 2.44, 2.41, 2.55, 2.93, 3.08, 3.11, 3.06, 'p'),
(212, 2.37, 2.14, 1.92, 1.75, 1.58, 1.51, 1.94, 1.84, 1.76, 'p'),
(214, 2.73, 2.57, 3.08, 2.62, 2.91, 2.71, 2.39, 2.42, 2.73, 'p'),
(215, 3.15, 2.90, 2.80, 3.17, 2.39, 3.01, 3.22, 2.75, 3.14, 'p'),
(216, 2.52, 3.02, 3.21, 3.17, 3.13, 3.38, 3.25, 3.29, 3.35, 'p'),
(217, 1.48, 1.35, 1.15, 1.24, 1.32, 0.95, 1.24, 1.04, 1.16, 'p'),
(218, 2.52, 2.61, 2.59, 2.77, 2.73, 2.70, 2.72, 2.71, 2.75, 'p'),
(219, 2.90, 2.91, 2.89, 3.01, 2.74, 2.71, 2.86, 2.95, 2.66, 'p'),
(220, 2.83, 2.78, 2.89, 2.77, 2.77, 2.69, 2.65, 2.84, 2.80, 'p'),
(221, 3.50, 3.81, 3.77, 3.78, 3.90, 3.80, 3.78, 3.70, 3.61, 'p'),
(222, 2.86, 3.06, 2.95, 3.07, 3.10, 2.67, 2.68, 2.94, 2.89, 'p'),
(223, 2.42, 2.87, 3.08, 3.02, 3.14, 3.67, 3.84, 3.55, 3.75, 'p'),
(224, 3.66, 3.98, 3.77, 3.65, 3.81, 3.77, 3.89, 3.63, 3.74, 'p'),
(232, 2.88, 3.04, 3.00, 3.24, 3.37, 2.69, 2.89, 2.89, 2.76, 'p'),
]
columns = [
'patient', 'basefev1',
'fev11h', 'fev12h', 'fev13h', 'fev14h',
'fev15h', 'fev16h', 'fev17h', 'fev18h',
'drug'
]
fev1 = pd.DataFrame(raw_data, columns=columns)
X = mt.create_design_matrix(fev1, 'drug')
Y = mt.create_response_matrix(fev1, [
'fev11h', 'fev12h', 'fev13h', 'fev14h',
'fev15h', 'fev16h', 'fev17h', 'fev18h'
])
L = mt.create_type_iii_hypothesis_contrast(X)
levels = np.array(np.arange(Y.shape[1])) + 1
degree = Y.shape[1] - 1
M = mt.create_orthopolynomial_contrasts(levels, degree)
results = mt.run_manova(X, Y, L, M)
print("L:\n", L)
print("M:\n", M)
print("b (aka beta hat):\n", results.b)
print("E (error SSCP):\n", results.E)
print("H (hypothesis SSCP):\n", results.H)
print(
f"Wilks's :\t{results.wilks_lambda.statistic}\n"
f"\tF stat:\t{results.wilks_lambda.F}\n"
f"\tNum DF:\t{results.wilks_lambda.df_n}\n"
f"\tDen DF:\t{results.wilks_lambda.df_d}\n"
f"\tPr > F:\t{results.wilks_lambda.p_value}"
)
print(
f"Pillai's :\t{results.pillais_trace.statistic}\n"
f"\tF stat:\t{results.pillais_trace.F}\n"
f"\tNum DF:\t{results.pillais_trace.df_n}\n"
f"\tDen DF:\t{results.pillais_trace.df_d}\n"
f"\tPr > F:\t{results.pillais_trace.p_value}"
)
print(
f"H-L Trace:\t{results.hotelling_lawley_trace.statistic}\n"
f"\tF stat:\t{results.hotelling_lawley_trace.F}\n"
f"\tNum DF:\t{results.hotelling_lawley_trace.df_n}\n"
f"\tDen DF:\t{results.hotelling_lawley_trace.df_d}\n"
f"\tPr > F:\t{results.hotelling_lawley_trace.p_value}"
)
print(
f"Roy's Root:\t{results.roys_largest_root.statistic}\n"
f"\tF stat:\t{results.roys_largest_root.F}\n"
f"\tNum DF:\t{results.roys_largest_root.df_n}\n"
f"\tDen DF:\t{results.roys_largest_root.df_d}\n"
f"\tPr > F:\t{results.roys_largest_root.p_value}"
)
The results, which correspond to the "hour*drug" output of SAS, are as follows:
L:
[[ 0. -1. 1. 0.]
[ 0. 0. -1. 1.]]
M:
[[-0.54006172 0.54006172 0.43082022 -0.28203804 0.14978617 0.06154575
0.01706972]
[-0.38575837 0.07715167 -0.30772873 0.52378493 -0.49215457 -0.30772873
-0.11948803]
[-0.23145502 -0.23145502 -0.43082022 0.12087344 0.36376642 0.55391171
0.35846409]
[-0.07715167 -0.38575837 -0.18463724 -0.36262033 0.32097037 -0.30772873
-0.59744015]
[ 0.07715167 -0.38575837 0.18463724 -0.36262033 -0.32097037 -0.30772873
0.59744015]
[ 0.23145502 -0.23145502 0.43082022 0.12087344 -0.36376642 0.55391171
-0.35846409]
[ 0.38575837 0.07715167 0.30772873 0.52378493 0.49215457 -0.30772873
0.11948803]
[ 0.54006172 0.54006172 -0.43082022 -0.28203804 -0.14978617 0.06154575
-0.01706972]]
b (aka beta hat):
[[2.4971875 2.47833333 2.41416667 2.3396875 2.2671875 2.219375
2.156875 2.14947917]
[0.9915625 0.93375 0.785 0.72197917 0.8015625 0.77520833
0.726875 0.72385417]
[1.1878125 1.14208333 1.15708333 1.0978125 0.97614583 0.85979167
0.81395833 0.8546875 ]
[0.3178125 0.4025 0.47208333 0.51989583 0.48947917 0.584375
0.61604167 0.5709375 ]]
E (error SSCP):
[[13.53008299 0.03729573 2.95577082 1.24673897 -0.207543 -0.13586661
0.89153884]
[ 0.03729573 3.3719559 0.63292572 0.65439853 0.09069842 0.79937725
-0.07909713]
[ 2.95577082 0.63292572 3.38663539 1.28422136 -0.23182537 0.68165963
0.46884371]
[ 1.24673897 0.65439853 1.28422136 2.65050559 1.08485336 -0.32215038
0.0957733 ]
[-0.207543 0.09069842 -0.23182537 1.08485336 3.32298731 -0.44344811
-0.7321604 ]
[-0.13586661 0.79937725 0.68165963 -0.32215038 -0.44344811 2.34518955
-0.2131844 ]
[ 0.89153884 -0.07909713 0.46884371 0.0957733 -0.7321604 -0.2131844
1.88290265]]
H (hypothesis SSCP):
[[ 5.08137267 -0.82869484 0.53019893 0.7140945 0.28970756 0.1449495
-0.40105116]
[-0.82869484 0.40594737 0.19859451 -0.00684876 -0.28693657 -0.05620211
0.14645219]
[ 0.53019893 0.19859451 0.355397 0.18989178 -0.22208473 -0.01915376
0.04346894]
[ 0.7140945 -0.00684876 0.18989178 0.14471854 -0.05630392 0.00718976
-0.02355585]
[ 0.28970756 -0.28693657 -0.22208473 -0.05630392 0.22867081 0.03708622
-0.09460136]
[ 0.1449495 -0.05620211 -0.01915376 0.00718976 0.03708622 0.00805041
-0.02118594]
[-0.40105116 0.14645219 0.04346894 -0.02355585 -0.09460136 -0.02118594
0.05590952]]
Wilks's : 0.5119131905306363
F stat: 3.578948797420759
Num DF: 14
Den DF: 126.0
Pr > F: 5.667506977740117e-05
Pillai's : 0.554909797309596
F stat: 3.5108265176304614
Num DF: 14.0
Den DF: 128.0
Pr > F: 7.191270361016422e-05
H-L Trace: 0.8229204275679236
F stat: 3.660676350971247
Num DF: 14
Den DF: 97.49520766773166
Pr > F: 6.871575915707377e-05
Roy's Root: 0.6083452743990918
F stat: 5.562013937363125
Num DF: 7
Den DF: 64
Pr > F: 4.939855356602607e-05
All four statistics are significant by any reasonable measure; the null hypothesis should be rejected.
The perform_box_m_test
function in MultivariaTeach can perform
Box's M test for the homogeneity of covariance matrices. It takes
To perform the test, use code such as:
BoxM = mt.perform_box_m_test(X, Y)
print(f"Box's M: {BoxM.statistic}; Chisq: {BoxM.chi2}; DF: {BoxM.df}; Pr > Chisq: {BoxM.p_value}")
In the case of Fisher's Iris data, this will return:
Box's M: 146.6632492125118; Chisq: 140.94304992349774; DF: 20.0; Pr > Chisq: 0.0
Note again that MultivariaTeach returns the result of Python's
computation; in reporting these results, you should indicate that
the calculated
To calculate Box's M Test, we first, with
$\mathbf{S}{\operatorname{pooled}}$ as the pooled sample
covariance matrix, $\mathbf{S}\ell$ the
The calculation is from page 310 of Johnson & Wichern's text.
Mauchly's test is used in repeated measures to test for equal
variances among the differences between all possible pairings of
the independent variables, aka "sphericity". To perform the test
with MultivariaTeach, you only need supply the mauchly
function; it will return
the same type of
To perform Mauchly's test:
Mauchly = mt.mauchly(X, Y)
print(f"Mauchly's W: {Mauchly.statistic}; Chisq: {Mauchly.chi2}; DF: {Mauchly.df}; Pr > Chisq: {Mauchly.p_value}")
When performed on the FEV1 data above, the result is:
Mauchly's W: 0.0654899236469594; Chisq: 181.13981973529806; DF: 27.0; Pr > Chisq: 0.0
The vanishingly-small
Let
If Mauchly's test indicates asphericity, you may wish to perform
univariate tests with the Greenhouse-Geisser correction. The
correction itself can be calculated with
mt.greenhouse_geisser_correction(Y, M)
, which will return the
value of
For example:
epsilon = mt.greenhouse_geisser_correction(X, Y, M)
print(f"Greenhouse-Geisser's epsilon: {epsilon}")
With the FEV1 data, this produces:
Greenhouse-Geisser's epsilon: 0.49705792649184893
With perfect sphericity,
MultivariaTeach calculates calculation_method="Original"
to the
greenhouse_geisser_correction()
function.
The function calculate_pooled_covariance_matrix(X, Y)
will
return
S_p = mt.calculate_pooled_covariance_matrix(X, Y)
print("S_pooled:\n", S_p)
The output:
S_pooled:
[[0.265 0.0927 0.1675 0.0384]
[0.0927 0.1154 0.0552 0.0327]
[0.1675 0.0552 0.1852 0.0427]
[0.0384 0.0327 0.0427 0.0419]]
For each of the four main MANOVA test statistics, you can directly
provide the
These functions would be of most interest if you are supplied not with data to analyze but instead with the above information and are asked to calculate the corresponding statistic. However, in such a circumstance, the math represented in the code may well be of more interest. In any case, the reader is strongly encouraged to examine the code directly should there be a need to directly perform these calculations.
The MultivariaTeach
package includes a suite of tests to ensure
its functionality. To run the tests, you'll need to have pytest
installed. If you don't have it installed, you can do so with the
following command:
pip install pytest
After installing pytest
, navigate to the directory containing
the MultivariaTeach
package, and run the following command:
pytest multivariateach
This will execute the included tests and provide you with the results.
Implementation of additional multivariate statistical methods would be most welcome. Please send an email with suggestions.
This project is licensed under the ISC License. The full license text can be found in the LICENSE file in the repository.
Thanks to ChatGPT, who has been a surprisingly good coach and even copilot. Though we've chased each other down rabbit holes, almost every rabbit hole has led to a deeper understanding.