rm(list=ls()) library(OpenMx) library(readr) library(haven) LongitudinalTobData <- readRDS("~/Desktop/TobaccoCohort/Data/LongitudinalTobData.rds") LongitudinalTobData$EverCigarette75[LongitudinalTobData$EverCigarette75==1] <- 0 LongitudinalTobData$EverCigarette75[LongitudinalTobData$EverCigarette75==2] <- 1 LongitudinalTobData$EverCigarette75F <- as_factor(LongitudinalTobData$EverCigarette75, levels="values", ordered=T) LongitudinalTobData$RegularCigarette75[LongitudinalTobData$RegularCigarette75==1] <- 0 LongitudinalTobData$RegularCigarette75[LongitudinalTobData$RegularCigarette75==2] <- 1 LongitudinalTobData$RegularCigarette75F <- as_factor(LongitudinalTobData$RegularCigarette75, levels="values", ordered=T) LongitudinalTobData <- subset(LongitudinalTobData, select=c("person_nb","family_nb","Twin","zyg", "sex","Education75", "age75", "EverCigarette75", "RegularCigarette75F", "RegularCigarette75", "CurrentCigarette75", "FormerQuantity75", "CurrentQuantity75", "EverCigarPipe75", "RegularCigarPipe75", "BeerQuant75", "WineQuant75", "SpiritsQuant75", "Binge75", "Status75","Quantity75")) LongitudinalTobDataWide <- reshape(LongitudinalTobData, v.names=c("person_nb","Education75", "age75", "EverCigarette75", "RegularCigarette75F", "RegularCigarette75", "CurrentCigarette75", "FormerQuantity75", "CurrentQuantity75", "EverCigarPipe75", "RegularCigarPipe75", "BeerQuant75", "WineQuant75", "SpiritsQuant75", "Binge75", "Status75","Quantity75"), timevar = 'Twin', idvar=c('family_nb'), direction="wide", sep="") # Select variables varsc <- c("Quantity75") nvc <- 1 ntvc <- nvc*2 conVars <- paste(varsc, c(rep(0,nvc), rep(1,nvc)), sep="") nth <- 1 varso <- c("RegularCigarette75F") nvo <- 1 ntvo <- nvo*2 ordVars <- paste(varso, c(rep(0,nvo), rep(1,nvo)), sep="") oc <- c(1,0) #1 ordinal 0 continuous Vars <- c("RegularCigarette75F", "Quantity75") nv <- nvc+nvo ntv <- nv*2 selVars <- paste(Vars, c(rep(0,nv), rep(1,nv)), sep="") DefVars <- c("sex", "RegularCigarette750", "RegularCigarette751") frMV <- c(F, T) # Select Data for Analysis MzData <- subset(LongitudinalTobDataWide, zyg==1, select=c(selVars, DefVars)) DzData <- subset(LongitudinalTobDataWide, zyg==2, select=c(selVars, DefVars)) MzData <- subset(MzData, !is.na(RegularCigarette750) & !is.na(RegularCigarette751)) DzData <- subset(DzData, !is.na(RegularCigarette750) & !is.na(RegularCigarette751)) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=ntv, ncol=1, free=c(F, T, F, T), values=c(0, 4, 0, 4), labels=c("meanDist","meanCig"), name="meanG" ) SexCov <- mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, labels="data.sex", name="SexCov") regCoefSex <- mxMatrix(type="Full", nrow=ntv, ncol=1, free=TRUE, values=.2, labels=c("SexbetaDist", "SexbetaCig", "SexbetaDist", "SexbetaCig"), name="regCoefSex") EverCov <- mxMatrix(type="Full", nrow=ntv, ncol=1, free=FALSE, labels=c("data.RegularCigarette750", "data.RegularCigarette751"), name="EverCov") regCoefEver <- mxMatrix(type="Full", nrow=1, ncol=ntv, free=c(F, T, F, T), values=c(0, .2, 0, .2), labels=c("NA","B", "NA","B"), name="regCoefEver") expMean <- mxAlgebra(expression=t(meanG+(regCoefSex%*%SexCov)+(regCoefEver%*%EverCov)), name="expMean") threG <- mxMatrix(type="Full", nrow=nth, ncol=ntvo, free=T, values=.01, labels=c("threshold", "threshold"), name="threG") # Create Matrices for Path Coefficients A <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(T, F, T), values=c(.5, 0,.5), labels=c("a1", "a2","a3"), name="A" ) C <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(T, F, T), values=c(.5, 0,.5), labels=c("c1", "c2", "c3"), name="C" ) E <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(T, F, T), values=c(.5, 0,.5), labels=c("e1", "e2", "e3"), name="E" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= A+C+E, name="V" ) covMZ <- mxAlgebra( expression= A+C, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%A+ C, 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" ) #constrain variance of binary variables matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="matUnv" ) matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="matOc" ) var1 <- mxConstraint( expression=diag2vec(matOc%&%V)==matUnv, name="Var1" ) # Putting MZ/DZ data into OpenMx objects 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, thresholds="threG", threshnames=ordVars) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( meanG, SexCov, regCoefSex, EverCov, regCoefEver, expMean, threG, A, C, E, covP, var1, matUnv, matOc ) 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") ) ciACE <- mxCI( c("a", "c", "e") ) # Build Model with Confidence Intervals modelACE <- mxModel( "oneACEc", modelMZ, modelDZ, multi, ciACE ) ################### ### Test Models ### ################### # Run ACE Model fitACE <- mxTryHardOrdinal( modelACE, intervals=F) sumACE <- summary( fitACE, vebose=T) sumACE fitACE$output$matrices$MZ.matUnv fitACE$output$matrices$MZ.matOc fitACE$output$status$code fitACE$output$algebras$MZ.V fitACE$output$matrices$MZ.A fitACE$output$matrices$MZ.C fitACE$output$matrices$MZ.E fitACE$output$matrices$MZ.threG fitACE$output$matrices$MZ.meanG fitACE$output$matrices$MZ.regCoefSex fitACE$output$matrices$MZ.regCoefEver fitACE$output$algebras$MZ.expMean ###################### ### Test Submodels ### ###################### # Run AE model modelAE <- mxModel( fitACE, name="oneAEc" ) modelAE <- omxSetParameters( modelAE, labels=c("c1", "c2", "c3"), free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=F ) # Run CE model modelCE <- mxModel( fitACE, name="oneCEc" ) modelCE <- omxSetParameters( modelCE, labels=c("a1", "a2", "a3"), free=FALSE, values=0 ) fitCE <- mxTryHard( modelCE, intervals=F ) # Run E model modelE <- mxModel( fitAE, name="oneEc" ) modelE <- omxSetParameters( modelE, labels=c("a1", "a2", "a3"), free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=F ) # Print Comparative Fit Statistics mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )