Login

News

doubt with binary variables

Posted on
No user picture. JuanJMV Joined: 07/20/2016
Forums
Behavioral Genetics Models
Hi everyone,

I am trying to do a joint analysis between 2 binary variables and 1 continuous variable I am having some problems to fix the variance for the 2 binary variables equal 1.

I think I should change something here

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" )

NVO=2
NV=3

However, I have done several changes and it doesn’t work since in some tries the matrix are not conformable and in others it doesn’t fix the variance to 1. How could I fix the variance for the two binary variables to 1?

Thank you so much in advance.
Juan J Madrid-Valero

Replied on Wed, 01/24/2018 - 10:35
Picture of user. AdminRobK Joined: 01/24/2014

Without the context of the full script, I can't really answer your question.
Log in or register to post comments
Replied on Wed, 01/24/2018 - 11:44
No user picture. JuanJMV Joined: 07/20/2016

Thanks Rob for your prompt response here you can find the full script:

# Select Continuous Variables
varsc <- c('SL') # list of continuous variables names
nvc <- 1 # 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('NE',"L") # list of ordinal variables names
nvo <- 2 # 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('SL','NE',"L") # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,1,1) # 0 for continuous, 1 for ordinal variables COMPROBAR
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Covariates for Analysis
covVars <- c('age')
nc <- 3 # number of covariates

# Select Data for Analysis
mzData <- subset(ordData, Zyg==1|Zyg==3, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2|Zyg==4|Zyg==5, c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$NE1 <- mxFactor(mzDataF$NE1, levels =c(1,2))
mzDataF$NE2 <- mxFactor(mzDataF$NE2, levels =c(1,2))
dzDataF$NE1 <- mxFactor(dzDataF$NE1, levels =c(1,2))
dzDataF$NE2 <- mxFactor(dzDataF$NE2, levels =c(1,2))
mzDataF$L1 <- mxFactor(mzDataF$L1, levels =c(1,2))
mzDataF$L2 <- mxFactor(mzDataF$L2, levels =c(1,2))
dzDataF$L1 <- mxFactor(dzDataF$L1, levels =c(1,2))
dzDataF$L2 <- mxFactor(dzDataF$L2, levels =c(1,2))

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

# Set Starting Values
frMV <- c(TRUE,FALSE,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(4,0,0) # start value for means
svPa <- .4 # start value for path coefficient
svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix
svPe <- .8 # 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.5 # start value for first threshold
svITh <- 1 # start value for increments
svTh <- 0 # start value for thresholds
lbTh <- 0 # 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" )
pathB <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=labBe, name="b" )

# 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(b%*%Age),t(b%*%Age)), 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(pathB, meanG, thinG, inc, threG, matI, invSD, matUnv, matOc,
pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE)
defs <- list(defAge)
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=F )
sumACE <- summary( fitACE )

Log in or register to post comments
Replied on Wed, 01/24/2018 - 13:42
Picture of user. AdminRobK Joined: 01/24/2014

At least in the script you posted, Oc%&%V is going to evaluate to a 1x1 matrix, but matUnv is 2x1, and the two sides of the comparator in an MxConstraint have to be of the same dimensions. A simple way to fix the variance of the two binary variables to 1 is to replace those three lines with:

matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
var1 <- mxConstraint(expression=rbind(V[2,2],V[3,3]) == Unv1, name="Var1")

You could have also just used two MxConstraint objects (one for each variance).

Since your analysis involves an MxConstraint, I recommend using SLSQP as the optimizer.

You could also avoid using an MxConstraint entirely if you change the parameterization of the model to direct variance components, and identify the scale of the ordinal variables by fixing their nonshared-environmental variances to 1.

Log in or register to post comments
Replied on Thu, 01/25/2018 - 06:31
No user picture. JuanJMV Joined: 07/20/2016

Thank you so much Rob,

I didn’t know why the matrix was ( expression=diag2vec(Oc%&%V)) 1x1.

Also, I have run the script and the results are pretty feasible but I got this 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 have tried tons of starting values but I have not been able to remove the “problem”
By the way I have used the optimizer that you told me and the function MxtryhardOrdinal with extra tries.

Juan J Madrid-Valero

Here the full message:
Retry limit reached

Solution found

Start values from best fit:
0.0873764505129859,-0.00302648728148331,0.00387006292591793,0.446376724057658,0.765174696854374,1.24795735809318,2.27769791884809,-0.336483766998107,-0.275906439822281,0.416526067571617,0.471557647293468,0.00245202168140257,0.202289271880198,0.352135547711807,0.0718736069163078,0.00290227641677533,0.00372693059844336,0.00118080806141676,3.1768820509539,-0.106534386756017,-0.147029899250398,0.76021560181588,0.40208133001671,0.716267708855051
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)

Log in or register to post comments
Replied on Thu, 01/25/2018 - 10:44
Picture of user. AdminRobK Joined: 01/24/2014

I didn’t know why the matrix was ( expression=diag2vec(Oc%&%V)) 1x1.

Oc is 1x3. Operator %&% is the quadratic-product operator (see the documentation for mxAlgebra()). Thus, Oc%&%V is equivalent to Oc%*%V%*%t(Oc), which is (1x3) times (3x3) times (3x1), which evaluates to a 1x1 product.

Also, I have run the script and the results are pretty feasible but I got this 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 have tried tons of starting values but I have not been able to remove the “problem”
By the way I have used the optimizer that you told me and the function MxtryhardOrdinal with extra tries.

