metals_asd.qmd

---
title: "ASD Urine Metals"
format: html
editor: source
---

# Packages

```{r load}
library(sas7bdat)
library(openxlsx)
library(dplyr)
library(reshape2)
library(ggplot2)
library(ggplot2)
library(ggrepel)
library(corrplot)
library(logbin)
library(epitools)
library(here)
```

# Load data files

```{r load}
# urine.metals <- readRDS(file=file.path(here(),'Data/urine_metals.RDS'))
# spec.grav <- readRDS(file=file.path(here(),'Data/spgrav.rds'))

urine.metals <- read.csv(file.path(here(), "../Data/urine_metals.csv"), header=T)

marbearl <- read.sas7bdat(file=file.path(here(),'../Data/e_m_covars_v5_20191115.sas7bdat'))
marb.cov <- read.sas7bdat(file=file.path(here(),'../Data/marbles_allcov_09mar20.sas7bdat'))

metals <- unique(urine.metals$Analyte.Code)
metals <- as.list(metals)
names(metals) <- unlist(metals)

samps.earli <- read.csv(file.path(here(),'../../maturine_12_15_2018.csv'))
samps.marb <- read.xlsx(file.path(here(),'../../IBC_list.xlsx'))
samps.earli <- samps.earli[!duplicated(samps.earli$Sample_ID),]
samps.marb <- samps.marb[!duplicated(samps.marb$ibc),]
```


# histograms of distribution

```{r}
#specific gravity adjusted values
pdf(file.path(here(),'../Plots/urine_metals_distributions_sg.pdf'))
for(metal in metals){
  sub <- urine.metals %>% filter(Analyte.Code==metal)
  p1 <- ggplot(sub, aes(x=Concentration_sg, fill=cohort)) +
    geom_histogram(alpha=0.6, position='identity') +
    theme_bw() +
    ggtitle(metal)
  print(p1)
}
dev.off()

#log specific gravity adjusted values
pdf(file.path(here(),'../Plots/urine_metals_distributions_sg_log_imp.pdf'))
for(metal in metals){
  sub <- urine.metals %>% filter(Analyte.Code==metal)
  p1 <- ggplot(sub, aes(x=Concentration_imp_sg_log2, fill=cohort)) +
    geom_histogram(alpha=0.6, position='identity') +
    theme_bw() +
    ggtitle(metal)
  print(p1)
}
dev.off()
```

# Merge in ASD data and remove multiple births

```{r dx_alg}
samples <- urine.metals[!duplicated(urine.metals$SID),]
samples <- samples[,c("SID","cohort","fam","child_id","Time")]
table(samples$cohort)
  # EARLI MARBLES 
  #   360     431
table(samples$fam %in% marbearl$famid) 
 #  TRUE 
 #   791

##### counting
e.count <- samples %>% filter(cohort=='EARLI')
length(unique(e.count$child_id))

m.count <- samples %>% filter(cohort=='MARBLES')
length(unique(m.count$child_id))
#####

marbearl.metal <- marbearl[marbearl$famid %in% samples$fam, ]

multiple <- marbearl.metal[duplicated(marbearl.metal$famid) | duplicated(marbearl.metal$famid, fromLast=T),]
multiple[multiple$project=='EARLI',c('coi_id','famid','n_pregs','Dx_alg','BirthMonth','BirthYear')]

multiple[multiple$project=='MARBLES',c('project','coi_id','famid','n_pregs','Dx_alg','BirthMonth','BirthYear')]
multiple.marb <- marbearl[marbearl$coi_id %in% urine.metals$child_id, ]
multiple.marb <- multiple.marb[duplicated(multiple.marb$famid) | duplicated(multiple.marb$famid, fromLast=T), ]
multiple.marb[multiple.marb$project=='MARBLES',c('project','coi_id','famid','n_pregs','Dx_alg','BirthMonth','BirthYear')]


#no multiples birhts
multi.birth.m <- multiple[multiple$project=='MARBLES',c('project','coi_id','famid','n_pregs','Dx_alg','BirthMonth','BirthYear')] %>%
  group_by(famid) %>% filter(min(n_pregs) == max(n_pregs))
multi.birth.e <- multiple[multiple$project=='EARLI',c('coi_id','famid','n_pregs','Dx_alg','BirthMonth','BirthYear')] %>%
  group_by(famid) %>% filter(min(n_pregs) == max(n_pregs))
nrow(multi.birth.e) #14
nrow(multi.birth.m) #10
multi.birth <- c(unique(as.character(multi.birth.m$famid)), unique(as.character(multi.birth.e$famid)))
marbearl.metal <- marbearl.metal[!marbearl.metal$famid %in% multi.birth,]

urine.nomulti <- urine.metals[!urine.metals$fam %in% multi.birth, ]
length(unique(urine.metals$SID)) #791
length(unique(urine.nomulti$SID)) #771

##### counting
urine.multi <- urine.metals[urine.metals$fam %in% multi.birth, ] %>% filter(Analyte.Code=="Lead")

e.count <- urine.multi %>% filter(cohort=='EARLI')
length(unique(e.count$child_id)) #7
nrow(e.count) #12

m.count <- urine.multi %>% filter(cohort=='MARBLES')
length(unique(m.count$child_id)) #5
nrow(m.count) #8
#####


#ASD outcomet
urine.nomulti$asd <- ifelse(urine.nomulti$cohort=='MARBLES',
                            marbearl.metal[match(urine.nomulti$child_id, marbearl.metal$coi_id),'Dx_alg'] , 
                            NA)
urine.nomulti$asd <- ifelse(urine.nomulti$cohort=='EARLI',
                            marbearl.metal[match(urine.nomulti$fam, marbearl.metal$famid),'Dx_alg'] , 
                            urine.nomulti$asd)
urine.nomulti <- urine.nomulti[!duplicated(paste0(urine.nomulti$SID,urine.nomulti$Analyte.Code)),]


urine.marbles <- urine.metals[urine.metals$cohort=='MARBLES',]
```

# Merge in gestational age at collection

