<jupyter_start><jupyter_text>Determining corrections to the number of positivesThe number of daily positives $P_i$, has to be scaled in relation to the number of daily tests $T_i$ as $(T_i/\bar{T})^{-\gamma}$, as the number of positives does not grow linearly with number of tests.Furthermore, variation with weekday remains, which should be subsequently corrected for.The following code determins from the (PCR) test data:- Measuring the scaling exponent of positives with respect to number of tests- Determining further weekday correctionsThus, the below code measures this exponent, as the current value of $\gamma = 0.7$ was determined several months ago (by Uffe from SSI), and may have changed. Also a measurement of this value for the AntiGen test is in order. The estimates should be done for both the entire period (e.g. 1st of September 2020 - 1st of March 2021) and possibly also monthly to test, if there are significant variations.The analysis is related to the observation, that the fraction of positives usually increases on Mondays (i.e. results from the weekend). However, one should **first** determine $\gamma$, and if any weekly variations remain, these should be accounted for **subsequently**.<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import scipy 
import datetime as dt
# from scipy.optimize import curve_fit, minimize
import scipy.stats
from scipy import stats
import iminuit
from iminuit import Minuit
from scipy.optimize import curve_fit

mpl.rcParams['font.size'] = 18      # Set the general plotting font size
plt.rc("font", family="serif")
#plt.rc("text", usetex=True)<jupyter_output><empty_output><jupyter_text>General settings:<jupyter_code>save_plots = True

day_0 = '2020-09-01'
fit_to_today = False
gamma_period = ['2020-09-01', '2021-03-05']

date_start = "2020-09-01"
date_end   = "2021-03-05"<jupyter_output><empty_output><jupyter_text>*** Fit to current test & positive numbers:  Read the data:<jupyter_code>data_posi = pd.read_csv("../Data-Epidemiologiske-Rapport-08032021-p9o8/Municipality_cases_time_series.csv", sep=';', thousands='.', index_col=0)
data_test = pd.read_csv("../Data-Epidemiologiske-Rapport-08032021-p9o8/Municipality_tested_persons_time_series.csv", sep=';', thousands='.', index_col=0)

# All of Denmark
if fit_to_today:
    date_end = data_test.index[-1]
    
Nposi_all = data_posi.loc[date_start:date_end].sum(axis=1)
eNposi_all = np.sqrt(Nposi_all)
Ntest_all = data_test.loc[date_start:date_end].sum(axis=1)
eNtest_all = np.sqrt(Ntest_all)
day  = np.arange(1,len(Nposi_all)+1)

# Number of positive tests and uncertainty, assuming fixed number of daily tests and power law correction:
Power = 0.56         # Power in power law for the increase in positives with increase in tests.
# SystError = 0.000   # To account for larger uncertainties on fraction positives than simply the statistical.

# TP 6. March: Are the above necessary? And the below scaling (yet)?
# TP 8. March: Better naming of....  P -> SP -> SCP

nAveDailyTests_all = Ntest_all.mean()
print("Average number of PCR tests: ", nAveDailyTests_all)
fPos_all = Nposi_all / Ntest_all
nPos_all  = Nposi_all * (Ntest_all / nAveDailyTests_all)**(-Power)
enPos_all = np.sqrt(Ntest_all*fPos_all*(1-fPos_all) ) * (Ntest_all / nAveDailyTests_all)**(-Power)<jupyter_output>Average number of PCR tests:  79765.63440860216<jupyter_text>Define ranges, fitting function, and model:<jupyter_code>Plot_StartDay = 1
Plot_EndDay   = 26

Fit_StartDay  = 0
Fit_EndDay    = 33
day1jan       = np.arange(Plot_StartDay, Plot_EndDay)
day1jan_fit   = np.arange(Fit_StartDay, Fit_EndDay)
Npoints = len(day1jan_fit)

day_extrapol = 20<jupyter_output><empty_output><jupyter_text>Comparison of geometric mean with scaled positivesFirst we define a likelihood for comparingf geometric mean of scaled positives and scaled positives with variable scaling, $\gamma$. Then best fit and uncertainty on gamma within time-period is defined<jupyter_code>def model_bllh_og(N_obs, N_exp):
    p_cdf = np.zeros_like(N_exp)
    p_cdf[N_obs < N_exp]  = scipy.stats.poisson.logcdf( N_obs[N_obs < N_exp], N_exp[N_obs < N_exp] )
    p_cdf[N_obs >= N_exp] = scipy.stats.poisson.logsf(  N_obs[N_obs >= N_exp], N_exp[N_obs >= N_exp] )
    return - 2 * sum( (p_cdf) )

