# Trivariate Model

31 posts / 0 new
Offline
Joined: 07/20/2016 - 13:13
Trivariate Model

Hi Everyone,

I am trying to fit a trivariate model since I think is the best option instead two separately bivariate models but I am not sure.

I have two errors:

1)
Warning messages:

• thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" )

1: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'values' in mxMatrix(type = "Full", nrow = nth, ncol = ntvo, free = TRUE, values = svTh, lbound = lbTh, labels = labTh, name = "thinG")
2: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'lbound' in mxMatrix(type = "Full", nrow = nth, ncol = ntvo, free = TRUE, values = svTh, lbound = lbTh, labels = labTh, name = "thinG")

2)
*fitACE <- mxRun( modelACE, intervals=T )

Error: In model 'twoACEja' the name 'mean_AG_Ln' is used as a free parameter in 'twoACEja.meanG', 'MZ.meanG', and 'DZ.meanG' and as a fixed parameter in 'twoACEja.meanG', 'MZ.meanG', and 'DZ.meanG'

And here the full script

# Select Continuous Variables
varsc     <- c('AG','RU')                   # list of continuous variables names
nvc       <- 2                         # number of continuous variables
ntvc      <- nvc*2                     # number of total continuous variables
conVars   <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="")

# Select Ordinal Variables
nth       <- 1                         # number of thresholds
varso     <- c('SL')                   # list of ordinal variables names
nvo       <- 1                         # number of ordinal variables
ntvo      <- nvo*2                     # number of total ordinal variables
ordVars   <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="")
ordData   <- twinData
#for (i in ordVars) { ordData[,i] <- cut(twinData[,i],
#breaks=quantile(ordData[,i],(0:(nth+1))/(nth+1),na.rm=TRUE), labels=c(0:nth)) }

