<jupyter_start><jupyter_text>Covid-19: From model prediction to model predictive control JPGMs tryout of the QALY calculation*Original code by <NAME>. Modified by <NAME> in consultation with the BIOMATH research unit headed by prof. <NAME>.*Copyright (c) 2020 by BIOMATH, Ghent University. All Rights Reserved. This notebook contains a demo calculation of the number of *Quality adjusted life years* lost due to COVID-19. The calculation is performed ex-post simulation. The calculation is a simple manipulation of the number of predicted deaths. For this demo, the deterministic national-level model was used. All details with regard to the model can be found in our preprint *A deterministic, age-stratified, extended SEIRD model for assessing the effect of non-pharmaceutical interventions on SARS-CoV-2 spread in Belgium*, which can be found here: https://www.medrxiv.org/content/10.1101/2020.07.17.20156034v2<jupyter_code>import gc
import sys, getopt
import ujson as json
import random
import datetime
import argparse
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.cm as cm<jupyter_output><empty_output><jupyter_text>Load the covid 19 custom development code<jupyter_code>from covid19model.models import models
from covid19model.data import mobility, sciensano, model_parameters, VOC
from covid19model.models.time_dependant_parameter_fncs import ramp_fun
from covid19model.visualization.output import _apply_tick_locator 
# OPTIONAL: Load the "autoreload" extension so that package code can change
%load_ext autoreload
# OPTIONAL: always reload modules so that as you change code in src, it gets loaded
%autoreload 2
# --------------------------
# Define simulation settings
# --------------------------

N=10
k=1
# Start and end of simulation
start_sim = '2020-09-01'
end_sim = '2022-03-01'
# Confidence level used to visualise model fit
conf_int = 0.05
# Path where figures and results should be stored
fig_path = '../../results/calibrations/COVID19_SEIRD/national/others/WAVE2/'
# Path where MCMC samples should be saved
samples_path = '../../data/interim/model_parameters/COVID19_SEIRD/calibrations/national/'<jupyter_output><empty_output><jupyter_text>Model initialization<jupyter_code># -----------------------
# Load samples dictionary
# -----------------------

from covid19model.models.utils import load_samples_dict
samples_dict = load_samples_dict(samples_path+'BE_WAVE2_stratified_vacc_R0_COMP_EFF_2021-08-31.json', wave=2)
warmup = int(samples_dict['warmup'])

# ---------
# Load data
# ---------

# Time-integrated contact matrices
initN, Nc_all = model_parameters.get_integrated_willem2012_interaction_matrices()
levels = initN.size
# Sciensano public data
df_sciensano = sciensano.get_sciensano_COVID19_data(update=False)
# Sciensano mortality data
df_sciensano_mortality =sciensano.get_mortality_data()
# Google Mobility data
df_google = mobility.get_google_mobility_data(update=False)
# Serological data
df_sero_herzog, df_sero_sciensano = sciensano.get_serological_data()
# VOC data
df_VOC_abc = VOC.get_abc_data()
# Start of data collection
start_data = df_sciensano.idxmin()
# Start of calibration warmup and beta
start_calibration = samples_dict['start_calibration']
# Last datapoint used to calibrate warmup and beta
end_calibration = samples_dict['end_calibration']

# --------------
# Initial states
# --------------

# Model initial condition on September 1st
with open('../../data/interim/model_parameters/COVID19_SEIRD/calibrations/national/initial_states_2020-09-01.json', 'r') as fp:
    initial_states = json.load(fp)
# Correct size of initial states
entries_to_remove = ('S_v', 'E_v', 'I_v', 'A_v', 'M_v', 'C_v', 'C_icurec_v', 'ICU_v', 'R_v')
for i in entries_to_remove:
    initial_states.pop(i, None)
for key, value in initial_states.items():
    initial_states[key] = np.concatenate((np.expand_dims(initial_states[key],axis=1),np.ones([9,2])),axis=1) 

# ---------------------------
# Time-dependant VOC function
# ---------------------------

from covid19model.models.time_dependant_parameter_fncs import make_VOC_function
VOC_function = make_VOC_function()