Which version of OpenMx are you using? As of v2.8.3, status Red should occur less frequently in ordinal analyses than it used to. Nonetheless, status Red is sometimes unavoidable in ordinal analyses, in that it can still occur even when the optimizer has found a local minimum. You could try using CSOLNP instead of SLSQP. If that still doesn't help, then it might be worth reparameterizing to eliminate the MxConstraint, which I alluded to in another post.

Log in or register to post comments
Replied on Fri, 01/26/2018 - 06:12
No user picture. JuanJMV Joined: 07/20/2016

Thanks Rob,

I am using R version 3.4.3 and I have tried with the CSOLNP optimizer but it doesn't work either. I am going to try to change the parameterization of the model to direct variance components. I am not sure how I should do it. Could you please give me some guidance or where I could find information about that?

Juan J Madrid-Valero

Log in or register to post comments
Replied on Fri, 01/26/2018 - 11:23
Picture of user. AdminRobK Joined: 01/24/2014

I am using R version 3.4.3

What about your OpenMx version?

I am going to try to change the parameterization of the model to direct variance components. I am not sure how I should do it. Could you please give me some guidance or where I could find information about that?

Let's see...one way, involving minimal changes to your script, would be as follows. First, remove the line that creates pathE. Then, define covE as:

covE <- mxMatrix(
type="Symm",nrow=nv,
free=c(T,T,T,F,T,F),
values=c(0.8,0,0,1,0,1),
labels=c("e_1_1", "e_2_1", "e_3_1", "e_2_2", "e_3_2", "e_3_3"),
lbound=c(0.0001,rep(NA,5)),
name="E")

Then, remove the lines that define matUnv, matOc, and var1. Finally, remove pathE, matUnv, matOc, and var1 from the definition of pars. What these changes do is identify the latent scale of the ordinal variables by fixing their nonshared-environmental variance, rather than their total variance, to 1.0. As the MxConstraint is no longer required, use the CSOLNP optimizer.

While I'm thinking about it--consider removing the lower bounds for pathA and pathE. You can put them back in if the optimizer keeps reaching really crazy solutions.

Log in or register to post comments
Replied on Sat, 01/27/2018 - 10:33
No user picture. JuanJMV Joined: 07/20/2016

Sorry! I am using OpenMx version: 2.8.3 [GIT v2.8.3] ; R version: R version 3.4.3 (2017-11-30).

I have made all the changes and now It works and I don't have a red mx status but I have this message:
Retry limit reached

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.086752181736907,-0.00393385954682701,0.00457737378618429,0.481260848285712,0.526771126020844,1.19685352917708,2.29257177590332,-0.392479607148631,-0.315671611477633,0.782950339036923,0.501766203589714,-0.266923949247171,10.0657799256565,-0.474366335537433,-0.565466421429395,0.490967007947968

the resuts are the same that with the mx status red.

Thank you so much!!!

Log in or register to post comments
Replied on Fri, 02/02/2018 - 10:32
No user picture. JuanJMV Joined: 07/20/2016

Do you think that these results are reliable?

by the way is it possible to know the influences shared by the three phenotypes (all together), not only two by two phenotypes?

and also I would like to know what is the easiest way to compare two correlations in the saturated model

I have tried creating these objects:

CICORMZF <-mxAlgebra(expression =cov2cor(MZf.covMZf), name="CORMZF")
CICORDZF <-mxAlgebra(expression =cov2cor(DZf.covDZf), name="CORDZF")
CICORMZM <-mxAlgebra(expression =cov2cor(MZm.covMZm), name="CORMZM")
CICORDZM <-mxAlgebra(expression =cov2cor(DZm.covDZm), name="CORDZM")
CICORDZO <-mxAlgebra(expression =cov2cor(DZo.covDZo), name="CORDZO")
rAcompare <- mxAlgebra(cbind(CORDZF[1,2] - CORDZM[1,2], name="rAcompare")

and adding them to the model but it doesn't work.

It is a 5group saturated model and I would like to compare the correlations between males and females.
I have thought to equal the correlations in the model using omxsetparameter but in the model I don't have the correlations for the five groups only in the submodels so I am not sure how to proceed.

Thanks!!!

Log in or register to post comments
Replied on Wed, 02/07/2018 - 10:45
Picture of user. AdminRobK Joined: 01/24/2014

Do you think that these results are reliable?

How do you mean?

by the way is it possible to know the influences shared by the three phenotypes (all together), not only two by two phenotypes?

Yes, that can be done. Look into the Common-Pathway and Independent-Pathway models. Both of them do what you're describing, but in different ways.

and also I would like to know what is the easiest way to compare two correlations in the saturated model

I have tried creating these objects:

CICORMZF <-mxAlgebra(expression =cov2cor(MZf.covMZf), name="CORMZF")
CICORDZF <-mxAlgebra(expression =cov2cor(DZf.covDZf), name="CORDZF")
CICORMZM <-mxAlgebra(expression =cov2cor(MZm.covMZm), name="CORMZM")
CICORDZM <-mxAlgebra(expression =cov2cor(DZm.covDZm), name="CORDZM")
CICORDZO <-mxAlgebra(expression =cov2cor(DZo.covDZo), name="CORDZO")
rAcompare <- mxAlgebra(cbind(CORDZF[1,2] - CORDZM[1,2], name="rAcompare")

and adding them to the model but it doesn't work.

It is a 5group saturated model and I would like to compare the correlations between males and females.
I have thought to equal the correlations in the model using omxsetparameter but in the model I don't have the correlations for the five groups only in the submodels so I am not sure how to proceed.

One thing wrong I notice immediately is that there's an unbalanced parenthesis in

rAcompare <- mxAlgebra(cbind(CORDZF[1,2] - CORDZM[1,2], name="rAcompare")

. I don't think you need the column-bind, so try changing it to

rAcompare <- mxAlgebra(CORDZF[1,2] - CORDZM[1,2], name="rAcompare")

. Otherwise, that block on syntax looks OK to me, but it's hard to say without seeing it in the context of the full script, and without any specifics beyond "it doesn't work." Since you're referencing the covariance matrices by "modelname-dot-objectname", I assume you're placing these algebras in the "container" MxModel, right?

Log in or register to post comments
Replied on Wed, 02/07/2018 - 11:18
No user picture. JuanJMV Joined: 07/20/2016

Sorry! When I said if these results are reliable I wanted to say if with this message:
"Retry limit reached

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.086752181736907,-0.00393385954682701,0.00457737378618429,0.481260848285712,0...."

I could be sure that the soltuion is "good" I mean since retry limit have been reached I don't know if I should try other starting values or something else (although there is a solution)

Here the 5 group script I am working with:

# Set Starting Values
svMe <- 0 # start value for means
svVa <- 1 # start value for variance
svVas <- diag(svVa,ntv,ntv) # assign start values to diagonal of matrix
lbVa <- .0001 # start value for lower bounds
lbVas <- diag(lbVa,ntv,ntv) # assign lower bounds values to diagonal of matrix

# Create Labels
labMeMZf <- c("mMZf1","mMZf2") # labels for means for MZ twins
labMeDZf <- c("mDZf1","mDZf2") # labels for means for DZ twins
labMeMZm <- c("mMZm1","mMZm2") # labels for means for MZ twins
labMeDZm <- c("mDZm1","mDZm2") # labels for means for DZ twins
labMeDZo <- c("mDZo1","mDZo2") # labels for means for DZ twins
labCvMZf <- c("vMZf1","cMZf21","vMZf2")# labels for (co)variances for MZ twins
labCvDZf <- c("vDZf1","cDZf21","vDZf2")# labels for (co)variances for DZ twins
labCvMZm <- c("vMZm1","cMZm21","vMZm2")# labels for (co)variances for MZ twins
labCvDZm <- c("vDZm1","cDZm21","vDZm2")# labels for (co)variances for DZ twins
labCvDZo <- c("vDZo1","cDZo21","vDZo2")# labels for (co)variances for DZ twins
labVaMZf <- c("vMZf1","vMZf2") # labels for (co)variances for MZ twins
labVaDZf <- c("vDZf1","vDZf2") # labels for (co)variances for DZ twins
labVaMZm <- c("vMZm1","vMZm2") # labels for (co)variances for MZ twins
labVaDZm <- c("vDZm1","vDZm2") # labels for (co)variances for DZ twins
labVaDZo <- c("vDZo1","vDZo2") # labels for (co)variances for DZ twins

# ------------------------------------------------------------------------------
# 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" )
pathBf <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="bf11", name="bf" )
pathBm <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="bm11", name="bm" )

# Create Algebra for expected Mean Matrices
meanMZf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZf, name="meanMZf" )
meanDZf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZf, name="meanDZf" )
meanMZm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZm, name="meanMZm" )
meanDZm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZm, name="meanDZm" )
meanDZo <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZo, name="meanDZo" )
expMeanMZf<- mxAlgebra( expression= meanMZf + cbind(bf%*%Age,bf%*%Age), name="expMeanMZf" )
expMeanDZf<- mxAlgebra( expression= meanDZf + cbind(bf%*%Age,bf%*%Age), name="expMeanDZf" )
expMeanMZm<- mxAlgebra( expression= meanMZm + cbind(bm%*%Age,bm%*%Age), name="expMeanMZm" )
expMeanDZm<- mxAlgebra( expression= meanDZm + cbind(bm%*%Age,bm%*%Age), name="expMeanDZm" )
expMeanDZo<- mxAlgebra( expression= meanDZo + cbind(bf%*%Age,bm%*%Age), name="expMeanDZo" )

# Create Algebra for expected Variance/Covariance Matrices
covMZf <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvMZf, name="covMZf" )
covDZf <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvDZf, name="covDZf" )
covMZm <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvMZm, name="covMZm" )
covDZm <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvDZm, name="covDZm" )
covDZo <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvDZo, name="covDZo" )

# Create Algebra for Maximum Likelihood Estimates of Twin Correlations
matI <- mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="I" )
corMZf <- mxAlgebra( solve(sqrt(I*covMZf)) %&% covMZf, name="corMZf" )
corDZf <- mxAlgebra( solve(sqrt(I*covDZf)) %&% covDZf, name="corDZf" )
corMZm <- mxAlgebra( solve(sqrt(I*covMZm)) %&% covMZm, name="corMZm" )
corDZm <- mxAlgebra( solve(sqrt(I*covDZm)) %&% covDZm, name="corDZm" )
corDZo <- mxAlgebra( solve(sqrt(I*covDZo)) %&% covDZo, name="corDZo" )

