Definition

One-way ANOVA, also known as single-factor ANOVA or single-factor factorial model, allows us to study whether there are significant differences between the mean of a continuous random variable at the different levels of another qualitative variable or factor, when the data are not paired. It is an extension of the independent t-tests for more than two groups.

The hypotheses tested in a one-way ANOVA are:

• $H_0$: There are no differences between the means of the different groups: ${\mu }_1={\mu }_2...={\mu }_k=\mu$

• $H_1$: At least one pair of means are significantly different from each other.

Another way to state the hypotheses of an ANOVA is as follows: if $\mu$ is considered as the expected value for any observation of the population (the mean of all observations without taking into account the different levels), and ${\alpha }_i$ is the effect introduced by level $i$. The mean of a certain level (${\mu }_i$) can be defined as:

• $H_0$: No level introduces an effect on the total mean: ${\alpha }_1={\alpha }_2=...{\alpha }_k=0$

• $H_1$: At least one level introduces an effect that shifts its mean: Some ${\alpha }_i\ne 0$
  
  

As mentioned before, the difference between means is detected through the study of the variance between groups and within groups. To achieve this, ANOVA requires a decomposition of the variance based on the following idea:

which is equivalent to:

which is equivalent to:

To be able to calculate the different variances, first the Sums of Squares (SS) must be obtained:

• Total sum of squares (TSS): measures the total variability of the data, it is defined as the sum of the squares of the differences of each observation with respect to the general mean of all observations. The degrees of freedom of the total sum of squares are equal to the total number of observations minus one (N-1).

• Sum of squares due to treatment (SST): measures the variability in the data associated with the effect of the factor on the mean (the difference of the means between the different levels or groups). It is obtained as the sum of the squares of the deviations of the mean of each group with respect to the general mean, weighting each squared difference by the number of observations of each group. The corresponding degrees of freedom are equal to the number of levels of the factor minus one (k-1).

• Sum of squares of errors (SSE): measures the variability within each level, that is, the variability that is not due to the qualitative variable. It is calculated as the sum of the squares of the deviations of each observation with respect to the mean of the level to which it belongs. The degrees of freedom assigned to the sum of residual squares are equivalent to the difference between the total degrees of freedom and the degrees of freedom of the factor, or what is the same (N-k). In statistics, the term error or residual is used since it is considered that this is the variability that the data show due to measurement errors. From a biological point of view, it makes more sense to call it sum of squares within groups since it is known that the observed variability is not only due to measurement errors, but to the many factors that are not controlled and that affect biological processes.

Once the Sum of Squares is decomposed, the decomposition of the variance can be obtained by dividing the Sum of Squares by the respective degrees of freedom. Strictly speaking, the quotient between the Sum of Squares and its corresponding degrees of freedom is called Mean Sum of Squares and can be used as an estimator of the variance:

• ${\hat{S}}_T^2=\frac{TSS}{N-1}=$ Total Mean Squares = Total Quasi-variance (total sample variance)
  
  

• ${\hat{S}}_t^2=\frac{SST}{k-1}=$ Factor Mean Squares = Between-group variance (variance between the means of the different levels)
  
  

• ${\hat{S}}_E^2=\frac{SSE}{N-k}=$ Error Mean Squares = Within-group variance (variance within the levels, known as residual or error variance)
  
  

Note: ANOVA is defined as analysis of variance, but in a strict sense, it is an analysis of the Mean Sum of Squares.

Once the estimation of the variance is decomposed, the $F_{ratio}$ statistic is obtained by dividing the between-group variance by the within-group variance:

Since by definition the $F_{ratio}$ statistic follows a F Fisher-Snedecor distribution with $k-1$ and $N-t$ degrees of freedom, the probability of obtaining values equal to or more extreme than those observed (p-values) can be known.

Conditions for one-way ANOVA with independent data

The conditions for a one-way ANOVA with independent data are:

• Independence: the observations must be random and independent of each other. The groups (levels of the factor) must be independent of each other.

