Introduction ¶
- We will briefly talk in this lesson about the Scientific Library for Python (
SciPy), which is the scientific toolbox for Python. - We will get a brief overview of the statistics subpackage and we will use it to perform many statistical calculations, including calculations of probabilities, probability distributions, and confidence intervals.
- We will also perform hypothesis testing on a real-world dataset, and we will be able to state conclusions that go beyond the sample that we have.
In this lesson we will show how to use SciPy Library to perform
- Statistics
- Probabilities
- Hypothesis testing
- Performing statistical tests
Introduction to SciPy ¶
SciPy
- is a tool for doing scientific computing in Python.
- is a Python-based ecosystem that is an open source software for math, science, and engineering.
- contains various toolboxes dedicated to common issues in scientific computing.
Let us first import the libraries
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from numpy import sqrt, arange
from scipy import stats
%matplotlib inline
In this lesson we use a dataset containing information about Portuguese students from two public schools. This is a real world dataset that was collected in order to study alcohol consumption in young people and its effects on students' academic performance. The dataset was built from two sources: school reports and questionnaires.
Attributes:
- 1 school - student's school (binary: 'GP' - Gabriel Pereira or 'MS' - Mousinho da Silveira)
- 2 sex - student's sex (binary: 'F' - female or 'M' - male)
- 3 age - student's age (numeric: from 15 to 22)
- 4 address - student's home address type (binary: 'U' - urban or 'R' - rural)
- 5 famsize - family size (binary: 'LE3' - less or equal to 3 or 'GT3' - greater than 3)
- 6 Pstatus - parent's cohabitation status (binary: 'T' - living together or 'A' - apart)
- 7 Medu - mother's education (numeric: 0 - none, 1 - primary education (4th grade), 2 - 5th to 9th grade, 3 -secondary education or 4 - higher education)
- 8 Fedu - father's education (numeric: 0 - none, 1 - primary education (4th grade), 2 - 5th to 9th grade, 3 -secondary education or 4 - higher education)
- 9 Mjob - mother's job (nominal: 'teacher', 'health' care related, civil 'services' (e.g. administrative or police), 'at_home' or 'other')
- 10 Fjob - father's job (nominal: 'teacher', 'health' care related, civil 'services' (e.g. administrative or police), 'at_home' or 'other')
- 11 reason - reason to choose this school (nominal: close to 'home', school 'reputation', 'course' preference or 'other')
- 12 guardian - student's guardian (nominal: 'mother', 'father' or 'other')
- 13 traveltime - home to school travel time (numeric: 1 <15 min., 2 - 15 to 30 min., 3 - 30 min. to 1 hour, or 4 - >1 hour)
- 14 studytime - weekly study time (numeric: 1 - <2 hours, 2 - 2 to 5 hours, 3 - 5 to 10 hours, or 4 - >10 hours)
- 15 failures - number of past class failures (numeric: n if 1<=n<3, else 4)
- 16 schoolsup - extra educational support (binary: yes or no)
- 17 famsup - family educational support (binary: yes or no)
- 18 paid - extra paid classes within the course subject (Math or Portuguese) (binary: yes or no)
- 19 activities - extra-curricular activities (binary: yes or no)
- 20 nursery - attended nursery school (binary: yes or no)
- 21 higher - wants to take higher education (binary: yes or no)
- 22 internet - Internet access at home (binary: yes or no)
- 23 romantic - with a romantic relationship (binary: yes or no)
- 24 famrel - quality of family relationships (numeric: from 1 - very bad to 5 - excellent)
- 25 freetime - free time after school (numeric: from 1 - very low to 5 - very high)
- 26 goout - going out with friends (numeric: from 1 - very low to 5 - very high)
- 27 Dalc - workday alcohol consumption (numeric: from 1 - very low to 5 - very high)
- 28 Walc - weekend alcohol consumption (numeric: from 1 - very low to 5 - very high)
- 29 health - current health status (numeric: from 1 - very bad to 5 - very good)
- 30 absences - number of school absences (numeric: from 0 to 93)
these grades are related with the course subject:
- 31 G1 - first period grade (numeric: from 0 to 20)
- 31 G2 - second period grade (numeric: from 0 to 20)
- 32 G3 - final grade (numeric: from 0 to 20, output target)
Importing the data into Python ¶
student = pd.read_csv("/Users/dhafermalouche/Documents/Teaching_2018_2019/Teaching_python/student.csv", sep=";")
A preview of the data
student.head()
Confidence Intervals ¶
G3 is the variable that gives the final grade (numeric: from 0 to 20, output target). We will show how to compute the confidence interval of the mean.
The sample size
sample_size = student.shape[0]
sample_size
Because we have a sample size that is much greater than 30, we can use the Central Limit Theorem to calculate confidence intervals. According to this theorem we can calculate a confidence interval for the mean using the normal distribution.
To get the confidence interval for the mean we need three numbers:
- Sample mean
- Standard error
- Confidence level
Formula for the standard error:
$$ SE = \frac{s}{\sqrt n} $$
sample_mean_grade = student['G3'].mean()
sample_mean_grade
std_error_grades = student['G3'].std()/sqrt(sample_size)
std_error_grades
stats.norm.interval(0.95, loc=sample_mean_grade, scale=std_error_grades)
Let us compare the confidence intervals of G3 between girls and boys
Girls data
SF = student.loc[(student.sex == 'F')]
SF.head()
Boys data
SM = student.loc[(student.sex == 'M')]
SM.head()
sample_mean_grade_F =SF['G3'].mean()
sample_mean_grade_F
sample_mean_grade_M =SM['G3'].mean()
sample_mean_grade_M
sample_size_F = SF.shape[0]
sample_size_F
sample_size_M = SM.shape[0]
sample_size_M
std_error_grades_F = SF['G3'].std()/sqrt(sample_size_F)
std_error_grades_F
std_error_grades_M = SM['G3'].std()/sqrt(sample_size_M)
std_error_grades_M
Girls Confidence Interval
stats.norm.interval(0.95, loc=sample_mean_grade_F, scale=std_error_grades_F)
Boys Confidence Interval
stats.norm.interval(0.95, loc=sample_mean_grade_M, scale=std_error_grades_M)
Let us show this graphicaly then
fig, axes = plt.subplots(1,2, figsize = (14,4))
sns.boxplot(x='sex', y='G3', data=student, ax=axes[0])
sns.pointplot(x='sex', y='G3', data=student, ax=axes[1]);
Let us try to answer to the question: Does alcohol consumption affect academic performance?
We can then draw the same graph as above by comparing the final grades between the five categories of drinker in the weekday. We will be using then Dalc variable
fig, axes = plt.subplots(1,2, figsize = (14,4))
sns.boxplot(x='Dalc', y='G3', data=student, ax=axes[0])
sns.pointplot(x='Dalc', y='G3', data=student, ax=axes[1]);
Hypothesis testing ¶
Some Statistical test from scipy.stats
- kurtosistest(a[, axis, nan_policy]) Tests whether a dataset has normal kurtosis
- normaltest(a[, axis, nan_policy]) Tests whether a sample differs from a normal distribution.
- skewtest(a[, axis, nan_policy]) Tests whether the skew is different from the normal distribution.
- pearsonr(x, y) Calculates a Pearson correlation coefficient and the p-value for testing non-correlation.
- ttest_1samp(a, popmean[, axis, nan_policy]) Calculates the T-test for the mean of ONE group of scores.
- ttest_1samp(a, popmean[, axis, nan_policy]) Calculates the T-test for the mean of ONE group of scores.
- ttest_ind(a, b[, axis, equal_var, nan_policy]) Calculates the T-test for the means of two independent samples of scores.
- ttest_ind_from_stats(mean1, std1, nobs1, ...) T-test for means of two independent samples from descriptive statistics.
- ttest_rel(a, b[, axis, nan_policy]) Calculates the T-test on TWO RELATED samples of scores, a and b.
- kstest(rvs, cdf[, args, N, alternative, mode]) Perform the Kolmogorov-Smirnov test for goodness of fit.
- chisquare(f_obs[, f_exp, ddof, axis]) Calculates a one-way chi square test.
- ansari(x, y) Perform the Ansari-Bradley test for equal scale parameters
- bartlett(*args) Perform Bartlett’s test for equal variances
- levene(*args, **kwds) Perform Levene test for equal variances.
- shapiro(x[, a, reta]) Perform the Shapiro-Wilk test for normality.
- anderson(x[, dist]) Anderson-Darling test for data coming from a particular distribution
- anderson_ksamp(samples[, midrank]) The Anderson-Darling test for k-samples.
Normality test ¶
Let's perform the Normality test whose Null Hypothesis is that the data is sampled from a Normal distribution. We will use a significance level of 5.0%
stats.kurtosistest(a=student['G3'])
Correlation test ¶
Let's perform the Correlation test whose Null Hypothesis is that the correlation coefficient is equal to zero. We will use a significance level of 5.0%
stats.pearsonr(student['G1'],student['G2'])
Test of Equal variances ¶
Let's perform the Bartetts's test whose Null Hypothesis is that the variances are equal. We will use a significance level of 5.0%
student.groupby('sex')['G3'].var()
grades_girls = student['G3'][student['sex']=='F']
grades_boys = student['G3'][student['sex']=='M']
stats.bartlett(grades_girls, grades_boys)
According to the test we cannot reject the Null hypothesis of equal variances, so we will use assume that the two groups are samples from a population with the same variances. This information will be useful in our next test.
Let's make again a graph showing the confidence intervals of the mean of the final grade between Boys and Girls
fig, axes = plt.subplots(1,2, figsize = (14,4))
sns.boxplot(x='sex', y='G3', data=student, ax=axes[0])
sns.pointplot(x='sex', y='G3', data=student, ax=axes[1]);
Comparing means ¶
The visualizations suggest there is a difference between the means of the final grade of the two groups. Now we will perform a formal statistical test to confirm the hypothesis that girls have higher final scores than the boys.
Null Hypothesis: for both groups (Boys and Girls) the population means of the final grade are equal.
Alternative Hypothesis: The population means of the final grades are different.
A common test to apply in for these cases is the two-sample t-test, which is used to determine if two population means are equal.
The assumptions of this test are:
- Independent samples: we will assume that the method for collecting the data assured that the answers given by the students are independent.
- Large enough sample size or observations come from a normally-distributed population: this assumption is required if we are working with small samples (less than 30), since in the smaller group we have 166 observations we can say that we have a "large enough" sample.
- Variances are equal
In addition this test have two versions: one assuming equal variances and the other assuming unequal variances.
stats.ttest_ind(grades_girls, grades_boys, equal_var=True)
Since we got such a low p-value we can reject the Null hypothesis of equal means for the two groups at a level of significance of 5%.
There is then a statistical significant difference between the grades in the two analyzed groups. The results suggest that the girls have higher grades than boys.
Do male teenagers drink more than female teenagers?
Let's first create two categorical variables Dalc_c and Walc_c from the variables Dalc and Walc by considering three categories of drinker:
- Low when the scores are 1 or 2
- Med when the scores are equal to 3
- High when the scores are 4 or 5
Dalc_cat = {
1: 'Low',
2: 'Low',
3: 'Med',
4: 'High',
5: 'High'
}
student['Dalcc'] = student['Dalc'].map(Dalc_cat)
student['Dalcc'].head()
student['Dalcc'].value_counts()
Make it as an ordered variable
student['Dalcc'] = student['Dalcc'].astype(dtype='category',
categories=['Low', 'Med', 'High'],
ordered=True)
student['Dalcc'].value_counts(sort=False)
gender_Dalc_table = pd.crosstab(student['Dalcc'], student['sex'])
gender_Dalc_table
fig, axes = plt.subplots(1,2, figsize = (14,4))
student['sex'].value_counts().plot(kind='bar', ax=axes[0], title='Gender')
student['Dalcc'].value_counts().plot(kind='bar', ax=axes[1], title='Alcohol Consumption Level');
And the plot of the cross-table
Given the alcohol consumption
100*(gender_Dalc_table.T/gender_Dalc_table.apply(sum, axis=1)).T
And Given the sex variable
100*(gender_Dalc_table/gender_Dalc_table.apply(sum, axis=0)).T
fig, axes = plt.subplots(1,2, figsize = (14,4))
gender_Dalc_table.plot(kind='bar', stacked=True, ax=axes[0]);
(100*(gender_Dalc_table.T/gender_Dalc_table.apply(sum, axis=1)).T).plot(kind='bar', stacked=True, ax=axes[1]);
fig, axes = plt.subplots(1,2, figsize = (14,4))
(gender_Dalc_table.T).plot(kind='bar', stacked=True, ax=axes[0]);
(100*(gender_Dalc_table/gender_Dalc_table.apply(sum, axis=0)).T).plot(kind='bar', stacked=True, ax=axes[1]);
Chi-square test of independence of variables in a contingency table.
This function computes the chi-square statistic and p-value for the hypothesis test of independence of the observed frequencies in the contingency table.
chi_stat, p_value, dof, expected = stats.chi2_contingency(gender_Dalc_table)
chi_stat
p_value
dof
One-Way Anova ¶
The one-way analysis of variance (ANOVA) is used to determine whether there are any statistically significant differences between the means of independent groups.
The Tukey Test (or Tukey procedure), also called Tukey’s Honest Significant Difference test, is a post-hoc test based on the studentized range distribution.
An ANOVA test can tell you if your results are significant overall, but it will not tell you exactly where those differences lie. After you have run an ANOVA and found significant results, then you can run Tukey’s HSD to find out which specific groups’s means (compared with each other) are different. The test compares all possible pairs of means.
As an example we will compare the final grades between the group of Alcohol consumption on the whole sample and we will do the same on the sub-samples of Girls and Boys.
Let us first import the libraries
import pandas as pd
import scipy.stats as stats
import researchpy as rp
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.multicomp import MultiComparison
We will now compute the confidence interval of the variable G3 and the confidence intervals by groups
rp.summary_cont(student['G3'])
rp.summary_cont(student['G3'].groupby(student['Dalcc']))
This is a first method to perform an ANOVA
stats.f_oneway(student['G3'][student['Dalcc'] == 'High'],
student['G3'][student['Dalcc'] == 'Med'],
student['G3'][student['Dalcc'] == 'Low'])
Second method
We need to create a new variable with type object
student['Dalcc'].head()
student['Dalcc2'] = student['Dalcc'].astype(dtype='object')
student['Dalcc2'].head()
results = ols('G3 ~ C(Dalcc2)', data=student).fit()
results.summary()
Let's us now compare the One-way ANOVA results between the girls and the boys sub-sample. To do so we will first consider these two sub-samples and perform two one-way ANOVA. We will then compare the coefficients and measure the difference of the impact of the consumption of Alcohol between Girls and Boys teenagers.
Girls sub-samples
student_girls=student[student['sex']=='F']
student_girls.head()
Let's then perform the One way-ANOVA
results_girls = ols('G3 ~ C(Dalcc2)', data=student_girls).fit()
results_girls.summary()
Boys sub-samples
student_boys=student[student['sex']=='M']
student_boys.head()
Let's then perform the One way-ANOVA
results_boys = ols('G3 ~ C(Dalcc2)', data=student_boys).fit()
results_boys.summary()
Comparison of the coefficients
I will draw in a same figure a plot with the confidence intervals of the coefficients of the Boys and Girls models. We will first create a dataframe containing the estimation of the models.
x=['Boys','Girls']
coef_df = pd.DataFrame()
for i, mod in enumerate([results_boys, results_girls]):
err_series = mod.params - mod.conf_int()[0]
coef_df = coef_df.append(pd.DataFrame({'coef': mod.params.values[1:],
'err': err_series.values[1:],
'varname': err_series.index.values[1:],
'model': x[i]
})
)
coef_df
## marker to use
marker_list = 'so'
width=0.05
## 5 covariates in total
base_x = pd.np.arange(2) - 0.2
base_x
from matplotlib import pyplot as plt
from matplotlib.lines import Line2D
%matplotlib inline
fig, ax = plt.subplots(figsize=(8, 5))
for i, mod in enumerate(coef_df.model.unique()):
mod_df = coef_df[coef_df.model == mod]
mod_df = mod_df.set_index('varname').reindex(coef_df['varname'].unique())
## offset x posistions
X = base_x + width*i
ax.bar(X, mod_df['coef'],
color='none',yerr=mod_df['err'])
## remove axis labels
ax.set_ylabel('')
ax.set_xlabel('')
ax.scatter(x=X,
marker=marker_list[i], s=120,
y=mod_df['coef'], color='black')
ax.axhline(y=0, linestyle='--', color='black', linewidth=4)
ax.xaxis.set_ticks_position('none')
_ = ax.set_xticklabels(['','','Low','','','','Med'],
rotation=0, fontsize=16)
## finally, build customized legend
legend_elements = [Line2D([0], [0], marker=m,
label=x[i],
color = 'k',
markersize=10)
for i, m in enumerate(marker_list)
]
_ = ax.legend(handles=legend_elements, loc=2,
prop={'size': 15}, labelspacing=1.2)
Tukey’s HSD Post-hoc comparison ¶
The overall model was significant, now to test which groups differ. Deciding which groups to compare should be theory driven. There are a few different techniques that can be used. Each of these techniques have different ways of controlling for familywise error rate. We will show here how to perform one them:
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.multicomp import MultiComparison
mc_girls = MultiComparison(student_girls['G3'], student_girls['Dalcc2'])
mc_results_girls = mc_girls.tukeyhsd()
print(mc_results_girls)
mc_boys = MultiComparison(student_boys['G3'], student_boys['Dalcc2'])
mc_results_boys = mc_boys.tukeyhsd()
print(mc_results_girls)