trychecking if summary stats may be added to the output

flu_stats_unpaired_npa.R

@@ -3,39 +3,47 @@ getwd()
 setwd("~/git/mosaic_2020/")
 getwd()
 ############################################################
-# TASK: unpaired (time) analysis of mediators: NPA 
+# TASK: unpaired (time) analysis of mediators:
+# sample type: NPA
+# data: Flu positive adult patients
+# group: obesity
 ############################################################
+
+my_sample_type = "npa"
+
 #=============
 # 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]
 
+#=========================================================
+# 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_adults_lf, sam_df_adults_clean
-, serum_adults_lf, serum_df_adults_clean) 
+rm(sam_lf, sam_wf
+, serum_lf, serum_wf) 
 rm(colnames_sam_df, expected_rows_sam_lf
 , colnames_serum_df, expected_rows_serum_lf)
 
 rm(pivot_cols)
 
-my_sample_type = "npa"
 #=============
 # 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)
-
-#%%========================================================
-# data assignment for stats
-wf = npa_df_adults_clean[npa_df_adults_clean$flustat == 1,]
-lf = npa_adults_lf[npa_adults_lf$flustat == 1,]
-#%%========================================================
+########################################################################
+# sanity checks
 table(lf$timepoint)
-lf$timepoint = paste0("t", lf$timepoint)
 
-if (table(lf$flustat) == table(npa_adults_lf$flustat)[[2]]){
+if (table(lf$flustat) == table(npa_lf$flustat)[[2]]){
   cat("Analysing Flu positive patients for:", my_sample_type)
 }else{
   cat("FAIL: problem with subsetting data for:", my_sample_type)
@@ -59,6 +67,9 @@ 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 
@@ -68,8 +79,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 +89,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
+#=========================
+# Merge1: 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)
@@ -100,6 +113,39 @@ 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)