What if there are many outcomes or many variables?

Start here

Multivariate or multivariable: which word describes your analysis?

Strictly, a multivariate analysis has two or more outcome (dependent) variables. A model with one outcome and several predictors is multivariable. Many papers use the first word when they mean the second.

In exam questions and textbooks, “multivariate analysis” is also used as a heading for all methods that handle many variables at once, including those with no outcome variable at all, such as cluster analysis and factor analysis. This chapter covers that whole group. When you write your own methods section, use “multivariable” for one outcome with several predictors, and name the actual method.

Two types of question

Which questions do multivariate methods answer?

Start by asking whether your variables fall into two sets, explanatory and outcome, or form one set in which you look for structure.

Sample size

How many patients do you need for the number of variables?

Every multivariate method estimates how each variable relates to every other variable. The number of those relationships grows much faster than the number of variables.

• The count grows quickly. With p variables, the covariance matrix contains p(p + 1) ÷ 2 separate values to estimate: 36 for 8 variables and 1,830 for 60 variables.
• Fewer patients than variables makes some methods impossible. With n patients, a sample covariance matrix has rank of at most n − 1. In the PLS teaching example below, 60 protein levels in 30 patients give a covariance matrix of rank 29, so it cannot be inverted. MANOVA, discriminant analysis and ordinary least squares regression all need that inverse. PLS regression and PCA do not.
• MANOVA and discriminant analysis: the within-group matrix W can be inverted only when the number of patients minus the number of groups is at least the number of variables. Having only just enough patients gives very unstable results.
• Factor analysis: MacCallum and colleagues showed that the idea of one minimum sample size, or one minimum ratio of patients to variables, is a misconception. The number needed depends on the communalities of the variables and on how many variables measure each factor. Nieser and colleagues, whose paper appears later in this chapter, quote a simulation result that for “moderate loadings (0.50), 2 or 3 factors, and 8 or more indicators, a sample size of 150–200 yields minimal bias”.
• Logistic and Cox models are limited by the number of events per predictor. See the correlation and regression chapter.

Choose the variables from the clinical question before the study. Selecting them because they were significant one at a time ignores the correlation between them, which is the reason for using a multivariate method.

Several outcomes

Several correlated outcomes in three groups: why not run one ANOVA per outcome?

Two problems. Each separate test has its own 5% false-positive rate, so the chance of at least one false positive rises. And separate tests ignore the correlation between the outcomes, which can hide a real difference.

What MANOVA does

Multivariate analysis of variance (MANOVA) compares the groups on all the outcomes at once. ANOVA splits the total sum of squares into a part between groups and a part within groups. MANOVA does the same with matrices that also hold the sums of cross-products, so the correlation between outcomes is included. With two groups, MANOVA is Hotelling's T², the multivariate form of the t-test.

Wilks' λ = |W| ÷ |T|

W

within-group matrix: sums of squares and cross-products of each patient around the group mean

T

total matrix: the same, around the overall mean

|W|, |T|

determinants, a single number summarising the spread in all the outcomes together

λ

the share of the total spread that is not explained by the groups. 1 = no group difference, towards 0 = large difference

Worked example: rise in heart rate and MAP with three drugs

1. Separate one-way ANOVAs: heart rate F(2, 27) = 1.57, p = 0.226. MAP F(2, 27) = 2.39, p = 0.111.
2. Within-group matrix W = [1685.0 1115.3 ; 1115.3 1338.1], determinant 1,010,804.4.
3. Total matrix T = [1881.2 1086.8 ; 1086.8 1574.7], determinant 1,781,191.4.
4. Wilks' λ = 1,010,804.4 ÷ 1,781,191.4 = 0.5675.
5. With two outcomes the conversion to F is exact: F = [(1 − √λ) ÷ √λ] × [(N − g − 1) ÷ (g − 1)] = (0.2467 ÷ 0.7533) × (26 ÷ 2) = 4.26, with 4 and 52 degrees of freedom.
6. p = 0.0047. Pillai's trace gives the same conclusion: 0.458, F(4, 54) = 4.01, p = 0.0064.