# ---------------------------------------------------------
# Time-dependant vaccination function and sampling function
# ---------------------------------------------------------

from covid19model.models.time_dependant_parameter_fncs import  make_vaccination_function
vacc_strategy = make_vaccination_function(df_sciensano).stratified_vaccination_strategy
from covid19model.models.utils import draw_fcn_WAVE2_stratified_vacc as draw_fcn_WAVE2

# --------------------------------------
# Time-dependant social contact function
# --------------------------------------

# Extract build contact matrix function
from covid19model.models.time_dependant_parameter_fncs import make_contact_matrix_function, delayed_ramp_fun, ramp_fun
contact_matrix_4prev = make_contact_matrix_function(df_google, Nc_all)
policies_WAVE2 = make_contact_matrix_function(df_google, Nc_all).policies_WAVE2_no_relaxation
    
# -----------------------------------
# Time-dependant seasonality function
# -----------------------------------

def seasonality(t,states,param,amplitude, peak_shift):
    maxdate = pd.Timedelta(days=peak_shift) + pd.to_datetime('2021-01-01')
    t = (t - pd.to_datetime(maxdate))/pd.Timedelta(days=1)/365
    return param*(1+amplitude*np.cos( 2*np.pi*(t)))

# ---------------------------------------------------
# Function to add poisson draws and sampling function
# ---------------------------------------------------

from covid19model.models.utils import output_to_visuals

# --------------------
# Initialize the model
# --------------------

# Load the model parameters dictionary
params = model_parameters.get_COVID19_SEIRD_parameters(vaccination=True)
params['dICUrec'] = [9.9, 3.4, 8.4, 6.6, 8.2, 10.1, 11.5, 15.2, 13.3]
# Add "size dummy" for vaccination stratification
params.update({'doses': np.zeros([3,3])})
# Correct size of other parameters
params.update({'e_s': np.array([[0, 0.5, 0.8],[0, 0.5, 0.8],[0, 0.3, 0.75]])}) # rows = VOC, columns = # no. doses
params.update({'e_h': np.array([[0,0.78,0.92],[0,0.78,0.92],[0,0.75,0.94]])})
params.pop('e_a')
params.update({'e_i': np.array([[0,0.5,0.5],[0,0.5,0.5],[0,0.5,0.5]])})  
params.update({'d_vacc': 31*36})
params.update({'N_vacc': np.zeros([9,3])})
# Social policies
params.update({'l': 21, 'prev_schools': 0, 'prev_work': 0.5, 'prev_rest_lockdown': 0.5, 'prev_rest_relaxation':0.5, 'prev_home': 0.5})
# Seasonality
params.update({'amplitude': 0.1, 'peak_shift': 0})
# Vaccination
params.update(
    {'initN': initN, 'vacc_order': np.array(range(9))[::-1], 'daily_first_dose': 55000,
     'refusal': 0.2*np.ones([9,2]), 'delay_immunity': 14, 'delay_doses': 3*7, 'stop_idx': 0}
)

model = models.COVID19_SEIRD_stratified_vacc(initial_states, params,
                                                time_dependent_parameters={'beta': seasonality,
                                                                           'Nc': policies_WAVE2,
                                                                           'N_vacc': vacc_strategy,
                                                                           'alpha':VOC_function})<jupyter_output><empty_output><jupyter_text>Perform simulation<jupyter_code># -------------------
# Perform simulations
# -------------------

