
source("utils.R")
source("config.R")
load("compute_SEMs.RData")

###############################################################################################
##### 1. Load data ############################################################################

print("Loading data...")
print(paste(cohort_name,"cohort"))
load(file_name)
i <- intersect(rownames(samples),colnames(betas))
b <- betas[,i]
s <- samples[i,]
check <- all.equal(rownames(s),colnames(b))
if(!check){stop("Rownames in sample table do not match with colnames in beta values matrix")}

###############################################################################################
##### 2. Remove X/Y chromosomes, cross reactive probes and probes containing SNPs #############

print("Removing X/Y chromosomes, cross reactive probes and probes containing SNPs...")
ann <- ann[which(ann$CHR %in% 1:22),]
i <- intersect(rownames(b),rownames(ann))
b <- b[i,]
ann <- ann[i,]
check <- all.equal(rownames(ann),rownames(b))
if(!check){stop("Rownames in annotation file do not match with rownames in beta values matrix")}

###############################################################################################
##### 3. Remove probes and samples for low call rate ##########################################

print("Probe and sample filtering for low call rate...")
probes_call_rate <- apply(b,1,function(x) sum(!is.na(x))/length(x))
b <- b[which(probes_call_rate > .95),]
samples_call_rate <- apply(b,2,function(x) sum(!is.na(x))/length(x))
s <- s[which(samples_call_rate > .95),]

i <- intersect(rownames(s),colnames(b))
b <- b[,i]
s <- s[i,]
ann <- ann[rownames(b),]
check1 <- all.equal(rownames(s),colnames(b))
check2 <- all.equal(rownames(b),rownames(ann))

if(!check1){stop("Rownames in sample table do not match with colnames in beta values matrix")}
if(!check2){stop("Rownames in annotation file do not match with rownames in beta values matrix")}

###############################################################################################
##### 4. Imputation ###########################################################################

ann <- ann[order(ann$CHR,ann$MAPINFO),]
b <- b[rownames(ann),]
for(i in 1:ncol(b)){
b[,i][which(b[,i] %in% c(0,1))] <- NA}

print("Imputation of missing data...")
cmd <- list()
for(i in 1:22){
cmd[[i]] <- b[which(ann$CHR == i),]}

no_cores <- detectCores() - 1
cl <- makeCluster(no_cores)
imputed_betas <- parLapply(cl,cmd,impute)
b_imp <- list.stack(imputed_betas)
rownames(b_imp) <- rownames(b)
b <- b_imp
rm(b_imp)

###############################################################################################
##### 5. test for bi-trimodality ##############################################################

print("Excluding probes with bi-trimodal distribution...")
peaks <- parRapply(cl,b,peaks)
npeaks <- unlist(lapply(peaks,length))
b <- b[which(npeaks==1),]
cmd <- b

###############################################################################################
##### 6. compute WBC ##########################################################################

print("Houseman WBC estimation...")
dnam <- t(b)
i <- intersect(colnames(dnam),rownames(model[["coefficients"]]))
dnam <- dnam[,i]
model[["coefficients"]] <- model[["coefficients"]][i,]
model[["df.residual"]] <- model[["df.residual"]][i]
model[["pvals"]] <- model[["pvals"]][i]
idx <- order(model$pvals)[1:500]
coefs <- coef(model)[idx,]

dnam <- dnam[,which(colnames(dnam) %in% rownames(coefs))]
coefs <- coefs[colnames(dnam),]
tmp <- re.match("^(.+)_(.+)$", rownames(dnam))
samples <- data.frame(row.names=rownames(dnam),chip=tmp[,2],chip.pos=tmp[,3])

for (i in 1:ncol(dnam)) {
    y <- dnam[,i]
    model <- lmer(y ~ (1 | chip) + (1 | chip.pos),
                  data=samples, REML=FALSE, na.action=na.exclude)
    dnam[,i] <- fixef(model)["(Intercept)"] + resid(model, type="response")}

types <- c("Monocytes","B","CD4T","NK","CD8T","Eosinophils","Neutrophils")
# Build projection matrix
L <- matrix(0, length(types), ncol(coefs),
            dimnames=list(types, colnames(coefs)))
