tidy scripts and using age as criteria for adults

This commit is contained in:
Tanushree Tunstall 2020-10-29 13:10:39 +00:00
parent 49c18d390f
commit 1e6e9eaef1
3 changed files with 110 additions and 115 deletions

View file

@ -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
#******************