nv <- 3 # number of variables per twin ntv <- nv*2 # number of variables per pair nth <- 6 # number of max thresholds # CREATE LABELS & START VALUES as objects(to ease specification)- 6 th LabThMZM <-c('Tmzm_11','Tmzm_21','Tmzm_31','Tmzm_41','Tmzm_51','Tmzm_61','Mmzm_12','Mmzm_22','Mmzm_32','Mmzm_42', 'Mmzm_52','Mmzm_62','Fmzm_13','Fmzm_23','Fmzm_33','Fmzm_43','Fmzm_53','Fmzm_63') # THs for 3 variables for a twin individual (mzm) LabCorMZM <-c('r21m','r31m','rMZ1m','MZxtm','MZxt2m','r32m','MZxtm','rMZ2m','MZxt3m','MZxt2m','MZxt3m','rMZ3m', 'r21m','r31m','r32m') LabThMZF <-c('Tmzf_11','Tmzf_21','Tmzf_31','Tmzf_41','Tmzf_51','Tmzf_61','Mmzf_12','Mmzf_22','Mmzf_32','Mmzf_42', 'Mmzf_52','Mmzf_62','Fmzf_13','Fmzf_23','Fmzf_33','Fmzf_43','Fmzf_53','Fmzf_63') # THs for 3 variables for a twin individual (mzf) LabCorMZF <-c('r21f','r31f','rMZ1f','MZxtf','MZxt2f','r32f','MZxtf','rMZ2f','MZxt3f','MZxt2f','MZxt3f','rMZ3f', 'r21f','r31f','r32f') LabThDZM <-c('Tdzm_11','Tdzm_21','Tdzm_31','Tdzm_41','Tdzm_51','Tdzm_61','Mdzm_12','Mdzm_22','Mdzm_32','Mdzm_42', 'Mdzm_52','Mdzm_62','Fdzm_13','Fdzm_23','Fdzm_33','Fdzm_43','Fdzm_53','Fdzm_63') # THs for 3 variables for a twin individual (dzm) LabCorDZM <-c('r21m','r31m','rDZ1m','DZxtm','DZxt2m','r32m','DZxtm','rDZ2m','DZxt3m','DZxt2m','DZxt3m','rDZ3m', 'r21m','r31m','r32m') LabThDZF <-c('Tdzf_11','Tdzf_21','Tdzf_31','Tdzf_41','Tdzf_51','Tdzf_61','Mdzf_12','Mdzf_22','Mdzf_32','Mdzf_42', 'Mdzf_52','Mdzf_62','Fdzf_13','Fdzf_23','Fdzf_33','Fdzf_43','Fdzf_53','Fdzf_63') # THs for 3 variables for a twin individual (dzf) LabCorDZF <-c('r21f','r31f','rDZ1f','DZxtf','DZxt2f','r32f','DZxtf','rDZ2f','DZxt3f','DZxt2f','DZxt3f','rDZ3f', 'r21f','r31f','r32f') LabThDZO <-c('Tdzo_11','Tdzo_21','Tdzo_31','Tdzo_41','Tdzo_51','Tdzo_61','Mdzo_12','Mdzo_22','Mdzo_32','Mdzo_42', 'Mdzo_52','Mdzo_62','Fdzo_13','Fdzo_23','Fdzo_33','Fdzo_43','Fdzo_53','Fdzo_63') # THs for 3 variables for a twin individual (dzo) LabCorDZO <-c('r21o','r31o','rDZ1o','DZxto','DZxt2o','r32o','DZxto','rDZ2o','DZxt3o','DZxt2o','DZxt3o','rDZ3o', 'r21o','r31o','r32o') #LabCov <-c('BageThT07', 'BageThT07','BageThT07', 'BageThT07','BageThT07', 'BageThT07', # 'BageThM07', 'BageThM07','BageThM07', 'BageThM07','BageThM07', 'BageThM07', # 'BageThF07', 'BageThF07','BageThF07', 'BageThF07','BageThF07', 'BageThF07') ThPat <-c(T,T,T,T,T,T) StCorMZ <-c(.3,.3,.8,.5,.5,.3,.5,.8,.5,.5,.5,.8,.3,.3,.3) StCorDZ <-c(.3,.3,.4,.2,.2,.3,.2,.4,.2,.2,.2,.4,.3,.3,.3) Thlbound <-c(-3,.001,.001,.001,.001,.001) StTHmzm <-c(.6,rep(.2,nth-1),.9,rep(.2,nth-1),1.2,rep(.2,nth-1)) #th order: twin,mother,father StTHmzf <-c(.6,rep(.2,nth-1),1.0,rep(.2,nth-1),1.4,rep(.2,nth-1)) StTHdzm <-c(.6,rep(.2,nth-1),.9,rep(.2,nth-1),1.5,rep(.2,nth-1)) StTHdzf <-c(.5,rep(.2,nth-1),.9,rep(.2,nth-1),.9,rep(.2,nth-1)) StTHdzo <-c(.7,rep(.2,nth-1),1.2,rep(.2,nth-1),1.4,rep(.2,nth-1)) # 1) Fits a constrained Polychoric correlation model # TH same across twins but different across zyg groups # There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric # ------------------------------------------------------------------------------------------------------------------------------ # Matrix & Algebra for expected means (SND), Thresholds, and correlations Mean <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" ) inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low") Tmzm <-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=ThPat, values=StTHmzm, lbound=Thlbound, ubound=3, labels=LabThMZM, name="ThMZM") ThresMZM <-mxAlgebra( expression= Low%*%ThMZM, name="expThresMZM") CorMZM <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZM, lbound=-.99, ubound=.99, name="expCorMZM") Tmzf <-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=ThPat, values=StTHmzf, lbound=Thlbound, ubound=3, labels=LabThMZF, name="ThMZF") ThresMZF <-mxAlgebra( expression= Low%*%ThMZF, name="expThresMZF") CorMZF <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZF, lbound=-.99, ubound=.99, name="expCorMZF") Tdzm <-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=ThPat, values=StTHdzm, lbound=Thlbound, ubound=3, labels=LabThDZM, name="ThDZM") ThresDZM <-mxAlgebra( expression= Low%*%ThDZM, name="expThresDZM") CorDZM <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZM, lbound=-.99, ubound=.99, name="expCorDZM") Tdzf <-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=ThPat, values=StTHdzf, lbound=Thlbound, ubound=3, labels=LabThDZF, name="ThDZF") ThresDZF <-mxAlgebra( expression= Low%*%ThDZF, name="expThresDZF") CorDZF <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZF, lbound=-.99, ubound=.99, name="expCorDZF") Tdzo <-mxMatrix( type="Full", nrow=nth, ncol=ntv, free=ThPat, values=StTHdzo, lbound=Thlbound, ubound=3, labels=LabThDZO, name="ThDZO") ThresDZO <-mxAlgebra( expression= Low%*%ThDZO, name="expThresDZO") CorDZO <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZO, lbound=-.99, ubound=.99, name="expCorDZO") # Data objects for Multiple Groups dataMZM <- mxData( observed=mzmData, type="raw" ) dataMZF <- mxData( observed=mzfData, type="raw" ) dataDZM <- mxData( observed=dzmData, type="raw" ) dataDZF <- mxData( observed=dzfData, type="raw" ) dataDZO <- mxData( observed=dzoData, type="raw" ) # Define the expectation expFunctionMZM <- mxExpectationNormal( covariance="expCorMZM", means="M", dimnames=selVars, thresholds="expThresMZM" ) expFunctionMZF <- mxExpectationNormal( covariance="expCorMZF", means="M", dimnames=selVars, thresholds="expThresMZF" ) expFunctionDZM <- mxExpectationNormal( covariance="expCorDZM", means="M", dimnames=selVars, thresholds="expThresDZM" ) expFunctionDZF <- mxExpectationNormal( covariance="expCorDZF", means="M", dimnames=selVars, thresholds="expThresDZF" ) expFunctionDZO <- mxExpectationNormal( covariance="expCorDZO", means="M", dimnames=selVars, thresholds="expThresDZO" ) # Choose a fit function fitfunction <- mxFitFunctionML() modelMZM <- mxModel( Mean, Tmzm, inc, ThresMZM, CorMZM, dataMZM, expFunctionMZM, fitfunction, name="MZM" ) modelMZF <- mxModel( Mean, Tmzf, inc, ThresMZF, CorMZF, dataMZF, expFunctionMZF, fitfunction, name="MZF" ) modelDZM <- mxModel( Mean, Tdzm, inc, ThresDZM, CorDZM, dataDZM, expFunctionDZM, fitfunction, name="DZM" ) modelDZF <- mxModel( Mean, Tdzf, inc, ThresDZF, CorDZF, dataDZF, expFunctionDZF, fitfunction, name="DZF" ) modelDZO <- mxModel( Mean, Tdzo, inc, ThresDZO, CorDZO, dataDZO, expFunctionDZO, fitfunction, name="DZO" ) SatModel <- mxModel( "Sat", modelMZM, modelMZF, modelDZM, modelDZF, modelDZO) # ------------------------------------------------------------------------------------------------------------------------------- # 1) RUN Saturated Model SatFit <- mxRun(SatModel, intervals=F) summary(SatFit)