You are here

doubt with binary variables

17 posts / 0 new
Last post
JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
doubt with binary variables

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
full script?

Without the context of the full script, I can't really answer your question.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thanks Rob for your prompt

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 )
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
modification

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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much Rob,

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)

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
answers

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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thanks Rob,

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
reparameterization
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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Sorry! I am using OpenMx

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!!!

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Do you think that these

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!!!

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
more answers
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?

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Sorry! When I said if these

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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
thresholds?

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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Yes sorry, I am working with

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)

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
correlations
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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Hi Rob! thank you so much

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!!!

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
CIs
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