• Normal distribution of each of the levels or groups: the variable under study must be normally distributed in each of the groups, this condition being less strict the larger the size of each group. The best way to verify normality is to study the residuals of each observation with respect to the mean of the group to which they belong.

  • Although ANOVA is quite robust even when there is a certain lack of normality, if the symmetry is very pronounced and the size of each group is not very large, the non-parametric Kruskal-Wallis H test can be used instead. It is often recommended to stick with ANOVA unless the lack of normality is very extreme.

  • Atypically extreme data can completely invalidate the conclusions of an ANOVA. If extreme residuals are observed, it is necessary to study in detail which observations they belong to, and it is advisable to recalculate the ANOVA without them and compare the results obtained.

• Same variance between groups (homoscedasticity): the variance within the groups must be approximately equal in all of them. This is because in the null hypothesis it is considered that all observations come from the same population, so they have the same mean and also the same variance.

  • This condition is more important the smaller the size of the groups.

  • ANOVA is quite robust to the lack of homoscedasticity if the design is balanced (same number of observations per group).

  • In unbalanced designs, the lack of homoscedasticity has a greater impact. If the smaller groups are the ones with the highest variance, the number of false positives will increase. If, on the other hand, the larger groups have a higher variance, the false negatives will increase.

  • When homoscedasticity cannot be accepted, a heteroscedastic ANOVA can be used, which uses the Welch correction (Welch test).

In the book Handbook of Biological Statistics it is considered highly advisable to use balanced designs. This being the case, they consider ANOVA reliable as long as the number of observations per group is not less than 10 and the standard deviation does not vary more than 3 times between groups. For unbalanced models, they recommend examining homoscedasticity in detail, if the variances of the groups are not very similar, it is better to use Welch's ANOVA.

Multiple comparison of means. POST-HOC contrasts

If an ANOVA is significant, it means that at least two of the compared means are significantly different from each other, but it does not determine which ones. To identify them, the means of all groups must be compared two by two using a t-test or another test that allows comparing 2 groups. This is known as post-hoc analysis. Due to the inflation of the type I error, the more comparisons are made, the greater the probability of finding significant differences. For example, for a significance level $\alpha$ = 0.05, 5 significant comparisons are expected out of every 100 just by chance.

To avoid this problem, the significance level can be adjusted based on the number of comparisons (significance correction). If no correction is made, the possibility of false positives (type I error) increases, but if the corrections are very strict, differences that are actually significant may be considered non-significant (type II error). The need for correction or not, and of what type, must be studied carefully in each case. The main post-hoc comparison methods are described below:

Effect size ${\eta }^2$

The effect size of an ANOVA is the value that allows measuring how much variance in the dependent variable is the result of the influence of the independent qualitative variable, or what is the same, how much the independent variable (factor) affects the dependent variable. In ANOVA, the most used measure of effect size is ${\eta }^2$ which is defined as:

The most used classification levels for the effect size are:

• 0.01 = small

• 0.06 = medium

• 0.14 = large

Power of one-way ANOVA

Power tests allow determining the probability of finding significant differences between the means for a given $\alpha$ indicating the size of the groups, or calculating the size that the groups must have to be able to detect with a certain probability a difference in the means if it exists. In those cases where you want to know the size that the samples must have, it is necessary to know (either from previous experiments or from pilot samples) an estimate of the population variance.

Results of an ANOVA

When communicating the results of an ANOVA, it is convenient to indicate the value obtained for the $F$ statistic, the degrees of freedom, the p-value, and the effect size (${\eta }^2$).

It may happen that an ANOVA analysis indicates that there are significant differences between the means and then the t-tests do not find any comparison that is significant. This is not contradictory, it is simply because they are two different techniques.

Libraries

The libraries used in this example are:

# ==============================================================================
# Libraries
# ==============================================================================
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['font.size'] = 8
import seaborn as sns
import pingouin as pg
from statsmodels.graphics.factorplots import interaction_plot

Data

The following table contains a sample of randomly selected league players.

