Udacity AB testing Online Experiment Design and Analysis - Final Project in Python

import plotly.express as px
import pandas as pd

df = pd.DataFrame([
    dict(Diversion="cookie diversion", Start='2021-01-01', End='2021-01-02', Phase="unique user cookie course overview page"),
    dict(Diversion="cookie diversion", Start='2021-01-04', End='2021-01-05', Phase="unique user cookie course overview page"),
    dict(Diversion="cookie diversion", Start='2021-01-05', End='2021-01-06', Phase="unique user cookie course overview page"),
    dict(Diversion="cookie diversion", Start='2021-01-06', End='2021-01-07', Phase="unique user cookie course overview page"),
    dict(Diversion="both diversions", Start='2021-01-06', End='2021-01-07', Phase="'Start free trial' button click"),
    dict(Diversion="user-id diversion", Start='2021-01-05', End='2021-02-20', Phase="user-id login"),
    dict(Diversion="user-id diversion", Start='2021-01-07', End='2021-01-21', Phase="checkout & free trial"),
    dict(Diversion="user-id diversion", Start='2021-01-21', End='2021-02-20', Phase="enrolled past 14 days"),
    dict(Diversion="user-id diversion", Start='2021-02-20', End='2021-02-28', Phase="unenrolled"),
])

# fig = px.timeline(df, x_start="Start", x_end="End", y="Phase", color="Diversion",
#                   color_discrete_sequence=px.colors.qualitative.Vivid,
#                  title='Exemplary phases of the Udacity user signup process' )
# fig.show()

df.head()

# Baseline Values

# udacity baseline google sheet 
google_sheet_url = 'https://docs.google.com/spreadsheets/d/1MYNUtC47Pg8hdoCjOXaHqF-thheGpUshrFA21BAJnNc/export?format=csv&gid=0'

import pandas as pd
df=pd.read_csv(google_sheet_url, names=['observation', 'baseline'], encoding="utf-8")

#df.iloc[2,0] = 'Enrollments'
df['shortname'] = ['cookies',
              'clicks',
              'enrollments',
              'p_click_through',
              'p_enroll_given_click',
               'p_payment_given_enroll',
               'p_payment_given_click'
              ]
df['metric'] = ['',
              '',
              '',
              '',
              'gross_conversion',
               'retention',
               'net_conversion',
              ]
df['dmin'] = [3000,240,50,0.01,0.01,0.01,0.0075]
df['potent_metric'] = [False,False,False,False,True,True,True]
df = df.set_index('shortname')
#df['n'] = 0

df

df['reference_value'] = ['cookie', 'click', 'enroll', 'cookie', 'click', 'enroll', 'click']
df.loc[df['reference_value'] == 'cookie', 'n'] = 5000
df.loc[df['reference_value'] == 'click', 'n'] = df.loc['cookies','n'] * df.loc['p_click_through','baseline'] 
df.loc[df['reference_value'] == 'enroll', 'n'] = df.loc['clicks','n'] * df.loc['p_enroll_given_click','baseline']

import numpy as np

def StandardError(p,n):
    stdev = np.sqrt(p*(1-p)/n)
    return stdev

df['SE'] = df[df['potent_metric']==True].apply( lambda row: StandardError(row['baseline'], row['n']), axis=1)

df

from scipy.stats import norm

#Inputs: required alpha value
#Returns: z-score for given alpha
def z_score(alpha):
    return norm.ppf(alpha)

# Inputs p-baseline conversion rate which is our estimated p and dmin detectable change
# Returns
def sds(p,d):
    sd1=np.sqrt(2*p*(1-p))
    sd2=np.sqrt(p*(1-p)+(p+d)*(1-(p+d)))
    sds=[sd1,sd2]
    return sds

# In: sd1-sd for the baseline, 
#     sd2-sd for the expected change,
#     alpha,
#     beta,
#     dmin,
#     p-baseline estimate p
# Returns: the minimum sample size required per group according to metric denominator
def sampSize(sds,alpha,beta,d):
    size=(z_score(1-alpha/2)*sds[0]+z_score(1-beta)*sds[1])**2/d**2
    return size

