# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
# Load libraries and functions
# %load_ext autoreload
# %autoreload 2
# %matplotlib inline
RANDOM_STATE = 42  # Pseudo-random state

from utils import *
sns.set_palette("tab10") # Default seaborn theme

# Extra libraries for this notebook

# +
# Upload dataset
fn_vae_data = glob.glob('./Updated*.pkl')
latest_fn_vae_data = max(fn_vae_data, key=os.path.getctime)

print("Loading... ",latest_fn_vae_data)
with open(latest_fn_vae_data, "rb") as f:
    vae_data_main = pickle.load(f)
print("Done")

# -

# # Daily prevalence of HARTI

# +
# Define function

def plot_regression_line(ts_infections, ts_total_people, a, b, c):
    data = pd.DataFrame(ts_infections * a/ ts_total_people).reset_index()
    data.columns = ['date', 'value']
    # Fit a polyfit model
    x = np.arange(data['date'].size) 
    slope, intercept = np.polyfit(x, data['value'], 1)
    fit = np.polyfit(x, data['value'], 1) # first degree polynomial
    fit_fn = np.poly1d(fit)
    # Plot scatter plot and regression model
    plt.figure(figsize=(10, 5))
    plt.plot(data['date'], data['value'], 'bo', ms=2, marker='x')
    plt.plot(data['date'], data['value'].rolling(window=30, center=True).mean())
    plt.plot(data['date'], fit_fn(x), 'k-')
#     fit = np.polyfit(x, data['value'], 9) # 9th degree polynomial
#     fit_fn = np.poly1d(fit)
#     plt.plot(data['date'], fit_fn(x), 'k-', color='green')
    # Add data decomposition to find a trend
    ts_ratio = (ts_infections / ts_total_people)
    decomposition = sm.tsa.seasonal_decompose(pd.DataFrame(ts_ratio), freq=365, model='additive')
    plt.plot(decomposition.trend*a, color='red')
    plt.grid(linestyle='dotted')
    plt.ylabel('Daily prevalence, %')
    plt.ylim(b, c)
    print('Slope (change per day): ', slope)
    print('Slope (change per year): ', slope*365)
    print('Slope (change per study period): ', slope*365*8)
    plt.legend(('one day value', 'smoothed monthly average', 'fitted linear model',
#                 '9th-degree polynomial',
                'annual trend'),
               fontsize=12, loc='upper right')
    
def get_data_for_regression(ts_infections, ts_total_people, a):
    data = pd.DataFrame(ts_infections * a/ ts_total_people).reset_index()
    data.columns = ['date', 'value']

    x = np.arange(data['date'].size) 
    return x, data['value']


# -

# ### Prevalence of VA-HARTI

# +
# Plot prevalence of VA-HARTI
ts_infections = vae_data_main.groupby('date').sum()['vap']
ts_total_people = vae_data_main.groupby('date').count()['ID']
ts_ratio = (ts_infections / ts_total_people)

plot_regression_line(ts_infections, ts_total_people, 100, -2, 70)
plt.title('Daily prevalence of VA-HARTI in ICU patients in 2011-2020')
plt.tight_layout()
plt.savefig('./pictures/prevalence_vap.pdf', dpi=600)

# +
# Linear regression model

# Define data
a = pd.DataFrame(ts_ratio)
a1 = a.index.values.astype('datetime64[D]')
x = a1.astype('int')
x = sm.add_constant(x)
y = a[0].values

# Fit OLS model
model = sm.OLS(y, x)
results = model.fit()
print(results.summary())

###### Check the quality of linear model
# Plot residuals
x = a1.astype('int')
sns.residplot(x, y)
###### Seemingly uniformal distribution of residuals -> can use linear model ######

# Plot QQ plot
res = results.resid # residuals
sm.qqplot(res, stats.t, fit=True, line='45')
##### Lightly tailed QQ plot -> can use linear model ######

# Get goodness-of-fit for polynomial regression
x = sm.add_constant(x)
polynomial_features= PolynomialFeatures(degree=9)
xp = polynomial_features.fit_transform(x)
model = sm.OLS(y, xp).fit()
#model.summary()
print('\nR squared for 9th-degree polynomial: ', model.rsquared)