L[,"(Intercept)"] <- 1
for (r in rownames(L)) {
    L[r,which(colnames(L) == r)] <- 1
}
L[c("Eosinophils", "Neutrophils"),"Granulocytes"] <- 1
L[c("Monocytes", "B", "CD4T", "NK", "CD8T"),"PBMC"] <- 1
# Project coefficients
coefs <- tcrossprod(coefs, L)
# Predict WBC differentials
wbc.predictions <- matrix(NA, nrow(dnam), ncol(coefs),
                          dimnames=list(rownames(dnam), colnames(coefs)))
A <- diag(ncol(coefs))
b <- rep(0, ncol(coefs))
for (i in 1:nrow(wbc.predictions)) {
    idx <- which(!is.na(dnam[i,]))
    D <- crossprod(coefs[idx,])
    d <- crossprod(coefs[idx,], dnam[i,idx])
    tmp <- try(solve.QP(D, d, A, b)$solution, silent=TRUE)
    if (!inherits(tmp, "try-error")) {
        wbc.predictions[i,] <- tmp
    }}

wbc.predictions <- as.data.frame(wbc.predictions)
wbc.predictions$Gran <- wbc.predictions$Eosinophils + wbc.predictions$Neutrophils
wbc.predictions <- wbc.predictions[,c(1:5,8)]

s <- cbind(s,wbc.predictions)
b <- cmd

###############################################################################################
##### 7. Epigenetic Age #######################################################################

print("Computing epigenetic age...")
## Horvath AA
datMethUsed <- as.matrix(b[which(rownames(b) %in% datClock$CpGmarker),])
datMethClock <- cbind(rownames(datMethUsed),datMethUsed)
colnames(datMethClock)[1] <- "ProbeID"
datMethClock <- data.frame(datMethClock)
for(i in 2:ncol(datMethClock)){datMethClock[,i] <- as.numeric(as.character(datMethClock[,i]))}
datMethClock[,1] <- as.character(datMethClock[,1])
nSamples <- ncol(datMethClock) - 1
datClock <- datClock[which(datClock$CpGmarker %in% c("(Intercept)",datMethClock[,1])),]
m <- match(datMethClock[,1],datClock$CpGmarker[2:nrow(datClock)])
datMethClock <- datMethClock[order(m),]
all.equal(datMethClock[,1],datClock$CpGmarker[2:nrow(datClock)])
datMethClock <- as.matrix(datMethClock[,-1])

HorvathAge <- rep(NA,length=nSamples)
for(i in 1:nSamples){
x <- datMethClock[,i]
rm <- names(x[which(is.na(x))])
if(length(rm)>0){x <- x[which(!is.na(x))]; datClocknew <- datClock[which(!datClock$CpGmarker %in% rm),]} else{datClocknew <- datClock}
HorvathAge[i]=as.numeric(anti.trafo(datClocknew$CoefficientTraining[1]+ t(x) %*% as.numeric(datClocknew$CoefficientTraining[-1])))}

HorvathAA <- HorvathAge - s$age
HorvathEEAA <- resid(glm(HorvathAA ~ s$age))
HorvathIEAA <- resid(glm(HorvathAA ~ s$age+s$Monocytes+s$B+s$CD4T+s$CD8T+s$NK+s$Gran))

## Hannum AA
datMethUsed <- as.matrix(b[which(rownames(b) %in% dat0$Name),])
datMethClock <- cbind(rownames(datMethUsed),datMethUsed)
colnames(datMethClock)[1] <- "ProbeID"
datMethClock <- data.frame(datMethClock)
for(i in 2:ncol(datMethClock)){datMethClock[,i] <- as.numeric(as.character(datMethClock[,i]))}
datMethClock[,1] <- as.character(datMethClock[,1])
nSamples <- ncol(datMethClock) - 1
datClock <- dat0[which(dat0$Name %in% datMethClock[,1]),]
datClock <- rbind(c("(Intercept)",rep(NA,length=5),0.695507258),datClock)
datClock$CoefficientHannum <- as.numeric(datClock$CoefficientHannum)
m <- match(datMethClock[,1],datClock$Name[2:nrow(datClock)])
datMethClock <- datMethClock[order(m),]
all.equal(datMethClock[,1],datClock$Name[2:nrow(datClock)])
rownames(datMethClock) <- datMethClock[,1]
datMethClock <- as.matrix(datMethClock[,-1])

