# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Data analyses in Python
#
# Imagine you have a `.csv` file (could also be `.txt`, `.log`, `.dat`, `.sav`, etc.) with some data in it and want to analyze that data. Your `worklfow` to do so will most likely consist of different steps: `loading the data`, `"cleaning" the data`, `exploring the data`, `analyzing the data` and `visualization of the data/results`. In the following, we'll briefly go through all of these steps.
#
# As you might be aware, the first step is always to load necessary `packages` and/or `functions`. Most of the time, it is not clear what exactly is needed along the `worklfow`. Hence, starting with respective `packages/functions` we're sure about is a good idea. This is most likely based on the data you want to analyze. As we want to have a look at some data in a `.csv` file, `numpy` is a good starting point.
# However, to already provide you with the full list of `packages` we are going to explore, here you::
#
# - [numpy](https://numpy.org/)
# - [pandas](https://pandas.pydata.org/)
# - [scipy's stats module](https://docs.scipy.org/doc/scipy/reference/stats.html)
# - [statsmodels](http://www.statsmodels.org/stable/index.html)
# - [seaborn](seaborn.pydata.org)
#
# Ok, nuff said, let's go to work:

import numpy as np

# Using `np.` + `tab` provides us with a nice overview of `numpy`'s functionality:

np.

# In case, we are not sure about how a certain function works or simply want to know more about it, we can use the `help` function:

help(np.array)

# Based on our goal, the `genfromtxt` function looks useful, as we initially need to load the data from the `.csv` file:

my_data_numpy = np.genfromtxt('data/brain_size.csv', delimiter=';')

# With that, we generated a variable called `my_data_numpy`. Now, let's check its `type`.

type(my_data_numpy)

# It is a `numpy.ndarray` and within it, the data is stored:

my_data_numpy

# As we saw in the `help` function above, `objects` and `data types` come with certain `functionality`:

my_data_numpy.

# We can, for example, check the `shape`, that is the `dimensionality`, of the data:

my_data_numpy.shape

# This returns not a `numpy.ndarray`, but a `tuple`:

type(my_data_numpy.shape)

# It is also possible to `concatenate functionality`. E.g., we could `transpose` our `numpy.ndarray` and check its resulting `shape` within one command:

my_data_numpy.transpose().shape

# Is it possible to `view only certain parts` of the data, i.e. the second row? Yes, using `slicing`.

my_data_numpy[1]

# The output is a `numpy.ndarray` again:

type(my_data_numpy[1])

# If we want to be more specific, it is also possible to only view one value, i.e. the fourth value of the second row:

my_data_numpy[1, 3]

# Now, the `data type` has changed to `numpy.float64`:

type(my_data_numpy[1, 3])

# However, getting more than one value, i.e. multiple values of the second row, results in a `numpy.ndarray` again:

my_data_numpy[1, 3:6]

# Let's look at our data again:

my_data_numpy

# Even though it's a small dataset, there's already a lot going on: different `data types`, different `columns`, etc. and apparently not everything is `"numpy"` compatible. `Numpy` is a great tool and very powerful, building the foundation for a lot of `python libraries`, but in a lot of cases you might want to prefer `packages` that build upon `numpy` and are intended for specific purposes. A good example for that is the amazing `pandas` library that should be your first address for everything `data wrangling`. In particular, this refers to `tabular data`, that is `multiple observations` or `samples` described by a set of different `attributes` or `features`. The data can than be seen as a `2D table`, or `matrix`, with `columns` giving the different `attributes` of the data, and rows the `observations`. Let's try `pandas`, but at first, we have to import it:

import pandas as pd

# Now we can check its functionality:

pd.

# `read_csv` looks helpful regarding loading the data:

my_data_pandas = pd.read_csv('data/brain_size.csv', delimiter=',')

# What do we have? A `pandas dataframe`:

type(my_data_pandas)

# Before, our data was in `np.ndarray` format:

type(my_data_numpy)

# How does our data look now?

my_data_pandas