def model_bllh(N_obs, N_exp):
    p_cdf = np.zeros_like(N_exp)
    p_cdf  = scipy.stats.poisson.logcdf( N_obs, N_exp )
    p_cdf[p_cdf >= np.log(0.5)] = scipy.stats.poisson.logsf(  N_obs[p_cdf >= np.log(0.5)], N_exp[p_cdf >= np.log(0.5)] )
    return - 2 * sum( (p_cdf) )

#Find likelihood of N_obs given N_exp with poisson statistics

def gamma_index(): #finds best fit gamma and uncertainty for period
    llh_list = []
    gamma_list = np.linspace(0,1.,101) #probe gamma in discrete values. 
    for gamma in gamma_list:
        nPos_all  = Nposi_all * (Ntest_all / nAveDailyTests_all)**(-gamma)
        enPos_all = np.sqrt(Nposi_all * (Ntest_all / nAveDailyTests_all)**(-gamma) )
    
        y = nPos_all[Plot_StartDay:Plot_EndDay]
        y_exp = []
        for i in range(3,len(y)-3):
            y_exp.append( (np.prod(y[i-3:i+4]))**(1/7) ) 
            #Geometric running average, predicted scaled positive test number
            
        llh = model_bllh( np.array(nPos_all[Plot_StartDay+4:Plot_EndDay-2].values), np.array(y_exp) )
        llh_list.append(llh)     
    min_llh = min(llh_list)

    
    p = np.polyfit(gamma_list[np.isclose(llh_list, min_llh, atol=10.)], 
               np.array(llh_list)[np.isclose(llh_list, min_llh, atol=10.)], 2)
    gamma = -p[1]/(2*p[0])
    dgamma = 1/np.sqrt(p[0]) #fit parabola fit intercept with 1 sigma
    
    return gamma, dgamma, llh_list, gamma_list, nPos_all, enPos_all, y_exp

#gamma, dgamma, llh_list, gamma_list, nPos_all, enPos_all, y_exp = gamma_index()<jupyter_output><empty_output><jupyter_text>Compare $\gamma (t)$ Compute $\gamma (t)$ as a linear function using locally computed geometric mean<jupyter_code>gamma_t = []
dgamma_t = []

days_between_points = 7
for i in np.arange(0,len(data_posi.loc[date_start:date_end].sum(axis=1)) - 28,days_between_points): 
    #loop over every days_between_points days from september to february
    Nposi_all = data_posi.loc[date_start:date_end].sum(axis=1)
    eNposi_all = np.sqrt(Nposi_all)
    Ntest_all = data_test.loc[date_start:date_end].sum(axis=1)
    eNtest_all = np.sqrt(Ntest_all)
    day  = np.arange(1,len(Nposi_all)+1)

    Plot_StartDay = i
    Plot_EndDay   = 28+i

    nAveDailyTests_all = Ntest_all.mean()
    fPos_all = Nposi_all / Ntest_all

    gamma, dgamma, llh_list, gamma_list, nPos_all, enPos_all, y_exp = gamma_index()
    gamma_t.append(gamma), dgamma_t.append(dgamma)<jupyter_output><empty_output><jupyter_text>Here we plot $\gamma(t)$<jupyter_code>plt.figure(figsize=(12,8))
plt.errorbar(days_between_points*np.arange(len(gamma_t)), gamma_t, yerr=dgamma_t, fmt='.', capsize=5)
days = days_between_points*np.arange(len(gamma_t)) 
plt.ylabel(r'Scaling Index [$\gamma$(t)]')#, plt.xlabel('Days from '+date_start)
plt.ylim(0,1)#, plt.legend() 

#fit linear function
def lin(x,a,b):
    return a*x+b
p, cov = curve_fit(lin, days, gamma_t, sigma=dgamma_t)
power = lin(np.arange(3,len(Nposi_all[date_start:date_end])-4), *p)

#Heavyside
#power = []
#for idx,i in enumerate(Ntest_all[date_start:date_end][3:-4].index):
#    if idx < 90:
#        power.append(0.78)
#    else:
#        power.append(0.57)
plt.plot(Ntest_all[date_start:date_end][3:-4].index, power, label=f'Linear Fit, $\gamma (t) = a \cdot t + b$')    