HannumAge <- rep(NA,length=nSamples)
for(i in 1:nSamples){
x <- datMethClock[,i]
rm <- names(x[which(is.na(x))])
if(length(rm)>0){x <- x[which(!is.na(x))]; datClocknew <- datClock[which(!datClock$Name %in% rm),]} else{datClocknew <- datClock}
HannumAge[i]=as.numeric(7.5 + t(x) %*% as.numeric(datClocknew$CoefficientHannum[-1]))}

HannumAA <- HannumAge - s$age
HannumEEAA <- resid(glm(HannumAA ~ s$age))
HannumIEAA <- resid(glm(HannumAA ~ s$age+s$Monocytes+s$B+s$CD4T+s$CD8T+s$NK+s$Gran))

## Levine
datMethUsed <- as.matrix(b[which(rownames(b) %in% lev$CpG),])
datMethClock <- cbind(rownames(datMethUsed),datMethUsed)
colnames(datMethClock)[1] <- "ProbeID"
datMethClock <- data.frame(datMethClock)
for(i in 2:ncol(datMethClock)){datMethClock[,i] <- as.numeric(as.character(datMethClock[,i]))}
datMethClock[,1] <- as.character(datMethClock[,1])
nSamples <- ncol(datMethClock) - 1
datClock <- lev[which(lev$CpG %in% datMethClock[,1]),]
datClock <- rbind(c("(Intercept)",60.664),datClock)
datClock$Weight <- as.numeric(datClock$Weight)
m <- match(datMethClock[,1],datClock$CpG[2:nrow(datClock)])
datMethClock <- datMethClock[order(m),]
all.equal(datMethClock[,1],datClock$CpG[2:nrow(datClock)])
datMethClock <- as.matrix(datMethClock[,-1])

LevineAge <- rep(NA,length=nSamples)
for(i in 1:nSamples){
x <- datMethClock[,i]
rm <- names(x[which(is.na(x))])
if(length(rm)>0){x <- x[which(!is.na(x))]; datClocknew <- datClock[which(!datClock$CpG %in% rm),]} else{datClocknew <- datClock}
#LevineAge[i]=as.numeric(anti.trafo(datClocknew$Weight[1]+ t(x) %*% as.numeric(datClocknew$Weight[-1])))}
LevineAge[i]=as.numeric(datClocknew$Weight[1]+ t(x) %*% as.numeric(datClocknew$Weight[-1]))}

LevineAA <- LevineAge - s$age
LevineEEAA <- resid(glm(LevineAA ~ s$age))
LevineIEAA <- resid(glm(LevineAA ~ s$age+s$Monocytes+s$B+s$CD4T+s$CD8T+s$NK+s$Gran))

### add epiclock to sample file
s$HorvathEEAA <- HorvathEEAA
s$HorvathIEAA <- HorvathIEAA
s$HannumEEAA <- HannumEEAA
s$HannumIEAA <- HannumIEAA
s$LevineEEAA <- LevineEEAA
s$LevineIEAA <- LevineIEAA

###############################################################################################
##### 8. Compute SEMs #########################################################################

print("Computing SEMs...")

bSEMs <- parRapply(cl,b,identify_SEMs)
bSEMs <- t(matrix(bSEMs,ncol(b),nrow(b)))
rownames(bSEMs) <- rownames(b)

sems <- parCapply(cl,bSEMs,Sum_abs)
log_sems <- log(sems+1)
sems_hyper <- parCapply(cl,bSEMs,Sum_hyper)
log_sems_hyper <- log(sems_hyper+1)
sems_hypo <- parCapply(cl,bSEMs,Sum_hypo)
log_sems_hypo <- log(sems_hypo+1)

s$log_sems <- log_sems
s$log_sems_hyper <- log_sems_hyper
s$log_sems_hypo <- log_sems_hypo