The drugs differ when both outcomes are considered together, although neither ANOVA reached 5%. The chance of at least one false positive in two separate tests at 5% would be 1 − 0.95² = 9.75% if the tests were independent.

The four test statistics and the assumptions

• IBM SPSS reports Pillai's trace, Wilks' lambda, Hotelling's trace and Roy's largest root, each with an approximate F.
• R's summary.manova() uses Pillai's trace by default. Its documentation notes that Wilks' statistic is the most popular in the literature and that the Pillai–Bartlett statistic is recommended by Hand and Taylor (1987).
• In a simulation study by Ateş and colleagues (2019), no statistic was best in every situation, and all four were affected when the covariance matrices differed or the data were not normal, particularly with unequal group sizes.
• Assumptions listed by IBM SPSS: the data are a random sample of vectors from a multivariate normal population, and the variance-covariance matrices are the same in all cells. Box's M tests the equality of the covariance matrices.
• After a significant MANOVA, SPSS runs post hoc tests for each outcome separately. Say in advance which follow-up analyses you will do.

Interactive

Change the correlation between two outcomes and watch the MANOVA p-value

The separate t-tests do not use the correlation between the outcomes. MANOVA does, so the same group differences can give very different MANOVA results.

Two sets of measurements

How strongly is one set of measurements related to another set?

Pearson r relates one measurement to one other. Canonical correlation relates a whole set to a whole set, by first combining each set into a weighted sum.

• First canonical correlation: 0.58. The second and third are 0.20 and 0.03.
• Test that all canonical correlations are zero: Wilks' λ = 0.636, F = 2.99 with 9 numerator degrees of freedom, p = 0.003.
• Reading it: the combination of body size measurements is moderately related to the combination of sleep study results. The weights show which variables contribute most to each combination, but they can be unstable when the variables within a set are strongly correlated with each other.
• Use: mainly exploratory, when both sets have several variables and neither is a single outcome. If one set is really a single outcome, the method becomes multiple regression and the canonical correlation becomes the multiple correlation R.

Known groups

Can a set of measurements sort patients into groups that are already known?

Discriminant analysis starts from patients whose group is known. It finds the weighted sum of the measurements that separates the groups best, then uses that score to classify new patients.

The discriminant function

D = b₀ + b₁x₁ + b₂x₂ + … + bₚxₚ

The weights make the ratio of between-group to within-group variation in D as large as possible. With g groups there are at most g − 1 functions. With two groups there is one.

D = −1.455 + 0.1613 × MCV − 1.9456 × red cell count

1. Group means: IDA MCV 72.96 fL, red cell count 4.45 × 10¹²/L. BTT MCV 63.68 fL, red cell count 5.38 × 10¹²/L.
2. Mean score (centroid): IDA +1.646, BTT −1.646. With equal group sizes and equal prior probabilities the cut-off is halfway, at 0.
3. A patient with MCV 66 fL and red cell count 5.0 × 10¹²/L: D = −1.455 + 10.646 − 9.728 = −0.54. Below 0, so classified as thalassaemia trait.
4. Standardised coefficients: MCV 0.75, red cell count −0.81. Both contribute, in opposite directions. Wilks' λ = 0.262.

• Classification table: actual group against predicted group. IBM SPSS calls it the summary table and notes it is sometimes called the confusion matrix.
• Leave-one-out classification: each patient is classified by a function built from all the other patients. Here it also gave 96%. Accuracy in the same patients used to build the function is not a test on new patients.
• Prior probabilities: SPSS offers “All groups equal” or “Compute from group sizes”. If thalassaemia trait is rarer than iron deficiency in your clinic, equal priors will over-call it.
• Assumptions listed by IBM SPSS: cases are independent, predictors have a multivariate normal distribution, and the within-group variance-covariance matrices are equal across groups. Box's M tests the last one and is sensitive to departures from multivariate normality.