# Select Variables for Analysis
vars      <- c('AG','RU','SL')              # list of variables names
nv        <- nvc+nvo                   # number of variables
ntv       <- nv*2                      # number of total variables
oc        <- c(0,0,1)                    # 0 for continuous, 1 for ordinal variables
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Covariates for Analysis
#ordData[,'age'] <- ordData[,'age']/100
#ordData   <- ordData[-which(is.na(ordData$age)),] covVars <- c('age', "Sex1" , "Sex2") nc <- 2 # number of covariates # number of covariates # Select Data for Analysis mzData <- subset(ordData, Zyg==1, c(selVars, covVars)) dzData <- subset(ordData, Zyg==2, c(selVars, covVars)) mzDataF <- mzData dzDataF <- dzData mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1)) mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1)) dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1)) dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1)) # Generate Descriptive Statistics apply(mzData,2,myMean) apply(dzData,2,myMean) # Set Starting Values frMV <- c(TRUE,FALSE) # free status for variables frCvD <- diag(frMV,ntv,ntv) # lower bounds for diagonal of covariance matrix frCvD[lower.tri(frCvD)] <- TRUE # lower bounds for below diagonal elements frCvD[upper.tri(frCvD)] <- TRUE # lower bounds for above diagonal elements frCv <- matrix(as.logical(frCvD),4) svMe <- c(1,1,0) # start value for means svPa <- .1 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .1 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix lbPa <- .0001 # start value for lower bounds lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements svLTh <- -1 # start value for first threshold svITh <- 1 # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds # Create Labels labMe <- paste("mean",vars,sep="_") labTh <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_") labBe <- labFull("beta",nc,1) # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12"), name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.Sex1", "data.Sex1","data.Sex1","data.Sex2","data.Sex2", "data.Sex2"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23","b21", "b22", "b23"), name="b2" ) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMe, name="meanG" ) expMean <- mxAlgebra( expression= meanG + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex , name="expMean" ) thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) threG <- mxAlgebra( expression= inc %*% thinG, name="threG" ) # Create Matrices for Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("a",nv), lbound=lbPaD, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), lbound=lbPaD, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" ) # Create Algebra for Variance Comptwonts covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), 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" ) # Create Algebra for Standardization matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Calculate genetic and environmental correlations corA <- mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="rA" ) #cov2cor() corC <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" ) corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" ) # Constrain Variance of Binary Variables matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" ) matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc" ) var1 <- mxConstraint( expression=diag2vec(Oc%&%V)==Unv1, name="Var1" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzDataF, type="raw" ) dataDZ <- mxData( observed=dzDataF, 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(pathB1,pathB2, meanG, thinG, inc, threG, matI, invSD, matUnv, matOc, pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE) defs <- list(defAge,defSex) modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Variance Components rowVC <- rep('VC',nv) colVC <- rep(c('A','C','E','SA','SC','SE'),each=nv) estVC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC)) # Create Confidence Interval Objects ciACE <- mxCI( "VC[1,1]" )# "VC[1,seq(1,3*nv,nv),(2,seq(1,3*nv,nv)),(2,seq(2,3*nv,nv)))]" ) # Build Model with Confidence Intervals modelACE <- mxModel( "twoACEja", pars, var1, modelMZ, modelDZ, multi, estVC, ciACE ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=T ) sumACE <- summary( fitACE ) # Compare with Saturated Model mxCompare( fit, fitACE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEsts(fitACE) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Run AE model modelAE <- mxModel( fitACE, name="twoAEja" ) modelAE <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=T ) mxCompare( fitACE, fitAE ) fitGofs(fitAE) fitEsts(fitAE) # Run CE model modelCE <- mxModel( fitACE, name="twoCEja" ) modelCE <- omxSetParameters( modelCE, labels=labLower("a",nv), free=FALSE, values=0 ) fitCE <- mxRun( modelCE, intervals=T ) mxCompare( fitACE, fitCE ) fitGofs(fitCE) fitEsts(fitCE) # Run E model modelE <- mxModel( fitAE, name="twoEja" ) modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) mxCompare( fitAE, fitE ) fitGofs(fitE) fitEsts(fitE) # Print Comparative Fit Statistics mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ) Thank you so much I have tried everything but I can not make work it Offline Joined: 07/20/2016 - 13:13 I have managed to make it I have managed to make it work but I'm not sure because I have to put 2 thresholds otherwise it does not work, but my ordinal variable is dichotomous. is it correct? Moreover, I am not sure if I have added the covariables properly. Here the full script: # Select Continuous Variables varsc <- c('AG','RU') # list of continuous variables names nvc <- 2 # number of continuous variables ntvc <- nvc*2 # number of total continuous variables conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="") # Select Ordinal Variables nth <- 2 # number of thresholds varso <- c('SL') # list of ordinal variables names nvo <- 1 # number of ordinal variables ntvo <- nvo*2 # number of total ordinal variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") ordData <- twinData # Select Variables for Analysis vars <- c('AG','RU','SL') # list of variables names nv <- nvc+nvo # number of variables ntv <- nv*2 # number of total variables oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis covVars <- c('age', "Sex1" , "Sex2") nc <- 3 # number of covariates # number of covariates # Select Data for Analysis mzData <- subset(ordData, Zyg==1, c(selVars, covVars)) dzData <- subset(ordData, Zyg==2, c(selVars, covVars)) mzDataF <- mzData dzDataF <- dzData mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1)) mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1)) dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1)) dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1)) # Set Starting Values frMV <- c(TRUE,FALSE) # free status for variables frCvD <- diag(frMV,ntv,ntv) # lower bounds for diagonal of covariance matrix frCvD[lower.tri(frCvD)] <- TRUE # lower bounds for below diagonal elements frCvD[upper.tri(frCvD)] <- TRUE # lower bounds for above diagonal elements frCv <- matrix(as.logical(frCvD),4) svMe <- c(1,1,0) # start value for means svPa <- .1 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .1 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix lbPa <- .0001 # start value for lower bounds lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements svLTh <- -1 # start value for first threshold svITh <- 1 # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds # Create Labels labMe <- paste("mean",vars,sep="_") labTh <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_") labBe <- labFull("beta",nc,1) # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12", "b13"), name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.Sex1", "data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23", "b21", "b22", "b23"), name="b2" ) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMe, name="meanG" ) expMean <- mxAlgebra( expression= meanG + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex , name="expMean" ) thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) threG <- mxAlgebra( expression= inc %*% thinG, name="threG" ) # Create Matrices for Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("a",nv), lbound=lbPaD, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), lbound=lbPaD, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" ) # Create Algebra for Variance Comptwonts covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), 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" ) # Create Algebra for Standardization matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Calculate genetic and environmental correlations corA <- mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="rA" ) #cov2cor() corC <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" ) corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" ) # Constrain Variance of Binary Variables matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" ) matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc" ) var1 <- mxConstraint( expression=diag2vec(Oc%&%V)==Unv1, name="Var1" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzDataF, type="raw" ) dataDZ <- mxData( observed=dzDataF, 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(pathB1,pathB2, meanG, thinG, inc, threG, matI, invSD, matUnv, matOc, pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE) defs <- list(defAge,defSex) modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Variance Components rowVC <- rep('VC',nv) colVC <- rep(c('A','C','E','SA','SC','SE'),each=nv) estVC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC)) # Create Confidence Interval Objects ciACE <- mxCI(c("VC[1,1]", "MZ.rA", "MZ.rC", "MZ.rE")) # Build Model with Confidence Intervals modelACE <- mxModel( "twoACEja", pars, var1, modelMZ, modelDZ, multi, estVC, ciACE ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE Model fitACE <- mxTryHard( modelACE, intervals=T, extraTries = 31 ) sumACE <- summary( fitACE ) # Compare with Saturated Model mxCompare( fit, fitACE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEsts(fitACE) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Run AE model modelAE <- mxModel( fitACE, name="twoAEja" ) modelAE <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0 ) fitAE <- mxTryHard( modelAE, intervals=T, extraTries = 31 ) mxCompare( fitACE, fitAE ) fitGofs(fitAE) fitEsts(fitAE) # Run CE model modelCE <- mxModel( fitACE, name="twoCEja" ) modelCE <- omxSetParameters( modelCE, labels=labLower("a",nv), free=FALSE, values=0 ) fitCE <- mxTryHard( modelCE, intervals=T, extraTries = 31 ) mxCompare( fitACE, fitCE ) fitGofs(fitCE) fitEsts(fitCE) # Run E model modelE <- mxModel( fitAE, name="twoEja" ) modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) mxCompare( fitAE, fitE ) fitGofs(fitE) fitEsts(fitE) # Print Comparative Fit Statistics mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ) Thank you so much again. Offline Joined: 03/01/2013 - 14:09 Only one threshold needed, check other issues Hi Juan You may be on the right track, but if the variables are binary then you should need at most one threshold in the threshold matrix. nth <- 1 I couldn't run your script because the variable names don't agree with those in the dataset. Also labFull function is missing. HTH Offline Joined: 01/24/2014 - 12:15 labFull() is from Hermine's labFull() is from Hermine's miFunctions.R . I also assume that Juan is using his own data, and not the 'twinData' object that comes with OpenMx. Offline Joined: 01/24/2014 - 12:15 ordinal variable The covariates look OK to me. If the ordinal variable is dichotomous, you should only need 1 threshold, which allows you to simplify your script a bit. You also need to either fix the one threshold, or the regression intercept for the dichotomous variable, in order to identify the model. I'd say the simplest modification to make would be to try changing thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" ) to thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, lbound=lbTh, labels=labTh, name="thinG" ) , changing expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars ) to expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars ) , and changing pars <- list(pathB1,pathB2, meanG, thinG, inc, threG, matI, invSD, matUnv, matOc, pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE) to pars <- list(pathB1,pathB2, meanG, thinG, matI, invSD, matUnv, matOc, pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE) Offline Joined: 07/20/2016 - 13:13 Thank you! I have made the Thank you! I have made the changes but I have a warning message in this line: thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, lbound=lbTh, labels=labTh, name="thinG" ) Warning message: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'lbound' in mxMatrix(type = "Full", nrow = 1, ncol = ntvo, free = FALSE, values = 0, lbound = lbTh, labels = labTh, name = "thinG") Thank you so much!! Offline Joined: 07/20/2016 - 13:13 Sorry I forgot to say that I Sorry I forgot to say that I use my own data Thanks Offline Joined: 01/24/2014 - 12:15 warning message Warning message: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'lbound' in mxMatrix(type = "Full", nrow = 1, ncol = ntvo, free = FALSE, values = 0, lbound = lbTh, labels = labTh, name = "thinG") That should be benign, since the threshold is fixed. In fact, because it's fixed, you don't need to give it labels or a lower bound. If you want, you could instead do thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinG" ) Offline Joined: 07/20/2016 - 13:13 Thank you so much it works Thank you so much it works perfect! Here you can see my results ACE A C E --------------------rA rC rE-----------Rph SL 0.34 0.35 0.31 AG 0.50 0.12 0.38 RB 0.48 0.22 0.30 SL&AG ---------------------------0.71(0.01,0.99) 0.84(0.121,1) -0.26(-0.56,0.27)------------- 0.38 SL&RB ---------------------------0.12(-0.45,0.67) 0.93(0.21,1) 0.05(-0.39,0.40)----------------0.32 AG&RB ---------------------------0.78(0.59,0.87) 0.98(0.77,1) 0.35(0.24,0.48)-----------------0.66 AE A E---------------rA rE-------------Rph SL 0.79 0.21 AG 0.62 0.38 RB 0.71 0.29 SL&AG-------------------0.61(0.46,0.75) -0.21(-0.27,0.07)-------0.37 SL&RB-------------------0.45(0.33,0.56) -0.08(-0.33,0.18)-------0.31 AG&RB-------------------0.81(0.77,0.86) 0.34(0.27,0.42)--------0.66 This is the comparison base comparison ep minus2LL df AIC diffLL diffdf p 1 twoACEja 27 11475.559 6073 -670.44092 NA NA NA 2 twoACEja twoAEja 21 11487.116 6079 -670.88425 11.556674 6 0.072621853 I am not sure what model should I chose and the criterion to choose one or the other Also, I need to know if A, C, E are significant in each model. I mean for example in the ACE model for sl C is 0.35 and for ag 0.12 I would like to know if C is significant for SL and AG or only for one of them. I think E is always significant because contain error. I have thought calculate CI and if it contais 0 is not significative but I am not sure. Finally, Can I compare the rA correlations in each trait (with statistic parameters)? I mean if the rA between 1-2 is bigger than 1-3… Thank you so much in advance I really appreciate it. Offline Joined: 01/24/2014 - 12:15 suggestions I am not sure what model should I chose and the criterion to choose one or the other Assuming you only want to report the results from a single "best" model, then I recommend selecting the model with the smallest AIC, and reporting confidence intervals for the quantities of interest. BTW, have you fit a saturated model? Is a CE model plausible for your data? However, I advocate reporting conclusions based on multiple models; see, e.g., my relevant publication for more information. Also, I need to know if A, C, E are significant in each model. I mean for example in the ACE model for sl C is 0.35 and for ag 0.12 I would like to know if C is significant for SL and AG or only for one of them. So just add the needed character strings to the reference argument to mxCI(). For instance, ciACE <- mxCI(c("VC[1,1]", "MZ.rA", "MZ.rC", "MZ.rE","MZ.A","MZ.C","MZ.E")) will request CIs for the elements of the additive-genetic, shared-environmental, and nonshared-environmental covariance matrices, in addition to that for which your script already requests CIs. I think E is always significant because contain error. It may or may not be, but inferences about E are usually not of interest. I have thought calculate CI and if it contais 0 is not significative but I am not sure. Correct. Finally, Can I compare the rA correlations in each trait (with statistic parameters)? I mean if the rA between 1-2 is bigger than 1-3… Yes. Create a suitable MxAlgebra and request confidence intervals for it, for instance, rAcompare <- mxAlgebra(cbind(rA[1,2] - rA[1,3], rA[1,2] - rA[2,3], rA[1,3] - rA[2,3]), name="rAcompare") , place it inside the MZ and DZ submodels, and add "MZ.rAcompare" to the vector of character strings provided for argument reference to mxCI() when creating object ciACE. This will give you confidence intervals for the differences in genetic correlations. Offline Joined: 07/20/2016 - 13:13 Once again thank you so much. Once again thank you so much. I have been working with the saturated model but I cannot make it works. I have tried all possible starting values but It doesn't find a solution. I don't know if I am doing something wrong here is the script.: # Select Continuous Variables varsc <- c('AG','RU') # list of continuous variables names nvc <- 2 # number of continuous variables ntvc <- nvc*2 # number of total continuous variables conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="") # Select Ordinal Variables nth <- 1 # number of thresholds varso <- c('SL') # list of ordinal variables names nvo <- 1 # number of ordinal variables ntvo <- nvo*2 # number of total ordinal variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") ordData <- twinData # Select Variables for Analysis vars <- c('AG','RU','SL') # list of variables names nv <- nvc+nvo # number of variables ntv <- nv*2 # number of total variables oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis covVars <- c('age', "Sex1" , "Sex2") nc <- 3 # number of covariates # number of covariates # Select Data for Analysis mzData <- subset(ordData, Zyg==1, c(selVars, covVars)) dzData <- subset(ordData, Zyg==2, c(selVars, covVars)) mzDataF <- mzData dzDataF <- dzData mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1)) mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1)) dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1)) dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1)) # Generate Descriptive Statistics apply(mzData,2,myMean) apply(dzData,2,myMean) #?apply(mzData,2,myCor) #?apply(dzData,2,myCor) # Set Starting Values frMV <- c(TRUE,FALSE) # free status for variables frCvD <- diag(frMV,ntv,ntv) # lower bounds for diagonal of covariance matrix frCvD[lower.tri(frCvD)] <- TRUE # lower bounds for below diagonal elements frCvD[upper.tri(frCvD)] <- TRUE # lower bounds for above diagonal elements frCv <- matrix(as.logical(frCvD),4) svMe <- c(1,1,0) # start value for means svVa <- c(0.5,0.5,0.5) # start value for variance lbVa <- .01 # start value for lower bounds svVas <- diag(svVa,ntv,ntv) # assign start values to diagonal of matrix lbVas <- diag(lbVa,ntv,ntv) # assign lower bounds to diagonal of matrix svLTh <- 0 # start value for first threshold svITh <- 1 # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds # Create Labels labMeMZ <- paste("meanMZ",selVars,sep="_") labMeDZ <- paste("meanDZ",selVars,sep="_") labMeZ <- paste("meanZ",selVars,sep="_") labCvMZ <- labLower("covMZ",ntv) labCvDZ <- labLower("covDZ",ntv) labCvZ <- labLower("covZ",ntv) labVaMZ <- labDiag("covMZ",ntv) labVaDZ <- labDiag("covDZ",ntv) labVaZ <- labDiag("covZ",ntv) labThMZ <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="_") labThDZ <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="_") labThZ <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="_") labBe <- labFull("beta",nc,1) # ------------------------------------------------------------------------------ # PREPARE MODEL # Saturated Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12", "b13"), name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.Sex1", "data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23", "b21", "b22", "b23"), name="b2" ) # Create Algebra for expected Mean & Threshold Matrices meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeMZ, name="meanMZ" ) meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZ, name="meanDZ" ) expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanMZ" ) expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanDZ" ) thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ" ) thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) threMZ <- mxAlgebra( expression= inc %*% thinMZ, name="threMZ" ) threDZ <- mxAlgebra( expression= inc %*% thinDZ, name="threDZ" ) # Create Algebra for expected Covariance Matrices covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" ) covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzDataF, type="raw" ) dataDZ <- mxData( observed=dzDataF, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="covMZ", means="expMeanMZ", dimnames=selVars, thresholds="thinMZ", threshnames=ordVars ) expDZ <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars, thresholds="thinDZ", threshnames=ordVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(pathB1,pathB2, inc ) defs <- list(defAge, defSex) modelMZ <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, thinMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, thinDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Confidence Interval Objects ciCor <- mxCI( c('MZ.covMZ','DZ.covDZ' )) #ciThre <- mxCI( c('MZ.threMZ','DZ.threDZ' )) # Build Saturated Model with Confidence Intervals model <- mxModel( "twoSATja", pars, modelMZ, modelDZ, multi, ciCor ) # ------------------------------------------------------------------------------ # RUN MODEL # Run Saturated Model fit <- mxTryHard( model, intervals=F, extraTries = 31 ) sum <- summary( fit ) round(fit$output$estimate,4) fitGofs(fit) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Test Significance of Covariates modelCov <- mxModel( fit, name="testCov" ) modelCov <- omxSetParameters( modelCov, label=labBe, free=FALSE, values=0 ) fitCov <- mxRun( modelCov ) mxCompare( fit, fitCov ) # Constrain expected Thresholds to be equal across twin order modelEMTVO <- mxModel( fit, name="eqThresTwin" ) svMe <- c(5,0); svVa <- c(2.5,1); svLTh <- 7; svITh <- 1; for (i in 1:nv) { modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labMeMZ[nv+i],labMeMZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeMZ[i] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labMeDZ[nv+i],labMeDZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeDZ[i] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labVaMZ[nv+i],labVaMZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaMZ[i] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labVaDZ[nv+i],labVaDZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaDZ[i] ) } modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThMZ[nvo*nth+1],labThMZ[1]), free=TRUE, values=svLTh, newlabels=labThMZ[1] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThDZ[nvo*nth+1],labThDZ[1]), free=TRUE, values=svLTh, newlabels=labThDZ[1] ) for (i in 2:nth) {for (j in seq(i,nvo*nth,nth)) { modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThMZ[nvo*nth+j],labThMZ[j]), free=TRUE, values=svITh, newlabels=labThMZ[j] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThDZ[nvo*nth+j],labThDZ[j]), free=TRUE, values=svITh, newlabels=labThDZ[j] ) }} fitEMTVO <- mxRun( modelEMTVO, intervals=F ) mxCompare( fit, subs <- list(fitCov, fitEMTVO) ) round(fitEMTVO$output$estimate,4) fitGofs(fitEMTVO) # Constrain expected Thresholds to be equal across twin order and zygosity modelEMTVZ <- mxModel( fitEMTVO, name="eqThresZyg" ) for (i in 1:nv) { modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labMeMZ[i],labMeDZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeZ[i] ) modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labVaMZ[i],labVaDZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaZ[i] ) } modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labThMZ[1],labThDZ[1]), free=TRUE, values=svLTh, newlabels=labThZ[1] ) for (i in 2:nth) {for (j in seq(i,nvo*nth,nth)) { modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labThMZ[j],labThDZ[j]), free=TRUE, values=svITh, newlabels=labThZ[j] ) }} fitEMTVZ <- mxRun( modelEMTVZ, intervals=F ) mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) ) round(fitEMTVZ$output$estimate,4) fitGofs(fitEMTVZ) Thanks!!! Offline Joined: 01/24/2014 - 12:15 specifically? I have been working with the saturated model but I cannot make it works. I have tried all possible starting values but It doesn't find a solution. So, what happens when you try to run it? Also, what do you get from mxVersion()? Offline Joined: 07/20/2016 - 13:13 This is the output: This is the output: Begin fit attempt 31 of at maximum 32 tries Lowest minimum so far: 12836.907293975 OpenMx status code 6 not in list of acceptable status codes, 0 OpenMx status code 6 not in list of acceptable status codes, 0 Not all eigenvalues of Hessian are greater than 0: 787961.505795225, 241712.973930221, 90801.8260774038, 37461.8509249169, 36108.3570058396, 20452.956895597, 16278.6263111995, 13514.5169971915, 12348.4836755664, 7865.0595383285, 7434.41254730365, 7351.21427400012, 6728.35738212829, 5648.20553748093, 4416.33288252297, 3674.66242643252, 3514.28595340669, 3224.37985610757, 3058.62507833283, 2824.51753405145, 2653.72023767685, 2290.52727017331, 2143.71360085118, 1845.17794007907, 1758.30696016368, 1351.05525266071, 1226.13780443413, 1173.91950744854, 867.931406722967, 807.202050148447, 615.874189695127, 502.241971723273, 418.459657188584, 364.354936896746, 36.5272816693136, 30.8297966297892, 25.5256524187929, 16.5921347221574, 12.4923808408779, 11.2477690283236, 5.48989588395323, 5.08592627903294, 4.88676405380493, 2.5967232615173, 2.08173667306169, 1.82647729161902, 1.52624731830076e-05, -8.33106327891148e-05 Begin fit attempt 32 of at maximum 32 tries Fit attempt worse than current best: 12836.908910494 vs 12836.907293975 Retry limit reached Start values from best fit: 0.10631708710754,0.0362139006101167,-0.0613654041393264,0.00569082513418116,-0.185581069609575,-0.344162428087438,1.06767187921709,-8.36632390547581,0.942344815595082,1.07312588323597,0.408909142438524,2.2005873612565,1.30673456966261,41.1843125633866,0.313146221049428,0.167776661842215,0.85467485365846,0.242455084108765,0.292201400080796,0.869137784815067,0.220708232633343,0.521510782740128,0.274754384528732,0.118846729885699,3.69004028161039,0.0849444409204678,0.103993475707172,1.09385843177197,-17.2859999506245,0.966720410717319,1.26641979849299,0.441386464784935,3.94902382148772,1.9668200201659,185.096253954615,0.164428802392192,0.108368824726336,1.08268432674425,0.219192741624143,0.173675315487782,2.13781627870702,0.227229733482567,0.578122306941328,0.197451965247457,0.0792066276546675,3.46936267863339,0.0862225357533247,0.088883979511605 And here my version: OpenMx version: 2.7.9 [GIT v2.7.9] R version: R version 3.3.2 (2016-10-31) Platform: x86_64-w64-mingw32 Default optimiser: NPSOL Thank you so much!!! Offline Joined: 01/24/2014 - 12:15 Because your ordinal variable Because your ordinal variable is dichotomous (has 2 levels), you need to fix its variance, and fix either its threshold or its regression intercept. Your script fixes the threshold, which is OK, but then the intercept needs to be free. Bearing in mind the value of selVars, look at meanMZ$free
meanDZ$free Your script only fixes SL's intercept for twin #2, and you fix the intercepts for twin #1's RU and twin #2's AG, both of which should be free. [Edit: to clarify, all the intercepts should be free if you're going to fix the threshold.] A similar situation holds for the diagonal elements of covMZ and covDZ: only the SL variances should be fixed (and most people fix them to 1.0, but that's arbitrary). Aside from that, are the start values you're using reasonable? You already have estimates of the regression coefficients, and of the MZ and DZ covariance-matrix elements, from your other models. You could use those as start values. You could also run cov() on your MZ and DZ data matrices, with argument use="pair". Since you have an ordinal variable, consider using mxTryHardOrdinal(). Also, be advised that parameterizing covariance matrices as a symmetric matrix of free parameters doesn't always work very well with with the TryHard functions, since the matrices are not guaranteed to be positive-definite at the start of an attempt. Offline Joined: 07/20/2016 - 13:13 Thanks!!! Thanks!!! I have made the changes but I am not sure If I have made the changes properly. Now I think the saturated model is OK here the output. Solution found Running final fit, for Hessian and/or standard errors and/or confidence intervals Warning messages generated from final fit for final fit for Hessian/SEs/CIs Start values from best fit: -0.0691057959162112,-0.0521420158570258,0.0384309413852813,-0.240706583661364,-0.298725394917852,-0.173318222954087,2.73579579403297,1.79440658210804,11.4497520581424,2.69020253759794,1.75369578879318,-0.782094367395943,1.04792129978035,0.40055361067745,2.73025605742978,1.60148107506046,61.4186624152812,0.303395019112009,0.163145240710836,1.15259082182226,0.230520829429149,0.288949826490561,1.08859911323993,0.215686628575476,0.517151654758388,0.332738963798313,0.130010886719768,4.50785881955327,0.123359536402215,0.110825292980583,24.1352589955029,0.478348813364634,2.97100594513919,1.95886486626441,-0.402130467155594,2.9284481924813,1.92691918020341,-0.311361676981213,1.25042143753879,0.436052231591234,0.254634500568326,0.120286524326224,0.663430914290249,0.156510657389065,0.103573383603742,0.0785470263867377,0.225754745967759,0.175168905853076,0.128880192532778,0.230800499171072,0.582082242296643,0.213125789106686,0.0840996230298986,0.191733399855721,0.0926517086525632,0.08013306011578... <truncated> But in the Test Significance of Covariates I have this warning: Running testCov with 58 parameters Warning message: In model 'testCov' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) > mxCompare( fit, fitCov ) base comparison ep minus2LL df AIC diffLL diffdf p 1 twoSATja 58 12559.294 6041 477.29376 NA NA NA 2 twoSATja testCov 58 12559.294 6041 477.29353 -0.00023061521 0 NA Warning message: In collectStatistics(nextBaseSummary, nextCompareSummary) : Model 'twoSATja' has the same degrees of freedom as testCov which means the models are not nested and should not be compared with the likelihood ratio test. Compare these models using the information criteria statistics > And in the # Constrain expected Thresholds to be equal across twin order Error in omxSetParameters(modelEMTVO, label = c(labThMZ[nvo * nth + 1], : The following labels are not present in the model (use 'strict' = FALSE to ignore): 't1MZ_SLEEP2' and 't1MZ_SLEEP1' The full script here: # Select Continuous Variables varsc <- c('AG','RU') # list of continuous variables names nvc <- 2 # number of continuous variables ntvc <- nvc*2 # number of total continuous variables conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="") # Select Ordinal Variables nth <- 1 # number of thresholds varso <- c('SL') # list of ordinal variables names nvo <- 1 # number of ordinal variables ntvo <- nvo*2 # number of total ordinal variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") ordData <- twinData # Select Variables for Analysis vars <- c('AG','RU','SL') # list of variables names nv <- nvc+nvo # number of variables ntv <- nv*2 # number of total variables oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis covVars <- c('age', "Sex1" , "Sex2") nc <- 3 # number of covariates # number of covariates # Select Data for Analysis mzData <- subset(ordData, Zyg==1, c(selVars, covVars)) dzData <- subset(ordData, Zyg==2, c(selVars, covVars)) mzDataF <- mzData dzDataF <- dzData mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1)) mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1)) dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1)) dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1)) # Generate Descriptive Statistics apply(mzData,2,myMean) apply(dzData,2,myMean) # Set Starting Values frMV <- c(TRUE,FALSE) # free status for variables frCvD <- diag(frMV,ntv,ntv) # lower bounds for diagonal of covariance matrix frCvD[lower.tri(frCvD)] <- TRUE # lower bounds for below diagonal elements frCvD[upper.tri(frCvD)] <- TRUE # lower bounds for above diagonal elements frCv <- matrix(as.logical(frCvD),4) svMe <- c(1,1,0) # start value for means svVa <- c(0.5,0.5,0.5) # start value for variance lbVa <- .01 # start value for lower bounds svVas <- diag(svVa,ntv,ntv) # assign start values to diagonal of matrix lbVas <- diag(lbVa,ntv,ntv) # assign lower bounds to diagonal of matrix svLTh <- 0 # start value for first threshold svITh <- 1 # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds # Create Labels labMeMZ <- paste("meanMZ",selVars,sep="_") labMeDZ <- paste("meanDZ",selVars,sep="_") labMeZ <- paste("meanZ",selVars,sep="_") labCvMZ <- labLower("covMZ",ntv) labCvDZ <- labLower("covDZ",ntv) labCvZ <- labLower("covZ",ntv) labVaMZ <- labDiag("covMZ",ntv) labVaDZ <- labDiag("covDZ",ntv) labVaZ <- labDiag("covZ",ntv) labThMZ <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="_") labThDZ <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="_") labThZ <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="_") labBe <- labFull("beta",nc,1) # ------------------------------------------------------------------------------ # PREPARE MODEL # Saturated Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12", "b13"), name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.Sex1", "data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23", "b21", "b22", "b23"), name="b2" ) # Create Algebra for expected Mean & Threshold Matrices meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZ, name="meanMZ" ) meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZ, name="meanDZ" ) expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanMZ" ) expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanDZ" ) thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ" ) thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) threMZ <- mxAlgebra( expression= inc %*% thinMZ, name="threMZ" ) threDZ <- mxAlgebra( expression= inc %*% thinDZ, name="threDZ" ) # Create Algebra for expected Covariance Matrices covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" ) covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzDataF, type="raw" ) dataDZ <- mxData( observed=dzDataF, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="covMZ", means="expMeanMZ", dimnames=selVars, thresholds="thinMZ", threshnames=ordVars ) expDZ <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars, thresholds="thinDZ", threshnames=ordVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(pathB1,pathB2, inc ) defs <- list(defAge, defSex) modelMZ <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, thinMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, thinDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Confidence Interval Objects ciCor <- mxCI( c('MZ.covMZ','DZ.covDZ' )) #ciThre <- mxCI( c('MZ.threMZ','DZ.threDZ' )) # Build Saturated Model with Confidence Intervals model <- mxModel( "twoSATja", pars, modelMZ, modelDZ, multi, ciCor ) # ------------------------------------------------------------------------------ # RUN MODEL # Run Saturated Model fit <- mxTryHardOrdinal( model, intervals=F, extraTries = 31 ) sum <- summary( fit ) round(fit$output$estimate,4) fitGofs(fit) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Test Significance of Covariates modelCov <- mxModel( fit, name="testCov" ) modelCov <- omxSetParameters( modelCov, label=labBe, free=FALSE, values=0 ) fitCov <- mxRun( modelCov ) mxCompare( fit, fitCov ) # Constrain expected Thresholds to be equal across twin order modelEMTVO <- mxModel( fit, name="eqThresTwin" ) svMe <- c(5,0); svVa <- c(2.5,1); svLTh <- 7; svITh <- 1; for (i in 1:nv) { modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labMeMZ[nv+i],labMeMZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeMZ[i] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labMeDZ[nv+i],labMeDZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeDZ[i] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labVaMZ[nv+i],labVaMZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaMZ[i] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labVaDZ[nv+i],labVaDZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaDZ[i] ) } modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThMZ[nvo*nth+1],labThMZ[1]), free=TRUE, values=svLTh, newlabels=labThMZ[1] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThDZ[nvo*nth+1],labThDZ[1]), free=TRUE, values=svLTh, newlabels=labThDZ[1] ) for (i in 2:nth) {for (j in seq(i,nvo*nth,nth)) { modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThMZ[nvo*nth+j],labThMZ[j]), free=TRUE, values=svITh, newlabels=labThMZ[j] ) modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThDZ[nvo*nth+j],labThDZ[j]), free=TRUE, values=svITh, newlabels=labThDZ[j] ) }} fitEMTVO <- mxRun( modelEMTVO, intervals=F ) mxCompare( fit, subs <- list(fitCov, fitEMTVO) ) round(fitEMTVO$output$estimate,4) fitGofs(fitEMTVO) # Constrain expected Thresholds to be equal across twin order and zygosity modelEMTVZ <- mxModel( fitEMTVO, name="eqThresZyg" ) for (i in 1:nv) { modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labMeMZ[i],labMeDZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeZ[i] ) modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labVaMZ[i],labVaDZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaZ[i] ) } modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labThMZ[1],labThDZ[1]), free=TRUE, values=svLTh, newlabels=labThZ[1] ) for (i in 2:nth) {for (j in seq(i,nvo*nth,nth)) { modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labThMZ[j],labThDZ[j]), free=TRUE, values=svITh, newlabels=labThZ[j] ) }} fitEMTVZ <- mxRun( modelEMTVZ, intervals=F ) mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) ) round(fitEMTVZ$output$estimate,4) fitGofs(fitEMTVZ) Thank you again and sorry because I am really new and I can not understand some things Offline Joined: 01/24/2014 - 12:15 suggestions Your script still does not address the issue I raised about the fixed vs. free status of the variances. You need to construct covMZ and covDZ differently. Try this: covfree <- matrix(TRUE,6,6) covfree[3,3] <- FALSE covfree[6,6] <- FALSE covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" ) covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" ) > mxCompare( fit, fitCov ) base comparison ep minus2LL df AIC diffLL diffdf p 1 twoSATja 58 12559.294 6041 477.29376 NA NA NA 2 twoSATja testCov 58 12559.294 6041 477.29353 -0.00023061521 0 NA Warning message: In collectStatistics(nextBaseSummary, nextCompareSummary) : Model 'twoSATja' has the same degrees of freedom as testCov which means the models are not nested and should not be compared with the likelihood ratio test. Compare these models using the information criteria statistics The relevant syntax is: modelCov <- mxModel( fit, name="testCov" ) modelCov <- omxSetParameters( modelCov, label=labBe, free=FALSE, values=0 ) fitCov <- mxRun( modelCov ) mxCompare( fit, fitCov ) I'm surprised it ran without omxSetParameters() throwing an error. Your script never actually uses labBe in the model. All the regression coefficients' labels start with "b", and not "beta_", as in labBe. As a result, the call to omxSetParameters() doesn't actually change anything, and fitCov is just fit, re-run from its solution. That's why it has the same degrees of freedom. Change the call to omxSetParameters() to something like this: modelCov <- omxSetParameters( modelCov, labels=c("b11","b12","b13","b21","b22","b23"), free=FALSE, values=0 ) And in the # Constrain expected Thresholds to be equal across twin order Error in omxSetParameters(modelEMTVO, label = c(labThMZ[nvo * nth + 1], : The following labels are not present in the model (use 'strict' = FALSE to ignore): 't1MZ_SLEEP2' and 't1MZ_SLEEP1' Your script doesn't give any labels to the matrix of thresholds. Try redefining them as: thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") ) thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") ) To constrain thresholds across twin order (I think this will work): modelEMTVO <- omxSetParameters(model=fit, labels=c("tauMZ1","tauMZ2","tauDZ1","tauDZ2"), newlabels=c("tauMZ","tauMZ","tauDZ","tauDZ"), name="eqThreshTwin") modelEMTVO <- omxAssignFirstParameters(modelEMTVO) To further constrain thresholds across zygosity: modelEMTVZ <- omxSetParameters(model=fitEMTVO, labels=c("tauMZ","tauDZ"), newlabels=c("tau","tau"), name="eqThresZyg") modelEMTVZ <- omxAssignFirstParameters(modelEMTVZ) Offline Joined: 07/20/2016 - 13:13 Thank you so much!!! Thank you so much!!! I have made the changes. Also I had to change the optimizer to : mxOption(NULL,"Default optimizer","SLSQP") because with the normal optimizer I always got the same error: Warning message: In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, : The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) I don't know if it is ok. Here you can see the results : base comparison ep minus2LL df AIC diffLL diffdf p 1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA 2 twoSATja testCov 54 11479.071 6045 -610.92935 6.5282698e+01 6 3.7769626e-12 3 twoSATja eqThreshTwin 58 11413.788 6041 -668.21204 4.0857321e-06 2 9.9999796e-01 4 twoSATja eqThresZyg 58 11413.788 6041 -668.21205 -1.0619260e-08 2 1.0000000e+00 And the full script # Select Continuous Variables varsc <- c('AG','RU') # list of continuous variables names nvc <- 2 # number of continuous variables ntvc <- nvc*2 # number of total continuous variables conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="") # Select Ordinal Variables nth <- 1 # number of thresholds varso <- c('SL') # list of ordinal variables names nvo <- 1 # number of ordinal variables ntvo <- nvo*2 # number of total ordinal variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") ordData <- twinData # Select Variables for Analysis vars <- c('AG','RU','SL') # list of variables names nv <- nvc+nvo # number of variables ntv <- nv*2 # number of total variables oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") #ordData <- ordData[-which(is.na(ordData$age)),]
covVars   <- c('age', "Sex1" , "Sex2")
nc        <- 3                         # number of covariates                    # number of covariates