# Even though we already have more information as in the `numpy array`, e.g., `headers`, `strings` and `indexes`, it still looks off. What's the problem? Well, we see that our data has a `;` as a `delimiter`, but we indicated `,` as delimiter when loading our data. Therefore, it is important to carefully check your data and beware of its specifics. Let's reload our data with the fitting `delimiter`:

my_data_pandas = pd.read_csv('data/brain_size.csv', delimiter=';')

# Investigating our `dataframe`, we see that it worked as expected this time:

my_data_pandas

# Thinking about our `numpy.ndarray` version, we see a difference in the `shape` of the data, which is related to the `header`:

my_data_pandas.shape

# What can we do with our `dataframe`:

my_data_pandas.

# For example we can and should rename `columns` with uninformative names:

my_data_pandas.rename(columns={'Unnamed: 0': 'sub-id'})

# That looks a bit more informative, doesn't it? Let's have a look at our columns again

my_data_pandas.columns

# Wait a minute, it's not `renamed`. Did we do something wrong? Let's check the respective functionality:

help(my_data_pandas.rename)

# Checking the functionality more in depth, a `dataframe` with the new `column names` is returned, but the old one `not automatically changed`. Hence, we have to do it again, this overwriting the original `dataframe`:

my_data_pandas=my_data_pandas.rename(columns={'Unnamed: 0': 'sub-id'})
my_data_pandas.columns

# Pandas also allows the easy and fast `exploration` of our data:

my_data_pandas.describe()

# Unfortunately, not all `columns` are there. But why is that? We need to investigate the `columns` more closely, beginning with one that was included:

type(my_data_pandas['sub-id'])

# The data in the `columns` is a `pandas series`, not a `dataframe` or `numpy.ndarray`, again with its own functionality. Nevertheless, it was included in our `numerical summary`. Let's check the first missing `column`:

type(my_data_pandas['Hair'])

# Well, that's not very informative on its own, as it's also a `pandas series`, but was not included. Maybe the `data type` is the problem? Luckily, the `pandas dataframe` object comes with a helpful functionality:

my_data_pandas.dtypes

# And a bit more closely using `indexing`:

type(my_data_pandas['Hair'][0])

# The data in `my_data_pandas['Hair']` has the `type str` and as you might have already guessed: it's rather hard to compute `summary statistics` from a `str`. We could re-code it, but given there are only two values, this might not be very useful for our current aim:

my_data_pandas['Hair'].unique()

# What about the other `missing columns`, e.g., `height`?

type(my_data_pandas['Height'][0])

# The `data type` is yet again `str`, but how many values do we have?

my_data_pandas['Height'].unique()

# Hm, we can see that `height` contains `numerical values`. However, the `data type` is `str`. Here it can be useful to change the `data type`, using `pandas dataframe` object functionality:

my_data_pandas['Height'].astype(float)

# And we're getting another `error`. This time, it's related to a `missing data point`, which needs to be addressed before the `conversion` is possible. We can simply use the `replace` functionality to `replace` the `missing data point` to `NaN`, which should as allow to do the `conversion`:

my_data_pandas['Height'] = my_data_pandas['Height'].replace('.', np.nan)

my_data_pandas['Height'] = my_data_pandas['Height'].astype(float)

# Let's check if the `column` is now included in the `summary`:

my_data_pandas.describe()

# Now, we can do the same for the `Weight` column, `concatenating` all necessary functions in one line:

my_data_pandas['Weight'] = my_data_pandas['Weight'].replace('.', np.nan).astype(float)

# Is `Weight` now included?

my_data_pandas.describe()

# We can also compute one statistical value for one column, for example the `mean` using `numpy`:

np.mean(my_data_pandas['Weight'])

# But the same is also possible using inbuilt `pandas data frame` functionality:

my_data_pandas['Weight'].mean()

# We can do the same for the standard deviation:

np.std(my_data_pandas['Weight'])

my_data_pandas['Weight'].std()

# Here we can see, the same `functionality` can lead to different `results`, potentially based on `different implementations`. Thus, always make sure to check every part of your code and re-run it to see if you get the same outputs. As you can see here, using a `jupyter notebook` for your analyses, this is comparably straightforward. Additionally, you can document each step of your workflow, from data loading, inspection, changes, etc. . While you should of course always use `version control` on your data, the format we've explored nicely allows to redo your analyses (excluding the `computational reproducibility` and `numerical instability` aspect). On top of that, you can document the executed steps so that your future self and everyone else knows what's going on. Enough chit-chat, now that we've loaded and inspected our data, as well as fixed some errors it's time to do some statistics. To show you a few nice `packages` that are out there, we will run different `analyses` with different `packages`. We will explore `pingouin`, `scipy`, `statsmodels` and `seaborn`.     

# <img src="https://github.com/raphaelvallat/pingouin/blob/master/docs/pictures/logo_pingouin.png?raw=true" height="300" width="700"/>
#
#
#
# ### _Pingouin is an open-source statistical package written in Python 3 and based mostly on Pandas and NumPy._
#
#
# - ANOVAs: one- and two-ways, repeated measures, mixed, ancova
# - Post-hocs tests and pairwise comparisons
# - Robust correlations
# - Partial correlation, repeated measures correlation and intraclass correlation
# - Linear/logistic regression and mediation analysis
# - Bayesian T-test and Pearson correlation
# - Tests for sphericity, normality and homoscedasticity
# - Effect sizes and power analysis
# - Parametric/bootstrapped confidence intervals around an effect size or a correlation coefficient
# - Circular statistics
# - Plotting: Bland-Altman plot, Q-Q plot, etc...
#
# **Pingouin is designed for users who want simple yet exhaustive statistical functions.**
#
#
# ##### **material scavenged from [10 minutes to Pingouin](https://pingouin-stats.org/index.html) and [the pingouin docs](https://pingouin-stats.org/api.html)
#
# Let's import the `package`:

import pingouin as pg

# ### Correlations
#
# "In the broadest sense correlation is any statistical association, though in common usage it most often refers to how close two variables are to having a linear relationship with each other" - [Wikipedia](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient)
#
# `Pingouin` supports a variety of [measures of correlation](https://pingouin-stats.org/generated/pingouin.corr.html#pingouin.corr). When talking about `correlation`, we commonly mean the `Pearson correlation coefficient`, also referred to as `Pearson's r`:
#
# <img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/93185aed3047ef42fa0f1b6e389a4e89a5654afa"/>
#
# Computing `Pearson's r` using `pingouin` is as easy as:

pearson_correlation = pg.corr(my_data_pandas['FSIQ'], my_data_pandas['VIQ'])
display(pearson_correlation)
cor_coeeficient = pearson_correlation['r']

# The output we got, is the `test summary`:
#
# - 'n' : Sample size (after NaN removal)
# - 'outliers' : number of outliers (only for 'shepherd' or 'skipped')
# - 'r' : Correlation coefficient
# - 'CI95' : [95% parametric confidence intervals](https://en.wikipedia.org/wiki/Confidence_interval)
# - 'r2' : [R-squared](https://en.wikipedia.org/wiki/Coefficient_of_determination)
# - 'adj_r2' : [Adjusted R-squared](https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2)
# - 'p-val' : one or two tailed p-value
# - 'BF10' : Bayes Factor of the alternative hypothesis (Pearson only)
# - 'power' : achieved power of the test (= 1 - type II error)

# What if we want to compute `pairwise correlations` between `columns` of a `dataframe`? With `pingouin` that's one line of code and we can even sort the results based on a `test statistic` of interest, e.g. `r2`:

pg.pairwise_corr(my_data_pandas, columns=['FSIQ', 'VIQ', 'Weight']).sort_values(by=['r2'], ascending=False)

# ### Before we calculate:  `Testing statistical premises`
#
# Statistical procedures can be classfied into either [`parametric`](https://en.wikipedia.org/wiki/Parametric_statistics) or `non parametric` prcedures, which require different necessary preconditions to be met in order to show consistent/robust results. 
# Generally people assume that their data follows a gaussian distribution, which allows for parametric tests to be run.
# Nevertheless it is essential to first test the distribution of your data to decide if the assumption of normally distributed data holds, if this is not the case we would have to switch to non parametric tests.