```{r models}
# are sample IDs all accounted for?
table(urine.nomulti$SID %in% samps.earli$Sample_ID | urine.nomulti$SID %in% samps.marb$ibc)
#  TRUE 
# 16962

urine.nomulti$GA_collect <- ifelse(urine.nomulti$cohort=="EARLI",
                                   samps.earli[match(urine.nomulti$SID, samps.earli$Sample_ID), "DCC_Weeks_After_Conception"],
                                   samps.marb[match(urine.nomulti$SID, samps.marb$ibc), "GA_sample_weeks"])

### zero GA incorrect, drop (3 samples)
who.missGA <- urine.nomulti %>% filter(GA_collect == 0)
urine.nomulti <- urine.nomulti %>% filter(GA_collect != 0)

##### counting
e.count <- urine.nomulti %>% filter(cohort=='EARLI') %>% filter(Analyte.Code=='Lead')
length(unique(e.count$child_id)) #180
nrow(e.count) #345

who.missGA <- who.missGA %>% filter(Analyte.Code=="Lead")
who.missGA$child_id %in% e.count$child_id

m.count <- urine.nomulti %>% filter(cohort=='MARBLES') %>% filter(Analyte.Code=='Lead')
length(unique(m.count$child_id)) #262
nrow(m.count) #423
#####

```


# EARLI

```{r metal_asd}
urine.earli <- urine.nomulti[urine.nomulti$cohort=='EARLI',]

### outliers (z score > 5 sd) eliminate
urine.earli$Con_sg_orig <- urine.earli$Concentration_sg
urine.earli$Con_imp_sg_orig <- urine.earli$Concentration_imp_sg
urine.earli$Con_imp_sg_log2_orig <- urine.earli$Concentration_imp_sg_log2
urine.earli <- urine.earli %>% group_by(Analyte.Code) %>%
  mutate(Z_sg=(Concentration_sg-mean(Concentration_sg, na.rm=T))/sd(Concentration_sg, na.rm=T),
         Z_imp_sg=(Concentration_imp_sg-mean(Concentration_imp_sg, na.rm=T))/sd(Concentration_imp_sg, na.rm=T),
         Z_imp_sg_log2=(Concentration_imp_sg_log2-mean(Concentration_imp_sg_log2, na.rm=T))/sd(Concentration_imp_sg_log2, na.rm=T)) %>%
  data.frame
urine.earli$Concentration_sg <- ifelse(urine.earli$Z_sg>5, NA, urine.earli$Concentration_sg)
urine.earli$Concentration_imp_sg <- ifelse(urine.earli$Z_imp_sg>5, NA, urine.earli$Concentration_imp_sg)
urine.earli$Concentration_imp_sg_log2 <- ifelse(urine.earli$Z_imp_sg_log2>5, NA, urine.earli$Concentration_imp_sg_log2)

urine.earli %>% group_by(Analyte.Code) %>%
  summarize(sum(is.na(Concentration_sg))) %>% data.frame
urine.earli %>% group_by(Analyte.Code) %>%
  summarize(sum(is.na(Concentration_imp_sg))) %>% data.frame
urine.earli %>% group_by(Analyte.Code) %>%
  summarize(sum(is.na(Concentration_imp_sg_log2))) %>% data.frame

### split into timepoints

#counting
earli.samples <- urine.earli %>% filter(Analyte.Code=="Antimony")
earli.counting <- earli.samples %>% group_by(fam) %>% summarize(n = n())
table(earli.counting$n)
 #  1   2 
 # 15 165

earli.single <- earli.counting %>% filter(n==1)
earli.two <- earli.counting %>% filter(n==2)

# get for those with two samples
fams.two <- unique(earli.two$fam)

pick <- lapply(fams.two, FUN=function(f){
  urine.earli %>% filter(fam==f) %>% filter(GA_collect == min(GA_collect))
})

pick.last <- lapply(fams.two, FUN=function(f){
  urine.earli %>% filter(fam==f) %>% filter(GA_collect == max(GA_collect))
})

metals.earli <- do.call(rbind, pick)
metals.earli.last <- do.call(rbind, pick.last)

summary(metals.earli$GA_collect)
summary(metals.earli.last$GA_collect)

# separate those with only one sample
fams.one <- unique(earli.single$fam)
earli.samples[earli.samples$fam %in% fams.one, ]

pick <- lapply(fams.one, FUN=function(f){
  urine.earli %>% filter(fam==f) %>% filter(GA_collect < 28)
})

pick.last <- lapply(fams.one, FUN=function(f){
  urine.earli %>% filter(fam==f) %>% filter(GA_collect >= 28)
})

metals.earli.addon <- do.call(rbind, pick)
metals.earli.last.addon <- do.call(rbind, pick.last)

metals.earli <- rbind(metals.earli, metals.earli.addon)
metals.earli.last <- rbind(metals.earli.last, metals.earli.last.addon)

summary(metals.earli$GA_collect)
summary(metals.earli.last$GA_collect)

unique(metals.earli$child_id) %>% length()
#173
unique(metals.earli.last$child_id)  %>% length()
#172


### no overlap GA range

# get for those with two samples
fams.two <- unique(earli.two$fam)

pick <- lapply(fams.two, FUN=function(f){
  urine.earli %>% filter(fam==f) %>% filter(GA_collect < 28 & GA_collect == min(GA_collect))
})

pick.last <- lapply(fams.two, FUN=function(f){
  urine.earli %>% filter(fam==f) %>% filter(GA_collect >= 28 & GA_collect == max(GA_collect))
})

metals.earli.strict <- do.call(rbind, pick)
metals.earli.last.strict <- do.call(rbind, pick.last)
metals.earli.strict <- rbind(metals.earli.strict, metals.earli.addon)
metals.earli.last.strict <- rbind(metals.earli.last.strict, metals.earli.last.addon)

unique(metals.earli.strict$child_id) %>% length()
#161
unique(metals.earli.last.strict$child_id)  %>% length()
#162

metals.earli.strict <- metals.earli %>% filter(GA_collect < 28)
metals.earli.last.strict <- metals.earli.last %>% filter(GA_collect >= 28)
unique(metals.earli.strict$child_id) %>% length() #161
unique(metals.earli.last.strict$child_id) %>% length() #162
```

```{r earli pc corr}
# get wide data set
concentrations.earli <- lapply(metals, FUN=function(x){
  values <- metals.earli %>% filter(Analyte.Code==x)
  concentrations <- values$Con_sg_orig
  names(concentrations) <- values$fam
  concentrations
})
concentrations.earli.last <- lapply(metals, FUN=function(x){
  values <- metals.earli.last %>% filter(Analyte.Code==x)
  concentrations <- values$Con_sg_orig
  names(concentrations) <- values$fam
  concentrations
})

wide.earli <- do.call(cbind.data.frame, concentrations.earli)
wide.earli <- wide.earli[,!colnames(wide.earli) %in% c('Beryllium','Platinum','Tungsten','Uranium','Vanadium')]
rownames(wide.earli) <- names(concentrations.earli[[1]])

wide.earli.last <- do.call(cbind.data.frame, concentrations.earli.last)
wide.earli.last <- wide.earli.last[,!colnames(wide.earli.last) %in% c('Beryllium','Platinum','Tungsten','Uranium','Vanadium')]
rownames(wide.earli.last) <- names(concentrations.earli.last[[1]])


# calculate cors
M.e <- cor(wide.earli, method='spearman', use='pairwise.complete.obs')
Ml.e <- cor(wide.earli.last, method='spearman', use='pairwise.complete.obs')

pdf(file.path(here(),'../Plots/EARLI_metal_correlation.pdf'), height=14, width=14)
  corrplot(M.e, method='ellipse', type='upper', addCoef.col='black')
  corrplot(Ml.e, method='ellipse', type='upper', addCoef.col='black')
dev.off()
```


```{r EARLI descriptives}
### sample characteristics
earl.cov <- read.csv(file.path(here(),"../Data Extraction Codebook 20150203_export_20150203_1422987480_wide_20150204_1423070501.csv"))

get_covars <- function(samps){
  samples <- samps[,c('SID','fam','asd')] %>% unique()
  demo <- marbearl.metal[match(samples$fam, marbearl.metal$famid),]
  demo <- demo[,c('famid','coi_id','coi_gender','maternal_age','GA_deliv_wks','MomEdu','BSRC_group_Sib_36mos')]
  demo$MomEdu <- demo$MomEdu > 3 #college degree or not
  
  covs <- earl.cov[match(samples$fam, earl.cov$Family_ID),]
  covs2 <- marbearl[match(samples$fam, marbearl$famid),]
  
  pd <- merge(samps, demo, by.x='fam', by.y='famid', all.x=T)
  pd$asd <- pd$BSRC_group_Sib_36mos
  pd$asd <- factor(pd$asd, levels=c('TD','Non-TD','ASD'))
  pd$ever_smoke <- covs[match(pd$fam, covs$Family_ID), 'cigarettes_1']
  pd$reg_smoke <- covs[match(pd$fam, covs$Family_ID), 'regularly_2']
  pd$preg_smoke <- covs[match(pd$fam, covs$Family_ID), 'EARLI_HB_smoke']
  pd$site <- covs2$site
  
  pd
}

metals.earli.T2 <- get_covars(metals.earli)
metals.earli.T3 <- get_covars(metals.earli.last)

earli.T2.summ <- metals.earli.T2 %>% group_by(Analyte.Code) %>%
  summarize(mean=mean(Concentration_sg, na.rm=T), sd=sd(Concentration_sg, na.rm=T), median=median(Concentration_sg, na.rm=T),
            IQR=IQR(Concentration_sg, na.rm=T),  above_lod=sum(Concentration>LOD, na.rm=T), n=n(), pct_over_LOD=round(above_lod/n*100,1)) %>% data.frame()

earli.T3.summ <- metals.earli.T3 %>% group_by(Analyte.Code) %>%
  summarize(mean=mean(Concentration_sg, na.rm=T), sd=sd(Concentration_sg, na.rm=T), median=median(Concentration_sg, na.rm=T),
            IQR=IQR(Concentration_sg, na.rm=T), above_lod=sum(Concentration>LOD, na.rm=T), n=n(), pct_over_LOD=round(above_lod/n*100,1)) %>% data.frame()

rounding <- function(x){
  x$mean <- signif(x$mean, 3)
  x$sd <- signif(x$sd, 3)
  x$median <- signif(x$median, 3)
  x$IQR <- signif(x$IQR, 3)
  
  x
}

### count missing
metals.earli.T2$period <- 'early'
metals.earli.T3$period <- 'late'
count.emiss <- rbind(metals.earli.T2 %>% filter(Analyte.Code=="Lead"), metals.earli.T3 %>% filter(Analyte.Code=="Lead"))
count.emiss <- count.emiss %>% filter(!duplicated(count.emiss$child_id))
is.na(count.emiss$asd) %>% sum() #0
is.na(count.emiss$MomEdu) %>% sum() #3
is.na(count.emiss$coi_gender) %>% sum() #0
is.na(count.emiss$maternal_age) %>% sum() #0
is.na(count.emiss$GA_collect) %>% sum() #0

count.emiss <- count.emiss %>% filter(!is.na(MomEdu))
nrow(count.emiss) #177

count.emiss2 <- rbind(metals.earli.T2 %>% filter(Analyte.Code=="Lead"), metals.earli.T3 %>% filter(Analyte.Code=="Lead")) %>% filter(!is.na(MomEdu))
table(count.emiss2$period)
# early  late 
#   170   171
nrow(count.emiss2) #341

### descriptive tables
library(compareGroups)

gather_descriptives <- function(pd){
  descriptives <- pd[!duplicated(pd$fam),]
  descriptives <- descriptives[!is.na(descriptives$asd) & !is.na(descriptives$MomEdu &
                              !is.na(descriptives$coi_gender)) & !is.na(descriptives$maternal_age)
                              & !is.na(descriptives$GA_deliv_wks),]
  
  
  descriptives$MomEdu <- ifelse(descriptives$MomEdu, "College Degree", "No College Degree")
  descriptives$ever_smoke <- factor(descriptives$ever_smoke)
  descriptives$preg_smoke <- factor(descriptives$preg_smoke)
  
  restab <- compareGroups(asd ~ MomEdu + coi_gender + maternal_age + GA_collect + GA_deliv_wks + ever_smoke + preg_smoke,
                          data=descriptives)
  createTable(restab)
}

tab.earli.T2 <- gather_descriptives(metals.earli.T2)
tab.earli.T3 <- gather_descriptives(metals.earli.T3)

export2word(tab.earli.T2, file=file.path(here(),'../Tables/EARLI_early_table1_20Sept2023.docx'))
export2word(tab.earli.T3, file=file.path(here(),'../Tables/EARLI_late_table1_20Sept2023.docx'))


### bivariate tests
lod_percents <- read.csv(file.path(here(),"../../Metals QC/over_LOD_cohort.csv"))
lod.earli <- lod_percents[lod_percents$cohort=='EARLI',]
lod.earli$mod <- ifelse(lod.earli$pct_over_LOD>70, 'continuous', 'binary')
mod.type <- lod.earli$mod
names(mod.type) <- lod.earli$Analyte.Code

IQRs.T2 <- metals.earli.T2 %>% group_by(Analyte.Code) %>% summarize(IQR(Concentration_sg, na.rm=T)) %>% data.frame
IQRs.T2 <- IQRs.T2[,2]
names(IQRs.T2) <- metals

IQRs.T3 <- metals.earli.T3 %>% group_by(Analyte.Code) %>% summarize(IQR(Concentration_sg, na.rm=T)) %>% data.frame
IQRs.T3 <- IQRs.T3[,2]
names(IQRs.T3) <- metals

# pdf(file.path(here(),"../Plots/earli_firstT_violin.pdf"))
# for(met in metals){
#   p1 = ggplot(metals.earli.T2 %>% filter(Analyte.Code==met), aes(x=factor(asd), y=Concentration_sg)) +
#   geom_violin() +geom_jitter(shape=16, position=position_jitter(0.2)) +
#     ggtitle(met) + xlab('') +
#     theme_bw() +
#     theme(axis.text.x=element_text(size=16, face="bold"))
#   print(p1)
# }
# dev.off()

# pdf(file.path(path,"earli_firstT_voilin_smoke.pdf"))
# for(met in metals){
#   p2 = ggplot(metals.earli.T2 %>% filter(Analyte.Code==met), aes(x=factor(ever_smoke), y=Concentration_sg)) +
#   geom_violin() +geom_jitter(shape=16, position=position_jitter(0.2)) +
#     ggtitle(met) + xlab('') +
#     theme_bw() +
#     theme(axis.text.x=element_text(size=16, face="bold"))
#   print(p2)
# }
# dev.off()
```