# bonferroni correction is needed to control the family-wise error rate as both metrics are not independent, 
# and both metrics share the Clicks number therefore the chance for a false positive is certainly higher. 
# it is assumed that one significant test of a single metric is sufficient to reject the Null Hypothesis
alpha = 0.05
numtests = 2 # retention metric excluded
bonferroni_individual_alpha = alpha / numtests

beta = 0.2

df['samplesize']=sampSize(sds(df.baseline[3:7],df.dmin[3:7]),bonferroni_individual_alpha, beta, df.dmin[3:7]) # only probab.
df

df.loc[df.metric == 'gross_conversion', 'required_pageviews'] = df.loc[df.metric == 'gross_conversion', 'samplesize'] / df.loc['p_click_through', 'baseline'] * 2 # 2 groups A + B
df.loc[df.metric == 'gross_conversion', 'required_pageviews'] = df.loc[df.metric == 'gross_conversion', 'required_pageviews'].astype(int)

df.loc[df.metric == 'retention', 'required_pageviews'] = df.loc[df.metric == 'retention', 'samplesize'] / df.loc['p_click_through', 'baseline'] / df.loc[df.metric == 'gross_conversion', 'baseline'][0] * 2 # 2 groups A + B
df.loc[df.metric == 'retention', 'required_pageviews'] = df.loc[df.metric == 'retention', 'required_pageviews'].astype(int)

df.loc[df.metric == 'net_conversion', 'required_pageviews'] = df.loc[df.metric == 'net_conversion', 'samplesize'] / df.loc['p_click_through', 'baseline'] * 2 # 2 groups A + B
df.loc[df.metric == 'net_conversion', 'required_pageviews'] = df.loc[df.metric == 'net_conversion', 'required_pageviews'].astype(int)

df['duration'] = df.apply( lambda row: row.required_pageviews / df.loc['cookies', 'baseline'], axis=1)

determine_duration = df[df.potent_metric==True & (df.metric!='retention')][['metric', 'required_pageviews', 'duration']]
determine_duration

print(f'the test has to run at least {int(np.ceil(max(determine_duration.duration)))} days assuming constant website traffic')

# https://stackoverflow.com/questions/19611729/
def load_from_gspreadsheet(sheet_name, key):
    url = 'https://docs.google.com/spreadsheets/d/{key}/gviz/tq?tqx=out:csv&sheet={sheet_name}&headers=1'.format(
        key=key, sheet_name=sheet_name.replace(' ', '%20'))

    #log.info('Loading google spreadsheet from {}'.format(url))

    df = pd.read_csv(url)
    return df.drop([col for col in df.columns if col.startswith('Unnamed')], axis=1)

cont = load_from_gspreadsheet('Control', '1Mu5u9GrybDdska-ljPXyBjTpdZIUev_6i7t4LRDfXM8')
cont.head()

exp = load_from_gspreadsheet('Experiment', gsheet_results_key)
exp.head()

# proportion checking

# H0: p = 0.5 # equal click proportions 

cont.Pageviews.sum()

exp.Pageviews.sum()

p = 0.5

stdev = np.sqrt(p*(1-p)/(cont.Pageviews.sum() + exp.Pageviews.sum()))
stdev

# zscore
alpha = 0.05
z= 1.96

# margin of error
margin_of_error = stdev * z
margin_of_error 

# confidence interval construction

min_ci = p - margin_of_error
max_ci = p + margin_of_error

print('confidence interval:')
print(round(min_ci,4))
print(p)
print(round(max_ci,4))

cont_cookie_prop = round(cont.Pageviews.sum() / (cont.Pageviews.sum() + exp.Pageviews.sum()),4)
if(cont_cookie_prop > min_ci and cont_cookie_prop < max_ci):
    print(f'sanity check passed: {cont_cookie_prop} is inside the confidence interval')

# proportion checking

# H0: p = 0.5 # equal click proportions 

p = 0.5

stdev = np.sqrt(p*(1-p)/(cont.Clicks.sum() + exp.Clicks.sum()))
stdev

# zscore
alpha = 0.05
z= 1.96

# margin of error
margin_of_error = stdev * z

# confidence interval construction

min_ci = p - margin_of_error
max_ci = p + margin_of_error