# Create Data Objects for Multiple Groups
dataMZf <- mxData( observed=mzfData, type="raw" )
dataDZf <- mxData( observed=dzfData, type="raw" )
dataMZm <- mxData( observed=mzmData, type="raw" )
dataDZm <- mxData( observed=dzmData, type="raw" )
dataDZo <- mxData( observed=dzoData, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZf <- mxExpectationNormal( covariance="covMZf", means="expMeanMZf", dimnames=selVars )
expDZf <- mxExpectationNormal( covariance="covDZf", means="expMeanDZf", dimnames=selVars )
expMZm <- mxExpectationNormal( covariance="covMZm", means="expMeanMZm", dimnames=selVars )
expDZm <- mxExpectationNormal( covariance="covDZm", means="expMeanDZm", dimnames=selVars )
expDZo <- mxExpectationNormal( covariance="covDZo", means="expMeanDZo", dimnames=selVars )
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars <- list(pathBf, pathBm)
defs <- list(defAge)
modelMZf <- mxModel( name="MZf", pars, defs, meanMZf, expMeanMZf, covMZf, matI, corMZf, dataMZf, expMZf, funML )
modelDZf <- mxModel( name="DZf", pars, defs, meanDZf, expMeanDZf, covDZf, matI, corDZf, dataDZf, expDZf, funML )
modelMZm <- mxModel( name="MZm", pars, defs, meanMZm, expMeanMZm, covMZm, matI, corMZm, dataMZm, expMZm, funML )
modelDZm <- mxModel( name="DZm", pars, defs, meanDZm, expMeanDZm, covDZm, matI, corDZm, dataDZm, expDZm, funML )
modelDZo <- mxModel( name="DZo", pars, defs, meanDZo, expMeanDZo, covDZo, matI, corDZo, dataDZo, expDZo, funML )
multi <- mxFitFunctionMultigroup( c("MZf","DZf","MZm","DZm","DZo") )

# Create Confidence Interval Objects
ciCov <- mxCI( c('MZf.covMZf','DZf.covDZf','MZm.covMZm','DZm.covDZm','DZo.covDZo' ))
ciMean <- mxCI( c('MZf.expMeanMZf','DZf.expMeanDZf','MZm.expMeanMZm','DZm.expMeanDZm','DZo.expMeanDZo' ))
CICORMZF <-mxAlgebra(expression =cov2cor(MZf.covMZf), name="CORMZF")
CICORDZF <-mxAlgebra(expression =cov2cor(DZf.covDZf), name="CORDZF")
CICORMZM <-mxAlgebra(expression =cov2cor(MZm.covMZm), name="CORMZM")
CICORDZM <-mxAlgebra(expression =cov2cor(DZm.covDZm), name="CORDZM")
CICORDZO <-mxAlgebra(expression =cov2cor(DZo.covDZo), name="CORDZO")
ci2 <- mxCI( c('CORMZF','CORDZF',"CORMZM","CORDZM","CORDZO") )

# Build Saturated Model with Confidence Intervals
model <- mxModel( "oneSAT5ca", pars, modelMZf, modelDZf, modelMZm, modelDZm, modelDZo, multi, ciCov, ciMean,CICORMZF,CICORDZF,CICORMZM,CICORDZM,CICORDZO,ci2)

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

# Run Saturated Mode
fit <- mxRun( model, intervals=T )
sum <- summary( fit )
round(fit$output$estimate,4)
fitGofs(fit)

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

# Test Significance of Covariate
modelCov <- mxModel( fit, name="testCov" )
modelCov <- omxSetParameters( modelCov, label=c("bf11","bm11"), free=FALSE, values=0 )
fitCov <- mxRun( modelCov )
mxCompare( fit, fitCov )

# Constrain expected Means to be equal across twin order
modelEMO <- mxModel( fit, name="eqMeansTwin" )
modelEMO <- omxSetParameters( modelEMO, label=labMeMZf, free=TRUE, values=svMe, newlabels='mMZf' )
modelEMO <- omxSetParameters( modelEMO, label=labMeDZf, free=TRUE, values=svMe, newlabels='mDZf' )
modelEMO <- omxSetParameters( modelEMO, label=labMeMZm, free=TRUE, values=svMe, newlabels='mMZm' )
modelEMO <- omxSetParameters( modelEMO, label=labMeDZm, free=TRUE, values=svMe, newlabels='mDZm' )
fitEMO <- mxRun( modelEMO, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMO) )
round(fitEMO$output$estimate,4)
fitGofs(fitEMO)

# Constrain expected Means and Variances to be equal across twin order
modelEMVO <- mxModel( fitEMO, name="eqMVarsTwin" )
modelEMVO <- omxSetParameters( modelEMVO, label=labVaMZf, free=TRUE, values=svVa, newlabels='vMZf' )
modelEMVO <- omxSetParameters( modelEMVO, label=labVaDZf, free=TRUE, values=svVa, newlabels='vDZf' )
modelEMVO <- omxSetParameters( modelEMVO, label=labVaMZm, free=TRUE, values=svVa, newlabels='vMZm' )
modelEMVO <- omxSetParameters( modelEMVO, label=labVaDZm, free=TRUE, values=svVa, newlabels='vDZm' )
fitEMVO <- mxRun( modelEMVO, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMO, fitEMVO) )
round(fitEMVO$output$estimate,4)
fitGofs(fitEMVO)

