###### Residulals for lipids

Residuals=function(x)rstandard(lm(x~Age+Sex+SmokingStatus+Experiment_batch+Education+BMI+LipidMedication+ApoE4carrierStatus,data=data))

######## Inverse normal transformation using R package library("RNOmni")
y=rankNorm(Residuals)

########## Heritability Estimation

OpenMx scripts were adapted from the OpenMx website https://protect-au.mimecast.com/s/W0C9CYW8ngs0zQ5sjAlaq?domain=openmx.ssri.psu.edu

############## Association analysis
## Association of Gene expression with lipids using twin data
library(nlme)
library("rcompanion")
library("caret")
lme.out=lme(Lipid~GeneExpression,random=~1|FamilyID,na.action=na.exclude,method="ML",data=data)  ## Lipids and Gene expressions - preprocessed

## To estimate the variance explained by the significantly assiated probes we used the following function.

r2=nagelkerke(lme.out) 

### For lipids associated with multiple gene expression probes we first obtained a reduced list of probes using penalised regression
model_glm <- caret::train(y ~.,
                         data = train_data,
                         method = "glmnet",
                         preProcess = c ("scale", "center"),
                         trControl = trainControl(method = "repeatedcv", 
                                                  number = 3, 
                                                  repeats = 5, 
                                                  savePredictions = FALSE, 
                                                  verboseIter = FALSE))

### The gene expression probes in from the final model were selected and refit using lme as in above to estimate the variance explained.

######## Association of methylation vs lipids

rxy=cor(lipid,methylation) ## correlation between lipid and methylation

lme.out=lme(lipid~methylation,random=~1|FamilyID,na.action=na.exclude,method="ML",data=dat)

######## Association of methylation vs gene expression

lme.out=lme(gene-expression~methylation,random=~1|FamilyID,na.action=na.exclude,method="ML",data=dat)