# Select Data for Analysis
mzData    <- subset(ordData, Zyg==1, c(selVars, covVars))
dzData    <- subset(ordData, Zyg==2, c(selVars, covVars))
mzDataF   <- mzData
dzDataF   <- dzData
mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1))
mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1))
dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1))
dzDataF$SL2 <- mxFactor(dzDataF$SLP2, levels =c(0,1))

# Generate Descriptive Statistics
apply(mzData,2,myMean)
apply(dzData,2,myMean)

# Set Starting Values
frMV      <- c(TRUE,FALSE)             # free status for variables
frCvD     <- diag(frMV,ntv,ntv)        # lower bounds for diagonal of covariance matrix
frCvD[lower.tri(frCvD)] <- TRUE        # lower bounds for below diagonal elements
frCvD[upper.tri(frCvD)] <- TRUE        # lower bounds for above diagonal elements
frCv      <- matrix(as.logical(frCvD),4)
svMe      <- c(1,1,0)                   # start value for means
svVa      <- c(0.5,0.5,0.5)                    # start value for variance
lbVa      <- .01                     # start value for lower bounds
svVas     <- diag(svVa,ntv,ntv)        # assign start values to diagonal of matrix
lbVas     <- diag(lbVa,ntv,ntv)        # assign lower bounds to diagonal of matrix
svLTh     <- 0                      # start value for first threshold
svITh     <- 1                         # start value for increments
svTh      <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv)     # start value for thresholds
lbTh      <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv)     # lower bounds for thresholds