sems_per_probe <- parRapply(cl,bSEMs,Sum_abs)
names(sems_per_probe) <- rownames(b)
sems_per_probe_hyper <- parRapply(cl,bSEMs,Sum_hyper)
names(sems_per_probe_hyper) <- rownames(b)
sems_per_probe_hypo <- parRapply(cl,bSEMs,Sum_hypo)
names(sems_per_probe_hypo) <- rownames(b)

s <- cbind(s,samples)

###############################################################################################
##### 9. Compute SEMs WBC Residuals ###########################################################

wbc_b <- parRapply(cl,b,resWBC,data=s)
wbc_b <- t(matrix(wbc_b,ncol(b),nrow(b)))

bSEMs <- parRapply(cl,wbc_b,identify_SEMs)
bSEMs <- t(matrix(bSEMs,ncol(wbc_b),nrow(wbc_b)))
rownames(bSEMs) <- rownames(b)

sems_wbc <- parCapply(cl,bSEMs,Sum_abs)
log_sems_wbc <- log(sems_wbc+1)
s$log_sems_wbc <- log_sems_wbc

###############################################################################################
##### 10. technical variables residuals #######################################################

print("Computing technical variables residuals...")
#tech_b <- parRapply(cl,wbc_b,resTechnical,data=s)
#tech_b <- t(matrix(tech_b,ncol(wbc_b),nrow(wbc_b)))

tech_b <- NULL
for(i in 1:108){
start <- (i-1)*5431 + 1
end <- start + 5430
print(paste("step",i,"start",start,"end",end))
tmp <- parRapply(cl,wbc_b[start:end,],resTechnical,data=s)
tmp <- t(matrix(tmp,ncol(wbc_b),1000))
tech_b <- rbind(tech_b,tmp)}

bSEMs_tech <- parRapply(cl,tech_b,identify_SEMs)
bSEMs_tech <- t(matrix(bSEMs_tech,ncol(tech_b),nrow(tech_b)))

sems <- parCapply(cl,bSEMs_tech,Sum_abs)
log_sems <- log(sems+1)
sems_hyper <- parCapply(cl,bSEMs_tech,Sum_hyper)
log_sems_hyper <- log(sems_hyper+1)
sems_hypo <- parCapply(cl,bSEMs_tech,Sum_hypo)
log_sems_hypo <- log(sems_hypo+1)

s$log_sems_tech <- log_sems
s$log_sems_hyper_tech <- log_sems_hyper
s$log_sems_hypo_tech <- log_sems_hypo

###############################################################################################
stopCluster(cl)
print("Saving preliminary results...")

###############################################################################################
##### 11. Descriptives ########################################################################

covar <- summary(s[,c("age","sex","bmi","smoking","alcohol","education","phys.activ")])
wbc <- summary(s[,c("Monocytes","B","CD4T","NK","CD8T","Gran")])
AA <- summary(s[,c("HorvathEEAA","HorvathIEAA","HannumEEAA","HannumIEAA","LevineEEAA","LevineIEAA")])
SEMs <- summary(s[,c("log_sems","log_sems_hyper","log_sems_hypo","log_sems_wbc",
"sems_ezh2","sems_ezh2_hyper","sems_ezh2_hypo","log_sems_tech",
"log_sems_hyper_tech","log_sems_hypo_tech")])

descr_bmi <- tapply(s$log_sems,s$bmi,summary_sd)
descr_smo <- tapply(s$log_sems,s$smoking,summary_sd)
descr_alc <- tapply(s$log_sems,s$alcohol,summary_sd)
descr_phy <- tapply(s$log_sems,s$phys.activ,summary_sd)
descr_edu <- tapply(s$log_sems,s$education,summary_sd)
descr_age_sems_aa <- round(cor(cbind(s[c("age","Monocytes","B","CD4T","NK","CD8T","Gran",
"HorvathEEAA","HorvathIEAA","HannumEEAA","HannumIEAA","LevineEEAA","LevineIEAA",
"log_sems","log_sems_hyper","log_sems_hypo","log_sems_wbc",
"sems_ezh2","sems_ezh2_hyper","sems_ezh2_hypo","log_sems_tech",
"log_sems_hyper_tech","log_sems_hypo_tech")])),digits=2)


