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"))