# +
# print results
ts_ratio_2011 = ts_ratio[((ts_ratio.index >= '2011') & (ts_ratio.index < '2012'))]
ts_ratio_2020 = ts_ratio[((ts_ratio.index >= '2020') & (ts_ratio.index < '2021'))]

print('Mean prevalence: ', ts_ratio.mean())
print('Mean prevalence 2011: ', ts_ratio_2011.mean())
print('Mean prevalence 2020: ', ts_ratio_2020.mean())
print('Lowest prevalence, 2016: ', ts_ratio[((ts_ratio.index >= '2016') & (ts_ratio.index < '2017'))].mean())
print('Highest prevalence, 2012: ', ts_ratio[((ts_ratio.index >= '2012') & (ts_ratio.index < '2013'))].mean())

# Welch’s t-test, which does not assume equal population variance
print('Welch ttest p-value: ', ttest_ind(ts_ratio_2011, ts_ratio_2020, equal_var=False)[1])
# -

# ### Prevalence of NVA-HARTI

# +
# Plot prevalence of VA-HARTI
ts_infections = vae_data_main.groupby('date').sum()['non_vap_resp_hai']
ts_total_people = vae_data_main.groupby('date').count()['ID']
ts_ratio = (ts_infections / ts_total_people)

plot_regression_line(ts_infections, ts_total_people, 100, -2, 70)
plt.title('Daily prevalence of NVA-HARTI in ICU patients in 2011-2020')
plt.tight_layout()
plt.savefig('./pictures/prevalence_nva.pdf', dpi=600)

# +
# Linear regression model

# Define data
a = pd.DataFrame(ts_ratio)
a1 = a.index.values.astype('datetime64[D]')
x = a1.astype('int')
x = sm.add_constant(x)
y = a[0].values

# Fit OLS model
model = sm.OLS(y, x)
results = model.fit()
print(results.summary())

###### Check the quality of linear model
# Plot residuals
x = a1.astype('int')
sns.residplot(x, y)
###### Non-uniformal distribution of residuals -> can't use linear model!!! ######

# Plot QQ plot
res = results.resid # residuals
sm.qqplot(res, stats.t, fit=True, line='45')
##### Disturbed QQ plot suggests right-skewed distribution-> can't use linear model!!! ######

# Get goodness-of-fit for polynomial regression
x = sm.add_constant(x)
polynomial_features= PolynomialFeatures(degree=9)
xp = polynomial_features.fit_transform(x)
model = sm.OLS(y, xp).fit()
#model.summary()
print('\nR squared for 9th-degree polynomial: ', model.rsquared)

# +
# print results
ts_ratio_2011 = ts_ratio[((ts_ratio.index >= '2011') & (ts_ratio.index < '2012'))]
ts_ratio_2020 = ts_ratio[((ts_ratio.index >= '2020') & (ts_ratio.index < '2021'))]

print('Mean prevalence: ', ts_ratio.mean())
print('Mean prevalence 2011: ', ts_ratio_2011.mean())
print('Mean prevalence 2021: ', ts_ratio_2020.mean())
print('Lowest prevalence, 2016: ', ts_ratio[((ts_ratio.index >= '2016') & (ts_ratio.index < '2017'))].mean())
print('Highest prevalence, 2017: ', ts_ratio[((ts_ratio.index >= '2017') & (ts_ratio.index < '2018'))].mean())

# Welch’s t-test, which does not assume equal population variance
print('Welch ttest p-value: ', ttest_ind(ts_ratio_2011, ts_ratio_2020, equal_var=False)[1])
# -

# # Cumulative proportional incidence of HARTI

# +
# Overall cumulative incidence

ir = vae_data_main[['ID_subid', 'infection_respiratory']].groupby('ID_subid').max().sum()
cases = vae_data_main[['ID_subid']].groupby('ID_subid').max().shape[0]
print('Overall incidence all HARTI: ', (ir/cases).values)
print('\n95% CI all HARTI: ', ci(ir, cases)[0].values, ci(ir, cases)[1].values)

vap = vae_data_main[['ID_subid', 'vap']].groupby('ID_subid').max().sum()
print('\nOverall incidence VA-HARTI: ', (vap/cases).values)
print('\n95% CI, VA-HARTI: ', ci(vap, cases)[0].values, ci(vap, cases)[1].values)

