Real Statistics Multivariate Functions

The following is a summary of all the multivariate functions provided in the Real Statistics Resource Pack.

These functions are organized into the following categories:

• Multivariate normal distribution
• Mahalanobis Distance and Outliers
• Hotelling T²
• One-way MANOVA
• Two-way MANOVA
• Permutational MANOVA
• Box’s Test
• Extraction
• Factor Analysis
• Cluster Analysis
• Jenks Natural Breaks
• Discriminant Analysis
• Correspondence Analysis
• Gaussian Mixture Models

Multivariate Normal Distribution

BNORMSDIST(x, y, r, cum, don) = the cdf of the standard bivariate normal distribution at x and y with correlation coefficient r if cum = TRUE and the pdf if cum = FALSE

BNORMDIST(x1, x2, m1, m2, s1, s2, r, cum, don) = the cdf of the standard bivariate normal distribution at at x1 and x2 with means m1 and m2, standard deviations s1 and s2 and correlation coefficient r if cum = TRUE and the pdf if cum = FALSE

BNORMSRECT(x1, x2, y1, y2, r, cum, don) = P(x1 < x < x2 and y1 < y < y2) for the standard bivariate normal distribution with correlation coefficient r.

BNORMRECT(x1, x2, y1, y2, mx, my, sx, sy, r, cum, don) = P(x1 < x < x2 and y1 < y < y2) for the bivariate normal distribution with means mx and my, standard deviations sx and sy and correlation coefficient r.

MNORMRAND(R1, R2): generates an array containing a random k × 1 vector from a k-dimensional multivariate normal distribution where R1 is the range that contains the k × 1 means vector and R2 is the range that contains the k × k covariance matrix

BNORMRAND(m1, m2, s1, s2, r): generates a random 2 × 1 vector from a bivariate normal distribution with means m1, m2, standard deviations s1, s2 and correlation coefficient r.

MNORMDIST(R0, Rm, Rc, cum, iter) = the cdf of the multivariate normal distribution at R0 if cum = TRUE and the pdf if cum = FALSE.

MNORMRECT(Rlo, Rhi, Rm, Rc, iter) = P(A < X < B) for the multivariate normal distribution where A and B are the column vectors corresponding to Rlo and Rhi.

MSKEWTEST(R1, lab): Mardia’s skewness test for multivariate normality; returns a column range with the values skewness, chi-square statistic, df and p-value, plus corrected statistic and p-value for small samples

MKURTTEST(R1, lab): Mardia’s kurtosis test for multivariate normality; returns a column range with the values kurtosis, z-statistic, and p-value

If lab = TRUE then an extra column of labels is appended to the results (defaults to FALSE). If don = TRUE the Donnelly algorithm is used; otherwise (default), the Genz algorithm is used.

R0, Rm, Rlo, and Rhi are k × 1 arrays or cell ranges with all the elements in Rlo less than or equal to the corresponding element in Rhi. Rm represents the mean vector and Rc is the k × k covariance matrix. iter = the number of samples used in the GHK algorithm (default 200).

Hyper-ellipse

The following functions relate to a hyper-ellipse:

CONF_MNORM(R1, alpha, iter): array function that returns a column array with the lengths of the axes in the 1-alpha confidence ellipse based on the data in R1

CONF_BNORM(R1, R2, lab, alpha): returns a column array with the lengths of the axes in the 1-alpha confidence ellipse based on the data in R1 and R2, plus theta (the angle at which the ellipse is tilted) and the percentage of data elements inside the ellipse. If lab = TRUE then a column of labels is appended to the output (default FALSE).

CONF_VOLUME(R1, alpha, iter): returns the volume (area if R1 only has two columns) of the 1-alpha confidence ellipse based on the data in R1.

The default for alpha is .05. iter = # of iterations in calculating eigenvalues (default 100).

Mahalanobis Distance and Outliers

MDistSq(R1, R2) = the Mahalanobis distance squared between the 1 × k row vector R2 and the mean vector of the sample contained in the n × k range R1

QuadForm(R1, R2, R3) = the Mahalanobis distance squared between the 1 × k row vector R1 and the 1 × k mean vector based on the 1 × k population covariance matrix in R3

MOUTLIERS(R1, alpha): when alpha = 0 or is omitted, then returns an n  × 2 array whose first column contains the Mahalanobis distance squared of each vector in R1 (i.e. the rows of the n  × k array R1) and the center of the hyper-ellipse defined by the data elements in R1 and whose second column contains the corresponding p-values based on a chi-square test. When alpha > 0 then only those values for which p-value  < alpha are returned.

Hotelling T²

T2TEST(R1, R2, type) = the p-value for Hotelling’s T² test
HotellingT2 (R1, R2, type) = T² statistic
HotellingF (R1, R2, type) = F statistic
Hotellingdf (R1, R2, type) = df₂ parameter

R1 and R2 are the ranges of the two samples, except in the one-sample case where R1 is the sample and R2 is a column vector representing the population mean column vector being tested.

The following are the acceptable values for type:

0 = one sample
1 = the samples have paired values from the same population
2 = the samples are from populations with the same covariance matrices
3 = the samples are from populations with different covariance matrices

Note that the df₁ parameter is always equal to the number of columns in R1.

Array functions

There are also the following array functions:

Hotelling(R1, R2, type, lab) – outputs a 5 × 1 column range with the values T-square, df1, df2, F, and p-value for Hotelling T2 test for the data in ranges R1 and R2.

T2_RepeatedMeasures(R1, lab) – outputs a 5 × 1 column range with the values T-square, df1, df2, F, and p-value for Hotelling T2 test for repeated measures design based on the data in ranges R1 (w/o headings).

RepMeasTransform(a,b) – outputs a ab × ab–1 array representing the transformation matrix in the case where there are two within-subjects factors A, with a levels, and B, with b levels.

If lab = TRUE then the output is a 5 × 2 range where the first column contains labels.

Power and sample size

The following functions provide support for statistical power and sample size:

Hotel1_POWER(d, n, k, ttype, alpha, iter, prec) = the statistical power for a one-sample or paired Hotelling’s T-square test where the sample size is n, the number of dependent variables is k and the effect size is d, where d = the Mahalanobis distance if ttype = 1 and d = the noncentrality parameter if ttype = 0.

Hotel2_POWER(d, n1, n2, k, ttype, alpha, iter, prec) = the statistical power for a two-sample Hotelling’s T-square test where the sample sizes are n1 and n2 and the other arguments are as for Hotel1_POWER.

Hotel1_SIZE(d, k, pow, alpha, iter, prec) = the minimum sample size required to obtain the statistical power pow for a one-sample or paired Hotelling’s T-square test where the number of dependent variables is k and the Mahalanobis effect size is d.

Hotel2_SIZE(d, k, pow, alpha, iter, prec) = the minimum size for sample 1 required to obtain statistical power of pow for a two-sample Hotelling’s T-square test where the number of dependent variables is k and the Mahalanobis effect size is d. nratio = the size of sample 2 divided by the size of sample 1 (default 1); if nratio is negative then the size of sample 2 is assumed to be –nratio.

alpha is the significance level (default .05), iter = the maximum number of iterations used in calculating the answer (default 1000) up to a precision of prec (default 0.000000001); the default for pow is .80.

One-way MANOVA

Where R1 is a data range in standard form (without column headings), we have the following functions:

Wilk’s Lambda Test

MANOVA_WilksLambda(R1) = Λ
MANOVA_Wilksdf1(R1) = df₁     
MANOVA_Wilksdf2(R1) = df₂
MANOVA_WilksF(R1) = F                    
MANOVA_WilksTest(R1) = p-value

Hotelling-Lawley Trace Test

MANOVA_HotelTrace(R1) = 
MANOVA_Hoteldf1(R1) = df₁    
MANOVA_Hoteldf2(R1) = df₂
MANOVA_HotelF(R1) = F                    
MANOVA_HotelTest(R1) = p-value

Pillai-Bartlett Trace Test

MANOVA_PillaiTrace(R1) = V
MANOVA_Pillaidf1(R1) = df₁      
MANOVA_Pillaidf2(R1) = df₂
MANOVA_PillaiF(R1) = F                     
MANOVA_PillaiTest(R1) = p-value