Comparison

Discriminant analysis or logistic regression?

Both can predict group membership from several measurements. They differ in their assumptions and in what they report.

In some clinical comparisons the two methods gave similar accuracy. Ebrahimzadeh and colleagues found similar accuracy for predicting unwanted pregnancy and recommended logistic regression when the interpretability of the results matters. Choose before you see the results, from the type of predictors and the question.

Unknown groups

Can the data find patient groups that nobody defined in advance?

Cluster analysis puts patients who are close to each other in the measured variables into the same group. There is no outcome variable and no known group label.

Step 1: measure how far apart two patients are

distance = √[(z₁ᵢ − z₁ⱼ)² + (z₂ᵢ − z₂ⱼ)² + …]

This is the Euclidean distance between patients i and j. Standardise each variable first (z = (x − mean) ÷ SD). IBM SPSS advises standardising when variables have large differences in scale, because BMI in the 20s and HbA1c in single figures would otherwise count unequally.

Step 2a: hierarchical clustering and the dendrogram

Hierarchical agglomerative clustering starts with every patient as a separate cluster and joins the two closest clusters, step by step, until only one is left. The linkage method decides the distance between two clusters: nearest neighbour (single linkage), furthest neighbour (complete linkage), between-groups average, centroid, or Ward's method, which joins the pair of clusters that increases the within-cluster sum of squares least.

Step 2b: k-means

1. Choose k, the number of clusters.
2. Pick k starting centres, for example k patients at random.
3. Assign every patient to the nearest centre.
4. Move each centre to the mean of its patients.
5. Repeat steps 3 and 4 until no patient changes cluster.

The result can depend on the starting centres. The documentation of R's kmeans() says that trying several random starts (nstart > 1) is often recommended.

Step 3: how many clusters?

• Silhouette width of a patient: close to 1, well inside its cluster. Around 0, between two clusters. Negative, probably in the wrong cluster.
• Cluster profiles: describe each cluster by the means of the variables. The names (“high BMI”) come from the researcher, not from the software.
• Differences between clusters are built in. IBM SPSS notes that the F tests in its k-means ANOVA table are only descriptive and their probabilities should not be interpreted, because the procedure forms groups that differ.
• k-means assumes k is right. IBM SPSS warns that if you choose an inappropriate number of clusters or omit important variables, the results may be misleading.
• Check before you believe the clusters: whether the data have any tendency to form clusters, and whether the clusters are stable when the data are resampled. The INSPIRED study in the “In a real paper” section did both.

Interactive

Step through k-means and watch the clusters and centres change

Each press shows one step: patients join the nearest centre, then the centres move to the mean of their patients.

Comparison

Hierarchical clustering or k-means?

Many correlated variables

Can many correlated measurements be replaced by a few?

If two measurements rise and fall together, one well-chosen weighted sum can keep most of the information in both. Principal component analysis (PCA) finds these sums.

The idea with two variables

Waist and hip circumference in 40 made-up adults: correlation 0.78. Draw any line through the middle of the scatter and project each patient onto it. The spread of the projected positions is the variance kept by that line. The first principal component (PC1) is the line that keeps the most variance. PC2 is at right angles to it and keeps what is left. Jolliffe and Cadima describe PCA as creating new uncorrelated variables that successively maximise variance, found by solving an eigenvalue problem.

The formula and the worked example

Covariance matrix (cm²): variance of waist 105.2, variance of hip 40.9, covariance 51.1. Total variance 146.1.

eigenvalues λ = (s₁₁ + s₂₂) ÷ 2 ± √[((s₁₁ − s₂₂) ÷ 2)² + s₁₂²]