#This creates weekly lines
now = dt.date(year=int(date_end[:4]), month = int(date_end[5:7]), day = int( date_end[8:] ))
then = now + dt.timedelta(days=-(len(data_posi.loc[date_start:date_end].sum(axis=1)))+14 )
dato_start = "2021-09-01"
for x in np.arange(0,165,7): #weeks
    plt.axvline( x, ls='--', color='grey', alpha=0.2)
#for x in np.arange(0,150,28): #month 
#    plt.axvline( x, ls='--', color='k', alpha=0.5)
    
# Plot monthly lines on first date
d_now = dt.date(year=int(date_start[:4]), month = int(date_start[5:7]), day = int( date_start[8:] ))
for i in range(len(data_posi.loc[date_start:date_end].sum(axis=1))):
    d_now = d_now + dt.timedelta(days=1)
    if d_now.day == 1:
        plt.axvline( i-12, ls='-', color='k', alpha=0.5)
    
times = []; N_weeks = int(len(data_posi.loc[date_start:date_end].sum(axis=1))/14)
for i in range(N_weeks):
    times.append(then + dt.timedelta(days=14*i))
    
plt.xticks(ticks=[14*i for i in range(N_weeks)], 
           labels=[str(time.year)+'-'+str(time.month)+'-'+str(time.day) for time in times], rotation=30 );

plt.savefig('gamma_t_2.png', dpi=300)
plt.figure(figsize=(20,10))
plt.errorbar(Nposi_all.index, Nposi_all[date_start:date_end], 
             yerr = eNposi_all[date_start:date_end], fmt='.', capsize=3, label='Positive Tests')

#includes date on x-label
ax = plt.gca()
temp = ax.xaxis.get_ticklabels()
temp = list(set(temp) - set(temp[::14]))
for label in temp:
    label.set_visible(False)
for x in np.arange(0,len(Nposi_all[date_start:date_end]),14): #weeks
    plt.axvline( x, ls='--', color='grey', alpha=0.5)
    
    
mean_arr = []
for i in np.arange(3,len(Nposi_all[date_start:date_end])-4): 
    mean_arr.append( np.mean(Ntest_all[date_start:date_end][i-3:i+4]) )
SP = Nposi_all[3:-4] * (Ntest_all[date_start:date_end][3:-4]/mean_arr)**(-np.asarray(power) )


plt.errorbar(SP.index, SP[date_start:date_end], 
             yerr = np.sqrt(SP[date_start:date_end]), fmt='.', capsize=3, label='Scaled Positives')
plt.ylabel(r'$N_{Pos} $')
plt.xticks(rotation = 40); plt.legend();

plt.savefig('Scaled_pos.png', dpi=300)<jupyter_output><empty_output><jupyter_text>Given $\gamma (t)$ compute the locally scaled positive test rate and plot residuals.Week_coor compares scaled positive to total positives to show weekly bias.<jupyter_code># THIS IS THE OLD THING - NOT TO BE USED!!!
"""
def week_corr(dato_start = "2020-08-28", dato_end   = "2020-12-01", power=0.56):
    Nposi_all = data_posi.loc[dato_start:dato_end].sum(axis=1)
    eNposi_all = np.sqrt(Nposi_all)
    Ntest_all = data_test.loc[dato_start:dato_end].sum(axis=1)
    eNtest_all = np.sqrt(Ntest_all)
    mean_arr = []
    
    #Calculate mean of tested in week
    for i in np.arange(3,len(Nposi_all)-4): 
        mean_arr.append( np.mean(Ntest_all[i-3:i+4]) )
    
    #Scale Positive tested
    SP = Nposi_all[3:-4] * (Ntest_all[3:-4]/mean_arr)**(-power)
    
    #Fractional correction by weekday
    n_days = (len(np.arange(3,len(Nposi_all)-4)) % 7)
    nweeks = int((len(np.arange(3,len(Nposi_all)-4))-n_days) / 7)
    if n_days == 0: 
        SP_week = np.reshape(SP.values, (nweeks, 7))
        Nposi_week = np.reshape(Nposi_all[3:-4].values, (nweeks, 7))
    else:
        SP_week = np.reshape(SP.values[:-(len(np.arange(3,len(Nposi_all)-4)) % 7)], (nweeks, 7))
        Nposi_week = np.reshape(Nposi_all[3:-4].values[:-(len(np.arange(3,len(Nposi_all)-4)) % 7)], (nweeks, 7))
        
    mean, spread = np.mean((SP_week-Nposi_week)/SP_week, axis=0), np.std((SP_week-Nposi_week)/SP_week, axis=0)/np.sqrt(nweeks-1)
    
    return mean, spread

mean, spread = week_corr(dato_start = "2020-09-01", dato_end   = "2020-12-15", power=0.8)
"""
def fractional_dev(mean, spread, date_start, dato_end):
    roll = dt.date(int(date_start[:4]), int(date_start[5:7]), int(date_start[8:]) ).weekday() + 3   
    
    plt.figure(figsize=(10, 7))
    plt.errorbar(np.arange(len(mean)), np.roll(mean, roll), yerr=np.roll(spread,roll), 
                 fmt='.', capsize=10,  color='b', label=date_start+' to '+dato_end)
    #plt.errorbar(np.arange(len(mean)), mean, yerr=spread, fmt='.', capsize=10,  color='r')
    
    plt.legend()
    plt.ylabel('Fractional Deviation')
    plt.axhline(0, ls='--', color='k')
    plt.xticks([0,1,2,3,4,5,6], ['Monday','Tuesday','Wednesday','Thursday','Friday','Saturday','Sunday'], rotation=25);
    