# Create Labels
labMeMZ   <- paste("meanMZ",selVars,sep="_")
labMeDZ   <- paste("meanDZ",selVars,sep="_")
labMeZ    <- paste("meanZ",selVars,sep="_")
labCvMZ   <- labLower("covMZ",ntv)
labCvDZ   <- labLower("covDZ",ntv)
labCvZ    <- labLower("covZ",ntv)
labVaMZ   <- labDiag("covMZ",ntv)
labVaZ    <- labDiag("covZ",ntv)

labThMZ   <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="_")
labThDZ   <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="_")
labThZ    <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="_")
labBe     <- labFull("beta",nc,1)

# ------------------------------------------------------------------------------
# PREPARE MODEL

# Saturated Model
# Create Matrices for Covariates and linear Regression Coefficients
defAge    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" )
pathB1     <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12", "b13"), name="b1" )
defSex    <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.Sex1", "data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2"), name="Sex" )
pathB2     <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23", "b21", "b22", "b23"), name="b2" )

# Create Algebra for expected Mean & Threshold Matrices
meanMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZ, name="meanMZ" )
meanDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZ, name="meanDZ" )
expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanDZ" )
thinMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
inc       <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" )
threMZ    <- mxAlgebra( expression= inc %*% thinMZ, name="threMZ" )
threDZ    <- mxAlgebra( expression= inc %*% thinDZ, name="threDZ" )