Roy’s Largest Root

MANOVA_RoyRoot(R1, b) = largest eigenvalue λ_p of HE^-1 if b = TRUE (default) and =  if b = FALSE

Array functions

The following are array functions that output a 5 × 1 array consisting of the specified Manova statistic, df₁, df₂, F statistic, and p-value.

MANOVA_Wilks(R1)
MANOVA_Hotel(R1)
MANOVA_Pillai(R1)

The following are array functions that calculate one of the SSCP matrices:

MANOVA_T(R1) = T                               
MANOVA_H(R1) = H                       
MANOVA_E(R1) = E

Power and sample size

The following functions provide support for statistical power and sample size:

MANOVA_POWER(f n, k, g, ttype, alpha, iter, prec) = the statistical power for one-way MANOVA where the sample size is n, the number of dependent variables is k, the number of groups is g and the effect size is f, where f = the partial eta-square effect size if ttype = 1, f = eta-square if ttype = 2 and f = Pillai’s V if ttype = 3.

MANOVA_SIZE(f, k, g, pow, ttype, alpha, iter, prec) = the minimum sample size to obtain statistical power of pow for one-way MANOVA where f, k, g and ttype are as for MANOVA_POWER.

alpha is the significance level (default .05), iter = the maximum number of iterations used in calculating the answer (default 1000) up to a precision of prec (default 0.000000001), the default for pow is .80.

Two-way MANOVA

Extract2Cov(R1, s1, s2) = covariance matrix for all the elements in R1 that match s1 and s2 in the first two columns; if s1 = “” then the pooled covariance matrix is returned.

COV2Pooled(R1) = pooled covariance matrix based on the data in R1.

Extract2Rows(R1, s1, s2, head) = an array containing all the rows in R1 that match s1 and s2 in the first two columns; if head = TRUE (default) then it is assumed that the first row in R1 contains column headings (which are included in the output)

MANOVA2_Tot(R1) = the Total (T) matrix based on the data in R1

MANOVA2_Res(R1) = the Residual (E) matrix (aka error matrix or within-subjects matrix) based on the data in R1

MANOVA2_Row(R1) = the H (hypothesis) matrix for the first factor (which we call the row factor or factor A) based on the data in R1

MANOVA2_Col(R1) = the H matrix for the second factor (which we call the column factor or factor B) based on the data in R1

MANOVA2_Int(R1) = the H matrix for the interaction between the two factors (which we call the interaction factor or factor AB) based on the data in R1

Array functions

In addition, there are the following array functions that return a 6 × 1 column array containing the values: the appropriate test statistic (as described below), df1, df2, F-stat, p-value, and partial eta-square effect size. If lab = TRUE (default FALSE) then an extra column of labels is appended to the output.

MANOVA2_WilksRow(R1, lab): uses Wilks test statistic on the column 1 factor

MANOVA2_WilksCol(R1, lab): uses Wilks test statistic on the column 2 factor

MANOVA2_WilksInt(R1, lab): uses Wilks test statistic on the interaction between the two factors

MANOVA2_PillaiRow(R1, lab): uses Pillai test statistic on the column 1 factor

MANOVA2_PillaiCol(R1, lab): uses Pillai test statistic on the column 2 factor

MANOVA2_PillaiInt(R1, lab): uses Pillai test statistic on the interaction between the two factors

MANOVA2_HotelRow(R1, lab): uses Hotelling test statistic on the column 1 factor

MANOVA2_HotelCol(R1, lab): uses Hotelling test statistic on the column 2 factor

MANOVA2_HotelInt(R1, lab): uses Hotelling test statistic on the interaction between the two factors

There are the following non-array functions. If b = TRUE (default) the largest eigenvalue, lambda, for the appropriate H covariance matrix is returned; otherwise 1/(1+lambda) is returned.

MANOVA2_RoyRow(R1, b) = Roy Root for the column 1 factor

MANOVA2_RoyCol(R1, b) = Roy Root for the column 2 factor

