Multivariate Regression: Interpretation and APA Reporting

The Basics

When to Use It

Use it when you have two or more related continuous outcomes and the same set of predictors for all of them, measured on the same people, and you want to know which predictors matter for the outcomes as a set, while keeping every outcome's own coefficients. For example:

Do sleep quality, social support and exercise predict how depressed, anxious and stressed adults feel, taken together?
Do hours of screen time and parental monitoring predict adolescents' scores on both the internalizing and externalizing scales of a behavior checklist?
Do working-memory span and processing speed predict children's reading fluency and math fluency test scores?

If the outcomes answer separate questions, separate multiple regressions say the same thing more simply: the coefficients are identical. If you think one outcome drives another, or the outcomes reflect a latent trait, you want path analysis or another structural equation model instead.

Vocabulary You Need

Definition Multivariate Outcome (

\mathbf{Y}

)

Two or more outcomes, such as depression, anxiety and stress scores, that a single analysis treats as a set rather than one at a time. Each person contributes a score on every outcome, so the data form a table with one column per outcome. Analyzing them together lets a model use how the outcomes vary together, and lets one test ask whether a predictor or group matters for any of them. It only makes sense when the outcomes belong together conceptually; outcomes that answer separate questions are better analyzed separately. Glossary →

Definition Partial Regression Coefficient (

b_j

)

In a regression with several predictors, the change in the predicted outcome for a one-unit rise in one predictor when the other predictors are held at the same values. It's in the outcome's units per unit of the predictor, and it describes only the part of that predictor that doesn't overlap with the others, so it usually differs from the slope the predictor would get on its own. “Holding constant” is statistical, not experimental: it doesn't show what would happen if you changed the predictor in real life. In a model with a product term, a predictor's coefficient is its slope where the other predictor in the product equals zero, and the product's coefficient is how much that slope differs between people one unit apart on the other. Glossary →

Definition Residual Correlation (

r_e

)

The correlation between two outcomes' residuals from the same model: whether people who score higher than predicted on one outcome also tend to score higher than predicted on the other. It's the part of the outcomes' overlap the predictors don't predict, often from a shared cause left out of the model or from measures that tap a common trait. It can differ a lot from the raw correlation between the outcomes, which also includes whatever the predictors explain. Glossary →

Definition Hotelling's T² (

T^2

)

A squared t statistic generalized to several outcomes at once. It adds up how far each outcome's estimate is from zero, in standard errors, after weighting by the inverse of the outcomes' covariance, so outcomes that carry overlapping information aren't counted twice. With a single outcome it is just t². It can be converted to an exact F statistic, and it underlies the multivariate tests of any effect with one degree of freedom, such as two groups or one continuous predictor. Named after Harold Hotelling (1931). Glossary →

Definition Pillai's Trace (

V

)

The sum, across the dimensions on which the groups can differ, of the share of variation each dimension explains, where a dimension is a weighted blend of the outcomes chosen to separate the groups (the best blend, then the next best uncorrelated with it, and so on). It runs from 0 (the group means are identical on every outcome) to s, the smaller of the number of outcomes and the effect's degrees of freedom (with groups, the number of groups minus 1; for one continuous predictor, 1), so it isn't a percentage of all the variation. Of the four usual MANOVA statistics it holds its false-positive rate best when the covariance matrices differ or the data aren't normal, which is why it is a common default. Also called the Pillai-Bartlett trace. Glossary →

Assumptions

1. Independent Observations The people (or other units) in your study don't influence each other.

What It Means: Each person's data comes from a separate, unconnected source, so knowing one person's scores tells you nothing about anyone else's. In designs that measure people more than once, their own scores are meant to be related; it's different people who must be independent. Methods built for clustered data, such as multilevel models and GEE, also allow members of the same group to be related, as long as the analysis names the grouping; then it's the groups that must be independent.
Applies To: Observations
How to Check: By design, not by a test. Did each person contribute their own data, with no clusters such as roommates, couples, therapy groups or classrooms whose members influence each other? If there are clusters, does the analysis model every one of them?
How Much It Matters: A lot. If people came in clusters (couples, coworkers, patients of the same therapist), their residuals are more alike than the model expects. Standard errors, especially for predictors that differ between clusters, usually come out too small, and so do the p-values, univariate and multivariate.

If It's Violated

Alternative Test Multilevel Regression (HLM). A multilevel model with the outcomes stacked, which estimates the same coefficients while accounting for the clustering.
Design Change Sample people independently, one per household or group.

2. Linear Relationship The relationship between your variables follows a straight line, not a curve.

What It Means: As one variable goes up, the other tends to change by about the same amount across its whole range. A curved relationship, one that levels off or a U shape, breaks this.
Applies To: Relationship
How to Check: Look at a scatterplot of each pair of variables the analysis relates: a roughly oval cloud tilted one way is fine; a clear bend, a ceiling or a U is not. In a regression with several predictors, plot the residuals against the predicted values instead: a flat, even band is fine; a curve means the model is missing a bend. A component-plus-residual plot for each predictor (R: car::crPlots) shows whether one of them bends once the others are held constant.
How Much It Matters: A lot, for what the coefficients mean. Each coefficient is a straight-line slope for one outcome; if depression drops steeply as sleep quality improves from poor to fair and then levels off, one slope misdescribes both ends. Check each outcome's residuals against its predicted values, not just the raw scatterplots.

If It's Violated

Correction Add a squared term for the predictor that bends, centered first. It then enters every outcome's equation, so its joint test asks whether the bend matters for any outcome.
Transformation Log-transform a predictor whose effect levels off, such as income. Its coefficients are then per unit of the log.

3. Equal Spread of Residuals Your scores scatter about equally around the predictions, whatever the prediction.

What It Means: The residuals have about the same spread at every predicted value, so the model misses by about as much for people predicted to score low as for people predicted to score high. The opposite, heteroscedasticity, often looks like a funnel. If the model's form is right, unequal spread doesn't bias the coefficients, though it makes them less precise than they could be; it does make their standard errors, and so the p-values and confidence intervals, unreliable.
Applies To: Residuals
How to Check: Plot the residuals against the predicted values (in the SPSS REGRESSION procedure, *ZRESID against *ZPRED; for MIXED, save them with /SAVE=PRED RESID and plot those). A band of even width is fine; a funnel or fan that widens or narrows is not. Don't use a test such as Breusch–Pagan as a gatekeeper: with small samples it misses real problems, and with large ones it flags trivial ones.
How Much It Matters: Moderate. The multivariate version asks that the residuals' whole covariance matrix, each outcome's spread and the correlations between them, stay the same at every predicted value. The coefficients stay unbiased if it fails, but the standard errors and the joint tests can be off. In the distress study each outcome's residuals spread about as much at low, middle and high predicted values (SDs within 0.3 points of each other across the thirds), an even band; a Breusch–Pagan test, read as a description rather than a gate, agrees.

