# Select Variables for Analysis selVars <- c('cohenpss_score_A','iv_A','cohenpss_score_B','iv_B') #Selecting these columns for analysis aceVars <- c("A1pss","C1pss","E1pss","A1iv","C1iv","E1iv", "A2pss","C2pss","E2pss","A2iv","C2iv","E2iv") #creating latent variables #"Dummy" latent variables, for nuisance covariates: dummyVars <- c("dummyAge_A","dummySex_A","dummyAge_B","dummySex_B") # Select Data for Analysis: Creating two subsets for MZ and DZ to generate descriptive stats for each group type mzData <- subset(df, zyg_A==1, selVars) dzData <- subset(df, zyg_A==0, selVars) # Generate Descriptive Statistics colMeans(mzData,na.rm=TRUE) #mean ~11 colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") #cov ~45 cov(dzData,use="complete") # Path objects for Multiple Groups manifestVars=selVars latentVars=c(aceVars,dummyVars) # variances of latent variables latVariances <- mxPath( from=aceVars, arrows=2, free=FALSE, values=1 ) # means of latent variables... AND SPECIFYING SINGLE=HEADED ARROWS FROM TRIANGLE (WHICH IS FIXED AT 1) TO EACH OF THE LATENT VARIABLES, FIXED AT 0 latMeans <- mxPath( from="one", to=aceVars, arrows=1, free=FALSE, values=0 ) # means of observed variables obsIntercepts <- mxPath( from="one", to=selVars, arrows=1, free=c(T,F,T,F), values=c(10,0,10,0), labels=c("intrcptPss","intrcptIv","intrcptPss","intrcptIv") ) #The mean cohen_pss_score is around 11 obsSlopes1 <- mxPath( from="one", to=dummyVars, free=F, labels=c("data.Age_A","data.Sex_A","data.Age_B","data.Sex_B") ) obsSlopes2 <- mxPath( from="dummyAge_A", to=c('cohenpss_score_A','iv_A'), free=T, values=0.1, labels=c("beta_pss_age","beta_iv_age") ) obsSlopes3 <- mxPath( from="dummySex_A", to=c('cohenpss_score_A','iv_A'), free=T, values=0.1, labels=c("beta_pss_sex","beta_iv_sex") ) obsSlopes4 <- mxPath( from="dummyAge_B", to=c('cohenpss_score_B','iv_B'), free=T, values=0.1, labels=c("beta_pss_age","beta_iv_age") ) obsSlopes5 <- mxPath( from="dummySex_B", to=c('cohenpss_score_B','iv_B'), free=T, values=0.1, labels=c("beta_pss_sex","beta_iv_sex") ) #Thresholds for ordinal IV: Tau <- mxThreshold( vars=c("iv_A","iv_B"),nThresh=4,free=T,values=c(-1.5,-0.5,0.5,1.5),labels=paste0("tau",1:4)) # path coefficients for twin 1 pathAceT1_pss <- mxPath( from=c("A1pss","C1pss","E1pss"), to="cohenpss_score_A", arrows=1, free=TRUE, values=.5, label=c("apss","cpss","epss") ) pathAceT1_iv <- mxPath( from=c("A1iv","C1iv","E1iv"), to="iv_A", arrows=1, free=TRUE, values=.5, label=c("aiv","civ","eiv") ) # path coefficients for twin 2 pathAceT2_pss <- mxPath( from=c("A2pss","C2pss","E2pss"), to="cohenpss_score_B", arrows=1, free=TRUE, values=.5, label=c("apss","cpss","epss") ) pathAceT2_iv <- mxPath( from=c("A2iv","C2iv","E2iv"), to="iv_B", arrows=1, free=TRUE, values=.5, label=c("aiv","civ","eiv") ) #Identifying constraint: idcons <- mxConstraint( (aiv^2) + (civ^2) + (eiv^2) - 1 == 0, name="identifyingConstraint") #Regression path of pss onto iv: regPaths <- mxPath( from=c("iv_A","iv_B"),to=c('cohenpss_score_A','cohenpss_score_B'), arrows=1, free=T, values=0.1, labels=c("b","b")) #The common environmental factors are the same for both twins (by definition), so setting the fixed correlation of 1 for C1 and C2 # covariance between C1 & C2 covC1C2 <- mxPath( from=c("C1pss","C1iv"), to=c("C2pss","C2iv"), arrows=2, free=FALSE, values=1 ) # covariance between A1 & A2 in MZ twins covA1A2_MZ <- mxPath( from=c("A1pss","A1iv"), to=c("A2pss","A2iv"), arrows=2, free=FALSE, values=1 ) # covariance between A1 & A2 in DZ twins covA1A2_DZ <- mxPath( from=c("A1pss","A1iv"), to=c("A2pss","A2iv"), arrows=2, free=FALSE, values=.5 ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Combine Groups paths <- list( latVariances, latMeans, obsIntercepts, obsSlopes1, obsSlopes2, obsSlopes3, obsSlopes4, obsSlopes5, regPaths, Tau, pathAceT1_pss, pathAceT1_iv, pathAceT2_pss, pathAceT2_iv , covC1C2 ) modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars, latentVars=aceVars, paths, covA1A2_MZ, dataMZ, idcons ) modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars, latentVars=aceVars, paths, covA1A2_DZ, dataDZ ) obj <- mxFitFunctionMultigroup(c("MZ", "DZ")) modelACE <- mxModel(model="ACE", modelMZ, modelDZ, obj ) # Run Model fitACE <- mxRun(modelACE) sumACE <- summary(fitACE)