Early vs Late Preg

```{r}
#remove ones with only single time
table(metals.earli.T2$fam %in% metals.earli.T3$fam)
both <- intersect(unique(metals.earli.T2$fam), unique(metals.earli.T3$fam))

earli.epreg <- metals.earli.T2[metals.earli.T2$fam %in% both, ]
earli.lpreg <- metals.earli.T3[metals.earli.T3$fam %in% both, ]
dim(earli.epreg)
dim(earli.lpreg)
identical(earli.epreg$fam, earli.lpreg$fam)

length(unique(earli.epreg$fam))
length(unique(earli.lpreg$fam))


cor(earli.epreg$Concentration , earli.lpreg$Concentration)
cor(earli.epreg$Concentration_imp_sg, earli.lpreg$Concentration_imp_sg, use="complete.obs")
time.cors <- sapply(metals, FUN=function(x){
  epreg = earli.epreg %>% filter(Analyte.Code==x)
  lpreg = earli.lpreg %>% filter(Analyte.Code==x)
  raw_concentration = round(cor(epreg$Concentration, lpreg$Concentration, 
                                   use="complete.obs", method='spearman'),2)
  imputed_sg_corrected = round(cor(epreg$Concentration_imp_sg, lpreg$Concentration_imp_sg, 
                                   use="complete.obs", method='spearman'),2)
  data.frame(cbind(raw_concentration, imputed_sg_corrected))
})
write.csv(t(time.cors), file.path(here(),"/../Tables/earli_time_cor.csv"))

plot.tm <- function(met){
  m.epreg <- earli.epreg %>% filter(Analyte.Code==met)
  m.lpreg <- earli.lpreg %>% filter(Analyte.Code==met)
  dat <- data.frame(epreg=m.epreg$Con_sg_orig, lpreg=m.lpreg$Con_sg_orig)
  
  ggplot(dat, aes(epreg, lpreg)) +
    geom_point() +
    theme_bw() +
    xlab(paste0("Early Pregnancy ",met)) +
    ylab(paste0("Late Pregnancy ",met))
}

plot.tm('Lead')
plot.tm('Mercury')
plot.tm('Cadmium')
```

# MARBLES