If It's Violated

Robust Estimation Use heteroscedasticity-consistent (HC3) standard errors for each outcome's coefficients, which don't assume equal spread (Long & Ervin, 2000). For the joint tests, bootstrap them: resample people with replacement and refit.
Transformation Log-transform a right-skewed outcome whose spread grows with its level, such as reaction time. Its coefficients are then on the log scale.

4. Multivariate Normality The variables are roughly bell-shaped together (within each group, when there are groups), or the sample is big enough that it matters less.

What It Means: The variables jointly follow a multivariate normal distribution: each is roughly normal, and every weighted combination of them is too, so the scatterplot of any two is an oval cloud. When the design has groups, as in MANOVA, it applies within each group; in a multivariate regression, to the outcomes' residuals; in factor analysis and SEM, it applies to the observed variables that maximum likelihood is fitted to (for outcomes, after allowing for any predictors the model takes as given). P-values and standard errors rely on it; estimates are usually less affected, and robust standard errors or a scaled test statistic can take its place.
Applies To: Residuals
How to Check: Look at a histogram of each variable (or each outcome's residuals, in a model with predictors) and a scatterplot of each pair, within each group if there are groups, for strong skew, curves or stray points. Multivariate outliers show up as large squared Mahalanobis distances from the means (each group's own, if there are groups). Don't rely on a normality test as a gatekeeper: with small samples it misses real problems, and with large ones it flags trivial ones.
How Much It Matters: Little for the coefficients and the separate t tests with 250 people, and not much for the joint tests either: they work on coefficient estimates, which are close to normal in samples this size. What matters more is a person who is unusual in combination, such as very depressed but barely anxious, which only a multivariate distance shows. Our largest squared Mahalanobis distance among the residuals is 13.32 (using the residual covariance matrix E/(n − k − 1); R's cov() divides by n − 1 instead and gives 13.48), under a common cutoff of 16.27 for three outcomes (χ² at p = .001).

If It's Violated

Robust Estimation Resample people with replacement and refit, for confidence intervals and joint tests that don't assume a residual shape.
Proceed with Caveat Check any case with a large Mahalanobis distance against the raw records, and report the result with and without it.
Transformation Transform a strongly skewed outcome, then interpret its coefficients on the new scale.

5. No Severe Multicollinearity Your predictors aren't so closely related that their effects can't be told apart.

What It Means: Each coefficient is estimated from the part of its predictor that doesn't overlap with the others. When predictors are highly correlated, that part is small, so the coefficients become unstable and their standard errors large, even though the model as a whole may predict well. Perfect overlap, such as including a total score and all of its parts, makes the model impossible to estimate.
Applies To: Predictors
How to Check: Look at the correlations among the predictors and at each predictor's variance inflation factor (VIF). A VIF near 1 means little overlap. Common cutoffs of 5 or 10 are rules of thumb, not tests (O'Brien, 2007); the real question is whether the coefficients you care about are still precise enough to interpret.
How Much It Matters: It depends on the question, exactly as in multiple regression: overlapping predictors inflate every outcome's standard errors for those predictors, and the joint tests along with them. Our VIFs are 1.05 to 1.07, so sleep, support and exercise overlap only slightly.

If It's Violated

Design Change If two predictors measure the same thing, average them into one score or keep the one that fits your question.
Proceed with Caveat If you only care about the predictors jointly, report each outcome's R² and don't lean on individual coefficients.

Limitations

This Is Not Your Test If…

You have only one outcome Multiple Linear Regression →
Your outcomes answer separate questions Multiple Linear Regression →
Your only predictor is a set of groups MANOVA →
You're comparing groups, adjusting for a covariate MANCOVA →
One outcome also predicts another Path Analysis →
Your outcomes are one measure repeated over time Growth-Curve Modeling →
Your people are nested in groups, like clinics Multilevel Regression (HLM) →

What It Doesn't Tell You

Which outcomes. The joint test says support predicts depression, anxiety and stress as a set; it doesn't say which. Here the joint relationship is carried mostly by depression: support's coefficients for anxiety and stress point the same way but are smaller and uncertain. Still, a test that support's coefficient is equal for all three outcomes wasn't significant (p = .075), so the data don't show that support relates to depression more than to the others.
Why. “Holding constant” is arithmetic, not an experiment. Nothing here shows that more support would lower anyone's depression; depressed people may withdraw from others, and anything left out, such as health or money worries, could drive both.
Whether it matters. Significance doesn't make the relationship large: support lines up with about 4% of the variation in the mix of outcomes it predicts best, and adults one SD (0.75 points) apart in support differ by about half a point of depression on a 0–21 scale.
That exercise doesn't matter. Its joint test wasn't significant (p = .126), yet on depression alone its coefficient just reached significance (p = .047). A test spread across three outcomes can miss a relationship confined to one; the result is inconclusive, not evidence that exercise is unrelated to distress.

Common Mistakes

Don'tRun a multivariate regression just because you measured several outcomes.

InsteadUse it when the outcomes answer one question together, such as distress measured three ways. Otherwise run separate regressions (Huang, 2020).
Don'tExpect the multivariate model to change the coefficients or their standard errors.

InsteadKnow they equal the separate regressions. What's new: joint tests, residual correlations, and comparisons of a predictor's coefficients across outcomes.
Don'tTreat a significant joint test as protection for every outcome's t test.

InsteadSay which outcomes the predictor relates to, and correct those follow-up tests if you test several.
Don'tRead a non-significant joint test as “no effect on any outcome.”

InsteadLook at each outcome's coefficient and interval. A real relationship with one outcome can be diluted across outcomes it doesn't touch.
Don'tSwitch to a structural equation model only to get the same estimates.

InsteadMove to path analysis when you have a theory about how the outcomes affect each other, or latent variables.
Don'tAdd up the partial η² values as if they split the variation between them.

InsteadRemember each has its own denominator, so they can sum past 1. Report η² if you want shares of the total.
Don'tReport partial η² under the label η².

InsteadSay which one you report. SPSS prints partial η², which is never smaller than η² and is usually larger, especially with other effects or repeated measures.
Don'tRead a p-value as the chance that the null hypothesis is true.

InsteadRead it as: if the null were true, how often data at least this extreme would turn up.

Result Decoder