print('confidence interval:')
print(round(min_ci,4))
print(p)
print(round(max_ci,4))

cont_click_prop = cont.Clicks.sum() / (cont.Clicks.sum() + exp.Clicks.sum())
if(cont_click_prop > min_ci and cont_click_prop < max_ci):
    print('sanity check passed - '+str(cont_click_prop)+' is inside the confidence interval')

# checking differences in probability-type-numbers

# H0: d = 0

# work with Pandas Series
ctp_cont = cont.Clicks.sum() / cont.Pageviews.sum()
ctp_exp = exp.Clicks.sum() / exp.Pageviews.sum()

d_hat = ctp_exp - ctp_cont 
d_hat

# gross conversion across both groups
p_pooled_hat = (cont.Clicks.sum() + exp.Clicks.sum()) / (cont.Pageviews.sum() + exp.Pageviews.sum())
round(p_pooled_hat,4)

# sample standard deviations

# SE_pooled
SE_pool = np.sqrt(p_pooled_hat*(1-p_pooled_hat) * (1/cont.Pageviews.sum() + 1/exp.Pageviews.sum()) )
SE_pool

import scipy.stats as st

# default one tail
# two tails needed
p_value = 0.975  # (1-.05/2)
z = st.norm.ppf(p_value)
z

# margin of error
margin_of_error = SE_pool * z
margin_of_error 

# confidence interval construction
# The expected difference d between the two groups is zero. 


min_ci = 0 - margin_of_error
max_ci = 0 + margin_of_error

print('confidence interval:')
print(round(min_ci,4))
print(0)
print(round(max_ci,4))

if(d_hat > min_ci and d_hat < max_ci):
    print(f'sanity check passed: d_hat {d_hat} is inside the confidence interval')

# disregard cookies which have not completed 14 days trial phase
cont_non_trial = cont.loc[~cont.Enrollments.isnull()]
exp_non_trial = exp.loc[~exp.Enrollments.isnull()]

print(len(exp))
print(len(exp_non_trial))

# gross conversions per group
cont_gc = cont_non_trial.Enrollments.sum() / cont_non_trial.Clicks.sum()
exp_gc = exp_non_trial.Enrollments.sum() / exp_non_trial.Clicks.sum()

# gross conversion across both groups
pooled_gc_hat = (cont_non_trial.Enrollments.sum() + exp_non_trial.Enrollments.sum()) / (cont_non_trial.Clicks.sum() + exp_non_trial.Clicks.sum())
pooled_gc_hat

# SE_pooled
SE_gc_pool = np.sqrt(pooled_gc_hat*(1-pooled_gc_hat) * (1/cont_non_trial.Clicks.sum() + 1/exp_non_trial.Clicks.sum()) )
SE_gc_pool

d_hat = exp_gc - cont_gc
d_hat

# bonferroni correction
alpha = 0.05
numtests = 2 # retention metric excluded
bonferroni_individual_alpha = alpha / numtests

# margin of error
p_value = 1 - (bonferroni_individual_alpha / 2) # two tails
#print(p_value)

z = st.norm.ppf(p_value)
#print(z)
m = SE_gc_pool * z
m

# confidence interval construction

min_ci = d_hat - m
max_ci = d_hat + m
print('confidence interval:')
print(round(min_ci,4))
print(round(d_hat,4))
print(round(max_ci,4))

if(0 >= min_ci and 0 <= max_ci):
    print('Null hypothesis cannot be rejected as zero is inside the confidence interval')
else:
    print('Null hypothesis can be rejected as zero is outside the confidence interval')
    print('as per experiment design there is a significant reason to doubt the Null hypothesis and') 
    print('conclude that there is indeed a change in the data between control and treatment groups')

# Gross Conversion
dmin_gc = df.loc[df['metric'] == 'gross_conversion', 'dmin'][0]

if (abs(d_hat) > dmin_gc):
    print(f'This observed change {round(abs(d_hat),4)} is practically signficant, as larger than {dmin_gc} and considered worth the implementation')
else:
    print(f'This observed change {round(abs(d_hat),4)} is not practically signficant, as smaller than {dmin_gc} and not considered worth the implementation')