print('\n1) Simulating COVID-19 SEIRD '+str(N)+' times')
start_sim = start_calibration
out = model.sim(end_sim,start_date=start_sim,warmup=warmup,N=N,draw_fcn=draw_fcn_WAVE2,samples=samples_dict)
simtime, df_2plot = output_to_visuals(out, ['H_in', 'H_tot', 'ICU', 'R'], N, k, LL = conf_int/2, UL = 1 - conf_int/2)
deaths_hospital = df_sciensano_mortality.xs(key='all', level="age_class", drop_level=True)['hospital','cumsum']
fig,ax=plt.subplots(figsize=(12,4))
ax.plot(simtime,df_2plot['H_in','mean'], linewidth=1.5, color='blue')
ax.fill_between(simtime,df_2plot['H_in','LL'],df_2plot['H_in','UL'],color='blue',alpha=0.2)
ax.scatter(df_sciensano.index, df_sciensano['H_in'], color='black', alpha=0.4, linestyle='None', facecolors='none', s=60, linewidth=2)
ax.set_xlim(['2020-09-01','2022-03-01'])<jupyter_output><empty_output><jupyter_text>Visualize the total number of patients in the hospital and the predicted number of deaths in the hospital<jupyter_code>print('2) Visualizing fit on deaths')

dates = ['2021-01-01','2021-03-01']

fig,axes = plt.subplots(nrows=len(dates),ncols=1,figsize=(12,4*len(dates)),sharex=True)
for idx,date in enumerate(dates):
    data_sciensano = []
    for jdx,age_group in enumerate(df_sciensano_mortality.index.get_level_values(0).unique().values[1:]):
        data_sciensano.append(df_sciensano_mortality.xs(key=age_group, level="age_class", drop_level=True).loc[dates[idx]]['hospital','cumsum'])
    
    axes[idx].scatter(df_sciensano_mortality.index.get_level_values(0).unique().values[1:],out['D'].sum(dim='doses').mean(dim='draws').loc[dict(time=date)],color='black',marker='v',zorder=1)
    yerr = np.zeros([2,len(out['D'].quantile(dim='draws',q=0.975).loc[dict(time=date)].values)])
    yerr[0,:] = out['D'].sum(dim='doses').mean(dim='draws').loc[dict(time=date)] - out['D'].sum(dim='doses').quantile(dim='draws',q=0.025).loc[dict(time=date)].values
    yerr[1,:] = out['D'].sum(dim='doses').quantile(dim='draws',q=0.975).loc[dict(time=date)].values - out['D'].sum(dim='doses').mean(dim='draws').loc[dict(time=date)]
    axes[idx].errorbar(x=df_sciensano_mortality.index.get_level_values(0).unique().values[1:],
                       y=out['D'].sum(dim='doses').mean(dim='draws').loc[dict(time=date)],
                       yerr=yerr,
                       color = 'black', fmt = '--v', zorder=1, linewidth=1, ecolor='black', elinewidth=1, capsize=5)
    axes[idx].bar(df_sciensano_mortality.index.get_level_values(0).unique().values[1:],data_sciensano,width=1,alpha=0.7,zorder=0)
    axes[idx].set_title('Cumulative hospital deaths on '+date)
    axes[idx].grid(False)
plt.show()<jupyter_output>2) Visualizing fit on deaths<jupyter_text>Comorbidity QALY calculations<jupyter_code># Obtain the number of QALYs that would be lost if a person of a given age dies
from covid19model.models.QALY import *<jupyter_output><empty_output><jupyter_text>Define the comorbidity profile of the studied population<jupyter_code># Implement the comorbidity distribution for coronavirus patients, depending on the disease outcome
# For now, the estimates by Harrison et. al 2020 are applied uniformly in every group
# Ideally, we'd use data by sciensano for this
######################################################################################

columns = [['R','R','R','R','D','D','D','D',],['0','1','2','3+','0','1','2','3+']]
tuples = list(zip(*columns))
columns = pd.MultiIndex.from_tuples(tuples, names=["population", "CCI"])
comorbidity_distribution = pd.DataFrame(index=pd.IntervalIndex.from_tuples([(0,12),(12,18),(18,25),(25,35),(35,45),(45,55),(55,65),(65,75),(75,85),(85,120)], closed='left'), columns=columns)
comorbidity_distribution.index.name = 'age_group'
for idx,key in enumerate(comorbidity_distribution.index):
    comorbidity_distribution.loc[key]['R'] = [0.614, 0.174, 0.076, 0.137]
    comorbidity_distribution.loc[key]['D'] = [0.173, 0.161, 0.144, 0.522]