#plt.savefig('Weekly_Oscilation.png', dpi=300)
# fractional_dev(mean, spread, "2020-09-01", date_end)

# plt.savefig('Fractional_deviation_wo_gamma.png', dpi=300)
def week_corr_sp_sys(dato_start = "2020-08-28", dato_end   = "2020-12-01", power=0.56):
    Nposi_all = data_posi.loc[dato_start:dato_end].sum(axis=1)
    eNposi_all = np.sqrt(Nposi_all)
    Ntest_all = data_test.loc[dato_start:dato_end].sum(axis=1)
    eNtest_all = np.sqrt(Ntest_all)
    mean_arr = []
    
    # Calculate mean of tested in week:
    for i in np.arange(3,len(Nposi_all)-4): 
        mean_arr.append( np.mean(Ntest_all[i-3:i+4]) )
    
    # Scale Positive tested in week:
    SP = Nposi_all[3:-4] * (Ntest_all[3:-4]/mean_arr)**(-power)
    SP_exp = np.zeros_like(SP)
    
    # Geometric mean of the Scaled Positives:
    for i in np.arange(3,len(SP)-4): 
        SP_exp[i] = (np.prod(SP[i-3:i+4]))**(1/7)
    
    # Fractional correction by weekday:
    n_days = (len(np.arange(3,len(SP_exp)-4)) % 7)
    nweeks = int((len(np.arange(3,len(SP_exp)-4))-n_days) / 7)
    if n_days == 0: 
        SP_week = np.reshape(SP_exp[3:-4], (nweeks, 7))
        Nposi_week = np.reshape(SP[3:-4].values, (nweeks, 7))
    else:
        SP_week = np.reshape(SP_exp[3:-4][:-(len(np.arange(3,len(SP)-4)) % 7)], (nweeks, 7))
        Nposi_week = np.reshape(SP[3:-4].values[:-(len(np.arange(3,len(SP)-4)) % 7)], (nweeks, 7))
        
    mean, spread = np.mean((SP_week-Nposi_week)/SP_week, axis=0), np.std((SP_week-Nposi_week)/SP_week, axis=0)/np.sqrt(nweeks-1)
    
    return mean, spread, SP, SP_exp, Nposi_all

mean1, spread1, SP, SP_exp, Nposi_all = week_corr_sp_sys(dato_start = date_start, dato_end = date_end)

fractional_dev(mean1, spread1, "2020-09-01", date_end)
plt.savefig('Fractional_deviation_wo_gamma.png', dpi=300)
mean, spread, SP, SP_exp,Nposi_all = week_corr_sp_sys(dato_start = '2021-01-01', dato_end   = date_end)

fractional_dev(mean, spread, "2021-01-01", date_end)
mean, spread, SP, SP_exp,Nposi_all = week_corr_sp_sys(dato_start = '2020-09-01', dato_end   = '2020-12-15')

fractional_dev(mean, spread, "2021-09-01", '2020-12-15')
mean
spread
np.sum( (mean/spread)**2)
stats.chi2.sf( 24, 7)  
np.sum( (mean-np.roll(mean1, -3))**2/(spread**2) )
(mean-np.roll(mean1, -3))**2/(spread**2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15188984122523536, 'fraction_numerical': 0.035689194885485456, 'mean_word_length': 4.394088669950739, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 40, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '287540', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6027, 'n_tokens_neox': 5591, 'n_words': 1133, 'path': 'MiscCode/MeasuringTestScalingExponent.ipynb', 'repo_name': 'troelspetersen/Covid19fitting', 'size': '14234'}