Contents

path2 = Path(r'./../data/')


country = ['Brazil', 'br']
country = ['all countries', 'all']


load_results = True  # load file with results if exists
if not load_results:
    # Warning: the bootstrap process takes hours!
    from bootstrap_stat import bootstrap_stat as bp

tqdm.pandas(desc='Boot', leave=False)
pd.set_option('display.precision', 3)
plt.rcParams.update({'font.size': 14, 'xtick.labelsize': 12,
                     'ytick.labelsize': 12})
sns.set_style('whitegrid', rc={'xtick.bottom': True, 'xtick.top': True,
                               'ytick.left': True,
              'ytick.right': True, 'xtick.direction': 'in',
                               'ytick.direction': 'in'})
Slice = pd.IndexSlice

# Set of colors that is unambiguous both to colorblinds and non-colorblinds
# https://www.cta-observatory.org/wp-content/uploads/2020/10/CTA_ColourBlindness_BestPractices-1.pdf
# https://jfly.uni-koeln.de/color/index.html
#matplotlib.style.use('seaborn-colorblind')
matplotlib.style.use('tableau-colorblind10')
colors = sns.color_palette()
display(colors)

colorbase = 'k'
color2019 = colors[0]
color2020 = colors[1]
colorn = colors[5]
colorvar = colors[4]
colorind = colors[6]

# Statistics
alpha = 0.05
two_tailed = True
# confidence interval
estimate = np.mean
boot = 100000  # number of bootstrap samples
seed = 0

# other variables:
years = ['2019', '2020']
quarters = ['First', 'Second', 'Third', 'Fourth']
months = ['January', 'February', 'March', 'April', 'May', 'June', 'July',
          'August', 'September', 'October', 'November', 'December']
freqs = ['d', 'w', 'm', 'q']
freqs_s = ['day', 'week', 'month', 'quarter']
# feature (dependent variable)
variables, units = ['distance', 'duration'], ['km', 'min']
vars_d = [v + '_d' for v in variables]
# period for analysis
freq, freq_s = freqs[1], freqs_s[1]

%watermark
%watermark --iversions

dfs, dfs_no0, nathletes, nruns = dict(), dict(), dict(), dict()
print('Country: {}'. format(country[0]))
for year in years:
    for f, frq in enumerate(freqs):
        # load dataset and select country from start to speedup process
        fname = path2 / 'run_ww_{}_{}.parquet'.format(year, frq)
        df = pd.read_parquet(fname, columns=['datetime', 'athlete', 'distance',
                                             'duration', 'country'])
        if country[1] != 'all':
            df = df[df['country'] == country[0]]
        df = df.set_index('datetime', drop=True)
        #df.index = df.index.normalize()
        # correct dates
        if df.index[0].is_leap_year and frq == 'd':
            # if leap year, remove Feb 29, average data of Feb 28 and 29
            feb28 = pd.Timestamp(unit='D', year=df.index[0].year, month=2, day=28)
            feb29 = pd.Timestamp(unit='D', year=df.index[0].year, month=2, day=29)
            nbs = df.select_dtypes(include='number').columns.to_list()
            df.loc[df.index==feb28, nbs] = (df.loc[df.index==feb28, nbs].values +
                                            df.loc[df.index==feb29, nbs].values)/2
            df = df.drop(index=feb29)  
        
        display(df.iloc[[0, -1]].style.set_caption(
                'Year: {}, Frequency: {} (first and last rows)'.format(year, freqs_s[f])))
        # create dictionary will all the data (including no runs data)
        dfs[year, frq] = df
        # create dictionary will all the data (excluding no runs data)
        dfs_no0[year, frq] = df[df[variables[0]] > 0]
        # number of athletes and running activities
        nathletes[year, frq] = dfs_no0[year, frq].groupby(level='datetime').size()
        x = dfs_no0[year, frq].groupby('athlete').size()
        nruns[year, frq] = x[x > 0]  # removes athletes with zero runs

