<jupyter_start><jupyter_text>XGBoost for ILI-onset prediction License© Evidation Health, Inc. (2020), All rights reserved. Provided pursuant to a non-exclusive, non-transferable license to Vector Institute (“Licensee”) for the term of, and solely for use in performing the Institute’s obligations under, the Data Use and Non-Disclosure Agreement between Vector Institute and Evidation Health, Inc. dated April 8, 2020 (“Data Use Agreement”).  The Evidation Code shall be considered Evidation Confidential Information under the Data Use Agreement.Evidation reserves all right, title, ownership and interest in and to the Evidation Code existing prior to and after the Effective Date, or created or generated by Evidation at any time, subject to the license granted to Licensee, including Derivative Works, whether created by Evidation or Licensee, of Evidation Code, unless otherwise specified in the Data Use Agreement. At Evidation’s request, which may be made at any time, Licensee shall promptly deliver to Evidation a copy of all Evidation Code in source code form, as it exists at the time of the request, along with the source code and related documentation for any derivative works. Licensee shall reproduce all copyright notices and other proprietary markings or legends contained within or on the Evidation Code on any copies made.Licensee shall not knowingly infringe upon the intellectual property rights of any third party when making Derivative Works to the Evidation Code.Licensee shall not use open source code for development of or in any authorized derivative work of Evidation Code in any manner that would subject the Evidation Code to open source distribution.  To extent there is any open source code in the Evidation Code and there exists a conflict between the this license and any applicable license to open source technology, the provisions of the open source license shall be followed, but only to the minimum extent reasonably necessary to comply with the applicable open source license.<jupyter_code>%load_ext autoreload
%autoreload 2
import pandas as pd; print('Pandas version:', pd.__version__)
import numpy as np; print('Numpy version:', np.__version__)
from multiprocessing import Pool
from time import time
import os
USE_CORES = os.cpu_count() - 1
print('Using cores=', USE_CORES)

pd.set_option('max.rows', 100)
pd.set_option('max.columns', 300)
pd.set_option('mode.chained_assignment', 'raise')
pd.set_option('display.float_format', lambda x: '%.4f' % x) #supress scientific notation

#pandarallel.initialize()<jupyter_output><empty_output><jupyter_text>Set DATA and OUTPUT paths here<jupyter_code># replace argparse args in notebook

@dataclass
class Args():
    hdf_file = '/SET/PATH/HERE/file.hdf' # HDF file with keys `activity` and `survey`
    test_csv = None # Can provide `test_participants.csv`
    output_dir = '/SET/PATH/HERE/OUTPUT_DIR' # Provide output directory path
    use_cores = None # If None, will use N-1 cores 

args = Args()
idx = pd.IndexSlice<jupyter_output><empty_output><jupyter_text>Lancet-RHR method* Using `measurement` & `mask` dataframes and columns:    * `heart_rate__resting_heart_rate` for RHR    * `sleep__asleep__sum` for Sleep    * `active_fitbit__sum` for Wear time* Data density requirements:    * 7-day window for weekly RHR calculation    * Minimum **4** valid RHR values required in window    * Minimum **100** valid-RHR-value days of past data for baseline calculation    * Minimum **1000** minutes of wear-time per valid day* Score & Threshold Rule:    * z-score = (Weekly Avg. RHR/Sleep  - Baseline Avg. RHR/Sleep)/(Baseline Std. Dev. RHR/Sleep)     * **z-RHR > 0.5 AND z-Sleep > -0.5**<jupyter_code>if not (os.path.exists(args.output_dir)):
        os.mkdir(args.output_dir)

OUTPUT_PATH = os.path.join(args.output_dir, 'lancet_rhr-testset_results.csv')
        
### Number of cores for multiprocessing
USE_CORES = args.use_cores
if USE_CORES is None:
    USE_CORES = os.cpu_count() - 1

### Default keys
activity_key = 'activity'
labels_key = 'survey'