Significance and Effect Size Guidelines

Benchmarks Significance (p Value)

p < .001
Very strong evidence against the null hypothesis.
.001 ≤ p < .01
Strong evidence against the null hypothesis.
.01 ≤ p < .05
Statistically significant at α = .05. Moderate evidence against the null hypothesis.
.05 ≤ p < .10
Not significant. Weak evidence; sometimes called a “trend.” Treat it as inconclusive.
p ≥ .10
Not significant. Little evidence against the null, which is not the same as evidence of no effect.

These evidence labels are one common convention for describing p in graded language rather than as significant or not; the boundaries are rough guides, not rules, and other authors draw them differently. Muff et al. (2022) argue for this kind of language of evidence over a single cutoff. α is chosen before the analysis, and p never says how big or important an effect is.

More about statistical significance ↓

Benchmarks Effect Size (Multivariate Partial

\eta^2

)

η²ₚ < .01 Below Small
Smaller than the conventional small benchmark. Whether an effect this size matters depends on the setting.
.01 ≤ η²ₚ < .06 Small
A sliver of the remaining variation lines up with the effect. Hard to see without measuring it.
.06 ≤ η²ₚ < .14 Medium
A noticeable share of the remaining variation lines up with the effect.
η²ₚ ≥ .14 Large
A large share of the remaining variation lines up with the effect.

For one predictor's joint test, partial η² equals Pillai's trace. These benchmarks are Cohen's (1988) for η² in group comparisons. His conventions for regression and for sets of outcomes use f² = .02, .15 and .35 (partial η² ≈ .02, .13 and .26); by those, support (.037) is small and sleep quality (.108) small-to-medium. Read either set loosely, and look at each outcome's coefficient in its own units too.

More about effect size ↓

How It's Reported

Equation Form

Example

social support: Pillai's V = η²_p = .037, F(3, 244) = 3.15, p = .026; on depression: b = −0.71, 95% CI [−1.18, −0.24], t(246) = −2.99, p = .003

Table 1. Means, Standard Deviations and Correlations

Example

	$M$	$SD$	1	2	3	4	5
1. Depression	9.00	2.99
2. Anxiety	7.50	2.51	.65
3. Stress	10.00	3.01	.69	.69
4. Sleep quality	6.00	1.50	−.35	−.27	−.33
5. Social support	4.50	0.75	−.25	−.14	−.15	.20
6. Exercise (hours a week)	3.00	1.00	−.20	−.09	−.09	.20	.15

Note. N = 250. Depression, anxiety and stress are scored 0–21; sleep quality 1–10; social support 1–7. Every correlation is significant at p < .05, two-sided, except exercise with anxiety (p = .147) and with stress (p = .136).

Table 2. Regression Coefficients for Each Outcome

Example

	$b$	$SE$	95% CI	$t$	$p$
Depression: Sleep quality	−0.58	0.12	[−0.82, −0.34]	−4.82	< .001
Depression: Social support	−0.71	0.24	[−1.18, −0.24]	−2.99	.003
Depression: Exercise	−0.36	0.18	[−0.71, −0.01]	−2.00	.047
Anxiety: Sleep quality	−0.42	0.11	[−0.63, −0.21]	−3.96	< .001
Anxiety: Social support	−0.28	0.21	[−0.69, 0.14]	−1.31	.191
Anxiety: Exercise	−0.08	0.16	[−0.39, 0.23]	−0.48	.632
Stress: Sleep quality	−0.63	0.12	[−0.87, −0.38]	−5.06	< .001
Stress: Social support	−0.36	0.25	[−0.84, 0.13]	−1.45	.149
Stress: Exercise	−0.06	0.18	[−0.42, 0.31]	−0.32	.752

Note. Each outcome's coefficients are exactly those of its own multiple regression on the three predictors; df = 246 for every t. R² = .17 for depression, .08 for anxiety and .12 for stress. Intercepts are left out. Social support's coefficient for depression is the one worked through in the Napkin Math.

Table 3. Multivariate Tests for Each Predictor

Example

	Pillai's $V$	$F$	$df_1$	$df_2$	$p$
Sleep quality	.108	9.82	3	244	< .001
Social support	.037	3.15	3	244	.026
Exercise (hours a week)	.023	1.93	3	244	.126

Note. Each row tests one predictor's coefficients on depression, anxiety and stress together, holding the other two predictors constant. With one predictor per test, Pillai's trace, Wilks' Λ (= 1 − V), Hotelling's trace and Roy's largest root all give this same exact F and p, and multivariate partial η² equals V.

Full Text (APA Style)

Example

Adults averaged 9.00 on depression (SD = 2.99), 7.50 on anxiety (SD = 2.51) and 10.00 on stress (SD = 3.01), each on a 0–21 scale, and the three correlated .65 to .69. Holding sleep quality and exercise constant, social support had a small joint relationship with the three outcomes, η²_p = .04. A multivariate regression showed it was statistically significant, Pillai's V = .037, F(3, 244) = 3.15, p = .026. The joint relationship was carried mostly by depression: adults one point higher in support on the 1–7 scale scored 0.71 points lower on depression, 95% CI [−1.18, −0.24], t(246) = −2.99, p = .003, while support's coefficients for anxiety (b = −0.28, p = .191) and stress (b = −0.36, p = .149) were smaller and not significant. A test that support's coefficient was equal for all three outcomes was not significant, F(2, 245) = 2.61, p = .075. Better sleep quality predicted lower scores on all three outcomes, Pillai's V = .108, F(3, 244) = 9.82, p < .001. Exercise's joint test was not significant, V = .023, F(3, 244) = 1.93, p = .126. The predictors together predicted 17% of the variation in depression, 8% in anxiety and 12% in stress, and the outcomes' residuals still correlated .61 to .66.

Give every outcome's coefficients in one table like Table 2 (b, SE, 95% CI, t, p), and say that each outcome's are the same as its own regression's.
Report each predictor's joint test like Table 3: the statistic (Pillai's V is the usual choice), F, both degrees of freedom and p. With one predictor per test all four statistics agree, so say which you used and move on.
Report the residual correlations between the outcomes. They're the part of the model separate regressions don't give.
Write “holding the other predictors constant” and “predicted” or “went with,” not “reduced,” unless the design supports a causal claim.

Reading Each Piece

The equation form above, piece by piece, in this order: how big, how certain, then how surprising.