```{r asd_metal_marb}
urine.marbles <- urine.nomulti[urine.nomulti$cohort=='MARBLES',]

### outliers (z score > 5 sd) eliminate
urine.marbles$Con_sg_orig <- urine.marbles$Concentration_sg
urine.marbles$Con_imp_sg_orig <- urine.marbles$Concentration_imp_sg
urine.marbles$Con_imp_sg_log2_orig <- urine.marbles$Concentration_imp_sg_log2
urine.marbles <- urine.marbles %>% group_by(Analyte.Code) %>%
  mutate(Z_sg=(Concentration_sg-mean(Concentration_sg, na.rm=T))/sd(Concentration_sg, na.rm=T),
         Z_imp_sg=(Concentration_imp_sg-mean(Concentration_imp_sg, na.rm=T))/sd(Concentration_imp_sg, na.rm=T),
         Z_imp_sg_log2=(Concentration_imp_sg_log2-mean(Concentration_imp_sg_log2, na.rm=T))/sd(Concentration_imp_sg_log2, na.rm=T)) %>%
  data.frame
urine.marbles$Concentration_sg <- ifelse(urine.marbles$Z_sg>5, NA, urine.marbles$Concentration_sg)
urine.marbles$Concentration_imp_sg <- ifelse(urine.marbles$Z_imp_sg>5, NA, urine.marbles$Concentration_imp_sg)
urine.marbles$Concentration_imp_sg_log2 <- ifelse(urine.marbles$Z_imp_sg_log2>5, NA, urine.marbles$Concentration_imp_sg_log2)

urine.marbles %>% group_by(Analyte.Code) %>%
  summarize(sum(is.na(Concentration_sg))) %>% data.frame
urine.marbles %>% group_by(Analyte.Code) %>%
  summarize(sum(is.na(Concentration_imp_sg))) %>% data.frame
urine.marbles %>% group_by(Analyte.Code) %>%
  summarize(sum(is.na(Concentration_imp_sg_log2))) %>% data.frame

### split into timepoints
urine.marblesT2 <- urine.marbles[urine.marbles$Time=='T2',]
urine.marblesT3 <- urine.marbles[urine.marbles$Time=='T3',]

marbles.counting <- data.frame(fam=unique(urine.marbles$fam))
marbles.counting$T2 <- marbles.counting$fam %in% urine.marblesT2$fam
marbles.counting$T3 <- marbles.counting$fam %in% urine.marblesT3$fam
marbles.counting$N_samp <- marbles.counting$T2 + marbles.counting$T3
table(marbles.counting$N_samp)


### pick earliest/latest measure
coi <- unique(urine.marbles$child_id)

pick <- lapply(coi, FUN=function(id){
  if(id %in% urine.marblesT2$child_id){
    urine.marblesT2[urine.marblesT2$child_id==id,]
  }else if(id %in% urine.marblesT3$child_id){
    urine.marblesT3[urine.marblesT3$child_id==id,]
  }
})

pick.last <- lapply(coi, FUN=function(id){
  if(id %in% urine.marblesT3$child_id){
    urine.marblesT3[urine.marblesT3$child_id==id,]
  }else if(id %in% urine.marblesT2$child_id){
    urine.marblesT2[urine.marblesT2$child_id==id,]
  }
})

metals.marbles <- do.call(rbind, pick)
metals.marbles.last <- do.call(rbind, pick.last)

### separate those with only one sample
onlyone.marbles <- marbles.counting[marbles.counting$N_samp==1,]
onlyone.marbles.early <- onlyone.marbles[!onlyone.marbles$T3,]
onlyone.marbles.late <- onlyone.marbles[onlyone.marbles$T3,]

metals.marbles <- metals.marbles[!metals.marbles$fam %in% onlyone.marbles.late$fam, ]
metals.marbles.last <- metals.marbles.last[!metals.marbles.last$fam %in% onlyone.marbles.early$fam, ]


### split into timepoints

#counting
marbles.samples <- urine.marbles %>% filter(Analyte.Code=="Antimony")
marbles.counting <- marbles.samples %>% group_by(child_id) %>% summarize(n = n())
table(marbles.counting$n)
#   1   2 
# 101 161

marbles.single <- marbles.counting %>% filter(n==1)
marbles.two <- marbles.counting %>% filter(n==2)

# get for those with two samples
coi.two <- unique(marbles.two$child_id)

pick <- lapply(coi.two, FUN=function(c){
  urine.marbles %>% filter(child_id==c) %>% filter(GA_collect == min(GA_collect))
})

pick.last <- lapply(coi.two, FUN=function(c){
  urine.marbles %>% filter(child_id==c) %>% filter(GA_collect == max(GA_collect))
})

metals.marbles <- do.call(rbind, pick)
metals.marbles.last <- do.call(rbind, pick.last)

summary(metals.marbles$GA_collect)
summary(metals.marbles.last$GA_collect)

# separate those with only one sample
coi.one <- unique(marbles.single$child_id)
marbles.samples[marbles.samples$child_id %in% coi.one, "GA_collect"] %>% summary()
marbles.samples[marbles.samples$child_id %in% coi.one, "Time"] %>% table()

pick <- lapply(coi.one, FUN=function(c){
  urine.marbles %>% filter(child_id==c) %>% filter(GA_collect < 28)
})

pick.last <- lapply(coi.one, FUN=function(c){
  urine.marbles %>% filter(child_id==c) %>% filter(GA_collect >= 28)
})

metals.marbles.addon <- do.call(rbind, pick)
metals.marbles.last.addon <- do.call(rbind, pick.last)

metals.marbles <- rbind(metals.marbles, metals.marbles.addon)
metals.marbles.last <- rbind(metals.marbles.last, metals.marbles.last.addon)

summary(metals.marbles$GA_collect)
summary(metals.marbles.last$GA_collect)

## for siblings, pick one randomly
sibs1 <- metals.marbles %>% filter(Analyte.Code=="Antimony") %>% group_by(fam) %>% summarize(n=n())
sibs2 <- metals.marbles.last %>% filter(Analyte.Code=="Antimony") %>% group_by(fam) %>% summarize(n=n())
table(sibs1$n)
table(sibs2$n)
fams1 <- sibs1[sibs1$n>1,]$fam
fams2 <- sibs2[sibs2$n>1,]$fam
fams1[fams1 %in% fams2]
fams2[fams2 %in% fams1]

check1 <- metals.marbles %>% filter(Analyte.Code=="Antimony") %>% filter(fam %in% fams1) %>% select(fam,child_id,asd) 
check1 <- check1[!duplicated(check1$child_id),]
check1[order(check1$child_id),]
check2 <- metals.marbles.last %>% filter(Analyte.Code=="Antimony") %>% filter(fam %in% fams2) %>% select(fam,child_id,asd)
check2 <- check2[!duplicated(check2$child_id),]
check2[order(check2$child_id),]

# select late preg first
set.seed(1453)
keep.late <- check2 %>% group_by(fam) %>% sample_n(1)
drop.late <- check2 %>% filter(!child_id %in% keep.late$child_id)

# select early preg, if relateds were also in late, choose same as the one chosen in that set
keep.early <- check1 %>% group_by(fam) %>% 
  filter(ifelse(fam %in% keep.late$fam, child_id %in% keep.late$child_id, T))
keep.early <- keep.early %>% group_by(fam) %>% sample_n(1)
drop.early <- check1 %>% filter(!child_id %in% keep.early$child_id)

nrow(drop.early) + nrow(drop.late)

# filtering after picking procedure
count.presibfilter1 <- metals.marbles %>% filter(Analyte.Code=="Lead")
count.presibfilter2 <- metals.marbles.last %>% filter(Analyte.Code=="Lead")

metals.marbles <- metals.marbles %>% filter(!child_id %in% drop.early$child_id)
metals.marbles.last <- metals.marbles.last %>% filter(!child_id %in% drop.late$child_id)

# count up whats here
marbles.samples <- rbind(metals.marbles %>% filter(Analyte.Code=="Antimony"), metals.marbles.last %>% filter(Analyte.Code=="Antimony"))
marbles.counting <- marbles.samples %>% group_by(child_id) %>% summarize(n = n())
table(marbles.counting$n)
 # 1   2 
 # 97 154 

marbles.presamples <- rbind(count.presibfilter1, count.presibfilter2)
marbles.precounting <- marbles.presamples %>% group_by(child_id) %>% summarize(n = n())
table(marbles.precounting$n)
#   1   2 
# 101 161

length(unique(marbles.presamples$child_id)) #262
length(unique(marbles.samples$child_id)) #251

length(unique(marbles.presamples$SID)) #423
length(unique(marbles.samples$SID)) #405


### no GA overlap

metals.marbles.strict <- metals.marbles %>% filter(GA_collect < 28)
metals.marbles.last.strict <- metals.marbles.last %>% filter(GA_collect >= 28)

length(unique(metals.marbles$child_id)) #167
length(unique(metals.marbles.strict$child_id)) #167

length(unique(metals.marbles.last$child_id)) #238
length(unique(metals.marbles.last.strict$child_id)) #238

``` 

