Multivariate cholesky model

Posted on Fri, 09/14/2018 - 03:56

sealx017 Joined: 08/17/2018

Forums

Hi,

I have a dataset of twins on 4 time-points. I can fit the Multivariate cholesky decompositon model with ACE components. It considers unstrcutured covariance matrices for A, C and E component. I was wondering if its possible to restrict it to a special structure? Like, say compound symmetric or AR1.

Thanks

Replied on Fri, 09/14/2018 - 12:05

AdminRobK Joined: 01/24/2014

yes; reparameterize

I was wondering if its possible to restrict it to a special structure? Like, say compound symmetric or AR1

Yes, that's possible, but it would be difficult to do under the Cholesky parameterization. It could be done with Cholesky, but you'd need to use MxConstraints on the off-diagonal elements of the _A_, _C_, and _E_ correlation matrices.

I recommend instead using a correlated-factors paramterization. Under that parameterization, the _A_, _C_, and _E_ covariance matrices would each be equal to a matrix of loadings times a factor correlation matrix times the transpose of the matrix of loadings. There would be a loadings matrix and a factor correlation matrix for each of _A_, _C_, and _E_. Each phenotype would load only onto its own _A_ (or whatever) factor; the loadings matrix could just be a diagonal matrix. The factor correlation matrix would have the structure you want. For compound symmetry, the factor correlation matrix could be an MxMatrix of `type="Stand"` and the same parameter label for all the off-diagonals. For AR1, I would suggest making a 1x1 matrix for the first-order correlation, and then define the factor correlation matrix as an MxAlgebra assembled, via `cbind()` and `rbind()`, from the first-order correlation raised to the appropriate power per cell (with ones on the diagonal of course). Loehlin (1996) and Neale, Røysamb, & Jacobson (2006) may be useful references.

Replied on Fri, 09/14/2018 - 15:48

sealx017 Joined: 08/17/2018

Can you guide me how to addd

Can you guide me how to addd constraints in Cholesky Model? This is how I am currently specifying.



pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                       values=svPa, labels=labLower("a",nv), lbound=lbPa, name="a" )

pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                       values=svPa, labels=labLower("c",nv), lbound=lbPa, name="c" )

pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                       values=svPa, labels=labLower("e",nv), lbound=lbPa, name="e" )

# Matrices generated to hold A, C, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" )

Replied on Fri, 09/14/2018 - 16:53

AdminRobK Joined: 01/24/2014

MxConstraints

Can you guide me how to addd constraints in Cholesky Model? This is how I am currently specifying.

OK, first make _A_, _C_, and _E_ correlation matrices:
corrA <- mxAlgebra(cov2cor(A), name="corA") corrC <- mxAlgebra(cov2cor(C), name="corC") corrE <- mxAlgebra(cov2cor(E), name="corE")
Now, you need to impose the desired structure via MxConstraints. To impose a compound-symmetric correlation structure on _A_:
Acon1 <- mxConstraint(corA[2,1] == corA[3,1]) Acon2 <- mxConstraint(corA[3,1] == corA[4,1]) Acon3 <- mxConstraint(corA[4,1] == corA[3,2]) Acon4 <- mxConstraint(corA[3,2] == corA[4,2]) Acon5 <- mxConstraint(corA[4,2] == corA[4,3])
Do likewise for _C_ and _E_. You can give the MxConstraints names if you want, via argument `name`. You might be able to come up with a way to do all this using fewer `mxConstraint()` statements, perhaps with the `vechs()` function. At any rate, you need to put these new objects into either the MZ or DZ group model.

For AR1, try:
#First, establish what the first-order correlation will be: Acon1 <- mxConstraint(corA[1,2] == corA[2,3]) Acon2 <- mxConstraint(corA[2,3] == corA[3,4]) #Now, impose the decay with increasing lag: Acon3 <- mxConstraint(corA[1,3] == corA[1,2]^2) Acon4 <- mxConstraint(corA[1,4] == corA[1,2]^3) Acon5 <- mxConstraint(corA[2,4] == corA[1,2]^2)
And likewise for _C_ and _E_.

