Analysis of Variance (ANOVA)

How does it work?

Statistical methods are based on several fundamental concepts, the most central of which is to consider the information available (in the form of data) resulting from a random process.

As such, the data represent a random sample of a totally or conceptually accessible population.

Then, statistical inference allows to infer the properties of a population based on the observed sample. This includes deriving estimates and testing hypotheses.

Source: luminousmen

Hypothesis testing

• In general (scientific) hypotheses can be translated into a set of (non-overlapping idealized) statistical hypotheses:

• In a hypothesis test, the statement being tested is called the null hypothesis \(\color{#373895}{H_0}\). A hypothesis test is designed to assess the strength of the evidence against the null hypothesis.
• The alternative hypothesis \(\color{#373895}{H_a}\) is the statement we hope or suspect to be true instead of \(\color{#373895}{H_0}\).
• Each hypothesis excludes the other, so that one can exclude one in favor of the other using the data.
• Example: a drug represses the progression of cancer

Hypothesis testing

Hypothesis testing

What are p-values?

In statistics, the p-value is defined as the probability of observing a test statistic that is at least as extreme as actually observed, assuming that the null hypothesis \(H_0\) is true.

Informally, a p-value can be understood as a measure of plausibility of the null hypothesis given the data. A small p-value indicates strong evidence against \(H_0\).

When the p-value is small enough (i.e., smaller than the significance level \(\alpha\)), the test based on the null and alternative hypotheses is considered significant, meaning we reject the null hypothesis in favor of the alternative. This is generally what we want because it “verifies” our (research) hypothesis.

When the p-value is not small enough, with the available data, we cannot reject the null hypothesis, so nothing can be concluded. 🤔

The obtained p-value summarizes the incompatibility between the data and the model constructed under the set of assumptions.

👋 From the British Medical Journal.

How to understand p-values?

👋 If you want to know more have a look here.

P-values may be controversial

P-values have been misused many times because understanding what they mean is not intuitive!

👋 If you want to know more have a look here.

Quick review: Normal distribution

Two-sample location tests

In practice, we often encounter problems where our goal is to compare the means (or locations) of two samples. For example,

1. A scientist is interested in comparing the vaccine efficacy of the Pfizer-BioNTech and the Moderna vaccine.

2. A bank wants to know which of two proposed plans will most increase the use of its credit cards.

3. A psychologist wants to compare male and female college students’ impression on a selected webpage.

We will discuss three two-sample location tests:

1. Two independent sample Student’s t-test

2. Two independent sample Welch’s t-test

3. Two independent sample Mann-Whitney-Wilcoxon test

Two independent sample Student’s t-test

This test considers the following assumed model for group A and B

where \(g=A,B\), \(i=1,...,n_{g}\), \(\varepsilon_{i(g)} \overset{iid}{\sim} N(0,\color{#eb078e}{\sigma^{2}}\)\()\) and \(\sum n_{g}\delta_{g} =0\).

📝 \(\color{#6A5ACD}{n_A}\) \(=\) sample size of group A, \(\color{#6A5ACD}{\mu_{A} = \mu + \delta_A}\) \(=\) population mean of group A, \(\color{#06bcf1}{n_B}\) and \(\color{#06bcf1}{\mu_{B} = \mu + \delta_B}\) are similarly defined for group B.

Hypotheses:

Test statistic’s distribution under \(H_0\):

\[ \color{#b4b4b4}{T = \frac{(\overline{X}_{A}-\overline{X}_{B})-\mu_0}{s_{p}\sqrt{n_{A}^{-1}+n_{B}^{-1}}} \ {\underset{H_0}{\sim}} \ \text{Student}(n_{A}+n_{B}-2).} \]

Discussion — Student’s t-test

Python function:

from scipy import stats

stats.ttest_ind(a = ..., b = ..., alternative = ..., equal_var = True)

This test strongly relies on the assumed absence of outliers. If outliers appear to be present the Mann-Whitney-Wilcoxon test (see later) is (probably) a better option.

For moderate and small sample sizes, the sample distribution should be at least approximately normal with no strong skewness to ensure the reliability of the test.

In practice, the assumption of equal variance is hard to verify so we recommend to avoid this test in practice.

Two independent sample Welch’s t-test

This test considers the following assumed model for group A and B

where \(g=A,B\), \(i=1,...,n_{g}\), \(\varepsilon_{i(g)} \overset{iid}{\sim} N(0,\color{#eb078e}{\sigma_g^{2}}\)\()\) and \(\sum n_{g}\delta_{g} =0\).

📝 \(\color{#6A5ACD}{n_A}\) \(=\) sample size of group A, \(\color{#6A5ACD}{\mu_{A} = \mu + \delta_A}\) \(=\) population mean of group A, \(\color{#06bcf1}{n_B}\) and \(\color{#06bcf1}{\mu_{B} = \mu + \delta_B}\) are similarly defined for group B.

Hypotheses:

Test statistic’s distribution under \(H_0\):

\[ \color{#b4b4b4}{T = \frac{(\overline{X}_{A}-\overline{X}_{B})-\mu_0}{\sqrt{s^2_A/n_{A} + s^2_B/n_{B}}} \ {\underset{H_0}{\sim}} \ \text{Student}(df).} \]

Discussion — Welch’s t-test

Python function:

from scipy import stats

stats.ttest_ind(a = ..., b = ..., alternative = ..., equal_var = False)

This test strongly relies on the assumed absence of outliers. If outliers appear to be present the Mann-Whitney-Wilcoxon test (see later) is (probably) a better option.

For moderate and small sample sizes, the sample distribution should be at least approximately normal with no strong skewness to ensure the reliability of the test.

This test does not require the variances of the two groups to be equal. If the variances of the two groups are the same (which is rather unlikely in practice), the Welch’s t-test loses a little bit of power compared to the Student’s t-test.

The computation of \(df\) (i.e. the degrees of freedom of the distribution under the null) is beyond the scope of this class.

Mann-Whitney-Wilcoxon test

This test considers the following assumed model for group A and B

where \(g=A,B\), \(i=1,...,n_{g}\), \(\varepsilon_{i(g)} \overset{iid}{\sim} N(0,\color{#eb078e}{\sigma^{2}}\)\()\) and \(\sum n_{g}\delta_{g} =0\).

📝 \(\color{#6A5ACD}{n_A}\) \(=\) sample size of group A, \(\color{#6A5ACD}{\theta_{A} = \theta + \delta_A}\) \(=\) population mean of group A, \(\color{#06bcf1}{n_B}\) and \(\color{#06bcf1}{\theta_{B} = \theta + \delta_B}\) are similarly defined for group B.

Hypotheses:

Test statistic’s distribution under \(H_0\):

\[ \color{#b4b4b4}{Z = \frac{\sum_{i=1}^{n_{B}}R_{i(g)}-[n_{B}(n_{A}+n_{B}+1)/2]}{\sqrt{n_{A}n_{B}(n_{A}+n_{B}+1)/12}},} \] where \(\color{#b4b4b4}{R_{i(g)}}\) denotes the global rank of the \(\color{#b4b4b4}{i}\)-th observation of group \(\color{#b4b4b4}{g}\).

Discussion – Mann–Whitney–Wilcoxon

Python function:

from scipy import stats

stats.mannwhitneyu(a = ..., b = ...)

This test is “robust” in the sense that (unlike the t-tests) it is not overly affected by outliers.

For the Mann–Whitney–Wilcoxon test to be comparable to the t-tests (i.e. testing for the mean) we need to assume symmetric distributions and equality in variances .

Then, we have \(\color{#6A5ACD}{\theta_A = \mu_A}\) and \(\color{#06bcf1}{\theta_B = \mu_B}\).

Compared to the t-tests, the Mann–Whitney–Wilcoxon test is less powerful if their requirements (Gaussian and possibly same variances) are met.

The distribution of this method under the null is complicated and can be obtained by different methods (e.g. exact, asymptotic normal, …).
The details are beyond the scope of this class.

Comparing diets A and B

Graph

Is there a problem in doing many tests?

Are jelly beans causing acne? Maybe… but why only green ones? 🤨

Source: xkcd

Are jelly beans causing acne?

Source: xkcd

Maybe a specific color?

Source: xkcd

Maybe a specific color?

Source: xkcd

And finally…

Source: xkcd

👋 If you want to know more about this comic, have a look here.

Multiple testing can be dangerous!

Possible solutions

Possible solutions

Multiple-sample location tests

Fisher’s one-way ANOVA

This test considers the following assumed model for G groups \[X_{i(g)} = {\color{#eb078e}{\mu_{g}}} + \varepsilon_{i(g)} = \mu + {\color{#eb078e}{\delta_{g}}} + \varepsilon_{i(g)},\] where \(g=1,\ldots, G\), \(i=1,...,n_{g}\), \(\varepsilon_{i(g)} \overset{iid}{\sim} N(0,{\color{#eb078e}{\sigma^{2}}})\) and \(\sum n_{g}\delta_{g}=0\).

📝 \(n_i =\) sample size of group \(i\), \(\mu_i = \mu + \delta_i =\) population mean of group \(i\), \(i=1,\ldots, G\).

Hypotheses: \[H_0: {\color{#6A5ACD}{\mu_1}} {\color{#eb078e}{=}} {\color{#06bcf1}{\mu_2}} {\color{#eb078e}{=}} \ldots {\color{#eb078e}{=}} {\color{#8bb174}{\mu_G}} \ \ \ \ \text{and} \ \ \ \ H_a: \mu_i {\color{#eb078e}{\neq}} \mu_j \ \ \text{for at least one pair of} \ \ (i,j).\]

Test statistic’s distribution under \(H_0\):

\(\color{#b4b4b4}{F = \frac{N s^2_{\overline{X}}}{s_p^2} \ {\underset{H_0}{\sim}} \ \text{Fisher}(G-1, N-G),}\) where \(\small \color{#b4b4b4}{s^2_{\overline{X}} = \frac{1}{G-1} \sum_{g=1}^G \frac{n_g}{N} (\overline{X}_g - \overline{X})^2}\),

\(\small \color{#b4b4b4}{s_p^2 = \frac{1}{N-G} \sum_{g=1}^G (n_g-1)s_g^2}\), \(\small \color{#b4b4b4}{N = \sum_{g=1}^G n_g}\), and \(\small \color{#b4b4b4}{\overline{X} = \frac{1}{N} \sum_{g=1}^G n_g \overline{X}_g}\).

Discussion - Fisher’s one-way ANOVA

Python function:

from scipy import stats

stats.f_oneway(group_A, group_B, group_C)

Welch’s one-way ANOVA

This test considers the following assumed model for G groups \[X_{i(g)} = {\color{#eb078e}{\mu_{g}}} + \varepsilon_{i(g)} = \mu + {\color{#eb078e}{\delta_{g}}} + \varepsilon_{i(g)},\] where \(g=1,\ldots, G\), \(i=1,...,n_{g}\), \(\varepsilon_{i(g)} \overset{iid}{\sim} N(0,{\color{#eb078e}{\sigma_g^{2}}})\) and \(\sum n_{g}\delta_{g}=0\).

Hypotheses: \[H_0: {\color{#6A5ACD}{\mu_1}} {\color{#eb078e}{=}} {\color{#06bcf1}{\mu_2}} {\color{#eb078e}{=}} \ldots {\color{#eb078e}{=}} {\color{#8bb174}{\mu_G}} \ \ \ \ \text{and} \ \ \ \ H_a: \mu_i {\color{#eb078e}{\neq}} \mu_j \ \ \text{for at least one pair of} \ \ (i,j).\]

Test statistic’s distribution under \(H_0\): \[\color{#b4b4b4}{F^* = \frac{s^{*^2}_{\overline{X}}}{1+\frac{2(G-2)}{3\Delta}} \ {\underset{H_0}{\sim}} \ \text{Fisher}(G-1, \Delta),}\] where \(\color{#b4b4b4}{s^{*^2}_{\overline{X}} = \frac{1}{G-1} \sum_{g=1}^G w_g (\overline{X}_g - \overline{X}^*)^2}\), \(\color{#b4b4b4}{\Delta = [\frac{3}{G^2-1} \sum_{g=1}^G \frac{1}{n_g} (1-\frac{w_g}{\sum_{g=1}^G w_g})]^{-1}}\), \(\color{#b4b4b4}{w_g = \frac{n_g}{s_g^2}}\),

and \(\small \color{#b4b4b4}{\overline{X}^* = \sum_{g=1}^G \frac{w_g\overline{X}_g}{\sum_{g=1}^G w_g}}\).

Discussion - Welch’s one-way ANOVA

Python function:

from statsmodels.stats.oneway import anova_oneway

anova_oneway(data, groups=groups, use_var='unequal', welch_correction=True)

This test strongly relies on the assumed absence of outliers. If outliers appear to be present the Kruskal-Wallis test (see later) is (probably) a better option.

For moderate and small sample sizes, the sample distribution should be at least approximately normal with no strong skewness to ensure the reliability of the test.

This test does not require the variances of the groups to be equal. If the variances of all the groups are the same (which is rather unlikely in practice), the Welch’s one-way ANOVA loses a little bit of power compared to the Fisher’s one-way ANOVA.

Kruskal-Wallis test

This test considers the following assumed model for G groups \[X_{i(g)} = {\color{#eb078e}{\theta_{g}}} + \varepsilon_{i(g)} = \theta + {\color{#eb078e}{\delta_{g}}} + \varepsilon_{i(g)},\] where \(g=1,\ldots, G\), \(i=1,...,n_{g}\), \(\varepsilon_{i(g)} \overset{iid}{\sim} N(0,{\color{#eb078e}{\sigma^{2}}})\) and \(\sum n_{g}\delta_{g}=0\).

📝 \(n_i =\) sample size of group \(i\), \(\theta_i = \theta + \delta_i =\) population location of group \(i\), \(i=1,\ldots, G\).

Hypotheses: \[H_0: {\color{#6A5ACD}{\theta_1}} {\color{#eb078e}{=}} {\color{#06bcf1}{\theta_2}} {\color{#eb078e}{=}} \ldots {\color{#eb078e}{=}} {\color{#8bb174}{\theta_G}} \ \ \ \ \text{and} \ \ \ \ H_a: \theta_i {\color{#eb078e}{\neq}} \theta_j \ \ \text{for at least one pair of} \ \ (i,j).\]

\(\color{#b4b4b4}{\text{Test statistic's distribution under } H_0: \quad H = \frac{\frac{12}{N(N+1)} \sum_{g=1}^G \frac{\overline{R}_g}{n_g}-3(N-1)}{1-\frac{\sum_{v=1}^V{t_v^3 - t_v}}{N^3 - N}} \ {\underset{H_0}{\sim}} \ \mathcal{X}(G-1)}\)

\(\color{#b4b4b4}{\text{where } \overline{R}_g = \frac{1}{n_g} \sum_{i=1}^{n_g} R_{i(g)} \ \text{with } R_{i(g)} \ \text{denoting the global rank of the } i^{th} \ \text{observation of group } g,}\)

\(\color{#b4b4b4}{V \ \text{is the number of different values/levels in } X \ \text{and } t_v \ \text{denotes the number of times a given}}\)

\(\color{#b4b4b4}{\text{value/level occurred in } X.}\)

Discussion - Kruskal-Wallis test

Python function:

from scipy import stats

stats.kruskal(group_A, group_B, group_C)

• This test is “robust” in the sense that (unlike the one-way ANOVA) it is not overly affected by outliers.
• For the Kruskal-Wallis test to be comparable to the one-way ANOVAs (i.e. testing for the mean) we need to assume: (1) The distributions are symmetric, (2) the variances are the same. Then, we have \(\color{#eb078e}{\theta_i = \mu_i,\hspace{0.1cm} i=1,\ldots,G}\).
• Compared to the one-way ANOVAs, the Kruskal-Wallis test is less powerful if their requirements (Gaussian and possibly same variances) are met.

Exercise: Comparing diets A, B and C

Graph

Exercise: Comparing diets A, B and C

Import

# Import data
import pandas as pd
diet = pd.read_csv("https://raw.githubusercontent.com/ELSTE-Master/Data-Science/main/Data/diet.csv")
diet["weight.loss"] = diet["initial.weight"] - diet["final.weight"]

# Variable of interest
dietA = diet["weight.loss"][diet["diet.type"]=="A"]
dietB = diet["weight.loss"][diet["diet.type"]=="B"]
dietC = diet["weight.loss"][diet["diet.type"]=="C"]

# Create data frame
dat = pd.DataFrame({
    "response": list(dietA) + list(dietB) + list(dietC),
    "groups": (["A"] * len(dietA)) + (["B"] * len(dietB)) + (["C"] * len(dietC))
})