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

from pandas import read_csv, read_excel, crosstab, DataFrame
from statistics import mean, median
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import fisher_exact
from helpers.plots import labels, theme, show_corelations

theme()

# ## Load cleaned data

g = read_csv('data/cleaned_answers.csv')
g.head(2)

# ## Demographic statistics

scouts = g[g.is_polish_scout]
not_scouts = g[g.is_polish_scout == False]

h_mean = mean(scouts['total_score_without_images'])
nh_mean = mean(not_scouts['total_score_without_images'])
g_mean = mean(g['total_score_without_images'])
h_mean, nh_mean, g_mean

g.place_of_residence.value_counts()

# ### Number of valid participants

participants_n = len(g)
participants_n

# ### Number of women

g.is_women.describe()

len(g[g.is_women])

# Majority participants were women

f'{round((len(g[g.is_women])/len(g))*100, 2)}% of participants were women.'

# ### Number of scouts

g.is_polish_scout.describe()

len(g[g.is_polish_scout])

# Only about 1/3 of participants were scouts

f'{round((len(g[g.is_polish_scout])/len(g))*100, 2)}% of participants were scouts.'

# ### Age of participants

f'Mean age of all participants was: {round(mean(g.age), 2)} years.'

f'Median age of all participants was: {median(g.age)} years.'

# Mean and median age of participants are comparable so I expect little skew in the age distribution.

s = g.age.std()

# Assuming that age distribution is normal I would expect 95% of participants to be in this age group:

mean(g.age) - 2 * s, mean(g.age) + 2 * s

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
sns.distplot(g.age, kde=False, ax=ax1, hist_kws={'rwidth': 0.5})
sns.distplot(g.age, hist=False, ax=ax2, kde_kws={'bw': .5});

g.age.describe()

# ### Age distribution: control group vs Polish population

# Those data were used in my master thesis. The data are from 2017 as 2018 data were not available at the time of writing my master thesis. Now I have decided to use 2018 data as they become available.
#
# GUS 2017 and 2018 data comparison is available in the [GUS_17_and_18](GUS_17_and_18.ipynb) notebook.
#
# The data is prepared in [Demographics_data](Demographics_data.ipynb) notebook.

# **Data 2018 source:** [Central Statistical Office in Poland, *Ludność. Stan i struktura oraz ruch naturalny w przekroju terytorialnym w 2018 r. Stan w dniu 31 XII*](https://stat.gov.pl/obszary-tematyczne/ludnosc/ludnosc/ludnosc-stan-i-struktura-oraz-ruch-naturalny-w-przekroju-terytorialnym-w-2018-r-stan-w-dniu-31-xii,6,25.html)
#

gus_2018 = read_csv('data/gus_2018_cleaned.csv')
gus_2018.head()

gus_age = np.repeat(gus_2018.age, gus_2018.total)
gus_age

sns.kdeplot(g[~g.is_polish_scout].age, bw=.5, label="Grupa kontrolna", shade=True)
sns.kdeplot(gus_age, bw=.5, label="Populacja wg GUS", shade=True);

# Age distribution for all participants was not normal - Shapiro test

stats.shapiro(g.age)

# ### Age distribution: scouts vs control group

theme()
sns.kdeplot(g[g.is_polish_scout == False].age, bw=.5, label="Grupa kontrolna", shade=True)
sns.kdeplot(g[g.is_polish_scout].age, bw=.5, label="Harcerze", shade=True)
sns.rugplot(g['age'], color='.2', height=0.01)
labels(title='Rozkład wieku pomiędzy badanymi grupami', x='Wiek', y='Gęstość', legend=True)
theme()

# ks test sprawdza czy dystrybucje są takie same, im bliżej zera tym większa szansa, że są podobne, im mniejsze pvalue tym bardziej istotne statystycznej
#
# `ks_2samp` p-value is computed more precisely since SciPy 1.3.0 ([stats.ks_2samp changes log](https://docs.scipy.org/doc/scipy/reference/release.1.3.0.html#scipy-stats-improvements)) - see difference in p-value below:

# old
stats.ks_2samp(g[g.is_polish_scout].age, g[g.is_polish_scout == False].age, mode='asymp')

# new
stats.ks_2samp(g[g.is_polish_scout].age, g[g.is_polish_scout == False].age)

# z-score = miara efektu, czyli jak ważne

stats.skewtest(g[g.is_polish_scout].age)

stats.skewtest(g[g.is_polish_scout == False].age)

# `scipy.stats.skew`
#
# >For normally distributed data, the skewness should be about 0. A skewness value > 0 means that there is more weight in the left tail of the distribution. The function skewtest can be used to determine if the skewness value is close enough to 0, statistically speaking.
#
# > Returns: The skewness of values along an axis, returning 0 where all values are equal.
#
# Source: [SciPy docs](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.skew.html)

stats.skew(g[g.is_polish_scout].age)

stats.shapiro(g[g.is_polish_scout].age)

stats.shapiro(g[g.is_polish_scout == False].age)

# ### Sex distribution comparison

# Female to male ratio from GUS 2018 data

t = gus_2018[['males', 'females']].apply(sum)
gus_female_male_ratio = t['females'] / t['males']
print(f'GUS women to men ratio: {round(gus_female_male_ratio, 2)}')

# Female to male ratio according to ZHP 2018 data.
#
# The data is prepared in [Demographics_data > ZHP](Demographics_data.ipynb#ZHP) notebook.
#
# Data: http://hib.zhp.pl/zhp-w-liczbach/
#
# In my Master Thesis I used data from 2017 (`female_male_ratio_zhp = 1.44`) as 2018 was not available at that time.

zhp_in_numbers = read_csv('data/zhp_2018_cleaned.csv', index_col=0)

zhp_in_numbers.tail(1)

zhp_members_total = zhp_in_numbers["Razem:"]
zhp_members_total.head(4)

zhp_females = zhp_members_total.loc["Wędr K"] + zhp_members_total.loc["Wędr/Inst K"]
zhp_males = zhp_members_total.loc["Wędr M"] + zhp_members_total.loc["Wędr/Inst M"]
zhp_female_male_ratio = zhp_females / zhp_males
print(f'ZHP women to men ratio: {round(zhp_female_male_ratio, 2)}')

# Sex distribution: scouts vs control group

contingency_table = crosstab(g.sex, g.is_polish_scout)
contingency_table

is_scout_sex = contingency_table.to_dict()

# odds_ratio, pvalue
fisher_exact(contingency_table)

# pl: Kobieta -> eng: Woman
#
# pl: Mężczyzna -> eng: Man

ct = DataFrame(contingency_table.unstack(), columns=['count']).reset_index()
ct

# pl: Grupa kontrolna -> eng: Contol group
#
# pl: Harcerze -> eng: Scouts

dictionary_is_polish_scout = {
    False: 'Grupa kontrolna',
    True: 'Harcerze'
}
ct.is_polish_scout = ct.is_polish_scout.replace(dictionary_is_polish_scout)
ct.is_polish_scout

gr = sns.catplot(
    x="sex", y="count", hue="is_polish_scout", data=ct,
    kind="bar", palette="muted", legend=False, height=8, aspect=14 / 8
)
gr.ax.legend(loc='upper right', frameon=False)
labels(title='Rozkład płci pomiędzy badanymi grupami', x=' ', y='Liczna respondentów', legend=False);

is_scout_sex[True]['Kobieta']
no_scout_sex_ratio_w_to_m = (is_scout_sex[False]['Kobieta']) / (is_scout_sex[False]['Mężczyzna'])
scout_sex_ratio_w_to_m = (is_scout_sex[True]['Kobieta']) / (is_scout_sex[True]['Mężczyzna'])
print(
    f'Control group women to men ratio: {round(no_scout_sex_ratio_w_to_m, 2)}\n'
    f'Scouts women to men ratio: {round(scout_sex_ratio_w_to_m, 2)}'
)

gr = sns.FacetGrid(g, row="sex", col="is_polish_scout", margin_titles=True)
gr.map(sns.kdeplot, "age", color="steelblue", shade=True, bw=0.45);
axes = gr.axes.flatten()
axes[0].set_title("Control\ngroup")
axes[1].set_title("Scouts");

