<jupyter_start><jupyter_text>Using Random EMA to check End-of-Day: Exploratory Data Analysis- This notebook is dedicated to understanding End-of-Day EMA using Random EMA- For every Random EMA where the response is 'Yes', check to see    + What is the fraction where the user clicked correct hour    + What is the fraction where the user clicked plus/minus one hour- We see that for _aggregated data_, the fraction where Random EMA and EOD agree is:    + 0.413 for current hour only    + 0.712 for plus/minus hour- If we assume normal and use midpoint of hour then we can convert to an estimated variance of a normal:    + \Phi(30/\sigma) - \Phi(-30/\sigma) = 2\Phi(30/sigma) - 1 = 0.413 -> sigma = 30/ \Phi^{-1} ( (0.413+1)/2 ) = 55 minutes    + \sigma = 60/ \Phi^{-1} ( (0.712+1)/2 ) = 56.47 minutes         * Here we use 60 since it was pm one hour so 90 seemed excessive    + So both suggest a prior with variance 55 minutes (or Inverse Wishart with mean 55<jupyter_code>import pandas as pd
import numpy as np
import datetime as datetime
import os
os.getcwd()
dir = "../final-data"
# different hour windows in eod_ema
keys = ['8to9', '9to10', '10to11', '11to12','12to13','13to14','14to15','15to16','16to17','17to18','18to19','19to20']

random_accptresponse = ['1 - 19 Minutes', '20 - 39 Minutes', '40 - 59 Minutes', 
                    '60 - 79 Minutes', '80 - 100 Minutes']
random_dictionary = {'1 - 19 Minutes': 10, 
                     '20 - 39 Minutes': 30, 
                     '40 - 59 Minutes':50,
                     '60 - 79 Minutes':70, 
                     '80 - 100 Minutes':90 } 
# read data
random_ema = pd.read_csv(os.path.join(os.path.realpath(dir), 'random-ema-final.csv'))
eod_ema = pd.read_csv(os.path.join(os.path.realpath(dir), 'eod-ema-final.csv'))
# Make a list of all random EMA event-times between 8AM and 8PM
# Throw away observations for 'when_smoke' is nan or 
# 'More than 30 minutes' to ensure we can calculate a meaningful 
# quantity.
days_smoked = {}
for index, row in random_ema.iterrows():
    try:
        time = datetime.datetime.strptime(row['date'], '%m/%d/%y %H:%M')
    except:
        time = datetime.datetime.strptime(row['date'], '%Y-%m-%d %H:%M:%S')
    if row['when_smoke'] in random_accptresponse:
        time = time - datetime.timedelta(minutes=random_dictionary[row['when_smoke']])
    date = (time.year, time.month, time.day, time.hour)
    if row['participant_id'] not in days_smoked:
        days_smoked[row['participant_id']] = set()
    if 8 <= date[3] < 20 and row['when_smoke'] in random_accptresponse:        
        days_smoked[row['participant_id']].add(date)
# Construct a list of id + dates for eod
# Use to look 
eod_dates = []
for irow in range(0,eod_ema.shape[0]):
    row = eod_ema.iloc[irow]
    if row['status'] == "MISSED":
        continue
    try:
        time = datetime.datetime.strptime(row['date'], '%m/%d/%Y %H:%M')
    except:
        time = datetime.datetime.strptime(row['date'], '%Y-%m-%d %H:%M:%S')
    if time.hour  == 0 or time.hour == 1:
        date = np.array([row['participant_id'], time.year, time.month, time.day-1])
        date = np.append(date, np.array(row[keys]))
    else:
        date = np.array([row['participant_id'], time.year, time.month, time.day])
        date = np.append(date, np.array(row[keys]))
    eod_dates.append(date)
    
eod_dates = np.asarray(eod_dates)
# For participants with both Random EMA and EOD measurements,
# on days when you give both, we ask whether they agree,
# up to the current hour, or +- 1 hour in either direction.
# The +-1 is max/min by 8AM and 8PM respectively.
matching_counts = []
max_iloc = 15; min_iloc = 4
for id in set(days_smoked.keys()) & set(eod_dates[:,0]):
    eod_dates_id = np.where(eod_dates[:,0] == id) 
    eod_dates_subset = eod_dates[eod_dates_id[0],:]
    total_count_id = 0
    hour_count_id_true = 0
    twohour_count_id_true = 0
    if days_smoked[id] == set():
        continue
    for ec_time in days_smoked[id]:
        row_iloc = np.where((eod_dates_subset[:,1:4] == ec_time[0:3]).all(axis=1))[0]
        if not row_iloc.size > 0:
            continue
        total_count_id+=1
        row = eod_dates_subset[row_iloc][0]
        ec_iloc = range(8,20).index(ec_time[3])+4
        if row[ec_iloc]==1:
            hour_count_id_true+=1
        if any(row[range(max(min_iloc, ec_iloc-1), min(max_iloc, ec_iloc+1)+1)] == 1):
            twohour_count_id_true+=1
    if total_count_id > 0:
        matching_counts.append(np.array([total_count_id, hour_count_id_true, twohour_count_id_true], dtype='f'))

matching_counts = np.asarray(matching_counts)

fraction_per_id_onehour = np.divide(matching_counts[:,1],matching_counts[:,0])
fraction_per_id_twohour = np.divide(matching_counts[:,2],matching_counts[:,0])

aggregate_matching_counts = np.sum(matching_counts, axis=0)

aggregate_frac_onehour = aggregate_matching_counts[1]/aggregate_matching_counts[0]
aggregate_frac_twohour = aggregate_matching_counts[2]/aggregate_matching_counts[0]

print('Current hour only:')
print('Aggregated data, Fraction agreement between EC and EOD: %s' % (np.round(aggregate_frac_onehour,3)))
print('Mean of Fraction agreement across indidivuals: %s' % (np.round(np.mean(fraction_per_id_onehour),3)))
print('Standard deviation of Fraction agreement across indidivuals: %s' %  (np.round(np.std(fraction_per_id_onehour),3)))
print()
print('Plus-minus one hour:')
print('Aggregated data, Fraction agreement between EC and EOD: %s' % (np.round(aggregate_frac_twohour,3)))
print('Mean of Fraction agreement across indidivuals: %s' % (np.round(np.mean(fraction_per_id_twohour),3)))
print('Standard deviation of Fraction agreement across indidivuals: %s' %  (np.round(np.std(fraction_per_id_twohour),3)))
# Compute an anova decomposition using the bernoulli likelihood
# This will test if there are significant differences across
# individuals.

llik_onehour = 0; llik_twohour = 0
for i in range(0, fraction_per_id_onehour.size):
    num_ones_onehour = matching_counts[i,1]
    num_zeros_onehour = matching_counts[i,0] - matching_counts[i,1]
    if num_ones_onehour > 0.0:
        llik_onehour += np.multiply(num_ones_onehour, np.log(fraction_per_id_onehour[i]))
    if num_zeros_onehour > 0.0:
        llik_onehour += np.multiply(num_zeros_onehour, np.log(1-fraction_per_id_onehour[i]))
    num_ones_twohour = matching_counts[i,2]
    num_zeros_twohour = matching_counts[i,0] - matching_counts[i,2]
    if num_ones_twohour > 0.0:
        llik_twohour += np.multiply(num_ones_twohour, np.log(fraction_per_id_twohour[i]))
    if num_zeros_twohour > 0.0:
        llik_twohour += np.multiply(num_zeros_twohour, np.log(1-fraction_per_id_twohour[i]))

agg_num_ones = aggregate_matching_counts[1]
agg_num_zeros = aggregate_matching_counts[0] - aggregate_matching_counts[1]
agg_llik_onehour = agg_num_ones*np.log(aggregate_frac_onehour)+agg_num_zeros*np.log(1-aggregate_frac_onehour)

D_onehour = -2*agg_llik_onehour + 2*llik_onehour

agg_num_ones_twohour = aggregate_matching_counts[2]
agg_num_zeros_twohour = aggregate_matching_counts[0] - aggregate_matching_counts[2]
agg_llik_twohour = agg_num_ones_twohour*np.log(aggregate_frac_twohour)+agg_num_zeros_twohour*np.log(1-aggregate_frac_twohour)

D_twohour = -2*agg_llik_twohour + 2*llik_twohour

from scipy.stats import chi2
n = aggregate_matching_counts[0]
k = matching_counts.shape[0]
df = k-1

print('ANOVA p-value for current hour: %s' % (1-chi2.cdf(D_onehour, df)))
print('ANOVA p-value for plus-minus one hour: %s' % (1-chi2.cdf(D_twohour, df)))<jupyter_output>ANOVA p-value for current hour: 1.306049392657549e-05
ANOVA p-value for plus-minus one hour: 0.00012842743519747568
{'fraction_non_alphanumeric': 0.1279843444227006, 'fraction_numerical': 0.037573385518591, 'mean_word_length': 4.141515761234071, 'pattern_counts': {'":': 1, '<': 6, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '34940', 'license': 'CECILL-B', 'n_tokens_mistral': 2973, 'n_tokens_neox': 2746, 'n_words': 638, 'path': 'exploratory_data_analysis/random_eod.ipynb', 'repo_name': 'wdempsey/sense2stop-lvm', 'size': '7665'}