If you take this route, I don't recommend the on-load default optmizer, CSOLNP. Switch to SLSQP (or perhaps NPSOL if your build of OpenMx supports it), e.g. mxOption(NULL,"Default optimizer","SLSQP").

Replied on Fri, 09/14/2018 - 16:58

sealx017 Joined: 08/17/2018

That is great!!

Can we also add multiple equalities like,

Acon1 <- mxConstraint(corA[1,2] == corA[2,3]==corC[1,2])

Replied on Fri, 09/14/2018 - 17:11

AdminRobK Joined: 01/24/2014

That would be nice

That is great!!

Can we also add multiple equalities like,

Acon1 <- mxConstraint(corA[1,2] == corA[2,3]==corC[1,2])

That would be nice, but I'm pretty sure the answer is "no" (though I admit that I've never tried). What happens at runtime is that OpenMx converts the `expression` you provided to `mxConstraint()` into a standard form, which has a matrix-valued function on the LHS of the comparator, and on the RHS, a matrix of the same dimensions containing all zeroes. So, you could use fewer `mxConstraint()` statements (and thereby make fewer MxConstraint objects) if you did something like this:
uno <- mxMatrix(type="Unit",nrow=6,ncol=1,name="Uno") vechsA <- mxAlgebra(vechs(A), name="vsA") Acon <- mxConstraint(vsA[1,1]%x%Uno == vsA)

Replied on Fri, 09/14/2018 - 17:08

sealx017 Joined: 08/17/2018

Matrices declared to store a,

  
  # Matrices declared to store a, c, and e Path Coefficients

  pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                         values=svPa, labels=labLower("a",nv), lbound=lbPa, name="a" )

  pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                         values=svPa, labels=labLower("c",nv), lbound=lbPa, name="c" )

  pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                         values=svPa, labels=labLower("e",nv), lbound=lbPa, name="e" )

# Matrices generated to hold A, C, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) corrA <- mxAlgebra(cov2cor(A), name="corA") corrC <- mxAlgebra(cov2cor(C), name="corC") corrE <- mxAlgebra(cov2cor(E), name="corE") Acon1 <- mxConstraint(corA[2,1] == corA[3,1]) Acon2 <- mxConstraint(corA[3,1] == corA[4,1]) Acon3 <- mxConstraint(corA[4,1] == corA[3,2]) Acon4 <- mxConstraint(corA[3,2] == corA[4,2]) Acon5 <- mxConstraint(corA[4,2] == corA[4,3])
I used the above code and was looking at results, but I still get an unstructured correlation matrix.

> A<-CholAceFit$A$result > C<-CholAceFit$C$result > E<-CholAceFit$E$result > A [,1] [,2] [,3] [,4] [1,] 0.3078549 0.2373371 0.1750760 0.1582570 [2,] 0.2373371 0.3636521 0.2733723 0.2153947 [3,] 0.1750760 0.2733723 0.2760969 0.2189010 [4,] 0.1582570 0.2153947 0.2189010 0.2465859 > cov2cor(A) [,1] [,2] [,3] [,4] [1,] 1.0000000 0.7093319 0.6005134 0.5743887 [2,] 0.7093319 1.0000000 0.8627416 0.7192964 [3,] 0.6005134 0.8627416 1.0000000 0.8389436 [4,] 0.5743887 0.7192964 0.8389436 1.0000000

Replied on Fri, 09/14/2018 - 17:16

AdminRobK Joined: 01/24/2014

Well, something's clearly not

Well, something's clearly not right. Please post your full script.

Edit: Also, what's your full output from `mxVersion()`?

Replied on Fri, 09/14/2018 - 18:24

sealx017 Joined: 08/17/2018

require(OpenMx)