```{r marbles pc corr}
# get wide data set
concentrations.marbles <- lapply(metals, FUN=function(x){
  values <- metals.marbles %>% filter(Analyte.Code==x)
  concentrations <- values$Con_sg_orig
  names(concentrations) <- values$child_id
  concentrations
})
concentrations.marbles.last <- lapply(metals, FUN=function(x){
  values <- metals.marbles.last  %>% filter(Analyte.Code==x)
  concentrations <- values$Con_sg_orig
  names(concentrations) <- values$child_id
  concentrations
})

wide.marbles <- do.call(cbind.data.frame, concentrations.marbles)
wide.marbles <- wide.marbles[,!colnames(wide.marbles) %in% c('Beryllium','Platinum','Tungsten','Uranium','Vanadium')]
rownames(wide.marbles) <- names(concentrations.marbles[[1]])

wide.marbles.last <- do.call(cbind.data.frame, concentrations.marbles.last)
wide.marbles.last <- wide.marbles.last[,!colnames(wide.marbles.last) %in% c('Beryllium','Platinum','Tungsten','Uranium','Vanadium')]
rownames(wide.marbles.last) <- names(concentrations.marbles.last[[1]])

# calculate cors
M.m <- cor(wide.marbles, method='spearman', use='pairwise.complete.obs')
Ml.m <- cor(wide.marbles.last, method='spearman', use='pairwise.complete.obs')

pdf(file.path(here(),'../Plots/MARBLES_metal_correlation.pdf'), height=14, width=14)
  corrplot(M.m, method='ellipse', type='upper', addCoef.col='black')
  corrplot(Ml.m, method='ellipse', type='upper', addCoef.col='black')
dev.off()

M.both <- M.e
M.both[lower.tri(M.both)] <- M.m[lower.tri(M.m)]

Ml.both <- Ml.e
Ml.both[lower.tri(Ml.both)] <- Ml.m[lower.tri(Ml.m)]

colnames(M.both) <- c('Sb', 'As', 'Ba', 'Cd', 'Cs', 'Cr', 'Co', 'Cu', 'Pb', 'Mn', 'Hg', 'Mo', 'Ni', 'Se', 'Tl', 'Sn', 'Zn')
rownames(M.both) <- c('Sb', 'As', 'Ba', 'Cd', 'Cs', 'Cr', 'Co', 'Cu', 'Pb', 'Mn', 'Hg', 'Mo', 'Ni', 'Se', 'Tl', 'Sn', 'Zn')
colnames(Ml.both) <- c('Sb', 'As', 'Ba', 'Cd', 'Cs', 'Cr', 'Co', 'Cu', 'Pb', 'Mn', 'Hg', 'Mo', 'Ni', 'Se', 'Tl', 'Sn', 'Zn')
rownames(Ml.both) <- c('Sb', 'As', 'Ba', 'Cd', 'Cs', 'Cr', 'Co', 'Cu', 'Pb', 'Mn', 'Hg', 'Mo', 'Ni', 'Se', 'Tl', 'Sn', 'Zn')

pdf(file.path(here(),'../Plots/metal_correlation_earlier.pdf'), height=7.5, width=7.5)
  corrplot(M.both, method='ellipse', addCoef.col='black', diag=F, tl.pos='d', tl.col='red', cl.pos='b', cl.cex=0.9, number.cex=0.8, tl.cex=1.0)
  text(18, 8, "EARLI", cex=2,  srt=90)
  text(0, 8, "MARBLES", cex=2, srt=90)
  text(9, -1.8, "Spearman Correlation", cex=1.1)
dev.off()

svg(file.path(here(),'../Plots/metal_correlation_later.svg'), height=14, width=14)
  corrplot(Ml.both, method='ellipse', addCoef.col='black', diag=F, tl.pos='d', tl.col='red', cl.pos='b', cl.cex=1.3, number.cex=1.4, tl.cex=2.5)
  text(18, 8, "EARLI", cex=3,  srt=90)
  text(0, 8, "MARBLES", cex=3, srt=90)
  text(9, -1.8, "Spearman Correlation", cex=1.6)
dev.off()
```

