#G*E model with quadratic effects for A, C and E #No OSDZ pairs rm(list = ls(all = TRUE)) # Load Library require(OpenMx) require(psych) require(foreign) setwd("D:/Odense/GEmodels") indata <- read.table ("BMIdata.dat", header=T, na=-99,dec = ".") indata2 <- subset(indata, age1 !="-99" & age2 !="-99" ) #create dataset without missing values for the "definition variable" indata2$age1 <- indata2$age1-50 indata2$age2 <- indata2$age2-50 head(indata2) # Select Variables for Analysis Vars <- c('bmi1','bmi2','age1','age2') mzData <- subset(indata2, zyg==1 & sp1==1 & sp2==1, select=Vars) dzData <- subset(indata2, zyg==2 & sp1==1 & sp2==1, select=Vars) #head(mzData) # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, labels="a11", values=30, lbound=-10, name="a" ) pathC <- mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, labels="c11", values=10, lbound=-10, name="c" ) pathE <- mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, labels="e11", values=20, lbound=-10, name="e" ) # GE-interaction part of path coefficients pathAI <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0.1, label="aI11", name="aI" ) pathCI <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0, label="cI11", name="cI" ) pathEI <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0.1, label="eI11", name="eI" ) # Matrices generated to hold moderated A, C, and E computed Variance Components covAI <- mxAlgebra( expression=(a+ age1%*%aI) %*% t(a+ age2%*%aI), name="AI" ) covCI <- mxAlgebra( expression=(c+ age1%*%cI) %*% t(c+ age2%*%cI), name="CI" ) covEI <- mxAlgebra( expression=(e+ age1%*%eI) %*% t(e+ age2%*%eI), name="EI" ) # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins meanGMZ <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=c(24, 24), labels="bmi_MeanMZ", name="MeanMZ" ) meanGDZ <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=c(24, 24), labels="bmi_MeanDZ", name="MeanDZ" ) # Matrix for moderating/interacting variable (linear age effect and quadratic age effect) def1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1"), name="age1") def2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age2"), name="age2") # Matrices declared to store linear and quadratic Coefficients for covariate pathB <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1, label="b11mz", name="bz" ) mean1MZ <- mxAlgebra( expression= bz%*%age1, name="Age1RMZ") mean2MZ <- mxAlgebra( expression= bz%*%age2, name="Age2RMZ") mean1DZ <- mxAlgebra( expression= bz%*%age1, name="Age1RDZ") mean2DZ <- mxAlgebra( expression= bz%*%age2, name="Age2RDZ") #different means for MZ and DZ twins meanGIMZ <- mxAlgebra( expression= cbind((MeanMZ + Age1RMZ),(MeanMZ + Age2RMZ)), name="expMeanMZ") meanGIDZ <- mxAlgebra( expression= cbind((MeanDZ + Age1RDZ),(MeanDZ + Age2RDZ)), name="expMeanDZ") covMZ <- mxAlgebra( expression= rbind( cbind(AI+CI+EI, AI+CI), cbind(AI+CI, AI+CI+EI)), name="expCovMZ") covDZ <- mxAlgebra( expression= rbind( cbind(AI+CI+EI, 0.5%x%AI+CI), cbind(0.5%x%AI+CI ,AI+CI+EI)), name="expCovDZ" ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=c("bmi1","bmi2")) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=c("bmi1","bmi2") ) funML <- mxFitFunctionML() # Combine Groups pars <- list( pathA, pathC, pathE, pathAI, pathCI, pathEI) pars2 <- list(covAI, covCI, covEI, def1, def2) modelMZ <- mxModel( "MZ", pars, pars2, meanGMZ, meanGIMZ, mean1MZ, mean2MZ, pathB, covMZ, dataMZ, expMZ, funML) modelDZ <- mxModel( "DZ", pars, pars2, meanGDZ, meanGIDZ, mean1DZ, mean2DZ, pathB, covDZ, dataDZ, expDZ, funML) minus2ll <- mxAlgebra( MZ.objective+ DZ.objective, name="minus2sumloglikelihood" ) obj <- mxFitFunctionAlgebra( "minus2sumloglikelihood" ) ci <- mxCI(c("aI11","cI11","eI11","b11mz")) model <- mxModel( "Model", pars, modelMZ, modelDZ, minus2ll, obj, ci ) modelFit <- mxTryHard( model, scale=0.05, intervals=F) modelSummary <- summary(modelFit) modelSummary #Model without C effect #modelNoC <- omxSetParameters(model, label="cI11", free=F, values=0) #modelNoC <- omxSetParameters(modelNoC, label="c11", free=F, values=0) #modelFitNoC <- mxTryHard( modelNoC, scale=0.05, intervals=F) #modelSummary <- summary(modelFitNoC) #modelSummary #mxCompare(modelFit, modelFitNoC)