position = ["OF", "IF", "IF", "OF", "IF", "IF", "OF", "OF", "IF", "IF", "OF",
            "OF", "IF", "OF", "IF", "IF", "IF", "OF", "IF", "OF", "IF", "OF",
            "IF", "OF", "IF", "DH", "IF", "IF", "IF", "OF", "IF", "IF", "IF",
            "IF", "OF", "IF", "OF", "IF", "IF", "IF", "IF", "OF", "OF", "IF",
            "OF", "OF", "IF", "IF", "OF", "OF", "IF", "OF", "OF", "OF", "IF",
            "DH", "OF", "OF", "OF", "IF", "IF", "IF", "IF", "OF", "IF", "IF",
            "OF", "IF", "IF", "IF", "OF", "IF", "IF", "OF", "IF", "IF", "IF",
            "IF", "IF", "IF", "OF", "DH", "OF", "OF", "IF", "IF", "IF", "OF",
            "IF", "OF", "IF", "IF", "IF", "IF", "OF", "OF", "OF", "DH", "OF",
            "IF", "IF", "OF", "OF", "C", "IF", "OF", "OF", "IF", "OF", "IF",
            "IF", "IF", "OF", "C", "OF", "IF", "C", "OF", "IF", "DH", "C", "OF",
            "OF", "IF", "C", "IF", "IF", "IF", "IF", "IF", "IF", "OF", "C", "IF",
            "OF", "OF", "IF", "OF", "IF", "OF", "DH", "C", "IF", "OF", "IF",
            "IF", "OF", "IF", "OF", "IF", "C", "IF", "IF", "OF", "IF", "IF",
            "IF", "OF", "OF", "OF", "IF", "IF", "C", "IF", "C", "C", "OF", "OF",
            "OF", "IF", "OF", "IF", "C", "DH", "DH", "C", "OF", "IF", "OF", "IF",
            "IF", "IF", "C", "IF", "OF", "DH", "IF", "IF", "IF", "OF", "OF", "C",
            "OF", "OF", "IF", "IF", "OF", "OF", "OF", "OF", "OF", "OF", "IF",
            "IF", "DH", "OF", "IF", "IF", "OF", "IF", "IF", "IF", "IF", "OF",
            "IF", "C", "IF", "IF", "C", "IF", "OF", "IF", "DH", "C", "OF", "C",
            "IF", "IF", "OF", "C", "IF", "IF", "IF", "C", "C", "C", "OF", "OF",
            "IF", "IF", "IF", "IF", "OF", "OF", "C", "IF", "IF", "OF", "C", "OF",
            "OF", "OF", "OF", "OF", "OF", "OF", "OF", "OF", "OF", "OF", "C",
            "IF", "DH", "IF", "C", "DH", "C", "IF", "C", "OF", "C", "C", "IF",
            "OF", "IF", "IF", "IF", "IF", "IF", "IF", "IF", "IF", "OF", "OF",
            "OF", "IF", "OF", "OF", "IF", "IF", "IF", "OF", "C", "IF", "IF",
            "IF", "IF", "OF", "OF", "IF", "OF", "IF", "OF", "OF", "OF", "IF",
            "OF", "OF", "IF", "OF", "IF", "C", "IF", "IF", "C", "DH", "OF", "IF",
            "C", "C", "IF", "C", "IF", "OF", "C", "C", "OF"]