1. (105.2 + 40.9) ÷ 2 = 73.0, and √[32.2² + 51.1²] = 60.4.
2. λ₁ = 133.5 cm², λ₂ = 12.6 cm². The eigenvalues add up to the total variance.
3. Variance explained: PC1 133.5 ÷ 146.1 = 91.4%, PC2 8.6%.
4. Loadings (eigenvector weights) of PC1: 0.875 for waist, 0.483 for hip. PC1 score = 0.875 × (waist − 92.6) + 0.483 × (hip − 101.05).

With more variables: the scree plot

With p variables there are p components. Plot the eigenvalues in order. Components before the point where the curve flattens carry most of the variance. The scree plot below comes from the 8-item questionnaire used in the factor analysis example: the first two components explain 41.8% and 18.9% of the variance.

• Covariance or correlation matrix: PCA on the covariance matrix lets variables with large variances dominate. IBM SPSS describes the correlation matrix as useful when variables are measured on different scales. In R, prcomp() does not scale by default (scale. = FALSE).
• Uses: reducing many correlated variables before another analysis, building a single index, showing patients on a two-dimensional plot, dietary patterns.
• Limits: components are weighted sums chosen for variance, not for clinical meaning, and PCA only captures straight-line relationships.

Interactive

Turn a line through the waist and hip data and watch the variance it keeps

PCA searches all directions for the one along which patients are most spread out. Here you do the search by hand.

Hidden factors

Are the items of a questionnaire driven by a few hidden factors?

Fatigue and low mood cannot be measured directly. A questionnaire measures them through items. Factor analysis asks how many hidden (latent) factors are needed to explain the correlations between the items, and which items belong to each factor.

Exploratory or confirmatory

Tavakol and Wetzel describe CFA as a theory-driven approach that tests how well the data fit a proposed model, in which the researcher must identify the factor model before analysing the data.

The steps of an exploratory factor analysis, with the worked example

Teaching example: 8 items scored 0 to 10 by 240 made-up patients after discharge. Items 1–4 describe tiredness, items 5–8 mood.

1. Are the items correlated enough? The Kaiser-Meyer-Olkin (KMO) measure of sampling adequacy compares the correlations with the partial correlations: KMO = 0.82. The psych package documentation quotes Kaiser's labels: above 0.9 marvelous, 0.80s meritorious, 0.70s middling, 0.60s mediocre, 0.50s miserable, below 0.5 unacceptable. Bartlett's test of sphericity tests whether the correlation matrix is an identity matrix (all correlations zero): χ² = 561.3, df 28, p < 0.001.
2. How many factors? Two eigenvalues are above 1 (3.34 and 1.51), and parallel analysis also gives 2 (see the scree plot above). jamovi uses parallel analysis as its default rule.
3. Extract the factors. IBM SPSS offers principal components, unweighted least squares, generalized least squares, maximum likelihood, principal axis factoring, alpha factoring and image factoring. This example uses minimum residual, the default in jamovi and in R's psych package.
4. Rotate. Turning the factor axes does not change how well the factors explain the correlations, but it can make each item load clearly on one factor. Varimax is an orthogonal rotation (factors stay uncorrelated). Direct oblimin and promax are oblique rotations (factors may correlate). For oblique rotations SPSS displays the pattern, structure and factor correlation matrices.
5. Name the factors from the items with large loadings, here factor 1 “physical fatigue” and factor 2 “low mood”.

• Loading: Tavakol and Wetzel describe it as the correlation between the item and the factor. That reading holds for uncorrelated factors such as these.
• Communality: the sum of an item's squared loadings, the share of its variance explained by the factors. Item 1: 0.79² + 0.12² = 0.64.
• Variance explained: factor 1 2.06 (25.7%), factor 2 1.79 (22.3%), together 48.0% of the variance of the 8 items. Before rotation the split was 35.4% and 12.6%, with the same total.
• Oblique solution: with oblimin the two factors correlate 0.43. Tavakol and Wetzel note that most constructs in social science research are assumed to be inter-related and should use an oblique rotation.