require(OpenMx)

  require(psych)

  setwd("C:\\Users\\Slsouvik\\Desktop\\openmx")

  source("myFunctions.R")

  source("GenEpiHelperFunctions.R")
  # --------------------------------------------------------------------

  # PREPARE DATA
  nmz<-mz_proportion*nfam

  ndz<-nfam-nmz
  nv        <- ntrait     # number of variables

  Vars      <- paste("X",c(1:nv),sep="")

  ntv       <- nv*2    # number of total variables

  selVars   <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")

  # Select Data for Analysis
  #Y1<-as.data.frame(scale(Y_mz,center=FALSE,scale = TRUE))

  #Y2<-as.data.frame(scale(Y_dz,center=FALSE,scale = TRUE))

  Y1<-as.data.frame((Y_mz))

  Y2<-as.data.frame(Y_dz)

  Y1<-as.data.frame(cbind(Y1[seq(1,(2*nmz),2),],Y1[seq(2,(2*nmz),2),]))

  Y2<-as.data.frame(cbind(Y2[seq(1,(2*ndz),2),],Y2[seq(2,(2*ndz),2),]))

  Y<-rbind(Y1,Y2)

  Y<-cbind(Y,c(rep(1,nmz),rep(2,ndz)))
  colnames(Y)[1:(2*ntrait)]<-selVars

  colnames(Y)[(2*ntrait+1)]<-"zyg2"
  # Select random subset to reduce time to fit examples

  testData <- Y
  mzData    <- subset(testData, zyg2==1, c(selVars))

  dzData    <- subset(testData, zyg2==2, c(selVars))

  describe(mzData, skew=F)

  describe(dzData, skew=F)

  dim(mzData)

  dim(dzData)
  #mzData    <- subset(testData, zyg2==1, c(selVars, varscov))

  #dzData    <- subset(testData, zyg2==2, c(selVars, varscov))
  # Set Starting Values

  svTh      <- .8                        # start value for thresholds

  svPa      <- .4                        # start value for path coefficient

  svPe      <- .8                        # start value for path coefficient for e

  lbPa      <- .0001                     # start value for lower bounds
  # Saturated Model

  # Set Starting Values

  svMe      <- c(0,0,0,0)            # start value for means

  svVa      <- valDiag(ntv,1.0)          # start values for variances

  lbVa      <- valODiag(ntv,.0001,-10)
  # -------|---------|---------|---------|---------|---------|---------|

  # Algebra for expected Mean Matrices in MZ & DZ twins

  meanMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,

                         values=svMe, labels=labFull("meMZ",1,ntv), name="expMeanMZ" )

  meanDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,

                         values=svMe, labels=labFull("meDZ",1,ntv), name="expMeanDZ" )
  # Matrices declared to store a, c, and e Path Coefficients

  pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                         values=svPa, labels=labLower("a",nv), lbound=lbPa, name="a" )

  pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                         values=svPa, labels=labLower("c",nv), lbound=lbPa, name="c" )

  pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

                         values=svPa, labels=labLower("e",nv), lbound=lbPa, name="e" )
  # Matrices generated to hold A, C, and E computed Variance Components

  covA      <- mxAlgebra( expression=a %*% t(a), name="A" )

  covC      <- mxAlgebra( expression=c %*% t(c), name="C" )

  covE      <- mxAlgebra( expression=e %*% t(e), name="E" )

  corrA <- mxAlgebra(cov2cor(A), name="corA")

  corrC <- mxAlgebra(cov2cor(C), name="corC")

  corrE <- mxAlgebra(cov2cor(E), name="corE")

  Acon1 <- mxConstraint(corA[2,1] == corA[3,1])

  Acon2 <- mxConstraint(corA[3,1] == corA[4,1])

  Acon3 <- mxConstraint(corA[4,1] == corA[3,2])

  Acon4 <- mxConstraint(corA[3,2] == corA[4,2])

  Acon5 <- mxConstraint(corA[4,2] == corA[4,3])
  # Algebra to compute total variances and standard deviations (diagonal only)

  covP      <- mxAlgebra( expression=A+C+E, name="V" )

  matI      <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")

  invSD     <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")

  var1     <- mxConstraint( expression=diag2vec(V)==1, name="Var1" )
  # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins

  meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,

                         values=svMe, labels=labFull("me",1,nv), name="expMean" )
  covMZ     <- mxAlgebra( expression= rbind( cbind(V, A+C), cbind(A+C, V)), name="expCovMZ" )

  covDZ     <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovDZ" )
  # Data objects for Multiple Groups

  dataMZ    <- mxData( observed=mzData, type="raw" )

  dataDZ    <- mxData( observed=dzData, type="raw" )
  # Objective objects for Multiple Groups

  objMZ     <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars )

  objDZ     <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars )
  # Combine Groups
  pars      <- list( pathA,pathC, pathE, covA, covC, covE, covP, matI, invSD, meanG )
  modelMZ   <- mxModel( pars, covMZ, dataMZ, objMZ, name="MZ" )

  modelDZ   <- mxModel( pars, covDZ, dataDZ, objDZ, name="DZ" )

  #obj       <- mxFitFunctionMultigroup(c("MZ","DZ"))

  #CholAceModel  <- mxModel( "CholACE", modelMZ, modelDZ, obj )
  minus2ll  <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )

  obj       <- mxAlgebraObjective( "m2LL" )

  CholAceModel  <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj )

  # ------------------------------------------------------------------------------

  # RUN GENETIC MODEL
  # Run Cholesky Decomposition ACE model

  CholAceFit    <- mxRun(CholAceModel)
  # ACE Covariance Matrices & Proportions of Variance Matrices

  ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V")

  ACEcovLabels <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE")

A<-CholAceFit$A$result C<-CholAceFit$C$result E<-CholAceFit$E$result

This is the code.

> mxVersion()
OpenMx version: 2.9.9 [GIT v2.9.9]
R version: R version 3.3.1 (2016-06-21)
Platform: x86_64-w64-mingw32
Default optimiser: CSOLNP
NPSOL-enabled?: Yes
OpenMP-enabled?: No

Replied on Sat, 09/15/2018 - 15:11

AdminRobK Joined: 01/24/2014

objects need to go in models

Like I said, the new objects have to actually go into an MxModel object. For instance,



modelMZ   <- mxModel( pars, covMZ, dataMZ, objMZ, corrA, corrC, corrE, Acon1, Acon2, Acon3, Acon4, Acon5, name="MZ" )

. I have a few additional remarks about the script as well. First, is there a reason you're not correcting for the main effects of age and sex? Second, I don't see any syntax in your script where the dataset is loaded. Are you just restoring an R workspace at the start of the session? If so, I strongly discourage doing that.

Replied on Sat, 09/15/2018 - 17:25

sealx017 Joined: 08/17/2018

I am getting this error,

I am getting this error,
do you think it is due to the data? I mean dataset does not suit the used model and hence estimates are getting biased and wrong.

Running CholACE with 34 parameters