# ### [<NAME> normality  test](https://pingouin-stats.org/generated/pingouin.normality.html#pingouin.normality)
#
# Standard procedure for testing `normal distribution` tests if the `distribution` of your data `deviates significantly` from a `normal distribution`.
# The function we're using returns the following information:
#
# - W  : Test statistic
#
# - p : float
#     P-value
#
# - normal : boolean
#     True if data comes from a normal distribution.

pg.normality(my_data_pandas['Height'], alpha=.05)

# ### [Henze-Zirkler multivariate normality test](https://pingouin-stats.org/generated/pingouin.multivariate_normality.html#pingouin.multivariate_normality)
#
# The same procedure, but for [multivariate normal distributions](https://en.wikipedia.org/wiki/Multivariate_normal_distribution).  

pg.multivariate_normality(my_data_pandas[['Height', 'Weight','VIQ']], alpha=.05)


# ### [Testing for homoscedasticity](https://pingouin-stats.org/generated/pingouin.homoscedasticity.html?highlight=homoscedasticity#pingouin.homoscedasticity)
#
#
# "In statistics, a sequence or a vector of random variables is homoscedastic /ˌhoʊmoʊskəˈdæstɪk/ if all random variables in the sequence or vector have the same finite variance." - Wikipedia
#
# returns:
#
# equal_var : boolean True if data have equal variance.
#
# p : float P-value.
#
# Note: This function first tests if the data are normally distributed using the Shapiro-Wilk test. If yes, then the homogeneity of variances is measured using the Bartlett test. If the data are not normally distributed, the Levene test, which is less sensitive to departure from normality, is used.

pg.homoscedasticity(my_data_pandas[['VIQ', 'FSIQ']], alpha=.05)


# ## Parametric tests
# ## Student's t-test: the simplest statistical test
#
# ### 1-sample t-test: testing the value of a population mean
#
# 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).
#
#
# `pingouin.ttest` returns the T_statistic, the p-value, the [degrees of freedom](https://en.wikipedia.org/wiki/Degrees_of_freedom_(statistics), the [Cohen d effect size](https://en.wikiversity.org/wiki/Cohen%27s_d), the achieved [power](https://en.wikipedia.org/wiki/Power_(statistics%29) of the test ( = 1 - type II error (beta) = [P(Reject H0|H1 is true)](https://deliveroo.engineering/2018/12/07/monte-carlo-power-analysis.html)), and the [Bayes Factor](https://en.wikipedia.org/wiki/Bayes_factor) of the alternative hypothesis
#
#
#
#

pg.ttest(my_data_pandas['VIQ'],0)

# ### 2-sample t-test: testing for difference across populations
#
# We have seen above that the mean VIQ in the dark hair and light hair populations
# were different. To test if this is significant, we do a 2-sample t-test:

light_viq = my_data_pandas[my_data_pandas['Hair'] == 'light']['VIQ']
dark_viq = my_data_pandas[my_data_pandas['Hair'] == 'dark']['VIQ']

pg.ttest(light_viq, dark_viq)

# ### Plot achieved power of a paired T-test
#
# Plot the curve of achieved power given the effect size (Cohen d) and the sample size of a paired T-test.

# +
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='ticks', context='notebook', font_scale=1.2)

d = 0.5  # Fixed effect size
n = np.arange(5, 80, 5)  # Incrementing sample size

# Compute the achieved power
pwr = pg.power_ttest(d=d, n=n, contrast='paired', tail='two-sided')

# Start the plot
plt.plot(n, pwr, 'ko-.')
plt.axhline(0.8, color='r', ls=':')
plt.xlabel('Sample size')
plt.ylabel('Power (1 - type II error)')
plt.title('Achieved power of a paired T-test')
sns.despine()
# -

