#Load dataset from github
url <- "https://github.com/dayoonkwon/BioAge/blob/master/data/NHANES4.rda?raw=true"
download.file(url, destfile = "NHANES4.rda", mode = "wb")
load("NHANES4.rda")
library(dplyr)
library(RNHANES)
library(ggplot2)
library(stringr)
library(survey)
NHANES4 <- NHANES4 %>%filter(year %in% c(2015, 2017))

#Renaming subject ID variable as they have different names in 2 datasets
NHANES4 <- NHANES4 %>%
  mutate(sampleID = str_remove(sampleID, "^\\d+_"))
NHANES4 <- NHANES4 %>%rename(SEQN = sampleID)
NHANES4 <- NHANES4 %>%mutate(SEQN = as.character(SEQN))



#DEMOGRAPHIC DATA
DEMO1516 <- nhanes_load_data("DEMO_I", "2015-2016")
DEMO1718 <- nhanes_load_data("DEMO_J", "2017-2018")
SMQ1516 <- nhanes_load_data("SMQ_I", "2015-2016")
SMQ1718 <- nhanes_load_data("SMQ_J", "2017-2018")
DEMO <- bind_rows( DEMO1516, DEMO1718)
SMQ <- bind_rows( SMQ1516, SMQ1718)
SMQ <- SMQ %>%mutate(SEQN = as.character(SEQN))
SMQ <- SMQ %>% select(SEQN, SMQ020)
SMQ <- SMQ %>%mutate(SMQ020 = case_when(SMQ020 == 1 ~ "YES",SMQ020 %in% c(2, 7, 9) ~ "NO",
    TRUE ~ as.character(SMQ020)))
DEMO <- DEMO %>%mutate(SEQN = as.character(SEQN))

NHANES4 <- left_join(DEMO, NHANES4 %>% select(SEQN, wave, phenoage0, phenoage_advance0), 
                     by = "SEQN")
NHANES4 <- NHANES4[,c("SEQN","SDDSRVYR","RIAGENDR","RIDAGEYR","RIDRETH1","DMDEDUC2","DMDMARTL","INDFMPIR","SDMVPSU","SDMVSTRA","phenoage0","phenoage_advance0")]

NHANES4 <- NHANES4 %>% rename(WAVE = SDDSRVYR, GENDER = RIAGENDR, AGE = RIDAGEYR, RACE = RIDRETH1, EDU = DMDEDUC2, MARITAL = DMDMARTL, `POVERTY_RATIO` = INDFMPIR)
NHANES4 <- NHANES4 %>%left_join(SMQ, by = "SEQN")

#removing missing data for bioage
NHANES4 <- NHANES4 %>%filter(!is.na(phenoage0))

#Discretizing demographic variables
NHANES4 <- NHANES4 %>%mutate(
  GENDER = factor(GENDER, levels = c(1, 2), labels = c("Male", "Female")),
  AGE = case_when(AGE <= 50 ~ "50 and below",AGE > 50 ~ "Above 50"),
  RACE = case_when(RACE == 3 ~ "White",RACE == 4 ~ "Black",TRUE ~ "Others"),
  EDU = case_when(EDU == 1 ~ "No Highschool",EDU %in% c(2, 3, 4) ~ "High School",EDU == 5 ~ "College Graduate"),
  MARITAL = case_when(MARITAL %in% c(1, 6) ~ "Mar or Coh",TRUE ~ "Single"),
  POVERTY_RATIO = case_when(POVERTY_RATIO <= 1.5 ~ "Low",POVERTY_RATIO > 1.5 & POVERTY_RATIO<= 3.0 ~ "Medium",POVERTY_RATIO > 3.0 ~ "High"))



#Blood metal data

PBCD1516 <- nhanes_load_data("PBCD_I", "2015-2016")
PBCD1718 <- nhanes_load_data("PBCD_J", "2017-2018")
PBCD1718 <- PBCD1718 %>%mutate(SEQN = as.character(SEQN))
PBCD1516 <- PBCD1516 %>%mutate(SEQN = as.character(SEQN))
PBCD1718 <- PBCD1718 %>%left_join(DEMO %>% select(SEQN, WTMEC2YR), by = "SEQN")
PBCD1718 <- PBCD1718%>%rename( WTSH2YR = WTMEC2YR)
PBCD <- bind_rows( PBCD1516, PBCD1718)

PBCD.vars <- PBCD[,c("SEQN","WTSH2YR","LBXBPB","LBXBCD","LBXBMN")]

PBCD.vars <- PBCD.vars %>%mutate(LBXBPB = LBXBPB * 10)

PBCD.vars <- PBCD.vars %>% mutate(CD = log(LBXBCD),PB = log(LBXBPB), MN = log(LBXBMN))

PBCD.vars <- PBCD.vars %>%filter(!is.na(MN))


PBCD.vars <- PBCD.vars %>%mutate(SEQN = as.character(SEQN))

BNHANES4 <- inner_join(NHANES4, PBCD.vars, by = "SEQN")

#quartiles
BNHANES4 <- BNHANES4 %>%
  mutate(
    MN_quartile = ntile(MN, 4),
    CD_quartile = ntile(CD, 4),
    PB_quartile = ntile(PB, 4))

# Adjust the weights by dividing by 2
BNHANES4 <- BNHANES4 %>%mutate(WTSH2YR= WTSH2YR / 2)

# Create a survey design object
svy_design_BN <- svydesign(
  id = ~SDMVPSU, 
  strata = ~SDMVSTRA, 
  weights = ~WTSH2YR, 
  data = BNHANES4, 
  nest = TRUE)

# Perform weighted linear regression for MN quartiles 
regression_results <- list()
regression_results$MN <- svyglm(phenoage_advance0 ~ factor(MN_quartile)  + GENDER + AGE + RACE + EDU +SMQ020+ MARITAL + POVERTY_RATIO, design = svy_design_BN)
summary_MN <- summary(regression_results$MN)

# Perform weighted linear regression for CD quartiles 
regression_results$CD <- svyglm(phenoage_advance0 ~ factor(CD_quartile)  + GENDER + AGE + RACE + EDU +SMQ020+ MARITAL + POVERTY_RATIO, design = svy_design_BN)
summary_CD <- summary(regression_results$CD)

# Perform weighted linear regression for PB quartiles
regression_results$PB <- svyglm(phenoage_advance0 ~ factor(PB_quartile)  + GENDER + AGE+ RACE + EDU +SMQ020+ MARITAL + POVERTY_RATIO, design = svy_design_BN)
summary_PB <- summary(regression_results$PB)

# Print summaries
print("Summary for MN:")
print(summary_MN)
print("Summary for CD:")
print(summary_CD)
print("Summary for PB:")
print(summary_PB)