MANOVA2_RoyInt(R1, b) = Roy Root for the interaction factor

Non-array functions

The following array function returns a 5 × 1 column array containing the values: M-statistic df1, df2, F-stat, and p-value. If lab = TRUE (default FALSE) then an extra column of labels is appended to the output.

MANOVA2_BOX(R1, lab): performs a Box M test for equality of the covariance matrices that correspond to the interactions between the two factors in R1.

Power and sample size

The following functions provide support for statistical power and sample size:

MANOVA2Row_POWER(f, n, k, g1, g2, ttype, alpha, iter, prec) = the statistical power of the first factor (aka the row factor) for two-way MANOVA where the sample size is n.

MANOVA2Col_POWER(f, n, k, g1, g2, ttype, alpha, iter, prec) = the statistical power of the second factor (aka the column factor) for two-way MANOVA

MANOVA2Int_POWER(f, n, k, g1, g2, ttype, alpha, iter, prec) = the statistical power of the interaction between the two factors for two-way MANOVA

MANOVA2Row_SIZE(f, k, g1, g2, pow, ttype, alpha, iter, prec) = the minimum sample size to obtain statistical power of pow of the first factor for two-way MANOVA.

MANOVA2Col_SIZE(f, k, g1, g2, pow, ttype, alpha, iter, prec) = the minimum sample size to obtain statistical power of pow of the second factor for two-way MANOVA.

MANOVA2Int_SIZE(f, k, g1, g2, pow, ttype, alpha, iter, prec) = the minimum sample size to obtain statistical power of pow for the interaction between the two factors for two-way MANOVA.

The number of dependent variables is k, the number of groups in the first factor is g1, the number of groups in the second factor is g2 and the effect size is f, where f = the partial eta-square effect size if ttype = 1, f = eta-square if ttype = 2 and f = Pillai’s V if ttype = 3.

alpha is the significance level (default .05), iter = the maximum number of iterations used in calculating the answer (default 1000) up to a precision of prec (default 0.000000001); the default for pow is .80.

Permutational MANOVA

In the following, R1 contains one-way MANOVA data in standard format (without row or column headings), and grps = the number of groups.

PseudoF(R1, grps) = the pseudo F statistic

PERMANOVA(R1, grps, iter) = the p-value for Permutational MANOVA based on iter many iterations (default 1,000).

Box’s Test

In the following, R1 is a range in standard formal if k = 0 (or is omitted) and R1 is in covariance matrix format if k > 0 where k = the number of dependent variables (i.e. all the covariance matrices are k × k).

BOXTEST(R1, k, b) = p-value
BOXM(R1, k) = M statistic
BOXF(R1, k) = F statistic
BOXdf1(R1, k) = df₁
BOXdf2(R1, k) = df₂
BOX(R1, k) – see below

There is also the following array function

BOX(R1, k, lab): outputs a 5 × 1 range consisting of the entries BOXM(R1, k), BOXdf1(R1, k), BOXdf2(R1, k), BOXF(R1, k) and BOXTEST(R1, k) in that order. If lab = TRUE (default FALSE), then a column of labels is appended to the output.

When b is TRUE (or omitted) BOXTEST(R1, k, b) performs the F test, while when it is FALSE the chi-square test is performed instead.

Where R1 and k are as described above, the following is an array function

COVPooled(R1, k) = k × k pooled covariance matrix

Extraction

The following are array functions where range R1 contains data in standard form (with or without column headings) and s is a string that presumably specifies a group (in the first column of R1):

ExtractRows(R1, s, b) = the array that contains all the elements in range R1 for the group labeled s. If b is set to TRUE, the first row of R1 is included in the output (presumably column headings) even if there is no match. If b is omitted it defaults to TRUE.

ExtractCov(R1, s) = the covariance matrix for the group labeled s, based on the data in range R1. If s is the empty string “” or is omitted then the pooled covariance matrix of R1 is returned.

For the following array function range R1 contains data in standard form (with column headings) and s is a string that specifies a column heading (i.e. an entry in the first row of R1) that selects one of the dependent variables.

