# ---------------------------------------------------------------------------------------------------------------------- # Script: oneACEvc7.R # Author: Sarah Medland - Hermine Maes # Date: 29 02 2020 - 05 22 2021 # # Twin Univariate ACE model to estimate causes of variation across multiple groups, including siblings # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # Load Libraries & Options rm(list=ls()) getwd() ## YOU WILL NEED TO CHANGE THIS TO YOUR DIRECTORY setwd("/home/hmaes/2021/day2/practical2") library(OpenMx) library(psych) source("miFunctions.R") mxOption(NULL,"Default optimizer","NPSOL") # if you do not have access to NPSOL-enabled version of OpenMx, uncomment the statement below: #mxOption(NULL,"Default optimizer","SLSQP") # Create Output filename <- "oneACEvc7" sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE) # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA # Load Data tsDataS <- read.table("tsDataS7.txt", header=T) dim(tsDataS) describe(tsDataS, skew=F) # Select Variables for Analysis nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- c('Twin1','Twin2') # list of variables names of observed variables covVars <- c('age1','age2','sex1','sex2') # list of variables names of covariates # Select Data for Analysis mzData <- subset(tsDataS, zyg==1, c(selVars,covVars)) dzData <- subset(tsDataS, zyg==2, c(selVars,covVars)) cov(mzData[,selVars],use="complete") cov(dzData[,selVars],use="complete") # Set Starting Values svBe <- .01 # start value for regressions svMu <- 0 # start value for means svVa <- .2 # start value for path coefficient svVe <- .5 # start value for path coefficient for e # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL 1: univariate twin model + covariates # see oneACEvc_1cov.R for single model script # ---------------------------------------------------------------------------------------------------------------------- # Create Matrices for Covariates and linear Regression Coefficients betaS <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, labels="betaS", name="bS" ) betaA <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, labels="betaA", name="bA" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex" ) defAge <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, labels=c("data.age1","data.age2"), name="Age" ) # Create Algebra for expected Mean Matrices intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMu, labels="interC", name="intercept" ) expMean <- mxAlgebra( expression = intercept + bS*Sex + bA*Age, name="expMean" ) # Create Matrices for Variance Components covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svVa, label="VA11", name="VA" ) covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svVa, label="VC11", name="VC" ) covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svVe, label="VE11", name="VE" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= VA+VC+VE, name="V" ) covMZ <- mxAlgebra( expression= VA+VC, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%VA+ VC, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups defs <- list( defAge, defSex) pars <- list( intercept, betaS, betaA, covA, covC, covE, covP ) modelMZ <- mxModel( pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Unstandardized & Standardized Variance Components rowUS <- rep('US',nv) colUS <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv) estUS <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="US", dimnames=list(rowUS,colUS) ) # Create Confidence Interval Objects ciACE <- mxCI( "US[1,1:6]" ) # Build Model with Confidence Intervals modelACE1 <- mxModel( "oneACEvc_1cov", pars, modelMZ, modelDZ, multi, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL 1 # Run ACE Model fitACE1 <- mxRun( modelACE1, intervals=T ) sumACE1 <- summary( fitACE1 ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE1) fitEstCIs(fitACE1,colUS) # this will give an error if you set intervals=F, but still print the parameter estimates # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL 2: univariate twin model + covariates + ADD SIBLING # see oneACEvc_2sib.R for single model script # ---------------------------------------------------------------------------------------------------------------------- # Select Variables for Analysis nv <- 1 # number of variables ntv <- nv*3 # number of total variables selVars <- c('Twin1','Twin2','Sib') # list of variables names of observed variables covVars <- c('age1','age2','age3','sex1','sex2','sex3') # list of variables names of covariates # Select Data for Analysis mzData <- subset(tsDataS, zyg==1, c(selVars,covVars)) dzData <- subset(tsDataS, zyg==2, c(selVars,covVars)) # Create Matrices for Covariates and linear Regression Coefficients defSex <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, labels=c("data.sex1","data.sex2","data.sex3"), name="Sex" ) defAge <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, labels=c("data.age1","data.age2","data.age3"), name="Age" ) # Create Algebra for expected Mean Matrices intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMu, labels="interC", name="intercept" ) expMean <- mxAlgebra( expression = intercept + bS*Sex + bA*Age, name="expMean" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins & extra sibling expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ, cDZ), cbind(t(cMZ), V, cDZ), cbind(t(cDZ), cDZ, V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ, cDZ), cbind(t(cDZ), V, cDZ), cbind(t(cDZ), cDZ, V)), name="expCovDZ" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) # Create Model Objects for Multiple Groups defs <- list( defAge, defSex) pars <- list( intercept, betaS, betaA, covA, covC, covE, covP, covDZ ) modelMZ <- mxModel( pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Build Model with Confidence Intervals modelACE2 <- mxModel( "oneACEvc_2sib", pars, modelMZ, modelDZ, multi, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL 2 # Run ACE Model fitACE2 <- mxRun( modelACE2, intervals=T ) sumACE2 <- summary( fitACE2 ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE2) fitEstCIs(fitACE2,colUS) # this will give an error if you set intervals=F, but still print the parameter estimates # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL 3: univariate twin model + covariates + SIBLING - ALTERNATE PARAMETERIZATION # see oneACEvc_3alt.R for single model script # ---------------------------------------------------------------------------------------------------------------------- # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins & extra sibling relAmz <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=FALSE, values=c(1, 1,.5,1,.5,1), name="rAmz" ) relAdz <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=FALSE, values=c(1,.5,.5,1,.5,1), name="rAdz" ) relC <- mxMatrix( type="Unit", nrow=ntv, ncol=ntv, free=FALSE, name="rC" ) relE <- mxMatrix( type="Iden", nrow=ntv, ncol=ntv, free=FALSE, name="rE" ) expCovMZ <- mxAlgebra( expression= VA%x%rAmz + VC%x%rC + VE%x%rE, name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= VA%x%rAdz + VC%x%rC + VE%x%rE, name="expCovDZ" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups defs <- list( defAge, defSex) pars <- list( intercept, betaS, betaA, covA, covC, covE, covP, relAmz, relAdz, relC, relE ) modelMZ <- mxModel( pars, defs, expMean, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, expMean, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Build Model with Confidence Intervals modelACE3 <- mxModel( "oneACEvc_3alt", pars, modelMZ, modelDZ, multi, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL 3 # Run ACE Model fitACE3 <- mxRun( modelACE3, intervals=T ) sumACE3 <- summary( fitACE3 ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE3) fitEstCIs(fitACE3,colUS) # this will give an error if you set intervals=F, but still print the parameter estimates # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL 4: univariate twin model + covariates + SIBLING - ALTERNATE PARAMETERIZATION with DEFINITION VARIABLES # see oneACEvc_4def.R for single model script # ---------------------------------------------------------------------------------------------------------------------- # Create Algebra for expected Variance/Covariance Matrices in twins & extra sibling relA <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=FALSE, labels=c("data.relT","data.relS","data.relS"), name="rA" ) expCovTW <- mxAlgebra( expression= VA%x%rA + VC%x%rC + VE%x%rE, name="expCovTW" ) # Create Data Objects for Single Group dataTW <- mxData( observed=tsDataS, type="raw" ) # Create Expectation Objects for Single Group expTW <- mxExpectationNormal( covariance="expCovTW", means="expMean", dimnames=selVars ) # Create Model Objects for Single Group defs <- list( defAge, defSex, relA) pars <- list( intercept, betaS, betaA, covA, covC, covE, covP, relC, relE ) # Build Model with Confidence Intervals modelACE4 <- mxModel( "oneACEvc_4def", pars, defs, expMean, expCovTW, dataTW, expTW, funML, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL 4 # Run ACE Model fitACE4 <- mxRun( modelACE4, intervals=T ) sumACE4 <- summary( fitACE4 ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE4) fitEstCIs(fitACE4,colUS) # this will give an error if you set intervals=F, but still print the parameter estimates # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL 5: univariate twin model + covariates + SIBLING - ALTERNATE SPEC with DEFINITION VARs for ACTUAL RELATEDNESS # see oneACEvc_5rel.R for single model script # ---------------------------------------------------------------------------------------------------------------------- # Create Algebra for expected Variance/Covariance Matrices in twins & extra sibling relA <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=FALSE, labels=c("data.rel12","data.rel13","data.rel23"), name="rA" ) # Create Model Objects for Single Group defs <- list( defAge, defSex, relA) pars <- list( intercept, betaS, betaA, covA, covC, covE, covP, relC, relE ) # Build Model with Confidence Intervals modelACE5 <- mxModel( "oneACEvc_5rel", pars, defs, expMean, expCovTW, dataTW, expTW, funML, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL 5 # Run ACE Model fitACE5 <- mxRun( modelACE5, intervals=T ) sumACE5 <- summary( fitACE5 ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE5) fitEstCIs(fitACE5,colUS) # this will give an error if you set intervals=F, but still print the parameter estimates # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL 6: univariate twin model + covariates + SIBLING - ALT SPEC with DEF VARs for ACTUAL RELATEDNESS for DZs only # see oneACEvc_6dzs.R for single model script # ---------------------------------------------------------------------------------------------------------------------- # Create Data Objects for DZ Group dataDZ <- mxData( observed=tsDataS[tsDataS$zyg==2,], type="raw" ) # Create Model Objects for DZ Group defs <- list( defAge, defSex, relA) pars <- list( intercept, betaS, betaA, covA, covC, covE, covP, relC, relE ) # Build Model with Confidence Intervals modelACE6 <- mxModel( "oneACEvc_6dzs", pars, defs, expMean, expCovTW, dataDZ, expTW, funML, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL 6 # Run ACE Model fitACE6 <- mxRun( modelACE6, intervals=T ) sumACE6 <- summary( fitACE6 ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE6) fitEstCIs(fitACE6,colUS) # this will give an error if you set intervals=F, but still print the parameter estimates # ---------------------------------------------------------------------------------------------------------------------- sink()