From e3259f2f1764f316f27a2fa534d4d26c74418e2f Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Mon, 23 Nov 2020 18:16:19 +0000 Subject: [PATCH] logistic regression for outcome and meds --- boxplot_linear_na.R | 166 --------- boxplot_log_na.R | 170 --------- data_extraction_formatting.R | 5 +- data_extraction_formatting_clinical.R | 55 ++- data_extraction_formatting_non_asthmatics.R | 355 ------------------ flu_stats_contingency.R | 5 +- flu_stats_unpaired_na_npa.R | 381 ------------------- flu_stats_unpaired_na_sam.R | 383 -------------------- flu_stats_unpaired_na_serum.R | 376 ------------------- logistic_regression.R | 273 +++++++++++--- 10 files changed, 265 insertions(+), 1904 deletions(-) delete mode 100644 boxplot_linear_na.R delete mode 100755 boxplot_log_na.R delete mode 100755 data_extraction_formatting_non_asthmatics.R delete mode 100755 flu_stats_unpaired_na_npa.R delete mode 100755 flu_stats_unpaired_na_sam.R delete mode 100755 flu_stats_unpaired_na_serum.R mode change 100755 => 100644 logistic_regression.R diff --git a/boxplot_linear_na.R b/boxplot_linear_na.R deleted file mode 100644 index 12f883e..0000000 --- a/boxplot_linear_na.R +++ /dev/null @@ -1,166 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd("~/git/mosaic_2020/") -getwd() -############################################################ -# TASK: boxplots at T1 -# FIXME: currently not rendering, problem with NAs for stats? -############################################################ -my_samples = "npa_sam_serum" -#============= -# Input -#============= -#source("data_extraction_formatting_non_asthmatics.R") -source("plot_data_na.R") - -# check: adult variable and age variable discrepancy! -metadata_all$mosaic[metadata_all$adult==1 & metadata_all$age<=18] - -#============= -# Output: -#============= -outfile_bp = paste0("boxplots_linear_NA_", my_samples, ".pdf") -output_boxplot = paste0(outdir_plots, outfile_bp); output_boxplot - -#=============================== -# data assignment for plots -#================================ -#----------- -# npa -#----------- -##wf_fp_npa = npa_wf[npa_wf$flustat == 1,] -#lf_fp_npa = npa_lf[npa_lf$flustat == 1,] -#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",] - -#----------- -# sam -#----------- -##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",] - -#----------- -# serum -#----------- -##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) -#lf_fp_serum$timepoint = as.factor(lf_fp_serum$timepoint) -#lf_fp_serum$obesity = as.factor(lf_fp_serum$obesity) - -######################################################################## -cat("Output plots will be in:", output_boxplot) -pdf(output_boxplot, width = 20, height = 15) - -#======================================================================= -# NPA -#======================================================================= -if (is.factor(lf_fp_npa$timepoint) && is.factor(lf_fp_npa$timepoint)){ - cat ("PASS: required groups are factors") -} -#------------------------------------------ -title_npa_linear = "NPA (Linear)" -#----------------------------------------- -bxp_npa_linear <- ggboxplot(lf_fp_npa, x = "timepoint", y = "value", - color = "obesity", palette = c("#00BFC4", "#F8766D")) + - facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+ - #scale_y_log10() + - theme(axis.text.x = element_text(size = 15) - , axis.text.y = element_text(size = 15 - , angle = 0 - , hjust = 1 - , vjust = 0) - , axis.title.x = element_text(size = 15) - , axis.title.y = element_text(size = 15) - , plot.title = element_text(size = 20, hjust = 0.5) - , strip.text.x = element_text(size = 15, colour = "black") - , legend.title = element_text(color = "black", size = 20) - , legend.text = element_text(size = 15) - , legend.direction = "horizontal") + - labs(title = title_npa_linear - , x = "" - , y = "Levels") - -bxp_npa_linear - -#======================================================================= -# SAM -#======================================================================= -if (is.factor(lf_fp_sam$timepoint) && is.factor(lf_fp_sam$timepoint)){ - cat ("PASS: required groups are factors") -} - -#------------------------------------------ -title_sam_linear = "SAM (Linear)" -#----------------------------------------- -bxp_sam_linear <- ggboxplot(lf_fp_sam, x = "timepoint", y = "value", - color = "obesity", palette = c("#00BFC4", "#F8766D")) + - facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+ - #scale_y_log10() + - theme(axis.text.x = element_text(size = 15) - , axis.text.y = element_text(size = 15 - , angle = 0 - , hjust = 1 - , vjust = 0) - , axis.title.x = element_text(size = 15) - , axis.title.y = element_text(size = 15) - , plot.title = element_text(size = 20, hjust = 0.5) - , strip.text.x = element_text(size = 15, colour = "black") - , legend.title = element_text(color = "black", size = 20) - , legend.text = element_text(size = 15) - , legend.direction = "horizontal") + - labs(title = title_sam_linear - , x = "" - , y = "Levels") - -bxp_sam_linear - -#======================================================================= -# SERUM -#======================================================================= -if (is.factor(lf_fp_serum$timepoint) && is.factor(lf_fp_serum$timepoint)){ - cat ("PASS: required groups are factors") -} - -table(lf_fp_serum$mediator) - -#------------------------------------------ -title_serum_linear = "SERUM (Linear)" -#----------------------------------------- -bxp_serum_linear <- ggboxplot(lf_fp_serum, x = "timepoint", y = "value", - color = "obesity", palette = c("#00BFC4", "#F8766D")) + - facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+ - #scale_y_log10() + - theme(axis.text.x = element_text(size = 15) - , axis.text.y = element_text(size = 15 - , angle = 0 - , hjust = 1 - , vjust = 0) - , axis.title.x = element_text(size = 15) - , axis.title.y = element_text(size = 15) - , plot.title = element_text(size = 20, hjust = 0.5) - , strip.text.x = element_text(size = 15, colour = "black") - , legend.title = element_text(color = "black", size = 20) - , legend.text = element_text(size = 15) - , legend.direction = "horizontal") + - labs(title = title_serum_linear - , x = "" - , y = "Levels") - -bxp_serum_linear - -#========================================================================== -dev.off() -############################################################################ diff --git a/boxplot_log_na.R b/boxplot_log_na.R deleted file mode 100755 index 7eff01c..0000000 --- a/boxplot_log_na.R +++ /dev/null @@ -1,170 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd("~/git/mosaic_2020/") -getwd() -############################################################ -# TASK: boxplots at T1 -# FIXME: currently not rendering, problem with NAs for stats? -############################################################ -my_samples = "npa_sam_serum" -#============= -# Input -#============= -#source("data_extraction_formatting_non_asthmatics.R") -source("plot_data_na.R") - -# check: adult variable and age variable discrepancy! -metadata_all$mosaic[metadata_all$adult==1 & metadata_all$age<=18] - -#============= -# Output: -#============= -outfile_bp_log = paste0("boxplots_log_NA_", my_samples, ".pdf") -output_boxplot_log = paste0(outdir_plots, outfile_bp_log); output_boxplot_log - -#=============================== -# data assignment for plots -#================================ -#----------- -# npa -#----------- -#wf_fp_npa = npa_wf[npa_wf$flustat == 1,] -lf_fp_npa = npa_lf[npa_lf$flustat == 1,] -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",] - -#----------- -# sam -#----------- -#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",] - -#----------- -# serum -#----------- -#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) -lf_fp_serum$timepoint = as.factor(lf_fp_serum$timepoint) -lf_fp_serum$obesity = as.factor(lf_fp_serum$obesity) - -######################################################################## -cat("Output plots will be in:", output_boxplot_log) -pdf(output_boxplot_log, width = 20, height = 15) - -#======================================================================= -# NPA -#======================================================================= -if (is.factor(lf_fp_npa$timepoint) && is.factor(lf_fp_npa$timepoint)){ - cat ("PASS: required groups are factors") -} - -#------------------------------------ -title_npa_log = "NPA (Log)" -#----------------------------------- - -bxp_npa_log <- ggboxplot(lf_fp_npa, x = "timepoint", y = "value", - color = "obesity", palette = c("#00BFC4", "#F8766D")) + - facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = F)+ - scale_y_log10() + - theme(axis.text.x = element_text(size = 15) - , axis.text.y = element_text(size = 15 - , angle = 0 - , hjust = 1 - , vjust = 0) - , axis.title.x = element_text(size = 15) - , axis.title.y = element_text(size = 15) - , plot.title = element_text(size = 20, hjust = 0.5) - , strip.text.x = element_text(size = 15, colour = "black") - , legend.title = element_text(color = "black", size = 20) - , legend.text = element_text(size = 15) - , legend.direction = "horizontal") + - labs(title = title_npa_log - , x = "" - , y = "Levels (Log)") - -bxp_npa_log - -#======================================================================= -# SAM -#======================================================================= -if (is.factor(lf_fp_sam$timepoint) && is.factor(lf_fp_sam$timepoint)){ - cat ("PASS: required groups are factors") -} - -#------------------------------------ -title_sam_log = "SAM (Log)" -#----------------------------------- - -bxp_sam_log <- ggboxplot(lf_fp_sam, x = "timepoint", y = "value", - color = "obesity", palette = c("#00BFC4", "#F8766D")) + - facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+ - scale_y_log10() + - theme(axis.text.x = element_text(size = 15) - , axis.text.y = element_text(size = 15 - , angle = 0 - , hjust = 1 - , vjust = 0) - , axis.title.x = element_text(size = 15) - , axis.title.y = element_text(size = 15) - , plot.title = element_text(size = 20, hjust = 0.5) - , strip.text.x = element_text(size = 15, colour = "black") - , legend.title = element_text(color = "black", size = 20) - , legend.text = element_text(size = 15) - , legend.direction = "horizontal") + - labs(title = title_sam_log - , x = "" - , y = "Levels (Log)") - -bxp_sam_log - -#======================================================================= -# SERUM -#======================================================================= -if (is.factor(lf_fp_serum$timepoint) && is.factor(lf_fp_serum$timepoint)){ - cat ("PASS: required groups are factors") -} - -table(lf_fp_serum$mediator) - -#------------------------------------ -title_serum_log = "SERUM (Log)" -#----------------------------------- - -bxp_serum_log <- ggboxplot(lf_fp_serum, x = "timepoint", y = "value", - color = "obesity", palette = c("#00BFC4", "#F8766D")) + - facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+ - scale_y_log10() + - theme(axis.text.x = element_text(size = 15) - , axis.text.y = element_text(size = 15 - , angle = 0 - , hjust = 1 - , vjust = 0) - , axis.title.x = element_text(size = 15) - , axis.title.y = element_text(size = 15) - , plot.title = element_text(size = 20, hjust = 0.5) - , strip.text.x = element_text(size = 15, colour = "black") - , legend.title = element_text(color = "black", size = 20) - , legend.text = element_text(size = 15) - , legend.direction = "horizontal") + - labs(title = title_serum_log - , x = "" - , y = "Levels (Log)") - -bxp_serum_log - -#========================================================================== -dev.off() -############################################################################ diff --git a/data_extraction_formatting.R b/data_extraction_formatting.R index cdc65a7..ad6a4e3 100755 --- a/data_extraction_formatting.R +++ b/data_extraction_formatting.R @@ -351,4 +351,7 @@ rm(sam_df_adults, serum_df_adults, npa_df_adults) # rm df rm(sam_df, serum_df, npa_df) -rm(colnames_check_f, fp_adults) +rm(colnames_check_f + #, fp_adults + #, fp_adults_na) +) diff --git a/data_extraction_formatting_clinical.R b/data_extraction_formatting_clinical.R index acd5ce1..0aef946 100755 --- a/data_extraction_formatting_clinical.R +++ b/data_extraction_formatting_clinical.R @@ -172,6 +172,12 @@ table(clinical_df$asthma, clinical_df$age_bins) #1 11 8 12 9 sum(table(clinical_df$asthma, clinical_df$age_bins)) == nrow(clinical_df) +table(clinical_df$age_int) +clinical_df = subset(clinical_df, select = -c(age_int)) +table(clinical_df$age_int) + +class(clinical_df$age_bins) +clinical_df$age_bins #=========================== # O2 saturation binning @@ -214,7 +220,6 @@ tot_onset2ini onset_initial_bin = cut(clinical_df$onset_2_initial, c(min_in, 4, max_in)) clinical_df$onset_initial_bin = onset_initial_bin - sum(table(clinical_df$onset_initial_bin)) == tot_onset2ini #======================= @@ -357,10 +362,6 @@ clinical_df_ics = merge(clinical_df, clinical_ics, by = merging_cols, all = T); colnames(clinical_df_ics) -# change colname of logistic_outcome -c1 = which(colnames(clinical_df_ics) == "logistic_outcome") -colnames(clinical_df_ics)[c1] <- "t1_resp_recoded" - if (nrow(clinical_df_ics) == nrow(clinical_df) & nrow(clinical_ics)){ cat("\nPASS: No. of rows match, nrow =", nrow(clinical_df_ics) , "\nChecking ncols...") @@ -376,14 +377,33 @@ if (nrow(clinical_df_ics) == nrow(clinical_df) & nrow(clinical_ics)){ , "\nExpected nrows:", nrow(fp_adults)) } -# change the factor vars to integers -str(clinical_df_ics) -factor_vars = lapply(clinical_df_ics, class) == "factor" -table(factor_vars) -clinical_df_ics[, factor_vars] <- lapply(clinical_df_ics[, factor_vars], as.integer) -table(factor_vars) +#========================= +# add binary outcome for T1 resp score +#========================= +table(clinical_df_ics$T1_resp_score) -str(clinical_df_ics) +clinical_df_ics$t1_resp_recoded = ifelse(clinical_df_ics$T1_resp_score <3, 0, 1) +table(clinical_df_ics$t1_resp_recoded) +#table(clinical_df_ics$steroid) +table(clinical_df_ics$steroid_ics) + +#========================= +# change the factor vars to integers +#========================= +#str(clinical_df_ics) +#factor_vars = lapply(clinical_df_ics, class) == "factor" +#table(factor_vars) + +#clinical_df_ics[, factor_vars] <- lapply(clinical_df_ics[, factor_vars], as.integer) +#table(factor_vars) + +#str(clinical_df_ics) + +#========================= +# remove cols +#========================= + +clinical_df_ics = subset(clinical_df_ics, select = -c(onset_2_initial)) #====================== # writing output file @@ -392,10 +412,17 @@ outfile_name_reg = "clinical_df_recoded.csv" outfile_reg = paste0(outdir, outfile_name_reg) cat("\nWriting clinical file for regression:", outfile_reg) - #write.csv(clinical_df_ics, file = outfile_reg) + +#========================= +# clinical_df_ics: without asthma +#========================= +clinical_df_ics_na = clinical_df_ics[clinical_df_ics$asthma == 0,] + ################################################################ rm(age_bins, max_age_interval, max_in, min_in , o2_sat_bin, onset_initial_bin, tot_o2 + , n_text_code, n1, n2, tot_onset2ini, infile_ics , tot_onset2ini, meta_data_cols - , clinical_df) + , clinical_df, clinical_ics) +################################################################ diff --git a/data_extraction_formatting_non_asthmatics.R b/data_extraction_formatting_non_asthmatics.R deleted file mode 100755 index af1912f..0000000 --- a/data_extraction_formatting_non_asthmatics.R +++ /dev/null @@ -1,355 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd('~/git/mosaic_2020/') -getwd() -######################################################################## -# TASK: Extract relevant columns from mosaic adults data -# sam -# serum -# npa -######################################################################## -#==================== -# Input: source data -#==================== -source("read_data.R") - -#============================ -# Data to use: Important step -#============================ -# select df to use -my_data = fp_adults_na - -# clear unnecessary variables -rm(all_df, adult_df, fp_adults) - -######################################################################## -if ( sum(table(my_data$asthma)) && sum(table(my_data$asthma, my_data$ia_exac_copd)) && sum(table(my_data$obesity)) == nrow(my_data) ){ - cat("PASS: fp adults WITHOUT asthma extracted sucessfully") -}else{ - cat("FAIL: numbers mismatch. Please check") - quit() -} - -#========= -# sam -#========= -sam_regex = regex(".*_sam[1-3]{1}$", ignore_case = T) -sam_cols_i = str_extract(colnames(my_data), sam_regex) # not boolean -#sam_cols_b = colnames(my_data)%in%sam_cols_i # boolean - -sam_cols = colnames(my_data)[colnames(my_data)%in%sam_cols_i] -sam_cols - -# this contains log columns + daysamp_samXX: omitting these -sam_regex_log_days = regex("log|day.*_sam[1-3]{1}$", ignore_case = T, perl = T) -sam_cols_to_omit = sam_cols[grepl(sam_regex_log_days, sam_cols)]; sam_cols_to_omit -sam_cols_to_omit - -sam_cols_clean = sam_cols[!sam_cols%in%sam_cols_to_omit]; sam_cols_clean -length(sam_cols_clean) - -if( length(sam_cols_clean) == length(sam_cols) - length(sam_cols_to_omit) ){ - cat("PASS: clean cols extracted" - , "\nNo. of clean SAM cols to extract:", length(sam_cols_clean)) -}else{ - cat("FAIL: length mismatch. Aborting further cols extraction" - , "Expected length:", length(sam_cols) - length(sam_cols_to_omit) - , "Got:", length(sam_cols_clean) ) -} - -sam_cols_to_extract = c(meta_data_cols, sam_cols_clean) - -cat("Extracting SAM cols + metadata_cols") - -if ( length(sam_cols_to_extract) == length(meta_data_cols) + length(sam_cols_clean) ){ - cat("Extracing", length(sam_cols_to_extract), "columns for sam") - sam_df = my_data[, sam_cols_to_extract] -}else{ - cat("FAIL: length mismatch" - , "Expeceted to extract:", length(meta_data_cols) + length(sam_cols_clean), "columns" - , "Got:", length(sam_cols_to_extract)) -} - -colnames_sam_df = colnames(sam_df); colnames_sam_df - -#========= -# serum -#========= -serum_regex = regex(".*_serum[1-3]{1}$", ignore_case = T) -serum_cols_i = str_extract(colnames(my_data), serum_regex) # not boolean -#serum_cols_b = colnames(my_data)%in%serum_cols_i # boolean - -serum_cols = colnames(my_data)[colnames(my_data)%in%serum_cols_i] - -# this contains log columns + dayserump_serumXX: omitting these -serum_regex_log_days = regex("log|day.*_serum[1-3]{1}$", ignore_case = T, perl = T) -serum_cols_to_omit = serum_cols[grepl(serum_regex_log_days, serum_cols)]; serum_cols_to_omit - -serum_cols_clean = serum_cols[!serum_cols%in%serum_cols_to_omit]; serum_cols_clean -length(serum_cols_clean) - -if( length(serum_cols_clean) == length(serum_cols) - length(serum_cols_to_omit) ){ - cat("PASS: clean cols extracted" - , "\nNo. of clean serum cols to extract:", length(serum_cols_clean)) -}else{ - cat("FAIL: length mismatch. Aborting further cols extraction" - , "Expected length:", length(serum_cols) - length(serum_cols_to_omit) - , "Got:", length(serum_cols_clean) ) -} - -serum_cols_to_extract = c(meta_data_cols, serum_cols_clean) - -cat("Extracting SERUM cols + metadata_cols") - -if ( length(serum_cols_to_extract) == length(meta_data_cols) + length(serum_cols_clean) ){ - cat("Extracing", length(serum_cols_to_extract), "columns for serum") - serum_df = my_data[, serum_cols_to_extract] -}else{ - cat("FAIL: length mismatch" - , "Expeceted to extract:", length(meta_data_cols) + length(serum_cols_clean), "columns" - , "Got:", length(serum_cols_to_extract)) -} - -colnames_serum_df = colnames(serum_df); colnames_serum_df - -#========= -# npa -#========= -npa_regex = regex(".*_npa[1-3]{1}$", ignore_case = T) -npa_cols_i = str_extract(colnames(my_data), npa_regex) # not boolean -#npa_cols_b = colnames(my_data)%in%npa_cols_i # boolean - -npa_cols = colnames(my_data)[colnames(my_data)%in%npa_cols_i] - -# this contains log columns + daynpap_npaXX: omitting these -npa_regex_log_days = regex("log|day|vl_samptime|ct.*_npa[1-3]{1}$", ignore_case = T, perl = T) -npa_cols_to_omit = npa_cols[grepl(npa_regex_log_days, npa_cols)]; npa_cols_to_omit - -npa_cols_clean = npa_cols[!npa_cols%in%npa_cols_to_omit]; npa_cols_clean -length(npa_cols_clean) - -if( length(npa_cols_clean) == length(npa_cols) - length(npa_cols_to_omit) ){ - cat("PASS: clean cols extracted" - , "\nNo. of clean npa cols to extract:", length(npa_cols_clean)) -}else{ - cat("FAIL: length mismatch. Aborting further cols extraction" - , "Expected length:", length(npa_cols) - length(npa_cols_to_omit) - , "Got:", length(npa_cols_clean) ) -} - -npa_cols_to_extract = c(meta_data_cols, npa_cols_clean) - -cat("Extracting NPA cols + metadata_cols") - -if ( length(npa_cols_to_extract) == length(meta_data_cols) + length(npa_cols_clean) ){ - cat("Extracing", length(npa_cols_to_extract), "columns for npa") - npa_df = my_data[, npa_cols_to_extract] -}else{ - cat("FAIL: length mismatch" - , "Expeceted to extract:", length(meta_data_cols) + length(npa_cols_clean), "columns" - , "Got:", length(npa_cols_to_extract)) -} - -colnames_npa_df = colnames(npa_df); colnames_npa_df - -#============== -# quick checks -#============== -colnames_check = as.data.frame(cbind(colnames_sam_df, colnames_serum_df, colnames_npa_df)) -tail(colnames_check) # gives a warning message due to differing no. of rows for cbind! - -# put NA where a match doesn't exist -# unmatched lengths -#colnames_check[117,1] <- NA -#colnames_check[117,2] <- NA - -if ( ncol(sam_df) == ncol(serum_df) ){ - start = ncol(sam_df)+1 - extra_cols = start:ncol(npa_df) -} - -colnames_check_f = colnames_check -tail(colnames_check_f) - -for (i in extra_cols){ - for (j in 1:2) { - cat("\ni:", i - ,"\nj:", j) - colnames_check_f[i,j] <- NA - #colnames_check_f[i, j]< - NA - - } -} -tail(colnames_check_f) - -########################################################################## -# LF data -########################################################################## -cols_to_omit = c("adult" - #, "obese2" - #, "height", "height_unit", "weight" - #, "weight_unit", "visual_est_bmi", "bmi_rating" - ) - -pivot_cols = meta_data_cols -# subselect pivot_cols -pivot_cols = meta_data_cols[!meta_data_cols%in%cols_to_omit];pivot_cols -ncols_omitted = table(meta_data_cols%in%cols_to_omit)[[2]] -ncols_omitted - -#============== -# lf data: sam -#============== -str(sam_df) -table(sam_df$obesity) -#table(sam_df$obese2) - -#sam_df_adults = sam_df[sam_df$adult == 1,] # resolved at source and only dealing wit age as adult -sam_df_adults = sam_df - -wf_cols = colnames(sam_df_adults)[!colnames(sam_df_adults)%in%cols_to_omit] -sam_wf = sam_df_adults[wf_cols] - -if (length(pivot_cols) == length(meta_data_cols) - ncols_omitted){ - cat("PASS: pivot cols successfully extracted") -}else{ - cat("FAIL: length mismatch! pivot cols could not be extracted" - , "\nExpected length:", length(meta_data_cols) - ncols_omitted - , "\nGot:",length(pivot_cols) ) - quit() -} - -expected_rows_sam_lf = nrow(sam_wf) * (length(sam_wf) - length(pivot_cols)); expected_rows_sam_lf - -# using regex: -sam_lf = sam_wf %>% - tidyr::pivot_longer(-all_of(pivot_cols) - , names_to = c("mediator", "sample_type", "timepoint") - , names_pattern = "(.*)_(.*)([1-3]{1})" - , values_to = "value") - -if ( - (nrow(sam_lf) == expected_rows_sam_lf) & (sum(table(is.na(sam_lf$mediator))) == expected_rows_sam_lf) - ) { - cat(paste0("PASS: long format data has correct no. of rows and NA in mediator:" - , "\nNo. of rows: ", nrow(sam_lf) - , "\nNo. of cols: ", ncol(sam_lf))) -} else{ - cat(paste0("FAIL:long format data has unexpected no. of rows or NAs in mediator" - , "\nExpected no. of rows: ", expected_rows_sam_lf - , "\nGot: ", nrow(sam_lf) - , "\ncheck expected rows calculation!")) - quit() -} - -#library(data.table) -#foo = sam_df_adults[1:10] -#long <- melt(setDT(sam_df_adults), id.vars = pivot_cols, variable.name = "levels") - -#============== -# lf data: serum -#============== -str(serum_df) -table(serum_df$obesity) -#table(serum_df$obese2) - -#serum_df_adults = serum_df[serum_df$adult == 1,] # extract based on age -serum_df_adults = serum_df - -wf_cols = colnames(serum_df_adults)[!colnames(serum_df_adults)%in%cols_to_omit] -serum_wf = serum_df_adults[wf_cols] - -if (length(pivot_cols) == length(meta_data_cols) - ncols_omitted){ - cat("PASS: pivot cols successfully extracted") -}else{ - cat("FAIL: length mismatch! pivot cols could not be extracted" - , "\nExpected length:", length(meta_data_cols) - ncols_omitted - , "\nGot:",length(pivot_cols) ) - quit() -} - -expected_rows_serum_lf = nrow(serum_wf) * (length(serum_wf) - length(pivot_cols)); expected_rows_serum_lf - -# using regex: -serum_lf = serum_wf %>% - tidyr::pivot_longer(-all_of(pivot_cols) - , names_to = c("mediator", "sample_type", "timepoint") - , names_pattern = "(.*)_(.*)([1-3]{1})" - , values_to = "value") - -if ( - (nrow(serum_lf) == expected_rows_serum_lf) & (sum(table(is.na(serum_lf$mediator))) == expected_rows_serum_lf) -) { - cat(paste0("PASS: long format data has correct no. of rows and NA in mediator:" - , "\nNo. of rows: ", nrow(serum_lf) - , "\nNo. of cols: ", ncol(serum_lf))) -} else{ - cat(paste0("FAIL:long format data has unexpected no. of rows or NAs in mediator" - , "\nExpected no. of rows: ", expected_rows_serum_lf - , "\nGot: ", nrow(serum_lf) - , "\ncheck expected rows calculation!")) - quit() -} - -#============== -# lf data: npa -#============== -str(npa_df) -table(npa_df$obesity) -#table(npa_df$obese2) - -#npa_df_adults = npa_df[npa_df$adult == 1,] # extract based on age -npa_df_adults = npa_df - -wf_cols = colnames(npa_df_adults)[!colnames(npa_df_adults)%in%cols_to_omit] -npa_wf = npa_df_adults[wf_cols] - -if (length(pivot_cols) == length(meta_data_cols) - ncols_omitted){ - cat("PASS: pivot cols successfully extracted") -}else{ - cat("FAIL: length mismatch! pivot cols could not be extracted" - , "\nExpected length:", length(meta_data_cols) - ncols_omitted - , "\nGot:",length(pivot_cols) ) - quit() -} - -expected_rows_npa_lf = nrow(npa_wf) * (length(npa_wf) - length(pivot_cols)); expected_rows_npa_lf - -# using regex: -npa_lf = npa_wf %>% - tidyr::pivot_longer(-all_of(pivot_cols) - , names_to = c("mediator", "sample_type", "timepoint") - , names_pattern = "(.*)_(.*)([1-3]{1})" - , values_to = "value") - -if ( - (nrow(npa_lf) == expected_rows_npa_lf) & (sum(table(is.na(npa_lf$mediator))) == expected_rows_npa_lf) -) { - cat(paste0("PASS: long format data has correct no. of rows and NA in mediator:" - , "\nNo. of rows: ", nrow(npa_lf) - , "\nNo. of cols: ", ncol(npa_lf))) -} else{ - cat(paste0("FAIL:long format data has unexpected no. of rows or NAs in mediator" - , "\nExpected no. of rows: ", expected_rows_npa_lf - , "\nGot: ", nrow(npa_lf) - , "\ncheck expected rows calculation!")) - quit() -} - -############################################################################### -# remove unnecessary variables -rm(sam_regex, sam_regex_log_days, sam_cols, sam_cols_clean, sam_cols_i, sam_cols_to_extract, sam_cols_to_omit) -rm(serum_regex, serum_regex_log_days, serum_cols, serum_cols_clean, serum_cols_i, serum_cols_to_extract, serum_cols_to_omit) -rm(npa_regex, npa_regex_log_days, npa_cols, npa_cols_clean, npa_cols_i, npa_cols_to_extract, npa_cols_to_omit) -rm(my_data) -rm(colnames_check) -rm(i, j - #, expected_cols - , start, wf_cols, extra_cols, cols_to_omit) - -# rm not_clean dfs -rm(sam_df_adults, serum_df_adults, npa_df_adults) - -# rm df -rm(sam_df, serum_df, npa_df) -rm(colnames_check_f, fp_adults_na) diff --git a/flu_stats_contingency.R b/flu_stats_contingency.R index a191912..b1d8648 100755 --- a/flu_stats_contingency.R +++ b/flu_stats_contingency.R @@ -214,7 +214,7 @@ comb_stats_categ_df_f = comb_stats_categ_df[order(comb_stats_categ_df$p_signif # write output file #****************** cat("Chisq and fishers test results in:", outfile_clin_categ) -#write.csv(comb_stats_categ_df_f, outfile_clin_categ, row.names = FALSE) +write.csv(comb_stats_categ_df_f, outfile_clin_categ, row.names = FALSE) #================== #0 date not recorded @@ -222,7 +222,8 @@ cat("Chisq and fishers test results in:", outfile_clin_categ) #-2 n/a specified by clinician #-3 unknown specified by - +######################################################################## +# checks! chisq_test(table(clinical_df_ics$obesity, clinical_df_ics$smoking)) chisq_test(table(clinical_df_ics$obesity, clinical_df_ics$max_resp_score)) diff --git a/flu_stats_unpaired_na_npa.R b/flu_stats_unpaired_na_npa.R deleted file mode 100755 index 3f5d3de..0000000 --- a/flu_stats_unpaired_na_npa.R +++ /dev/null @@ -1,381 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd("~/git/mosaic_2020/") -getwd() -############################################################ -# 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_non_asthmatics.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_NA_", my_sample_type, ".csv") -flu_stats_time_unpaired_na = 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) -lf = lf[!lf$mediator == "vitd",] -######################################################################## -# clear variables -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) - -# sanity checks -table(lf$timepoint) - -######################################################################## -# 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)) -#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" - #, data = lf_t1 - , data = lf_t1_comp - , paired = FALSE - , p.adjust.method = my_adjust_method) - -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) - -#---------------------------------------- -# 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 -#=================================== -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) - stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t1) - , "\nncol:", ncol(stats_un_t1)) -}else{ - nf = n_t1$mediator[!n_t1$mediator%in%stats_un_t1$mediator] - stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, 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)) -} - -# add bonferroni adjustment as well -stats_un_t1$p_adj_bonferroni = p.adjust(stats_un_t1$p, method = "bonferroni") - -rm(n_t1) -rm(lf_t1_comp) - -#============== -# unpaired: t2 -#============== -lf_t2 = lf[lf$timepoint == "t2",] -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 - , 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) - -#---------------------------------------- -# 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 -#================================== -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) - stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t2) - , "\nncol:", ncol(stats_un_t2)) -}else{ - nf = n_t2$mediator[!n_t2$mediator%in%stats_un_t2$mediator] - stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, 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_t2) - , "\nncol:", ncol(stats_un_t2)) -} - -# add bonferroni adjustment as well -stats_un_t2$p_adj_bonferroni = p.adjust(stats_un_t2$p, method = "bonferroni") - -rm(n_t2) -rm(lf_t2_comp) - -#============== -# unpaired: t3 -#============== -lf_t3 = lf[lf$timepoint == "t3",] -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 - , 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) - -#---------------------------------------- -# 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 -#================================== -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) - stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t3) - , "\nncol:", ncol(stats_un_t3)) -}else{ - nf = n_t3$mediator[!n_t3$mediator%in%stats_un_t3$mediator] - stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, 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_t3) - , "\nncol:", ncol(stats_un_t3)) -} - -# add bonferroni adjustment as well -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 - -if ( all.equal(nrow(stats_un_t1), nrow(stats_un_t2), nrow(stats_un_t3)) && - all.equal(ncol(stats_un_t1), ncol(stats_un_t2), ncol(stats_un_t3)) ) { - expected_rows = nrow(stats_un_t1) * n_dfs - expected_cols = ncol(stats_un_t1) - print("PASS: expected_rows and cols variables generated for downstream sanity checks") -}else{ - cat("FAIL: dfs have different no. of rows and cols" - , "\nCheck harcoded value of n_dfs" - , "\nexpected_rows and cols could not be generated") - quit() -} - -if ( all.equal(colnames(stats_un_t1), colnames(stats_un_t2), colnames(stats_un_t3)) ){ - print("PASS: colnames match. Rbind the 3 dfs...") - combined_unpaired_stats = rbind(stats_un_t1, stats_un_t2, stats_un_t3) -} else{ - cat("FAIL: cannot combined dfs. Colnames don't match!") - quit() -} - -if ( nrow(combined_unpaired_stats) == expected_rows && ncol(combined_unpaired_stats) == expected_cols ){ - cat("PASS: combined_df has expected dimension" - , "\nNo. of rows in combined_df:", nrow(combined_unpaired_stats) - , "\nNo. of cols in combined_df:", ncol(combined_unpaired_stats) ) -}else{ - cat("FAIL: combined_df dimension mismatch") - quit() -} - -####################################################################### -#================= -# formatting df -#================= -# delete: unnecessary column -combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.)) - -# add sample_type -cat("Adding sample type info as a column", my_sample_type, "...") -combined_unpaired_stats$sample_type = my_sample_type - -# add: reflect stats method correctly i.e paired or unpaired -# incase there are NA due to LLODs, the gsub won't work! -#combined_unpaired_stats$method = gsub("Wilcoxon", "Wilcoxon_unpaired", combined_unpaired_stats$method) -combined_unpaired_stats$method = "wilcoxon unpaired" -combined_unpaired_stats$method - -# add an extra column for padjust_signif: my_adjust_method -combined_unpaired_stats$padjust_signif = combined_unpaired_stats$p.adj -# add appropriate symbols for padjust_signif: my_adjust_method -combined_unpaired_stats = dplyr::mutate(combined_unpaired_stats, padjust_signif = case_when(padjust_signif == 0.05 ~ "." - , padjust_signif <=0.0001 ~ '****' - , padjust_signif <=0.001 ~ '***' - , padjust_signif <=0.01 ~ '**' - , padjust_signif <0.05 ~ '*' - , TRUE ~ 'ns')) -# add an extra column for p_bon_signif -combined_unpaired_stats$p_bon_signif = combined_unpaired_stats$p_adj_bonferroni -# add appropriate symbols for p_bon_signif -combined_unpaired_stats = dplyr::mutate(combined_unpaired_stats, p_bon_signif = case_when(p_bon_signif == 0.05 ~ "." - , p_bon_signif <=0.0001 ~ '****' - , p_bon_signif <=0.001 ~ '***' - , p_bon_signif <=0.01 ~ '**' - , p_bon_signif <0.05 ~ '*' - , TRUE ~ 'ns')) -# reorder columns -print("preparing to reorder columns...") -colnames(combined_unpaired_stats) -my_col_order2 = c("mediator" - , "timepoint" - , "sample_type" - , "n_obs" - , "n_obs_complete" - , "group1" - , "group2" - , "method" - , "p" - , "p.format" - , "p.signif" - , "p.adj" - , "padjust_signif" - , "p_adj_bonferroni" - , "p_bon_signif") - -if( length(my_col_order2) == ncol(combined_unpaired_stats) && (all(my_col_order2%in%colnames(combined_unpaired_stats))) ){ - print("PASS: Reordering columns...") - combined_unpaired_stats_f = combined_unpaired_stats[, my_col_order2] - print("Successful: column reordering") - print("formatted df called:'combined_unpaired_stats_f'") - cat('\nformatted df has the following dimensions\n') - print(dim(combined_unpaired_stats_f )) -} else{ - cat(paste0("FAIL:Cannot reorder columns, length mismatch" - , "\nExpected column order for: ", ncol(combined_unpaired_stats) - , "\nGot:", length(my_col_order2))) - quit() -} -# assign nice column names like replace "." with "_" -colnames(combined_unpaired_stats_f) = c("mediator" - , "timepoint" - , "sample_type" - , "n_obs" - , "n_obs_complete" - , "group1" - , "group2" - , "method" - , "p" - , "p_format" - , "p_signif" - , paste0("p_adj_fdr_", my_adjust_method) - , paste0("p_", my_adjust_method, "_signif") - , "p_adj_bonferroni" - , "p_bon_signif") - -colnames(combined_unpaired_stats_f) - -#--------------- -# quick summary -#--------------- -# count how many meds are significant -n_sig = length(combined_unpaired_stats_f$mediator[combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)]) -sig_meds = combined_unpaired_stats_f[(combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)),] - -sig_meds$med_time = paste0(sig_meds$mediator, "@", sig_meds$timepoint) - -cat("\nTotal no. of statistically significant mediators in", toupper(my_sample_type) - , "are:", n_sig - , "\nThese are:", sig_meds$med_time) - -######################################################################## -#****************** -# write output file -#****************** -cat("\nUNpaired stats for groups will be:", flu_stats_time_unpaired_na) -#write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired_na, row.names = FALSE) diff --git a/flu_stats_unpaired_na_sam.R b/flu_stats_unpaired_na_sam.R deleted file mode 100755 index b90ac0d..0000000 --- a/flu_stats_unpaired_na_sam.R +++ /dev/null @@ -1,383 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd("~/git/mosaic_2020/") -getwd() -############################################################ -# 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_non_asthmatics.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_NA_", my_sample_type, ".csv") -flu_stats_time_unpaired_na = paste0(outdir_stats, outfile_name) - -#=============================== -# data assignment for stats -#================================ -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 -, 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)) - -######################################################################## -# 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)) -#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" - #, data = lf_t1 - , data = lf_t1_comp - , paired = FALSE - , p.adjust.method = my_adjust_method) - -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) - -#---------------------------------------- -# 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 -#=================================== -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) - stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t1) - , "\nncol:", ncol(stats_un_t1)) -}else{ - nf = n_t1$mediator[!n_t1$mediator%in%stats_un_t1$mediator] - stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, 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)) -} - -# add bonferroni adjustment as well -stats_un_t1$p_adj_bonferroni = p.adjust(stats_un_t1$p, method = "bonferroni") - -rm(n_t1) -rm(lf_t1_comp) - -#============== -# unpaired: t2 -#============== -lf_t2 = lf[lf$timepoint == "t2",] -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 - , 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) - -#---------------------------------------- -# 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 -#=================================== -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) - stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t2) - , "\nncol:", ncol(stats_un_t2)) -}else{ - nf = n_t2$mediator[!n_t2$mediator%in%stats_un_t2$mediator] - stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, 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_t2) - , "\nncol:", ncol(stats_un_t2)) -} - -# add bonferroni adjustment as well -stats_un_t2$p_adj_bonferroni = p.adjust(stats_un_t2$p, method = "bonferroni") - -rm(n_t2) -rm(lf_t2_comp) - -#============== -# unpaired: t3 -#============== -lf_t3 = lf[lf$timepoint == "t3",] -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 - , 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) - -#---------------------------------------- -# 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 -#================================== -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) - stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t3) - , "\nncol:", ncol(stats_un_t3)) -}else{ - nf = n_t3$mediator[!n_t3$mediator%in%stats_un_t3$mediator] - stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, 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_t3) - , "\nncol:", ncol(stats_un_t3)) -} - -# check: satisfied!!!! -# FIXME: supply the col name automatically? -wilcox.test(wf$ifna2a_sam3[wf$obesity == 1], wf$ifna2a_sam3[wf$obesity == 0]) - -# add bonferroni adjustment as well -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 - -if ( all.equal(nrow(stats_un_t1), nrow(stats_un_t2), nrow(stats_un_t3)) && - all.equal(ncol(stats_un_t1), ncol(stats_un_t2), ncol(stats_un_t3)) ) { - expected_rows = nrow(stats_un_t1) * n_dfs - expected_cols = ncol(stats_un_t1) - print("PASS: expected_rows and cols variables generated for downstream sanity checks") -}else{ - cat("FAIL: dfs have different no. of rows and cols" - , "\nCheck harcoded value of n_dfs" - , "\nexpected_rows and cols could not be generated") - quit() -} - -if ( all.equal(colnames(stats_un_t1), colnames(stats_un_t2), colnames(stats_un_t3)) ){ - print("PASS: colnames match. Rbind the 3 dfs...") - combined_unpaired_stats = rbind(stats_un_t1, stats_un_t2, stats_un_t3) -} else{ - cat("FAIL: cannot combined dfs. Colnames don't match!") - quit() -} - -if ( nrow(combined_unpaired_stats) == expected_rows && ncol(combined_unpaired_stats) == expected_cols ){ - cat("PASS: combined_df has expected dimension" - , "\nNo. of rows in combined_df:", nrow(combined_unpaired_stats) - , "\nNo. of cols in combined_df:", ncol(combined_unpaired_stats) ) -}else{ - cat("FAIL: combined_df dimension mismatch") - quit() -} - -####################################################################### -#================= -# formatting df -#================= -# delete: unnecessary column -combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.)) - -# add sample_type -cat("Adding sample type info as a column", my_sample_type, "...") -combined_unpaired_stats$sample_type = my_sample_type - -# add: reflect stats method correctly i.e paired or unpaired -# incase there are NA due to LLODs, the gsub won't work! -#combined_unpaired_stats$method = gsub("Wilcoxon", "Wilcoxon_unpaired", combined_unpaired_stats$method) -combined_unpaired_stats$method = "wilcoxon unpaired" -combined_unpaired_stats$method - -# add an extra column for padjust_signif: my_adjust_method -combined_unpaired_stats$padjust_signif = combined_unpaired_stats$p.adj -# add appropriate symbols for padjust_signif: my_adjust_method -combined_unpaired_stats = dplyr::mutate(combined_unpaired_stats, padjust_signif = case_when(padjust_signif == 0.05 ~ "." - , padjust_signif <=0.0001 ~ '****' - , padjust_signif <=0.001 ~ '***' - , padjust_signif <=0.01 ~ '**' - , padjust_signif <0.05 ~ '*' - , TRUE ~ 'ns')) -# add an extra column for p_bon_signif -combined_unpaired_stats$p_bon_signif = combined_unpaired_stats$p_adj_bonferroni -# add appropriate symbols for p_bon_signif -combined_unpaired_stats = dplyr::mutate(combined_unpaired_stats, p_bon_signif = case_when(p_bon_signif == 0.05 ~ "." - , p_bon_signif <=0.0001 ~ '****' - , p_bon_signif <=0.001 ~ '***' - , p_bon_signif <=0.01 ~ '**' - , p_bon_signif <0.05 ~ '*' - , TRUE ~ 'ns')) -# reorder columns -print("preparing to reorder columns...") -colnames(combined_unpaired_stats) -my_col_order2 = c("mediator" - , "timepoint" - , "sample_type" - , "n_obs" - , "n_obs_complete" - , "group1" - , "group2" - , "method" - , "p" - , "p.format" - , "p.signif" - , "p.adj" - , "padjust_signif" - , "p_adj_bonferroni" - , "p_bon_signif") - -if( length(my_col_order2) == ncol(combined_unpaired_stats) && (all(my_col_order2%in%colnames(combined_unpaired_stats))) ){ - print("PASS: Reordering columns...") - combined_unpaired_stats_f = combined_unpaired_stats[, my_col_order2] - print("Successful: column reordering") - print("formatted df called:'combined_unpaired_stats_f'") - cat('\nformatted df has the following dimensions\n') - print(dim(combined_unpaired_stats_f )) -} else{ - cat(paste0("FAIL:Cannot reorder columns, length mismatch" - , "\nExpected column order for: ", ncol(combined_unpaired_stats) - , "\nGot:", length(my_col_order2))) - quit() -} -# assign nice column names like replace "." with "_" -colnames(combined_unpaired_stats_f) = c("mediator" - , "timepoint" - , "sample_type" - , "n_obs" - , "n_obs_complete" - , "group1" - , "group2" - , "method" - , "p" - , "p_format" - , "p_signif" - , paste0("p_adj_fdr_", my_adjust_method) - , paste0("p_", my_adjust_method, "_signif") - , "p_adj_bonferroni" - , "p_bon_signif") - -colnames(combined_unpaired_stats_f) - -#--------------- -# quick summary -#--------------- -# count how many meds are significant -n_sig = length(combined_unpaired_stats_f$mediator[combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)]) -sig_meds = combined_unpaired_stats_f[(combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)),] - -sig_meds$med_time = paste0(sig_meds$mediator, "@", sig_meds$timepoint) - -cat("\nTotal no. of statistically significant mediators in", toupper(my_sample_type) - , "are:", n_sig - , "\nThese are:", sig_meds$med_time) - -####################################################################### -#****************** -# write output file -#****************** -cat("\nUNpaired stats for groups will be:", flu_stats_time_unpaired_na) -#write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired_na, row.names = FALSE) diff --git a/flu_stats_unpaired_na_serum.R b/flu_stats_unpaired_na_serum.R deleted file mode 100755 index 97649c3..0000000 --- a/flu_stats_unpaired_na_serum.R +++ /dev/null @@ -1,376 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd("~/git/mosaic_2020/") -getwd() -############################################################ -# TASK: unpaired (time) analysis of mediators: serum -############################################################ -my_sample_type = "serum" - -#============= -# Input -#============= -source("data_extraction_formatting_non_asthmatics.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_NA_", my_sample_type, ".csv") -flu_stats_time_unpaired_na = paste0(outdir_stats, outfile_name) - -#=============================== -# data assignment for stats -#================================ -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) - -######################################################################## -# 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)) -#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" - #, data = lf_t1 - , data = lf_t1_comp - , paired = FALSE - , p.adjust.method = my_adjust_method) - -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) - -#---------------------------------------- -# 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 -#================================== -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) - stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t1) - , "\nncol:", ncol(stats_un_t1)) -}else{ - nf = n_t1$mediator[!n_t1$mediator%in%stats_un_t1$mediator] - stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, 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)) -} - -# add bonferroni adjustment as well -stats_un_t1$p_adj_bonferroni = p.adjust(stats_un_t1$p, method = "bonferroni") - -rm(n_t1) -rm(lf_t1_comp) - -#============== -# unpaired: t2 -#============== -lf_t2 = lf[lf$timepoint == "t2",] -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 - , 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) - -#---------------------------------------- -# 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 -#================================== -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) - stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t2) - , "\nncol:", ncol(stats_un_t2)) -}else{ - nf = n_t2$mediator[!n_t2$mediator%in%stats_un_t2$mediator] - stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, 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_t2) - , "\nncol:", ncol(stats_un_t2)) -} - -# add bonferroni adjustment as well -stats_un_t2$p_adj_bonferroni = p.adjust(stats_un_t2$p, method = "bonferroni") - -rm(n_t2) -rm(lf_t2_comp) - -#============== -# unpaired: t3 -#============== -lf_t3 = lf[lf$timepoint == "t3",] -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 - , 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) - -#---------------------------------------- -# 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 -#================================== -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) - stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all = T) - cat("\nsuccessfull merge:" - , "\nnrow:", nrow(stats_un_t3) - , "\nncol:", ncol(stats_un_t3)) -}else{ - nf = n_t3$mediator[!n_t3$mediator%in%stats_un_t3$mediator] - stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, 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_t3) - , "\nncol:", ncol(stats_un_t3)) -} - -# add bonferroni adjustment as well -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 - -if ( all.equal(nrow(stats_un_t1), nrow(stats_un_t2), nrow(stats_un_t3)) && - all.equal(ncol(stats_un_t1), ncol(stats_un_t2), ncol(stats_un_t3)) ) { - expected_rows = nrow(stats_un_t1) * n_dfs - expected_cols = ncol(stats_un_t1) - print("PASS: expected_rows and cols variables generated for downstream sanity checks") -}else{ - cat("FAIL: dfs have different no. of rows and cols" - , "\nCheck harcoded value of n_dfs" - , "\nexpected_rows and cols could not be generated") - quit() -} - -if ( all.equal(colnames(stats_un_t1), colnames(stats_un_t2), colnames(stats_un_t3)) ){ - print("PASS: colnames match. Rbind the 3 dfs...") - combined_unpaired_stats = rbind(stats_un_t1, stats_un_t2, stats_un_t3) -} else{ - cat("FAIL: cannot combined dfs. Colnames don't match!") - quit() -} - -if ( nrow(combined_unpaired_stats) == expected_rows && ncol(combined_unpaired_stats) == expected_cols ){ - cat("PASS: combined_df has expected dimension" - , "\nNo. of rows in combined_df:", nrow(combined_unpaired_stats) - , "\nNo. of cols in combined_df:", ncol(combined_unpaired_stats) ) -}else{ - cat("FAIL: combined_df dimension mismatch") - quit() -} - -####################################################################### -#================= -# formatting df -#================= -# delete: unnecessary column -combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.)) - -# add sample_type -cat("Adding sample type info as a column", my_sample_type, "...") -combined_unpaired_stats$sample_type = my_sample_type - -# add: reflect stats method correctly i.e paired or unpaired -# incase there are NA due to LLODs, the gsub won't work! -#combined_unpaired_stats$method = gsub("Wilcoxon", "Wilcoxon_unpaired", combined_unpaired_stats$method) -combined_unpaired_stats$method = "wilcoxon unpaired" -combined_unpaired_stats$method - -# add an extra column for padjust_signif: my_adjust_method -combined_unpaired_stats$padjust_signif = combined_unpaired_stats$p.adj -# add appropriate symbols for padjust_signif: my_adjust_method -combined_unpaired_stats = dplyr::mutate(combined_unpaired_stats, padjust_signif = case_when(padjust_signif == 0.05 ~ "." - , padjust_signif <=0.0001 ~ '****' - , padjust_signif <=0.001 ~ '***' - , padjust_signif <=0.01 ~ '**' - , padjust_signif <0.05 ~ '*' - , TRUE ~ 'ns')) -# add an extra column for p_bon_signif -combined_unpaired_stats$p_bon_signif = combined_unpaired_stats$p_adj_bonferroni -# add appropriate symbols for p_bon_signif -combined_unpaired_stats = dplyr::mutate(combined_unpaired_stats, p_bon_signif = case_when(p_bon_signif == 0.05 ~ "." - , p_bon_signif <=0.0001 ~ '****' - , p_bon_signif <=0.001 ~ '***' - , p_bon_signif <=0.01 ~ '**' - , p_bon_signif <0.05 ~ '*' - , TRUE ~ 'ns')) -# reorder columns -print("preparing to reorder columns...") -colnames(combined_unpaired_stats) -my_col_order2 = c("mediator" - , "timepoint" - , "sample_type" - , "n_obs" - , "n_obs_complete" - , "group1" - , "group2" - , "method" - , "p" - , "p.format" - , "p.signif" - , "p.adj" - , "padjust_signif" - , "p_adj_bonferroni" - , "p_bon_signif") - -if( length(my_col_order2) == ncol(combined_unpaired_stats) && (all(my_col_order2%in%colnames(combined_unpaired_stats))) ){ - print("PASS: Reordering columns...") - combined_unpaired_stats_f = combined_unpaired_stats[, my_col_order2] - print("Successful: column reordering") - print("formatted df called:'combined_unpaired_stats_f'") - cat('\nformatted df has the following dimensions\n') - print(dim(combined_unpaired_stats_f )) -} else{ - cat(paste0("FAIL:Cannot reorder columns, length mismatch" - , "\nExpected column order for: ", ncol(combined_unpaired_stats) - , "\nGot:", length(my_col_order2))) - quit() -} -# assign nice column names like replace "." with "_" -colnames(combined_unpaired_stats_f) = c("mediator" - , "timepoint" - , "sample_type" - , "n_obs" - , "n_obs_complete" - , "group1" - , "group2" - , "method" - , "p" - , "p_format" - , "p_signif" - , paste0("p_adj_fdr_", my_adjust_method) - , paste0("p_", my_adjust_method, "_signif") - , "p_adj_bonferroni" - , "p_bon_signif") - -colnames(combined_unpaired_stats_f) -#--------------- -# quick summary -#--------------- -# count how many meds are significant -n_sig = length(combined_unpaired_stats_f$mediator[combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)]) -sig_meds = combined_unpaired_stats_f[(combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)),] - -sig_meds$med_time = paste0(sig_meds$mediator, "@", sig_meds$timepoint) - -cat("\nTotal no. of statistically significant mediators in", toupper(my_sample_type) - , "are:", n_sig - , "\nThese are:", sig_meds$med_time) - - -######################################################################## -#****************** -# write output file -#****************** -cat("\nUNpaired stats for groups will be:", flu_stats_time_unpaired_na) -#write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired_na, row.names = FALSE) diff --git a/logistic_regression.R b/logistic_regression.R old mode 100755 new mode 100644 index d483858..e7f7a8f --- a/logistic_regression.R +++ b/logistic_regression.R @@ -18,80 +18,241 @@ getwd() #==================== source("data_extraction_formatting_clinical.R") -rm(fp_adults, metadata_all) +# quick sanity checks +table(clinical_df_ics$ia_exac_copd==1 & clinical_df_ics$asthma == 1) +table(fp_adults$ia_exac_copd==1 & fp_adults$asthma == 1) +table(fp_adults_na$ia_exac_copd==1 & fp_adults_na$asthma == 1) + +table(clinical_df_ics$asthma) + +#-------------------- +# Data reassignment +#-------------------- +my_data = clinical_df_ics +my_data_na = clinical_df_ics_na + +table(my_data$ia_exac_copd==1 & my_data$asthma == 1) +table(my_data_na$ia_exac_copd==1 & my_data_na$asthma == 1) + +# clear variables +#rm(fp_adults, fp_adults_na) -######################################################################## -my_data = reg_data ######################################################################### -# check factor of each column -lapply(my_data, class) -character_vars <- lapply(my_data, class) == "character" -character_vars -table(character_vars) +if ( names(which(lapply(my_data, class) == "character")) == "mosaic" ){ + cat("Character class for 1 column only:", "mosaic") +}else{ + cat("More than one character class detected: Resolve!") + quit() +} -factor_vars <- lapply(my_data, class) == "factor" -table(factor_vars) +#============================ +# Identifying column types: Reg data +#=========================== +cols_to_omit = c("mosaic", "flustat", "onset_2_initial", "ia_exac_copd") -my_data[, character_vars] <- lapply(my_data[, character_vars], as.factor) -factor_vars <- lapply(my_data, class) == "factor" -factor_vars +my_reg_data = my_data[!colnames(my_data)%in%cols_to_omit] + +my_vars = colnames(my_reg_data) +my_vars + +lapply(my_reg_data, class) +numerical_vars = c("age" + , "vl_pfu_ul_npa1" + , "los" + , "onset2final" + , "onsfindeath" + , "o2_sat_admis") + +my_reg_data[numerical_vars] <- lapply(my_reg_data[numerical_vars], as.numeric) + +my_reg_params = my_vars + +na_count = sapply(my_reg_data, function(x) sum(is.na(x)));na_count +names(na_count)[na_count>0] + +vars_to_factor = my_vars[!my_vars%in%numerical_vars] + +# convert to factor +lapply(my_reg_data, class) +my_reg_data[vars_to_factor] <- lapply(my_reg_data[vars_to_factor], as.factor) +factor_vars <- colnames(my_reg_data)[lapply(my_reg_data, class) == "factor"] table(factor_vars) # check again -lapply(my_data, class) +lapply(my_reg_data, class) -table(my_data$ethnicity) -my_data$ethnicity = as.factor(my_data$ethnicity) -class(my_data$ethnicity) -colnames(my_data) -reg_param = c("age" - , "age_bins" - #, "death" # outcome - , "asthma" - , "obesity" - , "gender" - , "los" - , "o2_sat_admis" - #, "logistic_outcome" - #, "steroid_ics" - , "ethnicity" - , "smoking" - , "sfluv" - , "h1n1v" - , "ia_cxr" - , "max_resp_score" - , "T1_resp_score" - , "com_noasthma" - , "onset_initial_bin") +# all parasm for reg +my_reg_params = c("age" + , "vl_pfu_ul_npa1" + , "los" + , "onset2final" + #, "onsfindeath" + #, "o2_sat_admis" + , "death" + , "obesity" + , "sfluv" + , "h1n1v" + , "gender" + , "asthma" + , "ethnicity" + , "smoking" + , "ia_cxr" + , "max_resp_score" + , "T1_resp_score" + , "com_noasthma" + , "T2_resp_score" + , "inresp_sev" + , "steroid" + , "age_bins" + , "o2_sat_bin" + , "onset_initial_bin" + , "steroid_ics" + , "t1_resp_recoded") -for(i in reg_param) { - # print (i) - p_form = as.formula(paste("death ~ ", i ,sep = "")) - model_reg = glm(p_form , family = binomial, data = my_data) +#================= +# reg data prepare +#================= +pv1 = "death" +pv2 = "t1_resp_recoded" + +#reg_params1 = factor_vars[!factor_vars%in%pv1] +#reg_params_mixed = my_vars[!my_vars%in%pv1] + +######################################################################## +#================= +# outcome2 +#================= +#----------------------------- +# outcome: death + obesity +# data: fp adults +#----------------------------- +my_reg_params1 = my_reg_params[!my_reg_params%in%c("death", "obesity")] + +for(i in my_reg_params1) { + #print (i) + p_form = as.formula(paste("death ~ obesity + ", i ,sep = "")) + print(p_form) + model_reg = glm(p_form , family = binomial, data = my_reg_data) print(summary(model_reg)) print(exp(cbind(OR = coef(model_reg), confint(model_reg)))) #print (PseudoR2(model_reg)) + print(nobs(model_reg)) + cat("=================================================================================\n") +} + +#----------------------------- +# outcome: death +# data: fp adults +#----------------------------- +my_reg_params1v2 = my_reg_params[!my_reg_params%in%c("death")] + +for(i in my_reg_params1v2) { + #print (i) + p_form = as.formula(paste("death ~ ", i ,sep = "")) + print(p_form) + model_reg = glm(p_form , family = binomial, data = my_reg_data) + print(summary(model_reg)) + print(exp(cbind(OR = coef(model_reg), confint(model_reg)))) + #print (PseudoR2(model_reg)) + print(nobs(model_reg)) + cat("=================================================================================\n") +} +######################################################################## +#================= +# outcome2 +#================= +#----------------------------- +# outcome: t1_resp_recoded + obesity +# data: fp adults +#----------------------------- +my_reg_params2 = my_reg_params[!my_reg_params%in%c("death" + , "obesity" + , "t1_resp_recoded" + , "T1_resp_score")] + +for(i in my_reg_params2) { + #print (i) + p_form = as.formula(paste("t1_resp_recoded ~ obesity + ", i ,sep = "")) + print(p_form) + model_reg = glm(p_form , family = binomial, data = my_reg_data) + print(summary(model_reg)) + print(exp(cbind(OR = coef(model_reg), confint(model_reg)))) + #print (PseudoR2(model_reg)) + print(nobs(model_reg)) cat("=================================================================================\n") } -full_mod = glm(death ~ asthma + - gender + - age_bins + - los + - #ethnicity + - onset_initial_bin + - o2_sat_bin + - com_noasthma + +#----------------------------- +# outcome: t1_resp_recoded +# data: fp adults +#----------------------------- +my_reg_params2v2 = my_reg_params[!my_reg_params%in%c("death" + #, "obesity" + , "t1_resp_recoded" + , "T1_resp_score")] + +for(i in my_reg_params2v2) { + #print (i) + p_form = as.formula(paste("t1_resp_recoded ~ ", i ,sep = "")) + print(p_form) + model_reg = glm(p_form , family = binomial, data = my_reg_data) + print(summary(model_reg)) + print(exp(cbind(OR = coef(model_reg), confint(model_reg)))) + #print (PseudoR2(model_reg)) + print(nobs(model_reg)) + cat("=================================================================================\n") +} + +######################################################################## +# Full model +######################################################################## + +full_mod = glm(death ~ obesity + + age + + #age_bins + obesity + - #ia_cxr + - smoking + - #sfluv + - #h1n1v - max_resp_score + - T1_resp_score + - , family = "binomial", data = my_data) + asthma + + t1_resp_recoded + + #ia_cxr + , family = "binomial", data = my_reg_data) summary(full_mod) + + + +######################################################################## +# mediators +######################################################################## +sig_npa_cols = c("mosaic", sig_npa_cols) + +my_med_sig = fp_adults[, sig_npa_cols] + +my_reg_data_med = merge(clinical_df_ics, my_med_sig + , by = intersect(names(clinical_df_ics), names(my_med_sig)) +) + +#my_reg_params_meds = c(my_reg_params, sig_npa_cols) +my_reg_params_meds = colnames(my_reg_data_med) +my_reg_params_meds1 = my_reg_params_meds[!my_reg_params_meds%in%c("mosaic", "flustat" + , "onset_2_initial" + , "onsfindeath" + , "ia_exac_copd" + , "death" + , "obesity")] + + + +for(i in my_reg_params_meds1) { + #print (i) + p_form = as.formula(paste("death ~ obesity + ", i ,sep = "")) + print(p_form) + model_reg = glm(p_form , family = binomial, data = my_reg_data_med) + print(summary(model_reg)) + print(exp(cbind(OR = coef(model_reg), confint(model_reg)))) + #print (PseudoR2(model_reg)) + print(nobs(model_reg)) + cat("=================================================================================\n") +} \ No newline at end of file