batting = [0.359, 0.34, 0.33, 0.341, 0.366, 0.333, 0.37, 0.331, 0.381, 0.332,
         0.365, 0.345, 0.313, 0.325, 0.327, 0.337, 0.336, 0.291, 0.34, 0.31,
         0.365, 0.356, 0.35, 0.39, 0.388, 0.345, 0.27, 0.306, 0.393, 0.331,
         0.365, 0.369, 0.342, 0.329, 0.376, 0.414, 0.327, 0.354, 0.321, 0.37,
         0.313, 0.341, 0.325, 0.312, 0.346, 0.34, 0.401, 0.372, 0.352, 0.354,
         0.341, 0.365, 0.333, 0.378, 0.385, 0.287, 0.303, 0.334, 0.359, 0.352,
         0.321, 0.323, 0.302, 0.349, 0.32, 0.356, 0.34, 0.393, 0.288, 0.339,
         0.388, 0.283, 0.311, 0.401, 0.353, 0.42, 0.393, 0.347, 0.424, 0.378,
         0.346, 0.355, 0.322, 0.341, 0.306, 0.329, 0.271, 0.32, 0.308, 0.322,
         0.388, 0.351, 0.341, 0.31, 0.393, 0.411, 0.323, 0.37, 0.364, 0.321,
         0.351, 0.329, 0.327, 0.402, 0.32, 0.353, 0.319, 0.319, 0.343, 0.288,
         0.32, 0.338, 0.322, 0.303, 0.356, 0.303, 0.351, 0.325, 0.325, 0.361,
         0.375, 0.341, 0.383, 0.328, 0.3, 0.277, 0.359, 0.358, 0.381, 0.324,
         0.293, 0.324, 0.329, 0.294, 0.32, 0.361, 0.347, 0.317, 0.316, 0.342,
         0.368, 0.319, 0.317, 0.302, 0.321, 0.336, 0.347, 0.279, 0.309, 0.358,
         0.318, 0.342, 0.299, 0.332, 0.349, 0.387, 0.335, 0.358, 0.312, 0.307,
         0.28, 0.344, 0.314, 0.24, 0.331, 0.357, 0.346, 0.351, 0.293, 0.308,
         0.374, 0.362, 0.294, 0.314, 0.374, 0.315, 0.324, 0.382, 0.353, 0.305,
         0.338, 0.366, 0.357, 0.326, 0.332, 0.323, 0.306, 0.31, 0.31, 0.333,
         0.34, 0.4, 0.389, 0.308, 0.411, 0.278, 0.326, 0.335, 0.316, 0.371,
         0.314, 0.384, 0.379, 0.32, 0.395, 0.347, 0.307, 0.326, 0.316, 0.341,
         0.308, 0.327, 0.337, 0.36, 0.32, 0.372, 0.306, 0.305, 0.347, 0.281,
         0.281, 0.296, 0.306, 0.343, 0.378, 0.393, 0.337, 0.327, 0.336, 0.32,
         0.381, 0.306, 0.358, 0.311, 0.284, 0.364, 0.315, 0.342, 0.367, 0.307,
         0.351, 0.372, 0.304, 0.296, 0.332, 0.312, 0.437, 0.295, 0.316, 0.298,
         0.302, 0.342, 0.364, 0.304, 0.295, 0.305, 0.359, 0.335, 0.338, 0.341,
         0.3, 0.378, 0.412, 0.273, 0.308, 0.309, 0.263, 0.291, 0.359, 0.352,
         0.262, 0.274, 0.334, 0.343, 0.267, 0.321, 0.3, 0.327, 0.313, 0.316,
         0.337, 0.268, 0.342, 0.292, 0.39, 0.332, 0.315, 0.298, 0.298, 0.331,
         0.361, 0.272, 0.287, 0.34, 0.317, 0.327, 0.354, 0.317, 0.311, 0.174,
         0.302, 0.302, 0.291, 0.29, 0.268, 0.352, 0.341, 0.265, 0.307, 0.36,
         0.305, 0.254, 0.279, 0.321, 0.305, 0.35, 0.308, 0.326, 0.219, 0.23,
         0.322, 0.405, 0.321, 0.291, 0.312, 0.357, 0.324]

         
data = pd.DataFrame({'position': position, 'batting': batting})
data.head(4)

Number of groups, observations per group and distribution of observations

The number of groups and the number of observations per group are identified to determine if it is a balanced model. The mean and standard deviation of each group are also calculated.

# Number of observations per group
# ==============================================================================
data.groupby('position').size()

position
C      39
DH     14
IF    154
OF    120
dtype: int64

# Mean and standard deviation per group
# ==============================================================================
data.groupby('position').agg(['mean', 'std'])

Since the number of observations per group is not constant, it is an unbalanced model. It is important to keep this in mind when checking the conditions of normality and homoscedasticity.