# ### Place distributipon

table = crosstab(g.place_of_residence, g.is_polish_scout)
table

t = DataFrame(table.unstack(), columns=['count']).reset_index()


def licz_procent(wiersz):
    return 100 * wiersz['count'] / len(g[g.is_polish_scout == wiersz.is_polish_scout])


t = t.assign(procent_grupy=t.apply(licz_procent, axis='columns'))
t

dict_place_of_residence = {
    'Miasto do 50 tys. mieszkańców': 'Miasto\ndo 50 tys.\nmieszkańców',
    'Miasto od 50 do 100 tys. mieszkańców': 'Miasto\nod 50 do 100 tys.\nmieszkańców',
    'Miasto powyżej 100 tys. mieszkańców': 'Miasto powyżej\n100 tys.\nmieszkańców'
}
t.place_of_residence = t.place_of_residence.replace(dict_place_of_residence)
t.place_of_residence

dict_is_polish_scout = {
    False: 'Grupa kontrolna',
    True: 'Harcerze'
}
t.is_polish_scout = t.is_polish_scout.replace(dict_is_polish_scout)
t.is_polish_scout

t.groupby(t.is_polish_scout).apply
t

# replace column name from 'is_polish_scout' to ' '
t.rename(columns={'is_polish_scout': 'group'}, inplace=True)

graph = sns.catplot(
    x="place_of_residence", y="procent_grupy", hue="group", data=t,
    kind="bar", palette="muted", aspect=3/2, height=10,
    legend=False
)
graph.ax.legend(loc='upper left', frameon=False)
labels(
    title='Miejsce zamieszkania w badanych grupach',
    x=' ', y='Procent respondentów w badanych grupach',
    legend=False
);

# +
counts = DataFrame(
    g[~g.is_polish_scout].place_of_residence.value_counts().reset_index()
).rename({'index': 'place', 'place_of_residence': 'total'}, axis='columns')
scout_counts = DataFrame(
    g[g.is_polish_scout].place_of_residence.value_counts().reset_index()
).rename({'index': 'place', 'place_of_residence': 'total'}, axis='columns')

control = {
    'miasto': sum(counts[counts.place.str.startswith('Miasto')].total),
    'wieś': counts[counts.place == 'Wieś'].total.iloc[0]
}
control
# -

scouts_miejsce = {
    'miasto': sum(scout_counts[scout_counts.place.str.startswith('Miasto')].total),
    'wieś': scout_counts[scout_counts.place == 'Wieś'].total.iloc[0]
}
scouts_miejsce

control_sum = control['miasto'] + control['wieś']
100 * control['miasto'] / control_sum

scouts_sum = scouts_miejsce['miasto'] + scouts_miejsce['wieś']
100 * scouts_miejsce['miasto'] / scouts_sum

total_gus = gus_2018[['city_total', 'rural_total']].apply(sum)
total_gus

# +
contingency_table = [
    # all populacja
    [control['miasto'], total_gus.loc['city_total']],   # miasto
    [control['wieś'], total_gus.loc['rural_total']],  # wieś
]

print('Miasto / wieś')
print('Kontrolna', control['miasto'] / control['wieś'])
print('Harcerze', scouts_miejsce['miasto'] / scouts_miejsce['wieś'])
print('GUS', total_gus.loc['city_total'] / total_gus.loc['rural_total'])

stats.fisher_exact(contingency_table), stats.chi2_contingency(contingency_table)
# -

(control['miasto'] / control['wieś']) / (total_gus.loc['city_total'] / total_gus.loc['rural_total'])

# ### Demographics corelations

show_corelations(g, only=[
    'age', 'is_women', 'size_of_residence_place', 'is_polish_scout'
])
stats.spearmanr(g.is_polish_scout, g.age);   # średniomaly efekt

stats.spearmanr(g.is_polish_scout, g.size_of_residence_place)

stats.spearmanr(g.is_polish_scout, g.is_women)

StarcoderdataJupyterScriptsDedupFiltered