dfs['2019', 'd'].info(memory_usage='deep')

dfs['2019', freq]

dfs['2020', freq]

def correct_date(y, freq):
    """
    Correct dates and adjust data of last days to avoid last week with < 7 days
     or of the last day of February if leap year.
    """
    # numerical columns
    cols_num = y.select_dtypes(include='number').columns.to_list()
    ndays_lastweek = 7 + y.index[-1].dayofyear % 7
    if freq == '7d' and ndays_lastweek > 7:
        ts = pd.Timestamp(y.index[-1].date() - pd.to_timedelta(ndays_lastweek-7, unit='D'))
        y.index = y.index.where(y.index <= ts, ts)    
        grouper = pd.Grouper(axis=0, freq=freq)
        y = y.groupby(grouper).mean()
        # correct the divisor if the last week doesn't have 7 days
        y.loc[y.index[-1], cols_num] = y.loc[y.index[-1], cols_num] * (7 / ndays_lastweek)
    elif freq == 'd' and y.index.unique()[59].day == 29:  # leap year
        grouper = pd.Grouper(axis=0, freq=freq)
        y = y.groupby(grouper).mean()
        # February 28 as the average February 28 and 29
        y.loc[y.index[58], cols_num] = y.loc[y.index[58:60], cols_num].mean(axis=0)  
        y = y.drop(index=y.index[59])
    elif freq == 'm' and y.index.unique()[59].day == 29:  # just resample data
        grouper = pd.Grouper(axis=0, freq=freq)
        y = y.groupby(grouper).mean()
    return y

c19idx = correct_date(c19idx, '7d' if freq == 'w' else freq)
c19idx

for year in years:
    print('Year: {}'.format(year))
    print('  Number of days:', dfs[year, 'd'].index.value_counts().size)
    print('  Number of athletes:', nruns[year, 'd'].size)
    print('  Number of activities:', dfs[year, 'd'].shape[0])
    print('  Number of actual running activities:', nathletes[year, 'd'].sum())

nathlete = nruns['2019', 'd'].size
nactivity = nathletes['2019', 'd'].sum() + nathletes['2020', 'd'].sum()

print('Total number of actual running activities in 2019 and 2020:', nactivity)

print('Median and interquartile range of number of days the athletes ran each year:')
for year in years:
    print('In {}: {:.0f} ({:.0f}, {:.0f}) days'.
          format(year, *np.percentile(nruns[year, 'd'], [50, 25, 75]).round()))

fig, axs = plt.subplots(3, 1, squeeze=True, sharex=True, figsize=(12, 7))
    
for i, year in enumerate(years):
    y = nathletes[year, 'd'].values
    axs[0].plot(np.linspace(0, 365, 366), np.append(y, y[-1]),
                color=colors[i], lw=2, label=year, drawstyle='steps-post')
    y = nathletes[year, 'w'].values
    axs[1].plot(np.linspace(0, 365, 53), np.append(y, y[-1]),
                color=colors[i], lw=3, label=year, drawstyle='steps-post')
    x = nathletes[year, 'm'].index.day_of_year.values
    y = nathletes[year, 'm'].values
    axs[2].plot(np.append([0], x), np.append(y, y[-1]),
                color=colors[i], lw=3, label=year, drawstyle='steps-post')

axs[0].set_xlim([0, 365])
axs[0].legend(['2019', '2020'])
axs[0].set_ylabel('Per day')
axs[1].set_ylabel('Per week')
axs[2].set_ylabel('Per month')
axs[2].set_xlabel('Month')
axs[2].set_xticks(np.append([0], x[:-1]))
axs[2].set_xticklabels(months, rotation=0, ha='center')
offset = matplotlib.transforms.ScaledTranslation(34/72, 0, fig.dpi_scale_trans)
for label in axs[2].get_xticklabels():
    label.set_transform(label.get_transform() + offset)