# Constrain expected Means and Variances to be equal across twin order and zygosity
modelEMVZ <- mxModel( fitEMVO, name="eqMVarsZyg" )
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("mMZf","mDZf"), free=TRUE, values=svMe, newlabels='mZf' )
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("vMZf","vDZf"), free=TRUE, values=svVa, newlabels='vZf' )
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("mMZm","mDZm"), free=TRUE, values=svMe, newlabels='mZm' )
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("vMZm","vDZm"), free=TRUE, values=svVa, newlabels='vZm' )
fitEMVZ <- mxRun( modelEMVZ, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMO, fitEMVO, fitEMVZ) )
round(fitEMVZ$output$estimate,4)
fitGofs(fitEMVZ)

# Constrain expected Means and Variances to be equal across twin order and zygosity and SS/OS
modelEMVP <- mxModel( fitEMVZ, name="eqMVarsSos" )
modelEMVP <- omxSetParameters( modelEMVP, label=c("mZf","mDZo1"), free=TRUE, values=svMe, newlabels='mZf' )
modelEMVP <- omxSetParameters( modelEMVP, label=c("vZf","vDZo1"), free=TRUE, values=svVa, newlabels='vZf' )
modelEMVP <- omxSetParameters( modelEMVP, label=c("mZm","mDZo2"), free=TRUE, values=svMe, newlabels='mZm' )
modelEMVP <- omxSetParameters( modelEMVP, label=c("vZm","vDZo2"), free=TRUE, values=svVa, newlabels='vZm' )
fitEMVP <- mxRun( modelEMVP, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMO, fitEMVO, fitEMVZ, fitEMVP) )
round(fitEMVP$output$estimate,4)
fitGofs(fitEMVP)

# Constrain expected Means and Variances to be equal across twin order and zygosity and SS/OS and sex
modelEMVS <- mxModel( fitEMVP, name="eqMVarsSex" )
modelEMVS <- omxSetParameters( modelEMVS, label=c("mZf","mZm"), free=TRUE, values=svMe, newlabels='mZ' )
modelEMVS <- omxSetParameters( modelEMVS, label=c("vZf","vZm"), free=TRUE, values=svVa, newlabels='vZ' )
fitEMVS <- mxRun( modelEMVS, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMO, fitEMVO, fitEMVZ, fitEMVP, fitEMVS) )
round(fitEMVS$output$estimate,4)
fitGofs(fitEMVS)

mxEval(cov2cor(MZf.covMZf), fit, T)
mxEval(cov2cor(DZf.covDZf), fit, T)
mxEval(cov2cor(MZm.covMZm), fit, T)
mxEval(cov2cor(DZm.covDZm), fit, T)
mxEval(cov2cor(DZo.covDZo), fit, T)

I would like to add to the model the correlations for each group (MZM,MZF,DZM,DZF,DZO) I have this correlations in the submodels but not in the model. I mean when I run summary(fit) this values are only in the Confident invervals.
Therefore I can't use omxsetparemateres to equal the correlations and test the model (for example corMZF equal to corMZM)

I think I should create five matrix for the correlations but I don't know how to do this. Also I tried adding CICORMZF, CICORMZM... to the model and submodels but I only have the correlations in the submodels.

Log in or register to post comments
Replied on Wed, 02/07/2018 - 12:38
Picture of user. AdminRobK Joined: 01/24/2014

Aren't you working with two binary phenotypes? If so, the saturated model still needs to have thresholds, and constraints imposed to identify the scale of the latent liability.
Log in or register to post comments
Replied on Wed, 02/07/2018 - 13:02
No user picture. JuanJMV Joined: 07/20/2016

Yes sorry, I am working with two binary variables and one continous. This script is for the continous variable. I am testing if there are sex differences (in the continous variable). It is a 5 group saturated univariate model but I don't know how to add the correlations for the 5 groups to the model (actually I have them in the submodels)

Log in or register to post comments
Replied on Thu, 02/08/2018 - 11:45
Picture of user. AdminRobK Joined: 01/24/2014

I could be sure that the soltuion is "good" I mean since retry limit have been reached I don't know if I should try other starting values or something else (although there is a solution)

What was the warning generated during the final fit? Did you get status Red or not? If you're really concerned, you could increase the number of extra tries. You could also try my Nelder-Mead optimizer, but I doubt it would work any better than CSOLNP.

I would like to add to the model the correlations for each group (MZM,MZF,DZM,DZF,DZO) I have this correlations in the submodels but not in the model. I mean when I run summary(fit) this values are only in the Confident invervals.
Therefore I can't use omxsetparemateres to equal the correlations and test the model (for example corMZF equal to corMZM)

I think I should create five matrix for the correlations but I don't know how to do this. Also I tried adding CICORMZF, CICORMZM... to the model and submodels but I only have the correlations in the submodels.

Ah, I see. The reason for why that doesn't work is unrelated to whether the correlations are in the submodel or in the container (so, you don't need to create new objects for the correlation matrices; the ones you have in the submodels already will do). The reason for why that doesn't work is that the correlations are not statistical parameters of the MxModel object. Parameters are elements of MxMatrices (or values of MxPaths), but the correlations are elements of MxAlgebras. That is, they are not parameters themselves, but are instead functions of parameters (informally called "algebras"). omxSetParameters() and the 'free parameters' table in summary() output only deal with parameters, not algebras. In contrast, confidence intervals can be calculated for algebras. MxConstraints can also be imposed on algebras.

