#import pandas and use magic function for plotting
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()

data.describe()

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

# check out one column with df['column name'] or df.column_name
#data['Gender']
data.Gender

# make a groupby object on the dataframe
data_grouped = data.groupby("Gender")


# look at mean weight per group in data_grouped
data_grouped['Height'].min()

# import plotting from pandas
from pandas import plotting

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

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

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


# import seaborn plotting package
import seaborn as sns

# use the kdeplot() function on FSIQ, PIQ and VIQ
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);

# cumulative = True


"""
Explanation: Plotting and Comparing Distributions
End of explanation
"""


# import stats from scipy 
from scipy import stats

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


# run a 1 sample t-test on VIQ against 0
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
"""


# set up two series for female and male VIQ
female = data[data.Gender == 'Female']['VIQ']
male = data[data.Gender == 'Male']['VIQ']


# use stats.ttest_ind() on the two series
stats.ttest_ind(female, male)

"""
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
"""


# run a 2-sample t-test on FSIQ and PIQ
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
"""


# use the ttest_rel() to do a paired t-test on FSIQ and PIQ
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
"""


# 1-sample ttest on the difference
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
"""


# for a non-parametric version of this test, use stats.wilcoxon()
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
import numpy as np
# create range of 20 x values over -5 to 5 using np.linspace()
x = np.linspace(-5,5,20)
# set random seed to 1
np.random.seed(1)
# normal distributed noise of y = -5 + 3x + 4 * np.random.normal(size=x.shape)
y = -5 + 3*x + 4 * np.random.normal(size=x.shape)
# Create a DataFrame named sim_data containing all the relevant variables
fake_data = pd.DataFrame({"x": x, "y": y})

# import ols from statsmodels.formmula.api
from statsmodels.formula.api import ols

# Specify an OLS model (y ~ x) and fit it using the .fit() method 
model = ols('y ~ x', fake_data).fit()

# Inspect the results of the model using model.summary()
model.summary()
print(model.params)

# Retrieve the model params, note tab completion


"""
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
"""


# remind ourselves of our data using head()
data.head()

# We can write a comparison between VIQ of male and female 
# (VIQ ~ Gender) using a linear model:
model = ols("VIQ ~ Gender", data).fit()

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

# view results of the model using .summary()
model.summary()

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


# read in our iris data https://www.scipy-lectures.org/_downloads/iris.csv


# use seaborns pairplot() function and hue='name', use ; to remove output


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


# let's make a model that takes into account this potential interaction of name (sepal_width ~ name + petal_length)


# view the results of the model


"""
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
"""


# remind ourselves of our iris data
iris = pd.read_csv("https://www.scipy-lectures.org/_downloads/iris.csv")

iris.head()

# let's visualize the difference species and their sepal_length in 
# a box plot using seaborns box_plot()
sns.boxplot(x="name", y="sepal_length", data=iris);

# we need to split our data to run an anova using stats.f_oneway()
# make a variable for setosa, versicolor and virginica on sepal_length
setosa = iris[iris['name'] == 'setosa']['sepal_length']
versicolor = iris[iris['name'] == 'versicolor']['sepal_length']
virginica = iris[iris['name'] == 'virginica']['sepal_length']

# use stats.f_oneway() on the series you made for each species and 
# get an f_value and p_value
f_value, p_value = stats.f_oneway(setosa, versicolor, virginica)

# print the results
print(f_value, p_value)

# another way to do it that prevents having to split your data
# import sm from statsmodels.api and import ols from statsmodels.formula.api
import statsmodels.api as sm
from statsmodels.formula.api import ols

# specify the model (sepal_length ~ name) and fit
model = ols('sepal_length ~ name', iris).fit()

# make an aov table using sm.stats.anova_lm() using type=2 since we don't have an interaction term
aov_table = sm.stats.anova_lm(model, typ=2)

# print the aov_table
print(aov_table)

model = ols('FSIQ ~ Gender + MRI_Count + Gender:MRI_Count', data).fit()

aov_table = sm.stats.anova_lm(model, typ=1)


print(aov_table)

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