for ax in axs:
    ax.tick_params(axis='both', which='major', labelsize=13)
plt.suptitle('Number of athletes running per period (out of a total of {:,d})'
             .format(nruns['2019', 'd'].size), fontsize=18)
plt.tight_layout(h_pad=.5)
plt.show()

nruns['2019', 'd'].describe()

365 / nruns['2019', 'd'].describe()

nruns['2020', 'd'].describe()

365 / nruns['2020', 'd'].describe()

print('Parameters:')
print(' years =', years)
print(' frequency =', freq)
print(' variables =', variables)
print(' estimate =', estimate)

grouper = pd.Grouper(axis=0, freq='y', sort=True)
for year in years:
    print(year)
    display(dfs[year, freq][['distance', 'duration']
                           ].groupby(grouper).describe().stack(level=0))
    display(dfs[year, freq][['distance', 'duration']
                           ].groupby(grouper).sum()) 

grouper = pd.Grouper(axis=0, freq='y', sort=True)
for year in years:
    print(year)
    display(dfs_no0[year, freq][['distance', 'duration']
                               ].groupby(grouper).describe().stack(level=0))    
    display(dfs_no0[year, freq][['distance', 'duration']
                               ].groupby(grouper).sum())

for year in years:
    y = dfs[year, freq][['duration', 'distance']].groupby(dfs[year, freq].index).agg('sum')
    y['pace'] = y['duration'] / y['distance']
    y['nathletes'] = nathletes[year, freq]
    y.loc[:, ['duration', 'distance']] = y.loc[:, ['duration', 'distance']
                                              ].div(y['nathletes'], axis=0)
    offset = 30 if freq == 'm' else 0
    y.index = dfs['2019', freq].index.dayofyear.unique() - offset 
    y.index.rename('Period', inplace=True)
    last = y.iloc[-1]
    last.name = 366
    y = pd.concat([y, last.to_frame().transpose()])  # for plotting with steps
    y = y[['nathletes', 'duration', 'distance', 'pace']]
    y.rename(columns={'nathletes': 'N athletes', 'distance': 'Distance (km)',
                      'duration': 'Duration (min)', 'pace': 'Pace (min/km)'}, inplace=True)
    if year == '2019':
        data = y.copy(deep=True)
    else:
        data = pd.concat([data, y], axis=1, keys=years, names=['year', 'var'])
        
display_df(data, subset=data.columns)

fig, axs = plt.subplots(4, 1, sharex=True, figsize=(12, 8))
cols = data['2019'].columns
precision = [0, 1, 1, 2]
for a, ax in enumerate(axs):
    for y, year in enumerate(years):
        ax.plot(np.linspace(0, 365, data[year].shape[0]), data[year][cols[a]],
                    color=colors[y], lw=3, label=year, drawstyle='steps-post')
        msd = '{:g}±{:g}'.format(np.round(data[year][cols[a]].iloc[:-1].mean(), precision[a]),
                                 np.round(data[year][cols[a]].iloc[:-1].std(), precision[a]))
        ax.text(.98, .85-y/7, msd, ha='right', c=colors[y],
                fontsize=14, transform=ax.transAxes)
    ax.margins(x=0, y=0.3)
    ax.tick_params(axis='both', which='major', labelsize=13)
    ax.set_ylabel(cols[a])

axs[0].text(.97, 1.03, 'Mean±SD', ha='right', c='k',
            fontsize=14, transform=axs[0].transAxes)
axs[0].set_xlim([0, 365])
axs[0].legend(['2019', '2020'], loc='lower left', frameon=True, framealpha=.5, fontsize=14)
axs[3].set_xlabel('Month')
x = nathletes['2019', 'm'].index.day_of_year.values
axs[3].set_xticks(np.append([0], x[:-1]))
axs[3].set_xticklabels(months, rotation=0, ha='center')
offset = matplotlib.transforms.ScaledTranslation(34/72, 0, fig.dpi_scale_trans)
for label in axs[3].get_xticklabels():
    label.set_transform(label.get_transform() + offset)