# ### Non parametric tests:
#
#
# Unlike the parametric test these do not require the assumption of normal distributions.
#
# "`Mann-Whitney U Test` (= Wilcoxon rank-sum test). It is the non-parametric version of the independent T-test.
# Mwu tests the hypothesis that data in x and y are samples from continuous distributions with equal medians. The test assumes that x and y are independent. This test corrects for ties and by default uses a continuity correction." - [mwu-function](https://pingouin-stats.org/generated/pingouin.mwu.html#pingouin.mwu)
#
# Test summary
#
# - 'W-val' : W-value
# - 'p-val' : p-value
# - 'RBC'   : matched pairs rank-biserial correlation (effect size)
# - 'CLES'  : common language effect size

pg.mwu(light_viq, dark_viq)

# "`Wilcoxon signed-rank test` is the non-parametric version of the paired T-test.
#
# The Wilcoxon signed-rank test tests the null hypothesis that two related paired samples come from the same distribution. A continuity correction is applied by default." - [wilcoxon - func](https://pingouin-stats.org/generated/pingouin.wilcoxon.html#pingouin.wilcoxon)
#

pg.wilcoxon(light_viq, dark_viq, tail='two-sided')

# ### `scipy.stats` - Hypothesis testing: comparing two groups
#
# For simple [statistical tests](https://en.wikipedia.org/wiki/Statistical_hypothesis_testing), it is also possible to use the `scipy.stats` sub-module of [`scipy`](http://docs.scipy.org/doc/).

from scipy import stats

# ### 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](https://en.wikipedia.org/wiki/Student%27s_t-test), and the [p-value](https://en.wikipedia.org/wiki/P-value) (see the function's help):

stats.ttest_1samp(my_data_pandas['VIQ'], 100)

# With a p-value of 10^-28 we can claim that the population mean for the IQ (VIQ measure) is not 0.

# ### 2-sample t-test: testing for difference across populations
#
# We have seen above that the mean VIQ in the dark hair and light hair populations
# were different. To test if this is significant, we do a 2-sample t-test
# with `scipy.stats.ttest_ind`:

light_viq = my_data_pandas[my_data_pandas['Hair'] == 'light']['VIQ']
dark_viq = my_data_pandas[my_data_pandas['Hair'] == 'dark']['VIQ']
stats.ttest_ind(light_viq, dark_viq)

# ## Paired tests: repeated measurements on the same indivuals
#
# PIQ, VIQ, and FSIQ give 3 measures of IQ. Let us test if FISQ and PIQ are significantly different. We can use a 2 sample test:

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

# The problem with this approach is that it forgets that there are links
# between observations: FSIQ and PIQ are measured on the same individuals.
#
# Thus the variance due to inter-subject variability is confounding, and
# can be removed, using a "paired test", or ["repeated measures test"](https://en.wikipedia.org/wiki/Repeated_measures_design):

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

# This is equivalent to a 1-sample test on the difference::

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

# T-tests assume Gaussian errors. We can use a [Wilcoxon signed-rank test](https://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test), that relaxes this assumption:

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

# **Note:** The corresponding test in the non paired case is the [Mann–Whitney U test](https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U), `scipy.stats.mannwhitneyu`.

# + [markdown] solution2="hidden" solution2_first=true
# ### Exercise 2
#
# * Test the difference between weights in people with dark and light hair.
# * Use non-parametric statistics to test the difference between VIQ in people with dark and light hair.

# + solution2="hidden"
light_weight = my_data_pandas[my_data_pandas['Hair'] == 'light']['Weight']
dark_weight = my_data_pandas[my_data_pandas['Hair'] == 'dark']['Weight']
stats.ttest_ind(light_weight, dark_weight, nan_policy='omit')

# + solution2="hidden"
stats.mannwhitneyu(light_viq, dark_viq)

# + [markdown] solution2="hidden"
# **Conclusion**: we find that the data does not support the hypothesis that people with dark and light hair have different VIQ.

# +
# Create solution here
# -