```{r marbles descriptives}
### sample characteristics
get_covarsM <- function(samps){
  samples <- samps[,c('SID','child_id','asd')] %>% unique()
  demo <- marbearl.metal[match(samples$child_id, marbearl.metal$coi_id),]
  demo <- demo[,c('famid','coi_id','coi_gender','maternal_age','GA_deliv_wks','MomEdu','BSRC_group_Sib_36mos')]
  demo$MomEdu <- demo$MomEdu > 3 #college degree or not
  
  marb.cov <- marb.cov[match(samples$child_id, marb.cov$COI_ID),]
  
  pd <- merge(samps, demo, by.x='child_id', by.y='coi_id', all.x=T)
  #pd$asd <- pd$BSRC_group_Sib_36mos
  pd$asd <- ifelse(pd$asd==0, 'TD',
                                 ifelse(pd$asd==1, 'ASD', 'Non-TD'))
  pd$asd <- factor(pd$asd, levels=c('TD','Non-TD','ASD'))
  pd$ever_smoke <- marb.cov[match(pd$child_id, marb.cov$COI_ID), 'MomEverSmoked']
  pd$reg_smoke <- marb.cov[match(pd$child_id, marb.cov$COI_ID), 'MomSmokedRegularly']
  pd$preg_smoke <- marb.cov[match(pd$child_id, marb.cov$COI_ID), 'SmokeYN_Pregnancy']
  pd$maternal_age <- marb.cov[match(pd$child_id, marb.cov$COI_ID),'MomAgeYr']

  pd
}

metals.marbles.T2 <- get_covarsM(metals.marbles)
metals.marbles.T3 <- get_covarsM(metals.marbles.last)

marbles.T2.summ <- metals.marbles.T2 %>% group_by(Analyte.Code) %>%
  summarize(mean=mean(Concentration_sg, na.rm=T), sd=sd(Concentration_sg, na.rm=T), median=median(Concentration_sg, na.rm=T),
            IQR=IQR(Concentration_sg, na.rm=T),  above_lod=sum(Concentration>LOD, na.rm=T), n=n(), pct_over_LOD=round(above_lod/n*100,1)) %>% data.frame()

marbles.T3.summ <- metals.marbles.T3 %>% group_by(Analyte.Code) %>%
  summarize(mean=mean(Concentration_sg, na.rm=T), sd=sd(Concentration_sg, na.rm=T), median=median(Concentration_sg, na.rm=T),
            IQR=IQR(Concentration_sg, na.rm=T),  above_lod=sum(Concentration>LOD, na.rm=T), n=n(), pct_over_LOD=round(above_lod/n*100,1)) %>% data.frame()


### metals descriptives tables

earli.T2.summ <- rounding(earli.T2.summ)
earli.T3.summ <- rounding(earli.T3.summ)
marbles.T2.summ <- rounding(marbles.T2.summ)
marbles.T3.summ <- rounding(marbles.T3.summ)

earli.T2.summ$cohort <- 'EARLI'
earli.T3.summ$cohort <- 'EARLI'
marbles.T2.summ$cohort <- 'MARBLES'
marbles.T3.summ$cohort <- 'MARBLES'

earli.T2.summ$Time <- 'earlier'
earli.T3.summ$Time <- 'later'
marbles.T2.summ$Time <- 'earlier'
marbles.T3.summ$Time <- 'later'


T2.summ <- rbind(earli.T2.summ[,c(1,9,2:8)], marbles.T2.summ[,c(1,9,2:8)])
T3.summ <- rbind(earli.T3.summ[,c(1,9,2:8)], marbles.T3.summ[,c(1,9,2:8)])

T2.summ <- T2.summ[order(T2.summ$Analyte.Code, T2.summ$cohort), ]
T3.summ <- T3.summ[order(T3.summ$Analyte.Code, T3.summ$cohort), ]

write.csv(T2.summ, file=file.path(here(),"../Tables/metal_descriptives_T2.csv"), row.names=F, quote=F)
write.csv(T3.summ, file=file.path(here(),"../Tables/metal_descriptives_T3.csv"), row.names=F, quote=F)


### count missing
metals.marbles.T2$period <- 'early'
metals.marbles.T3$period <- 'late'
count.mmiss <- rbind(metals.marbles.T2 %>% filter(Analyte.Code=="Lead"), metals.marbles.T3 %>% filter(Analyte.Code=="Lead"))
count.mmiss <- count.mmiss %>% filter(!duplicated(count.mmiss$child_id))
is.na(count.mmiss$asd) %>% sum() #11
is.na(count.mmiss$MomEdu) %>% sum() #0
is.na(count.mmiss$coi_gender) %>% sum() #0
is.na(count.mmiss$maternal_age) %>% sum() #0
is.na(count.mmiss$GA_collect) %>% sum() #0

count.mmiss <- count.mmiss %>% filter(!is.na(asd))
nrow(count.mmiss) #240

count.mmiss2 <- rbind(metals.marbles.T2 %>% filter(Analyte.Code=="Lead"), metals.marbles.T3 %>% filter(Analyte.Code=="Lead")) %>% filter(!is.na(asd))
table(count.mmiss2$period)
# early  late 
#   158   231 
nrow(count.mmiss2) #389


### descriptive tables
library(compareGroups)

gather_descriptivesM <- function(pd){
  descriptives <- pd[!duplicated(pd$child_id),]
  descriptives <- descriptives[!is.na(descriptives$asd) & !is.na(descriptives$MomEdu &
                              !is.na(descriptives$coi_gender)) & 
                              !is.na(descriptives$maternal_age) &
                              !is.na(descriptives$GA_deliv_wks),]
  
  descriptives$MomEdu <- ifelse(descriptives$MomEdu, "College Degree", "No College Degree")
  descriptives$ever_smoke <- ifelse(is.nan(descriptives$ever_smoke), NA, descriptives$ever_smoke)
  descriptives$ever_smoke <- factor(descriptives$ever_smoke)
  descriptives$preg_smoke <- ifelse(is.nan(descriptives$preg_smoke), NA, descriptives$preg_smoke)
  descriptives$preg_smoke <- factor(descriptives$preg_smoke)
  
  
  restab <- compareGroups(asd ~ MomEdu + coi_gender + maternal_age + GA_collect + GA_deliv_wks + ever_smoke + preg_smoke,
                          data=descriptives)
  createTable(restab)
}

tab.marbles.T2 <- gather_descriptivesM(metals.marbles.T2)
tab.marbles.T3 <- gather_descriptivesM(metals.marbles.T3)

export2word(tab.marbles.T2, file=file.path(here(),'../Tables/marbles_early_table1_20Sept2023.docx'))
export2word(tab.marbles.T3, file=file.path(here(),'../Tables/marbles_late_table1_20Sept2023.docx'))



IQRs.mT2 <- metals.marbles.T2 %>% group_by(Analyte.Code) %>% summarize(IQR(Concentration_sg, na.rm=T)) %>% data.frame
IQRs.mT2 <- IQRs.mT2[,2]
names(IQRs.mT2) <- metals

IQRs.mT3 <- metals.marbles.T3 %>% group_by(Analyte.Code) %>% summarize(IQR(Concentration_sg, na.rm=T)) %>% data.frame
IQRs.mT3 <- IQRs.mT3[,2]
names(IQRs.mT3) <- metals

# pdf(file.path(here(),"../Plots/marbles_firstT_violin.pdf"))
# for(met in metals){
#   p = ggplot(metals.marbles.T2 %>% filter(Analyte.Code==met), aes(x=factor(asd), y=Concentration_sg)) +
#   geom_violin() +geom_jitter(shape=16, position=position_jitter(0.2)) +
#     ggtitle(met) + xlab('') +
#     theme_bw() +
#     theme(axis.text.x=element_text(size=16, face="bold"))
#   print(p)
# }
# dev.off()
```