Graphical analysis

Two of the most useful representations before performing an ANOVA are the Box-Plot and Violin-Plot graphs.

fig, ax = plt.subplots(1, 1, figsize=(8, 4))
sns.boxplot(x="position", y="batting", data=data, ax=ax)
sns.swarmplot(x="position", y="batting", data=data, color='black', alpha = 0.5, ax=ax);

This type of representation allows for a preliminary identification of asymmetries, outliers, or differences in variances. In this case, the 4 groups seem to follow a symmetrical distribution. At the IF level, some extreme values are detected that will need to be studied in detail in case they need to be removed. The size of the boxes is similar for all levels, so there are no signs of lack of homoscedasticity.

Verify conditions for an ANOVA

Independence

The groups (categorical variable) and players within each group are independent of each other since a random sampling of players from the entire league (not just from the same team) has been done.

Normal distribution of observations

The quantitative variable must be normally distributed in each of the groups. The study of normality can be done graphically (qqplot) or with hypothesis tests.

# qqplot graphs
# ==============================================================================
fig, axs = plt.subplots(2, 2, figsize=(8, 7))
pg.qqplot(data.loc[data.position=='OF', 'batting'], dist='norm', ax=axs[0,0])
axs[0,0].set_title('OF')
pg.qqplot(data.loc[data.position=='IF', 'batting'], dist='norm', ax=axs[0,1])
axs[0,1].set_title('IF')
pg.qqplot(data.loc[data.position=='DH', 'batting'], dist='norm', ax=axs[1,0])
axs[1,0].set_title('DH')
pg.qqplot(data.loc[data.position=='C', 'batting'], dist='norm', ax=axs[1,1])
axs[1,1].set_title('C')
plt.tight_layout()

# Shapiro-Wilk normality test
# ==============================================================================
pg.normality(data=data, dv='batting', group='position')

Neither the graphical analysis nor the hypothesis tests show evidence of lack of normality.

Constant variance between groups (homoscedasticity)

Since there is an IF group that is on the limit to accept that it is normally distributed, the Levene test is more appropriate than the Bartlett test.

# Homoscedasticity test
# ==============================================================================
pg.homoscedasticity(data=data, dv='batting', group='position', method='levene')

According to the Levene test, there is no significant evidence to indicate a lack of homoscedasticity.

ANOVA Test

# One-way ANOVA test
# ==============================================================================
pg.anova(data=data, dv='batting', between='position', detailed=True)

A p-value greater than 0.1 is very weak evidence against the null hypothesis that all groups have the same mean. The eta-squared value (${\eta }^2$) is 0.018, which can be considered a small effect size.

Post-hoc multiple comparisons

In this case, the ANOVA was not significant, so it makes no sense to perform pairwise comparisons. However, for educational purposes, it is shown how to obtain the TukeyHSD tests.

# Post-hoc Tukey test
# ==============================================================================
pg.pairwise_tukey(data=data, dv='batting', between='position').round(3)

As expected, no significant difference is found between any pair of means.

Conclusion

In the study carried out, a small effect size was observed and the ANOVA inference techniques have not found statistical significance to reject that the means are equal among all groups.

Definition

The two-way analysis of variance, also known as a factorial plan with two factors, is used to study the relationship between a quantitative dependent variable and two qualitative independent variables (factors), each with several levels. The two-way ANOVA allows studying how each of the factors influences the dependent variable by themselves (additive model) as well as the influence of the combinations that can occur between them (model with interaction).

Suppose we want to study the effect of a drug on blood pressure (quantitative dependent variable) depending on the patient's sex (levels: male, female) and age group (levels: child, adult, elderly).

The simple effect of the factors consists of studying how the effect of the drug varies depending on sex without differentiating by age, as well as studying how the effect of the drug varies depending on age without taking sex into account. The effect of the double interaction consists of studying whether the influence of one of the factors varies depending on the levels of the other factor. That is, if the influence of the sex factor on the activity of the drug is different depending on the age of the patient or, what is the same, if the activity of the drug for a certain age is different depending on whether you are a man or a woman.

