You are here

Trivariate Model

31 posts / 0 new
Last post
JuanJMV's picture
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

JuanJMV's picture
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.

AdminNeale's picture
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

AdminRobK's picture
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.

AdminRobK's picture
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)
JuanJMV's picture
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!!

JuanJMV's picture
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

AdminRobK's picture
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" )
JuanJMV's picture
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.

AdminRobK's picture
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.

JuanJMV's picture
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!!!

AdminRobK's picture
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()?

JuanJMV's picture
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!!!

AdminRobK's picture
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.

JuanJMV's picture
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

AdminRobK's picture
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)
JuanJMV's picture
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)
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", 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!

AdminRobK's picture
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?

JuanJMV's picture
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

> mxCompare( fit, fitADE )
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.

AdminRobK's picture
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.

AdminRobK's picture
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).

JuanJMV's picture
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).

I had thought change this:
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.

AdminRobK's picture
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]
I had thought change this:
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.

JuanJMV's picture
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!!!

JuanJMV's picture
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.

Thank you so much any insight you could give me about this or materials to learn more would be really helpful

AdminRobK's picture
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).

JuanJMV's picture
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!

AdminRobK's picture
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?

JuanJMV's picture
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

Thanks in advance

AdminRobK's picture
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.

JuanJMV's picture
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!