Validate results wolf again

data/QALY_model/interim/long_COVID/average_QALY_losses_per_age_group_SMR100.csv

@@ -0,0 +1,41 @@
+hospitalisation,age_group,mean,sd,lower,upper
+Non-hospitalised,"[0, 12)",,,,
+Non-hospitalised,"[12, 18)",,,,
+Non-hospitalised,"[18, 25)",,,,
+Non-hospitalised,"[25, 35)",,,,
+Non-hospitalised,"[35, 45)",,,,
+Non-hospitalised,"[45, 55)",,,,
+Non-hospitalised,"[55, 65)",,,,
+Non-hospitalised,"[65, 75)",,,,
+Non-hospitalised,"[75, 85)",,,,
+Non-hospitalised,"[85, 120)",,,,
+Hospitalised (no IC),"[0, 12)",,,,
+Hospitalised (no IC),"[12, 18)",,,,
+Hospitalised (no IC),"[18, 25)",,,,
+Hospitalised (no IC),"[25, 35)",,,,
+Hospitalised (no IC),"[35, 45)",,,,
+Hospitalised (no IC),"[45, 55)",,,,
+Hospitalised (no IC),"[55, 65)",,,,
+Hospitalised (no IC),"[65, 75)",,,,
+Hospitalised (no IC),"[75, 85)",,,,
+Hospitalised (no IC),"[85, 120)",,,,
+Hospitalised (IC),"[0, 12)",,,,
+Hospitalised (IC),"[12, 18)",,,,
+Hospitalised (IC),"[18, 25)",,,,
+Hospitalised (IC),"[25, 35)",,,,
+Hospitalised (IC),"[35, 45)",,,,
+Hospitalised (IC),"[45, 55)",,,,
+Hospitalised (IC),"[55, 65)",,,,
+Hospitalised (IC),"[65, 75)",,,,
+Hospitalised (IC),"[75, 85)",,,,
+Hospitalised (IC),"[85, 120)",,,,
+Non-hospitalised (no AD),"[0, 12)",,,,
+Non-hospitalised (no AD),"[12, 18)",,,,
+Non-hospitalised (no AD),"[18, 25)",,,,
+Non-hospitalised (no AD),"[25, 35)",,,,
+Non-hospitalised (no AD),"[35, 45)",,,,
+Non-hospitalised (no AD),"[45, 55)",,,,
+Non-hospitalised (no AD),"[55, 65)",,,,
+Non-hospitalised (no AD),"[65, 75)",,,,
+Non-hospitalised (no AD),"[75, 85)",,,,
+Non-hospitalised (no AD),"[85, 120)",,,,

data/QALY_model/interim/long_COVID/average_QALY_losses_per_age_group_SMR150.csv