# # `statsmodels` - use "formulas" to specify statistical models in Python
#
# Use `statsmodels` to perform linear models, multiple factors or analysis of variance.
#
#
# ## A simple linear regression
#
# Given two set of observations, `x` and `y`, we want to test the hypothesis that `y` is a linear function of `x`.
#
# In other terms:
#
# $y = x * coef + intercept + e$
#
# where $e$ is observation noise. We will use the [statsmodels](http://statsmodels.sourceforge.net) module to:
#
# 1. Fit a linear model. We will use the simplest strategy, [ordinary least squares](https://en.wikipedia.org/wiki/Ordinary_least_squares) (OLS).
# 2. Test that $coef$ is non zero.
#
# First, we generate simulated data according to the model. Then we specify an OLS model and fit it:

from statsmodels.formula.api import ols
model = ols("FSIQ ~ VIQ", my_data_pandas).fit()

# **Note:** For more about "formulas" for statistics in Python, see the [statsmodels documentation](http://statsmodels.sourceforge.net/stable/example_formulas.html).

# We can inspect the various statistics derived from the fit:

print(model.summary())

# ### Terminology
#
# Statsmodels uses a statistical terminology: the `y` variable in statsmodels is called *endogenous* while the `x` variable is called *exogenous*. This is discussed in more detail [here](http://statsmodels.sourceforge.net/devel/endog_exog.html). To simplify, `y` (endogenous) is the value you are trying to predict, while `x` (exogenous) represents the features you are using to make the prediction.

# + [markdown] solution2="hidden" solution2_first=true
# ### Exercise 3
#
# Retrieve the estimated parameters from the model above.  
# **Hint**: use tab-completion to find the relevant attribute.

# + solution2="hidden"
model.params

# +
# Create solution here
# -

# ## Categorical variables: comparing groups or multiple categories

model = ols("VIQ ~ Hair + 1", my_data_pandas).fit()
print(model.summary())

# ### Tips on specifying model
#  
# ***Forcing categorical*** - the 'Hair' is automatically detected as a categorical variable, and thus each of its different values is treated as different entities.
#
# An integer column can be forced to be treated as categorical using:
#
# ```python
# model = ols('VIQ ~ C(Hair)', my_data_pandas).fit()
# ```
#
# ***Intercept***: We can remove the intercept using `- 1` in the formula, or force the use of an intercept using `+ 1`.

# ### 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:

data_fisq = pd.DataFrame({'iq': my_data_pandas['FSIQ'], 'type': 'fsiq'})
data_piq = pd.DataFrame({'iq': my_data_pandas['PIQ'], 'type': 'piq'})
data_long = pd.concat((data_fisq, data_piq))
print(data_long[::8])

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

# We can see that we retrieve the same values for t-test and corresponding p-values for the effect of the type of IQ than the previous t-test:

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

# ## Multiple Regression: including multiple factors
#
# Consider a linear model explaining a variable `z` (the dependent
# variable) with 2 variables `x` and `y`:
#
# $z = x \, c_1 + y \, c_2 + i + e$
#
# Such a model can be seen in 3D as fitting a plane to a cloud of (`x`,
# `y`, `z`) points (see the following figure).

# +
from mpl_toolkits.mplot3d import Axes3D
x = np.linspace(-5, 5, 21)

# We generate a 2D grid
X, Y = np.meshgrid(x, x)

# To get reproducable values, provide a seed value
np.random.seed(1)

# Z is the elevation of this 2D grid
Z = -5 + 3*X - 0.5*Y + 8 * np.random.normal(size=X.shape)

# Plot the data
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_trisurf(my_data_pandas['VIQ'].to_numpy(), my_data_pandas['PIQ'].to_numpy(), 
                       my_data_pandas['FSIQ'].to_numpy(), cmap=plt.cm.plasma)
ax.view_init(20, -120)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
# -

model = ols('FSIQ ~ VIQ + PIQ', my_data_pandas).fit()
print(model.summary())