η²_p = .04	Holding sleep quality and exercise constant, support lines up with about 4% of the variation in the weighted mix of depression, anxiety and stress it predicts best, a mix that here is almost all depression. That's small by the benchmarks above; whether it matters is a question about mental health, not statistics. Benchmarks ↑
b = −0.71, 95% CI [−1.18, −0.24]	Among adults with the same sleep quality and exercise, those one point apart in social support differ by 0.71 points of depression on the 0–21 scale, plausibly anywhere from about a quarter of a point to 1.2 points. One point of support is more than its SD of 0.75.
t(246) = −2.99, p = .003	Support's coefficient for depression is 2.99 standard errors below zero; if support's population coefficient for depression were 0, holding sleep and exercise constant, a result at least this far from zero, in either direction, would turn up about 3 times in 1,000 (if the model's assumptions hold). The 246 is n − k − 1: 250 adults minus 3 predictors and the intercept.
Pillai's V = .037	Holding sleep quality and exercise constant, support lines up with about 3.7% of the variation in the weighted mix of depression, anxiety and stress it predicts best, a mix that here is almost all depression. That isn't a share of the three outcomes' variation added up, which is smaller (about 1.8%). It could reach 1; 0 would mean support predicts none of the three.
F(3, 244) = 3.15	Support's joint test, as an F of 3.15. The 3 and 244 are its degrees of freedom: 3 outcomes, and 250 adults minus 3 predictors minus 3 outcomes. Pillai's, Wilks', Hotelling's and Roy's statistics all give this same F here.
p = .026	If support's population coefficients were 0 for all three outcomes, holding sleep and exercise constant, a joint result at least this strong would turn up about 2.6% of the time (if the model's assumptions hold). That's below .05, so it's statistically significant. Benchmarks ↑

Calculations

The Running Example

Example

The distress study: 250 adults, three related outcomes (depression, anxiety and stress, each 0–21) and the same three predictors for each: sleep quality, social support and weekly exercise.

Null Hypothesis:: support's population coefficients for depression, anxiety and stress are all 0
Alternative:: at least one of the three is not 0

This is the null for social support's joint test, with sleep quality and exercise held constant: among adults like these, support predicts none of the three outcomes. Sleep quality and exercise get joint tests of their own (Table 3), and every coefficient also gets its own two-sided t test, whose null is that one population coefficient is 0.

What Your Data Should Look Like

One row per person, one numeric column for each outcome and one for each predictor. Every outcome is predicted from the same predictors, and anyone missing any value, outcome or predictor, is usually dropped from the whole analysis.

The summary below shows depression only; Table 1 under How It's Reported gives all three outcomes.

Example

participant	depression	anxiety	stress	sleep_quality	social_support	exercise_hours
1	7	4	8	5.6	4.8	2.9
2	9	10	12	6.0	4.1	2.6
3	6	7	7	8.0	5.6	2.8
4	8	7	14	7.4	4.5	3.7
5	10	7	7	4.2	4.0	1.7
6	10	12	12	8.2	4.5	2.6
… 244 more rows

	M	SD	1	2	3
1. Depression	9.00	2.99
2. Sleep quality	6.00	1.50	−.35
3. Social support	4.50	0.75	−.25	.20
4. Exercise (hours a week)	3.00	1.00	−.20	.20	.15

N = 250.

Correlations are Pearson r.

Simulated data.

Download all 250 rows (data.csv)

Napkin Math

The whole test, one idea at a time. Under each heading, a line says why the step matters, and the green note reads the equation in plain English with the distress study's numbers. Collapse all to read just the headings: they tell the whole story on their own.

1. Fit the Same Regression to Each Outcome

A multivariate regression predicts all three outcomes from the same predictors at once. The coefficients it finds are exactly what three separate regressions give, one per outcome.

$\mathbf{B} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{Y}$ = $\mathbf{Y} = \begin{bmatrix} \text{depression} & \text{anxiety} & \text{stress} \end{bmatrix}, \quad \mathbf{X} = \begin{bmatrix} 1 & \text{sleep} & \text{support} & \text{exercise} \end{bmatrix}$

$b_{\text{dep}}$ = −0.71 points
$b_{\text{anx}}$ = −0.28 points
$b_{\text{str}}$ = −0.36 points
$R^2_{\text{dep}}$ = .17
$R^2_{\text{anx}}$ = .08
$R^2_{\text{str}}$ = .12

$\mathbf{Y}$ holds our 250 adults' three scores, one column per outcome, and $\mathbf{X}$ their three predictors plus a column of 1s. Software solves for $\mathbf{B}$ by least squares, the same way each separate regression does, giving a 4 × 3 table of intercepts and coefficients. Support's column entries are −0.71 for depression, −0.28 for anxiety and −0.36 for stress: points per point of support, holding sleep and exercise constant. Each outcome's column, and its $R^2$ , matches its own regression to the last decimal.
2. Give a Range of Plausible Values for Support's Coefficient on Depression

Each coefficient comes with a standard error, the same one its own regression reports. A confidence interval turns it into a range of population values that fit the data.

$b_j \pm t^{*}_{n-k-1} \times SE_{b_j}$ = $-0.71 \pm 1.970 \times 0.238$

= [−1.18, −0.24] points per point of support

$b_j$ is support's −0.71 points of depression per point of support and $SE_{b_j}$ is 0.238. $t^{*}_{246}$ = 1.970 holds the middle 95% of a t distribution with 246 df, so the margin is 1.970 × 0.238 = 0.47, and the range runs from −1.18 to −0.24 points.
3. Test Support on Each Outcome Separately

Each outcome gets the ordinary t test for support's coefficient, the coefficient divided by its standard error. These are the three tests three separate regressions would print.

$t = \dfrac{b_j}{SE_{b_j}}$ = $t_{\text{anxiety}} = \dfrac{-0.28}{0.210}, \quad t_{\text{stress}} = \dfrac{-0.36}{0.247}, \quad t_{\text{depression}} = \dfrac{-0.71}{0.238}$

= −2.99

$p_{\text{dep}}$ = .003
$t_{\text{anx}}$ = −1.31
$t_{\text{str}}$ = −1.45

For depression, −0.71 ÷ 0.238 = −2.99, $p$ = .003. For anxiety, −0.28 ÷ 0.210 = −1.31, $p$ = .191, and for stress, −0.36 ÷ 0.247 = −1.45, $p$ = .149. All three point the same way, but only depression's is significant on its own.
4. Measure How Much the Outcomes' Leftovers Go Together

Here is what separate regressions never show. After the predictors have done their work, people who are more depressed than predicted also tend to be more anxious and more stressed than predicted. That overlap decides how much three t tests really add up to.

