# Open libraries require(OpenMx) require(psych) # Reads data data <- read.csv ('data.csv', header=T) #Make the integer variables ordered factors data$ebent1 <-mxFactor(data$ebent1, levels=c(0:1) ) data$ebent2 <-mxFactor(data$ebent2, levels=c(0:1) ) data$Sex1 <-mxFactor(data$Sex1, levels=c(0:1) ) data$Sex2 <-mxFactor(data$Sex2, levels=c(0:1) ) data$partner1 <-mxFactor(data$partner1, levels=c(0:1) ) data$partner2 <-mxFactor(data$partner2, levels=c(0:1) ) #uncomment the below to include PHU as categorical variable instead of continuous-ordinal data$PHU1 <- data$PHU1 -1 # this changes the PHU categories from 1:3 to 0:2 data$PHU2 <- data$PHU2 -1 data$PHU1 <-mxFactor(data$PHU1, levels=c(0:2) ) data$PHU2 <-mxFactor(data$PHU2, levels=c(0:2) ) selVars <- c('ebent1','ebent2') useVars <- c('ebent1','ebent2', 'zyg', 'Sex1', 'Sex2','Age1', 'Age2', 'PHU1', 'PHU2', 'partner1', 'partner2') # Select Data for Analysis mzData <- subset(data, zyg==1, useVars) dzData <- subset(data, zyg==2, useVars) describe(mzData) describe(dzData) nv <- 1 # number of variables per twin ntv <- nv*2 # number of variables per pair nth <- 1 # number of max thresholds (I changed this to 1 since I only have 2 categories (0,1)) # 1) Specify Saturated Model (Polychoric correlations) with covariate effects on the thresholds # --------------------------------------------------------------------------------------- # Define definition variables to hold the Covariates obsAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.Age1","data.Age2"), name="Age") obsSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.Sex1","data.Sex2"), name="Sex") obsPHU <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.PHU1","data.PHU2"), name="PHU") obspartner <- mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.partner1","data.partner2"), name="partner") # Matrix & Algebra for expected means (SND), Cor, Thresholds, Covariates MeanL <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" ) betaA <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.1, labels=c('BaTH'), name="BageTH" ) betaS <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.1, labels=c('BsTH'), name="BsexTH" ) betaPHU <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.1, labels=c('BpTH'), name="BphuTH" ) betapart <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.1, labels=c('BparTH'), name="BpartTH" ) Inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="L") Tmz <-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=c(.8, 1), lbound=c(-2,.001), ubound=2, labels=c("Tmz11","imz11","Tmz12","imz12"), name="ThMZ") ThresMZ <-mxAlgebra( expression= L%*%ThMZ + BageTH%x%Age + BsexTH%x%Sex + BphuTH%x%PHU + BpartTH%x%partner, name="expThresMZ") CorMZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=.6, lbound=-.99, ubound=.99, name="expCorMZ") Tdz <-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=c(.8, 1), lbound=c(-2,.001), ubound=2, labels=c("Tdz11","idz11","Tdz12","idz12"), name="ThDZ") ThresDZ <-mxAlgebra( expression= L%*%ThDZ + BageTH%x%Age + BsexTH%x%Sex + BphuTH%x%PHU + BpartTH%x%partner, name="expThresDZ") CorDZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=.3, lbound=-.99, ubound=.99, name="expCorDZ") # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="expCorMZ", means="M", dimnames=selVars, thresholds="expThresMZ" ) objDZ <- mxExpectationNormal( covariance="expCorDZ", means="M", dimnames=selVars, thresholds="expThresDZ" ) fitFunction <- mxFitFunctionML() # Combine Groups modelMZ <- mxModel( obsAge, obsSex, obsPHU, obspartner, MeanL, betaA, betaS, betaPHU, betapart, Inc, Tmz, ThresMZ, CorMZ, dataMZ, objMZ, fitFunction, name="MZ" ) modelDZ <- mxModel( obsAge, obsSex, obsPHU, obspartner, MeanL, betaA, betaS, betaPHU, betapart, Inc, Tdz, ThresDZ, CorDZ, dataDZ, objDZ, fitFunction, name="DZ" ) minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) Conf <- mxCI (c ('MZ.expCorMZ[2,1]', 'DZ.expCorDZ[2,1]', 'MZ.BageTH[1,1]','MZ.BsexTH[1,1]', 'MZ.BphuTH[1,1]', 'MZ.BpartTH[1,1]') ) SatModel <- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf) # ----------------------------------------------------------------------------------------------- # 1) RUN Saturated Model SatFit <- mxRun(SatModel, intervals=TRUE) (SatSumm <- summary(SatFit))