Hi Everyone,
I am trying to fit a trivariate model since I think is the best option instead two separately bivariate models but I am not sure.
I have two errors:
1)
Warning messages:
- thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" )
1: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'values' in mxMatrix(type = "Full", nrow = nth, ncol = ntvo, free = TRUE, values = svTh, lbound = lbTh, labels = labTh, name = "thinG")
2: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'lbound' in mxMatrix(type = "Full", nrow = nth, ncol = ntvo, free = TRUE, values = svTh, lbound = lbTh, labels = labTh, name = "thinG")
2)
*fitACE <- mxRun( modelACE, intervals=T )
Error: In model 'twoACEja' the name 'mean_AG_Ln' is used as a free parameter in 'twoACEja.meanG', 'MZ.meanG', and 'DZ.meanG' and as a fixed parameter in 'twoACEja.meanG', 'MZ.meanG', and 'DZ.meanG'
And here the full script
# Select Continuous Variables varsc <- c('AG','RU') # list of continuous variables names nvc <- 2 # number of continuous variables ntvc <- nvc*2 # number of total continuous variables conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="") # Select Ordinal Variables nth <- 1 # number of thresholds varso <- c('SL') # list of ordinal variables names nvo <- 1 # number of ordinal variables ntvo <- nvo*2 # number of total ordinal variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") ordData <- twinData #for (i in ordVars) { ordData[,i] <- cut(twinData[,i], #breaks=quantile(ordData[,i],(0:(nth+1))/(nth+1),na.rm=TRUE), labels=c(0:nth)) } # Select Variables for Analysis vars <- c('AG','RU','SL') # list of variables names nv <- nvc+nvo # number of variables ntv <- nv*2 # number of total variables oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis #ordData[,'age'] <- ordData[,'age']/100 #ordData <- ordData[-which(is.na(ordData$age)),] covVars <- c('age', "Sex1" , "Sex2") nc <- 2 # number of covariates # number of covariates # Select Data for Analysis mzData <- subset(ordData, Zyg==1, c(selVars, covVars)) dzData <- subset(ordData, Zyg==2, c(selVars, covVars)) mzDataF <- mzData dzDataF <- dzData mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1)) mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1)) dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1)) dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1)) # Generate Descriptive Statistics apply(mzData,2,myMean) apply(dzData,2,myMean) # Set Starting Values frMV <- c(TRUE,FALSE) # free status for variables frCvD <- diag(frMV,ntv,ntv) # lower bounds for diagonal of covariance matrix frCvD[lower.tri(frCvD)] <- TRUE # lower bounds for below diagonal elements frCvD[upper.tri(frCvD)] <- TRUE # lower bounds for above diagonal elements frCv <- matrix(as.logical(frCvD),4) svMe <- c(1,1,0) # start value for means svPa <- .1 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .1 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix lbPa <- .0001 # start value for lower bounds lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements svLTh <- -1 # start value for first threshold svITh <- 1 # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds # Create Labels labMe <- paste("mean",vars,sep="_") labTh <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_") labBe <- labFull("beta",nc,1) # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12"), name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.Sex1", "data.Sex1","data.Sex1","data.Sex2","data.Sex2", "data.Sex2"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23","b21", "b22", "b23"), name="b2" ) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMe, name="meanG" ) expMean <- mxAlgebra( expression= meanG + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex , name="expMean" ) thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) threG <- mxAlgebra( expression= inc %*% thinG, name="threG" ) # Create Matrices for Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("a",nv), lbound=lbPaD, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), lbound=lbPaD, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" ) # Create Algebra for Variance Comptwonts covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= A+C+E, name="V" ) covMZ <- mxAlgebra( expression= A+C, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) # Create Algebra for Standardization matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Calculate genetic and environmental correlations corA <- mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="rA" ) #cov2cor() corC <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" ) corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" ) # Constrain Variance of Binary Variables matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" ) matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc" ) var1 <- mxConstraint( expression=diag2vec(Oc%&%V)==Unv1, name="Var1" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzDataF, type="raw" ) dataDZ <- mxData( observed=dzDataF, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(pathB1,pathB2, meanG, thinG, inc, threG, matI, invSD, matUnv, matOc, pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE) defs <- list(defAge,defSex) modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Variance Components rowVC <- rep('VC',nv) colVC <- rep(c('A','C','E','SA','SC','SE'),each=nv) estVC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC)) # Create Confidence Interval Objects ciACE <- mxCI( "VC[1,1]" )# "VC[1,seq(1,3*nv,nv),(2,seq(1,3*nv,nv)),(2,seq(2,3*nv,nv)))]" ) # Build Model with Confidence Intervals modelACE <- mxModel( "twoACEja", pars, var1, modelMZ, modelDZ, multi, estVC, ciACE ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=T ) sumACE <- summary( fitACE ) # Compare with Saturated Model mxCompare( fit, fitACE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEsts(fitACE) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Run AE model modelAE <- mxModel( fitACE, name="twoAEja" ) modelAE <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=T ) mxCompare( fitACE, fitAE ) fitGofs(fitAE) fitEsts(fitAE) # Run CE model modelCE <- mxModel( fitACE, name="twoCEja" ) modelCE <- omxSetParameters( modelCE, labels=labLower("a",nv), free=FALSE, values=0 ) fitCE <- mxRun( modelCE, intervals=T ) mxCompare( fitACE, fitCE ) fitGofs(fitCE) fitEsts(fitCE) # Run E model modelE <- mxModel( fitAE, name="twoEja" ) modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) mxCompare( fitAE, fitE ) fitGofs(fitE) fitEsts(fitE) # Print Comparative Fit Statistics mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )
Thank you so much I have tried everything but I can not make work it
Printer-friendly version
I have managed to make it work but I'm not sure because I have to put 2 thresholds otherwise it does not work, but my ordinal variable is dichotomous. is it correct?
Moreover, I am not sure if I have added the covariables properly.
Here the full script:
Thank you so much again.
Hi Juan
You may be on the right track, but if the variables are binary then you should need at most one threshold in the threshold matrix.
I couldn't run your script because the variable names don't agree with those in the dataset. Also labFull function is missing.
HTH
labFull()is from Hermine's miFunctions.R . I also assume that Juan is using his own data, and not the 'twinData' object that comes with OpenMx.The covariates look OK to me.
If the ordinal variable is dichotomous, you should only need 1 threshold, which allows you to simplify your script a bit. You also need to either fix the one threshold, or the regression intercept for the dichotomous variable, in order to identify the model. I'd say the simplest modification to make would be to try changing
to
, changing
to
, and changing
to
Thank you! I have made the changes but I have a warning message in this line:
thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, lbound=lbTh, labels=labTh, name="thinG" )
Warning message:
data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'lbound' in mxMatrix(type = "Full", nrow = 1, ncol = ntvo, free = FALSE, values = 0, lbound = lbTh, labels = labTh, name = "thinG")
Thank you so much!!
Sorry I forgot to say that I use my own data
Thanks
That should be benign, since the threshold is fixed. In fact, because it's fixed, you don't need to give it labels or a lower bound. If you want, you could instead do
Thank you so much it works perfect!
Here you can see my results
ACE
A C E --------------------rA rC rE-----------Rph
SL 0.34 0.35 0.31
AG 0.50 0.12 0.38
RB 0.48 0.22 0.30
SL&AG ---------------------------0.71(0.01,0.99) 0.84(0.121,1) -0.26(-0.56,0.27)------------- 0.38
SL&RB ---------------------------0.12(-0.45,0.67) 0.93(0.21,1) 0.05(-0.39,0.40)----------------0.32
AG&RB ---------------------------0.78(0.59,0.87) 0.98(0.77,1) 0.35(0.24,0.48)-----------------0.66
AE
A E---------------rA rE-------------Rph
SL 0.79 0.21
AG 0.62 0.38
RB 0.71 0.29
SL&AG-------------------0.61(0.46,0.75) -0.21(-0.27,0.07)-------0.37
SL&RB-------------------0.45(0.33,0.56) -0.08(-0.33,0.18)-------0.31
AG&RB-------------------0.81(0.77,0.86) 0.34(0.27,0.42)--------0.66
This is the comparison
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoACEja 27 11475.559 6073 -670.44092 NA NA NA
2 twoACEja twoAEja 21 11487.116 6079 -670.88425 11.556674 6 0.072621853
I am not sure what model should I chose and the criterion to choose one or the other
Also, I need to know if A, C, E are significant in each model. I mean for example in the ACE model for sl C is 0.35 and for ag 0.12 I would like to know if C is significant for SL and AG or only for one of them. I think E is always significant because contain error. I have thought calculate CI and if it contais 0 is not significative but I am not sure.
Finally, Can I compare the rA correlations in each trait (with statistic parameters)? I mean if the rA between 1-2 is bigger than 1-3…
Thank you so much in advance I really appreciate it.
Assuming you only want to report the results from a single "best" model, then I recommend selecting the model with the smallest AIC, and reporting confidence intervals for the quantities of interest. BTW, have you fit a saturated model? Is a CE model plausible for your data?
However, I advocate reporting conclusions based on multiple models; see, e.g., my relevant publication for more information.
So just add the needed character strings to the
referenceargument tomxCI(). For instance,will request CIs for the elements of the additive-genetic, shared-environmental, and nonshared-environmental covariance matrices, in addition to that for which your script already requests CIs.
It may or may not be, but inferences about E are usually not of interest.
Correct.
Yes. Create a suitable MxAlgebra and request confidence intervals for it, for instance,
, place it inside the MZ and DZ submodels, and add "MZ.rAcompare" to the vector of character strings provided for argument
referencetomxCI()when creating objectciACE. This will give you confidence intervals for the differences in genetic correlations.Once again thank you so much.
I have been working with the saturated model but I cannot make it works. I have tried all possible starting values but It doesn't find a solution. I don't know if I am doing something wrong here is the script.:
Thanks!!!
So, what happens when you try to run it? Also, what do you get from
mxVersion()?This is the output:
Begin fit attempt 31 of at maximum 32 tries
Lowest minimum so far: 12836.907293975
OpenMx status code 6 not in list of acceptable status codes, 0
OpenMx status code 6 not in list of acceptable status codes, 0
Not all eigenvalues of Hessian are greater than 0: 787961.505795225, 241712.973930221, 90801.8260774038, 37461.8509249169, 36108.3570058396, 20452.956895597, 16278.6263111995, 13514.5169971915, 12348.4836755664, 7865.0595383285, 7434.41254730365, 7351.21427400012, 6728.35738212829, 5648.20553748093, 4416.33288252297, 3674.66242643252, 3514.28595340669, 3224.37985610757, 3058.62507833283, 2824.51753405145, 2653.72023767685, 2290.52727017331, 2143.71360085118, 1845.17794007907, 1758.30696016368, 1351.05525266071, 1226.13780443413, 1173.91950744854, 867.931406722967, 807.202050148447, 615.874189695127, 502.241971723273, 418.459657188584, 364.354936896746, 36.5272816693136, 30.8297966297892, 25.5256524187929, 16.5921347221574, 12.4923808408779, 11.2477690283236, 5.48989588395323, 5.08592627903294, 4.88676405380493, 2.5967232615173, 2.08173667306169, 1.82647729161902, 1.52624731830076e-05, -8.33106327891148e-05
Begin fit attempt 32 of at maximum 32 tries
Fit attempt worse than current best: 12836.908910494 vs 12836.907293975
Retry limit reached
Start values from best fit:
0.10631708710754,0.0362139006101167,-0.0613654041393264,0.00569082513418116,-0.185581069609575,-0.344162428087438,1.06767187921709,-8.36632390547581,0.942344815595082,1.07312588323597,0.408909142438524,2.2005873612565,1.30673456966261,41.1843125633866,0.313146221049428,0.167776661842215,0.85467485365846,0.242455084108765,0.292201400080796,0.869137784815067,0.220708232633343,0.521510782740128,0.274754384528732,0.118846729885699,3.69004028161039,0.0849444409204678,0.103993475707172,1.09385843177197,-17.2859999506245,0.966720410717319,1.26641979849299,0.441386464784935,3.94902382148772,1.9668200201659,185.096253954615,0.164428802392192,0.108368824726336,1.08268432674425,0.219192741624143,0.173675315487782,2.13781627870702,0.227229733482567,0.578122306941328,0.197451965247457,0.0792066276546675,3.46936267863339,0.0862225357533247,0.088883979511605
And here my version:
OpenMx version: 2.7.9 [GIT v2.7.9]
R version: R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32
Default optimiser: NPSOL
Thank you so much!!!
Because your ordinal variable is dichotomous (has 2 levels), you need to fix its variance, and fix either its threshold or its regression intercept. Your script fixes the threshold, which is OK, but then the intercept needs to be free. Bearing in mind the value of
selVars, look atYour script only fixes SL's intercept for twin #2, and you fix the intercepts for twin #1's RU and twin #2's AG, both of which should be free. [Edit: to clarify, all the intercepts should be free if you're going to fix the threshold.] A similar situation holds for the diagonal elements of
covMZandcovDZ: only the SL variances should be fixed (and most people fix them to 1.0, but that's arbitrary).Aside from that, are the start values you're using reasonable? You already have estimates of the regression coefficients, and of the MZ and DZ covariance-matrix elements, from your other models. You could use those as start values. You could also run
cov()on your MZ and DZ data matrices, with argumentuse="pair".Since you have an ordinal variable, consider using
mxTryHardOrdinal(). Also, be advised that parameterizing covariance matrices as a symmetric matrix of free parameters doesn't always work very well with with the TryHard functions, since the matrices are not guaranteed to be positive-definite at the start of an attempt.Thanks!!!
I have made the changes but I am not sure If I have made the changes properly.
Now I think the saturated model is OK here the output.
Solution found
Running final fit, for Hessian and/or standard errors and/or confidence intervals
Warning messages generated from final fit for final fit for Hessian/SEs/CIs
Start values from best fit:
But in the Test Significance of Covariates I have this warning:
Running testCov with 58 parameters
Warning message:
In model 'testCov' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)
> mxCompare( fit, fitCov )
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 58 12559.294 6041 477.29376 NA NA NA
2 twoSATja testCov 58 12559.294 6041 477.29353 -0.00023061521 0 NA
Warning message:
In collectStatistics(nextBaseSummary, nextCompareSummary) :
Model 'twoSATja' has the same degrees of freedom as testCov which means the models are not nested and should not be compared with the likelihood ratio test. Compare these models using the information criteria statistics
>
And in the # Constrain expected Thresholds to be equal across twin order
Error in omxSetParameters(modelEMTVO, label = c(labThMZ[nvo * nth + 1], :
The following labels are not present in the model (use 'strict' = FALSE to ignore): 't1MZ_SLEEP2' and 't1MZ_SLEEP1'
The full script here:
Thank you again and sorry because I am really new and I can not understand some things
Your script still does not address the issue I raised about the fixed vs. free status of the variances. You need to construct
covMZandcovDZdifferently. Try this:The relevant syntax is:
I'm surprised it ran without
omxSetParameters()throwing an error. Your script never actually useslabBein the model. All the regression coefficients' labels start with "b", and not "beta_", as inlabBe. As a result, the call toomxSetParameters()doesn't actually change anything, andfitCovis justfit, re-run from its solution. That's why it has the same degrees of freedom. Change the call toomxSetParameters()to something like this:Your script doesn't give any labels to the matrix of thresholds. Try redefining them as:
To constrain thresholds across twin order (I think this will work):
To further constrain thresholds across zygosity:
Thank you so much!!!
I have made the changes. Also I had to change the optimizer to : mxOption(NULL,"Default optimizer","SLSQP") because with the normal optimizer I always got the same error:
Warning message:
In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, :
The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)
I don't know if it is ok.
Here you can see the results :
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA
2 twoSATja testCov 54 11479.071 6045 -610.92935 6.5282698e+01 6 3.7769626e-12
3 twoSATja eqThreshTwin 58 11413.788 6041 -668.21204 4.0857321e-06 2 9.9999796e-01
4 twoSATja eqThresZyg 58 11413.788 6041 -668.21205 -1.0619260e-08 2 1.0000000e+00
And the full script
I really appreciate it!
Something is still not quite right here. The fourth model, "eqThresZyg", should have one less free parameter than the third, "eqThreshTwin". Maybe run
omxGetParameters()on the model named "eqThresZyg" before you run it, to check whether my suggested syntax did what it was supposed to?You were right I had a mistake in my script:
This line was before mxModel
modelEMTVZ <- omxSetParameters(model=fitEMTVO, labels=c("tauMZ","tauDZ"), newlabels=c("tau","tau"), name="eqThresZyg")
Thanks!!
Here you can see the correct results saturated model and saturated vs ADE and ACE
> mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) )
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA
2 twoSATja testCov 54 11479.071 6045 -610.92935 6.5282698e+01 6 3.7769626e-12
3 twoSATja eqThreshTwin 58 11413.788 6041 -668.21205 1.7364073e-08 2 9.9999999e-01
4 twoSATja eqThresZyg 57 11413.788 6042 -670.21205 4.1691237e-09 3 1.0000000e+00
>
> mxCompare( fit, fitACE )
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA
2 twoSATja twoACEja 27 11475.559 6073 -670.44092 61.771126 34 0.0024790555
> mxCompare( fit, fitADE )
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja 60 11413.788 6039 -664.21205 NA NA NA
2 twoSATja twoADEja 27 11474.970 6073 -671.02971 61.182339 34 0.0028790296
And here ADE vs AE and ACE vs AE
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoADEja 27 11474.970 6073 -671.02971 NA NA NA
2 twoADEja twoAEja 21 11487.115 6079 -670.88451 12.145205 6 0.05880676
1 twoACEja 27 11475.559 6073 -670.44092 NA NA NA
2 twoACEja twoAEja 21 11487.116 6079 -670.88425 11.556674 6 0.072621853
So, what model is the best? Should I explain all as you mentioned before?
Also, I would like to know if it is possible:
1-Can I get CI for phenotypic correlations? I have been trying to get but I don't know how since phenotypic correlations are not in VC I mean I have got SA,SD and SC since they are in VC but I don't know how get the CI's phenotypic correlations since to get rPH y use cov2cor.
2-As you told me in order to know if A, C or D are significant I have calculated CI's for A C and D but for example in the ADE model for one variable D(unstandardized) is 0.05(0,0.15) I mean is weird that I don't have any negative value since I think 0.05 should be non-significant. I have a big sample but I don't know if I am doing something wrong.
3- Also, I would like to know where and how should I get the MZ and DZ correlations in the saturated model? with covariables? or in the ACE or ADE model? I mean I need the cross twin and cross trait correlations
Thank you so much. I really need it and I am learning a lot.
cov2cor()works in MxAlgebras, too. You can make MxAlgebras for the correlation matrices, e.g.and put them into the MZ and DZ submodels, respectively. Then, request confidence intervals for those two algebras; in your case, I think that would entail adding "MZ.corMZ" and "DZ.corDZ" to the character vector provided for argument
referencetomxCI(). You could also instead request CIs only for specific elements of the correlation matrices (e.g., "MZ.corMZ[1,2]"), since requesting CIs for all the elements is computationally inefficient, since a correlation matrix is symmetric and has a fixed diagonal. This approach should work for all of the models.Does "D" here refer to a path coefficient, or to a variance component? If it's a variance component, then you should not expect to see a negative lower confidence limit, since the Cholesky parameterization prevents variance components from being negative (they are sums of squares of diagonal elements of the lower-triangular Cholesky factor). The fact that OpenMx's profile-likelihood CIs respect constraints on the parameters, and can be asymmetric when necessary, is among their advantages.
Something I just noticed about the script in post #17: the dichotomous phenotype has free thresholds and free regression intercepts (in
meanMZandmeanDZ). One of those two things has to be fixed to identify the model. The "status RED" warning was probably a result of this unidentification. I suggest just fixing the threshold to zero (which would make constraining it equal across zygosity and twin order pointless).Thank you so much!
1-OK! but using this I have the correlations for MZ1 and MZ2 I mean this:
mxAlgebra 'corDZ'
$formula: cov2cor(covDZ)
$result:
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00000000 0.67553492 0.35730803 0.34876643 0.41805389 0.35158189
[2,] 0.67553492 1.00000000 0.30300303 0.33465014 0.46008228 0.26706220
[3,] 0.35730803 0.30300303 1.00000000 0.25237671 0.30618880 0.38199174
[4,] 0.34876643 0.33465014 0.25237671 1.00000000 0.64877935 0.32860511
[5,] 0.41805389 0.46008228 0.30618880 0.64877935 1.00000000 0.28502398
[6,] 0.35158189 0.26706220 0.38199174 0.32860511 0.28502398 1.00000000
dimnames: NULL
How can I get the total value for MZ and DZ?
2-I am refering to this value:
VC 0.7930 0.6586 0.5146 0.2176 -0.0449 0.0040 0.5668 0.1360 -0.0552 0.5027 0.8786 1.1104 0.1380 -0.0600 0.0087
VC 0.6586 0.5482 0.4233 -0.0449 0.0471 -0.1371 0.1360 0.2326 -0.0004 0.8786 0.6621 1.4809 -0.0600 0.0569 -0.4796
VC 0.5146 0.4233 0.3468 0.0040 -0.1371 0.4905 -0.0552 -0.0004 0.1627 1.1104 1.4809 0.3468 0.0087 -0.4796 0.4905
SE SE SE
VC 0.3593 0.1814 -0.1190
VC 0.1814 0.2809 -0.0013
VC -0.1190 -0.0013 0.1627
and I calculated its CI in order to know if it is significative am I wrong?
3-Also, I would like to know If I can fit a AE model for one variable and ACE for the other 2 variables. is it possible? I should fix C to 0 for one variable? how can I do this only for one variable?
4- for this:
Something I just noticed about the script in post #17: the dichotomous phenotype has free thresholds and free regression intercepts (in meanMZ and meanDZ). One of those two things has to be fixed to identify the model. The "status RED" warning was probably a result of this unidentification. I suggest just fixing the threshold to zero (which would make constraining it equal across zygosity and twin order pointless).
I had thought change this:
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
to :
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
But I am not sure if it would be correct.
Thank you so much again in adavance, I really appreciate your help.
Sorry, I don't understand the question?
That sounds good to me. The element you bolded is the dominance-genetic variance component for the second phenotype, so (given the parameterization you're using) its CI should not have a negative lower limit.
Yes, that's possible. The easiest way to do it would be to fix to zero the appropriate path coefficients in the common-environmental Cholesky factor (named "c" in your script), using
omxSetParameters(), in such a way that the phenotype has no common-environmental variance, nor common-environmental covariance with the other traits. I guess you could instead use MxConstraints to force the phenotype's variance and covariance elements to equal zero in the MxAlgebra named "C", but that is a clumsy approach. Be warned that if you drop the effect of C for only one phenotype, the MxAlgebra named "C" will no longer be a valid (positive-definite) covariance matrix, so the results of other MxAlgebras that use it may be rendered uninterpretable.That should work.
Thanks!
1-Sorry I meant the correlations for MZ and DZ I mean the correlations that we would use in the Falconer's formula
2-So, if its CI should not have a negative lower limit how can I know if D (or other parameter) is significant?
3-I have made the change to this :
modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0)
And It works !! :)
Thank you so much!!!
Hi again,
I have been working in the saturated model, but if I do this:
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
then Taumz1, Taumz2, TauDZ1 and Taudz2 are not in the model therefore the saturated model and eqThresTwin and eqThresZyg models have the same degrees of freedom so what changes should I do in these models?
Also, I have been thinking in de ACE/AE model about what you told me “I guess you could instead use MxConstraints to force the phenotype's variance and covariance elements to equal zero in the MxAlgebra named "C", but that is a clumsy approach. Be warned that if you drop the effect of C for only one phenotype, the MxAlgebra named "C" will no longer be a valid (positive-definite) covariance matrix, so the results of other MxAlgebras that use it may be rendered uninterpretable.” But, I am not sure how I should do a ACE/AE model I thought
modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0) but I'm not sure if it's both practically and theoretically correct
I made this decision because I think C in one trait is not significant 0.35(0,0.57)
But as I mentioned earlier. I'm not sure if I'm calculating well the intervals to know if it is significant or not. In order to do this I calculated the intervals for VC (VC[1,1], [2,2]…) I mean the unstandardized value. I don’t know if it is correct.
Thank you so much any insight you could give me about this or materials to learn more would be really helpful
?? I'm still confused as to which correlations you might want that aren't in the 'corMZ' and 'corDZ' MxAlgebras.
If the CI for a quantity contains zero, then that means that you cannot reject the null hypothesis that the quantity equals zero, i.e. it does not differ significantly from zero.
If you want to drop C from the third phenotype, then that looks reasonable to me.
Are the "eqThresTwin" and "eqThresZyg" models actually of interest? In any event, you could either (1) fix the intercepts of the ordinal phenotype to zero, and free the thresholds, or (2) leave the thresholds fixed, parameterize the ordinal phenotype's intercepts, in the saturated model, in terms of an MZ1 intercept, an MZ2 intercept, a DZ1 intercept, and a DZ2 intercept, and finally fit an analogous pair of additional models that reduce the four free intercept parameters to two (to equate across twin order) and reduce the two free parameters to one (to equate across zygosity).
Hi!
1- Sorry but my English is bad and my knowledge about Open Mx too ☹ If I use
corMZ <- mxAlgebra(cov2cor(covMZ), name="corMZ")
corDZ <- mxAlgebra(cov2cor(covDZ), name="corDZ")
Then I got this matrix:
dimnames: NULL
mxAlgebra 'expCovDZ'
$formula: rbind(cbind(V, cDZ), cbind(t(cDZ), V))
$result:
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.60428543 0.76958186 0.48280016 0.59664363 0.41088273 0.40671628
[2,] 0.76958186 0.84383543 0.29644796 0.41088273 0.38962653 0.26024648
[3,] 0.48280016 0.29644796 1.00000003 0.40671628 0.26024648 0.52439786
[4,] 0.59664363 0.41088273 0.40671628 1.60428543 0.76958186 0.48280016
[5,] 0.41088273 0.38962653 0.26024648 0.76958186 0.84383543 0.29644796
[6,] 0.40671628 0.26024648 0.52439786 0.48280016 0.29644796 1.00000003
I would like to know the correlation for MZ and DZ ( rMZ and rDZ) I am not sure if I should get from the saturated model or ACE model and with or without covariables. For example, rMZ could be 0.40 and for DZ 0.20. It should be a heritability about 40%. I don’t know if I am expressing properly
2- So, If I have 0 in one interval the parameter is non-significant, right? I expected for example intervals such as (-0.425, 0.100) I mean not exactly the zero
3- If I do :modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0)
Then the mode does not calculate rC between the other 2 parameters is it correct?
4- I am going to keep on working in the saturated model.
Thank you so much!
I'd say, get it from the saturated model, with covariates.
Right. If you really want, you could do a likelihood-ratio test by comparing a model in which the quantity of interest is free versus a model in which it is fixed to the null value (which is usually zero). You can usually fix the quantity to the null value via judicious use of
omxSetParameters(); occasionally, you might need to use an MxConstraint.I think so? Could you post the output that doesn't look right to you?
Thank you so much! yes of course:
> fitACE$algebras
Thanks in advance
What's happening is that, internally,
sqrt(I * C)isn't invertible, so no result is returned for MxAlgebra 'rC'. Try something likemxEval(cov2cor(C[1:2,1:2]),fitACE,T)if you want to see the shared-environmental correlation for the first two phenotypes.It works perfectly as always!
Thank you so much!! I really appreciate it!