@@ -0,0 +1,41 @@
+hospitalisation,age_group,mean,sd,lower,upper
+Non-hospitalised,"[0, 12)",0.668462419110261,0.2582313828457331,0.2036324547925269,1.2266842634252282
+Non-hospitalised,"[12, 18)",0.5798944237428355,0.22187795956279144,0.18094233731503914,1.0596337319588378
+Non-hospitalised,"[18, 25)",0.5190581364867876,0.1968431158317833,0.16537131828772747,0.9447050484446313
+Non-hospitalised,"[25, 35)",0.44763489426351,0.16736730072686784,0.1471094089416646,0.8092417221722648
+Non-hospitalised,"[35, 45)",0.37120795187222455,0.13572307883109985,0.12759033988719307,0.6642768756212402
+Non-hospitalised,"[45, 55)",0.3002411506332599,0.10627491577506192,0.10947793400313126,0.5297938114190412
+Non-hospitalised,"[55, 65)",0.23622082675172743,0.07971458883081076,0.0929959256939798,0.4084045703816341
+Non-hospitalised,"[65, 75)",0.17618743525018343,0.05486694992256693,0.07783840424384207,0.29459863667090785
+Non-hospitalised,"[75, 85)",0.12206526711280687,0.03255206915706073,0.06381993641723875,0.1920038620032939
+Non-hospitalised,"[85, 120)",0.08103127670002823,0.015751336483312568,0.05284996746867399,0.1150524706145426
+Hospitalised (no IC),"[0, 12)",1.8980894392621346,0.7267827620017653,0.6207356011360611,3.4240096409038068
+Hospitalised (no IC),"[12, 18)",1.6478459932072107,0.6300813814747115,0.546874783644301,2.9726944049100523
+Hospitalised (no IC),"[18, 25)",1.472820461853593,0.5617587044599202,0.49489406166793193,2.654800869057853
+Hospitalised (no IC),"[25, 35)",1.2634184870428975,0.47914425110921366,0.4321870750253125,2.272532748393311
+Hospitalised (no IC),"[35, 45)",1.0346805795675365,0.3878601694223998,0.3623500214659087,1.8514667391139488
+Hospitalised (no IC),"[45, 55)",0.8193694351258769,0.30129304663868833,0.2966104797693478,1.4529560947554594
+Hospitalised (no IC),"[55, 65)",0.6252626513069367,0.22328391548470894,0.23728421113389478,1.0937858217247964
+Hospitalised (no IC),"[65, 75)",0.44562604115188553,0.15161325700772058,0.18080980173221772,0.7631334746039069
+Hospitalised (no IC),"[75, 85)",0.28685503228783793,0.08900307321463005,0.13044941001272536,0.4732435233307565
+Hospitalised (no IC),"[85, 120)",0.16891493919738063,0.04328019853546922,0.09158887797754195,0.2612073485364238
+Hospitalised (IC),"[0, 12)",2.9237322020738303,0.9793705841671974,1.2107236313271368,4.999388921316399
+Hospitalised (IC),"[12, 18)",2.531638312858393,0.8497821429058505,1.0519529833728634,4.331883645375432
+Hospitalised (IC),"[18, 25)",2.257562381291156,0.757962284073033,0.9411960051611258,3.8628058096578908
+Hospitalised (IC),"[25, 35)",1.9298543316512762,0.6466155391861154,0.8090419079356348,3.2987897263962442
+Hospitalised (IC),"[35, 45)",1.5721185906351458,0.5232154765316691,0.6651014247066298,2.679412531598891
+Hospitalised (IC),"[45, 55)",1.2355267693612049,0.40597563326945363,0.5298650819811873,2.094388339582217
+Hospitalised (IC),"[55, 65)",0.9320835477121102,0.3003429071567427,0.40793600939386565,1.56709365027751
+Hospitalised (IC),"[65, 75)",0.6511568673099988,0.2034594772850809,0.2948885959200651,1.0807438310262218
+Hospitalised (IC),"[75, 85)",0.4027184979804148,0.1189921363368849,0.1936167990895866,0.6530541829319151
+Hospitalised (IC),"[85, 120)",0.21808592502322638,0.057238589088626377,0.11680411982667387,0.33711143404710353
+Non-hospitalised (no AD),"[0, 12)",0.06014137527980275,0.002961992425985527,0.05453426633193338,0.06582908835826626
+Non-hospitalised (no AD),"[12, 18)",0.06011715444004468,0.002961992425985527,0.054510045492175314,0.0658048675185082
+Non-hospitalised (no AD),"[18, 25)",0.06009237832769276,0.0029619924259855273,0.05448526937982339,0.06578009140615626
+Non-hospitalised (no AD),"[25, 35)",0.06005474395604065,0.0029619924259855278,0.054447635008171284,0.06574245703450415
+Non-hospitalised (no AD),"[35, 45)",0.06000375519059451,0.002961992425985528,0.054396646242725145,0.06569146826905803
+Non-hospitalised (no AD),"[45, 55)",0.059945480395673406,0.002961992425985528,0.05433837144780404,0.06563319347413692
+Non-hospitalised (no AD),"[55, 65)",0.05988345587040076,0.002961992425985525,0.054276346922531395,0.06557116894886428
+Non-hospitalised (no AD),"[65, 75)",0.05981520268326211,0.0029619924259855247,0.054208093735392746,0.06550291576172562
+Non-hospitalised (no AD),"[75, 85)",0.05974245826449364,0.0029619924191818985,0.05413534932950368,0.06543017132989262
+Non-hospitalised (no AD),"[85, 120)",0.05963967771162379,0.0029605970042578026,0.05403521032407272,0.06532471125642572

notebooks/analysis/woldmuyn_long_COVID.py

@@ -38,16 +38,14 @@
     ## Define simulation settings ##
     ################################
 
-    # Number of simulations
-    #N=50
     # Number of neg. binomial draws/ simulation
     K=10
     # Number of cpu's
     processes=int(mp.cpu_count()/2)
     # Number of age groups
     age_stratification_size=10
     # End of simulation
-    end_sim='2021-12-31'
+    end_sim='2021-07-01'
     # Confidence level used to visualise model fit
     conf_int=0.05
 
@@ -116,14 +114,12 @@
       ax.grid(False)
       ax.tick_params(axis='both', which='major', labelsize=10)
       ax.tick_params(axis='x', which='major', rotation=30)
-      ax.set_ylim([0,7000])
+      ax.set_ylim([0,6600])
     axes[0].legend(loc=2, framealpha=1, fontsize=10)
     axes[0].set_ylabel('QALYs lost per 100K inhab.', size=10)
 
     plt.tight_layout()