# Create Algebra for expected Covariance Matrices
covfree <- matrix(TRUE,6,6)
covfree[3,3] <- FALSE
covfree[6,6] <- FALSE
covMZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" )
covDZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" )
# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=mzDataF, type="raw" )
dataDZ    <- mxData( observed=dzDataF, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="covMZ", means="expMeanMZ", dimnames=selVars, thresholds="thinMZ", threshnames=ordVars )
expDZ     <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars, thresholds="thinDZ", threshnames=ordVars )
funML     <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars      <- list(pathB1,pathB2, inc )
defs      <- list(defAge, defSex)
modelMZ   <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, thinMZ,  dataMZ, expMZ, funML )
modelDZ   <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, thinDZ,  dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )

# Create Confidence Interval Objects
ciCor     <- mxCI( c('MZ.covMZ','DZ.covDZ' ))
#ciThre    <- mxCI( c('MZ.threMZ','DZ.threDZ' ))

# Build Saturated Model with Confidence Intervals
model     <- mxModel( "twoSATja", pars, modelMZ, modelDZ, multi, ciCor )

# ------------------------------------------------------------------------------
# RUN MODEL

# Run Saturated Model
fit       <- mxTryHardOrdinal( model, intervals=F, extraTries = 101, finetuneGradient = F,exhaustive = F)
sum       <- summary( fit )
round(fit$output$estimate,4)
fitGofs(fit)