###############################################################################################
##### 12. Regression analyses #################################################################

#### s$center <- s$sample.type
#### In the EPIC cohort 'center' indicate the center of recruitment.
#### Please remove the variable if it is not relevant for your cohort.
#### Add any other covariate if it is relevant for your cohort.

s$y <- s$log_sems

fit0 <- glm(y~age+sex,data=s)
fit0.smo.add <- glm(y~age+smoking+sex,data=s)
fit0.smo.int <- glm(y~age*smoking+sex,data=s)
fit0.bmi.add <- glm(y~age+bmi+sex,data=s)
fit0.bmi.int <- glm(y~age*bmi+sex,data=s)
fit0.alc.add <- glm(y~age+alcohol+sex,data=s)
fit0.alc.int <- glm(y~age*alcohol,data=s)
fit0.phy.add <- glm(y~age+phys.activ+sex,data=s)
fit0.phy.int <- glm(y~age*phys.activ,data=s)
fit0.edu.add <- glm(y~age+education+sex,data=s)
fit0.edu.int <- glm(y~age*education+sex,data=s)

fit1 <- glm(y~age+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.smo.int <- glm(y~age*smoking+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.bmi.int <- glm(y~age*bmi+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.alc.int <- glm(y~age*alcohol+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.phy.int <- glm(y~age*phys.activ+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.edu.int <- glm(y~age*education+sex+bmi+smoking+alcohol+phys.activ+education,data=s)

res_log_sems <- list(fit0,fit0.smo.add,fit0.smo.int,fit0.bmi.add,fit0.bmi.int,
     fit0.alc.add,fit0.alc.int,fit0.edu.add,fit0.edu.int,fit0.phy.add,fit0.phy.int,
     fit1,fit1.smo.int,fit1.bmi.int,fit1.alc.int,fit1.edu.int,fit1.phy.int)

s$y <- s$HorvathEEAA

fit0 <- glm(y~age+sex,data=s)
fit0.smo.add <- glm(y~age+smoking+sex,data=s)
fit0.smo.int <- glm(y~age*smoking+sex,data=s)
fit0.bmi.add <- glm(y~age+bmi+sex,data=s)
fit0.bmi.int <- glm(y~age*bmi+sex,data=s)
fit0.alc.add <- glm(y~age+alcohol+sex,data=s)
fit0.alc.int <- glm(y~age*alcohol,data=s)
fit0.phy.add <- glm(y~age+phys.activ+sex,data=s)
fit0.phy.int <- glm(y~age*phys.activ,data=s)
fit0.edu.add <- glm(y~age+education+sex,data=s)
fit0.edu.int <- glm(y~age*education+sex,data=s)

fit1 <- glm(y~age+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.smo.int <- glm(y~age*smoking+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.bmi.int <- glm(y~age*bmi+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.alc.int <- glm(y~age*alcohol+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.phy.int <- glm(y~age*phys.activ+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.edu.int <- glm(y~age*education+sex+bmi+smoking+alcohol+phys.activ+education,data=s)

res_HorvathEEAA <- list(fit0,fit0.smo.add,fit0.smo.int,fit0.bmi.add,fit0.bmi.int,
     fit0.alc.add,fit0.alc.int,fit0.edu.add,fit0.edu.int,fit0.phy.add,fit0.phy.int,
     fit1,fit1.smo.int,fit1.bmi.int,fit1.alc.int,fit1.edu.int,fit1.phy.int)

###############################################################################################
###############################################################################################

s$y <- s$HorvathIEAA

fit0 <- glm(y~age+sex,data=s)
fit0.smo.add <- glm(y~age+smoking+sex,data=s)
fit0.smo.int <- glm(y~age*smoking+sex,data=s)
fit0.bmi.add <- glm(y~age+bmi+sex,data=s)
fit0.bmi.int <- glm(y~age*bmi+sex,data=s)
fit0.alc.add <- glm(y~age+alcohol+sex,data=s)
fit0.alc.int <- glm(y~age*alcohol,data=s)
fit0.phy.add <- glm(y~age+phys.activ+sex,data=s)
fit0.phy.int <- glm(y~age*phys.activ,data=s)
fit0.edu.add <- glm(y~age+education+sex,data=s)
fit0.edu.int <- glm(y~age*education+sex,data=s)

fit1 <- glm(y~age+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.smo.int <- glm(y~age*smoking+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.bmi.int <- glm(y~age*bmi+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.alc.int <- glm(y~age*alcohol+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.phy.int <- glm(y~age*phys.activ+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.edu.int <- glm(y~age*education+sex+bmi+smoking+alcohol+phys.activ+education,data=s)

res_HorvathIEAA <- list(fit0,fit0.smo.add,fit0.smo.int,fit0.bmi.add,fit0.bmi.int,
     fit0.alc.add,fit0.alc.int,fit0.edu.add,fit0.edu.int,fit0.phy.add,fit0.phy.int,
     fit1,fit1.smo.int,fit1.bmi.int,fit1.alc.int,fit1.edu.int,fit1.phy.int)

###############################################################################################
###############################################################################################

s$y <- s$HannumEEAA

fit0 <- glm(y~age+sex,data=s)
fit0.smo.add <- glm(y~age+smoking+sex,data=s)
fit0.smo.int <- glm(y~age*smoking+sex,data=s)
fit0.bmi.add <- glm(y~age+bmi+sex,data=s)
fit0.bmi.int <- glm(y~age*bmi+sex,data=s)
fit0.alc.add <- glm(y~age+alcohol+sex,data=s)
fit0.alc.int <- glm(y~age*alcohol,data=s)
fit0.phy.add <- glm(y~age+phys.activ+sex,data=s)
fit0.phy.int <- glm(y~age*phys.activ,data=s)
fit0.edu.add <- glm(y~age+education+sex,data=s)
fit0.edu.int <- glm(y~age*education+sex,data=s)

fit1 <- glm(y~age+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.smo.int <- glm(y~age*smoking+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.bmi.int <- glm(y~age*bmi+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.alc.int <- glm(y~age*alcohol+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.phy.int <- glm(y~age*phys.activ+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.edu.int <- glm(y~age*education+sex+bmi+smoking+alcohol+phys.activ+education,data=s)

res_HannumEEAA <- list(fit0,fit0.smo.add,fit0.smo.int,fit0.bmi.add,fit0.bmi.int,
     fit0.alc.add,fit0.alc.int,fit0.edu.add,fit0.edu.int,fit0.phy.add,fit0.phy.int,
     fit1,fit1.smo.int,fit1.bmi.int,fit1.alc.int,fit1.edu.int,fit1.phy.int)

###############################################################################################
###############################################################################################

s$y <- s$HannumIEAA

fit0 <- glm(y~age+sex,data=s)
fit0.smo.add <- glm(y~age+smoking+sex,data=s)
fit0.smo.int <- glm(y~age*smoking+sex,data=s)
fit0.bmi.add <- glm(y~age+bmi+sex,data=s)
fit0.bmi.int <- glm(y~age*bmi+sex,data=s)
fit0.alc.add <- glm(y~age+alcohol+sex,data=s)
fit0.alc.int <- glm(y~age*alcohol,data=s)
fit0.phy.add <- glm(y~age+phys.activ+sex,data=s)
fit0.phy.int <- glm(y~age*phys.activ,data=s)
fit0.edu.add <- glm(y~age+education+sex,data=s)
fit0.edu.int <- glm(y~age*education+sex,data=s)

fit1 <- glm(y~age+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.smo.int <- glm(y~age*smoking+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.bmi.int <- glm(y~age*bmi+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.alc.int <- glm(y~age*alcohol+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.phy.int <- glm(y~age*phys.activ+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.edu.int <- glm(y~age*education+sex+bmi+smoking+alcohol+phys.activ+education,data=s)

res_HannumIEAA <- list(fit0,fit0.smo.add,fit0.smo.int,fit0.bmi.add,fit0.bmi.int,
     fit0.alc.add,fit0.alc.int,fit0.edu.add,fit0.edu.int,fit0.phy.add,fit0.phy.int,
     fit1,fit1.smo.int,fit1.bmi.int,fit1.alc.int,fit1.edu.int,fit1.phy.int)

###############################################################################################
###############################################################################################

s$y <- s$LevineEEAA

fit0 <- glm(y~age+sex,data=s)
fit0.smo.add <- glm(y~age+smoking+sex,data=s)
fit0.smo.int <- glm(y~age*smoking+sex,data=s)
fit0.bmi.add <- glm(y~age+bmi+sex,data=s)
fit0.bmi.int <- glm(y~age*bmi+sex,data=s)
fit0.alc.add <- glm(y~age+alcohol+sex,data=s)
fit0.alc.int <- glm(y~age*alcohol,data=s)
fit0.phy.add <- glm(y~age+phys.activ+sex,data=s)
fit0.phy.int <- glm(y~age*phys.activ,data=s)
fit0.edu.add <- glm(y~age+education+sex,data=s)
fit0.edu.int <- glm(y~age*education+sex,data=s)

fit1 <- glm(y~age+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.smo.int <- glm(y~age*smoking+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.bmi.int <- glm(y~age*bmi+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.alc.int <- glm(y~age*alcohol+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.phy.int <- glm(y~age*phys.activ+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.edu.int <- glm(y~age*education+sex+bmi+smoking+alcohol+phys.activ+education,data=s)

res_LevineEEAA <- list(fit0,fit0.smo.add,fit0.smo.int,fit0.bmi.add,fit0.bmi.int,
     fit0.alc.add,fit0.alc.int,fit0.edu.add,fit0.edu.int,fit0.phy.add,fit0.phy.int,
     fit1,fit1.smo.int,fit1.bmi.int,fit1.alc.int,fit1.edu.int,fit1.phy.int)

###############################################################################################
###############################################################################################

s$y <- s$LevineIEAA

fit0 <- glm(y~age+sex,data=s)
fit0.smo.add <- glm(y~age+smoking+sex,data=s)
fit0.smo.int <- glm(y~age*smoking+sex,data=s)
fit0.bmi.add <- glm(y~age+bmi+sex,data=s)
fit0.bmi.int <- glm(y~age*bmi+sex,data=s)
fit0.alc.add <- glm(y~age+alcohol+sex,data=s)
fit0.alc.int <- glm(y~age*alcohol,data=s)
fit0.phy.add <- glm(y~age+phys.activ+sex,data=s)
fit0.phy.int <- glm(y~age*phys.activ,data=s)
fit0.edu.add <- glm(y~age+education+sex,data=s)
fit0.edu.int <- glm(y~age*education+sex,data=s)

fit1 <- glm(y~age+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.smo.int <- glm(y~age*smoking+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.bmi.int <- glm(y~age*bmi+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.alc.int <- glm(y~age*alcohol+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.phy.int <- glm(y~age*phys.activ+sex+bmi+smoking+alcohol+phys.activ+education,data=s)
fit1.edu.int <- glm(y~age*education+sex+bmi+smoking+alcohol+phys.activ+education,data=s)

res_LevineIEAA <- list(fit0,fit0.smo.add,fit0.smo.int,fit0.bmi.add,fit0.bmi.int,
     fit0.alc.add,fit0.alc.int,fit0.edu.add,fit0.edu.int,fit0.phy.add,fit0.phy.int,
     fit1,fit1.smo.int,fit1.bmi.int,fit1.alc.int,fit1.edu.int,fit1.phy.int)

###############################################################################################
###############################################################################################


save(covar,wbc,AA,SEMs,descr_bmi,descr_smo,descr_alc,descr_phy,descr_edu,descr_age_sems_aa,
     res_log_sems,res_log_sems_wbc,res_log_sems_hyper,res_log_sems_hypo,res_sems_ezh2,
	 res_HorvathIEAA,res_HorvathEEAA,res_HannumIEAA,res_HannumEEAA,res_LevineIEAA,res_LevineEEAA,
	 sems_per_probe,sems_per_probe_hyper,sems_per_probe_hypo,
	 file=paste0("output/",cohort_name,"_Results.rda"))




