EN030120201045
In the last post we did the ANOVA analysis of the variance comparison in groups. It turned out, unexpectedly, that one of the conditions of the correct ANOVA ols model was not met.
ANOVA the one-factor models compatibility of variance 020120200957
The ANOVA model has three basic 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.
ANOVA ols did not pass the first Wilks Shapiro test for normal distribution of residual values and Levene test.
Therefore, ANOVA cannot be used.
Now I will repeat the whole procedure from the previous entry. I will then perform an analysis of variance using a nonparametric test to replace ANOVA.
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.
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 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 numpy as np
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)
import matplotlib.pyplot as plt
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 non normality of the probability distribution of residual values.
If ANOVA does not meet the condition for the normal distribution of residual values, a nonparametric test should be performed.
The Kruskal-Wallis H-test
The Kruskal-Wallis H test (called nonparametric ANOVA) tests the null hypothesis that the median population of all groups is equal. This is a non-parametric version of ANOVA. The test works on 2 or more independent samples that can be of different sizes. The test does not assume distribution normality. Sometimes it is considered as a nonparametric alternative to one-way analysis of variance between groups. Note that rejecting the null hypothesis does not indicate which group is different. Post hoc comparisons between groups are required to determine which groups are different. The method accepts array structures, but not DataFrames.
The null hypothesis – the equality of distribution functions in compared populations.
PKS.head(3)
We delete empty values from the data, otherwise the test will not work.
PKS = PKS.dropna(how='any')
PKS.isnull().sum()
import scipy.stats as ss
H, p = ss.kruskal(PKS['0-5'], PKS['6-12'], PKS['12-18'], PKS['18-24'])
print('p-value: ',p)
Interpretation: p-value is significant (p <0.05) which means the need to reject the null hypothesis that the median population of all groups is equal.
The median is not equal!
To find out what time periods (data groups) differ in terms of the median, we need to perform post hoc tests. Conover test
Conover’s test
In statistics, the Conover quadratic rank test is a nonparametric version of the parametric Levene test for equality of variances. The Conover rank squares test is the only test of equality of variances that seems not to be parametric.
import scikit_posthocs as sp
FF = PKS.columns
x = [PKS['0-5'], PKS['6-12'], PKS['12-18'], PKS['18-24']]
CT = sp.posthoc_conover(x, p_adjust = 'holm')
CT
sp.sign_table(CT)
P-values are replaced with asterisks: – p < 0.05, – p < 0.01, – p < 0.001..
pc = sp.posthoc_conover(x, val_col='values', group_col='groups')
heatmap_args = {'linewidths': 0.25, 'linecolor': '0.5', 'clip_on': False, 'square': True, 'cbar_ax_bbox': [0.80, 0.35, 0.04, 0.3]}
sp.sign_plot(pc, **heatmap_args)
From what could be determined only at night not from 0 to 5 am no contamination with PT08.S1 (CO) is different (smaller). This can also be seen on the graph.
import matplotlib.pyplot as plt
PKS.boxplot(column=['0-5', '6-12', '12-18', '18-24'], grid=False)
