#import pandas and use magic function
import pandas as pd
%matplotlib inline

# import our data using pandas read_csv() function where delimiter = ';', index_col = 0, na_values = '.'
data = pd.read_csv('https://www.scipy-lectures.org/_downloads/brain_size.csv', 
                   delimiter=';', index_col=0, na_values='.')

# check out our data using pandas df.head() function
data.head()

# how many observations do we have? use pandas df.shape attribute
data.shape

# check out one column with df['column name'] or df.column_nae
data.Gender.head()

# make a groupby object on the dataframe
groupby_gender = data.groupby('Gender')

# take the mean of the groupby object across all measures using the .mean() method
groupby_gender.mean()

# take a look at our data distributions and pair-wise correlations
from pandas import plotting
plotting.scatter_matrix(data[['Weight', 'Height', 'MRI_Count']]);

plotting.scatter_matrix(data[['PIQ', 'VIQ', 'FSIQ']]);

## Distributions and 

sns.kdeplot(data['FSIQ'])
sns.kdeplot(data['PIQ'])
sns.kdeplot(data['VIQ'])

sns.kdeplot(data['FSIQ'],cumulative=True)
sns.kdeplot(data['PIQ'],cumulative=True)
sns.kdeplot(data['VIQ'],cumulative=True)

"""
Explanation: Introduction to Statistics in Python
https://www.scipy-lectures.org/packages/statistics/index.html
Importing and Exploring our Data
End of explanation
"""


from scipy import stats

"""
Explanation: Hypothesis Testing
End of explanation
"""


# runa 1 sample t-test
stats.ttest_1samp(data['VIQ'], 0)   

"""
Explanation: 1-sample t-test: Testing the value of a population mean.
scipy.stats.ttest_1samp() tests if the population mean of data is likely to 
be equal to a given value (technically if observations are drawn from a Gaussian
distributions of given population mean). It returns the T statistic, and the 
p-value (see the function’s help)
End of explanation
"""


female_viq = data[data['Gender'] == 'Female']['VIQ']
male_viq = data[data['Gender'] == 'Male']['VIQ']
stats.ttest_ind(female_viq, male_viq)   



"""
Explanation: 2-sample t-test: testing for difference across populations
We have seen above that the mean VIQ in the male and female populations were different.
To test if this is significant, we do a 2-sample t-test with scipy.stats.ttest_ind()
End of explanation
"""


stats.ttest_ind(data['FSIQ'], data['PIQ'])  

"""
Explanation: Paired tests: repeated measurements on the same indivuals
The PIQ, VIQ and FSIQ are three different measures of IQ in the same individual. 
We can first look if FSIQ and PIQ are different using the 2-sample t-test.
End of explanation
"""


stats.ttest_rel(data['FSIQ'], data['PIQ'])   

"""
Explanation: However this doesn't account for individual differences contributing to variance in data.
We can use a paired t-test or repeated measures test to account for these individual differences.
End of explanation
"""


stats.ttest_1samp(data['FSIQ'] - data['PIQ'], 0)   

"""
Explanation: This is actually equivalent to doing a 1-sample t-test on the difference of the two measures.
End of explanation
"""


stats.wilcoxon(data['FSIQ'], data['PIQ'])   

"""
Explanation: These tests assume normality in the data. A non-parametric alternative is the Wilcoxian signed rank test
End of explanation
"""


# Let's simulate some data according to the model
import numpy as np

x = np.linspace(-5, 5, 20)
np.random.seed(1)

# normal distributed noise
y = -5 + 3*x + 4 * np.random.normal(size=x.shape)

# Create a data frame containing all the relevant variables
sim_data = pd.DataFrame({'x': x, 'y': y})

# Specify an OLS model and fit it
from statsmodels.formula.api import ols
model = ols("y ~ x", sim_data).fit()

# Inspect the results of the model fit
print(model.summary()) 

# Retrieve the model params, note tab completion
model.params

"""
Explanation: Note: 
The corresponding test in the non paired case is the Mann–Whitney U test, scipy.stats.mannwhitneyu(). 
Linear models, multiple factors, and analysis of variance
Given two set of observations, x and y, we want to test the hypothesis that y is a linear function of x.
We will use the statsmodels module to:

Fit a linear model. We will use the simplest strategy, ordinary least squares (OLS).
Test that coef is non zero.
End of explanation
"""


data.head()

# We can write a comparison between IQ of male and female using a linear model:

model = ols("VIQ ~ Gender", data).fit()

# ols automatically detects Gender as categorical
# model = ols('VIQ ~ C(Gender)', data).fit()

print(model.summary())

"""
Explanation: Categorical variables: comparing groups or multiple categories
End of explanation
"""


iq_melt = pd.melt(data, value_vars=['FSIQ', 'PIQ'], value_name='iq', var_name="type")

iq_melt.head()

model = ols("iq ~ type", iq_melt).fit()
print(model.summary())  



"""
Explanation: Link to t-tests between different FSIQ and PIQ
To compare different types of IQ, we need to create a “long-form” table, 
listing IQs, where the type of IQ is indicated by a categorical variable:
End of explanation
"""


stats.ttest_ind(data['VIQ'], data['PIQ']) 

"""
Explanation: Note that we can get the same results now of doing the independant t-tests for each pariing.
End of explanation
"""


iris = pd.read_csv('https://www.scipy-lectures.org/_downloads/iris.csv')

import seaborn as sns
sns.pairplot(iris, hue='name');

"""
Explanation: Multiple Regression: including multiple factors
End of explanation
"""


model = ols('sepal_width ~ name + petal_length', iris).fit()
print(model.summary())  


"""
Explanation: Sepal and petal size tend to be related: bigger flowers are bigger!
But is there in addition a systematic effect of species?
End of explanation
"""


data.head()

sns.boxplot('name', "sepal_length", data=iris)

setosa = iris[iris['name'] == 'setosa']['sepal_length']
versicolor = iris[iris['name'] == 'versicolor']['sepal_length']
virginica = iris[iris['name'] == 'virginica']['sepal_length']

f_value, p_value = stats.f_oneway(setosa, versicolor, virginica)

print(f_value, p_value)

data.head()

data = data.dropna(axis=0, how='any')

data.isna().sum()

import statsmodels.api as sm
from statsmodels.formula.api import ols
 
mod = ols('sepal_length ~ name',
                data=iris).fit();
                
aov_table = sm.stats.anova_lm(mod, typ=2);
print(aov_table);

"""
Explanation: Analysis of Variance (ANOVA)
1-Way Anova
End of explanation
"""