# ------------------------------------------------------------------------------
# RUN SUBMODELS

# Test Significance of Covariates
modelCov  <- mxModel( fit, name="testCov" )
modelCov  <- omxSetParameters( modelCov, labels=c("b11","b12","b13","b21","b22","b23"), free=FALSE, values=0 )
fitCov    <- mxTryHardOrdinal( modelCov, extraTries = 31 )
mxCompare( fit, fitCov )

# Constrain expected Thresholds to be equal across twin order
modelEMTVO  <- mxModel( fit, name="eqThresTwin" )
modelEMTVO <- omxSetParameters(model=fit, labels=c("tauMZ1","tauMZ2","tauDZ1","tauDZ2"), newlabels=c("tauMZ","tauMZ","tauDZ","tauDZ"), name="eqThreshTwin")
modelEMTVO <- omxAssignFirstParameters(modelEMTVO)
fitEMTVO    <- mxTryHardOrdinal( modelEMTVO, intervals=F, extraTries=10 )
mxCompare( fit, subs <- list(fitCov, fitEMTVO) )
round(fitEMTVO$output$estimate,4)
fitGofs(fitEMTVO)

# Constrain expected Thresholds to be equal across twin order and zygosity
modelEMTVZ <-  omxSetParameters(model=fitEMTVO, labels=c("tauMZ","tauDZ"), newlabels=c("tau","tau"), name="eqThresZyg")
modelEMTVZ <- omxAssignFirstParameters(modelEMTVZ)
modelEMTVZ <- mxModel( fitEMTVO, name="eqThresZyg" )
fitEMTVZ    <- mxTryHardOrdinal( modelEMTVZ, intervals=F, extraTries=10)
mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) )
round(fitEMTVZ $output$estimate,4)
fitGofs(fitEMTVZ)

I really appreciate it!

Offline
Joined: 01/24/2014 - 12:15
eqThreshTwin vs. eqThresZyg
Here you can see the results :

base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA
2 twoSATja testCov 54 11479.071 6045 -610.92935 6.5282698e+01 6 3.7769626e-12
3 twoSATja eqThreshTwin 58 11413.788 6041 -668.21204 4.0857321e-06 2 9.9999796e-01
4 twoSATja eqThresZyg 58 11413.788 6041 -668.21205 -1.0619260e-08 2 1.0000000e+00

Something is still not quite right here. The fourth model, "eqThresZyg", should have one less free parameter than the third, "eqThreshTwin". Maybe run omxGetParameters() on the model named "eqThresZyg" before you run it, to check whether my suggested syntax did what it was supposed to?

Offline
Joined: 07/20/2016 - 13:13
You were right I had a

You were right I had a mistake in my script:

This line was before mxModel
modelEMTVZ <- omxSetParameters(model=fitEMTVO, labels=c("tauMZ","tauDZ"), newlabels=c("tau","tau"), name="eqThresZyg")

Thanks!!
Here you can see the correct results saturated model and saturated vs ADE and ACE

> mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) )
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA
2 twoSATja testCov 54 11479.071 6045 -610.92935 6.5282698e+01 6 3.7769626e-12
3 twoSATja eqThreshTwin 58 11413.788 6041 -668.21205 1.7364073e-08 2 9.9999999e-01
4 twoSATja eqThresZyg 57 11413.788 6042 -670.21205 4.1691237e-09 3 1.0000000e+00

>

> mxCompare( fit, fitACE )
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA
2 twoSATja twoACEja 27 11475.559 6073 -670.44092 61.771126 34 0.0024790555