Error: The job for model 'CholACE' exited abnormally with the error message: fit is not finite (The continuous part of the model implied covariance (loc2) is not positive definite in data 'MZ.data' row 28. Detail:
covariance = matrix(c( # 8x8
0.480003, 0.320338, 0.324131, 1.9695e+019, 0.320001, 0.160337, 0.164131, 1.9695e+019
, 0.320338, 3.22184e+026, 3.25516e+031, 4.14316e+016, 0.160337, 3.22184e+026, 3.25516e+031, 4.14316e+016
, 0.324131, 3.25516e+031, 1.84611e+038, 1.41805e+018, 0.164131, 3.25516e+031, 1.84611e+038, 1.41805e+018
, 1.9695e+019, 4.14316e+016, 1.41805e+018, 2.42432e+039, 1.9695e+019, 4.14316e+016, 1.41805e+018, 2.42432e+039
, 0.320001, 0.160337, 0.164131, 1.9695e+019, 0.480003, 0.320338, 0.324131, 1.9695e+019
, 0.160337, 3.22184e+026, 3.25516e+031, 4.14316e+016, 0.320338, 3.22184e+026, 3.25516e+031, 4.14316e+016
, 0.164131, 3.25516e+031, 1.84611e+038, 1.41805e+018, 0.324131, 3.25516e+031, 1.84611e+038, 1.41805e+018
, 1.9695e+019, 4.14316e+016, 1.41805e+018, 2.424
In addition: Warning message:
In model 'CholACE' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard()

Replied on Sat, 09/15/2018 - 18:49

AdminRobK Joined: 01/24/2014

some tips

do you think it is due to the data? I mean dataset does not suit the used model and hence estimates are getting biased and wrong.

Possibly, but some of those elements of the model-expected covariance matrix just look ridiculously wrong. I think it's more likely that the system of nonlinear constraints your model imposes is making optimization difficult, and/or the start values in your script are poor. A few tips:

If you use this website's content-search feature for "start values", you can find plenty of discussion and advice about how to choose good start values. One of the diagnostic messages you're seeing is advising you that the start values aren't feasible. So, perhaps they can be improved.
The message also suggests `mxTryHard()`, which is a function that makes multiple attempts to fit a model to data, and randomly perturbs the start values between attempts. Try replacing `mxRun()` with `mxTryHard()`in your script. If you do so, it would be a good idea to at least skim the help page for `mxTryHard()`.
Did you switch optimizers to SLSQP or NPSOL as I suggested here?
And on a related note, did you put the MxConstraints into more than one `mxModel()` statement? If so, your best bet is to only put them into one. There are several Known Issues with redundant MxConstraints, and only using NPSOL as the optimizer substantially mitigates them.

I'll remind you that you're using the "Cholesky parameterization plus MxConstraints" approach to fit these structured-correlation-matrix models, which goes against my recommendation (to instead use a correlated-factor parameterization) in the first place. You are not the first user to encounter difficulties using such an approach.

Replied on Sat, 09/15/2018 - 21:08

sealx017 Joined: 08/17/2018

I am getting fine result on a

I am getting fine result on a simpler data. I guess for the other one, the starting values need to be more accurate.

Regarding this,

Acon1 <- mxConstraint(corA[1,2] == corA[2,3]==corC[1,2])

Do you think it is possible in OpenMx to allow such restriction in any way, i.e, correlation of A component is the same as that of C? like if I do

Acon1 <- mxConstraint(corA[1,2] == corA[2,3]
Ccon1 <- mxConstraint(corC[1,2] == corA[2,3]

Would it work?
I tried once, it didnt work.

Replied on Sun, 09/16/2018 - 15:00

AdminRobK Joined: 01/24/2014

That is indeed possible to do

That is indeed possible to do. Assuming you put the new MxConstraints into an MxModel, my guess for why "it didn't work" is that the optimizer is failing to find a solution that satisfies the constraints.

Replied on Sun, 09/16/2018 - 15:50

sealx017 Joined: 08/17/2018

I used this code,



Acon1 <- mxConstraint(corA[2,1] == corA[3,1])

  Acon2 <- mxConstraint(corA[3,1] == corA[4,1])

  Acon3 <- mxConstraint(corA[4,1] == corA[3,2])

  Acon4 <- mxConstraint(corA[3,2] == corA[4,2])

  Acon5 <- mxConstraint(corA[4,2] == corA[4,3])
  Ccon1 <- mxConstraint(corA[2,1] == corC[3,1])

  Ccon2 <- mxConstraint(corA[3,1] == corC[4,1])

  Ccon3 <- mxConstraint(corA[4,1] == corC[3,2])

  Ccon4 <- mxConstraint(corA[3,2] == corC[4,2])

  Ccon5 <- mxConstraint(corA[4,2] == corC[4,3])
  Econ1 <- mxConstraint(corA[2,1] == corE[3,1])

  Econ2 <- mxConstraint(corA[3,1] == corE[4,1])

  Econ3 <- mxConstraint(corA[4,1] == corE[3,2])

  Econ4 <- mxConstraint(corA[3,2] == corE[4,2])

  Econ5 <- mxConstraint(corA[4,2] == corE[4,3])

  # Algebra to compute total variances and standard deviations (diagonal only)

  covP      <- mxAlgebra( expression=A+C+E, name="V" )

  matI      <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")

  invSD     <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")

  var1     <- mxConstraint( expression=diag2vec(V)==1, name="Var1" )
  # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins

  meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,

                         values=svMe, labels=labFull("me",1,nv), name="expMean" )
  covMZ     <- mxAlgebra( expression= rbind( cbind(V, A+C), cbind(A+C, V)), name="expCovMZ" )

  covDZ     <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovDZ" )
  # Data objects for Multiple Groups

  dataMZ    <- mxData( observed=mzData, type="raw" )

  dataDZ    <- mxData( observed=dzData, type="raw" )
  # Objective objects for Multiple Groups

  objMZ     <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars )

  objDZ     <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars )
  # Combine Groups
  pars      <- list( pathA,pathC, pathE, covA, covC, covE, covP,corrA, corrC, corrE, Acon1, Acon2, Acon3, Acon4, Acon5,Ccon1, Ccon2, Ccon3, Ccon4, Ccon5,Econ1, Econ2, Econ3, Econ4, Econ5, matI, invSD, meanG )

  #pars      <- list( pathA,pathC, pathE, covA, covC, covE, covP, matI, invSD, meanG )

  #pars      <- list( pathA,pathC, pathE, covA, covC, covE, covP,corrA, corrC, corrE, Acon1, Acon2, Acon3, Acon4, Acon5, matI, invSD, meanG )
 # modelDZ   <- mxModel( pars,corrA, corrC, corrE, Acon1, Acon2, Acon3, Acon4, Acon5,Ccon1, Ccon2, Ccon3, Ccon4, Ccon5,Econ1, Econ2, Econ3, Econ4, Econ5, covDZ, dataDZ, objDZ, name="DZ" )

 # modelMZ   <- mxModel( pars, corrA, corrC, corrE, Acon1, Acon2, Acon3, Acon4, Acon5,Ccon1, Ccon2, Ccon3, Ccon4, Ccon5,Econ1, Econ2, Econ3, Econ4, Econ5,covMZ, dataMZ, objMZ, name="MZ" )
 # modelDZ   <- mxModel( pars, covDZ, dataDZ, objDZ,corrA, corrC, corrE, Acon1, Acon2, Acon3, Acon4, Acon5, name="DZ" )

 # modelMZ   <- mxModel( pars,covMZ, dataMZ, objMZ, corrA, corrC, corrE, Acon1, Acon2, Acon3, Acon4, Acon5,name="MZ" )
  modelDZ   <- mxModel( pars, covDZ, dataDZ, objDZ, name="DZ" )

  modelMZ   <- mxModel( pars, covMZ, dataMZ, objMZ, name="MZ" )
  #obj       <- mxFitFunctionMultigroup(c("MZ","DZ"))

  #CholAceModel  <- mxModel( "CholACE", modelMZ, modelDZ, obj )
  minus2ll  <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )

  obj       <- mxAlgebraObjective( "m2LL" )

  CholAceModel  <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj )

  mxOption(NULL,"Default optimizer","SLSQP")

  # ------------------------------------------------------------------------------

  # RUN GENETIC MODEL