$r_{e_1 e_2} = \dfrac{E_{12}}{\sqrt{E_{11}\,E_{22}}}$ = $\dfrac{992.62}{\sqrt{1852.32 \times 1437.44}}$

$r_{e,\text{dep,anx}}$ = .61
$r_{e,\text{dep,str}}$ = .65
$r_{e,\text{anx,str}}$ = .66

$E_{11}$ = 1852.32 and $E_{22}$ = 1437.44 are the squared residuals for depression and anxiety, added up over all 250 people, and $E_{12}$ = 992.62 adds up each person's depression residual times their anxiety residual. 992.62 ÷ √(1852.32 × 1437.44) = .61. Depression with stress is .65, and anxiety with stress .66.
5. Combine the Three Tests, Allowing for the Overlap

If the outcomes' leftovers were unrelated, the evidence from three t tests would simply add up as their squares. Because they're correlated, Hotelling's T² weighs each t against what the others already lead you to expect. When a predictor's t's are all sizable and point the same way, as support's do, they partly repeat each other and T² comes out smaller than the sum. When a predictor relates to one outcome but hardly to its correlated neighbors, T² can come out larger.

$T^2 = \mathbf{t}^\top \mathbf{R}_e^{-1}\, \mathbf{t}$ = $\begin{pmatrix} -2.99 & -1.31 & -1.45 \end{pmatrix} \begin{pmatrix} 1 & .61 & .65 \\ .61 & 1 & .66 \\ .65 & .66 & 1 \end{pmatrix}^{-1} \begin{pmatrix} -2.99 \\ -1.31 \\ -1.45 \end{pmatrix}$

= 9.54

$\textstyle\sum t^2$ = 12.76

$\mathbf{t}$ holds support's three t's and $\mathbf{R}_e$ the residual correlations from the step above; software inverts the matrix, which here means judging each t after allowing for what the other outcomes' t's already predict. Added up plainly, the squared t's come to 2.99² + 1.31² + 1.45² = 12.76. Allowing for the overlap gives $T^2$ = 9.54, smaller because support's three t's are all sizable and point the same way.
6. Find Support's Share of the Variation It Predicts Best

T² measures how strong support's joint signal is. Adding $df_e$ , the amount of leftover noise, gives the total, so Pillai's trace is the signal's share, a number between 0 and 1. It's the share of variation in the one weighted mix of the outcomes that the predictor predicts best, with the weights chosen from the data.

$V = \dfrac{T^2}{T^2 + df_e}$ = $\dfrac{9.54}{9.54 + 246}$

= 0.037

$\Lambda$ = .963

$T^2$ is our 9.54 and $df_e$ is 246, the error degrees of freedom of each outcome's regression. 9.54 ÷ 255.54 = .037. Wilks' Λ is just 1 − .037 = .963 here.
7. Turn That Share into an F Ratio

With one predictor, Pillai's trace converts to an F that is exact, not an approximation. Its degrees of freedom come from the number of outcomes and the error degrees of freedom.

$F = \dfrac{V}{1 - V} \times \dfrac{df_e - p + 1}{p}$ = $\dfrac{.0373}{1 - .0373} \times \dfrac{246 - 3 + 1}{3}$

= 3.15

$df_1$ = 3
$df_2$ = 244

$p$ is our 3 outcomes and $df_e$ = 246. .0373 ÷ .9627 = .0388, and 244 ÷ 3 = 81.33, so $F$ = .0388 × 81.33 = 3.15, on 3 and 244 degrees of freedom.
8. Measure How Strong Support's Joint Relationship Is

F grows with sample size; an effect size doesn't. For a single predictor, multivariate partial eta squared is Pillai's trace itself.

$\eta^2_p = \dfrac{V}{s} = V$ = $\dfrac{.037}{1}$

= .04

$s$ is 1, because one predictor can shift the three outcomes along only one direction. So $\eta^2_p$ = $V$ = .037, about .04.
9. Ask How Surprising This Would Be If Support Predicted None of the Three

Finally, ask how often a joint result at least this strong would turn up if support's population coefficients were 0 for depression, anxiety and stress alike.

$p = P(F_{df_1,\,df_2} \ge F)$ = $P(F_{3,\,244} \ge 3.15)$

= .026

$F_{3,\,244}$ is an F distribution with 3 and 244 df. The chance it lands at or above 3.15 is .026, about 1 in 39. Only large values count against the null, but they catch a relationship with any of the outcomes, in either direction.

Run It in Software

The code below reads the running example's data from data.csv; save it in the same folder you run the code from.

SPSS (version 30, as documented by IBM)

Analyze → General Linear Model → Multivariate.
Move depression, anxiety and stress into Dependent Variables, and sleep_quality, social_support and exercise_hours into Covariate(s). Continuous predictors go in Covariate(s); there are no fixed factors.
Options: check Estimates of effect size and Parameter estimates, then Continue and OK.

Or paste this syntax:

GLM depression anxiety stress WITH sleep_quality social_support exercise_hours
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PRINT=PARAMETER ETASQ
  /CRITERIA=ALPHA(.05)
  /DESIGN=sleep_quality social_support exercise_hours.

Reading the Output

Read Parameter Estimates for each outcome's coefficients (B, its 95% Confidence Interval, t and Sig.), then each predictor's rows in Multivariate Tests: Partial Eta Squared, Value, F and Sig. on Pillai's Trace.

Multivariate Tests (predictor rows; the Intercept rows are left out)

Effect		Value	F	Hypothesis df	Error df	Sig.	Partial Eta Squared
sleep_quality	Pillai's Trace	.108	9.815	3.000	244.000	.000	.108
	Wilks' Lambda	.892	9.815	3.000	244.000	.000	.108
	Hotelling's Trace	.121	9.815	3.000	244.000	.000	.108
	Roy's Largest Root	.121	9.815	3.000	244.000	.000	.108
social_support	Pillai's Trace	.037	3.155	3.000	244.000	.026	.037
	Wilks' Lambda	.963	3.155	3.000	244.000	.026	.037
	Hotelling's Trace	.039	3.155	3.000	244.000	.026	.037
	Roy's Largest Root	.039	3.155	3.000	244.000	.026	.037
exercise_hours	Pillai's Trace	.023	1.927	3.000	244.000	.126	.023
	Wilks' Lambda	.977	1.927	3.000	244.000	.126	.023
	Hotelling's Trace	.024	1.927	3.000	244.000	.126	.023
	Roy's Largest Root	.024	1.927	3.000	244.000	.126	.023

