# ------------------------------------------------------------------------------ # Title: Biv_WB_SC_a7 for forum OpenMx # Author: based on script Hermine Maes # Date: 07/10/2016 # Problem: ACE estimates do not match twin correlations # # Twin Bivariate Saturated model to estimate means and (co)variances across multiple groups # # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------| # Load Libraries & Options library(OpenMx) library(psych) source("miFunctions.R") #for now, I the use NPSOL optimizer # Create Output filename <- "Biv_WB_SC_a7" sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE) # ------------------------------------------------------------------------------ # PREPARE DATA # Load Data Data <- read.delim("U:/data/SC scale/September2016/SCscale/fortwins/fortwinsage7/persvar_a7_twins.dat", dec=",", na.strings="-99") describe(Data, skew=F) #matches output SPSS #WB7_tw1 3 2010 8.40 0.98 8.00 8.44 1.48 3.00 10.00 7.00 0.02 #WB7_tw2 4 2006 8.41 0.95 8.00 8.43 1.48 0.00 10.00 10.00 0.02 #SC7M_tw1 5 11906 3.67 3.21 3.20 3.30 2.67 0.00 16.00 16.00 0.03 #SC7M_tw2 6 11875 3.18 3.02 3.00 2.77 2.97 0.00 16.00 16.00 0.03 colnames(Data) colnames(Data)[1] <- "age7" colnames(Data)[3] <- "WB_tw1" colnames(Data)[4] <- "WB_tw2" colnames(Data)[5] <- "SC_tw1" colnames(Data)[6] <- "SC_tw2" colnames(Data) # Select Variables for Analysis Vars <- c('WB_tw','SC_tw') nv <- 2 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") #c('WB_tw1","SC_tw1',WB_tw2","SC_tw2) selVars # Select Data for Analysis mzmData <- subset(Data, zyg5==1, selVars) dzmData <- subset(Data, zyg5==2, selVars) mzfData <- subset(Data, zyg5==3, selVars) dzfData <- subset(Data, zyg5==4, selVars) dosmfData <- subset(Data, zyg5==5, selVars) # Generate Descriptive Statistics colMeans(mzmData,na.rm=TRUE) colMeans(dzmData,na.rm=TRUE) colMeans(mzfData,na.rm=TRUE) colMeans(dzfData,na.rm=TRUE) colMeans(dosmfData,na.rm=TRUE) cov(mzmData,use="complete") cov(dzmData,use="complete") cov(mzfData,use="complete") cov(dzfData,use="complete") cov(dosmfData,use="complete") cor(mzmData,use="complete") cor(dzmData,use="complete") cor(mzfData,use="complete") cor(dzfData,use="complete") cor(dosmfData,use="complete") # Set Starting Values svMe <- c(7.5,3) # start value for means svVa <- c(.7,4) # start value for variances svVas <- diag(svVa,ntv,ntv) # assign start values to diagonal of matrix # Create Labels labMeMZM <- paste("meanMZM",selVars,sep="_") labMeDZM <- paste("meanDZM",selVars,sep="_") labMeMZF <- paste("meanMZF",selVars,sep="_") labMeDZF <- paste("meanDZF",selVars,sep="_") labMeDOSMF <- paste("meanDOSMF",selVars,sep="_") labMeZ <- paste("meanZ",selVars,sep="_") labCvMZM <- labLower("covMZM",ntv) labCvDZM <- labLower("covDZM",ntv) labCvMZF <- labLower("covMZF",ntv) labCvDZF <- labLower("covDZF",ntv) labCvDOSMF <- labLower("covDOSMF",ntv) labCvZ <- labLower("covZ",ntv) labVaMZM <- labDiag("covMZM",ntv) labVaDZM <- labDiag("covDZM",ntv) labVaMZF <- labDiag("covMZF",ntv) labVaDZF <- labDiag("covDZF",ntv) labVaDOSMF <- labDiag("covDOSMF",ntv) labVaZ <- labDiag("covZ",ntv) # ------------------------------------------------------------------------------ # PREPARE MODEL # Saturated Model # Create Algebra for expected Mean Matrices meanMZM <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZM, name="meanMZM" ) meanDZM <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZM, name="meanDZM" ) meanMZF <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZF, name="meanMZF" ) meanDZF <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZF, name="meanDZF" ) meanDOSMF <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDOSMF, name="meanDOSMF" ) # Create Algebra for expected Variance/Covariance Matrices covMZM <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=labCvMZM, name="covMZM" ) #lbound=lbVas, covDZM <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=labCvDZM, name="covDZM" ) #lbound=lbVas, covMZF <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=labCvMZF, name="covMZF" ) #lbound=lbVas, covDZF <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=labCvDZF, name="covDZF" ) #lbound=lbVas, covDOSMF <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=labCvDOSMF, name="covDOSMF" ) #lbound=lbVas, # Create Algebra for Maximum Likelihood Estimates of Twin Correlations matI <- mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="I" ) corMZM <- mxAlgebra( solve(sqrt(I*covMZM)) %&% covMZM, name="corMZM" ) corDZM <- mxAlgebra( solve(sqrt(I*covDZM)) %&% covDZM, name="corDZM" ) corMZF <- mxAlgebra( solve(sqrt(I*covMZF)) %&% covMZF, name="corMZF" ) corDZF <- mxAlgebra( solve(sqrt(I*covDZF)) %&% covDZF, name="corDZF" ) corDOSMF <- mxAlgebra( solve(sqrt(I*covDOSMF)) %&% covDOSMF, name="corDOSMF" ) # Create Data Objects for Multiple Groups dataMZM <- mxData( observed=mzmData, type="raw" ) dataDZM <- mxData( observed=dzmData, type="raw" ) dataMZF <- mxData( observed=mzfData, type="raw" ) dataDZF <- mxData( observed=dzfData, type="raw" ) dataDOSMF <- mxData( observed=dosmfData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZM <- mxExpectationNormal( covariance="covMZM", means="meanMZM", dimnames=selVars ) expDZM <- mxExpectationNormal( covariance="covDZM", means="meanDZM", dimnames=selVars ) expMZF <- mxExpectationNormal( covariance="covMZF", means="meanMZF", dimnames=selVars ) expDZF <- mxExpectationNormal( covariance="covDZF", means="meanDZF", dimnames=selVars ) expDOSMF <- mxExpectationNormal( covariance="covDOSMF", means="meanDOSMF", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups modelMZM <- mxModel( name="MZM", meanMZM, covMZM, matI, corMZM, dataMZM, expMZM, funML ) modelDZM <- mxModel( name="DZM", meanDZM, covDZM, matI, corDZM, dataDZM, expDZM, funML ) modelMZF <- mxModel( name="MZF", meanMZF, covMZF, matI, corMZF, dataMZF, expMZF, funML ) modelDZF <- mxModel( name="DZF", meanDZF, covDZF, matI, corDZF, dataDZF, expDZF, funML ) modelDOSMF <- mxModel( name="DOSMF", meanDOSMF, covDOSMF, matI, corDOSMF, dataDOSMF, expDOSMF, funML ) multi <- mxFitFunctionMultigroup( c("MZM","DZM", "MZF", "DZF","DOSMF") ) # Create Confidence Interval Objects ciCov <- mxCI( c('MZM.covMZM','DZM.covDZM','MZF.covMZF','DZF.covDZF','DOSMF.covDOSMF') ) ciMean <- mxCI( c('MZM.meanMZM','DZM.meanDZM','MZF.meanMZF','DZF.meanDZF','DOSMF.meanDOSMF') ) ciCor <- mxCI( c('MZM.corMZM','DZM.corDZM','MZF.corMZF','DZF.corDZF','DOSMF.corDOSMF') ) # Build Saturated Model with Confidence Intervals model <- mxModel( "fiveSATc", modelMZM, modelDZM, modelMZF, modelDZF, modelDOSMF, multi, ciMean,ciCov, ciCor ) # ------------------------------------------------------------------------------ # RUN MODEL # Run Saturated Model fit <- mxTryHard( model, intervals=F ) sum <- summary( fit ) sum round(fit$output$estimate,4) fitGofs(fit) ExpCovMZm <- mxEval(MZM.corMZM, fit) ExpCovMZm ExpCovDZm <- mxEval(DZM.corDZM, fit) ExpCovDZm ExpCovMZf <- mxEval(MZF.corMZF, fit) ExpCovMZf ExpCovDZf <- mxEval(DZF.corDZF, fit) ExpCovDZf ExpCovDOSMF <- mxEval(DOSMF.corDOSMF, fit) ExpCovDOSMF fit$output$confidenceIntervals # ------------------------------------------------------------------------------ # ACE MODEL # -------|---------|---------|---------|---------|---------|---------|---------| # Load Libraries & Options # Create Output filename <- "twoACEc" sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE) # ------------------------------------------------------------------------------ # PREPARE DATA #------------------------------------------------------------------------------ # Load Data # Set Starting Values svPa <- .4 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .8 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix # Create Labels labMeM <- paste("meanM",Vars,sep="_") labMeF <- paste("meanF",Vars,sep="_") labMeO <- paste("meanO",Vars,sep="_") # ------------------------------------------------------------------------------ # PREPARE ACE MODEL # ------------------------------------------------------------------------------- # Create Algebra for expected Mean Matrices meanGm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeM, name="meanGm" ) meanGf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeF, name="meanGf" ) meanGo <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeO, name="meanGo" ) # Create Matrices for Path Coefficients pathAm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("am",nv), name="am" ) #lbound=lbPaD, pathCm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("cm",nv), name="cm" )#lbound=lbPaD, pathEm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("em",nv), name="em" )#lbound=lbPaD, pathAf <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("af",nv), name="af" ) #lbound=lbPaD, pathCf <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("cf",nv), name="cf" )#lbound=lbPaD, pathEf <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("ef",nv), name="ef" )#lbound=lbPaD, pathRg <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="rg11", lbound=0, ubound=1, name="rg" ) # Create Algebra for Variance Comptwonts covAm <- mxAlgebra( expression=am %*% t(am), name="Am" ) covCm <- mxAlgebra( expression=cm %*% t(cm), name="Cm" ) covEm <- mxAlgebra( expression=em %*% t(em), name="Em" ) covAf <- mxAlgebra( expression=af %*% t(af), name="Af" ) covCf <- mxAlgebra( expression=cf %*% t(cf), name="Cf" ) covEf <- mxAlgebra( expression=ef %*% t(ef), name="Ef" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covPm <- mxAlgebra( expression= Am+Cm+Em, name="Vm" ) covPf <- mxAlgebra( expression= Af+Cf+Ef, name="Vf" ) covMZm <- mxAlgebra( expression= Am+Cm, name="cMZm" ) covMZf <- mxAlgebra( expression= Af+Cf, name="cMZf" ) covDZm <- mxAlgebra( expression= 0.5%x%Am+ Cm, name="cDZm" ) covDZf <- mxAlgebra( expression= 0.5%x%Af+ Cf, name="cDZf" ) covDZfm <- mxAlgebra( expression= 0.5%*%rg%x%(af%*%t(am))+cf%*%t(cm), name="covDZfm" ) covDZmf <- mxAlgebra( expression= 0.5%*%rg%x%(am%*%t(af))+cm%*%t(cf), name="covDZmf" ) expCovMZm <- mxAlgebra( expression= rbind( cbind(Vm, cMZm), cbind(t(cMZm), Vm)), name="expCovMZm" ) expCovMZf <- mxAlgebra( expression= rbind( cbind(Vf, cMZf), cbind(t(cMZf), Vf)), name="expCovMZf" ) expCovDZm <- mxAlgebra( expression= rbind( cbind(Vm, cDZm), cbind(t(cDZm), Vm)), name="expCovDZm" ) expCovDZf <- mxAlgebra( expression= rbind( cbind(Vf, cDZf), cbind(t(cDZf), Vf)), name="expCovDZf" ) expCovDosmf <- mxAlgebra( expression= rbind( cbind(Vm, covDZmf), cbind(covDZfm, Vf)), name="expCovDosmf" ) # Create Algebra for Standardization matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSDm <- mxAlgebra( expression=solve(sqrt(I*Vm)), name="iSDm") invSDf <- mxAlgebra( expression=solve(sqrt(I*Vf)), name="iSDf") # Calculate genetic and environmental correlations corAm <- mxAlgebra( expression=solve(sqrt(I*Am))%&%Am, name ="rAm" ) #cov2cor() corCm <- mxAlgebra( expression=solve(sqrt(I*Cm))%&%Cm, name ="rCm" ) corEm <- mxAlgebra( expression=solve(sqrt(I*Em))%&%Em, name ="rEm" ) corAf <- mxAlgebra( expression=solve(sqrt(I*Af))%&%Af, name ="rAf" ) #cov2cor() corCf <- mxAlgebra( expression=solve(sqrt(I*Cf))%&%Cf, name ="rCf" ) corEf <- mxAlgebra( expression=solve(sqrt(I*Ef))%&%Ef, name ="rEf" ) # Create Data Objects for Multiple Groups dataMZM <- mxData( observed=mzmData, type="raw" ) dataDZM <- mxData( observed=dzmData, type="raw" ) dataMZF <- mxData( observed=mzfData, type="raw" ) dataDZF <- mxData( observed=dzfData, type="raw" ) dataDOSMF <- mxData( observed=dosmfData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZm <- mxExpectationNormal( covariance="expCovMZm", means="meanGm", dimnames=selVars) expMZf <- mxExpectationNormal( covariance="expCovMZf", means="meanGf", dimnames=selVars) expDZm <- mxExpectationNormal( covariance="expCovDZm", means="meanGm", dimnames=selVars) expDZf <- mxExpectationNormal( covariance="expCovDZf", means="meanGf", dimnames=selVars) expDOSMF <- mxExpectationNormal( covariance="expCovDosmf", means="meanGo", dimnames=selVars) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups parsm <- list(meanGm, matI, invSDm, #matUnv, matOc, pathAm, pathCm, pathEm, covAm, covCm, covEm, covPm, corAm, corCm, corEm) parsf <- list(meanGf, matI, invSDf,#matUnv, matOc, pathAf, pathCf, pathEf, covAf, covCf, covEf, covPf, corAf, corCf, corEf) parso <- list(meanGo , pathRg, covDZmf, covDZfm) modelMZm <- mxModel( name="MZM", parsm, covMZm, expCovMZm, dataMZM, expMZm, funML ) modelMZf <- mxModel( name="MZF", parsf, covMZf, expCovMZf, dataMZF, expMZf, funML ) modelDZm <- mxModel( name="DZM", parsm, covDZm, expCovDZm, dataDZM, expDZm, funML ) modelDZf <- mxModel( name="DZF", parsf, covDZf, expCovDZf, dataDZF, expDZf, funML ) modelDZo <- mxModel( name="DOSMF", parsf, parsm, parso, covDZfm, covDZmf, expCovDosmf, dataDOSMF, expDOSMF, funML ) multi <- mxFitFunctionMultigroup( c("MZM","DZM","MZF","DZF","DOSMF" ) ) # Create Algebra for Variance Components rowVC <- rep('VC',nv) colVCm <- rep(c('Am','Cm','Em','SAm','SCm','SEm'),each=nv) colVCf <- rep(c('Af','Cf','Ef','SAf','SCf','SEf'),each=nv) estVCm <- mxAlgebra( expression=cbind(Am,Cm,Em,Am/Vm,Cm/Vm,Em/Vm), name="VCm", dimnames=list(rowVC,colVCm)) estVCf <- mxAlgebra( expression=cbind(Af,Cf,Ef,Af/Vf,Cf/Vf,Ef/Vf), name="VCf", dimnames=list(rowVC,colVCf)) # Create Confidence Interval Objects ciACE <- mxCI( c("VCf[1,7:12]", "VCf[2,7:12]", "VCm[1,7:12]","VCm[2,7:12]" ) ) cicor <- mxCI( c('rAf', 'rCf', 'rEf', 'rAm', 'rCm', 'rEm')) # Build Model with Confidence Intervals modelACE <- mxModel( "twoACEc", parsm, parsf, parso, modelMZm, modelMZf, modelDZm, modelDZf, modelDZo, multi, estVCm, estVCf, ciACE,cicor ) #var1m, var1f, # ------------------------------------------------------------------------------ # RUN ACE model # ------------------------------------------------------------------------------ # order # ACE seperate boys and girls # ACE equal for boys & girls # AE seperate for boys and girls (so for WB, SC and crosstrait) # Run ACE Model fitACE <- mxTryHard( modelACE, intervals=T ) summary(fitACE) # Compare with Saturated Model mxCompare( fit, fitACE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEsts(fitACE) # ------------------------------------------------------------------------------ # RUN SUBMODELS # ------------------------------------------------------------------------------ # check equality a, c, e for boys and girls #Constrain a , c, & e to be equal for boys and girls. Doing so shows you if we should present a , c, & e seperate for b/g modelgender <- mxModel( fitACE, name = "ace equal boys and girls") modelgender <- omxSetParameters(modelgender, label= c("am_1_1","af_1_1"), newlabels= c("am_1_1"), values=2, free=T ) modelgender <- omxSetParameters(modelgender, label= c("am_2_1","af_2_1"), newlabels= c("am_2_1"), values=2, free=T ) modelgender <- omxSetParameters(modelgender, label= c("am_2_2","af_2_2"), newlabels= c("am_2_2"), values=2, free=T ) modelgender <- omxSetParameters(modelgender, label= c("cm_1_1","cf_1_1"), newlabels= c("cm_1_1"), values=2, free=T ) modelgender <- omxSetParameters(modelgender, label= c("cm_2_1","cf_2_1"), newlabels= c("cm_2_1"), values=2, free=T ) modelgender <- omxSetParameters(modelgender, label= c("cm_2_2","cf_2_2"), newlabels= c("cm_2_2"), values=2, free=T ) modelgender <- omxSetParameters(modelgender, label= c("em_1_1","ef_1_1"), newlabels= c("em_1_1"), values=2, free=T ) modelgender <- omxSetParameters(modelgender, label= c("em_2_1","ef_2_1"), newlabels= c("em_2_1"), values=2, free=T ) modelgender <- omxSetParameters(modelgender, label= c("em_2_2","ef_2_2"), newlabels= c("em_2_2"), values=2, free=T ) modelgenderfit <- mxRun(modelgender) summary(modelgenderfit) mxCompare(fitACE,modelgenderfit) # Run ACE for males, ACE wb AE cross ACE females modelACEAE <- mxModel( fitACE, name= "ACEAE") modelACEAE <- omxSetParameters( modelACEAE, labels= "cf_2_1", free=FALSE, values=0 ) fitACEAE <- mxRun(modelACEAE, intervals=F) mxCompare(fitACE, fitACEAE) summary(fitACEAE) # ------------------------------------------------------------------------------ sink() save.image(paste(filename,".Ri",sep=""))