# Run Cholesky Decomposition ACE model CholAceFit <- mxRun(CholAceModel)
It does gives me A, C and E unstructured. But if I keep only Acon1,..,Acon5 in pars, it gives me compound symmetric correlation structure of A.

Replied on Mon, 09/17/2018 - 10:34

AdminRobK Joined: 01/24/2014

MxConstraints on whole matrices

Hmm. That's odd. Try again, retaining `Acon1` thru `Acon5`, and including two new MxConstraints:



Ccon <- mxConstraint(corC == corA)

Econ <- mxConstraint(corE ==corA)

Replied on Mon, 09/17/2018 - 15:48

sealx017 Joined: 08/17/2018

I tried what you said. Only

I tried what you said. Only retained con1 through Acon5 and added last two constraints. It is weird it did not work.

Replied on Tue, 09/18/2018 - 11:37

AdminRobK Joined: 01/24/2014

output?

Could you post some post-mxRun() output, such as the estimated values of the 3 correlation matrices, and the summary() output? I might be able to diagnose what's going wrong.

My advice to adopt the correlated-factors parameterization still stands.

Replied on Tue, 09/18/2018 - 17:46

sealx017 Joined: 08/17/2018

This is when I specify

This is when I specify compound symmetry for all 3,

A
[,1] [,2] [,3] [,4]
[1,] 44.04476 18.84978 24.12863 20.55763
[2,] 18.84978 44.75724 25.76197 20.95504
[3,] 24.12863 25.76197 49.13861 23.88618
[4,] 20.55763 20.95504 23.88618 42.12631
> C
[,1] [,2] [,3] [,4]
[1,] 5.9019408291 4.1313542437 0.0002429391 6.765129377
[2,] 4.1313542437 3.0425636732 0.0002088668 4.935798747
[3,] 0.0002429391 0.0002088668 1.9565944167 0.000469937
[4,] 6.7651293767 4.9357987469 0.0004699370 8.020701197
> E
[,1] [,2] [,3] [,4]
[1,] 13.707064 7.331347 6.261522 5.814909
[2,] 7.331347 12.460888 6.820699 6.240138
[3,] 6.261522 6.820699 10.797324 5.555669
[4,] 5.814909 6.240138 5.555669 11.915593