Early vs Late Preg

```{r}
#remove ones with only single time
table(metals.marbles.T2$SID %in% metals.marbles.T3$SID)
only.one <- metals.marbles.T2$SID[metals.marbles.T2$SID %in% metals.marbles.T3$SID] %>% unique()

marbles.epreg <- metals.marbles.T2[!metals.marbles.T2$SID %in% only.one, ]
marbles.lpreg <- metals.marbles.T3[!metals.marbles.T3$SID %in% only.one, ]

both <- intersect(unique(marbles.epreg$child_id), unique(marbles.lpreg$child_id))


marbles.epreg <- metals.marbles.T2[metals.marbles.T2$child_id %in% both, ]
marbles.lpreg <- metals.marbles.T3[metals.marbles.T3$child_id %in% both, ]

dim(marbles.epreg)
dim(marbles.lpreg)

#remove those with no sg
table(is.na(marbles.epreg$sg))
table(is.na(marbles.lpreg$sg))
no.sg <- marbles.epreg$child_id[is.na(marbles.epreg$sg)] %>% unique

marbles.epreg <- marbles.epreg[!marbles.epreg$child_id %in% no.sg, ]
marbles.lpreg <- marbles.lpreg[!marbles.lpreg$child_id %in% no.sg, ]
dim(marbles.epreg)
dim(marbles.lpreg)
identical(marbles.epreg$child_id, marbles.lpreg$child_id)

length(unique(marbles.epreg$child_id))
length(unique(marbles.lpreg$child_id))


cor(marbles.epreg$Concentration_imp_sg, marbles.lpreg$Concentration_imp_sg, use="complete.obs")
time.cors <- sapply(metals, FUN=function(x){
  epreg = marbles.epreg %>% filter(Analyte.Code==x)
  lpreg = marbles.lpreg %>% filter(Analyte.Code==x)
  raw_concentration = round(cor(epreg$Concentration, lpreg$Concentration, 
                                   use="complete.obs", method='spearman'),2)
  imputed_sg_corrected = round(cor(epreg$Concentration_imp_sg, lpreg$Concentration_imp_sg, 
                                   use="complete.obs", method='spearman'),2)
  data.frame(cbind(raw_concentration, imputed_sg_corrected))
})

write.csv(t(time.cors), file.path(here(),'../Tables/marbles_time_cor.csv'))

# plot.tm <- function(met){
#   m.epreg <- marbles.epreg %>% filter(Analyte.Code==met)
#   m.lpreg <- marbles.lpreg %>% filter(Analyte.Code==met)
#   dat <- data.frame(epreg=m.epreg$Concentration_imp_sg, lpreg=m.lpreg$Concentration_imp_sg)
#   
#   ggplot(dat, aes(epreg, lpreg)) +
#     geom_point() +
#     theme_bw() +
#     xlab(paste0("Early Pregnancy ",met)) +
#     ylab(paste0("Late Pregnancy ",met))
# }
# 
# plot.tm('Lead')
# plot.tm('Mercury')
# plot.tm('Cadmium')
```


# LOD descriptives

```{r}
metals.earli.T2$Time <- "early"
metals.earli.T3$Time <- "late"
metals.marbles.T2$Time <- "early"
metals.marbles.T3$Time <- "late"


urine.all <- rbind(metals.earli.T2 %>% select(Analyte.Code, Concentration, cohort, LOD, Time), 
                   metals.earli.T3 %>% select(Analyte.Code, Concentration, cohort, LOD, Time), 
                   metals.marbles.T2 %>% select(Analyte.Code, Concentration, cohort, LOD, Time), 
                   metals.marbles.T3 %>% select(Analyte.Code, Concentration, cohort, LOD, Time))

LOD.descriptives <- urine.all %>% group_by(Analyte.Code, cohort, LOD, Time) %>% 
  summarize(n=n(), above_lod=sum(Concentration>LOD), 
            pct_over_LOD=round(above_lod/n*100,1)) %>%
  data.frame()
head(LOD.descriptives)


write.csv(LOD.descriptives, file=file.path(here(),'../Tables/LOD_descriptives.csv'), row.names = F)

```


# Get rid of metals with too low proportion above LOD

```{r}
filter_metals <- function(met){
  met[!met$Analyte.Code %in% c("Beryllium","Platinum","Tungsten","Uranium","Vanadium"),]
}

metals.earli.T2 <- filter_metals(metals.earli.T2)
metals.earli.T3 <- filter_metals(metals.earli.T3)

metals.marbles.T2 <- filter_metals(metals.marbles.T2)
metals.marbles.T3 <- filter_metals(metals.marbles.T3)

mod.type <- mod.type[!names(mod.type) %in% c("Beryllium","Platinum","Tungsten","Uranium","Vanadium")]
metals <- metals[!names(metals) %in% c("Beryllium","Platinum","Tungsten","Uranium","Vanadium")]

```

# How many outliers

```{r outlier}
out.et2 <- metals.earli.T2 %>% group_by(Analyte.Code) %>% summarize(sum(is.na(Concentration_imp_sg_log2)))
out.et3 <- metals.earli.T3 %>% group_by(Analyte.Code) %>% summarize(sum(is.na(Concentration_imp_sg_log2)))
out.mt2 <- metals.marbles.T2 %>% group_by(Analyte.Code) %>% summarize(sum(is.na(Concentration_imp_sg_log2)))
out.mt3 <- metals.marbles.T3 %>% group_by(Analyte.Code) %>% summarize(sum(is.na(Concentration_imp_sg_log2)))

outliers <- cbind(out.et2, out.et3[,2], out.mt2[,2], out.mt3[,2])
colnames(outliers) <- c('Metal', 'EARLI_TE', 'EARLI_TL', 'MARBLES_TE', 'MARBLES_TL')
write.csv(outliers, file=file.path(here(),'../Tables/outliers_excluded_count.csv'), row.names=F, quote=F)
```

# Save

```{r models}
save(metals.earli.T2, metals.earli.T3, metals.marbles.T2, metals.marbles.T3, mod.type, metals, IQRs.T2, IQRs.T3, IQRs.mT2, IQRs.mT3, file=file.path(here(),'../Data/dat_for_models_20Sept2023.RDA'))
```