nonvap = vae_data_main[['ID_subid', 'non_vap_resp_hai']].groupby('ID_subid').max().sum()
print('\nOverall incidence NVA-HARTI: ', (nonvap/cases).values)
print('\n95% CI NVA-HARTI: ', ci(nonvap, cases)[0].values, ci(nonvap, cases)[1].values)

# +
# Cumulative incidence by years
inc = {}
for col in ['vap', 'non_vap_resp_hai']:
    inc[col] = {}
    a = vae_data_main[['ID_subid', 'year', col]].groupby('ID_subid').max()
    x = a.groupby('year').sum()
    x = x.head(10)
    n = a.groupby('year').count()
    n = n.head(10)
    
    inc[col]['inc'] = (x / n).values.reshape(-1)
    inc[col]['ci_low'] = ci(x, n, alpha=0.05, method='normal')[0].values.reshape(-1)
    inc[col]['ci_up'] = ci(x, n, alpha=0.05, method='normal')[1].values.reshape(-1)

fig, ax = plt.subplots(1, figsize=(10,5.5))
for key, val in inc.items():
    val_df = pd.DataFrame(val) * 100
    val_df['inc'].plot(style='-o', linewidth=1.5, ax=ax)
    ax.fill_between(np.arange(val_df.shape[0]), val_df['ci_low'], val_df['ci_up'], alpha=.2)


plt.legend(labels = ['VA-HARTI *', 'NVA-HARTI'])
plt.title('Cumulative incidence of HARTI per 100 ICU admissions in 2011-2020')
plt.ylabel('Yearly cumulative incidence per 100 ICU admissions')
plt.xlabel('')

xlabels = ('2011', '2012', '2013', '2014', '2015', '2016', '2017', '2018', '2019', '2020')
# create an index for each tick position
xi = list(range(len(xlabels)))
plt.xticks(xi, xlabels)
ax.yaxis.set_ticks_position('both')
ax.tick_params(labeltop=False, labelright=True)
plt.minorticks_on()
plt.grid(linestyle='dotted', which='both', axis='y')
nv = round(((nonvap/cases).values *100).item(0), 1)
v = round(((vap/cases).values *100).item(0), 1)
plt.axhline(y=nv, color='black', linestyle='dotted')
plt.text(-0.2, nv - 2.8, 'mean NVA-HARTI\nincidence = %s'%(nv), fontsize=12)
plt.axhline(y=v, color='black', linestyle='dotted')
plt.text(-0.2, v - 2.8, 'mean VA-HARTI\nincidence = %s'%(v), fontsize=12)
plt.tight_layout()
plt.savefig('pictures/incidence_years.pdf', bbox_inches="tight", dpi=600)

# Adfuller test for stationarity
print('VA-HARTI, Adfuller test: ', adfuller(inc['vap']['inc'])[1])
print('NVA-HARTI, Adfuller test: ', adfuller(inc['non_vap_resp_hai']['inc'])[1])
# -

# ### Test trends in incidence by linear regression

# +
# Fit linear regression for VA-HARTI incidence
# Define data
x = [2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020]
x = sm.add_constant(x)
y = inc['vap']['inc']

# Fit OLS model
model = sm.OLS(y, x)
results = model.fit()
print(results.summary())

# +
# Fit linear regression for NVA-HARTI incidence
# Define data
x = [2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020]
x = sm.add_constant(x)
y = inc['non_vap_resp_hai']['inc']

# Fit OLS model
model = sm.OLS(y, x)
results = model.fit()
print(results.summary())
# -

#Print results
print("VA-HARTI, yearly incidence rate: ", inc['vap']['inc'])
print("\nVA-HARTI, 95% CI: ", inc['vap']['ci_low'], inc['vap']['ci_up'])
print("\nNVA-HARTI, yearly incidence rate: ", inc['non_vap_resp_hai']['inc'])
print("\nNVA-HARTI, 95% CI: ", inc['non_vap_resp_hai']['ci_low'], inc['non_vap_resp_hai']['ci_up'])

# # VA- and NVA-HARTI rates per patient-days

# +
# get first and last day in the ICU for each patient
max_date_per_subid = vae_data_main[['ID_subid', 'date']].groupby('ID_subid').max()['date']
min_date_per_subid = vae_data_main[['ID_subid', 'first_day_in_icu']].groupby('ID_subid').first()['first_day_in_icu']