plt.suptitle('Number of athletes and running volume per athlete at each {} (from only actual runs)'
             .format(freq_s), fontsize=18)
fig.align_ylabels(axs)
plt.tight_layout(h_pad=.05)
plt.show()

nathletesd = dict()
for frq in freqs:
    # 2019 total number of athletes
    Na = nruns['2019', 'd'].size
    # 2019 average number of athletes per period (day, week or month)
    #Na = estimate(nathletes['2019', frq])
    # 2019 average running volume per period (day, week or month)    
    Vr = estimate(dfs['2019', frq][variables].values, axis=0)
    
    dfs['2020', frq][vars_d] = (100 * (dfs['2020', frq][variables].values -
                                       dfs['2019', frq][variables].values) / Vr)

    nathletesd[frq] = pd.Series(100 * ((nathletes['2020', frq].values -
                                        nathletes['2019', frq].values
                                       ) / estimate(nathletes['2019', frq])),
                                index=nathletes['2020', frq].index)
    
display(dfs['2020', freq])

print('Relative difference in total number of athletes between years')
100 * (nruns['2020', 'd'].size - nruns['2019', 'd'].size) / nruns['2019', 'd'].size

print('Relative difference in number of athletes per {} between years'.format(freq_s))
display(nathletesd[freq].describe().to_frame().T)    
display_df(nathletesd[freq].groupby('datetime').agg(['mean']))

grouper = pd.Grouper(axis=0, freq='y', sort=True)
display(dfs['2020', freq][vars_d].groupby(grouper).describe().stack(level=0))    

effect_sizes = pd.DataFrame({'d': [.01, .2, .5, .8, 1.2, 2], 'effect size':
                             ['Very small', 'Small', 'Medium', 'Large', 'Very large', 'Huge']})
effect_sizes

n_periods = dfs['2020', freq].index.unique().size
alphaB = alpha / n_periods
alphaS = 1-(1-alpha)**(1/n_periods)
print('Number of comparisons: {}'.format(n_periods))
print('Alpha level: {:.2f}'.format(alpha))
print('Alpha level with Bonferroni correction: {:.5f}'.format(alphaB))
print('Alpha level with Sidak correction: {:.5f}'.format(alphaS))

data = pd.DataFrame([])
frq = 'd'
for year in years:
    y = dfs[year, frq][['athlete', 'distance', 'duration']].groupby('athlete').sum() * 7 /365
    y['pace'] = y['duration'] / y['distance']
    for col in y.columns:
        y.rename(columns={col: col+year}, inplace=True)
    data = pd.concat([data, y], axis=1)
    data = pd.concat([data, nruns[year, 'd'].to_frame(name='nruns'+year) * 7 / 365], axis=1) 
# replace nans with zeros; merging 2019 and 2020 creates nans in nruns2020 due to nonexistent data  
data['nruns2020'] = data['nruns2020'].fillna(0)
# replace nans with zeros; merging 2019 and 2020 creates nans in pace2020 due to no runs,
#  not replacing nans will use only actual runs to compute pace (makes sense), but it will
#  result in 2029 and 2020 having different number of athletes when computing averages.
#data['pace2020'] = data['pace2020'].fillna(0)
# relative differences
# normalize by global values
Vr = estimate(data[['distance2019', 'duration2019', 'pace2019', 'nruns2019']].values, axis=0)
# normalize by individual values
#Vr = data[['distance2019', 'duration2019', 'pace2019', 'nruns2019']].values

data[['distanced','durationd','paced','nrunsd']] = (100*(data[['distance2020','duration2020',
                                                               'pace2020','nruns2020']].values -
                                                         data[['distance2019','duration2019',
                                                               'pace2019','nruns2019']].values) / Vr)
display(data)