Conditions for two-way ANOVA for independent data

The necessary conditions for a two-way ANOVA to be valid, as well as the process to follow to perform it, are similar to the one-way ANOVA. The only differences are:

• Hypothesis: the two-way ANOVA with repetitions combines 3 null hypotheses: that the means of the observations grouped by one factor are equal, that the means of the observations grouped by the other factor are equal, and that there is no interaction between the two factors.

• It requires calculating the Sum of Squares and Mean Squares for both factors. It is common to find that one factor is called "treatment" and the other block.

It is important to keep in mind that the order in which the factors are multiplied does not affect the result of the ANOVA only if the size of the groups is equal (balanced model), otherwise it does matter. For this reason, it is recommended that the design be balanced.

In addition, the study of the interaction of the two factors is only possible if there are several observations for each of the combinations of the levels.

Effect size

In the case of ANOVA with two factors, the eta-squared (${\eta }^2$) effect size can be calculated for each of the two factors as well as for the interaction.

Libraries

The libraries used in this example are:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pingouin as pg

Data

strength = [15.29, 15.89, 16.02, 16.56, 15.46, 16.91, 16.99, 17.27, 16.85,
               16.35, 17.23, 17.81, 17.74, 18.02, 18.37, 12.07, 12.42, 12.73,
               13.02, 12.05, 12.92, 13.01, 12.21, 13.49, 14.01, 13.30, 12.82,
               12.49, 13.55, 14.53]

tempering = ["fast"] * 15 + ["slow"] * 15

thickness = ([8] * 5 + [16] * 5 + [24] * 5) * 2

data = pd.DataFrame({
    'tempering': tempering,
    'thickness': thickness,
    'strength': strength
})

data.head()

Descriptive and graphical analysis

First, Box-plot diagrams are generated to identify possible notable differences, asymmetries, outliers, and homogeneity of variance between the different levels. The mean and variance of each group are also calculated.

fig, axs = plt.subplots(1, 2, figsize=(10, 4))
axs[0].set_title('Strength vs tempering')
sns.boxplot(x="tempering", y="strength", data=data, ax=axs[0])
sns.swarmplot(x="tempering", y="strength", data=data, color='black',
              alpha = 0.5, ax=axs[0])

axs[1].set_title('Strength vs thickness')
sns.boxplot(x="thickness", y="strength", data=data, ax=axs[1])
sns.swarmplot(x="thickness", y="strength", data=data, color='black',
              alpha = 0.5, ax=axs[1]);

fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.set_title('Strength vs tempering and thickness')
sns.boxplot(x="tempering", y="strength", hue='thickness', data=data, ax=ax);

print('Mean strength and standard deviation by tempering')
data.groupby('tempering')['strength'].agg(['mean', 'std'])

Mean strength and standard deviation by tempering

print('Mean strength and standard deviation by thickness')
data.groupby('thickness')['strength'].agg(['mean', 'std'])

Mean strength and standard deviation by thickness

print('Mean strength and standard deviation by tempering and thickness')
data.groupby(['tempering', 'thickness'])['strength'].agg(['mean', 'std'])

Mean strength and standard deviation by tempering and thickness

From the graphical representation and the calculation of the means, it can be inferred that there is a difference in the strength achieved depending on the type of tempering. The strength seems to increase as the thickness of the sheet increases, although it is not clear that the difference in the means is significant. The distribution of the observations of each level seems symmetrical without the presence of outliers. A priori, it seems that the necessary conditions for an ANOVA are met.

Another way to graphically identify possible interactions between the two factors is through what are known as "interaction plots". If the lines that describe the data for each of the levels are parallel, it means that the behavior is similar regardless of the level of the factor, that is, there is no interaction.

fig, ax = plt.subplots(figsize=(6, 4))
fig = interaction_plot(
    x        = data.tempering,
    trace    = data.thickness,
    response = data.strength,
    ax       = ax,
)

fig, ax = plt.subplots(figsize=(6, 4))
fig = interaction_plot(
    x        = data.thickness,
    trace    = data.tempering,
    response = data.strength,
    ax       = ax,
)