comorbidity_distribution
# = inpt of Comorbidity-QALY calculation
# To be extracted per age group from the Sciensano comorbidity dataset<jupyter_output><empty_output><jupyter_text>Initialize the QALY_model object<jupyter_code>QALY_model = QALY_model(comorbidity_distribution)
QALY_model.QoL_Van_Wilder
# Number of 'chronic diseases' of the Belgian population per age group
# Must be matched with the comorbidity definition of the general population
QALY_model.QoL_df
# By comparing the comorbidities of the studied population with the general population,
# quality-of-life scores can be computed
# = impact of comorbidities on life quality
QALY_model.SMR_pop_df
# Standardized mortality ratios for the Belgian population can be computed based on the
# "relative risk of death" (Defined in CCI units!) and the comorbidity prevalence
# If the number of comorbidities is higher than average, then SMR > 1 and vice versa.
# A person with the average number of comorbidities has an SMR = 1 which results in an average life expectancy
# = impact of comorbidities on life expectancy
QALY_model.SMR_df
# By comparing the comorbidities of the studied population with the general population,
# standardized mortality ratios can be computed
# = impact of comorbidities on life expectancy<jupyter_output><empty_output><jupyter_text>Compute and visualize survival functions and life expectancy<jupyter_code>fig,ax = plt.subplots(figsize=(12,4))
SMR_comorbid=2
ax.set_title('Survival function ($S_x$)')
ax.set_ylabel('cumulative chance of survival (-)')
ax.set_xlabel('age (y)')
ax.plot(QALY_model.survival_function(SMR=1),'.',color='black')
ax.plot(QALY_model.survival_function(SMR=SMR_comorbid),'--',color='black')
ax.legend(['baseline: SMR = 1', 'comorbid: SMR = '+ str(SMR_comorbid)], bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=13)
ax.grid(False)
plt.show()