Interactive

Rotate the factor axes and watch the loadings table simplify

Rotation does not change the fit. It changes how the same solution is described, so that each item loads mainly on one factor.

Comparison

Principal component analysis or factor analysis: what is the difference?

The two are often run from the same menu and give similar-looking tables, but they answer different questions.

Tavakol and Wetzel note that PCA is often applied and described, in error, as exploratory factor analysis. If your aim is to measure a latent construct, use factor analysis with maximum likelihood or principal axis factoring.

Paths and hidden variables

Can you test a diagram of causes that includes hidden variables?

Ordinary regression gives each predictor's effect adjusted for the others. Path analysis goes further: you draw which variable affects which, and the analysis estimates every arrow, including effects that pass through another variable.

Worked example: a simple path analysis

1. Regress adherence on depression: a = −0.46 (p < 0.001).
2. Regress HbA1c on adherence and depression together: b = −0.36 (p < 0.001), direct effect c′ = 0.13 (p = 0.084).
3. Regress HbA1c on depression alone: total effect c = 0.29.
4. Indirect effect through adherence = a × b = (−0.463) × (−0.363) = 0.168. Check: c = c′ + a × b = 0.125 + 0.168 = 0.293. About 57% of the total effect passes through adherence.
5. Bootstrap 95% confidence interval for the indirect effect (2,000 resamples): 0.10 to 0.24. It excludes 0.

From path analysis to structural equation modelling (SEM)

• Beran and Violato describe SEM as a set of techniques that analyse the relationships of observed and latent variables while accounting for measurement error.
• SEM has two parts. The measurement model links each latent variable to the measured items that reflect it, as in confirmatory factor analysis. The structural model gives the relationships among the latent variables.
• In path analysis, as in the example above, every variable is measured directly.
• Fit is judged with indices such as the comparative fit index (CFI), the Tucker-Lewis index (TLI) and the root mean square error of approximation (RMSEA). jamovi's confirmatory factor analysis reports these three by default. Beran and Violato cite CFI of 0.90 or more and RMSEA below 0.05 as adequate. Whichever cut-offs you use, state them and cite their source.
• A model that fits well is consistent with the data, not proof that the arrows point the right way. Other diagrams can fit equally well. Draw the model from clinical knowledge and the study design before analysing.

More variables than patients

What does partial least squares regression add when predictors outnumber patients?

The correlation and regression chapter introduces PLS regression. Here is how it works, how the number of components is chosen, and why the error in the fitting data cannot be trusted.

PLS regression builds a few components from the predictors, like PCA, but chooses them to predict the outcome. IBM SPSS describes PLS as particularly useful when predictor variables are highly correlated or when the number of predictors exceeds the number of cases, and says it first extracts latent factors that explain as much of the covariance as possible between the independent and dependent variables. See the basic introduction and the PLS diagram in the correlation and regression chapter.

The algorithm, one component at a time (PLS1, one outcome)

1. Standardise every predictor and centre the outcome.
2. Give each predictor a weight proportional to its covariance with the outcome. The first component is the weighted sum of the predictors.
3. Regress the outcome on this component.
4. Remove from the predictors and from the outcome the part explained by the component.
5. Repeat on what is left to get the second component, which is uncorrelated with the first, and so on.
6. Choose the number of components by cross-validation, then predict the outcome from those components.

Worked example: eGFR from 60 serum proteins in 30 patients

Made-up teaching data. The 60 protein levels are strongly correlated. Their covariance matrix has rank 29, so ordinary regression cannot be fitted. PLS models with 1 to 10 components were checked by leave-one-out cross-validation: each patient is predicted by a model built from the other 29.