The first interaction plot seems to indicate that the increase in strength between the two types of tempering is proportional for the three thicknesses. In the second plot, a certain deviation is observed in the 24mm thickness. This could be due to simple variability or because there is an interaction between the thickness and tempering variables. These indications will be confirmed or discarded by the ANOVA.

ANOVA Test

# Two-way ANOVA test
# ==============================================================================
pg.anova(
    data     = data,
    dv       = 'strength',
    between  = ['tempering', 'thickness'],
    detailed = True
).round(4)

The analysis of variance confirms that there is a significant influence on the strength of the sheets by both factors (tempering and thickness) with large and medium ${\eta }^2$ effect sizes, respectively. However, no significant interaction between them is detected.

Libraries

The libraries used in this example are:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pingouin as pg

Data

subject = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
           21, 22, 23, 24, 25, 26, 27, 28, 29, 30]

sex = ["female", "male", "male", "female", "male", "male", "male", "female",
       "female", "male", "male", "male", "male", "female", "female", "female",
       "male", "female", "female", "male", "male", "female", "male", "male",
       "male", "male", "male", "male", "female", "male"]

age = ["adult", "adult", "adult", "adult", "adult", "adult", "young", "young",
       "adult", "young", "young", "adult", "young", "young", "young", "adult",
       "young", "adult", "young", "young", "young", "young", "adult", "young",
       "young", "young", "young", "young", "young", "adult"]

result = [7.1, 11.0, 5.8, 8.8, 8.6, 8.0, 3.0, 5.2, 3.4, 4.0, 5.3, 11.3, 4.6, 6.4,
          13.5, 4.7, 5.1, 7.3, 9.5, 5.4, 3.7, 6.2, 10.0, 1.7, 2.9, 3.2, 4.7, 4.9,
          9.8, 9.4]

data = pd.DataFrame({
    'subject': subject,
    'sex': sex,
    'age': age,
    'result': result
})

data.head()

Descriptive and graphical analysis

First, Box-plot diagrams are generated to identify possible notable differences, asymmetries, outliers, and homogeneity of variance between the different levels. The mean and variance of each group are also calculated.

fig, axs = plt.subplots(1, 2, figsize=(10, 4))
axs[0].set_title('Results vs age')
sns.boxplot(x="age", y="result", data=data, ax=axs[0])
sns.swarmplot(x="age", y="result", data=data, color='black',
              alpha = 0.5, ax=axs[0])

axs[1].set_title('Results vs sex')
sns.boxplot(x="sex", y="result", data=data, ax=axs[1])
sns.swarmplot(x="sex", y="result", data=data, color='black',
              alpha = 0.5, ax=axs[1]);

fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.set_title('Results vs sex and age')
sns.boxplot(x="age", y="result", hue='sex', data=data, ax=ax);

print('Mean results and standard deviation by age')
data.groupby('age')['result'].agg(['mean', 'std'])

Mean results and standard deviation by age

print('Mean results and standard deviation by sex')
data.groupby('sex')['result'].agg(['mean', 'std'])

Mean results and standard deviation by sex

print('Mean results and standard deviation by age and sex')
data.groupby(['age', 'sex'])['result'].agg(['mean', 'std'])

Mean results and standard deviation by age and sex

From the graphical representation and the calculation of the means, it can be inferred that there is a difference in the effect of the drug depending on age and also on sex. The effect seems to be greater in women than in men and in adults than in young people, although the significance will have to be confirmed with the ANOVA. The distribution of the observations of each level seems symmetrical with the presence of a single outlier. A priori, it seems that the necessary conditions for an ANOVA are met, although they will have to be confirmed by studying the residuals.

# Interaction plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 4))
fig = interaction_plot(
    x        = data.age,
    trace    = data.sex,
    response = data.result,
    ax       = ax,
)

A clear interaction between both factors is observed. The response to the drug is different between adults and young people, and has an inverse trend depending on sex. In women, the response is greater when they are young than when they are adults, and in men, it is greater when they are adults and lower when they are young. The ANOVA will confirm if the observed differences are significant.