### Define parameters
feature_name_rhr='heart_rate__resting_heart_rate'
feature_name_sleep='sleep__asleep__sum'
feature_name_wear='active_fitbit__sum'
measurement_col='measurement'
mask_col='mask'
window='7d'
min_week_valid_days=4
min_baseline_valid_days=100
min_wear_time_mins=1000
rhr_thresh=0.5
sleep_thresh=-0.5
label_col = 'ili'
output_label_col = 'label'<jupyter_output><empty_output><jupyter_text>Define functions and utils<jupyter_code>def lancet_rhr_method(x, feature_name_rhr='heart_rate__resting_heart_rate', feature_name_sleep='sleep__asleep__sum',
                      feature_name_wear='active_fitbit__sum', measurement_col='measurement', mask_col='mask',
                      window='7d', min_week_valid_days=4, min_baseline_valid_days=100, min_wear_time_mins=1000,
                      rhr_thresh=0.5, sleep_thresh=-0.5):
    """ 
    Input: single-participant time-series data-frame with RHR, Sleep, Wear-time columns
    Output: Lancet-RHR decision and scores per date in data-frame
    """
    
    x = x.sort_index(level=1).droplevel(0) #drop participant_id since rolling/expanding not supported for multi-index
    
    drop_rhr = ((~x[mask_col, feature_name_rhr]) | (x[measurement_col, feature_name_wear] < min_wear_time_mins)).rename('drop_rhr')
    drop_sleep = ((~x[mask_col, feature_name_sleep]) | (x[measurement_col, feature_name_wear] < min_wear_time_mins)).rename('drop_sleep')

    valid_rhr = x[measurement_col, feature_name_rhr].mask(drop_rhr).rename('valid_rhr')
    valid_sleep = x[measurement_col, feature_name_sleep].mask(drop_sleep).rename('valid_sleep')

    baseline_valid_days = (0 + ~drop_rhr).expanding().sum().shift(freq=window).rename('baseline_valid_days')
    baseline_rhr_mean = valid_rhr.expanding().mean().shift(freq=window).rename('baseline_rhr_mean') # .mask(baseline_valid_days < min_baseline_valid_days).rename('baseline_rhr_mean') # only those with 100 days
    baseline_rhr_stdv = valid_rhr.expanding().std().shift(freq=window).rename('baseline_rhr_stdv') # .mask(baseline_valid_days < min_baseline_valid_days).rename('baseline_rhr_stdv')

    baseline_sleep_mean = valid_sleep.expanding().mean().shift(freq=window).rename('baseline_sleep_mean') #.mask(baseline_valid_days < min_baseline_valid_days).rename('baseline_sleep_mean') # only those with 100 days
    baseline_sleep_stdv = valid_sleep.expanding().std().shift(freq=window).rename('baseline_sleep_stdv') #.mask(baseline_valid_days < min_baseline_valid_days).rename('baseline_sleep_stdv')


    weekly_valid_days = (0 + ~drop_rhr).rolling(window).sum().rename('weekly_valid_days')
    weekly_rhr_mean = valid_rhr.rolling(window).mean().rename('weekly_rhr_mean') #.mask(weekly_valid_days < min_week_valid_days).rename('weekly_rhr_mean')
    weekly_sleep_mean = valid_sleep.rolling(window).mean().rename('weekly_sleep_mean') #.mask(weekly_valid_days < min_week_valid_days).rename('weekly_sleep_mean')

    daily_z_rhr = ((weekly_rhr_mean - baseline_rhr_mean)/baseline_rhr_stdv).rename('daily_z_rhr')
    daily_z_sleep = ((weekly_sleep_mean - baseline_sleep_mean)/baseline_sleep_stdv).rename('daily_z_sleep')

    lancet_rhr_decision = (0 + ((daily_z_rhr > rhr_thresh) & (daily_z_sleep > sleep_thresh))).rename('lancet_rhr_decision')
    lancet_rhr_decision[daily_z_rhr.isna() | daily_z_sleep.isna() | baseline_valid_days < min_baseline_valid_days] = np.nan
    
    return pd.concat([lancet_rhr_decision, daily_z_rhr, daily_z_sleep, weekly_rhr_mean, weekly_sleep_mean, 
                      weekly_valid_days, baseline_valid_days, baseline_rhr_mean, baseline_rhr_stdv, 
                      baseline_sleep_mean, baseline_sleep_stdv, valid_rhr, drop_rhr, valid_sleep, drop_sleep],
                      join='outer', axis=1).loc[x.index,:] #return values for dates in original dataframe

