#!/usr/bin/Rscript
############################################################
# TASK: to check unpaired stats from WF data
############################################################
#=============
# Input
#=============
source("plot_data.R")
#%%========================================================
# data assignment for stats
wf_npa = wf_fp_npa
wf_sam = wf_fp_sam
wf_serum = wf_fp_serum
############################################################
# Numbers
table(fp_adults_ics$asthma, fp_adults_ics$obesity)

table(fp_adults_ics$obesity)


#=========
# t1
#=========
wilcox.test(wf_npa$eotaxin3_npa1[wf_npa$obesity == 1], wf_npa$eotaxin3_npa1[wf_npa$obesity == 0], paired = F)
length(wf_npa$eotaxin3_npa1) - sum(is.na(wf_npa$eotaxin3_npa1))

wilcox.test(wf_npa$il12p70_npa1[wf_npa$obesity == 1], wf_npa$il12p70_npa1[wf_npa$obesity == 0], paired = F)
length(wf_npa$il12p70_npa1) - sum(is.na(wf_npa$il12p70_npa1))

tapply(wf_npa$eotaxin3_npa1, wf_npa$obesity, mean, na.rm = T)
tapply(wf_npa$eotaxin3_npa1, wf_npa$obesity, median, na.rm = T)


npa_wf_a = npa_wf[npa_wf$asthma == 1,] 
npa_wf_na = npa_wf[npa_wf$asthma == 0,]

wilcox.test(npa_wf_a$il4_npa1[npa_wf_a$obesity == 1]
            , npa_wf_a$il4_npa1[npa_wf_a$obesity == 0], paired = F)

length(npa_wf_a$il4_npa1) - sum(is.na(npa_wf_a$il4_npa1))

wilcox.test(npa_wf_a$il13_npa1[npa_wf_a$obesity == 1]
            , npa_wf_a$il13_npa1[npa_wf_a$obesity == 0], paired = F)

wilcox.test(npa_wf_na$il4_npa1[npa_wf_na$obesity == 1]
            , npa_wf_na$il4_npa1[npa_wf_na$obesity == 0], paired = F)
wilcox.test(npa_wf_na$il13_npa1[npa_wf_na$obesity == 1]
            , npa_wf_na$il13_npa1[npa_wf_na$obesity == 0], paired = F)

#=========
# t2
#=========
wilcox.test(wf_npa$il4_npa2[wf_npa$obesity == 1], wf_npa$il4_npa2[wf_npa$obesity == 0], paired = F)
length(wf_npa$il4_npa2) - sum(is.na(wf_npa$il4_npa2))

tapply(wf_npa$il4_npa2, wf_npa$obesity, mean, na.rm = T)
tapply(wf_npa$il4_npa2, wf_npa$obesity, median, na.rm = T)

#=========
# t3
#=========
wilcox.test(wf_npa$ip10_npa3[wf_npa$obesity == 1], wf_npa$ip10_npa3[wf_npa$obesity == 0], paired = F)
length(wf_npa$ip10_npa3) - sum(is.na(wf_npa$ip10_npa3))
###################################################################
##########
# sam
##########
#=========
# t1
#=========
wilcox.test(wf_sam$eotaxin3_sam1[wf_sam$obesity == 1], wf_sam$eotaxin3_sam1[wf_sam$obesity == 0], paired = F)
length(wf_sam$eotaxin3_sam1) - sum(is.na(wf_sam$eotaxin3_sam1))

wilcox.test(wf_sam$il12p70_sam1[wf_sam$obesity == 1], wf_sam$il12p70_sam1[wf_sam$obesity == 0], paired = F)
length(wf_sam$il12p70_sam1) - sum(is.na(wf_sam$il12p70_sam1))

tapply(wf_sam$il13_sam1, wf_sam$obesity, mean, na.rm = T)
tapply(wf_sam$il13_sam1, wf_sam$obesity, median, na.rm = T)

#=========
# t2
#=========
wilcox.test(wf_sam$il5_sam2[wf_sam$obesity == 1], wf_sam$il5_sam2[wf_sam$obesity == 0], paired = F)
length(wf_sam$il5_sam2) - sum(is.na(wf_sam$il5_sam2))

tapply(wf_sam$il29_sam2, wf_sam$obesity, mean, na.rm = T)
tapply(wf_sam$il29_sam2, wf_sam$obesity, median, na.rm = T)

#=========
# t3
#=========
wilcox.test(wf_sam$tnfr1_sam3[wf_sam$obesity == 1], wf_sam$tnfr1_sam3[wf_sam$obesity == 0], paired = F)
length(wf_sam$tnfr1_sam3) - sum(is.na(wf_sam$tnfr1_sam3))

###################################################################
##########
# serum
##########

#=========
# t1
#=========
wilcox.test(wf_serum$eotaxin3_serum1[wf_serum$obesity == 1], wf_serum$eotaxin3_serum1[wf_serum$obesity == 0], paired = F)
length(wf_serum$eotaxin3_serum1) - sum(is.na(wf_serum$eotaxin3_serum1))

wilcox.test(wf_serum$neopterin_serum1[wf_serum$obesity == 1], wf_serum$neopterin_serum1[wf_serum$obesity == 0], paired = F)
length(wf_serum$neopterin_serum1) - sum(is.na(wf_serum$neopterin_serum1))

tapply(wf_serum$eotaxin3_serum1, wf_serum$obesity, mean, na.rm = T)
tapply(wf_serum$eotaxin3_serum1, wf_serum$obesity, median, na.rm = T)

#=========
# t2
#=========
wilcox.test(wf_serum$il8_2_serum2[wf_serum$obesity == 1], wf_serum$il8_2_serum2[wf_serum$obesity == 0], paired = F)
length(wf_serum$il8_2_serum2) - sum(is.na(wf_serum$il8_2_serum2))

#=========
# t3
#=========
wilcox.test(wf_serum$rantes_serum3[wf_serum$obesity == 1], wf_serum$rantes_serum3[wf_serum$obesity == 0], paired = F)
length(wf_serum$rantes_serum3) - sum(is.na(wf_serum$rantes_serum3))

tapply(wf_serum$itac_serum3, wf_serum$obesity, mean, na.rm = T)
tapply(wf_serum$itac_serum3, wf_serum$obesity, median, na.rm = T)