base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA
2 twoSATja twoADEja 27 11474.970 6073 -671.02971 61.182339 34 0.0028790296

And here ADE vs AE and ACE vs AE

base comparison ep minus2LL df AIC diffLL diffdf p
1 twoADEja 27 11474.970 6073 -671.02971 NA NA NA
2 twoADEja twoAEja 21 11487.115 6079 -670.88451 12.145205 6 0.05880676

base comparison    ep  minus2LL  df     AIC        diffLL diffdf            p


1 twoACEja 27 11475.559 6073 -670.44092 NA NA NA
2 twoACEja twoAEja 21 11487.116 6079 -670.88425 11.556674 6 0.072621853
So, what model is the best? Should I explain all as you mentioned before?

Also, I would like to know if it is possible:

1-Can I get CI for phenotypic correlations? I have been trying to get but I don't know how since phenotypic correlations are not in VC I mean I have got SA,SD and SC since they are in VC but I don't know how get the CI's phenotypic correlations since to get rPH y use cov2cor.
2-As you told me in order to know if A, C or D are significant I have calculated CI's for A C and D but for example in the ADE model for one variable D(unstandardized) is 0.05(0,0.15) I mean is weird that I don't have any negative value since I think 0.05 should be non-significant. I have a big sample but I don't know if I am doing something wrong.
3- Also, I would like to know where and how should I get the MZ and DZ correlations in the saturated model? with covariables? or in the ACE or ADE model? I mean I need the cross twin and cross trait correlations

Thank you so much. I really need it and I am learning a lot.

Offline
Joined: 01/24/2014 - 12:15
cov2cor() works in MxAlgebras

cov2cor() works in MxAlgebras, too. You can make MxAlgebras for the correlation matrices, e.g.

corMZ <- mxAlgebra(cov2cor(covMZ), name="corMZ")
corDZ <- mxAlgebra(cov2cor(covDZ), name="corDZ")

and put them into the MZ and DZ submodels, respectively. Then, request confidence intervals for those two algebras; in your case, I think that would entail adding "MZ.corMZ" and "DZ.corDZ" to the character vector provided for argument reference to mxCI(). You could also instead request CIs only for specific elements of the correlation matrices (e.g., "MZ.corMZ[1,2]"), since requesting CIs for all the elements is computationally inefficient, since a correlation matrix is symmetric and has a fixed diagonal. This approach should work for all of the models.

2-As you told me in order to know if A, C or D are significant I have calculated CI's for A C and D but for example in the ADE model for one variable D(unstandardized) is 0.05(0,0.15) I mean is weird that I don't have any negative value since I think 0.05 should be non-significant. I have a big sample but I don't know if I am doing something wrong.

Does "D" here refer to a path coefficient, or to a variance component? If it's a variance component, then you should not expect to see a negative lower confidence limit, since the Cholesky parameterization prevents variance components from being negative (they are sums of squares of diagonal elements of the lower-triangular Cholesky factor). The fact that OpenMx's profile-likelihood CIs respect constraints on the parameters, and can be asymmetric when necessary, is among their advantages.

Offline
Joined: 01/24/2014 - 12:15
Something I just noticed

Something I just noticed about the script in post #17: the dichotomous phenotype has free thresholds and free regression intercepts (in meanMZ and meanDZ). One of those two things has to be fixed to identify the model. The "status RED" warning was probably a result of this unidentification. I suggest just fixing the threshold to zero (which would make constraining it equal across zygosity and twin order pointless).

Offline
Joined: 07/20/2016 - 13:13
Thank you so much!

Thank you so much!
1-OK! but using this I have the correlations for MZ1 and MZ2 I mean this:

mxAlgebra 'corDZ'
$formula: cov2cor(covDZ)$result:
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00000000 0.67553492 0.35730803 0.34876643 0.41805389 0.35158189
[2,] 0.67553492 1.00000000 0.30300303 0.33465014 0.46008228 0.26706220
[3,] 0.35730803 0.30300303 1.00000000 0.25237671 0.30618880 0.38199174
[4,] 0.34876643 0.33465014 0.25237671 1.00000000 0.64877935 0.32860511
[5,] 0.41805389 0.46008228 0.30618880 0.64877935 1.00000000 0.28502398
[6,] 0.35158189 0.26706220 0.38199174 0.32860511 0.28502398 1.00000000
dimnames: NULL

How can I get the total value for MZ and DZ?
2-I am refering to this value:

    A      A      A       D       D       D       E       E       E     SA     SA     SA      SD      SD      SD


VC 0.7930 0.6586 0.5146 0.2176 -0.0449 0.0040 0.5668 0.1360 -0.0552 0.5027 0.8786 1.1104 0.1380 -0.0600 0.0087
VC 0.6586 0.5482 0.4233 -0.0449 0.0471 -0.1371 0.1360 0.2326 -0.0004 0.8786 0.6621 1.4809 -0.0600 0.0569 -0.4796
VC 0.5146 0.4233 0.3468 0.0040 -0.1371 0.4905 -0.0552 -0.0004 0.1627 1.1104 1.4809 0.3468 0.0087 -0.4796 0.4905
SE SE SE
VC 0.3593 0.1814 -0.1190
VC 0.1814 0.2809 -0.0013
VC -0.1190 -0.0013 0.1627

and I calculated its CI in order to know if it is significative am I wrong?

3-Also, I would like to know If I can fit a AE model for one variable and ACE for the other 2 variables. is it possible? I should fix C to 0 for one variable? how can I do this only for one variable?

4- for this:
Something I just noticed about the script in post #17: the dichotomous phenotype has free thresholds and free regression intercepts (in meanMZ and meanDZ). One of those two things has to be fixed to identify the model. The "status RED" warning was probably a result of this unidentification. I suggest just fixing the threshold to zero (which would make constraining it equal across zygosity and twin order pointless).

thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
to :

thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )

But I am not sure if it would be correct.

Thank you so much again in adavance, I really appreciate your help.

Offline
Joined: 01/24/2014 - 12:15
1-OK! but using this I have
1-OK! but using this I have the correlations for MZ1 and MZ2 I mean this:
[snip]
How can I get the total value for MZ and DZ?

Sorry, I don't understand the question?

2-I am refering to this value:
[snip]
and I calculated its CI in order to know if it is significative am I wrong?

That sounds good to me. The element you bolded is the dominance-genetic variance component for the second phenotype, so (given the parameterization you're using) its CI should not have a negative lower limit.

3-Also, I would like to know If I can fit a AE model for one variable and ACE for the other 2 variables. is it possible? I should fix C to 0 for one variable? how can I do this only for one variable?

Yes, that's possible. The easiest way to do it would be to fix to zero the appropriate path coefficients in the common-environmental Cholesky factor (named "c" in your script), using omxSetParameters(), in such a way that the phenotype has no common-environmental variance, nor common-environmental covariance with the other traits. I guess you could instead use MxConstraints to force the phenotype's variance and covariance elements to equal zero in the MxAlgebra named "C", but that is a clumsy approach. Be warned that if you drop the effect of C for only one phenotype, the MxAlgebra named "C" will no longer be a valid (positive-definite) covariance matrix, so the results of other MxAlgebras that use it may be rendered uninterpretable.

4- for this:
[snip]
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
to :

thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )

But I am not sure if it would be correct.

That should work.

Offline
Joined: 07/20/2016 - 13:13
Thanks!

Thanks!

1-Sorry I meant the correlations for MZ and DZ I mean the correlations that we would use in the Falconer's formula

2-So, if its CI should not have a negative lower limit how can I know if D (or other parameter) is significant?

3-I have made the change to this :
modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0)
And It works !! :)
Thank you so much!!!

Offline
Joined: 07/20/2016 - 13:13
Hi again,

Hi again,

I have been working in the saturated model, but if I do this:

thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )

then Taumz1, Taumz2, TauDZ1 and Taudz2 are not in the model therefore the saturated model and eqThresTwin and eqThresZyg models have the same degrees of freedom so what changes should I do in these models?