stat = ['count', 'mean', 'std', cohensd]
ys = data.agg(stat, axis=0).T
for var in tqdm(data.columns):
    ys.loc[var, ['ci_inf', 'ci_sup', 'pvalue']] = ci_asl(data[var],
                                                         alpha2=alpha/2, boot=10000)
    if var[-1] != 'd':
        ys.loc[var, ['cohensd', 'pvalue']] = ['-', '-']
display(ys)

data = pd.DataFrame([])
frq = 'd'
for year in years:
    y = dfs[year, frq][['athlete', 'distance', 'duration']].groupby('athlete').sum() * 7 /365
    y['pace'] = y['duration'] / y['distance']
    for col in y.columns:
        y.rename(columns={col: col+year}, inplace=True)
    data = pd.concat([data, y], axis=1)
    data = pd.concat([data, nruns[year, 'd'].to_frame(name='nruns'+year) * 7 / 365], axis=1) 
# replace nans with zeros; merging 2019 and 2020 creates nans in nruns2020 due to nonexistent data  
data['nruns2020'] = data['nruns2020'].fillna(0)
# replace nans with zeros; merging 2019 and 2020 creates nans in pace2020 due to no runs,
#  not replacing nans will use only actual runs to compute pace (makes sense), but it will
#  result in 2029 and 2020 having different number of athletes when computing averages.
data['pace2020'] = data['pace2020'].fillna(0)
# relative differences
# normalize by global values
Vr = estimate(data[['distance2019', 'duration2019', 'pace2019', 'nruns2019']].values, axis=0)
# normalize by individual values
#Vr = data[['distance2019', 'duration2019', 'pace2019', 'nruns2019']].values

data[['distanced','durationd','paced','nrunsd']] = (100*(data[['distance2020','duration2020',
                                                               'pace2020','nruns2020']].values -
                                                         data[['distance2019','duration2019',
                                                               'pace2019','nruns2019']].values) / Vr)

stat = ['count', 'mean', 'std', cohensd]
ys = data.agg(stat, axis=0).T
for var in tqdm(data.columns):
    ys.loc[var, ['ci_inf', 'ci_sup', 'pvalue']] = ci_asl(data[var],
                                                         alpha2=alpha/2, boot=10000)
    if var[-1] != 'd':
        ys.loc[var, ['cohensd', 'pvalue']] = ['-', '-']
display(ys)

for year in years:
    print('Year: {}'.format(year))
    print('  Number of days:', dfs[year, 'd'].index.value_counts().size)
    print('  Number of athletes with a least one session:', nruns[year, 'd'].size)
    print('  Number of activities:', dfs[year, 'd'].shape[0])
    print('  Number of actual running activities:', nathletes[year, 'd'].sum())

nathlete = nruns['2019', 'd'].size
nactivity = nathletes['2019', 'd'].sum() + nathletes['2020', 'd'].sum()

print('Total number of actual running activities in 2019 and 2020:', nactivity)
nd = 100 * (nruns['2020', 'd'].size - nruns['2019', 'd'].size) / nruns['2019', 'd'].size
print('Relative difference in number of athletes [%]: {}'.format(nd))

nwd = 100 * (nathletes['2020', 'w'].values - nathletes['2019', 'w'].values
            ) / nathletes['2019', 'w'].mean()
data = pd.DataFrame(data=np.c_[nathletes['2019', 'w'].values, nathletes['2020', 'w'].values, nwd],
                    columns=['nathletesw2019', 'nathletesw2020', 'nathleteswd'])
#display(data)

stat = ['count', 'mean', 'std', cohensd]
ys = data.agg(stat, axis=0).T
for var in tqdm(data.columns):
    ys.loc[var, ['ci_inf', 'ci_sup', 'pvalue']] = ci_asl(data[var],
                                                         alpha2=alpha/2, boot=1000)
    if var[-1] != 'd':
        ys.loc[var, ['cohensd', 'pvalue']] = ['-', '-']
display(ys)