# ## Post-hoc hypothesis testing: analysis of variance (ANOVA)
#
# In the above iris example, we wish to test if the petal length is different between versicolor and virginica, after removing the effect of sepal width. This can be formulated as testing the difference between the coefficient associated to versicolor and virginica in the linear model estimated above (it is an Analysis of Variance, [ANOVA](https://en.wikipedia.org/wiki/Analysis_of_variance). For this, we write a **vector of 'contrast'** on the parameters estimated with an [F-test](https://en.wikipedia.org/wiki/F-test):

print(model.f_test([0, 1, -1]))

# Is this difference significant?

# + [markdown] solution2="hidden" solution2_first=true
# ### Exercise 4
#
# Going back to the brain size + IQ data, test if the VIQ of people with dark and light hair are different after removing the effect of brain size, height, and weight.

# + solution2="hidden"
data = pd.read_csv('data/brain_size.csv', sep=';', na_values=".")
model = ols("VIQ ~ Hair + Height + Weight + MRI_Count", data).fit()
print(model.summary())

# +
# Create solution here
# -

# ### Throwback to pingouin and pandas
#
# Remember `pingouin`? As briefly outlined, it can also compute `ANOVA`s and other types of models fairly easy. For example, let's compare `VIQ` between `light` and `dark` `hair`ed participants. 

pg.anova(dv='VIQ', between='Hair', data=my_data_pandas,
               detailed=True)

# It gets even better, `pandas` actually support some `pingouin` functions directly as an in-built `dataframe method`: 

my_data_pandas.anova(dv='VIQ', between='Hair', detailed=True)

# # `seaborn` - use visualization for statistical exploration
#
# [Seaborn](http://stanford.edu/~mwaskom/software/seaborn/) combines simple statistical fits with plotting on `pandas dataframes`, `numpy arrays`, etc. .

# ## Pairplot: scatter matrices
#
# We can easily have an intuition on the interactions between continuous variables using `seaborn.pairplot` to display a scatter matrix:

import seaborn
seaborn.set()
seaborn.pairplot(my_data_pandas, vars=['FSIQ', 'PIQ', 'VIQ'], kind='reg')

# Categorical variables can be plotted as the hue:

seaborn.pairplot(my_data_pandas, vars=['FSIQ', 'VIQ', 'PIQ'], kind='reg', hue='Hair')

# ## lmplot: plotting a univariate regression
#
# A regression capturing the relation between one variable and another, e.g. `FSIQ` and `VIQ`, can be plotted using `seaborn.lmplot`:

seaborn.lmplot(y='FSIQ', x='VIQ', data=my_data_pandas)

# ### Robust regression
# Given that, in the above plot, there seems to be a couple of data points that are outside of the main cloud to the right, they might be outliers, not representative of the population, but driving the regression.
#
# To compute a regression that is less sensitive to outliers, one must use a [robust model](https://en.wikipedia.org/wiki/Robust_statistics). This is done in seaborn using ``robust=True`` in the plotting functions, or in `statsmodels` by replacing the use of the OLS by a "Robust Linear Model", `statsmodels.formula.api.rlm`.

# # Testing for interactions
#
# Does `FSIQ` increase more with `PIQ` for people with dark hair than with light hair?

seaborn.lmplot(y='FSIQ', x='PIQ', hue='Hair', data=my_data_pandas)

# The plot above is made of two different fits. We need to formulate a single model that tests for a variance of slope across the population. This is done via an ["interaction"](http://statsmodels.sourceforge.net/devel/example_formulas.html#multiplicative-interactions).

from statsmodels.formula.api import ols
result = ols(formula='FSIQ ~ PIQ + Hair + PIQ * Hair', data=my_data_pandas).fit()
print(result.summary())

# # Take home messages
#
# * Hypothesis testing and p-value give you the **significance** of an effect / difference
#
# * **Formulas** (with categorical variables) enable you to express rich links in your data
#
# * **Visualizing** your data and simple model fits matters!
#
# * **Conditioning** (adding factors that can explain all or part of the variation) is an important modeling aspect that changes the interpretation.

StarcoderdataJupyterScriptsDedupFiltered