Also, I have been thinking in de ACE/AE model about what you told me “I guess you could instead use MxConstraints to force the phenotype's variance and covariance elements to equal zero in the MxAlgebra named "C", but that is a clumsy approach. Be warned that if you drop the effect of C for only one phenotype, the MxAlgebra named "C" will no longer be a valid (positive-definite) covariance matrix, so the results of other MxAlgebras that use it may be rendered uninterpretable.” But, I am not sure how I should do a ACE/AE model I thought
modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0) but I'm not sure if it's both practically and theoretically correct

I made this decision because I think C in one trait is not significant 0.35(0,0.57)
But as I mentioned earlier. I'm not sure if I'm calculating well the intervals to know if it is significant or not. In order to do this I calculated the intervals for VC (VC[1,1], [2,2]…) I mean the unstandardized value. I don’t know if it is correct.

Offline
Joined: 01/24/2014 - 12:15
1-Sorry I meant the
1-Sorry I meant the correlations for MZ and DZ I mean the correlations that we would use in the Falconer's formula

?? I'm still confused as to which correlations you might want that aren't in the 'corMZ' and 'corDZ' MxAlgebras.

2-So, if its CI should not have a negative lower limit how can I know if D (or other parameter) is significant?
...
I made this decision because I think C in one trait is not significant 0.35(0,0.57)
But as I mentioned earlier. I'm not sure if I'm calculating well the intervals to know if it is significant or not. In order to do this I calculated the intervals for VC (VC[1,1], [2,2]…) I mean the unstandardized value. I don’t know if it is correct.

If the CI for a quantity contains zero, then that means that you cannot reject the null hypothesis that the quantity equals zero, i.e. it does not differ significantly from zero.

I thought
modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0) but I'm not sure if it's both practically and theoretically correct

If you want to drop C from the third phenotype, then that looks reasonable to me.

I have been working in the saturated model, but if I do this:

thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )

then Taumz1, Taumz2, TauDZ1 and Taudz2 are not in the model therefore the saturated model and eqThresTwin and eqThresZyg models have the same degrees of freedom so what changes should I do in these models?

Are the "eqThresTwin" and "eqThresZyg" models actually of interest? In any event, you could either (1) fix the intercepts of the ordinal phenotype to zero, and free the thresholds, or (2) leave the thresholds fixed, parameterize the ordinal phenotype's intercepts, in the saturated model, in terms of an MZ1 intercept, an MZ2 intercept, a DZ1 intercept, and a DZ2 intercept, and finally fit an analogous pair of additional models that reduce the four free intercept parameters to two (to equate across twin order) and reduce the two free parameters to one (to equate across zygosity).

Offline
Joined: 07/20/2016 - 13:13
Hi!

Hi!
1- Sorry but my English is bad and my knowledge about Open Mx too ☹ If I use

corMZ <- mxAlgebra(cov2cor(covMZ), name="corMZ")
corDZ <- mxAlgebra(cov2cor(covDZ), name="corDZ")

Then I got this matrix:

dimnames: NULL
mxAlgebra 'expCovDZ'
$formula: rbind(cbind(V, cDZ), cbind(t(cDZ), V))$result:
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.60428543 0.76958186 0.48280016 0.59664363 0.41088273 0.40671628
[2,] 0.76958186 0.84383543 0.29644796 0.41088273 0.38962653 0.26024648
[3,] 0.48280016 0.29644796 1.00000003 0.40671628 0.26024648 0.52439786
[4,] 0.59664363 0.41088273 0.40671628 1.60428543 0.76958186 0.48280016
[5,] 0.41088273 0.38962653 0.26024648 0.76958186 0.84383543 0.29644796
[6,] 0.40671628 0.26024648 0.52439786 0.48280016 0.29644796 1.00000003

I would like to know the correlation for MZ and DZ ( rMZ and rDZ) I am not sure if I should get from the saturated model or ACE model and with or without covariables. For example, rMZ could be 0.40 and for DZ 0.20. It should be a heritability about 40%. I don’t know if I am expressing properly

2- So, If I have 0 in one interval the parameter is non-significant, right? I expected for example intervals such as (-0.425, 0.100) I mean not exactly the zero

3- If I do :modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0)
Then the mode does not calculate rC between the other 2 parameters is it correct?

4- I am going to keep on working in the saturated model.
Thank you so much!

Offline
Joined: 01/24/2014 - 12:15
I would like to know the
I would like to know the correlation for MZ and DZ ( rMZ and rDZ) I am not sure if I should get from the saturated model or ACE model and with or without covariables.

I'd say, get it from the saturated model, with covariates.

2- So, If I have 0 in one interval the parameter is non-significant, right? I expected for example intervals such as (-0.425, 0.100) I mean not exactly the zero

Right. If you really want, you could do a likelihood-ratio test by comparing a model in which the quantity of interest is free versus a model in which it is fixed to the null value (which is usually zero). You can usually fix the quantity to the null value via judicious use of omxSetParameters(); occasionally, you might need to use an MxConstraint.

3- If I do :modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0)
Then the mode does not calculate rC between the other 2 parameters is it correct?

I think so? Could you post the output that doesn't look right to you?

Offline
Joined: 07/20/2016 - 13:13
Thank you so much! yes of

Thank you so much! yes of course:

> fitACE$algebras $iSD
mxAlgebra 'iSD'
$formula: solve(sqrt(I * V))$result:
[,1]      [,2]       [,3]
[1,] 0.79311073 0.0000000 0.00000000
[2,] 0.00000000 1.0937125 0.00000000
[3,] 0.00000000 0.0000000 0.99999978
dimnames: NULL

$A mxAlgebra 'A'$formula:  a %*% t(a)
$result: [,1] [,2] [,3] [1,] 0.81389313 0.46688500 0.54976880 [2,] 0.46688500 0.45165068 0.31551369 [3,] 0.54976880 0.31551369 0.77769329 dimnames: NULL$C
mxAlgebra 'C'
$formula: c %*% t(c)$result:
[,1]       [,2] [,3]
[1,] 0.16839824 0.15043721    0
[2,] 0.15043721 0.13439190    0
[3,] 0.00000000 0.00000000    0
dimnames: NULL

$E mxAlgebra 'E'$formula:  e %*% t(e)
$result: [,1] [,2] [,3] [1,] 0.607471487 0.141730178 -0.089903038 [2,] 0.141730178 0.249933037 -0.033470557 [3,] -0.089903038 -0.033470557 0.222307147 dimnames: NULL$V
mxAlgebra 'V'
$formula: A + C + E$result:
[,1]       [,2]       [,3]
[1,] 1.58976285 0.75905239 0.45986576
[2,] 0.75905239 0.83597562 0.28204313
[3,] 0.45986576 0.28204313 1.00000043
dimnames: NULL

$rA mxAlgebra 'rA'$formula:  solve(sqrt(I * A)) %&% A
$result: [,1] [,2] [,3] [1,] 1.00000000 0.77006058 0.69102255 [2,] 0.77006058 1.00000000 0.53236893 [3,] 0.69102255 0.53236893 1.00000000 dimnames: NULL$rC
mxAlgebra 'rC'
$formula: solve(sqrt(I * C)) %&% C$result: (not yet computed) <0 x 0 matrix>
dimnames: NULL

$rE mxAlgebra 'rE'$formula:  solve(sqrt(I * E)) %&% E
\$result:
[,1]        [,2]        [,3]
[1,]  1.00000000  0.36373706 -0.24464408
[2,]  0.36373706  1.00000000 -0.14199544
[3,] -0.24464408 -0.14199544  1.00000000
dimnames: NULL

Offline
Joined: 01/24/2014 - 12:15
looks OK

What's happening is that, internally, sqrt(I * C) isn't invertible, so no result is returned for MxAlgebra 'rC'. Try something like mxEval(cov2cor(C[1:2,1:2]),fitACE,T) if you want to see the shared-environmental correlation for the first two phenotypes.

Offline
Joined: 07/20/2016 - 13:13
It works perfectly as always!

It works perfectly as always!

Thank you so much!! I really appreciate it!