If you want to fit models in which the correlations are constrained across groups in some way, you basically have two options. One option is to reparameterize your MxModels so that the correlations are in fact parameters rather than algebras. For advice on doing that, see this other thread. The other option, which would involve less work, is to use MxConstraints. For instance, if you want to constrain the MZF and MZM correlations to be equal to one another, make an object like this,

MZconstraint <- mxConstraint(MZm.corMZm[1,2] == MZf.corMZf[1,2])

and put it in the container MxModel. I recommend using SLSQP when MxConstraints are involved.

You could instead request confidence intervals on the differences in the correlations. For instance,

MZdiff <- mxAlgebra(MZm.corMZm[1,2] - MZf.corMZf[1,2], name="MZdiff")
ciMZdiff <- mxCI("MZdiff")

, and put them into the container MxModel.

Log in or register to post comments
Replied on Mon, 02/12/2018 - 13:03
No user picture. JuanJMV Joined: 07/20/2016

Hi Rob! thank you so much again!

1-I think the warnings are for the CI since I have NA in some of them. I don't have mx status red now.
The message is "Retry limit reached

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.086752181736907,-0.00393385954682701,0.00457737378618429,0.48126084828...."

I have tried changing the optimizer and 100 extra tries but I have the same resuls and always retry limit is reached.

In this model we have fixed the nonshared-environmental variances to 1 but why not for A and C?

2-Where can I find information about Common-Pathway and Independent-Pathway models? I need to use another script right?

3- With regards to the correations for the five group univariate model it works perfect now! thank you!!!

Log in or register to post comments
Replied on Tue, 02/13/2018 - 13:54
Picture of user. AdminRobK Joined: 01/24/2014

1-I think the warnings are for the CI since I have NA in some of them. I don't have mx status red now.
The message is "Retry limit reached

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.086752181736907,-0.00393385954682701,0.00457737378618429,0.48126084828...."

I have tried changing the optimizer and 100 extra tries but I have the same resuls and always retry limit is reached.

If you run summary() with argument verbose=TRUE, you can see diagnostics for the confidence intervals. Also, mxTryHardOrdinal() runs with default argument exhaustive=TRUE, so by default it will always use up all of its extra tries.

2-Where can I find information about Common-Pathway and Independent-Pathway models? I need to use another script right?

You could always look around in the materials from the 2016 Boulder Workshop, e.g. here. See also:
Posthuma, D. (2009). Multivariate genetic analysis. In Kim, Y-K. (Ed.), Handbook of Behavior Genetics (pp. 47-59). New York: Springer Science+Business Media, LLC.
Neale, M.C., & Cardon, L.R. (1992). Methodology for genetic studies of twins and families. Dordrecht, The Netherlands: Kluwer Academic Publishers

Log in or register to post comments
Replied on Thu, 02/22/2018 - 12:12
No user picture. JuanJMV Joined: 07/20/2016

Hi Rob,

Thanks for your response. The diagnosis is "alpha level not reached infeasible start". I have tried with the optimizer "SLSQP" and removing lower bounds on the ACE path coefficients (as I saw in another of your post). But It does not work. It took more than 20 hours and I had to stop the analysis.

What should I do?

Thanks!!

Log in or register to post comments
Replied on Mon, 03/12/2018 - 11:30
Picture of user. AdminRobK Joined: 01/24/2014

Hi, Juan. Sorry for the delayed reply. I was at the Boulder Workshop last week.

I am not really sure that CSOLNP found a minimum of the fitfunction. The information matrix isn't positive-definite, and the gradient is asymmetric around every free-parameter estimate. How many extra tries did you run it with? Nonetheless, with ordinal data, I guess it's still possible that's close to a solution. Could you post your current script, so I know exactly what you're running?

Log in or register to post comments
Replied on Mon, 03/12/2018 - 13:48
No user picture. JuanJMV Joined: 07/20/2016

Hi! Don't worry! I wanted to go but.... I couldn't :(

Here the script. I think I also tried with 101 extraTries and It was the same result.

# Load Libraries & Options
library(OpenMx)
library(psych)
source("miFunctions.R")
mxOption(NULL,"Default optimizer","CSOLNP")
# ------------------------------------------------------------------------------
# PREPARE DATA
#filename <- "tri"
#sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE)
# Load Data
twinData <- read.table ("data.dat", sep= "\t", dec= ",", header=T, na.strings= "-9")

dim(twinData)
describe(twinData, skew=F)

# Select Continuous Variables
varsc <- c('SL') # list of continuous variables names
nvc <- 1 # 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('NE',"LB") # list of ordinal variables names
nvo <- 2 # 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('SL','NE',"LB") # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,1,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')
nc <- 3 # number of covariates

# Select Data for Analysis
mzData <- subset(ordData, Zyg==1|Zyg==3, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2|Zyg==4|Zyg==5, c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$NE1 <- mxFactor(mzDataF$NE1, levels =c(1,2))
mzDataF$NE2 <- mxFactor(mzDataF$NE2, levels =c(1,2))
dzDataF$NE1 <- mxFactor(dzDataF$NE1, levels =c(1,2))
dzDataF$NE2 <- mxFactor(dzDataF$NE2, levels =c(1,2))
mzDataF$LB1 <- mxFactor(mzDataF$LB1, levels =c(1,2))
mzDataF$LB2 <- mxFactor(mzDataF$LB2, levels =c(1,2))
dzDataF$LB1 <- mxFactor(dzDataF$LB1, levels =c(1,2))
dzDataF$LB2 <- mxFactor(dzDataF$LB2, levels =c(1,2))

# 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,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(3,1.5,1.5) # start value for means
svPa <- .3 # start value for path coefficient
svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix
svPe <- .4 # 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 <- 0.5 # start value for first threshold
svITh <- 0.2 # start value for increments
svTh <- -1 # start value for thresholds
lbTh <- -1 # 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" )
pathB <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=labBe, name="b" )

# 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(b%*%Age),t(b%*%Age)), 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), name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), name="e" )

# Create Algebra for Variance Comptwonts
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxMatrix(
type="Symm",nrow=nv,
free=c(T,T,T,F,T,F),
values=c(0.8,0,0,1,0,1),
labels=c("e_1_1", "e_2_1", "e_3_1", "e_2_2", "e_3_2", "e_3_3"),
lbound=c(0.0001,rep(NA,5)),
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

# 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(pathB, meanG, thinG, inc, threG, matI, invSD,
pathA, pathC, covA, covC, covE, covP, corA, corC, corE)
defs <- list(defAge)
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","MZ.rA", "MZ.rC", "MZ.rE","MZ.A","MZ.C","MZ.E" ))# "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, modelMZ, modelDZ, multi, estVC, ciACE )

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