With one predictor per effect, the four statistics give the same exact F and the same partial η² (Pillai's V). Sig. = .000 means p < .001.

Parameter Estimates (the Intercept rows are left out)

Dependent Variable	Parameter	B	Std. Error	t	Sig.	95% CI Lower	Upper	Partial Eta Squared
depression	sleep_quality	-.580	.120	-4.817	.000	-.817	-.343	.086
	social_support	-.714	.238	-2.992	.003	-1.183	-.244	.035
	exercise_hours	-.357	.179	-1.998	.047	-.709	-.005	.016
anxiety	sleep_quality	-.420	.106	-3.962	.000	-.629	-.211	.060
	social_support	-.275	.210	-1.311	.191	-.689	.138	.007
	exercise_hours	-.075	.157	-.480	.632	-.385	.234	.001
stress	sleep_quality	-.630	.124	-5.057	.000	-.875	-.384	.094
	social_support	-.357	.247	-1.446	.149	-.843	.129	.008
	exercise_hours	-.058	.185	-.316	.752	-.422	.306	.000

These are the same numbers Analyze → Regression → Linear gives for each outcome on its own. SPSS also prints Tests of Between-Subjects Effects (each predictor's univariate F, which is the square of its t here); the residual correlations come from saving the residuals (Save → Residuals: Unstandardized) and correlating them.

Output laid out from this page's computed results, not captured from a live run; labels and layout may differ slightly in your version.

R

distress <- read.csv("data.csv")
fit <- lm(cbind(depression, anxiety, stress) ~ sleep_quality + social_support + exercise_hours,
          data = distress)
summary(fit)            # one coefficient table per outcome, identical to three separate lm() fits
car::Anova(fit)         # each predictor's joint test (Pillai's trace by default)
round(cor(residuals(fit)), 2)

Response depression :

Call:
lm(formula = depression ~ sleep_quality + social_support + exercise_hours,
    data = distress)

Residuals:
    Min      1Q  Median      3Q     Max
-6.3368 -1.9288  0.0823  2.0146  6.3177

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)     16.7623     1.2057  13.902  < 2e-16 ***
sleep_quality   -0.5797     0.1204  -4.817 2.55e-06 ***
social_support  -0.7136     0.2385  -2.992  0.00305 **
exercise_hours  -0.3568     0.1786  -1.998  0.04680 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.744 on 246 degrees of freedom
Multiple R-squared:  0.1705,	Adjusted R-squared:  0.1604
F-statistic: 16.85 on 3 and 246 DF,  p-value: 5.48e-10


Response anxiety :

Call:
lm(formula = anxiety ~ sleep_quality + social_support + exercise_hours,
    data = distress)

Residuals:
    Min      1Q  Median      3Q     Max
-6.2197 -1.9616  0.2229  1.6476  7.5980

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)    11.48969    1.06214  10.817  < 2e-16 ***
sleep_quality  -0.42007    0.10602  -3.962 9.73e-05 ***
social_support -0.27537    0.21007  -1.311    0.191
exercise_hours -0.07548    0.15732  -0.480    0.632
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.417 on 246 degrees of freedom
Multiple R-squared:  0.08238,	Adjusted R-squared:  0.07119
F-statistic: 7.362 on 3 and 246 DF,  p-value: 9.602e-05


Response stress :

Call:
lm(formula = stress ~ sleep_quality + social_support + exercise_hours,
    data = distress)

Residuals:
   Min     1Q Median     3Q    Max
-7.156 -2.128  0.232  1.869  8.081

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)    15.56152    1.24718  12.477  < 2e-16 ***
sleep_quality  -0.62958    0.12449  -5.057 8.31e-07 ***
social_support -0.35674    0.24666  -1.446    0.149
exercise_hours -0.05834    0.18472  -0.316    0.752
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.838 on 246 degrees of freedom
Multiple R-squared:  0.1203,	Adjusted R-squared:  0.1096
F-statistic: 11.22 on 3 and 246 DF,  p-value: 6.345e-07



Type II MANOVA Tests: Pillai test statistic
               Df test stat approx F num Df den Df    Pr(>F)
sleep_quality   1  0.107686   9.8155      3    244 3.892e-06 ***
social_support  1  0.037337   3.1545      3    244   0.02552 *
exercise_hours  1  0.023147   1.9272      3    244   0.12573
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
           depression anxiety stress
depression       1.00    0.61   0.65
anxiety          0.61    1.00   0.66
stress           0.65    0.66   1.00