ExtractCol(R1, s) = an array that consists of all the data in R1 for the dependent variable identified by s, now organized by columns with one column for each group. You should highlight a range with g columns where g = the number of groups.

Factor Analysis

For a k × k range R1, the following is an array function:

VARIMAX(R1): Produces a k × k array containing the loading factor matrix after applying a Varimax rotation to the loading factor matrix contained in range R1.

For a k × k range R1 representing a correlation matrix, the following is an array function:

ExtractCommunalities(R1, iter, prec, eigen) = the 1 × k row vector with the communalities after convergence based on a precision value of prec (default = .00001), maximum number of iter iteration (default = 25) and maximum number of iterations eigen in finding eigenvalues/vectors (default = 100).

Cluster Analysis

The following are array functions:

CLUST(R1, k, R2) = m × 1 column range of cluster numbers 1, 2, …, k calculated by the k-means algorithm corresponding to the m n-tuple data elements in the m × n range R1 where R2 is an m × 1 column range containing the initial cluster number assignments. If R2 is omitted then the k-means++ algorithm is used to calculate the initial cluster number assignments.

ClustAnal(R1, k, iter) = the m × 1 column range of cluster numbers produced by CLUST(R1, k) with the lowest SS_E after iter iterations (here each of the iterations uses the k-means++ algorithm). If iter = 0, however, then the m × 1 column range of initial cluster numbers based on the k-means++ algorithm is returned.

CLUSTERS(R1, R2) = m × 1 column range of cluster numbers 1, 2, …, k corresponding to the centroid described in the n × k range R2 which is closest to the respective data element in the m × n range R1.

CENTROIDS(R1, R2) = n × k range defining the k centroids corresponding to the  m× 1 column range R2 of clusters for the respective data elements in the m × n range R1. Here k = the maximum value in R2.

CENTROIDS(R1, k) = n × k range defining the k initial centroids corresponding to the  data elements in the m × n range R1 based on the k-means++ algorithm

Non-array functions

The following are non-array functions:

CLUSTErr(R1, R2, k) = SS_E value for the data in the m × n range R1 based on the cluster assignment in the m × 1 column range which contains cluster number assignments based on k clusters 1, 2, …, k.

CENTROIDErr(R1, R2, k) = SS_E value for the data in the m × n range R1 based on the centroids described in the n × k range R2.

