Multivariate cholesky model

20 posts / 0 new

Last post

Fri, 09/14/2018 - 04:56

sealx017

Offline

Joined: 08/17/2018 - 21:55

Multivariate cholesky model

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

Fri, 09/14/2018 - 13:05

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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.

Fri, 09/14/2018 - 16:48

sealx017

Offline

Joined: 08/17/2018 - 21:55

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

Fri, 09/14/2018 - 17:53

(Reply to #3) #4

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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

Fri, 09/14/2018 - 17:58

sealx017

Offline

Joined: 08/17/2018 - 21:55

That is great!!

Can we also add multiple equalities like,

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

Fri, 09/14/2018 - 18:11

(Reply to #5) #6

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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)

Fri, 09/14/2018 - 18:08

sealx017

Offline

Joined: 08/17/2018 - 21:55

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

Fri, 09/14/2018 - 18:16

(Reply to #7) #8

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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()?

Fri, 09/14/2018 - 19:24

sealx017

Offline

Joined: 08/17/2018 - 21:55

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

Sat, 09/15/2018 - 16:11

(Reply to #9) #10

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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.

Sat, 09/15/2018 - 18:25

#11

sealx017

Offline

Joined: 08/17/2018 - 21:55

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

Sat, 09/15/2018 - 19:49

(Reply to #11) #12

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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.

Sat, 09/15/2018 - 22:08

#13

sealx017

Offline

Joined: 08/17/2018 - 21:55

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.

Sun, 09/16/2018 - 16:00

(Reply to #13) #14

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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.

Sun, 09/16/2018 - 16:50

#15

sealx017

Offline

Joined: 08/17/2018 - 21:55

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.

Mon, 09/17/2018 - 11:34

(Reply to #15) #16

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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)

Mon, 09/17/2018 - 16:48

#17

sealx017

Offline

Joined: 08/17/2018 - 21:55

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.

Tue, 09/18/2018 - 12:37

(Reply to #17) #18

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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.

Tue, 09/18/2018 - 18:46

#19

sealx017

Offline

Joined: 08/17/2018 - 21:55

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

Wed, 09/19/2018 - 11:46

(Reply to #19) #20

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

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?

Main menu

Navigation

You are here

Multivariate cholesky model

Main menu

User login

Navigation

You are here

Multivariate cholesky model

Search form