fig,ax = plt.subplots(figsize=(12,4))
ax.set_title('Life expectancy ($LE_x$)')
ax.set_ylabel('cumulative chance of survival (-)')
ax.set_xlabel('age (y)')
ax.plot(QALY_model.life_expectancy(SMR=1),'.',color='black')
ax.plot(QALY_model.life_expectancy(SMR=SMR_comorbid),'--',color='black')
ax.legend(['baseline: SMR = 1', 'comorbid: SMR = '+ str(SMR_comorbid)], bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=13)
ax.grid(False)
plt.show()<jupyter_output><empty_output><jupyter_text>Compute and visualize QALEs en QALYs of the comorbid population<jupyter_code>fig,ax = plt.subplots(figsize=(12,4))
SMR_method = 'convergent'
# No comorbidity
ax.plot(QALY_model.compute_QALE_x(population='BE', SMR_method='convergent'), color='black')
# Comorbidity
ax.plot(QALY_model.compute_QALE_x(population='R', SMR_method='convergent'), '--',color='black')
# Comorbidity
ax.plot(QALY_model.compute_QALE_x(population='D', SMR_method='convergent'), '.',color='black')
ax.set_xlabel('age (y)')
ax.set_ylabel('QALE upon death (-)')
ax.legend(['Belgium', 'COVID-19 \nRecovered in hospital','COVID-19 \nDeceased in hospital'], bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=13)
ax.set_title('Quality-adjusted life expectancy ($QALE_x$)')
ax.grid(False)
plt.show()
fig,ax = plt.subplots(figsize=(12,4))
SMR_method = 'convergent'
r=0.03
# No comorbidity
ax.plot(QALY_model.compute_QALY_x(population='BE', SMR_method=SMR_method,r=r), color='black')
# Comorbidity
ax.plot(QALY_model.compute_QALY_x(population='R', SMR_method=SMR_method,r=r), '--',color='black')
# Comorbidity
ax.plot(QALY_model.compute_QALY_x(population='D', SMR_method=SMR_method,r=r), '.',color='black')
ax.set_xlabel('age (year)')
ax.set_ylabel('QALYs lost upon death (-)')
ax.legend(['Belgium', 'COVID-19 \nRecovered in hospital','COVID-19 \nDeceased in hospital'], bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=13)
ax.set_title('Quality-adjusted life years remaining at age $x$ ($QALY_x$)')
ax.grid(False)
plt.show()<jupyter_output><empty_output><jupyter_text>Map QALYs to COVID-19 SEIQRD age-bins<jupyter_code># Manual approach
QALY_model.bin_QALY_x(QALY_model.compute_QALY_x(population='BE', SMR_method='convergent',r=r))
# Automatically loop over populations and make a pd.Dataframe
binned_QALY_df = QALY_model.build_binned_QALY_df()<jupyter_output><empty_output><jupyter_text>Perform the acute QALY calculation and add to the xarray output of the model<jupyter_code>out = QALY_model.append_acute_QALY_losses(out, binned_QALY_df)<jupyter_output><empty_output><jupyter_text>Visualize number of QALYs lost<jupyter_code>fig,ax = plt.subplots(figsize=(12,4))
ax.plot(out['time'],out['QALYs_mild'].mean(dim='draws').cumsum().sum(dim='Nc'))
ax.fill_between(out['time'].values,out['QALYs_mild'].quantile(dim='draws',q=1-0.05/2).cumsum().sum(dim='Nc'),out['QALYs_mild'].quantile(dim='draws',q=0.05/2).cumsum().sum(dim='Nc'),alpha=0.2)
ax.plot(out['time'],out['QALYs_cohort'].mean(dim='draws').cumsum().sum(dim='Nc'))
ax.fill_between(out['time'].values,out['QALYs_cohort'].quantile(dim='draws',q=1-0.05/2).cumsum().sum(dim='Nc'),out['QALYs_cohort'].quantile(dim='draws',q=0.05/2).cumsum().sum(dim='Nc'),alpha=0.2)
ax.plot(out['time'],out['QALYs_ICU'].mean(dim='draws').cumsum().sum(dim='Nc'))
ax.fill_between(out['time'].values,out['QALYs_ICU'].quantile(dim='draws',q=1-0.05/2).cumsum().sum(dim='Nc'),out['QALYs_ICU'].quantile(dim='draws',q=0.05/2).cumsum().sum(dim='Nc'),alpha=0.2)
ax = _apply_tick_locator(ax)
ax.legend(['Mild symptoms', 'Hospitalization - Cohort', 'Hospitalization - IC'], bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=13)
ax.set_title('Acute QALY losses')
Palette=cm.get_cmap('tab10_r', levels).colors
age_group=['0-9','10-19','20-29','30-39','40-59','50-59','60-69','70-79','80+']


fig, axs = plt.subplots(3,figsize=(12,12),sharex=True)

axs[0].stackplot(out['time'].values,out['R_C'].mean(dim='draws').values+out['R_ICU'].mean(dim='draws').values, labels=age_group, colors=Palette, alpha=0.8)
axs[0].set_ylabel('Recovered in hospital')
axs[0].legend(fancybox=True, frameon=True, framealpha=1, fontsize=15,title='Age Group', loc="upper left", bbox_to_anchor=(1,1))  
axs[0] = _apply_tick_locator(axs[0])
axs[0].grid(False)

axs[1].stackplot(simtime,out['QALYs_treatment'].mean(dim="draws"),linewidth=3, labels=age_group, colors=Palette, alpha=0.8)
axs[1].set_ylabel('QALYs gained')
axs[1] = _apply_tick_locator(axs[1])
axs[1].grid(False)

y=[]
for i in range(levels):
    y.append(out['QALYs_treatment'].mean(dim="draws")[i].cumsum())
axs[2].stackplot(simtime,y,linewidth=3, labels=age_group, colors=Palette, alpha=0.8)
axs[2].set_ylabel('Cumulative QALYs gained')
axs[2] = _apply_tick_locator(axs[2])
axs[2].grid(False)

plt.savefig('../../results/QALYs/QALYs_gained.png',bbox_inches='tight',dpi=300)
#Plot

Palette=cm.get_cmap('tab10_r', levels).colors
age_group=['0-9','10-19','20-29','30-39','40-59','50-59','60-69','70-79','80+']