-    plt.show()
-    fig.savefig('QALY_losses_per_age_group.pdf')
-    plt.close()
+    fig.savefig(result_folder+'QALY_losses_per_age_group.pdf')
 
     # Wolf's code starts here
     label_font = font_manager.FontProperties(family='CMU Sans Serif',
@@ -206,6 +202,6 @@
         lower = np.quantile(total,0.025)
         upper = np.quantile(total,0.975)
 
-        QALY_table[total_label]['Total'] = f'{mean:.0f}\n({lower:.0f};{upper:.0f})'
+        QALY_table[total_label]['Total'] = f'{mean:.0f}\n({lower:.0f}; {upper:.0f})'
 
     QALY_table.to_csv(os.path.join(result_folder,f'Long_COVID_summary_SMR{SMR*100:.0f}.csv'))

notebooks/preprocessing/woldmuyn_long_covid_average_QALY_loss.py

@@ -97,14 +97,6 @@
 # results #
 # ------- #
 
-# latex font
-label_font = font_manager.FontProperties(family='CMU Sans Serif',
-                                   style='normal', 
-                                   size=10)
-legend_font = font_manager.FontProperties(family='CMU Sans Serif',
-                                   style='normal', 
-                                   size=8)
-
 # result folders
 data_result_folder = '../../data/QALY_model/interim/long_COVID/'
 fig_result_folder = '../../results/QALY_model/direct_QALYs/prepocessing/'
@@ -247,7 +239,7 @@ def log_probability(theta, tau, x, y, bounds):
     )
     samplers.update({hospitalisation:sampler})
     # run sampler
-    sampler.run_mcmc(pos, 25000, progress=True)
+    sampler.run_mcmc(pos, 50000, progress=True)
     # extract chains
     flat_samples = sampler.get_chain(discard=5000, thin=100, flat=True)
     # print results
@@ -256,6 +248,8 @@ def log_probability(theta, tau, x, y, bounds):
     p_AD_summary['lower'][hospitalisation] = np.quantile(flat_samples,0.025,axis=0)
     p_AD_summary['upper'][hospitalisation] = np.quantile(flat_samples,0.975,axis=0)
 
+print(p_AD_summary)
+
 # visualise MCMC results    
 fig,axes = plt.subplots(nrows=1,ncols=3,figsize=(8.3,0.25*11.7),sharey=True)
 
@@ -305,8 +299,8 @@ def log_probability(theta, tau, x, y, bounds):
 ax.set_xlabel('Age (years)',size=10)
 ax.set_ylabel('QoL (-)',size=10)
 ax.set_ylim([0.5, 0.9])
-ax.grid(False)
 ax.tick_params(axis='both', which='major', labelsize=10)
+ax.grid(False)
 plt.tight_layout()
 fig.savefig(os.path.join(abs_dir,fig_result_folder,'QoL_Belgium_fit.pdf'))
 
@@ -323,7 +317,7 @@ def QALY_loss_func(t,tau,p_AD,age,QoL_after):
     beta = QoL_Belgium_func(age+t/12)-QoL_after
     return prevalence_func(t,tau,p_AD) * max(0,beta)
 
-draws = 1000
+draws = 500
 
 # Pre-allocate new multi index series with index=hospitalisation,age,draw
 multi_index = pd.MultiIndex.from_product([hospitalisation_groups+['Non-hospitalised (no AD)'],np.arange(draws),LE_table.index.values],names=['hospitalisation','draw','age'])
@@ -351,8 +345,8 @@ def QALY_loss_func(t,tau,p_AD,age,QoL_after):
     QoL_after = QoL_Belgium_func(age)-beta
     # integrate QALY_loss_func from 0 to LE (OPM: kan ook discreet) 
     QALY_loss = quad(QALY_loss_func,0,LE,args=(tau,p_AD,age,QoL_after))[0]/12 
-    average_QALY_losses.iloc[idx] = QALY_loss
-    
+    average_QALY_losses_per_age.iloc[idx] = QALY_loss
+
 print('\n(4.2) Bin average QALY loss per age to age groups\n')
 
 # bin data
@@ -394,13 +388,12 @@ def get_sd(x):
 
 print('\n(6) Visualise results\n')
 
-# Visualise results (TWA)
+# QALY loss per age group
 fig,axes = plt.subplots(nrows=1,ncols=3,figsize=(8.3,0.25*11.7),sharey=True)
-
 for ax,hospitalisation in zip(axes,hospitalisation_groups):