def lancet_rhr_parallel_apply(x):
    y = lancet_rhr_method(x[1])
    y.index = pd.MultiIndex.from_product([[x[0]], y.index])
    y.index.names = ['participant_id', 'date']
    return y

def prop_table(x, dropna=True):
    """ Unique values and count/percentage in pandas series"""
    tmp = (x.value_counts(sort=False, dropna=dropna).reset_index()
           .merge((100 * x.value_counts(sort=False, normalize=True, dropna=dropna)).round(2).reset_index(), on='index',
                  how='inner'))
    tmp.columns = [x.name, 'count', 'percent']
    tmp = tmp.sort_values('count', ascending=False)
    tot = x.notnull().sum() if dropna else len(x)
    return tmp.append(pd.DataFrame([['Total', tot, 100]], columns=tmp.columns), ignore_index=True)<jupyter_output><empty_output><jupyter_text>Load activity `measurement` and `mask` dataframe columns<jupyter_code>#Load processed-imputed activity
df = pd.read_hdf(args.hdf_file, activity_key)

print(df.shape)
print(df.columns.get_level_values(0).unique())<jupyter_output><empty_output><jupyter_text>Load survey dataframe (for ILI labels)<jupyter_code>#Load survey
mf = pd.read_hdf(args.hdf_file, labels_key)

print(mf.shape)
mf.columns.get_level_values(0).unique()<jupyter_output><empty_output><jupyter_text>Join activity and labels<jupyter_code>mf.columns = pd.MultiIndex.from_product([['labels'], mf.columns])
df = df.join(mf, how='left')
print(df.shape)
df.columns.get_level_values(0).unique()
pd.Series(df.index.get_level_values(1)).describe()<jupyter_output><empty_output><jupyter_text>Test with single participant<jupyter_code>x = df.loc[df.index.get_level_values(0)=='viAz0JRFPhIUyKb1',:].copy()#.droplevel(0)
print(x.shape)
x.loc[:,idx['measurement', ['heart_rate__resting_heart_rate', 'sleep__asleep__sum', 'active_fitbit__sum']]].head(10)
x.loc[:,idx['measurement', ['heart_rate__resting_heart_rate', 'sleep__asleep__sum', 'active_fitbit__sum']]].tail(10)
y = lancet_rhr_method(x)
print(y.shape)
y[['lancet_rhr_decision', 'daily_z_rhr', 'daily_z_sleep']].head(10)
y[['lancet_rhr_decision', 'daily_z_rhr', 'daily_z_sleep']].tail(10)<jupyter_output><empty_output><jupyter_text>Test with multiple participants and `parallel_apply`<jupyter_code>tmp = df.loc[idx[list(df.sample(5).index.get_level_values(0)),:],:]
print(tmp.shape)
tmp.sample(5)
print('Using cores=', USE_CORES)

p = Pool(USE_CORES)
stime = time()
foo = pd.concat(p.map(lancet_rhr_parallel_apply, tmp.groupby('participant_id')))
etime = time()
print(f'Time: {etime-stime:.2f}s')
print(foo.shape)
foo.head()
foo.sample(5)<jupyter_output><empty_output><jupyter_text>(MAIN) Run on whole dataset<jupyter_code>stime = time()
out = pd.concat(p.map(lancet_rhr_parallel_apply, df.groupby('participant_id')))
etime = time()
print(f'Time: {etime-stime:.2f}s')
print(out.shape)<jupyter_output><empty_output><jupyter_text>Test NaN/0/1 prediction distribution<jupyter_code>prop_table(out.lancet_rhr_decision, dropna=False)<jupyter_output><empty_output><jupyter_text>Merge with day-level ILI labels<jupyter_code>tmp = mf[[label_col]].rename(columns={label_col: output_label_col})
print(tmp.shape)
display(tmp.head())<jupyter_output><empty_output><jupyter_text>TANH the z-scores before saving result to deal with infinities<jupyter_code>out = tmp.join(out, how='left')
out['tanh_daily_z_rhr'] = out.loc[:, 'daily_z_rhr'].apply(np.tanh)
out['label'] = out['label'].astype(int)
out['score'] = out['lancet_rhr_decision']