fig, axs = plt.subplots(2,figsize=(12,8),sharex=True)

axs[0].stackplot(out['time'].values,out['D'].mean(dim='draws').values, labels=age_group, colors=Palette, alpha=0.8)
axs[0].legend(fancybox=True, frameon=True, framealpha=1, fontsize=15,title='Age Group', loc="upper left", bbox_to_anchor=(1,1))  
axs[0] = _apply_tick_locator(axs[0])
axs[0].set_ylabel('Deceased in hospital')
axs[0].grid(False)

axs[1].stackplot(simtime,out['QALYs_death'].mean(dim="draws"),linewidth=3, labels=age_group, colors=Palette, alpha=0.8)
axs[1].set_ylabel('QALYs lost')
axs[1] = _apply_tick_locator(axs[1])
axs[1].grid(False)

states = ['QALYs_death','QALYs_treatment']
title = ['QALYs lost to COVID-19 death', 'QALYs gained by treating COVID-19']
ylabel = ['lost QALYs','gained QALYs']

plt.savefig('../../results/QALYs/QALYs_lost.png',bbox_inches='tight',dpi=300)<jupyter_output><empty_output><jupyter_text>The effect of hospital care reduction vs. COVID-19+ In the folowing example, the loss of life due to a constant average percentual reduction in regular hospital is calculated and compared with the overall loss of life to Covid-19.+ The provided function can be used to investigate the effect of time-dependent changes in hospital care as well as reductions in the treatment of different disease categories.<jupyter_code>from covid19model.models.QALY import lost_QALYs_hospital_care

# Define percentual reduction range
reduction_range=np.linspace(0.10,0.99,100)
# Calculate QALYs lost due to hospital care reduction
lost_hospital=lost_QALYs_hospital_care(reduction_range)
# Calculate how much QALYs are lost over a time period delta_T
delta_T = (simtime[-1]-simtime[0])/pd.Timedelta(days=1)
QALY_indirect = lost_hospital*delta_T
# Extract direct QALY loss from model
QALY_direct_mean = np.cumsum(out['QALYs_treatment'].mean(dim="draws").sum(dim="Nc").values)[-1]
QALY_direct_lower = np.cumsum(out['QALYs_treatment'].quantile(conf_int/2,dim="draws").sum(dim="Nc").values)[-1]
QALY_direct_upper = np.cumsum(out['QALYs_treatment'].quantile(1-conf_int/2,dim="draws").sum(dim="Nc").values)[-1]
# Calculate ratio
ratio_mean = QALY_direct_mean/QALY_indirect
ratio_upper = QALY_direct_lower/QALY_indirect
ratio_lower = QALY_direct_upper/QALY_indirect
#ratio_data = sum(deaths*lost_QALY_pp)/QALY_indirect
from covid19model.visualization.output import _apply_tick_locator

# Initialize plot
fig,ax = plt.subplots(figsize=(12,6))

# Add lines
ax.plot(reduction_range*100,np.ones(100),color='red',alpha=0.8,linewidth=2)
ax.plot(reduction_range*100,ratio_mean,'--',color='black')
ax.fill_between(reduction_range*100,ratio_lower,ratio_upper,alpha=0.2,color='black')

# limit the number of ticks on the axis
ax = _apply_tick_locator(ax)

# labels and title
ax.set_xlabel('Reduction of hospital care (%)')
ax.set_ylabel('QALY ratio (direct/indirect)')
#plt.title('Comparison of qualitative life years lost\n directly vs. indirectly to SARS-CoV-2')

# save figure
plt.savefig('../../results/QALYs/QALY_ratio.pdf',dpi=300)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1820721659248035, 'fraction_numerical': 0.033810019888246994, 'mean_word_length': 5.353489771359808, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 60, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '290529', 'license': 'MIT', 'n_tokens_mistral': 8065, 'n_tokens_neox': 7457, 'n_words': 1385, 'path': 'notebooks/scratch/twallema-QALYS.ipynb', 'repo_name': 'UGentBiomath/COVID19-Model', 'size': '21141'}