####################################################################################
###                    ANALYSIS OF PHS VITAMIN K DATASET                         ###
###                THOMAS WILLIAMS UNI EDINBURGH FEBRUARY 2023                   ###
####################################################################################
#load pdplyr
library(dplyr)
library(tidyr)

#set working directory
setwd(""")

#import data for pre and post pandemic
VitK_overall_dataset <- read.table("vitK_rates_pre_post_pandemic.txt",
                                        header = T, sep = "\t")


#Fisher's exact test to look at rates pre and post pandemic
oral_none_pre_pandemic <- sum(VitK_overall_dataset$X2018[2:3]) + sum(VitK_overall_dataset$X2019[2:3])
IM_pre_pandemic <- VitK_overall_dataset$X2018[1] + VitK_overall_dataset$X2019[1]
oral_none_post_pandemic <- sum(VitK_overall_dataset$X2021[2:3])
IM_post_pandemic <- VitK_overall_dataset$X2021[1]


#Fisher's exact test to look at dataset
fisher_test_pre_post <- fisher.test(rbind(c(IM_pre_pandemic,IM_post_pandemic),c(oral_none_pre_pandemic,oral_none_post_pandemic)),
                                    alternative = "two.sided")

#Fisher's exact test to compare 2018 to other years
oral_none_2018 <- sum(VitK_overall_dataset$X2018[2:3])
IM_2018 <- VitK_overall_dataset$X2018[1]
oral_none_2019_2021 <-  sum(VitK_overall_dataset$X2019[2:3]) + sum(VitK_overall_dataset$X2020[2:3])+ sum(VitK_overall_dataset$X2021[2:3])
IM_2019_2021        <-  VitK_overall_dataset$X2019[1]+VitK_overall_dataset$X2020[1]+VitK_overall_dataset$X2021[1]

fisher_test_2018 <- fisher.test(rbind(c(oral_none_2018,oral_none_2019_2021),c( IM_2018,IM_2019_2021)),
                                    alternative = "two.sided")

#Fisher's exact test to compare 2019 to 2020/21
oral_none_2019 <- sum(VitK_overall_dataset$X2019[2:3])
IM_2019 <- VitK_overall_dataset$X2019[1]
oral_none_2020_2021 <-  sum(VitK_overall_dataset$X2020[2:3])+ sum(VitK_overall_dataset$X2021[2:3])
IM_2020_2021        <-VitK_overall_dataset$X2020[1]+VitK_overall_dataset$X2021[1]

fisher_test_2019 <- fisher.test(rbind(c(oral_none_2019,oral_none_2020_2021),c(IM_2019,IM_2020_2021)),
                                alternative = "two.sided")

#ethnicity
#import data
ethnicity_data <- read.table("VitaminK_Analysis_HV_ethnicity_aggregate.txt", 
                             header = T, sep = "\t")

#run chi squared test on the dataset
chisq.test(ethnicity_data$Oral_none,ethnicity_data$IM)

#English as a first language
English_first_language <- read.table("VitaminK_Analysis_english_first_language_for_testing.txt",
                                     header = T, sep = "\t")
#run Fisher's exact test
fisher.test(rbind(c(English_first_language$oral_none[1],English_first_language$oral_none[2]),
                  c(English_first_language$IM[1],English_first_language$IM[2])))


#SIMD data
#import data
SIMD_data <- read.table("VitaminK_Analysis_SES_for_testing.txt", 
                             header = T, sep = "\t")

#run chi squared test on the dataset
chisq.test(SIMD_data$oral_none,SIMD_data$IM)

#maternal age data
#import data
maternal_age_data <- read.table("VitaminK_Analysis_maternal_age_for_testing.txt", 
                        header = T, sep = "\t")

#run chi squared test on the dataset
chisq.test(maternal_age_data$oral_none,maternal_age_data$IM)