# YOU MAY ALSO VISIT THIS WEBSITE FOR THE TUTORIAL TO GO WITH THIS CODE # https://qimr.az1.qualtrics.com/jfe/form/SV_eqWgUTNJghdOBVA # ---------------------------------------------------------------------------------------------------------------------- # Program: oneSATca.R # Author: Hermine Maes # Website: https://hermine-maes.squarespace.com/openmx#/one/ # Date: 01 04 2018 # # Twin Univariate Saturated model to estimate means and (co)variances across multiple groups # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # PREPARE DATA # Setting working directory. # Use workshop server address. Put YOUR folder address here. Showing mine as an example ONLY. setwd('/home/elizabeth/2022') # Load Libraries & Options #rm(list=ls()) #library(OpenMx) #library(psych) #library(polycor) #source("miFunctions5272022.R") # Create Output #filename <- "oneSATca" #sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE) # ---------------------------------------------------------------------------------------------------------------------- # Open site-specific data # ODD NUMBER ROOMS- Open SimWtDataSite1D1.csv # EVEN NUMBER ROOMS- Open SimWtDataSite2D1.csv #Twins <- read.table(file="SimWtDataSite1D1.csv",header=T,sep=",") #Twins <- read.table(file="SimWtDataSite2D1.csv",header=T,sep=",") dim(Twins) describe(Twins[,1:5], skew=F) #################################################################### #### Use of Saturated Model to Test Assumptions-- ONE SITE #### #################################################################### # ---------------------------------------------------------------------------------------------------------------------- # Program: oneSATc.R # Author: Hermine Maes # Date: 10 22 2018 # Website: https://hermine-maes.squarespace.com/openmx#/one/ # Twin Univariate Saturated model to estimate means and (co)variances across multiple groups # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # Load Libraries & Options #rm(list=ls()) library(OpenMx) mxOption( NULL, "Default optimizer", "NPSOL" ) library(psych) library(polycor) source("miFunctions5272022.R") # Create Output--Commented out for now # filename <- "oneSATc" # sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE) # ---------------------------------------------------------------------------------------------------------------------- # Select Variables for Analysis vars <- 'WT' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="_T") # Select Data for Analysis mzData <- subset(Twins, ZYG==1, selVars) dzData <- subset(Twins, ZYG==2, selVars) # Calcuating means, variances, and covariances by twin member (Twin 1 or Twin2) and zygosity groups in OpenMx colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") # Set Starting Values svMe <- 50 # start value for means svVa <- 3 # start value for variance lbVa <- 0.00001 # lower bound for variance # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Algebra for expected Mean Matrices meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mMZ1","mMZ2"), name="meanMZ" ) meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mDZ1","mDZ2"), name="meanDZ" ) # Create Algebra for expected Variance/Covariance Matrices covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv), labels=c("vMZ1","cMZ21","vMZ2"), name="covMZ" ) covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv), labels=c("vDZ1","cDZ21","vDZ2"), name="covDZ" ) # 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="covMZ", means="meanMZ", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="covDZ", means="meanDZ", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups modelMZ <- mxModel( meanMZ, covMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( meanDZ, covDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Confidence Interval Objects ciCov <- mxCI( c('MZ.covMZ','DZ.covDZ') ) ciMean <- mxCI( c('MZ.meanMZ','DZ.meanDZ') ) # Build Saturated Model with Confidence Intervals modelSAT <- mxModel( "oneSATc", modelMZ, modelDZ, multi, ciCov, ciMean ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run Saturated Model fitSAT <- mxRun( modelSAT, intervals=F ) sumSAT <- summary( fitSAT ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitSAT) fitEsts(fitSAT) mxGetExpected( fitSAT, c("means","covariance") ) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Constrain expected Means to be equal across Twin Order modelEMO <- mxModel( fitSAT, name="oneEMOc" ) modelEMO <- omxSetParameters( modelEMO, label=c("mMZ1","mMZ2"), free=TRUE, values=svMe, newlabels='mMZ' ) modelEMO <- omxSetParameters( modelEMO, label=c("mDZ1","mDZ2"), free=TRUE, values=svMe, newlabels='mDZ' ) fitEMO <- mxRun( modelEMO, intervals=F ) fitGofs(fitEMO); fitEsts(fitEMO) # Constrain expected Means and Variances to be equal across Twin Order modelEMVO <- mxModel( fitEMO, name="oneEMVOc" ) modelEMVO <- omxSetParameters( modelEMVO, label=c("vMZ1","vMZ2"), free=TRUE, values=svVa, newlabels='vMZ' ) modelEMVO <- omxSetParameters( modelEMVO, label=c("vDZ1","vDZ2"), free=TRUE, values=svVa, newlabels='vDZ' ) fitEMVO <- mxRun( modelEMVO, intervals=F ) fitGofs(fitEMVO); fitEsts(fitEMVO) # Constrain expected Means and Variances to be equal across Twin Order and Zygosity modelEMVZ <- mxModel( fitEMVO, name="oneEMVZc" ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("mMZ","mDZ"), free=TRUE, values=svMe, newlabels='mZ' ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("vMZ","vDZ"), free=TRUE, values=svVa, newlabels='vZ' ) fitEMVZ <- mxRun( modelEMVZ, intervals=F ) fitGofs(fitEMVZ); fitEsts(fitEMVZ) # Print Comparative Fit Statistics mxCompare( fitSAT, subs <- list(fitEMO, fitEMVO, fitEMVZ) ) #### ACE Model 1- One Site Only #### # Set Starting Values svMe <- 50 # start value for means svPa <- .2 # start value for path coefficient svPe <- .5 # start value for path coefficient for e # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" ) # Create Matrices for Variance Components covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" ) covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" ) covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, 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="meanG", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( meanG, covA, covC, covE, covP ) modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Unstandardized and 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 modelACE <- mxModel( "oneACEvc", pars, modelMZ, modelDZ, multi, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=T ) sumACE <- summary( fitACE ) # Compare with Saturated Model if saturated model was fitted in same session and if saturated model prior to genetic model mxCompare( fitSAT, fitACE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEstCis(fitACE) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Run AE model modelAE <- mxModel( fitACE, name="oneAEvc" ) modelAE <- omxSetParameters( modelAE, labels="VC11", free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=T ) fitGofs(fitAE); fitEstCis(fitAE) # Run CE model modelCE <- mxModel( fitACE, name="oneCEvc" ) modelCE <- omxSetParameters( modelCE, labels="VA11", free=FALSE, values=0 ) modelCE <- omxSetParameters( modelCE, labels=c("VE11","VC11"), free=TRUE, values=.6 ) fitCE <- mxRun( modelCE, intervals=T ) fitGofs(fitCE); fitEstCis(fitCE) # Run E model modelE <- mxModel( fitAE, name="oneEvc" ) modelE <- omxSetParameters( modelE, labels="VA11", free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) fitGofs(fitE); fitEstCis(fitE) # Print Comparative Fit Statistics mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ) round(rbind(fitACE$US$result,fitAE$US$result,fitCE$US$result,fitE$US$result),4) #### FULL DATASET - ACE MODEL WITHOUT SITE ##### #### ACE Model1 #### Twins2<- read.table(file='SimWtDataPairD1.csv',header=T,sep=",") # Select Variables for Analysis vars <- 'WT' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="_T") # Select Data for Analysis mzData <- subset(Twins2, ZYG==1, selVars) dzData <- subset(Twins2, ZYG==2, selVars) # Set Starting Values svMe <- 50 # start value for means svPa <- .2 # start value for path coefficient svPe <- .5 # start value for path coefficient for e # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" ) # Create Matrices for Variance Components covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" ) covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" ) covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, 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="meanG", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( meanG, covA, covC, covE, covP ) modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Unstandardized and 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 modelACE <- mxModel( "oneACEvc", pars, modelMZ, modelDZ, multi, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=TRUE ) sumACE <- summary( fitACE ) # Compare with Saturated Model if saturated model was fitted in same session and if saturated model prior to genetic model #mxCompare( fitSAT, fitACE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEstCis(fitACE) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Run AE model modelAE <- mxModel( fitACE, name="oneAEvc" ) modelAE <- omxSetParameters( modelAE, labels="VC11", free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=T ) fitGofs(fitAE); fitEstCis(fitAE) # Run CE model modelCE <- mxModel( fitACE, name="oneCEvc" ) modelCE <- omxSetParameters( modelCE, labels="VA11", free=FALSE, values=0 ) modelCE <- omxSetParameters( modelCE, labels=c("VE11","VC11"), free=TRUE, values=.6 ) fitCE <- mxRun( modelCE, intervals=T ) fitGofs(fitCE); fitEstCis(fitCE) # Run E model modelE <- mxModel( fitAE, name="oneEvc" ) modelE <- omxSetParameters( modelE, labels="VA11", free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) fitGofs(fitE); fitEstCis(fitE) # Print Comparative Fit Statistics mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ) round(rbind(fitACE$US$result,fitAE$US$result,fitCE$US$result,fitE$US$result),4) #### FULL DATASET- ANALYSES WITH SITE ##### ### Saturated Model ### # Load Data Twins3<- read.table(file='SimWtDataPairD1.csv',header=T,sep=",") dim(Twins3) describe(Twins3[,1:5], skew=F) # Select Variables for Analysis vars <- 'WT' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="_T") # Select Covariates for Analysis. Make sure that are no observations that are missing on the covariate. Twins3a <- subset(Twins3, !is.na(Twins3$Site)) covVars <- 'Site' # Select Data for Analysis mzData <- subset(Twins3a, ZYG==1, c(selVars,covVars)) dzData <- subset(Twins3a, ZYG==2, c(selVars,covVars)) # Generate Descriptive Statistics colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") # Set Starting Values svBe <- 0.01 # start value for regressions svMe <- 60 # start value for means svVa <- 3 # start value for variance lbVa <- .0001 # lower bound for variance # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Matrices for Covariates and linear Regression Coefficients defL <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Site"), name="defL" ) pathBl <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="b11", name="bl" ) # Create Algebra for expected Mean Matrices meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mMZ1","mMZ2"), name="meanMZ" ) meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mDZ1","mDZ2"), name="meanDZ" ) expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(defL%*%bl,defL%*%bl), name="expMeanMZ" ) expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(defL%*%bl,defL%*%bl), name="expMeanDZ" ) # Create Algebra for expected Variance/Covariance Matrices covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv), labels=c("vMZ1","cMZ21","vMZ2"), name="covMZ" ) covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv), labels=c("vDZ1","cDZ21","vDZ2"), name="covDZ" ) # 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="covMZ", means="expMeanMZ", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( pathBl ) defs <- list( defL ) modelMZ <- mxModel( pars, defs, meanMZ, expMeanMZ, covMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, meanDZ, expMeanDZ, covDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Confidence Interval Objects ciCov <- mxCI( c('MZ.covMZ','DZ.covDZ') ) ciMean <- mxCI( c('MZ.meanMZ','DZ.meanDZ') ) # Build Saturated Model with Confidence Intervals modelSAT <- mxModel( "oneSATca", pars, modelMZ, modelDZ, multi, ciCov, ciMean ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run Saturated Model fitSAT <- mxRun( modelSAT, intervals=F ) sumSAT <- summary( fitSAT ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitSAT) fitEsts(fitSAT) mxGetExpected( fitSAT, c("means","covariance") ) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Test Significance of Covariate modelCOV <- mxModel( fitSAT, name="oneCOVca" ) modelCOV <- omxSetParameters( modelCOV, label="b11", free=FALSE, values=0 ) fitCOV <- mxRun( modelCOV ) # Constrain expected Means to be equal across Twin Order modelEMO <- mxModel( fitSAT, name="oneEMOca" ) modelEMO <- omxSetParameters( modelEMO, label=c("mMZ1","mMZ2"), free=TRUE, values=svMe, newlabels='mMZ' ) modelEMO <- omxSetParameters( modelEMO, label=c("mDZ1","mDZ2"), free=TRUE, values=svMe, newlabels='mDZ' ) fitEMO <- mxRun( modelEMO, intervals=F ) fitGofs(fitEMO); fitEsts(fitEMO) # Constrain expected Means and Variances to be equal across Twin Order modelEMVO <- mxModel( fitEMO, name="oneEMVOca" ) modelEMVO <- omxSetParameters( modelEMVO, label=c("vMZ1","vMZ2"), free=TRUE, values=svVa, newlabels='vMZ' ) modelEMVO <- omxSetParameters( modelEMVO, label=c("vDZ1","vDZ2"), free=TRUE, values=svVa, newlabels='vDZ' ) fitEMVO <- mxRun( modelEMVO, intervals=F ) fitGofs(fitEMVO); fitEsts(fitEMVO) # Constrain expected Means and Variances to be equal across Twin Order and Zygosity modelEMVZ <- mxModel( fitEMVO, name="oneEMVZca" ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("mMZ","mDZ"), free=TRUE, values=svMe, newlabels='mZ' ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("vMZ","vDZ"), free=TRUE, values=svVa, newlabels='vZ' ) fitEMVZ <- mxRun( modelEMVZ, intervals=F ) fitGofs(fitEMVZ); fitEsts(fitEMVZ) # Print Comparative Fit Statistics mxCompare( fitSAT, subs <- list(fitCOV, fitEMO, fitEMVO, fitEMVZ) ) # ---------------------------------------------------------------------------------------------------------------------- #sink() #save.image(paste(filename,".Ri",sep="")) # ---------------------------------------------------------------------------------------------------------------------- # Program: oneACEvca.R # Author: Hermine Maes # Date: 02 15 2020 # # Twin Univariate ACE model to estimate causes of variation across multiple groups # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # Create Output #filename <- "oneACEvca" #sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE) # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA # Load Data with both study site locations Twins3<- read.table(file='SimWtDataPairD1.csv',header=T,sep=",") dim(Twins3) describe(Twins3[,1:5], skew=F) # Select Variables for Analysis vars <- 'WT' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="_T") # Select Covariates for Analysis. Make sure that are no observations that are missing on the covariate. Twins3a <- subset(Twins3, !is.na(Twins3$Site)) covVars <- 'Site' # Select Data for Analysis mzData <- subset(Twins3a, ZYG==1, c(selVars,covVars)) dzData <- subset(Twins3a, ZYG==2, c(selVars,covVars)) # Generate Descriptive Statistics in OpenMx colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") # Set Starting Values svBe <- 0.01 # start value for regressions svMe <- 50 # start value for means svPa <- .2 # start value for path coefficient svPe <- .5 # start value for path coefficient for e # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Matrices for Covariates and linear Regression Coefficients defL <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Site"), name="defL" ) pathBl <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="b11", name="bl" ) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" ) expMean <- mxAlgebra( expression= meanG + cbind(defL%*%bl,defL%*%bl), name="expMeanG" ) # Create Matrices for Variance Components covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" ) covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" ) covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, 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="expMeanG", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanG", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( pathBl, meanG, covA, covC, covE, covP ) defs <- list( defL ) 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 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 modelACE <- mxModel( "oneACEvca", pars, modelMZ, modelDZ, multi, estUS, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=T ) sumACE <- summary( fitACE ) # Compare with Saturated Model if saturated model fitted in same session if saturated model prior to genetic model #mxCompare( fitSAT, fitACE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEstCis(fitACE) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Run AE model modelAE <- mxModel( fitACE, name="oneAEvca" ) modelAE <- omxSetParameters( modelAE, labels="VC11", free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=T ) fitGofs(fitAE); fitEstCis(fitAE) # Run CE model modelCE <- mxModel( fitACE, name="oneCEvca" ) modelCE <- omxSetParameters( modelCE, labels="VA11", free=FALSE, values=0 ) modelCE <- omxSetParameters( modelCE, labels=c("VE11","VC11"), free=TRUE, values=.6 ) fitCE <- mxRun( modelCE, intervals=T ) fitGofs(fitCE); fitEstCis(fitCE) # Run E model modelE <- mxModel( fitAE, name="oneEvca" ) modelE <- omxSetParameters( modelE, labels="VA11", free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) fitGofs(fitE); fitEstCis(fitE) # Print Comparative Fit Statistics mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ) round(rbind(fitACE$US$result,fitAE$US$result,fitCE$US$result,fitE$US$result),4) # ---------------------------------------------------------------------------------------------------------------------- #sink() #save.image(paste(filename,".Ri",sep=""))