Read each outcome's Estimate, t value and Pr(>|t|), then the car::Anova table for each predictor's joint test (test stat is Pillai's V), then the residual correlations. confint(fit) gives the 95% intervals. Base R's anova(fit) also runs Pillai tests, but they are sequential, so each predictor is adjusted only for those listed before it; car::Anova adjusts each for all the others. For all four statistics, use summary(car::Anova(fit)); with one predictor per test they agree.

To test whether support's coefficient is the same for all three outcomes, run car::linearHypothesis(fit, "social_support", P = cbind(c(1, -1, 0), c(1, 0, -1))): here Pillai = .021, F(2, 245) = 2.61, p = .075. For the predictors as a set, anova(fit, update(fit, . ~ 1)) gives Pillai = .20, F(9, 738) = 5.83, p < .001.

Output from a real run: R 4.6.1 (car 3.1.5), run 2026-09-29.

Python

import pandas as pd
import statsmodels.formula.api as smf
from statsmodels.multivariate.manova import MANOVA

distress = pd.read_csv("data.csv")
predictors = "sleep_quality + social_support + exercise_hours"

# Each predictor's joint test across the three outcomes
print(MANOVA.from_formula("depression + anxiety + stress ~ " + predictors, data=distress).mv_test())

# Coefficients: one ordinary regression per outcome gives exactly the same numbers
fits = {o: smf.ols(f"{o} ~ {predictors}", data=distress).fit() for o in ["depression", "anxiety", "stress"]}
print(fits["depression"].summary().tables[1])
print(pd.DataFrame({o: f.resid for o, f in fits.items()}).corr().round(2))

-------------------------------------------------------------
     sleep_quality      Value  Num DF  Den DF  F Value Pr > F
-------------------------------------------------------------
          Wilks' lambda 0.8923 3.0000 244.0000  9.8155 0.0000
         Pillai's trace 0.1077 3.0000 244.0000  9.8155 0.0000
 Hotelling-Lawley trace 0.1207 3.0000 244.0000  9.8155 0.0000
    Roy's greatest root 0.1207 3.0000 244.0000  9.8155 0.0000
-------------------------------------------------------------

-------------------------------------------------------------
     social_support     Value  Num DF  Den DF  F Value Pr > F
-------------------------------------------------------------
          Wilks' lambda 0.9627 3.0000 244.0000  3.1545 0.0255
         Pillai's trace 0.0373 3.0000 244.0000  3.1545 0.0255
 Hotelling-Lawley trace 0.0388 3.0000 244.0000  3.1545 0.0255
    Roy's greatest root 0.0388 3.0000 244.0000  3.1545 0.0255
-------------------------------------------------------------

-------------------------------------------------------------
     exercise_hours     Value  Num DF  Den DF  F Value Pr > F
-------------------------------------------------------------
          Wilks' lambda 0.9769 3.0000 244.0000  1.9272 0.1257
         Pillai's trace 0.0231 3.0000 244.0000  1.9272 0.1257
 Hotelling-Lawley trace 0.0237 3.0000 244.0000  1.9272 0.1257
    Roy's greatest root 0.0237 3.0000 244.0000  1.9272 0.1257
=============================================================

==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
Intercept         16.7623      1.206     13.902      0.000      14.387      19.137
sleep_quality     -0.5797      0.120     -4.817      0.000      -0.817      -0.343
social_support    -0.7136      0.238     -2.992      0.003      -1.183      -0.244
exercise_hours    -0.3568      0.179     -1.998      0.047      -0.709      -0.005
==================================================================================
            depression  anxiety  stress
depression        1.00     0.61    0.65
anxiety           0.61     1.00    0.66
stress            0.65     0.66    1.00

The Intercept block that mv_test() prints first is left out, and only depression's coefficient table is shown; anxiety's and stress's match Table 2. Despite its name, statsmodels' MANOVA class accepts continuous predictors and gives the same joint tests as R's car::Anova.

Output from a real run: Python 3.9.6 (statsmodels 0.14.6, pandas 2.3.3), run 2026-10-03; intercept tests trimmed.

Impact

Effect Size

A multivariate regression has effect sizes at two levels. Each coefficient b gives how one outcome's prediction changes with a predictor, in that outcome's own units, and each outcome's R² is the share of its variation the predictors predict, exactly as in separate regressions. For a predictor's joint relationship with the outcomes, multivariate partial eta squared is the share of variation, after the other predictors, in the weighted mix of the outcomes that the predictor predicts best. The weights are chosen from the data. It isn't a share of all the outcomes' variation added up, which is smaller. With one predictor per test it equals Pillai's trace.

\eta^2_p = V = \dfrac{T^2}{T^2 + df_e} = \dfrac{9.54}{9.54 + 246} = .04

Holding sleep quality and exercise constant, social support lines up with about 4% of the variation in its best-predicted mix of depression, anxiety and stress, a mix that here is almost all depression: small by the benchmarks. Its share of the three outcomes' leftover variation added up is about 2%. Sleep quality's η²_p is larger, .11, and exercise's .02.

Like any variance share, η²_p can't go below 0, and because the weights are fitted to the sample it comes out above 0 even when a predictor relates to nothing: with 3 outcomes and 246 error df it averages about .01 (3 ÷ 247) under the null.

A single joint number hides where the relationship is, so also give the coefficients. Adults one point higher in support on the 1–7 scale scored 0.71 points lower on depression, 0.28 lower on anxiety and 0.36 lower on stress, on 0–21 scales, holding the other predictors constant. Support's SD is 0.75, so adults one SD apart in support differ by about half a point of depression (0.71 × 0.75 = 0.53; standardized β = −.18). Benchmarks are conventions. Whether half a point of depression matters is a question about mental health, not statistics.

Benchmarks Multivariate Partial

\eta^2

at a Glance

Below Small

η²ₚ < .01

Small

.01 ≤ η²ₚ < .06

Medium

.06 ≤ η²ₚ < .14

Large

η²ₚ ≥ .14

Full benchmarks and caveats ↑

Statistical Significance

A multivariate regression answers two kinds of significance question. Each coefficient's t test asks: if this predictor's population coefficient for this outcome were 0, holding the others constant, how often would its coefficient land at least this far from zero, in either direction? Those are the same tests separate regressions give. Each predictor's joint test asks: if its coefficients were 0 for every outcome, how often would they, taken together, be at least this far from zero? Both answers assume the model's assumptions hold. Below the conventional α = .05, a result is called statistically significant.

If support's population coefficients were 0 for all three outcomes, holding sleep and exercise constant, a joint result at least as strong as ours would turn up about 2.6% of the time (p = .026). On its own, support's depression coefficient would land at least this far from zero, in either direction, about 0.3% of the time (p = .003); its anxiety and stress coefficients weren't significant.

The joint test isn't just the separate tests added up. Our three outcomes' leftovers correlate about .6, so the tests overlap. Support's three t's are all sizable and point the same way, so they partly repeat each other: squared and summed they come to 12.76, but T² is 9.54. Exercise shows the reverse. Its t is −2.00 on depression but near 0 on anxiety and stress, even though their leftovers go with depression's. That mismatch is informative in itself, so its T² (5.83) is larger than its summed squared t's (4.32). Exercise's joint test still wasn't significant (p = .126, against p = .047 on depression alone), because the joint test spreads its attention across three outcomes. Decide in advance which question you're asking.

Two more tests come from fitting the outcomes together. To ask whether the predictors as a set relate to the outcomes as a set, compare the model with an intercept-only model (in R, anova(fit, update(fit, . ~ 1))), the multivariate counterpart of each regression's overall F. Here Pillai's V = .20, F(9, 738) = 5.83, p < .001. To ask whether a predictor's coefficient differs across outcomes on the same scale, test that they're equal: for support, F(2, 245) = 2.61, p = .075, so the data can't tell its depression coefficient from its anxiety and stress ones.

p is not the probability that a predictor is useless, and it says nothing about how much it matters or on which outcome. That's what the coefficients, intervals and effect sizes are for. With one predictor per test, Pillai's trace, Wilks' Λ, Hotelling's trace and Roy's largest root all give the same p, so the usual MANOVA debate about which to report doesn't arise.

Benchmarks p at a Glance

Very Strong

p < .001

Strong

.001 ≤ p < .01

Significant

.01 ≤ p < .05

Weak

.05 ≤ p < .10

Little

p ≥ .10

Strength of evidence against the null hypothesis.

Full benchmarks and caveats ↑

Power and Sample Size

Power is the chance that a study detects an effect that's really there. A real effect doesn't guarantee a significant result: with too few people, it can easily come out non-significant. For a predictor's joint test, power depends on the size of its coefficient on every outcome, how strongly the outcomes' residuals correlate, how many outcomes and predictors there are, how many people you have, and α.

The pattern matters, not just the sizes. Correlated outcomes that all move the same way carry overlapping evidence, so adding a third outcome that repeats the first two buys little power, while each extra outcome costs a degree of freedom. There's no calculator on this page. Plan by simulation: write down the smallest coefficients that would matter for each outcome, plausible residual SDs and correlations, generate many datasets, and count how often the joint test comes out significant.

We simulated 4,000 studies from the fitted model, taking our coefficients and residual correlations as the truth. With 250 adults, support's joint test came out significant 74% of the time; it took about 300 adults to reach 80%. If support's true coefficients were half as large, 250 adults would detect them only 22% of the time, and it would take about 1,150.

Don't plan the next study from our observed effect alone: with a p of .026 it's a noisy estimate, and the true effect could easily be smaller. Also plan for the question you'll actually interpret. If what matters is support's coefficient on depression alone, the power you need is for that one coefficient's t test, as in multiple regression.

Researchers usually aim for at least 80% power. Plan it before collecting data, using the smallest effects that would matter.

Other Info

Follow-Up and Post Hoc Tests

Post Hoc Tests

Holm Correction: Use after a significant joint test, to say which outcomes the predictor relates to: take its coefficient's t test for each outcome (Table 2) and Holm-correct across the outcomes. For social support, depression holds up (Holm-adjusted p = .009); anxiety and stress don't (p = .299 for both).

Related Methods

General CasePath Analysis: a multivariate regression is a path model in which every predictor has an arrow to every outcome and the outcomes' residuals are free to correlate. That model is saturated (zero degrees of freedom), so it fits perfectly and gives the same coefficients. Path analysis adds value once you leave some paths out or let outcomes predict each other.
Special CaseMANCOVA: multivariate regression whose predictors are k − 1 group dummy variables plus a continuous covariate. The multivariate test of the dummies, over and above the covariate, is the MANCOVA.
Special CaseMANOVA: multivariate regression whose only predictors are k − 1 dummy variables for the groups. The multivariate test of the dummies is the MANOVA.
Special CaseMultiple Linear Regression: the one-outcome case. Fit it to each outcome separately and you get the same coefficients; what the multivariate model adds is tests that consider all the outcomes together.
AlternativeMultilevel Regression (HLM): An alternative if “Independent Observations” doesn't hold. A multilevel model with the outcomes stacked, which estimates the same coefficients while accounting for the clustering.

References

Hotelling, H. (1931). The generalization of Student's ratio. The Annals of Mathematical Statistics, 2(3), 360–378. https://doi.org/10.1214/aoms/1177732979
Pillai, K. C. S. (1955). Some new test criteria in multivariate analysis. The Annals of Mathematical Statistics, 26(1), 117–121. https://doi.org/10.1214/aoms/1177728599
“Three new test criteria are proposed for overall tests of hypotheses in multivariate analysis.” (Abstract)
Wilks, S. S. (1932). Certain generalizations in the analysis of variance. Biometrika, 24(3–4), 471–494. https://doi.org/10.1093/biomet/24.3-4.471
Friendly, M., Monette, G., & Fox, J. (2013). Elliptical insights: Understanding statistical methods through elliptical geometry. Statistical Science, 28(1), 1–39. https://doi.org/10.1214/12-STS402
“all tests of linear hypotheses (for null effects) can be represented in the form of a general linear test” (Section 5, p. 19)
Cohen, J., Cohen, P., West, S. G., & Aiken, L. S. (2003). Applied multiple regression/correlation analysis for the behavioral sciences (3rd ed.). Lawrence Erlbaum Associates. https://doi.org/10.4324/9780203774441
Huang, F. L. (2020). MANOVA: A procedure whose time has passed? Gifted Child Quarterly, 64(1), 56–60. https://doi.org/10.1177/0016986219887200
“Just because a researcher has multiple dependent variables of interest does not mean that a MANOVA should be used at all.” (Abstract)
Long, J. S., & Ervin, L. H. (2000). Using heteroscedasticity consistent standard errors in the linear regression model. The American Statistician, 54(3), 217–224. https://doi.org/10.1080/00031305.2000.10474549
Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum Associates.
Olson, C. L. (1974). Comparative robustness of six tests in multivariate analysis of variance. Journal of the American Statistical Association, 69(348), 894–908. https://doi.org/10.1080/01621459.1974.10480224
“the Pillai-Bartlett trace test may be recommended as the most robust of the MANOVA tests, with adequate power to detect true differences in a variety of situations” (Abstract, p. 894)
Cohen, J. (1973). Eta-squared and partial eta-squared in fixed factor ANOVA designs. Educational and Psychological Measurement, 33(1), 107–112. https://doi.org/10.1177/001316447303300111
Pierce, C. A., Block, R. A., & Aguinis, H. (2004). Cautionary note on reporting eta-squared values from multifactor ANOVA designs. Educational and Psychological Measurement, 64(6), 916–924. https://doi.org/10.1177/0013164404264848
Lakens, D. (2013). Calculating and reporting effect sizes to facilitate cumulative science: A practical primer for t-tests and ANOVAs. Frontiers in Psychology, 4, 863. https://doi.org/10.3389/fpsyg.2013.00863
“using these benchmarks when interpreting the η2p effect size in designs that include covariates or repeated measures is not consistent with the considerations upon which the benchmarks were based” (Eta-squared in between and within-subjects comparisons)
Muff, S., Nilsen, E. B., O'Hara, R. B., & Nater, C. R. (2022). Rewriting results sections in the language of evidence. Trends in Ecology & Evolution, 37(3), 203–210. https://doi.org/10.1016/j.tree.2021.10.009
“we suggest language of evidence that allows for a more nuanced approach to communicate scientific findings as a simple and intuitive alternative to statistical significance testing” (Abstract)
Williams, M. N., Gómez Grajales, C. A., & Kurkiewicz, D. (2013). Assumptions of multiple regression: Correcting two misconceptions. Practical Assessment, Research & Evaluation, 18(11), 1–14. https://doi.org/10.7275/55hn-wk47
“ordinary least squares estimates will be unbiased and consistent as long as the errors are independent, but will not be efficient” (p. 9)
Hayes, A. F., & Cai, L. (2007). Using heteroskedasticity-consistent standard error estimators in OLS regression: An introduction and software implementation. Behavior Research Methods, 39(4), 709–722. https://doi.org/10.3758/BF03192961
O'Brien, R. M. (2007). A caution regarding rules of thumb for variance inflation factors. Quality & Quantity, 41(5), 673–690. https://doi.org/10.1007/s11135-006-9018-6