<jupyter_start><jupyter_text>Data PreparationThis notebook is used to clean up, preprocess, and aggregate data necessary for STXGB models.Please refer to the `Method` section of the article and Supplementary Information document for more information.In order to run this notebook, you have to download the required datasets from this [address](https://drive.google.com/drive/u/1/folders/1laAZFCvsPLLaKDvg0isTMMr20kMe0x_r). Once downloaded, set the directory in which you have save the files as `input_directory` in the cell below.<jupyter_code>input_directory = './input/'
import pandas as pd
import numpy as np
import geopandas as gpd
import matplotlib
import matplotlib.pyplot as plt
import platform
from datetime import datetime, timedelta
#Set date string formatting based on operating system

if platform.system() == 'Windows':
    conversion_format = '%#m/%#d/%y'
else:
    conversion_format = '%-m/%-d/%y'
data = gpd.read_file( input_directory + 'Contiguous_US.geojson')

# Or alternatively:
# url='https://drive.google.com/file/d/1MVyLzzHl3hzno4o1rLZtI0peqIi23zsr/view?usp=sharing'
# url_counties='https://drive.google.com/uc?id=' + url.split('/')[-2]
# data = gpd.read_file(url_counties)
data['STATEFP'] = data.apply(lambda L: L.GEOID[:2], axis=1)
global number_counties 
number_counties = data.shape[0] #3103
data.sort_values(by='GEOID', inplace=True)<jupyter_output><empty_output><jupyter_text>Load COVID data and apply smoothing To alleviate inconsistencies in reporting COVID-19 cases, we apply a 7-day moving average to the case data published by JHU.<jupyter_code>def get_JH_covid_data(target, smooth):
    
    '''
    Parameters:
    --------------
        target: str
            the target variable: either 'case' or 'death'
            
        smooth: bool
            wether to smooth the data frame or not.
            The smoothing is done by using a 7-day rolling average   
    '''
    
    assert isinstance(smooth, bool), "Smooth must be a boolean variable!"
    
    base_url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/'
    
    
    if target == 'case':
        jh_data_url = base_url + 'csse_covid_19_time_series/time_series_covid19_confirmed_US.csv'
        cols_to_drp = ['UID', 'iso2', 'iso3', 'code3','Country_Region', 'Lat', 'Long_']

    elif target=='death':
        jh_data_url = base_url + 'csse_covid_19_time_series/time_series_covid19_deaths_US.csv'
        cols_to_drp = ['UID', 'iso2', 'iso3', 'code3','Country_Region', 'Lat', 'Long_','Population']
    else:
        print("invalid argument for target. Acceptable values are: 'case' or 'death'")
        return None

    jh_covid_df = pd.read_csv(jh_data_url)

    # preprocessing JH COVID data
    jh_covid_df.dropna(axis=0, how='any', inplace=True)

    jh_covid_df['FIPS'] = jh_covid_df['FIPS'].astype('int64')

    jh_covid_df.drop(columns=cols_to_drp, inplace=True)

    #Important: check to see the column index is adherent to the imported df

    first_date = datetime.strptime(jh_covid_df.columns[4], '%m/%d/%y').date()

    last_date = datetime.strptime(jh_covid_df.columns[-1], '%m/%d/%y').date()


    current_date = last_date

    previous_date = last_date - timedelta (days=1)


    while current_date > first_date:

        #For unix, replace # with - in the time format

        current_col = current_date.strftime(conversion_format) #replace # with - in Mac or Linux

        previous_col = previous_date.strftime(conversion_format)

        jh_covid_df[previous_col] = np.where(jh_covid_df[previous_col] > jh_covid_df[current_col], 
                                             jh_covid_df[current_col], jh_covid_df[previous_col])

        current_date = current_date - timedelta(days=1)

        previous_date = previous_date - timedelta(days=1)
        
    
    if smooth:
        jh_covid_df.iloc[:,4:] = jh_covid_df.iloc[:,4:].rolling(7,min_periods=1,axis=1).mean()


    return jh_covid_df
covid_df = get_JH_covid_data('case', True)<jupyter_output><empty_output><jupyter_text>Load Facebook Movement Range Data<jupyter_code>fb_mobility = pd.read_csv(input_directory + 'movement-range-2021-03-02.txt', sep="\t", dtype={'polygon_id':str})
fb_mobility_us = fb_mobility[fb_mobility['country']=='USA']
# Number of unique counties for which we have at least one day of data
len(fb_mobility_us['polygon_id'].unique())
# counting Counties in the contiguous US for which there is no data in FB mobility
contiguous_fips = set(data['GEOID']) # number of unique fips: 3103
mobility_fips = set(fb_mobility_us['polygon_id']) # number of unique fips: 2694
i = 0
missing_fips = []
for fips in contiguous_fips:
    if (fips in mobility_fips):
        i+=1
    else:
        missing_fips.append(fips)
# number of Counties in the contiguous US for which there is no data in FB mobility
len(missing_fips)
# Create a new dataframe as transpose of the above, with days as columns and counties as rows

relative_df = pd.DataFrame(columns=fb_mobility_us['ds'].unique(), index=data['GEOID'].unique(), dtype=float)
ratio_df = pd.DataFrame(columns=fb_mobility_us['ds'].unique(), index=data['GEOID'].unique(), dtype=float)
fb_mobility_us
idx_contiguous = fb_mobility_us.index[fb_mobility_us['polygon_id'].isin(contiguous_fips)]
fb_mobility_contiguous = fb_mobility_us.loc[idx_contiguous]
import time
start = time.time()

for index, row in fb_mobility_contiguous.iterrows():
    relative_df.loc[row['polygon_id']][row['ds']] = row['all_day_bing_tiles_visited_relative_change']
    ratio_df.loc[row['polygon_id']][row['ds']] = row['all_day_ratio_single_tile_users']
    

end = time.time()
print(end - start)
relative_df.head()
print(relative_df.shape , ratio_df.shape)<jupyter_output>(3103, 367) (3103, 367)<jupyter_text>Impute FB mobility dataframesThe two dataframes above have a lot of Nan values which should be imputed by state average<jupyter_code>temp_ratio_df = data[['GEOID', 'STATEFP']].merge(ratio_df, left_on='GEOID', right_index=True, how='left')
for col in temp_ratio_df.iloc[:,2:].columns:
    temp_ratio_df[col] = temp_ratio_df.groupby('STATEFP')[col].transform(lambda x: x.fillna(x.mean()))
temp_ratio_df.isna().sum().sum()
temp_relative_df = data[['GEOID', 'STATEFP']].merge(relative_df, left_on='GEOID', right_index=True, how='left')
for col in temp_relative_df.iloc[:,2:].columns:
    temp_relative_df[col] = temp_relative_df.groupby('STATEFP')[col].transform(lambda x: x.fillna(x.mean()))
temp_relative_df.isna().sum().sum()
relative_df_smooth = temp_relative_df.copy()
ratio_df_smooth = temp_ratio_df.copy()

relative_df_smooth.iloc[:,2:] = relative_df_smooth.iloc[:,2:].rolling(7,min_periods=1, axis=1).mean()
ratio_df_smooth.iloc[:,2:] = ratio_df_smooth.iloc[:,2:].rolling(7,min_periods=1, axis=1).mean()
relative_df_smooth.iloc[:,2:].columns
relative_df_smooth.shape<jupyter_output><empty_output><jupyter_text>Add Social Proximity to Cases<jupyter_code>SCI_df = pd.read_csv(input_directory + 'SCI_matrix.csv', dtype={'Unnamed: 0':str})
SCI_df.set_index('Unnamed: 0', inplace=True)
# Create normalized SCI. It is calculated by dividing all the columns of the sci_matrix by the sum of the rpw
# This would give us the second term in social proximity formula above

sci_matrix_normal = SCI_df.div(SCI_df.sum(axis=1), axis=0)
# set diagonal to zero
sci_matrix_normal.values[[np.arange(sci_matrix_normal.shape[0])]*2] = 0<jupyter_output>C:\Users\BehzadVahediTorghabe\anaconda3\lib\site-packages\ipykernel_launcher.py:2: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.<jupyter_text>**The matrix above is created for the entire US, but we are using contiguous US data here, therefore some rows and columns should be removed**<jupyter_code>to_drop=[]

for index in sci_matrix_normal.index:
    if not index in contiguous_fips:
        to_drop.append(index)
sci_matrix_normal.drop(to_drop, inplace=True)
sci_matrix_normal.drop(to_drop, axis=1, inplace=True)
sci_matrix_normal.shape<jupyter_output><empty_output><jupyter_text>Add SafeGraph mobility featuresUpdated based on forecast hub dates<jupyter_code>safegraph_mobility = pd.read_csv(input_directory + 'safegraph_mobility.csv', dtype={'county_fips':str})
safegraph_contiguous = safegraph_mobility[safegraph_mobility['county_fips'].isin(contiguous_fips)]
len(safegraph_contiguous['county_fips'].unique())
temp_cols = safegraph_contiguous.drop(['start_date', 'end_date', 'base_start', 'base_end'], axis=1)
safegraph_metrics = temp_cols.columns
len(safegraph_metrics)<jupyter_output><empty_output><jupyter_text>Load Temperature Data<jupyter_code>max_temp = pd.read_csv(input_directory + 'max_temp_df_2021.csv', dtype={'GEOID':str})
min_temp = pd.read_csv(input_directory + 'min_temp_df_2021.csv', dtype={'GEOID':str})<jupyter_output><empty_output><jupyter_text>Helper functions<jupyter_code># return FCI-normal table for the input date
# set path_to_fci to where FCI matrices are stored
def get_normal_fci(date):
    path_to_fci = './output/' + str(date.year) + '/'+ date.strftime('%m') + 
                '/FCI_normal/' + date.strftime('%Y-%m-%d') + '-FCI-normal.csv'
    fci_norm = pd.read_csv(path_to_fci, dtype={'Unnamed: 0':str})
    fci_norm.set_index('Unnamed: 0', inplace=True)
    
    to_drop=[]

    for index in fci_norm.index:
        if not index in contiguous_fips:
            to_drop.append(index)
            
    fci_norm.drop(to_drop, inplace=True)
    fci_norm.drop(to_drop, axis=1, inplace=True)
    return fci_norm
# calculates weekly average FPC using the end date and the start date of the week
# the input to this fuction should be of type datetime
def weekly_mean_FPC(end_date, start_date, logged=False):
    
    dates = [end_date]
    while end_date> start_date:
        end_date -= timedelta(days=1)
        dates.append(end_date)
    
    
    temp = data[['GEOID','FIPS', 'POPULATION']]
    
    
    for date in dates:
        # convert date to String
        date_str = date.strftime('%Y-%m-%d')
        
        temp = temp.merge(covid_df_contiguous[['FIPS', date_str]], on='FIPS', how='left')
        
        if logged:
            temp['inc_rate_' + date_str] = np.log(temp[date_str] / temp['POPULATION'] * 10000 + 1)
            
        else:
            temp['inc_rate_' + date_str] = temp[date_str] / temp['POPULATION'] * 10000
        
        normal_fci = get_normal_fci(date)
        normal_fci = normal_fci.merge(temp, left_index= True, right_on='GEOID')
        
        normal_fci['fpc_'+ date_str] = np.dot(normal_fci.iloc[:,:number_counties], normal_fci['inc_rate_' + date_str])
        
    
    normal_fci['mean_fpc'] = normal_fci.iloc[:,-len(dates):].mean(axis=1)
        
    return normal_fci[['GEOID','mean_fpc']]    
# calculates weekly average SPC
# the input to this fuction should be of type datetime
def weekly_mean_SPC(end_date, start_date, logged=False):
    
    dates = [end_date]
    while end_date> start_date:
        end_date -= timedelta(days=1)
        dates.append(end_date)
    
    
    temp = data[['GEOID','FIPS', 'POPULATION']]
    
    for date in dates:
        # convert date to String
        date_str = date.strftime('%Y-%m-%d')
        
        temp = temp.merge(covid_df_contiguous[['FIPS', date_str]], on='FIPS', how='left')
        
        if logged:
            temp['inc_rate_' + date_str] = np.log(temp[date_str] / temp['POPULATION'] * 10000 + 1)
            
        else:
            temp['inc_rate_' + date_str] = temp[date_str] / temp['POPULATION'] * 10000
        
        
        normal_sci = sci_matrix_normal.merge(temp, left_index= True, right_on='GEOID')
        
        normal_sci['spc_'+ date_str] = np.dot(normal_sci.iloc[:,:number_counties], normal_sci['inc_rate_' + date_str])
        
    
    normal_sci['mean_spc'] = normal_sci.iloc[:,-len(dates):].mean(axis=1)
        
    return normal_sci[['GEOID','mean_spc']]
# the input to this fuction should be of type datetime.
# returns a subset of FB movement range dfs based on the given week
def weekly_fb_mobility(end_date, start_date, df):
    
    dates = [end_date]
    while end_date> start_date:
        end_date -= timedelta(days=1)
        dates.append(end_date)
    
    dates_str=[]
    for date in dates:
        # convert date to String
        dates_str.append(date.strftime('%Y-%m-%d'))
    
    return df[['GEOID', *dates_str]]
# calculate slope features

from scipy.stats import linregress

def linear_reg(week_df):
    
    x = np.arange(1,(week_df.shape[1]),1)
    x = (x - np.mean(x))/ np.std(x)
    
    slopes=[]
    
    for index, row in week_df.iloc[:,1:].iterrows():
        y = row
        y = (y - np.mean(y))/ np.std(y)
        slopes.append(linregress(x, y)[0])
        
    week_df.loc[:,'slope'] = slopes
    return week_df[['GEOID','slope']]<jupyter_output><empty_output><jupyter_text>Final function to combine all features generated aboveThis function generates a dataframe and for a given date, will add the following features to the dataframe- incidence rate data- FB mobility data (ratio, relative)- SPC (facebook SCI and incidence rates)- SafeGraph mobility - FPC (FCI and incidence rate)For each period, there is a 5 week difference between the actual date (t) and the start of the second lag. For example if `T: Oct 1 (Sep 24 to Oct 1)`, then `T-1: Sep 10 to Sep 24`, and `T-2: August 27 to Sep 10`.Since the earliest day for which we have FB mobility data is March 1, the rearliest  (end) date for T will be April 5th.<jupyter_code>covid_df_contiguous = data[['FIPS','STATEFP','COUNTYFP','GEOID']].merge(covid_df, on='FIPS', how='left')
# Main function
# the input to the funtion is the end date for which we want to add data

from datetime import datetime, timedelta


def add_lagged_data(date):
    global data
    
    T_end = date
    T_start = T_end - timedelta(days=6)
    
    T_1_end = T_start - timedelta(days=1)
    T_1_start = T_1_end - timedelta(days=6)
    
    T_2_end = T_1_start - timedelta(days=1)
    T_2_start = T_2_end - timedelta(days=6)
    
    T_3_end = T_2_start - timedelta(days=1)
    T_3_start = T_3_end - timedelta(days=6)
    
    T_4_end = T_3_start - timedelta(days=1)
    T_4_start = T_4_end - timedelta(days=6)
    
    # These dates are used for cumulative cases (Saturday to Saturday)
    T_start_case = T_end - timedelta(days=7)
    T_1_start_case = T_1_end - timedelta(days=7)
    T_2_start_case = T_2_end - timedelta(days=7)
    T_3_start_case = T_3_end - timedelta(days=7)
    T_4_start_case = T_4_end - timedelta(days=7)
    

    dates = [T_end.strftime('%Y-%m-%d'), T_start.strftime('%Y-%m-%d'),
             T_1_end.strftime('%Y-%m-%d'), T_1_start.strftime('%Y-%m-%d'), 
             T_2_end.strftime('%Y-%m-%d'), T_2_start.strftime('%Y-%m-%d'),
             T_3_end.strftime('%Y-%m-%d'), T_3_start.strftime('%Y-%m-%d'), 
             T_4_end.strftime('%Y-%m-%d'), T_4_start.strftime('%Y-%m-%d')]
    
    
    dates_non_str = [T_end, T_start,
             T_1_end, T_1_start, 
             T_2_end, T_2_start,
             T_3_end, T_3_start, 
             T_4_end, T_4_start]
    
    dates_case = [T_end.strftime('%Y-%m-%d'), T_start_case.strftime('%Y-%m-%d'),
             T_1_end.strftime('%Y-%m-%d'), T_1_start_case.strftime('%Y-%m-%d'), 
             T_2_end.strftime('%Y-%m-%d'), T_2_start_case.strftime('%Y-%m-%d'),
             T_3_end.strftime('%Y-%m-%d'), T_3_start_case.strftime('%Y-%m-%d'), 
             T_4_end.strftime('%Y-%m-%d'), T_4_start_case.strftime('%Y-%m-%d')]
    
    
    dates_non_str_case = [T_end, T_start_case,
             T_1_end, T_1_start_case, 
             T_2_end, T_2_start_case,
             T_3_end, T_3_start_case, 
             T_4_end, T_4_start_case]

    
    temp = data.copy()
    
    temp['date_end_period'] = T_end.strftime('%Y-%m-%d')
    temp['date_start_period'] = T_start.strftime('%Y-%m-%d')
    temp['date_end_lag'] = T_1_end.strftime('%Y-%m-%d')
    temp['date_start_lag'] = T_4_start.strftime('%Y-%m-%d')
    
    time_periods = ['T_end', 'T_start', 'T_1_end', 'T_1_start','T_2_end','T_2_start',
                    'T_3_end','T_3_start','T_4_end','T_4_start']
    i = 0
    for period in time_periods:
        
        
        temp = temp.merge(covid_df_contiguous[['GEOID',dates_case[i]]], on='GEOID', how='left')
        temp.rename(columns={dates_case[i]:'case_'+ period}, inplace=True)

        
        temp['inc_rate_' + period] = temp['case_'+ period] / temp['POPULATION'] * 10000
        
        temp = temp.merge(relative_df_smooth[['GEOID',dates[i]]], on='GEOID', how='left')
        temp.rename(columns={dates[i]:'relative_'+ period}, inplace=True) 


        temp = temp.merge(ratio_df_smooth[['GEOID',dates[i]]], on='GEOID', how='left')
        temp.rename(columns={dates[i]:'ratio_'+ period}, inplace=True)
        
        # The same date is used as the input to weekly_mean_SPC function to calculate
        # SPC for that given date (instead of an average over a period)
        temp = temp.merge(weekly_mean_SPC(dates_non_str_case[i],dates_non_str_case[i], logged=False), on='GEOID', how='left')
        temp.rename(columns={'mean_spc':'SPC_'+ period}, inplace=True)

        # simiar to SPC, add FPC values
        temp = temp.merge(weekly_mean_FPC(dates_non_str_case[i],dates_non_str_case[i], logged=False), on='GEOID', how='left')
        temp.rename(columns={'mean_fpc':'FPC_'+ period}, inplace=True)
        
        
        # add logged SPC (defined as log(delta incidence rate)*sci/sum(sci))
        temp = temp.merge(weekly_mean_SPC(dates_non_str_case[i],dates_non_str_case[i], logged=True), on='GEOID', how='left')
        temp.rename(columns={'mean_spc':'SPC_logged_'+ period}, inplace=True)
        
        
        # add logged FPC (defined as log(delta incidence rate)*fci/sum(fci))
        temp = temp.merge(weekly_mean_FPC(dates_non_str_case[i],dates_non_str_case[i], logged=True), on='GEOID', how='left')
        temp.rename(columns={'mean_fpc':'FPC_logged_'+ period}, inplace=True)

        # add raw John Hopkins case data
        temp = temp.merge(jh_covid_df[['FIPS',dates_non_str_case[i].strftime('%#m/%#d/%y')]], on='FIPS', how='left')
        temp.rename(columns={dates_non_str_case[i].strftime('%#m/%#d/%y'):'case_JH_'+ period}, inplace=True)
        
        # add smoothed John Hopkins case data
        temp = temp.merge(covid_df_contiguous[['FIPS',dates_case[i]]], on='FIPS', how='left')
        temp.rename(columns={dates_case[i]:'case_JH_smoothed_'+ period}, inplace=True)
       
        
        
        i+=1
    
    times = ['T', 'T_1', 'T_2', 'T_3', 'T_4']
    
    j = 0
    for period in times:
        temp['DELTA_INC_RATE_' + period] = temp['inc_rate_' + time_periods[j]] - temp['inc_rate_' + time_periods[j+1]]
        
        temp['DELTA_REL_MOB_' + period] = temp['relative_' + time_periods[j]] - temp['relative_' + time_periods[j+1]]
        temp['DELTA_RATIO_MOB_' + period] = temp['ratio_' + time_periods[j]] - temp['ratio_' + time_periods[j+1]]
        
        temp['DELTA_SPC_' + period] = temp['SPC_' + time_periods[j]] - temp['SPC_' + time_periods[j+1]]
        temp['DELTA_FPC_' + period] = temp['FPC_' + time_periods[j]] - temp['FPC_' + time_periods[j+1]]
        temp['DELTA_SPC_LOGGED_' + period] = temp['SPC_logged_' + time_periods[j]] - temp['SPC_logged_' + time_periods[j+1]]
        temp['DELTA_FPC_LOGGED_' + period] = temp['FPC_logged_' + time_periods[j]] - temp['FPC_logged_' + time_periods[j+1]]
        
        temp['DELTA_CASE_JH_' + period] = temp['case_JH_'+ time_periods[j]] - temp['case_JH_'+ time_periods[j+1]]
        temp['DELTA_CASE_JH_SMOOTH_' + period] = temp['case_JH_smoothed_'+ time_periods[j]] - 
                                                 temp['case_JH_smoothed_'+ time_periods[j+1]]
        
        # mean incidence rate is calculated between Sunday and Saturday
        temp['MEAN_INC_RATE_' + period] = covid_df_contiguous[pd.date_range(start=dates[j+1], end=dates[j])
                                            .strftime('%Y-%m-%d')].mean(axis=1) / temp['POPULATION'] * 10000
        temp['MEAN_REL_MOB_' + period] = relative_df_smooth[pd.date_range(start=dates[j+1], end=dates[j])
                                            .strftime('%Y-%m-%d')].mean(axis=1)
        temp['MEAN_RATIO_MOB_' + period] = ratio_df_smooth[pd.date_range(start=dates[j+1], end=dates[j])
                                            .strftime('%Y-%m-%d')].mean(axis=1)

        
        # add Safegraph mobility features
        safegraph_data = safegraph_contiguous[safegraph_contiguous['end_date']==dates[j]][safegraph_metrics]
        temp = temp.merge(safegraph_data, left_on='GEOID', right_on='county_fips')
        
        rename_dict = dict()
        for col in safegraph_metrics[1:]:
            rename_dict[col] = col + '_' + period
            
        temp.rename(columns=rename_dict, inplace=True)
        
        
        # add MEAN_FPC
        temp = temp.merge(weekly_mean_FPC(dates_non_str_case[j], dates_non_str_case[j+1]), on='GEOID', how='left')
        temp.rename(columns={'mean_fpc':'MEAN_FPC_'+ period}, inplace=True)
        
        # add logged MEAN_FPC 
        temp = temp.merge(weekly_mean_FPC(dates_non_str_case[j], dates_non_str_case[j+1], logged=True), on='GEOID', how='left')
        temp.rename(columns={'mean_fpc':'MEAN_FPC_LOGGED_'+ period}, inplace=True)
        
        # add MEAN_SPC
        temp = temp.merge(weekly_mean_SPC(dates_non_str_case[j], dates_non_str_case[j+1]), on='GEOID', how='left')
        temp.rename(columns={'mean_spc':'MEAN_SPC_'+ period}, inplace=True)
        
        # add logged MEAN_SPC
        temp = temp.merge(weekly_mean_SPC(dates_non_str_case[j], dates_non_str_case[j+1], logged=True), on='GEOID', how='left')
        temp.rename(columns={'mean_spc':'MEAN_SPC_LOGGED_'+ period}, inplace=True)
        
        # add FB mobility slopes
        ratio_slope = linear_reg(weekly_fb_mobility(dates_non_str[j], dates_non_str[j+1], ratio_df_smooth))
        temp = temp.merge(ratio_slope, on='GEOID', how='left')
        temp.rename(columns={'slope':'SLOPE_RATIO_MOB_'+ period}, inplace=True)
        
        relative_slope = linear_reg(weekly_fb_mobility(dates_non_str[j], dates_non_str[j+1], relative_df_smooth))
        temp = temp.merge(relative_slope, on='GEOID', how='left')
        temp.rename(columns={'slope':'SLOPE_REL_MOB_'+ period}, inplace=True)

        
        
        # add temperature features
        # to update for the new dates, min and max temperature are used with one day offset
        adj_temp_date = (dates_non_str[j] + timedelta(days=1)).strftime('%Y-%m-%d')
        
        temp = temp.merge(max_temp[['GEOID', adj_temp_date]], on='GEOID', how='left')
        temp.rename(columns={adj_temp_date:'MAX_TEMP_'+ period}, inplace=True)
        
        temp = temp.merge(min_temp[['GEOID', adj_temp_date]], on='GEOID', how='left')
        temp.rename(columns={adj_temp_date:'MIN_TEMP_'+ period}, inplace=True)

        j += 2

    output_df = temp.copy()
    
    return output_df
import warnings
warnings.filterwarnings('ignore')
end_date = datetime(2021, 2, 27)
week_counter = 0
df_list = []

while end_date > datetime(2020, 3, 29):
    
    start_time = time.time()
    
    df_list.append(add_lagged_data(end_date))
    
    end_time = time.time()
    print('Feature generation for the week ending in {} finished in {} seconds'.format(end_date.strftime('%Y-%m-%d'),
                                                                                       round(end_time-start_time,1)))
    
    end_date -= timedelta(weeks=1)
    week_counter+=1
final_df = pd.concat(df_list, axis=0)
# check number of weeks for which we have features
final_df.shape[0]/3103<jupyter_output><empty_output><jupyter_text>Imputing NAs<jupyter_code># pd.set_option('max_rows', 200)
final_df.isna().sum().sum()
pd.set_option('max_rows', 400)
final_df.isna().sum()
data_to_save = final_df.copy()
na_columns = data_to_save.columns[data_to_save.isnull().any()]
# impute missing values by state average
for col in data_to_save[na_columns].columns:
    data_to_save[col] = data_to_save.groupby('STATEFP')[col].transform(lambda x: x.fillna(x.mean()))
data_to_save.isna().sum().sum()
times = ['T', 'T_1', 'T_2', 'T_3', 'T_4']

for period in times:
    data_to_save['LOG_DELTA_INC_RATE_' + period] = np.log(data_to_save['DELTA_INC_RATE_' + period] + 1)
    data_to_save['LOG_MEAN_INC_RATE_' + period] = np.log(data_to_save['MEAN_INC_RATE_' + period] + 1)
    data_to_save['LOG_DELTA_SPC_' + period] = np.log(data_to_save['DELTA_SPC_' + period] + 1)
    data_to_save['LOG_DELTA_FPC_' + period] = np.log(data_to_save['DELTA_FPC_' + period] + 1)
final_cols = [
'GEOID',
'NAME',
'State_Name',
'STATEFP', 
'COUNTYFP', 
'date_end_period',
'date_start_period',
'date_end_lag',
'date_start_lag',
'LOG_DELTA_INC_RATE_T',
'PCT_MALE',
'PCT_BLACK',
'PCT_HISPAN', 
'PCT_AMIND',
'PCT_RURAL',
'PCT_COL_DE' ,
'PCT_TRUMP_',
'MED_HOS_IN',
'POPULATION',
'DELTA_CASE_JH_T',
'DELTA_CASE_JH_SMOOTH_T'
]
additional_cols=[
'LOG_DELTA_INC_RATE_T_',
'DELTA_REL_MOB_T_',
'DELTA_RATIO_MOB_T_',
'DELTA_SPC_T_',
'DELTA_SPC_LOGGED_T_',
'DELTA_FPC_T_',
'DELTA_FPC_LOGGED_T_',
'LOG_MEAN_INC_RATE_T_',
'MEAN_REL_MOB_T_',
'MEAN_RATIO_MOB_T_',
'MEAN_FPC_T_',
'MEAN_SPC_T_',
'SLOPE_RATIO_MOB_T_',
'SLOPE_REL_MOB_T_',
'MAX_TEMP_T_',
'MIN_TEMP_T_',
'pct_completely_home_device_count_current_T_',
'pct_full_time_work_behavior_devices_current_T_',
'pct_part_time_work_behavior_devices_current_T_',
'pct_delivery_behavior_devices_current_T_',
'distance_traveled_from_home_current_T_',
'median_home_dwell_time_current_T_',
'pct_completely_home_device_count_baselined_T_',
'pct_full_time_work_behavior_devices_baselined_T_',
'pct_part_time_work_behavior_devices_baselined_T_',
'pct_delivery_behavior_devices_baselined_T_',
'distance_traveled_from_home_baselined_T_',
'median_home_dwell_time_baselined_T_',
'pct_completely_home_device_count_slope_T_',
'pct_full_time_work_behavior_devices_slope_T_',
'pct_part_time_work_behavior_devices_slope_T_',
'pct_delivery_behavior_devices_slope_T_',
'distance_traveled_from_home_slope_T_',
'median_home_dwell_time_slope_T_',
'DELTA_CASE_JH_T_',
'MEAN_SPC_LOGGED_T_',
'MEAN_FPC_LOGGED_T_'
]

for i in range(1,5):
    for col in additional_cols:
        final_cols.append(col+str(i))
    
output = data_to_save[final_cols]
output.to_csv('./output/all_features_updated_incidence.csv')<jupyter_output><empty_output><jupyter_text>Creating dataframes for 2, 3, and 4-week predictionsin this dataframe, the target variables is the the number of cumulative cases in 2, 3, and 4 weeks ahead, denoted by `LOG_DELTA_INC_RATE_T_14`, `LOG_DELTA_INC_RATE_T_21`, and `LOG_DELTA_INC_RATE_T_28` respectively.<jupyter_code># Main function
# the input to the funtion is the end date for which we want to add data

from datetime import datetime, timedelta


def add_lagged_y(date):
    global output, jh_covid_df
    
    T_end = date
    T_start = T_end - timedelta(days=7)
    
    T_start_period = (T_end - timedelta(days=6)).strftime('%Y-%m-%d')
    
    T_14 =  T_end + timedelta(days=7)
    T_21 =  T_end + timedelta(days=14)
    T_28 =  T_end + timedelta(days=21)
    

    dates_non_str = [T_end, T_start, T_14, T_21, T_28]
    
    dates = [item.strftime('%Y-%m-%d') for item in dates_non_str]
    
    dates_jh = [item.strftime('%#m/%#d/%y') for item in dates_non_str]
    
    
    periods = ['T_end', 'T_start', 'T_14', 'T_21', 'T_28']
    
    temp = output.loc[(output.date_end_period==dates[0]) & (output.date_start_period==T_start_period)].copy()
    
    temp['FIPS'] = temp['GEOID'].astype(int)
    
    #print('check 1 {}'.format(temp.shape))
    temp['target_date_2wk'] = T_14.strftime('%Y-%m-%d')
    temp['target_date_3wk'] = T_21.strftime('%Y-%m-%d')
    temp['target_date_4wk'] = T_28.strftime('%Y-%m-%d')
        
    temp = temp.merge(covid_df_contiguous[['GEOID',*dates]], on='GEOID', how='left')
    
    temp = temp.merge(jh_covid_df[['FIPS',*dates_jh]], on='FIPS', how='left')
    #print('check 2 {}'.format(temp.shape))

    for period, date in zip(periods, dates):
        temp['inc_rate_' + period] = temp[date] / temp['POPULATION'] * 10000


    for period, date in zip(periods[-3:], dates[-3:]):
        temp['DELTA_CASE_SMOOTHED_' + period] = temp[date] - temp[dates[1]]
        temp['DELTA_INC_RATE_' + period] = temp['inc_rate_' + period] - temp['inc_rate_T_start']
        temp['LOG_DELTA_INC_RATE_' + period] = np.log(temp['DELTA_INC_RATE_' + period] + 1)
    
    for period, date in zip(periods[-3:], dates_jh[-3:]):
        temp['DELTA_CASE_JH_' + period] = temp[date] - temp[dates_jh[1]]
    
    temp['DELTA_CASE_JH_T'] = temp[dates_jh[0]] - temp[dates_jh[1]]
        
    
    
    cols = ['target_date_2wk','LOG_DELTA_INC_RATE_T_14', 
            'target_date_3wk','LOG_DELTA_INC_RATE_T_21',
            'target_date_4wk','LOG_DELTA_INC_RATE_T_28' ]
    #print('check 3 {}'.format(temp.shape))
    
    return temp[[*output.columns,'DELTA_CASE_JH_T',
            'target_date_2wk','LOG_DELTA_INC_RATE_T_14', 'DELTA_CASE_SMOOTHED_T_14', 'DELTA_CASE_JH_T_14',
            'target_date_3wk','LOG_DELTA_INC_RATE_T_21', 'DELTA_CASE_SMOOTHED_T_21', 'DELTA_CASE_JH_T_21',
            'target_date_4wk','LOG_DELTA_INC_RATE_T_28', 'DELTA_CASE_SMOOTHED_T_28', 'DELTA_CASE_JH_T_28']]
    
    
   
end_date = datetime(2021, 2, 27)
df_lagged_list = []

while end_date > datetime(2020, 3, 29):
    
    start_time = time.time()
    
    df_lagged_list.append(add_lagged_y(end_date))
    
    end_time = time.time()
    print('Feature generation for the week ending in {} finished in {} seconds'.format(end_date.strftime('%Y-%m-%d'),
                                                                                       round(end_time-start_time,1)))
    end_date -= timedelta(weeks=1)
df_lagged = pd.concat(df_lagged_list, axis=0)
df_lagged.isna().sum().sum()
df_lagged.shape, df_lagged.shape[0]/3103
df_lagged.to_csv('./output/all_features.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13847734579563156, 'fraction_numerical': 0.014848695767458817, 'mean_word_length': 3.8391339214113875, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 54, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '605228', 'license': 'MIT', 'n_tokens_mistral': 11755, 'n_tokens_neox': 11139, 'n_words': 1984, 'path': 'Data Preparation.ipynb', 'repo_name': 'geohai/COVID19-STGXB', 'size': '30171'}