ANOVA Test

# Two-way ANOVA test
# ==============================================================================
pg.anova(
    data     = data,
    dv       = 'result',
    between  = ['sex', 'age'],
    detailed = True
).round(4)

The analysis of variance finds no significant differences in the effect of the drug between men and women (sex factor) but does find significant differences between young and old and between at least two groups of the combinations of sex and age, that is, there is significance for the interaction. The effect size ${\eta }^2$ is large for both age and the interaction of age and sex.

Note: in this case, the order in which the factors are multiplied does affect the results since the size of the groups is not equal.

Conditions for ANOVA with dependent variables

Sphericity

This is the only condition to be able to apply this type of analysis. Sphericity implies that the variance of the differences between all pairs of variables to be compared is equal. If the ANOVA is performed to compare the mean of a quantitative variable between three levels (A, B, C), sphericity will be accepted if the variance of the differences between A-B, A-C, B-C is the same.

Sphericity can be analyzed using Mauchly's test, whose null hypothesis is that sphericity exists. In practice, it is common for the sphericity condition not to be met, in which case, the Greenhouse-Geisser and Huynh-Feldt corrections can be applied to the p-values. Another alternative, if sphericity is not met, is the non-parametric Friedman test.

The rm_anova function of the pingouin library, if the correction = True argument is indicated, automatically applies Mauchly's sphericity test, if it finds evidence that the condition is not met, it corrects the p-values with the Greenhouse-Geisser method.

Libraries

The libraries used in this example are:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pingouin as pg
from statsmodels.graphics.factorplots import interaction_plot

Data

product = ["lettuce", "potatoes", "milk", "eggs", "bread", "cereal", "ground.beef",
            "tomato.soup", "laundry.detergent", "aspirin"]
store_A = [1.755, 2.655, 2.235, 0.975, 2.370, 4.695, 3.135, 0.930, 8.235, 6.690]
store_B = [1.78, 1.98, 1.69, 0.99, 1.70, 3.15, 1.88, 0.65, 5.99, 4.84]
store_C = [1.29, 1.99, 1.79, 0.69, 1.89, 2.99, 2.09, 0.65, 5.99, 4.99]
store_D = [1.29, 1.99, 1.59, 1.09, 1.89, 3.09, 2.49, 0.69, 6.99, 5.15]

data = pd.DataFrame({
    'product':product * 4,
    'store': np.repeat(['A', 'B', 'C', 'D'], 10),
    'price': store_A + store_B + store_C + store_D
})

data.head()

Descriptive and graphical analysis

fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.set_title('Prices vs store')
sns.boxplot(x="store", y="price", data=data, ax=ax)
sns.swarmplot(x="store", y="price", data=data, color='black', alpha=0.5, ax=ax);

print('Mean prices and standard deviation by store')
data.groupby('store')['price'].agg(['mean', 'std'])

Mean prices and standard deviation by store

ANOVA Test

# Repeated measures ANOVA test
# ==============================================================================
pg.rm_anova(
    data       = data,
    dv         = 'price',
    within     = 'store',
    subject    = 'product',
    detailed   = True,
    correction = 'auto'
).round(4)

Since sphericity is not met (sphericity = False), the corrected p-value must be used, which is in the p-GG-corr column. The analysis of variance shows significant evidence with a large effect size.

Multiple comparisons

Since the data are paired, Tukey HSD intervals cannot be calculated. Instead, paired t-tests can be used.

# Post-hoc pairwise t-test
# ==============================================================================
pg.pairwise_tests(
    dv      = 'price',
    within  = 'store',
    subject = 'product',
    padjust = 'holm',
    data    = data
)

Conclusion

The ANOVA analysis (including corrections given the lack of sphericity) finds significant differences in the price of food between at least 2 stores, with a large effect size. The subsequent pairwise comparison by t-student with holm significance correction identifies the differences between stores A-B, A-C, and A-D as significant, but not between B-C, B-D, and C-D.