# disregard cookies which have not completed 14 days trial phase
#cont_non_trial = cont.loc[~cont.Enrollments.isnull()]
#exp_non_trial = exp.loc[~exp.Enrollments.isnull()]

# net conversions per group
cont_nc = cont_non_trial.Payments.sum() / cont_non_trial.Clicks.sum()
exp_nc = exp_non_trial.Payments.sum() / exp_non_trial.Clicks.sum()

# gross conversion across both groups
pooled_nc_hat = (cont_non_trial.Payments.sum() + exp_non_trial.Payments.sum()) / (cont_non_trial.Clicks.sum() + exp_non_trial.Clicks.sum())
pooled_nc_hat

# sample standard deviations
#stdev_cont = (cont_gc*(1-cont_gc))**0.5
#stdev_exp = (exp_gc*(1-exp_gc))**0.5

# SE_pooled
nc_SE_pool = np.sqrt(pooled_nc_hat*(1-pooled_nc_hat) * (1/cont_non_trial.Clicks.sum() + 1/exp_non_trial.Clicks.sum()) )
nc_SE_pool

nc_d_hat = exp_nc - cont_nc
nc_d_hat

# bonferroni correction
alpha = 0.05
numtests = 2 # retention metric excluded
bonferroni_individual_alpha = alpha / numtests

# margin of error
p_value = 1 - (bonferroni_individual_alpha / 2) # two tails
#p_value = 0.975
print(p_value)
z = st.norm.ppf(p_value)
print(z)
m = nc_SE_pool * z
print(m)

# confidence interval construction

nc_min_ci = nc_d_hat - m
nc_max_ci = nc_d_hat + m
print(round(nc_min_ci,4))
print(round(nc_d_hat, 4))
print(round(nc_max_ci,4)) 

if(0 >= nc_min_ci and 0 <= nc_max_ci):
    print('The Null hypothesis cannot be rejected as zero is inside the confidence interval,')
    print('as per experiment design there is no significant reason to doubt the Null hypothesis and conclude that there is no difference')
else:
    print('Null hypothesis can be rejected as zero is outside the confidence interval')

# prepare diff series

cont_gc_daily = cont_non_trial.Enrollments / cont_non_trial.Clicks
exp_gc_daily = exp_non_trial.Enrollments / exp_non_trial.Clicks
diff_hat_gc_daily = exp_gc_daily - cont_gc_daily

cont_nc_daily = cont_non_trial.Payments / cont_non_trial.Clicks
exp_nc_daily = exp_non_trial.Payments / exp_non_trial.Clicks
diff_hat_nc_daily = exp_nc_daily - cont_nc_daily

# first approach: built-in statsmodels function sign_test

from statsmodels.sandbox.descstats import sign_test

# gross conversion
print(f'gc sign test p-value: {sign_test(diff_hat_gc_daily, 0)[1]}')

# net conversion
print(f'nc sign test p-value: {sign_test(diff_hat_nc_daily, 0)[1]}')

# second approach: calculate manually with which we are able to state on how many days there is a significant change
#                  instead of just receiving the result itself as in first approach

def signtest_diffseries(diff_series, alternative = "two_sided"):
    from scipy.stats import binom
    pos_diff_count = len([k for k in diff_series if k > 0])
    timeunits = len(diff_series)
    median_count = np.median(np.arange(0, timeunits + 1))
    
    if (alternative == "two_sided") and (pos_diff_count >= median_count):
      p_value = round(2 * (1 - binom.cdf((pos_diff_count - 1), timeunits, 0.5)),4)
    elif pos_diff_count < median_count:
      p_value = round(2 * binom.cdf(pos_diff_count, timeunits, 0.5),4)
    
    print(f"Number of positive time units: {pos_diff_count} out of a sample size of {timeunits} time units")
    if(p_value < 0.05):
        print(f"The p-value of {p_value} does indicate a significant change in the data. The Null hypothesis H0: d=0 has to be rejected")
    else:
        print(f"The p-value of {p_value} does not indicate a significant change in the data.")

signtest_diffseries(diff_hat_gc_daily)

signtest_diffseries(diff_hat_nc_daily)