# patient-days per quarter
patients_per_day = {}
index = np.arange(min_date_per_subid.min(), max_date_per_subid.max(), datetime.timedelta(days=1))

for i, day in enumerate(index):
    patients_per_day[day] = ((min_date_per_subid <= day) & (max_date_per_subid >= day)).sum()
    
patients_per_day = pd.DataFrame(pd.Series(patients_per_day))
month_from_year = (patients_per_day.index.year - patients_per_day.index.year.min()) * 12
month_from_month = patients_per_day.index.month
patients_per_day['quarter'] = ((month_from_month + month_from_year) // 3).astype('int')

# Select group data for VA- and NVA-HARTI
data_vap = vae_data_main.loc[(vae_data_main.group == 'VA-HARTI')]
data_nonvap = vae_data_main.loc[(vae_data_main.group == 'NVA-HARTI')]

dsets = {'vap': data_vap, 
         'nonvap': data_nonvap}

# VA-HARTI cases per quarter
vap_date = dsets['vap'][['ID_subid', 'harti_first_date']].groupby('ID_subid').first().values.astype('datetime64')
vap_per_day = {}
index = np.arange(min_date_per_subid.min(), max_date_per_subid.max(), datetime.timedelta(days=1))

for day in index:
    vap_per_day[day] = (day == vap_date).sum()
    
vap_per_day = pd.DataFrame(pd.Series(vap_per_day))
vap_per_day['quarter'] = ((month_from_month + month_from_year) // 3).astype('int')

# NVA-HARTI cases per quarter
nva_date = dsets['nonvap'][['ID_subid', 'harti_first_date']].groupby('ID_subid').first().values.astype('datetime64')

nva_per_day = {}
index = np.arange(min_date_per_subid.min(), max_date_per_subid.max(), datetime.timedelta(days=1))

for day in index:
    nva_per_day[day] = (day == nva_date).sum()
    
nva_per_day = pd.DataFrame(pd.Series(nva_per_day))
nva_per_day['quarter'] = ((month_from_month + month_from_year) // 3).astype('int')

# Incidence per quarter
inc_vap = []
for i in range(4):
    a = (vap_per_day.groupby('quarter').sum()[i::4][0] / patients_per_day.groupby('quarter').sum()[i::4][0]*1000)
    inc_vap.append(a)
inc_vap = pd.DataFrame(inc_vap)
inc_vap = inc_vap.T.apply(lambda x: (x.dropna()), axis=1)
inc_vap = inc_vap[:-1]

inc_nva = []
for i in range(4):
    a = (nva_per_day.groupby('quarter').sum()[i::4][0] / patients_per_day.groupby('quarter').sum()[i::4][0]*1000)
    inc_nva.append(a)
inc_nva = pd.DataFrame(inc_nva)
inc_nva = inc_nva.T.apply(lambda x: (x.dropna()), axis=1)
inc_nva = inc_nva[:-1]

# +
# Prepare long dataset
inc = pd.concat([inc_vap, inc_nva], axis=1)
inc.columns = ('VA-HARTI', 'NVA-HARTI')
inc['quarter'] = inc_vap.index
inc_long = pd.melt(inc, id_vars='quarter')

# fit linear regression
seaborn_grid = sns.lmplot(x='quarter', y='value', hue="variable", data=inc_long, robust=True, legend=False)
seaborn_grid.fig.set_figwidth(7)
plt.grid(linestyle='dotted')
plt.ylabel("Incidence rate per 1000 patient-days")
xlabels = ('2011', '2012', '2013', '2014', '2015', '2016', '2017', '2018', '2019', '2020')
# create an index for each tick position
xi = list(range(len(inc)))
xi = xi[0::4]
plt.xticks(xi, xlabels)
plt.xlabel("")
plt.ylim(0, 21)
plt.title("Incidence rate of VA- and NVA-HARTI per 1000 patient-days in the ICU. Robust regression")
plt.legend(ncol=2, mode="expand")

# Get p-value from linear regression, VA-HARTI
x = inc.quarter.values
x = sm.add_constant(x)
y = inc['VA-HARTI'].values
results_va = sm.OLS(y, x).fit()
pval_va = results_va.pvalues[1]
text_va = inc['VA-HARTI'].mean()
plt.text(10, 18, 'VA-HARTI mean rate = '+ "%.1f" % text_va + ' (+/-'+ "%.1f" % inc['VA-HARTI'].std() + '); p-value = '+ "%.4f" % pval_va, fontsize=10)

# Get p-value from linear regression, NVA-HARTI
x = inc.quarter.values
x = sm.add_constant(x)
y = inc['NVA-HARTI'].values
results_nva = sm.OLS(y, x).fit()
pval_nva = results_nva.pvalues[1]
text_nva = inc['NVA-HARTI'].mean()
plt.text(10, 16.5, 'NVA-HARTI mean rate = '+ "%.1f" % text_nva + ' (+/-'+ "%.1f" % inc['NVA-HARTI'].std() + '); p-value = '+ "%.4f" % pval_nva, fontsize=10)

plt.savefig('./pictures/incidence_days_VA_NVA.pdf', bbox_inches="tight", dpi=600)

# Adfuller test for stationarity
print('VA-HARTI per 1000 pts-days, Adfuller test: ', adfuller(inc['VA-HARTI'])[1])
print('NVA-HARTI per 1000 pts-days, Adfuller test: ', adfuller(inc['NVA-HARTI'])[1])

# -

# # VA-HARTI rates per ventilator-days

# +
# Prepare data
# a/b
# a numerator
harti_first_dates = vae_data_main.loc[vae_data_main.vap == 1, ['ID_subid', 'harti_first_date']].groupby('ID_subid').first()
quarters = ((harti_first_dates.harti_first_date.dt.year - harti_first_dates.harti_first_date.dt.year.min()) * 12 + harti_first_dates.harti_first_date.dt.month) // 3
quarters = quarters.astype(int)
numer = harti_first_dates.groupby(quarters).count()

# b denominator
quarters = ((vae_data_main.date.dt.year - vae_data_main.date.dt.year.min()) * 12 + vae_data_main.date.dt.month) // 3
quarters = quarters.astype(int)
denom = vae_data_main[['mech_vent']].groupby(quarters).sum()
denom = denom[:-1]

case_per_ventilator_days = (pd.DataFrame(numer.to_numpy().reshape(-1) / denom.to_numpy().reshape(-1)) * 1000)

# Prepare long dataset
case_per_ventilator_days['quarter'] = case_per_ventilator_days.index
inc_vap_long = pd.melt(case_per_ventilator_days, id_vars='quarter')

# +
# Plot incidence rate
# fit a nonparametric regression using a lowess smoother
seaborn_grid = sns.lmplot(x='quarter', y='value', data=inc_vap_long, robust=True)
seaborn_grid.fig.set_figwidth(7)
plt.grid(linestyle='dotted')
plt.ylabel("Incidence rate per 1000 ventilator-days")
plt.xlabel("")
xlabels = ('2011', '2012', '2013', '2014', '2015', '2016', '2017', '2018', '2019', '2020')
# create an index for each tick position
xi = list(range(len(inc)))
xi = xi[0::4]
plt.xticks(xi, xlabels)
plt.title("VA-HARTI, cases per 1000 ventilator-days in the ICU; quarter data; robust regression")
plt.ylim(0, 30)

# Get p-value from linear regression
x = case_per_ventilator_days.quarter.values
x = sm.add_constant(x)
y = case_per_ventilator_days[0].values
results = sm.OLS(y, x).fit()
pval = results.pvalues[1]
text = case_per_ventilator_days[0].mean()
plt.text(10, 27.5, 'Mean rate = '+ "%.1f" % text + ' (+/-'+ "%.1f" % case_per_ventilator_days[0].std() + '); p-value = '+ "%.2f" % pval, fontsize=12)


plt.savefig('./pictures/incidence_vent-days.pdf', bbox_inches="tight", dpi=600)

print('VA-HARTI per 1000 vent-days, Adfuller test: ', adfuller(inc_vap_long['value'])[1])

# -

# plot residuals
sns.residplot(x='quarter', y='value', data=inc_vap_long)
####### uniform distribution of residuals -> can use linear model

# Show number of vent-days
mech_vent_days = vae_data_main[['ID_subid', 'year', 'mech_vent']]
mech_vent_days = mech_vent_days.groupby('year').sum().astype('int')
mech_vent_days

# _____

StarcoderdataJupyterScriptsDedupFiltered