CLUST_Converge(R1, R2) = TRUE if the m × 1 column array  R2 of cluster numbers 1, 2, …, k calculated by the k-means algorithm for the m n-tuple data elements in the m × n range R1 has converged (i.e. one additional iteration of the algorithm does not result in any changes to the cluster assignments.

These functions also support non-Euclidean distances as well as weighted distances. This is implemented by adding the p (default 2) parameter and Rw weights range to the above functions (except for CENTROIDS). Click here for more details.

Jenks Natural Breaks

JENKS(R1,, lab, iter) – performs Jenks Natural Breaks optimization on the data in range R1 for k classes and outputs a k+1 × 3 range, whose first k rows contain the left and right endpoints of each of the classes followed by the number of elements in the class. The last row contains the total squared deviation of the k classes, followed by the squared deviation for the data in R1 and GVF.

Here k+1 = the number of rows in the highlighted output range. If iter is included then the algorithm uses iter iterations, otherwise, all possible partitions are tested.

JENKS(R1, k, lab, iter) – outputs a k × 1 column range, whose first k−1 rows contain the right endpoint of the first k−1 classes (the endpoint of the k^th category is the largest element in R1) and whose last element is the GVF.

GVF(R1, R2) = GVF for the right endpoints (breaks) in range R2 based on the data in column range R1.

If lab = TRUE (default = FALSE) then an extra column of labels is included.

Discriminant Analysis

LDACoeff(Rt, head): returns an array with the LDA coefficients for the (training) data in Rt consisting of one row for each independent variable whose columns consist of the name of the independent variable, the intercept coefficient, and a coefficient for each dependent variable. If head = TRUE (default), then the data in Rt contains column headings (corresponding to the names of the dependent variables); these headings are also appended to the output from this function.

LDAPredC(R0, Rc, Rp, lab): returns an array whose rows contain the probabilities for each vector in the data array R0 (which contains no row/column headings) using the LDA coefficient array in array Rc (without column headings) and prior probabilities in the column array Rp. A column of names of the independent variable that has the highest probability is also appended to the output.

LDAPred(R0, Rt, Rp, lab) = LDAPredC(R0, LDACoeff(Rt,FALSE), Rp, lab), i.e. the predictions for the vectors in R0 based on the LDA model defined by Rt and Rp.

QDAPred(R0, Rt, Rp, lab): returns an array whose rows contain the probabilities for each vector in the data array R0 (which contains no row/column headings) and the name of the independent variable with the highest probability as for the LDAPred function, except that the QDA model is used instead of the LDA model.

DAClassification(Rt, Rp, linear): returns a classification for the training data in Rt and priors in Rp. If linear = TRUE (default), then the classification table is based on an LDA model, while if linear = FALSE then a QDA model is used instead.

DASummary(R1): returns a summary of the classification table in range R1

If Rp is omitted then equally probably priors are used. If lab = TRUE (default FALSE) then column headings are added to the output.

Correspondence Analysis

In the following R1 is an array containing an m x n contingency table, R2 is an mm × n array containing supplementary row profiles and R3 is an m × nn array containing supplementary column profiles. None of these arrays include row or column headings or totals. We also assume that k = min(m, n) – 1.

CARowFactors(R1): returns an m × k array in which the i^th row contains the i^th row factor vector for R1

CAColFactors(R1): returns an n × k array in which the i^th row contains the i^th column factor vector for R1

CARowFactors(R1, R2): returns an mm × n array in which the i^th row contains the row factor vector for the i^th supplementary row profile in R2 based on the data in R1.

CAColFactors(R1, R3): returns an m × nn array in which the i^th row contains the factor vector for the i^th supplementary column profile in R3 based on the data in R1.

CAEigen(R1): returns an k × 1 column array with the CA eigenvalues for R1

Gaussian Mixture Models (GMM)

UniGMM(R1, k, lab, iter, prec, R0, kiter, nruns): outputs a GMM model which for each cluster 1, 2, …, k, estimates the univariate normal distribution for that cluster (consisting of the mean vector and standard deviation) as well as the probability that a data element in R1 best fits that particular cluster.

UniGMMClust(R1, R2, sorted): outputs an array whose first column contains the data in R1, and whose next k columns show the probability that each data element in R1 belongs to each one of the k clusters based on the GMM model described in R2. Finally, the last column of the output shows the cluster with the highest probability for each data vector in R1. If sorted = TRUE then the output is sorted by the cluster number in the last column; otherwise (default) the order is the same as in R1. 

MultiGMM(R1, k, lab, iter, prec, R0, kiter, nruns): outputs a GMM model which for each cluster 1, 2, …, k, estimates the multivariate normal distribution for that cluster (consisting of the mean vector and covariance matrix) as well as the probability that a data vector in R1 best fits that particular cluster.

MultiGMMClust(R1, R2, sorted): outputs an array whose first columns consist of the data in R1, and whose next k columns show the probability that each data vector in R1 belongs to each one of the k clusters based on the GMM model described in R2. Finally, the last column of the output shows the cluster with the highest probability for each data vector in R1. If sorted = TRUE then the output is sorted by the cluster number in the last column; otherwise (default) the order is the same as in R1.

References

Penn State University (2013) STAT 505: Applied multivariate statistical analysis (course notes)
https://online.stat.psu.edu/stat505/lesson/7

Rencher, A.C. (2002) Methods of multivariate analysis (2nd Ed). Wiley-Interscience, New York.
https://www.ipen.br/biblioteca/slr/cel/0241

Johnson, R. A. and Wichern, D. W. (2007) Applied multivariate statistical analysis. 6th Ed. Pearson.
https://www.webpages.uidaho.edu/~stevel/519/Applied%20Multivariate%20Statistical%20Analysis%20by%20Johnson%20and%20Wichern.pdf

Pituch, K. A. and Stevens, J. P. (2016) Applied multivariate statistical analysis for the social sciences. Routledge.