From a2dbecd03d2db2d60ec67362b7c3f1e2cafcab1d Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Fri, 30 Oct 2020 13:04:45 +0000 Subject: [PATCH] updated un[aired stats with n_obs_complete and removed mediator vitd from sam and npa --- boxplot.R | 17 ++++---- flu_stats_unpaired_npa.R | 71 ++++++++++++++++++++++++--------- flu_stats_unpaired_sam.R | 73 ++++++++++++++++++++++++---------- flu_stats_unpaired_serum.R | 67 ++++++++++++++++++++++--------- summary_stats_time_group_npa.R | 2 +- 5 files changed, 165 insertions(+), 65 deletions(-) diff --git a/boxplot.R b/boxplot.R index acf827d..c669c70 100644 --- a/boxplot.R +++ b/boxplot.R @@ -30,12 +30,22 @@ lf_fp_npa$timepoint = paste0("t", lf_fp_npa$timepoint) lf_fp_npa$timepoint = as.factor(lf_fp_npa$timepoint) lf_fp_npa$obesity = as.factor(lf_fp_npa$obesity) +table(lf_fp_npa$mediator) +head(lf_fp_npa$value[lf_fp_npa$mediator == "vitd"]) +lf_fp_npa = lf_fp_npa[!lf_fp_npa$mediator == "vitd",] + + #wf_fp_sam = samm_wf[samm_wf$flustat == 1,] lf_fp_sam = sam_lf[sam_lf$flustat == 1,] lf_fp_sam$timepoint = paste0("t", lf_fp_sam$timepoint) lf_fp_sam$timepoint = as.factor(lf_fp_sam$timepoint) lf_fp_sam$obesity = as.factor(lf_fp_sam$obesity) +table(lf_fp_sam$mediator) +head(lf_fp_sam$value[lf_fp_sam$mediator == "vitd"]) +lf_fp_sam = lf_fp_sam[!lf_fp_sam$mediator == "vitd",] + + #wf_fp_serum = serum_wf[serum_wf$flustat == 1,] lf_fp_serum = serum_lf[serum_lf$flustat == 1,] lf_fp_serum$timepoint = paste0("t", lf_fp_serum$timepoint) @@ -52,11 +62,6 @@ pdf(output_boxplot, width = 20, height = 15) if (is.factor(lf_fp_npa$timepoint) && is.factor(lf_fp_npa$timepoint)){ cat ("PASS: required groups are factors") } - -table(lf_fp_npa$mediator) -head(lf_fp_npa$value[lf_fp_npa$mediator == "vitd"]) -lf_fp_npa = lf_fp_npa[!lf_fp_npa$mediator == "vitd",] - #------------------------------------------ title_npa_linear = "Boxplot: NPA (Linear)" #----------------------------------------- @@ -114,8 +119,6 @@ bxp_npa_log if (is.factor(lf_fp_sam$timepoint) && is.factor(lf_fp_sam$timepoint)){ cat ("PASS: required groups are factors") } -table(lf_fp_sam$mediator) -lf_fp_sam = lf_fp_sam[!lf_fp_sam$mediator == "vitd",] #------------------------------------------ title_sam_linear = "Boxplot: sam (Linear)" diff --git a/flu_stats_unpaired_npa.R b/flu_stats_unpaired_npa.R index ff04e8e..ebd8647 100755 --- a/flu_stats_unpaired_npa.R +++ b/flu_stats_unpaired_npa.R @@ -30,7 +30,7 @@ flu_stats_time_unpaired = paste0(outdir_stats, outfile_name) wf = npa_wf[npa_wf$flustat == 1,] lf = npa_lf[npa_lf$flustat == 1,] lf$timepoint = paste0("t", lf$timepoint) - +lf = lf[!lf$mediator == "vitd",] ######################################################################## # clear variables rm(sam_lf, sam_wf @@ -43,12 +43,12 @@ rm(pivot_cols) # sanity checks table(lf$timepoint) -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) - quit() -} +#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) +# quit() +#} ######################################################################## # Unpaired stats at each timepoint b/w groups: wilcoxon UNpaired analysis # with correction @@ -58,6 +58,7 @@ my_adjust_method = "BH" #============== # unpaired: t1 +#============== lf_t1 = lf[lf$timepoint == "t1",] sum(is.na(lf_t1$value)) @@ -80,10 +81,20 @@ stats_un_t1$timepoint = "t1" stats_un_t1 = as.data.frame(stats_un_t1) class(stats_un_t1) -# calculate n_obs for each mediator -n_t1 = data.frame(table(lf_t1_comp$mediator)) -colnames(n_t1) = c("mediator", "n_obs") -n_t1$mediator = as.character(n_t1$mediator) +#---------------------------------------- +# calculate n_obs for each mediator: t1 +#---------------------------------------- +#n_t1 = data.frame(table(lf_t1_comp$mediator)) +n_t1_all = data.frame(table(lf_t1$mediator)) +colnames(n_t1_all) = c("mediator", "n_obs") +n_t1_all$mediator = as.character(n_t1_all$mediator) + +n_t1_comp = data.frame(table(lf_t1_comp$mediator)) +colnames(n_t1_comp) = c("mediator", "n_obs_complete") +n_t1_comp$mediator = as.character(n_t1_comp$mediator) + +merge_cols = intersect(names(n_t1_all), names(n_t1_comp)); merge_cols +n_t1= merge(n_t1_all, n_t1_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -130,10 +141,20 @@ stats_un_t2$timepoint = "t2" stats_un_t2 = as.data.frame(stats_un_t2) class(stats_un_t2) -# calculate n_obs for each mediator -n_t2 = data.frame(table(lf_t2_comp$mediator)) -colnames(n_t2) = c("mediator", "n_obs") -n_t2$mediator = as.character(n_t2$mediator) +#---------------------------------------- +# calculate n_obs for each mediator: t2 +#---------------------------------------- +#n_t2 = data.frame(table(lf_t2_comp$mediator)) +n_t2_all = data.frame(table(lf_t2$mediator)) +colnames(n_t2_all) = c("mediator", "n_obs") +n_t2_all$mediator = as.character(n_t2_all$mediator) + +n_t2_comp = data.frame(table(lf_t2_comp$mediator)) +colnames(n_t2_comp) = c("mediator", "n_obs_complete") +n_t2_comp$mediator = as.character(n_t2_comp$mediator) + +merge_cols = intersect(names(n_t2_all), names(n_t2_comp)); merge_cols +n_t2= merge(n_t2_all, n_t2_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -180,10 +201,20 @@ stats_un_t3$timepoint = "t3" stats_un_t3 = as.data.frame(stats_un_t3) class(stats_un_t3) -# calculate n_obs for each mediator -n_t3 = data.frame(table(lf_t3_comp$mediator)) -colnames(n_t3) = c("mediator", "n_obs") -n_t3$mediator = as.character(n_t3$mediator) +#---------------------------------------- +# calculate n_obs for each mediator: t3 +#---------------------------------------- +#n_t3 = data.frame(table(lf_t3_comp$mediator)) +n_t3_all = data.frame(table(lf_t3$mediator)) +colnames(n_t3_all) = c("mediator", "n_obs") +n_t3_all$mediator = as.character(n_t3_all$mediator) + +n_t3_comp = data.frame(table(lf_t3_comp$mediator)) +colnames(n_t3_comp) = c("mediator", "n_obs_complete") +n_t3_comp$mediator = as.character(n_t3_comp$mediator) + +merge_cols = intersect(names(n_t3_all), names(n_t3_comp)); merge_cols +n_t3= merge(n_t3_all, n_t3_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -291,6 +322,7 @@ my_col_order2 = c("mediator" , "timepoint" , "sample_type" , "n_obs" + , "n_obs_complete" , "group1" , "group2" , "method" @@ -320,6 +352,7 @@ colnames(combined_unpaired_stats_f) = c("mediator" , "timepoint" , "sample_type" , "n_obs" + , "n_obs_complete" , "group1" , "group2" , "method" diff --git a/flu_stats_unpaired_sam.R b/flu_stats_unpaired_sam.R index 2c0c67e..0797766 100755 --- a/flu_stats_unpaired_sam.R +++ b/flu_stats_unpaired_sam.R @@ -30,6 +30,7 @@ flu_stats_time_unpaired = paste0(outdir_stats, outfile_name) wf = sam_wf[sam_wf$flustat == 1,] lf = sam_lf[sam_lf$flustat == 1,] lf$timepoint = paste0("t", lf$timepoint) +lf = lf[!lf$mediator == "vitd",] ######################################################################## # clear variables rm(npa_lf, npa_wf @@ -43,12 +44,12 @@ rm(pivot_cols) table(lf$timepoint) length(unique(lf$mosaic)) -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() -} +#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 @@ -70,8 +71,8 @@ foo = lf_t1[which(is.na(lf_t1$value)),] lf_t1_comp = lf_t1[-which(is.na(lf_t1$value)),] stats_un_t1 = compare_means(value~obesity , group.by = "mediator" - #, data = lf_t1 - , data = lf_t1_comp + , data = lf_t1 + #, data = lf_t1_comp , paired = FALSE , p.adjust.method = my_adjust_method) @@ -82,10 +83,20 @@ stats_un_t1$timepoint = "t1" stats_un_t1 = as.data.frame(stats_un_t1) class(stats_un_t1) -# calculate n_obs for each mediator -n_t1 = data.frame(table(lf_t1_comp$mediator)) -colnames(n_t1) = c("mediator", "n_obs") -n_t1$mediator = as.character(n_t1$mediator) +#---------------------------------------- +# calculate n_obs for each mediator: t1 +#---------------------------------------- +#n_t1 = data.frame(table(lf_t1_comp$mediator)) +n_t1_all = data.frame(table(lf_t1$mediator)) +colnames(n_t1_all) = c("mediator", "n_obs") +n_t1_all$mediator = as.character(n_t1_all$mediator) + +n_t1_comp = data.frame(table(lf_t1_comp$mediator)) +colnames(n_t1_comp) = c("mediator", "n_obs_complete") +n_t1_comp$mediator = as.character(n_t1_comp$mediator) + +merge_cols = intersect(names(n_t1_all), names(n_t1_comp)); merge_cols +n_t1= merge(n_t1_all, n_t1_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -131,10 +142,20 @@ stats_un_t2$timepoint = "t2" stats_un_t2 = as.data.frame(stats_un_t2) class(stats_un_t2) -# calculate n_obs for each mediator -n_t2 = data.frame(table(lf_t2_comp$mediator)) -colnames(n_t2) = c("mediator", "n_obs") -n_t2$mediator = as.character(n_t2$mediator) +#---------------------------------------- +# calculate n_obs for each mediator: t2 +#---------------------------------------- +#n_t2 = data.frame(table(lf_t2_comp$mediator)) +n_t2_all = data.frame(table(lf_t2$mediator)) +colnames(n_t2_all) = c("mediator", "n_obs") +n_t2_all$mediator = as.character(n_t2_all$mediator) + +n_t2_comp = data.frame(table(lf_t2_comp$mediator)) +colnames(n_t2_comp) = c("mediator", "n_obs_complete") +n_t2_comp$mediator = as.character(n_t2_comp$mediator) + +merge_cols = intersect(names(n_t2_all), names(n_t2_comp)); merge_cols +n_t2= merge(n_t2_all, n_t2_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -180,10 +201,20 @@ stats_un_t3$timepoint = "t3" stats_un_t3 = as.data.frame(stats_un_t3) class(stats_un_t3) -# calculate n_obs for each mediator -n_t3 = data.frame(table(lf_t3_comp$mediator)) -colnames(n_t3) = c("mediator", "n_obs") -n_t3$mediator = as.character(n_t3$mediator) +#---------------------------------------- +# calculate n_obs for each mediator: t3 +#---------------------------------------- +#n_t3 = data.frame(table(lf_t3_comp$mediator)) +n_t3_all = data.frame(table(lf_t3$mediator)) +colnames(n_t3_all) = c("mediator", "n_obs") +n_t3_all$mediator = as.character(n_t3_all$mediator) + +n_t3_comp = data.frame(table(lf_t3_comp$mediator)) +colnames(n_t3_comp) = c("mediator", "n_obs_complete") +n_t3_comp$mediator = as.character(n_t3_comp$mediator) + +merge_cols = intersect(names(n_t3_all), names(n_t3_comp)); merge_cols +n_t3= merge(n_t3_all, n_t3_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -294,6 +325,7 @@ my_col_order2 = c("mediator" , "timepoint" , "sample_type" , "n_obs" + , "n_obs_complete" , "group1" , "group2" , "method" @@ -323,6 +355,7 @@ colnames(combined_unpaired_stats_f) = c("mediator" , "timepoint" , "sample_type" , "n_obs" + , "n_obs_complete" , "group1" , "group2" , "method" diff --git a/flu_stats_unpaired_serum.R b/flu_stats_unpaired_serum.R index b61e032..a069975 100755 --- a/flu_stats_unpaired_serum.R +++ b/flu_stats_unpaired_serum.R @@ -67,8 +67,8 @@ ci = which(is.na(lf_t1$value)) lf_t1_comp = lf_t1[-which(is.na(lf_t1$value)),] stats_un_t1 = compare_means(value~obesity , group.by = "mediator" - #, data = lf_t1 - , data = lf_t1_comp + , data = lf_t1 + #, data = lf_t1_comp , paired = FALSE , p.adjust.method = my_adjust_method) @@ -79,10 +79,20 @@ stats_un_t1$timepoint = "t1" stats_un_t1 = as.data.frame(stats_un_t1) class(stats_un_t1) -# calculate n_obs for each mediator -n_t1 = data.frame(table(lf_t1_comp$mediator)) -colnames(n_t1) = c("mediator", "n_obs") -n_t1$mediator = as.character(n_t1$mediator) +#---------------------------------------- +# calculate n_obs for each mediator: t1 +#---------------------------------------- +#n_t1 = data.frame(table(lf_t1_comp$mediator)) +n_t1_all = data.frame(table(lf_t1$mediator)) +colnames(n_t1_all) = c("mediator", "n_obs") +n_t1_all$mediator = as.character(n_t1_all$mediator) + +n_t1_comp = data.frame(table(lf_t1_comp$mediator)) +colnames(n_t1_comp) = c("mediator", "n_obs_complete") +n_t1_comp$mediator = as.character(n_t1_comp$mediator) + +merge_cols = intersect(names(n_t1_all), names(n_t1_comp)); merge_cols +n_t1= merge(n_t1_all, n_t1_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -119,8 +129,8 @@ lf_t2_comp = lf_t2[-which(is.na(lf_t2$value)),] stats_un_t2 = compare_means(value~obesity , group.by = "mediator" - #, data = lf_t2 - , data = lf_t2_comp + , data = lf_t2 + #, data = lf_t2_comp , paired = FALSE , p.adjust.method = my_adjust_method) # add timepoint and convert to df @@ -128,10 +138,20 @@ stats_un_t2$timepoint = "t2" stats_un_t2 = as.data.frame(stats_un_t2) class(stats_un_t2) -# calculate n_obs for each mediator -n_t2 = data.frame(table(lf_t2_comp$mediator)) -colnames(n_t2) = c("mediator", "n_obs") -n_t2$mediator = as.character(n_t2$mediator) +#---------------------------------------- +# calculate n_obs for each mediator: t2 +#---------------------------------------- +#n_t2 = data.frame(table(lf_t2_comp$mediator)) +n_t2_all = data.frame(table(lf_t2$mediator)) +colnames(n_t2_all) = c("mediator", "n_obs") +n_t2_all$mediator = as.character(n_t2_all$mediator) + +n_t2_comp = data.frame(table(lf_t2_comp$mediator)) +colnames(n_t2_comp) = c("mediator", "n_obs_complete") +n_t2_comp$mediator = as.character(n_t2_comp$mediator) + +merge_cols = intersect(names(n_t2_all), names(n_t2_comp)); merge_cols +n_t2= merge(n_t2_all, n_t2_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -168,8 +188,8 @@ lf_t3_comp = lf_t3[-which(is.na(lf_t3$value)),] stats_un_t3 = compare_means(value~obesity , group.by = "mediator" - #, data = lf_t3 - , data = lf_t3_comp + , data = lf_t3 + #, data = lf_t3_comp , paired = FALSE , p.adjust.method = my_adjust_method) @@ -178,11 +198,20 @@ stats_un_t3$timepoint = "t3" stats_un_t3 = as.data.frame(stats_un_t3) class(stats_un_t3) +#---------------------------------------- +# calculate n_obs for each mediator: t3 +#---------------------------------------- +#n_t3 = data.frame(table(lf_t3_comp$mediator)) +n_t3_all = data.frame(table(lf_t3$mediator)) +colnames(n_t3_all) = c("mediator", "n_obs") +n_t3_all$mediator = as.character(n_t3_all$mediator) -# calculate n_obs for each mediator -n_t3 = data.frame(table(lf_t3_comp$mediator)) -colnames(n_t3) = c("mediator", "n_obs") -n_t3$mediator = as.character(n_t3$mediator) +n_t3_comp = data.frame(table(lf_t3_comp$mediator)) +colnames(n_t3_comp) = c("mediator", "n_obs_complete") +n_t3_comp$mediator = as.character(n_t3_comp$mediator) + +merge_cols = intersect(names(n_t3_all), names(n_t3_comp)); merge_cols +n_t3 = merge(n_t3_all, n_t3_comp, by = merge_cols, all = T) #================================== # Merge: merge stats + n_obs df @@ -289,6 +318,7 @@ my_col_order2 = c("mediator" , "timepoint" , "sample_type" , "n_obs" + , "n_obs_complete" , "group1" , "group2" , "method" @@ -318,6 +348,7 @@ colnames(combined_unpaired_stats_f) = c("mediator" , "timepoint" , "sample_type" , "n_obs" + , "n_obs_complete" , "group1" , "group2" , "method" diff --git a/summary_stats_time_group_npa.R b/summary_stats_time_group_npa.R index 0a47227..cd2f4f6 100755 --- a/summary_stats_time_group_npa.R +++ b/summary_stats_time_group_npa.R @@ -31,7 +31,7 @@ wf = npa_wf[npa_wf$flustat == 1,] lf = npa_lf[npa_lf$flustat == 1,] lf$timepoint = paste0("t", lf$timepoint) lf = lf[!lf$mediator == "vitd",] - +######################################################################## #======================================================= # summary stats by timepoint and obesity: each mediator