-    mean = average_QALY_losses_summary.loc[hospitalisation]['mean']
-    lower = average_QALY_losses_summary.loc[hospitalisation]['lower']
-    upper = average_QALY_losses_summary.loc[hospitalisation]['upper']
+    mean = average_QALY_losses_per_age_summary.loc[hospitalisation]['mean']
+    lower = average_QALY_losses_per_age_summary.loc[hospitalisation]['lower']
+    upper = average_QALY_losses_per_age_summary.loc[hospitalisation]['upper']
     ax.plot(LE_table.index.values,mean,color='black',linewidth=1.5, linestyle='--')
     ax.plot(LE_table.index.values,lower,color='black',linewidth=1, linestyle='-')
     ax.plot(LE_table.index.values,upper,color='black',linewidth=1, linestyle='-')
@@ -413,42 +406,24 @@ def get_sd(x):
 plt.tight_layout()
 fig.savefig(os.path.join(abs_dir,fig_result_folder,'average_QALY_losses_per_age.pdf'))
 
-# Visualise results (WD)
-fig,axs = plt.subplots(1,3,figsize=(6,2.5),sharey=True,sharex=True)
-for ax,hospitalisation in zip(axs,hospitalisation_groups):
-    mean = average_QALY_losses_per_age_summary.loc[hospitalisation]['mean']
-    lower = average_QALY_losses_per_age_summary.loc[hospitalisation]['lower']
-    upper = average_QALY_losses_per_age_summary.loc[hospitalisation]['upper']
-    ax.plot(LE_table.index.values,mean,color=palette_colors[color_dict[hospitalisation]],linestyle='--',label=f'{hospitalisation}',linewidth=1)
-    ax.fill_between(LE_table.index.values,lower,upper,alpha=0.20, color=palette_colors[color_dict[hospitalisation]])
-
-    ax.grid(False)
-    ax.set_xlabel('Age when infected (years)',font=label_font)
-    ax.tick_params(axis='both', which='major', labelsize=8)
-    ax.set_title(hospitalisation,font=label_font)
-
-axs[0].set_ylabel('Average QALY loss',font=label_font)
-fig.tight_layout()
-fig.savefig(os.path.join(abs_dir,fig_result_folder,f'average_QALY_losses_per_age_SMR{SMR*100:.0f}.png'),dpi=600,bbox_inches='tight')
-
 # QALY losses due COVID death
 fig,ax = plt.subplots(figsize=(5,3))
-ax.plot(Life_table.compute_QALY_D_x(r=0,SMR=SMR),color=palette_colors['black'],label=r'$r=0\%$',linewidth=2)
-ax.plot(Life_table.compute_QALY_D_x(r=0.03,SMR=SMR),color=palette_colors['black'],linestyle=':',label=r'$r=3\%$',linewidth=2)
+ax.plot(Life_table.compute_QALY_D_x(r=0,SMR=SMR),color='black',label=r'$r=0\%$',linewidth=2)
+ax.plot(Life_table.compute_QALY_D_x(r=0.03,SMR=SMR),color='black',linestyle=':',label=r'$r=3\%$',linewidth=2)
 ax.grid(False)
-ax.set_xlabel('Age (years)',font=label_font)
-ax.set_ylabel(r'$QALY_D$',font=label_font)
+ax.set_xlabel('Age (years)')
+ax.set_ylabel(r'$QALY_D$')
 ax.tick_params(axis='both', which='major', labelsize=8)
-ax.legend(prop=legend_font)
+ax.legend()
 fig.tight_layout()
-fig.savefig(os.path.join(abs_dir,fig_result_folder,f'QALY_D_SMR{SMR*100:.0f}.png'),dpi=600,bbox_inches='tight')
+fig.savefig(os.path.join(abs_dir,fig_result_folder,f'QALY_D_SMR{SMR*100:.0f}.pdf'),dpi=600,bbox_inches='tight')
 
 # Life expectancy
 fig,ax = plt.subplots(figsize=(3,3))
 ax.plot(LE_table,'black',linewidth=1.5)
 ax.grid(False)
-ax.set_ylabel('Life expectancy (years)',font=label_font)
-ax.set_xlabel('Age (years)',font=label_font)
+ax.set_ylabel('Life expectancy (years)')
+ax.set_xlabel('Age (years)')
 ax.tick_params(axis='both', which='major', labelsize=8)
 fig.tight_layout()
-fig.savefig(os.path.join(abs_dir,fig_result_folder,f'LE_SMR{SMR*100:.0f}.png'),dpi=600,bbox_inches='tight')
+fig.savefig(os.path.join(abs_dir,fig_result_folder,f'LE_SMR{SMR*100:.0f}.pdf'),dpi=600,bbox_inches='tight')