Analysis of variance, ANOVA (statistical analysis of variance) – a statistical method used to study observations that depend on one or many factors acting simultaneously. This method explains with what probability the extracted factors may be the reason for differences between the observed group means. Analysis of variance was created in the 1920s by Ronald Fisher.
Variance analysis models can be divided into:
- one-factor models – the impact of each factor is considered separately, this class of issues is dealt with by one-way analysis of variance
- multifactorial models – the impact of various factors is considered together, this class of issues is dealt with by multifactorial analysis of variance.
Source: https://pl.wikipedia.org/wiki/Analiza_wariancji
We check if the insulin level in the blood of patients differs depending on age.
We analyze only one factor because we are interested in the quality of treatment (this is one factor) so we choose one-factor model of ANOVA.
Data source: https://www.kaggle.com/saurabh00007/diabetescsv
import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv('c:/1/diabetes.csv', usecols=['Insulin','Age'])
df.head(3)
We delete the missing data. The level of insulin in the blood cannot be zero, these records are removed as unsuitable for our tests.
import numpy as np
df['Insulin']=df['Insulin'].replace(0,np.nan)
df = df.dropna(how ='any')
df.isnull().sum()
We divide patients into four groups depending on their age.
df['Age_G']=pd.qcut(df['Age'], 4,labels=["young","adults","middle-aged","old"])
df.head()
Test ANOVA
import statsmodels.api as sm
from statsmodels.formula.api import ols
model2 = ols('Insulin ~ C(Age_G)', data=df).fit()
anova_table = sm.stats.anova_lm(model2, typ=2)
print(anova_table)
Interpretation: The P value obtained from ANOVA is significant (P <0.05), so we conclude that there is a statistically significant difference between blood insulin levels between the four age groups listed.
ANOVA have shown that there is a statistically significant difference in insulin levels between the age groups of patients. Unfortunately, ANOVA does not indicate which age groups are significantly different from each other.
To find out age group pairs with significantly different blood insulin levels, multiple pairs comparison analysis (post-hoc comparison) should be performed using Tukey’s HSD test.
Post hoc tests
Post-hoc tests should only be performed if p ANOVA is significant (i.e. greater than 0.05).
Post hoc (post hoc comparisons) – statistical tests performed after obtaining a significant F value after analyzing variance. They show which means differ statistically significantly. If the analysis of variance shows significant differences between the groups under consideration, further tests should be carried out to answer the question of which of the compared populations are responsible for rejecting the null hypothesis. We want to know which of the n averages are different and which are equal. For this purpose we use special post-hoc tests (after the fact), also called multiple comparison tests. The name of these tests results from the fact that we use them only after finding (by analysis of variance) the lack of equality between the means. These tests are also sometimes called homogeneous grouping tests, because after their application we can get medium groups. means belonging to the same group do not differ significantly, while means belonging to different groups will differ significantly (http://home.agh.edu.pl/~bartus/index.php?action=dydaktyka&subaction=statystyka&item=testy_post_hoc)
Tukey HSD test
from statsmodels.stats.multicomp import pairwise_tukeyhsd
m_comp = pairwise_tukeyhsd(endog=df['Insulin'], groups=df['Age_G'], alpha=0.05)
print(m_comp)
The Tukey HSD test does not give us information on average location and variance. Fortunately, we have other ways to suspect this. The position of the red vertical line is rather indicative. This confirms the results of the Tukey HSD test.
m_comp.plot_simultaneous() # Plot group confidence intervals
plt.vlines(x=160,ymin=-0.5,ymax=3.5, color="red")
The Tukey HSD test showed a statistically significant difference in blood insulin content between the “old” and “young” flus as well as “adults” and “old”
The above results from Tukey HSD suggest that, apart from the old-young and adults-old comparison, all other pair comparisons retain the null hypothesis, i.e. in the other groups there are no significant statistical differences in blood insulin content.
To better understand mutual similarities of variance (or their absence) we can perform a t test. Inside the matrix there are values of p-value coefficients.
T-test
Pair T-test for many comparisons of independent groups. It can be used after parametric ANOVA for pairwise comparisons.
import scikit_posthocs as sp
sp.posthoc_ttest(df, val_col='Insulin', group_col='Age_G', p_adjust='holm')
ttest = sp.posthoc_ttest(df, val_col='Insulin', group_col='Age_G', p_adjust='holm')
ttest
As you can easily see, there is a significant difference between the young -old and adults -old groups (high p-value> 0.05).
where ones indicate which pairs do not have the same average. P-values are replaced with asterisks: – p < 0.05, – p < 0.01, – p < 0.001..
sp.sign_table(ttest)
Checking compliance with ANOVA¶
Conditions:
- residual values have a normal distribution (Wilks Shapiro test)
- variances in groups are homogeneous (Levene or Bartlett test)
- observations are carried out independently of each other
Levene’s test Checking the homogeneity of variance
Null hypothesis: a group from the population have equal variances.
df.head()
DFF = df.reset_index()
KOT = pd.pivot_table(DFF,index='index',columns='Age_G', values='Insulin')
young=KOT['young'].dropna(how='any')
old=KOT['old'].dropna(how='any')
middleaged=KOT['middle-aged'].dropna(how='any')
adults=KOT['adults'].dropna(how='any')
import scipy.stats as stats
w,p = stats.levene(young,old,middleaged,adults)
print("Value: ",w)
print("p-value: ",p)
Because the P value is not significant (p> 0.05), we do not reject the null hypothesis – i.e. age groups have uniform variances.
Shapiro-Wilk test Checking normality of residue distribution
Null hypothesis: residual values are normally distributed.
import scipy.stats as stats
w, p = stats.shapiro(model2.resid)
print("Value: ",w)
print("p-value: ",p)
Because the P value is significant (p< 0.05), we reject the null hypothesis – that is, the residual values have a normal distribution.
At the end we can make a chart comparing age groups:
import matplotlib.pyplot as plt
KOT.boxplot(column=['young', 'old', 'middle-aged', 'adults'], grid=False)
Analysis of air pollution by PT08.S1 (CO) using ANOVA
- H0: Air pollution PT08.S1 (CO) does NOT differ significantly over four daily periods
- H1: Air pollution PT08.S1 (CO) significantly differs over four daily periods
We analyze only one factor because we are interested in the quality of treatment (this is one factor) so we choose one-factor model of ANOVA.
data source:
import pandas as pd
df2 = pd.read_csv('c:/TF/AirQ_filled.csv')
df2.head(3)
df2.dtypes
We process the Time column into a number format.
df2['Time'] = df2.Time.str.slice(0,2)
df2['Time']= df2['Time'].convert_objects(convert_numeric=True)
df2['Time'].dtype
We check the completeness of the data.
df2.isnull().sum()
I separate four daily periods.
df2['Periods'] = pd.qcut(df2.Time,4,labels=["0-5","6-12","12-18","18-24"])
df2.dtypes
I check daily times.
pd.pivot_table(df2, index='Periods', values='Time', aggfunc=['min', 'max'])
Now we create a pivot table, where the columns will have four daily periods.
PKS = pd.pivot_table(df2, index = 'Date', columns = 'Periods', values='PT08.S1(CO)')
PKS.head(4)
We accept research hypotheses for pollutants with the substance PT08.S1 (CO)
- H0: Air pollution PT08.S1 (CO) does NOT differ significantly over four daily periods
- H1: Air pollution PT08.S1 (CO) significantly differs over four daily periods
Test ANOVA
I change the name PT08.S1(CO) to PT08S1CO.
df2.rename(columns={'PT08.S1(CO)':'PT08S1CO'},inplace=True)
import statsmodels.api as sm
from statsmodels.formula.api import ols
model3 = ols('PT08S1CO ~ C(Periods)', data=df2).fit()
anova_table = sm.stats.anova_lm(model3, typ=2)
print(anova_table)
Interpretation: The P value obtained from ANOVA is significant (P <0.05), so we conclude that there are statistically significant differences in air pollution PT08.S1 (CO) over four daily periods.
ANOVA has shown that there is a statistically significant difference in the level of contamination and the null hypothesis has been rejected in favor of the alternative hypothesis. Unfortunately, ANOVA does not indicate which periods (possibly) differ from each other.
To check and detect statistically significant differences, multiple pairs comparison analysis (post-hoc comparison) should be performed using Tukey’s HSD test.
Tukey HSD test
from statsmodels.stats.multicomp import pairwise_tukeyhsd
m_comp = pairwise_tukeyhsd(endog=df2['PT08S1CO'], groups=df2['Periods'], alpha=0.05)
print(m_comp)
Turkey’s HDS test showed that only a comparison of groups 6-12 from 12-18 did not give statistically significant differences in PT08.S1 (CO) contamination.
Checking compliance with ANOVA
Conditions:
1.residual values have a normal distribution (Wilks Shapiro test)
- variances in groups are homogeneous (Levene or Bartlett test)
- observations are carried out independently of each other
Levene’s test Checking variance homogeneity
Null hypothesis: a group from the population have equal variances.
Now we create a pivot table, where there will be 4 daily periods in columns.
PKS = pd.pivot_table(df2, index = 'Date', columns = 'Periods', values='PT08S1CO')
PKS.head(4)
P01=PKS['0-5'].dropna(how='any')
P02=PKS['6-12'].dropna(how='any')
P03=PKS['12-18'].dropna(how='any')
P04=PKS['18-24'].dropna(how='any')
import scipy.stats as stats
w,p = stats.levene(P01,P02,P03,P04)
print("Value: ",w)
print("p-value: ",p)
Because the P value is significant (p< 0.05), we reject the null hypothesis – i.e. age groups have not uniform variances.
Shapiro-Wilk test Checking the normality of residue distribution
Null hypothesis: residual values are normally distributed.
import scipy.stats as stats
w, p = stats.shapiro(model3.resid)
print("Value: ",w)
print("p-value: ",np.round(p, decimals=2))
Because the P value of ZERO is significant because it is smaller than the confidence factor 0.05 (p> 0.05), we reject the null hypothesis – the residual values do not have a normal distribution.
We’ll take a closer look at the rest of the model. Let’s see how it looks in the chart.
x=model3.resid
title = "Residuals"
x_label = "level"
y_label = "probability"
def Dist1(x, ax, title, x_label, y_label):
x.plot.kde(ax=ax, legend=False)
ax.set_title(title, color='darkred', alpha=1)
ax.set_ylabel(y_label, color='grey', alpha=0.6)
ax.set_xlabel(x_label, color='grey', alpha=0.6)
fig, ax = plt.subplots(figsize=(6, 3))
Dist1(x, ax, title, x_label, y_label)
import scipy
scipy.stats.anderson(model3.resid, dist='norm')
This time check the normal distribution of residues using the QQ chart.
Q-Q plot
from statsmodels.graphics.gofplots import qqplot
from matplotlib import pyplot
qqplot(model3.resid, line='s')
pyplot.show()
The Shapiro-Wilk test again
made from scipy library.
from scipy.stats import shapiro
stat, p = shapiro(model3.resid)
print('Statistics=
K ^ 2 D’Agostino tests
The K ^ 2 D’Agostino test (Ralpha D’Agostino) calculates summary statistics from the data, namely kurtosis and skewness, to determine whether the distribution of the data deviates from the normal distribution. This is a simple and commonly used statistical test of normality.
– Skew (skew) is a quantification of how far the distribution is shifted left or right, a measure of asymmetry in the distribution.
– Kurtosis quantifies the distribution of the tail.
from scipy.stats import normaltest
stat, p = normaltest(model3.resid)
print('Statistics=
The Shapiro-Wilk test and D’Agostino’s K^2 test indicated a disturbance in the normality of the probability distribution of residual values.
I can’t comment on that.
It turned out that the distribution of residual values from the OLS ANOVA model does not have a normal distribution, which means that we cannot use parametric ANOVA. So we have to use nonparametric methods to reduce variance. I will deal with this in the next post.
Finally a chart comparing the groups.
import matplotlib.pyplot as plt
PKS.boxplot(column=['0-5', '6-12', '12-18', '18-24'], grid=False)