• Chosen model: 3 components. Cross-validated RMSE 3.25 mL/min/1.73 m² (SD of eGFR 13.25), cross-validated Q² = 0.94. R² in the fitting data 0.97.
• Why not PCA first? With one component, principal component regression had a cross-validated RMSE of 7.91, against 3.56 for PLS, because the first principal component was not chosen to predict eGFR.
• Variable importance in projection (VIP): IBM SPSS produces VIP scores by default, a summary of how much each predictor contributes to the components.
• PLS-DA: for a yes/no outcome, PLS-discriminant analysis uses PLS regression with a binary dummy outcome. It is commonly used in metabolomics, and Szymańska and colleagues list double cross-validation and permutation testing among the approaches used to validate it.
• Limits: components mix many predictors, so the model is better for prediction than for explaining the effect of one variable. Cross-validation inside one data set is not external validation.

Decide

Which multivariate method fits your question?

Answer the questions in order. Each answer links to its section.

In a real paper

Cluster analysis and factor analysis in published studies

How to run it

Running multivariate methods in SPSS, R and jamovi

Menu names are from the IBM SPSS Statistics documentation and the jamovi documentation. Function names and defaults are from the R and CRAN reference manuals.

Defaults that change your results

• prcomp() does not standardise unless you set scale. = TRUE.
• kmeans() uses one random start by default (nstart = 1) and the Hartigan-Wong algorithm, with iter.max = 10.
• hclust() uses complete linkage by default. Its documentation states that "ward.D2", not "ward.D", implements Ward's (1963) criterion.
• MASS::lda() uses the class proportions in your data as prior probabilities unless you give prior.
• psych::fa() defaults to minimum residual extraction and oblimin rotation. factanal() in base R uses maximum likelihood and varimax.
• psych::cortest.bartlett() assumes n = 100 if you do not give the sample size.
• jamovi Exploratory Factor Analysis defaults to parallel analysis, minimum residual extraction and oblimin rotation. Its Principal Component Analysis defaults to varimax.

R code for the teaching examples

# MANOVA
fit <- manova(cbind(hr_rise, map_rise) ~ drug, data = haemo)
summary(fit)                    # Pillai's trace
summary(fit, test = "Wilks")    # Wilks' lambda

# Discriminant analysis with equal priors and leave-one-out classification
library(MASS)
ld  <- lda(group ~ mcv + rbc, data = anaemia, prior = c(0.5, 0.5))
loo <- lda(group ~ mcv + rbc, data = anaemia, prior = c(0.5, 0.5), CV = TRUE)
table(anaemia$group, loo$class)

# Hierarchical clustering and k-means on standardised variables
z  <- scale(dm[, c("bmi", "hba1c")])
hc <- hclust(dist(z), method = "ward.D2")
plot(hc)
cl <- cutree(hc, k = 3)
km <- kmeans(z, centers = 3, nstart = 100)
library(cluster)
mean(silhouette(km$cluster, dist(z))[, "sil_width"])

# PCA
pc <- prcomp(body[, c("waist", "hip")])
pc$sdev^2       # eigenvalues
pc$rotation     # loadings

# Exploratory factor analysis
library(psych)
KMO(items)
cortest.bartlett(cor(items), n = nrow(items))
fa.parallel(items)
efa <- fa(items, nfactors = 2, fm = "minres", rotate = "varimax")
efa$loadings