# Run ACE Model
fitACE <- mxTryHardOrdinal( 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 <- mxTryHardOrdinal( modelAE, intervals=F )
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=F )
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=F )
mxCompare( fitAE, fitE )
fitGofs(fitE)
fitEsts(fitE)

# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )
#round(rbind(fitACE$VC$result,fitAE$VC$result,fitCE$VC$result,fitE$VC$result),4)
fitACE$algebras
mxEval(cov2cor(V), fitACE, T)
fitAE$algebras
mxEval(cov2cor(V), fitAE, T)

sink()
save.image(paste(filename,".Ri",sep=""))

Thank you so much

Log in or register to post comments
Replied on Fri, 03/16/2018 - 15:00
Picture of user. AdminRobK Joined: 01/24/2014

Do you already know for sure that you're going to report CIs from the ACE model? If not, you should first try to be sure you're getting to a decent maximum-likelihood solution for all four models, select the one with smallest AIC, and only then worry about confidence intervals. Try mxTryHardOrdinal() without intervals=T; stick with CSOLNP for now. You can subsequently use omxRunCI() on a fitted MxModel object to calculate CIs without redoing the primary optimization to find the MLE.

I don't see that the primary optimization should be too difficult. Granted, you're working with ordinal data, but your ordinal variables are binary with neither outcome being rare, and you don't have a whole lot of missing data.

The CI details table is so big that it can be hard to read. You could simplify it if you request CIs for a vector of MxAlgebra elements, e.g.,

mxCI( c("VC[1,1]", "VC[1,2]" #<--Continue in this fashion
))

and just skip the diagonal elements of the correlation matrices.

It would be a good idea to make your mxTryHardOrdinal() results reproducible. The best way to do so is probably to add a line like

set.seed(180316)

before each call to mxTryHardOrdinal(). Also, if you pass argument bestInitsOutput=TRUE to mxTryHardOrdinal(), you can see the start values that the program used to find its best solution, and then modify your script to actually use those start values.

Log in or register to post comments
Replied on Sun, 03/25/2018 - 13:13
No user picture. JuanJMV Joined: 07/20/2016

Hi Rob,