This is when I specify CS for only A,

[,1] [,2] [,3] [,4]
[1,] 42.59026 22.69322 23.48385 21.82458
[2,] 22.69322 45.94477 24.39115 22.66777
[3,] 23.48385 24.39115 49.20198 23.45751
[4,] 21.82458 22.66777 23.45751 42.49477
>
> C
[,1] [,2] [,3] [,4]
[1,] 7.8527242 1.2478234 0.6538513 5.8108463
[2,] 1.2478234 1.7890718 0.1040252 3.2639430
[3,] 0.6538513 0.1040252 1.3143163 0.4841337
[4,] 5.8108463 3.2639430 0.4841337 7.7436756
> E
[,1] [,2] [,3] [,4]
[1,] 13.683333 7.053824 6.238655 5.705981
[2,] 7.053824 12.407678 6.930252 6.152708
[3,] 6.238655 6.930252 10.851366 5.556251
[4,] 5.705981 6.152708 5.556251 11.896003

Replied on Wed, 09/19/2018 - 10:46

AdminRobK Joined: 01/24/2014

mxEval()

OK, but what do you get from mxEval(MZ.corA, CholAceFit, TRUE), and likewise for _C_ and _E_? Remember, the constraints are imposed on the correlation matrices, not the covariance matrices. And what about the summary() output?

News

Recent Posts

Multivariate cholesky model

yes; reparameterize

Can you guide me how to addd

MxConstraints

That is great!!

That would be nice

Matrices declared to store a,

Well, something's clearly not

require(OpenMx)

objects need to go in models

I am getting this error,

some tips

I am getting fine result on a

That is indeed possible to do

I used this code,

MxConstraints on whole matrices

I tried what you said. Only

output?

This is when I specify

mxEval()

News

Recent Posts