# Confirmatory factor analysis and a path model with an indirect effect
library(lavaan)
cfa_fit <- cfa(' fatigue =~ i1 + i2 + i3 + i4
                 mood    =~ i5 + i6 + i7 + i8 ', data = items)
fitMeasures(cfa_fit, c("cfi", "rmsea"))
path_fit <- sem(' adherence ~ a*depression
                  hba1c ~ b*adherence + depression
                  indirect := a*b ', data = dm, se = "bootstrap")

# PLS regression with leave-one-out cross-validation
library(pls)
pl <- plsr(egfr ~ ., data = proteins, ncomp = 10, scale = TRUE, validation = "LOO")
RMSEP(pl)

How to report it

What to write in your thesis or paper

Name the method, the software and version, every choice that changes the result, and the numbers a reader needs to check your conclusion.

Worked sentences using this chapter's teaching data

• MANOVA: “The rise in heart rate and in MAP were compared across the three groups with a one-way MANOVA. The groups differed on the two outcomes considered together (Wilks' λ = 0.57, F(4, 52) = 4.26, p = 0.005; Pillai's trace = 0.46, p = 0.006). One-way ANOVAs for each outcome gave p = 0.23 for heart rate and p = 0.11 for MAP.”
• Discriminant analysis: “Linear discriminant analysis with MCV and red cell count, using equal prior probabilities, gave D = −1.455 + 0.161 × MCV − 1.946 × red cell count. It correctly classified 48 of 50 patients (96%), and 96% with leave-one-out classification.”
• Cluster analysis: “BMI and HbA1c were standardised to z-scores. k-means clustering (100 random starts) was run for k = 2 to 6, and k = 3 was chosen because it gave the largest average silhouette width (0.67). Ward's hierarchical clustering gave the same three clusters.”
• Exploratory factor analysis: “Sampling adequacy was acceptable (KMO = 0.82; Bartlett's test χ²(28) = 561.3, p < 0.001). Two factors were retained by parallel analysis and eigenvalues greater than 1 (3.34 and 1.51). Factors were extracted by minimum residual and rotated by varimax. Loadings of 0.40 or more defined factor membership. The two factors explained 48.0% of the variance.”
• PLS regression: “PLS regression was fitted to 60 standardised protein levels. The number of components was chosen by leave-one-out cross-validation (3 components, RMSE 3.25 mL/min/1.73 m², Q² = 0.94). The model has not been validated in an external sample.”

Always report

• Why a multivariate method was needed, and whether it was planned in the protocol.
• The variables included and how they were chosen, standardisation, and how missing values were handled.
• For MANOVA and discriminant analysis: the test statistic named, and the check of equal covariance matrices.
• For cluster analysis: Clatworthy and colleagues found that published cluster analyses most often omitted the similarity measure and the computer program, and that the procedures used to choose the number of clusters and to validate them were often reported in too little detail to repeat.
• For factor analysis: KMO and Bartlett's test, the rule for the number of factors, extraction and rotation methods, the loading cut-off, the full loading table and the variance explained.
• For any classification or prediction model: accuracy by cross-validation or in a separate sample, not only in the fitting data.

By specialty

Multivariate problems by specialty

Common myths

Common mistaken beliefs

• ✗Myth: "Multivariate logistic regression" describes a model with one outcome and several predictors.

  ✓That model is multivariable. Multivariate means two or more outcome variables. In one review, 25 of 30 articles that said "multivariate" had used multivariable models (Hidalgo and Goodman, 2013).

• ✗Myth: One ANOVA per outcome is as good as a MANOVA.

  ✓Separate tests each carry their own false-positive rate and ignore the correlation between outcomes. In this chapter's example, two ANOVAs gave p = 0.23 and 0.11, and MANOVA gave p = 0.005.

• ✗Myth: PCA is a type of factor analysis.

  ✓PCA builds weighted sums that keep total variance. Factor analysis models hidden factors and separates common from unique variance. Tavakol and Wetzel note that PCA is often described, in error, as exploratory factor analysis.

• ✗Myth: Eigenvalues greater than 1 fix the number of factors.

  ✓It is one guide among several. Nieser and colleagues found three eigenvalues above 1 but kept two factors, because only one item loaded on the third.

• ✗Myth: Cluster analysis discovers real disease subtypes.

  ✓k-means always returns the number of clusters you ask for. Clusters need checks of cluster tendency, stability and outcomes, and replication in another sample.

• ✗Myth: Ten patients per variable always makes a factor analysis stable.

  ✓MacCallum and colleagues showed that no single ratio works. The sample needed depends on the communalities and on how many items measure each factor.

Why this matters to you

Five consequences of choosing the wrong multivariate method

1. Protocol sent back

   "Multivariate analysis will be done" does not name a method or an outcome. ICH E9 expects one primary variable and the handling of multiplicity to be stated in the protocol.

2. A wrong conclusion in the thesis

   Separate tests on correlated outcomes can report "no difference" when the outcomes together differ, or a false difference when many outcomes are tested.

3. Questions you cannot answer at the viva

   "Why PCA and not factor analysis?", "How did you choose k?" and "Which rotation, and why?" need answers written into your methods.

4. Misreading papers

   Clusters named by the authors can look like new diseases. Factor names are the authors' interpretation of the loadings. Check how both were chosen and validated.

5. Harm to patients

   A classification rule with 96% accuracy in the patients used to build it may do worse in new patients. Using it without validation can send patients down the wrong diagnostic path.

History

History of multivariate methods

1. 1901

   Karl Pearson

   Published "On lines and planes of closest fit to systems of points in space". Jolliffe and Cadima date the earliest literature on principal component analysis to Pearson and Hotelling.

2. 1904

   Charles Spearman

   Published "General intelligence, objectively determined and measured", the start of factor analysis.

3. 1921

   Sewall Wright

   Developed the method of path coefficients (path analysis) between 1914 and 1921.

4. 1933

   Harold Hotelling

   Published "Analysis of a complex of statistical variables into principal components" in two parts in the Journal of Educational Psychology.

5. 1936

   Ronald A. Fisher

   Published "The use of multiple measurements in taxonomic problems", the paper behind Fisher's linear discriminant.

6. 1936

   Harold Hotelling

   Published "Relations between two sets of variates", the reference R gives for canonical correlation.

7. 1958

   Henry Kaiser

   Defined the varimax criterion for analytic rotation in factor analysis.

8. 1960s

   Herman Wold

   First proposed partial least squares. Interest in chemistry expanded in the 1980s.

9. 1966

   Raymond Cattell

   Published "The scree test for the number of factors".

10. 1969

    Karl Jöreskog

    Published "A general approach to confirmatory maximum likelihood factor analysis".

11. 1987

    Peter Rousseeuw

    Published silhouettes as a graphical aid to the interpretation and validation of cluster analysis.

Guidelines

Several outcomes in guidelines

Reading a paper

Six checks on a multivariate analysis in any paper

• Does "multivariate" really mean two or more outcomes, or is it a multivariable model with one outcome?
• Were correlated outcomes analysed together, or with a correction for multiple testing, and is the test statistic named?
• For cluster analysis: were the variables standardised, the distance and algorithm named, the choice of k explained and the clusters checked for stability?
• For PCA or factor analysis: are KMO, Bartlett's test, the rule for the number of factors, the extraction and rotation methods and the full loading table reported?
• For discriminant analysis, PLS or any classification: is accuracy shown by cross-validation or in a separate sample, not only in the patients used to build the model?
• For path analysis or SEM: was the diagram specified before the analysis, are fit indices given, and do indirect effects have confidence intervals?

Why is it called that?

Word surgery

Exam answers

Model answer outlines

Past papers

TNMGRMU questions on multivariate methods

Questions from MD, MS, DM, MCh and Ph.D. papers that this chapter helps you answer. Tap a sitting to open the university's own question paper.

1. ×3Multivariate analysisPh.D. RM · Jul 2022Ph.D. RM · Jul 2014DM Cardiology · Oct 1998
2. ×1Principal component analysisDM Clinical Immunology and Rheumatology · Aug 2023
3. ×1(a) Explain what is meant by ‘Multivariate Analysis’. (b) Explain in detail Partial Least Square Regression. (c) Explain what is Meta Analysis, the steps of Meta Analysis, the methods and assumptions of Meta AnalysisPh.D. RM · Jul 2013