diff --git a/flu_stats_unpaired_npa.R b/flu_stats_unpaired_npa.R
index e8282ce..bcc0c81 100755
--- a/flu_stats_unpaired_npa.R
+++ b/flu_stats_unpaired_npa.R
@@ -8,7 +8,6 @@ getwd()
 # data: Flu positive adult patients
 # group: obesity
 ############################################################
-
 my_sample_type = "npa"
 
 #=============
@@ -19,12 +18,20 @@ source("data_extraction_formatting.R")
 # check: adult variable and age variable discrepancy!
 metadata_all$mosaic[metadata_all$adult==1 & metadata_all$age<=18]
 
-#=========================================================
+#=============
+# Output
+#=============
+outfile_name = paste0("flu_stats_time_unpaired_", my_sample_type, ".csv")
+flu_stats_time_unpaired = paste0(outdir_stats, outfile_name)
+
+#===============================
 # data assignment for stats
+#================================
 wf = npa_wf[npa_wf$flustat == 1,]
 lf = npa_lf[npa_lf$flustat == 1,]
 lf$timepoint = paste0("t", lf$timepoint)
-#=========================================================
+
+########################################################################
 # clear variables
 rm(sam_lf, sam_wf
 , serum_lf, serum_wf) 
@@ -33,13 +40,6 @@ rm(colnames_sam_df, expected_rows_sam_lf
 
 rm(pivot_cols)
 
-#=============
-# Output: unpaired analysis of time for npa
-#=============
-outfile_name = paste0("flu_stats_time_unpaired_", my_sample_type, ".csv")
-
-flu_stats_time_unpaired = paste0(outdir_stats, outfile_name)
-########################################################################
 # sanity checks
 table(lf$timepoint)
 
@@ -52,24 +52,20 @@ if (table(lf$flustat) == table(npa_lf$flustat)[[2]]){
 ########################################################################
 # Unpaired stats at each timepoint b/w groups: wilcoxon UNpaired analysis 
 # with correction
-#######################################################################
+########################################################################
 # with adjustment: fdr and BH are identical
 my_adjust_method = "BH"
 
 #==============
 # unpaired: t1
-#==============
 lf_t1 = lf[lf$timepoint == "t1",]
 sum(is.na(lf_t1$value))
 
-foo = lf_t1[which(is.na(lf_t1$value)),]
-ci = which(is.na(lf_t1$value))
-
+#foo = lf_t1[which(is.na(lf_t1$value)),]
+#ci = which(is.na(lf_t1$value))
 #lf_t1_comp = lf_t1[-ci,]
+
 lf_t1_comp = lf_t1[-which(is.na(lf_t1$value)),]
-#--------------------
-# unpaired stats: t1
-#--------------------
 stats_un_t1 = compare_means(value~obesity
                             , group.by = "mediator"
                             #, data = lf_t1 
@@ -89,9 +85,9 @@ n_t1 = data.frame(table(lf_t1_comp$mediator))
 colnames(n_t1) = c("mediator", "n_obs")
 n_t1$mediator = as.character(n_t1$mediator)
 
-#=========================
-# Merge1: merge stats + n_obs df
-#=========================
+#==================================
+# Merge: merge stats + n_obs df
+#===================================
 merging_cols = intersect(names(stats_un_t1), names(n_t1)); merging_cols
 if (all(n_t1$mediator%in%stats_un_t1$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -113,39 +109,6 @@ if (all(n_t1$mediator%in%stats_un_t1$mediator)) {
 # add bonferroni adjustment as well
 stats_un_t1$p_adj_bonferroni = p.adjust(stats_un_t1$p, method = "bonferroni")
 
-#--------------------
-# summary stats: t1
-#--------------------
-gp_stats_t1 = groupedstats::grouped_summary(
-  data = lf_t1_comp,
-  grouping.vars = c(mediator, obesity),
-  measures = value,
-  measures.type = "numeric")
-
-#====================================
-# Merge 2: Merge1 + summary_stats
-#====================================
-merging_cols2 = intersect(names(stats_un_t1), names(gp_stats_t1)); merging_cols2
-if (all(gp_stats_t1$mediator%in%stats_un_t1$mediator)) {
-  cat("PASS: merging summary stats with unpaired analysis on column/s:", merging_cols2)
-  stats_un_t1 = merge(stats_un_t1, gp_stats_t1, by = merging_cols2, all = T)
-  cat("\nsuccessfull merge:"
-      , "\nnrow:", nrow(stats_un_t1)
-      , "\nncol:", ncol(stats_un_t1))
-}else{
-  nf = gp_stats_t1$mediator[!gp_stats_t1$mediator%in%stats_un_t1$mediator]
-  stats_un_t1 = merge(stats_un_t1, gp_stats_t1, by = merging_cols2, all.y = T)
-  cat("\nMerged with caution:"
-      , "\nnrows mismatch:", nf
-      , "not found in stats possibly due to all obs being LLODs"
-      , "\nintroduced NAs for:", nf
-      , "\nnrow:", nrow(stats_un_t1)
-      , "\nncol:", ncol(stats_un_t1))
-}
-
-
-
-
 rm(n_t1)
 rm(lf_t1_comp)
 
@@ -161,8 +124,9 @@ stats_un_t2 = compare_means(value~obesity
                             , data = lf_t2_comp
                             , paired = FALSE
                             , p.adjust.method = my_adjust_method)
-stats_un_t2$timepoint = "t2"
 
+# add timepoint and convert to df
+stats_un_t2$timepoint = "t2"
 stats_un_t2 = as.data.frame(stats_un_t2)
 class(stats_un_t2)
 
@@ -171,7 +135,9 @@ n_t2 = data.frame(table(lf_t2_comp$mediator))
 colnames(n_t2) = c("mediator", "n_obs")
 n_t2$mediator = as.character(n_t2$mediator)
 
-# merge stats + n_obs df
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
 merging_cols = intersect(names(stats_un_t2), names(n_t2)); merging_cols
 if (all(n_t2$mediator%in%stats_un_t2$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -209,8 +175,8 @@ stats_un_t3 = compare_means(value~obesity
                             , paired = FALSE
                             , p.adjust.method = my_adjust_method)
 
+# add timepoint and convert to df
 stats_un_t3$timepoint = "t3"
-
 stats_un_t3 = as.data.frame(stats_un_t3)
 class(stats_un_t3)
 
@@ -219,7 +185,9 @@ n_t3 = data.frame(table(lf_t3_comp$mediator))
 colnames(n_t3) = c("mediator", "n_obs")
 n_t3$mediator = as.character(n_t3$mediator)
 
-# merge stats + n_obs df
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
 merging_cols = intersect(names(stats_un_t3), names(n_t3)); merging_cols
 if (all(n_t3$mediator%in%stats_un_t3$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -244,7 +212,7 @@ stats_un_t3$p_adj_bonferroni = p.adjust(stats_un_t3$p, method = "bonferroni")
 rm(n_t3)
 rm(lf_t3_comp)
 
-###########################################################################################
+########################################################################
 #==============
 # Rbind these dfs
 #==============
@@ -364,7 +332,7 @@ colnames(combined_unpaired_stats_f) = c("mediator"
                                         , "p_bon_signif")
 
 colnames(combined_unpaired_stats_f)
-
+########################################################################
 #******************
 # write output file
 #******************
diff --git a/flu_stats_unpaired_sam.R b/flu_stats_unpaired_sam.R
index 570c575..67be452 100755
--- a/flu_stats_unpaired_sam.R
+++ b/flu_stats_unpaired_sam.R
@@ -3,45 +3,54 @@ getwd()
 setwd("~/git/mosaic_2020/")
 getwd()
 ############################################################
-# TASK: unpaired (time) analysis of mediators: SAM
+# TASK: unpaired (time) analysis of mediators:
+# sample type: SAM
+# data: Flu positive adult patients
+# group: obesity
 ############################################################
+my_sample_type = "sam"
+
 #=============
 # Input
 #=============
 source("data_extraction_formatting.R")
 
-table(metadata_all$flustat[metadata_all$adult == 1])
+# check: adult variable and age variable discrepancy!
+metadata_all$mosaic[metadata_all$adult==1 & metadata_all$age<=18]
 
-# clear variables
-rm(npa_adults_lf, npa_df_adults_clean
-, serum_adults_lf, serum_df_adults_clean) 
-rm(colnames_npa_df, expected_rows_npa_lf
-, colnames_serum_df, expected_rows_serum_lf) 
-
-rm(pivot_cols)
-
-my_sample_type = "sam"
 #=============
 # Output: unpaired analysis of time for sam
 #=============
 outfile_name = paste0("flu_stats_time_unpaired_", my_sample_type, ".csv")
 flu_stats_time_unpaired = paste0(outdir_stats, outfile_name)
 
-#%%========================================================
+#===============================
 # data assignment for stats
-wf = sam_df_adults_clean[sam_df_adults_clean$flustat == 1,]
-lf = sam_adults_lf[sam_adults_lf$flustat == 1,]
- #%%========================================================
+#================================
+wf = sam_wf[sam_wf$flustat == 1,]
+lf = sam_lf[sam_lf$flustat == 1,]
+lf$timepoint = paste0("t", lf$timepoint)
+########################################################################
+# clear variables
+rm(npa_lf, npa_wf
+, serum_lf, serum_wf) 
+rm(colnames_npa_df, expected_rows_npa_lf
+, colnames_serum_df, expected_rows_serum_lf) 
+
+rm(pivot_cols)
+
+# sanity checks
 table(lf$timepoint)
 length(unique(lf$mosaic))
 lf$timepoint = paste0("t", lf$timepoint)
 
-if (table(lf$flustat) == table(sam_adults_lf$flustat)[[2]]){
+if (table(lf$flustat) == table(sam_lf$flustat)[[2]]){
   cat("Analysing Flu positive patients for:", my_sample_type)
 }else{
   cat("FAIL: problem with subsetting data for:", my_sample_type)
   quit()
 }
+
 ########################################################################
 # Unpaired stats at each timepoint b/w groups: wilcoxon UNpaired analysis 
 # with correction
@@ -55,10 +64,10 @@ my_adjust_method = "BH"
 lf_t1 = lf[lf$timepoint == "t1",]
 sum(is.na(lf_t1$value))
 
-foo = lf_t1[which(is.na(lf_t1$value)),]
-ci = which(is.na(lf_t1$value))
-
+#foo = lf_t1[which(is.na(lf_t1$value)),]
+#ci = which(is.na(lf_t1$value))
 #lf_t1_comp = lf_t1[-ci,]
+
 lf_t1_comp = lf_t1[-which(is.na(lf_t1$value)),]
 stats_un_t1 = compare_means(value~obesity
                             , group.by = "mediator"
@@ -69,8 +78,8 @@ stats_un_t1 = compare_means(value~obesity
 
 foo$mosaic[!unique(foo$mosaic)%in%unique(lf_t1_comp$mosaic)]
 
+# add timepoint and convert to df
 stats_un_t1$timepoint = "t1"
-
 stats_un_t1 = as.data.frame(stats_un_t1)
 class(stats_un_t1)
 
@@ -79,7 +88,9 @@ n_t1 = data.frame(table(lf_t1_comp$mediator))
 colnames(n_t1) = c("mediator", "n_obs")
 n_t1$mediator = as.character(n_t1$mediator)
 
-# merge stats + n_obs df
+#==================================
+# Merge: merge stats + n_obs df
+#===================================
 merging_cols = intersect(names(stats_un_t1), names(n_t1)); merging_cols
 if (all(n_t1$mediator%in%stats_un_t1$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -116,8 +127,8 @@ stats_un_t2 = compare_means(value~obesity
                             , data = lf_t2_comp
                             , paired = FALSE
                             , p.adjust.method = my_adjust_method)
+# add timepoint and convert to df                           
 stats_un_t2$timepoint = "t2"
-
 stats_un_t2 = as.data.frame(stats_un_t2)
 class(stats_un_t2)
 
@@ -126,7 +137,9 @@ n_t2 = data.frame(table(lf_t2_comp$mediator))
 colnames(n_t2) = c("mediator", "n_obs")
 n_t2$mediator = as.character(n_t2$mediator)
 
-# merge stats + n_obs df
+#==================================
+# Merge: merge stats + n_obs df
+#===================================
 merging_cols = intersect(names(stats_un_t2), names(n_t2)); merging_cols
 if (all(n_t2$mediator%in%stats_un_t2$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -163,9 +176,8 @@ stats_un_t3 = compare_means(value~obesity
                             , data = lf_t3_comp
                             , paired = FALSE
                             , p.adjust.method = my_adjust_method)
-
+# add timepoint and convert to df                            
 stats_un_t3$timepoint = "t3"
-
 stats_un_t3 = as.data.frame(stats_un_t3)
 class(stats_un_t3)
 
@@ -174,7 +186,9 @@ n_t3 = data.frame(table(lf_t3_comp$mediator))
 colnames(n_t3) = c("mediator", "n_obs")
 n_t3$mediator = as.character(n_t3$mediator)
 
-# merge stats + n_obs df
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
 merging_cols = intersect(names(stats_un_t3), names(n_t3)); merging_cols
 if (all(n_t3$mediator%in%stats_un_t3$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -202,10 +216,10 @@ stats_un_t3$p_adj_bonferroni = p.adjust(stats_un_t3$p, method = "bonferroni")
 
 rm(n_t3)
 rm(lf_t3_comp)
-
-#==============
+########################################################################
+#=================
 # Rbind these dfs
-#==============
+#=================
 str(stats_un_t1);str(stats_un_t2); str(stats_un_t3)
 
 n_dfs = 3
@@ -322,6 +336,7 @@ colnames(combined_unpaired_stats_f) = c("mediator"
                                         , "p_bon_signif")
 
 colnames(combined_unpaired_stats_f)
+#######################################################################
 #******************
 # write output file
 #******************
diff --git a/flu_stats_unpaired_serum.R b/flu_stats_unpaired_serum.R
index 0400963..42c9687 100755
--- a/flu_stats_unpaired_serum.R
+++ b/flu_stats_unpaired_serum.R
@@ -5,37 +5,43 @@ getwd()
 ############################################################
 # TASK: unpaired (time) analysis of mediators: serum
 ############################################################
+my_sample_type = "serum"
+
 #=============
 # Input
 #=============
 source("data_extraction_formatting.R")
 
-table(metadata_all$flustat[metadata_all$adult == 1])
+# check: adult variable and age variable discrepancy!
+metadata_all$mosaic[metadata_all$adult==1 & metadata_all$age<=18]
 
-# clear variables
-rm(sam_adults_lf, sam_df_adults_clean
-, npa_adults_lf, npa_df_adults_clean) 
-rm(colnames_sam_df, expected_rows_sam_lf
-, colnames_npa_df, expected_rows_npa_lf)
-
-rm(pivot_cols)
-
-my_sample_type = "serum"
 #=============
 # Output: unpaired analysis of time for serum
 #=============
 outfile_name = paste0("flu_stats_time_unpaired_", my_sample_type, ".csv")
 flu_stats_time_unpaired = paste0(outdir_stats, outfile_name)
 
-#%%========================================================
+#===============================
 # data assignment for stats
-wf = serum_df_adults_clean[serum_df_adults_clean$flustat == 1,]
-lf = serum_adults_lf[serum_adults_lf$flustat == 1,]
-#%%========================================================
+#================================
+wf = serum_wf[serum_wf$flustat == 1,]
+lf = serum_lf[serum_lf$flustat == 1,]
+lf$timepoint = paste0("t", lf$timepoint)
+
+########################################################################
+# clear variables
+rm(sam_lf, sam_wf
+, npa_lf, npa_wf) 
+rm(colnames_sam_df, expected_rows_sam_lf
+, colnames_npa_df, expected_rows_npa_lf)
+
+rm(pivot_cols)
+
+# sanity checks
 table(lf$timepoint)
 lf$timepoint = paste0("t", lf$timepoint)
 
-if (table(lf$flustat) == table(serum_adults_lf$flustat)[[2]]){
+if (table(lf$flustat) == table(serum_lf$flustat)[[2]]){
   cat("Analysing Flu positive patients for:", my_sample_type)
 }else{
   cat("FAIL: problem with subsetting data for:", my_sample_type)
@@ -43,7 +49,8 @@ if (table(lf$flustat) == table(serum_adults_lf$flustat)[[2]]){
 }
 
 ########################################################################
-# Unpaired stats at each timepoint b/w groups: wilcoxon UNpaired analysis with correction
+# Unpaired stats at each timepoint b/w groups: wilcoxon UNpaired analysis
+# with correction
 #######################################################################
 # with adjustment: fdr and BH are identical
 my_adjust_method = "BH"
@@ -56,8 +63,8 @@ sum(is.na(lf_t1$value))
 
 foo = lf_t1[which(is.na(lf_t1$value)),]
 ci = which(is.na(lf_t1$value))
-
 #lf_t1_comp = lf_t1[-ci,]
+
 lf_t1_comp = lf_t1[-which(is.na(lf_t1$value)),]
 stats_un_t1 = compare_means(value~obesity
                             , group.by = "mediator"
@@ -68,8 +75,8 @@ stats_un_t1 = compare_means(value~obesity
 
 foo$mosaic[!unique(foo$mosaic)%in%unique(lf_t1_comp$mosaic)]
 
+# add timepoint and convert to df
 stats_un_t1$timepoint = "t1"
-
 stats_un_t1 = as.data.frame(stats_un_t1)
 class(stats_un_t1)
 
@@ -78,7 +85,9 @@ n_t1 = data.frame(table(lf_t1_comp$mediator))
 colnames(n_t1) = c("mediator", "n_obs")
 n_t1$mediator = as.character(n_t1$mediator)
 
-# merge stats + n_obs df
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
 merging_cols = intersect(names(stats_un_t1), names(n_t1)); merging_cols
 if (all(n_t1$mediator%in%stats_un_t1$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -115,8 +124,8 @@ stats_un_t2 = compare_means(value~obesity
                             , data = lf_t2_comp
                             , paired = FALSE
                             , p.adjust.method = my_adjust_method)
+# add timepoint and convert to df
 stats_un_t2$timepoint = "t2"
-
 stats_un_t2 = as.data.frame(stats_un_t2)
 class(stats_un_t2)
 
@@ -125,7 +134,9 @@ n_t2 = data.frame(table(lf_t2_comp$mediator))
 colnames(n_t2) = c("mediator", "n_obs")
 n_t2$mediator = as.character(n_t2$mediator)
 
-# merge stats + n_obs df
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
 merging_cols = intersect(names(stats_un_t2), names(n_t2)); merging_cols
 if (all(n_t2$mediator%in%stats_un_t2$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -163,8 +174,8 @@ stats_un_t3 = compare_means(value~obesity
                             , paired = FALSE
                             , p.adjust.method = my_adjust_method)
 
+# add timepoint and convert to df
 stats_un_t3$timepoint = "t3"
-
 stats_un_t3 = as.data.frame(stats_un_t3)
 class(stats_un_t3)
 
@@ -174,7 +185,9 @@ n_t3 = data.frame(table(lf_t3_comp$mediator))
 colnames(n_t3) = c("mediator", "n_obs")
 n_t3$mediator = as.character(n_t3$mediator)
 
-# merge stats + n_obs df
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
 merging_cols = intersect(names(stats_un_t3), names(n_t3)); merging_cols
 if (all(n_t3$mediator%in%stats_un_t3$mediator)) {
   cat("PASS: merging stats and n_obs on column/s:", merging_cols)
@@ -198,7 +211,7 @@ stats_un_t3$p_adj_bonferroni = p.adjust(stats_un_t3$p, method = "bonferroni")
 
 rm(n_t3)
 rm(lf_t3_comp)
-
+########################################################################
 #==============
 # Rbind these dfs
 #==============
@@ -318,8 +331,7 @@ colnames(combined_unpaired_stats_f) = c("mediator"
                                         , "p_bon_signif")
 
 colnames(combined_unpaired_stats_f)
-
-
+########################################################################
 #******************
 # write output file
 #******************