diff --git a/plot_data_ns.R b/plot_data_ns.R
index e7b2eab..da156c1 100644
--- a/plot_data_ns.R
+++ b/plot_data_ns.R
@@ -55,4 +55,4 @@ lf_fp_serum$obesity = as.factor(lf_fp_serum$obesity)
 head(lf_fp_sam$value[lf_fp_serum$mediator == "vitd"])
 ########################################################################
 
-rm(npa_lf, npa_wf, sam_lf, sam_wf, serum_lf, serum_wf, fp_adults_ics)
+rm(npa_lf, npa_wf, sam_lf, sam_wf, serum_lf, serum_wf)
diff --git a/plot_data_s.R b/plot_data_s.R
index 137eabd..3ef273c 100644
--- a/plot_data_s.R
+++ b/plot_data_s.R
@@ -55,4 +55,4 @@ lf_fp_serum$obesity = as.factor(lf_fp_serum$obesity)
 head(lf_fp_sam$value[lf_fp_serum$mediator == "vitd"])
 ########################################################################
 
-rm(npa_lf, npa_wf, sam_lf, sam_wf, serum_lf, serum_wf, fp_adults_ics)
+rm(npa_lf, npa_wf, sam_lf, sam_wf, serum_lf, serum_wf)
diff --git a/read_data.R b/read_data.R
index c6be2f3..d30bf74 100755
--- a/read_data.R
+++ b/read_data.R
@@ -23,6 +23,12 @@ ifelse(!dir.exists(outdir_stats_na), dir.create(outdir_stats_na), FALSE)
 outdir_stats_ns = paste0(maindir, "output/stats/non_severe/")
 ifelse(!dir.exists(outdir_stats_ns), dir.create(outdir_stats_ns), FALSE)
 
+outdir_stats_s = paste0(maindir, "output/stats/severe/")
+ifelse(!dir.exists(outdir_stats_s), dir.create(outdir_stats_s), FALSE)
+
+outdir_stats_a = paste0(maindir, "output/stats/asthma/")
+ifelse(!dir.exists(outdir_stats_a), dir.create(outdir_stats_a), FALSE)
+
 outdir_plots = paste0(maindir, "output/plots/")
 ifelse(!dir.exists(outdir_plots), dir.create(outdir_plots), FALSE)
 ########################################################################
diff --git a/stats_check.R b/stats_check.R
index 7d01038..6b5a967 100755
--- a/stats_check.R
+++ b/stats_check.R
@@ -15,7 +15,15 @@ wf_serum = wf_fp_serum
 # Numbers
 table(fp_adults_ics$asthma, fp_adults_ics$obesity)
 
-table(fp_adults_ics$obesity)
+table(fp_adults_ics$obesity); table(fp_adults_ics$ia_cxr)
+
+table(fp_adults_ics$obesity, fp_adults_ics$ia_cxr)
+fisher.test(table(fp_adults_ics$obesity, fp_adults_ics$ia_cxr))
+
+table(fp_adults_ics$obesity, fp_adults_ics$com_noasthma)
+fisher.test(table(fp_adults_ics$obesity, fp_adults_ics$com_noasthma))
+
+chisq.test(table(fp_adults_ics$obesity, fp_adults_ics$com_noasthma))
 
 
 #=========