print(out.shape)
display(out.head())
print(f'RHR NaN frac={out.daily_z_rhr.isna().sum()/out.shape[0]:.2f}, Sleep NaN frac={out.daily_z_sleep.isna().sum()/out.shape[0]:.2f}')
print()
pd.concat([out.tanh_daily_z_rhr.describe(), out.daily_z_rhr.describe(), out.daily_z_sleep.describe()], axis=1)
col_order = ['label', 'score', 'tanh_daily_z_rhr', 'lancet_rhr_decision', 'daily_z_rhr', 'daily_z_sleep',
             'weekly_rhr_mean', 'weekly_sleep_mean', 'weekly_valid_days',
             'baseline_valid_days', 'baseline_rhr_mean', 'baseline_rhr_stdv',
             'baseline_sleep_mean', 'baseline_sleep_stdv', 'valid_rhr', 'drop_rhr', 'valid_sleep', 'drop_sleep']

out = out[col_order]
print(out.shape)
display(out.sample(10))<jupyter_output><empty_output><jupyter_text>Save output to `pickle` and results for test participants to `csv`<jupyter_code>out.to_pickle(OUTPUT_PATH.replace('testset_results.csv', 'full_output.pkl'))
### Prepare dataframe for saving results
dump = out.reset_index()[['participant_id', 'date', 'label', 'score', 'tanh_daily_z_rhr']]

## SAVE ALL RESULT CSV
dump.to_csv(OUTPUT_PATH.replace('testset', 'all'), index=None)

## Load test participants list & Save result for test set only
if args.test_csv is not None:
    test_idx = pd.read_csv(args.test_csv)

    dump[dump.participant_id.isin(test_idx.iloc[:,0])].to_csv(OUTPUT_PATH, index=None)<jupyter_output><empty_output><jupyter_text>Metrics and performance Lancet-RHR method makes a prediction 42% of the time with recall of 30% (precision 7%) for ILI days* Does not use imputed values* Min. 100 "valid" days for baseline stats (wear-time > 1000 mins and RHR available for day to be "valid")* Min. 4 "valid" days in a week for weekly stats* Rule: z-RHR > 0.5 SD and z-Sleep > -0.5 SD<jupyter_code>pd.crosstab(out.label, out.score, margins=True)
pd.crosstab(out.label, out.score, normalize='index')
from sklearn.metrics import roc_auc_score, roc_curve, precision_recall_curve, classification_report, average_precision_score, plot_precision_recall_curve
keep = out.lancet_rhr_decision.notna()
y_true = out.loc[keep, 'label']
y_pred = out.loc[keep, 'lancet_rhr_decision']
out_rhr_auroc = roc_auc_score(y_true, y_pred)
out_rhr_avgprec = average_precision_score(y_true, y_pred)
print(f'AUROC = {out_rhr_auroc:.3f}')
print(f'Avg. Precision = {out_rhr_avgprec:.3f}')
print(classification_report(y_true, y_pred))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12152676221117285, 'fraction_numerical': 0.009578824217607487, 'mean_word_length': 6.24417372881356, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 68, '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': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '521688', 'license': 'Apache-2.0', 'n_tokens_mistral': 4878, 'n_tokens_neox': 4673, 'n_words': 1009, 'path': 'nb_lancet_rhr_method.ipynb', 'repo_name': 'goldenberg-lab/DearWatchCOVID-19', 'size': '13676'}