Thank you so much, now the model estimates all the CI for the best fitted model. I have adjust the starting values and I have modified the line mxCI( c("VC[1,1]…. As you told me.

I had to use the optimizer “SLSQP” is it ok?

Just one more question I am not sure why we have to fix to 1 only nonshared-environmental variance and not for the other components?

Thank you so much again!!!

Log in or register to post comments
Replied on Mon, 03/26/2018 - 09:55
Picture of user. AdminRobK Joined: 01/24/2014

I had to use the optimizer “SLSQP” is it ok?

If it worked, it seems OK to me!

Just one more question I am not sure why we have to fix to 1 only nonshared-environmental variance and not for the other components?

The way OpenMx models ordinal variables is by assuming that the ordinal variable represents a latent, normally-distributed continuous variable that has been discretized into ordered categories by "cutting" at the thresholds (a binary variable will only have 2 categories and 1 threshold). Since the underlying continuous variable isn't observed, the model isn't identified without some additional constraint. With a binary variable, either the threshold or the underlying mean has to be fixed to some value; typically the variance is also fixed, to 1.0. However, in an ACE twin model, fixing the variance (i.e., the total phenotypic variance) to 1.0 requires an MxConstraint. I was concerned the MxConstraint was complicating the optimization of your model, so I suggested an alternative constraint, which is to fix the nonshared-environmental variance to 1.0. When the model is identified in that fashion, the values of the additive-genetic and shared-environmental variance components are interpretable as the factors by which those variance components are proportionately larger than (if >1.0) or smaller than (if <1.0) the nonshared-environmental variance.

Log in or register to post comments
Replied on Mon, 05/07/2018 - 08:02
Picture of user. Julia Joined: 03/29/2012

Hi Rob,

Just reading this great thread which has helped me a lot in setting up my model and I am wondering about omxRunCI. This sounds like a great function to calculate the CI's without refitting the whole model. However, I cannot find it. I have OpenMx version 2.6.9 and when I type ?omxRunCI, it says that no documentation is found. Could you please advise?
Thank you,
Julia

Log in or register to post comments
Replied on Mon, 05/07/2018 - 08:25
Picture of user. Julia Joined: 03/29/2012

Sorry, my bad! Found out that this function is only available with version 2.9.6. Just updated OpenMx, but can't figure out how to use omxRunCI(). The CI's are specified as follows.

ciPars <- mxCI( c('MZ.expThre', 'MZ.bAge', 'MZ.bSex')

And when I pass AceFit into omRunCI, I get the following:

> AceFitCI = omxRunCI(AceFit)
Error in imxConvertIdentifier(constraint@jac[1], modelname, namespace) :
no slot of name "jac" for this object of class "MxConstraint"

Could you please help me with this one?

Log in or register to post comments
Replied on Mon, 05/07/2018 - 09:53
Picture of user. AdminRobK Joined: 01/24/2014

Found out that this function is only available with version 2.9.6.

Actually, it's available in v2.7.18 or newer.


Error in imxConvertIdentifier(constraint@jac[1], modelname, namespace) :
no slot of name "jac" for this object of class "MxConstraint"

I am willing to bet a lot of money that the MxModel and other OpenMx objects were created with an older version of the package--something older than v2.7.x--and were saved to disk in an .RData file (possibly as part of R autosaving your workspace at the end of a session) and subsequently reloaded after you updated to v2.9.6. The MxConstraint object class was changed between version 2.6.x and version 2.7.4 by adding new slots to it, in order to allow users to provide analytic constraint Jacobians. You'll need to reconstruct the MxConstraint objects that are in your MxModel, and put them into the MxModel in some way. Probably the best thing to do is just re-run your whole script, since it's very possible that other object classes have had their signatures altered since v2.6.x. Alternately, you could reconstruct the MxConstraints only, and directly overwrite the 'constraints' slot of the model (or of the submodels, if that's where the MxConstraints go).

If you could run traceback() right after seeing the error message and post the output, that would be helpful, since it's possible some code changes could be made so that future OpenMx releases don't suffer from this lack of backward compatibility.

Log in or register to post comments
Replied on Mon, 05/07/2018 - 14:18
Picture of user. Julia Joined: 03/29/2012

Yes, the models were indeed built under the older 2.6.9 version of OpenMx. Actually now when I updated to 2.9.6 and tried to load workspace with my fitted models, I can't even run the summary of the models. Getting the following:

> summary(AceFit)
Error in imxConvertIdentifier(constraint@jac[1], modelname, namespace) :
no slot of name "jac" for this object of class "MxConstraint"
In addition: Warning messages:
1: package ‘OpenMx’ was built under R version 3.4.3
2: In fun(libname, pkgname) : bytecode version mismatch; using eval
> traceback()
17: imxConvertIdentifier(constraint@jac[1], modelname, namespace)
16: FUN(X[[i]], ...)
15: lapply(slot(model, slotName), convertFunction, model@name, namespace)
14: collectComponentsHelper(model, namespace, slotName, convertFunction)
13: collectComponents(model, namespace, "constraints", qualifyNamesConstraint)
12: imxFlattenModel(model, namespace, FALSE)
11: EvalInternal(expression, model, modelvariable, compute, show,
defvar.row, cache, cacheBack, 2L)
10: mxEval(diag2vec(ACE.V), model, compute = TRUE)
9: eval(expr, envir, enclos)
8: eval(substitute(mxEval(x, model, compute = TRUE), list(x = leftHandSide)))
7: FUN(X[[i]], ...)
6: lapply(X = X, FUN = FUN, ...)
5: sapply(constraints, calculateConstraintsHelper, model)
4: calculateConstraints(model, useSubmodels)
3: computeOptimizationStatistics(model, numStats, useSubmodels,
saturatedDoF, independenceDoF, retval)
2: summary.MxModel(AceFit)
1: summary(AceFit)

I can extract estimated parameters though. Is there a way to get the summary of the model fit or will I have to rerun all the models (it takes a couple of hours for each)...

Log in or register to post comments
Replied on Mon, 05/07/2018 - 15:53
Picture of user. AdminRobK Joined: 01/24/2014

Good idea, Joshua. I have a feeling, though, that the confidence-interval optimization will take the lion's share of the computational effort.
Log in or register to post comments
Replied on Tue, 11/12/2019 - 07:06
No user picture. yaning Joined: 04/12/2019

Hello, I want to run a G*E model (2 binary variables ) using the script from “JuanJMV”. Unfortunately, there seems to be some errors in my model

I just met the following Warning message when I add one of the Constraint lines
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)

# Constraint on variance of Binary variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
var_constraint <- mxConstraint( expression=diag2vec(V0)==Unv1, name="Var1" )
or
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
var_constraint <- mxConstraint(expression=rbind(V0[2,2],V0[1,1]) == Unv1, name="Var1")

Another errors occurred when I remove the Constraint lines and define covE as:
myVarE <- mxMatrix( type="Symm",nrow=nv, ncol=nv, free=c(F,T,F), values=c(1, 0.1, 1),
labels=c("e_1_1", "e_2_1", "e_2_2"), lbound=c(0.0001,rep(NA,2)),name="E0")

** Information matrix is not positive definite (not at a candidate optimum).
Be suspicious of these results. At minimum, do not trust the standard errors.

Any suggestions from you?
Thanks!

File attachments
code.txt
Log in or register to post comments
  Login

News