Hello,
I am trying to do a univariate analysis with a binary variable. I am not sure how to treat siblings
I have created different groups MZ, DZ, SIB (which I think is the same as introduce SIB in the DZ group).
I have also checked if DZ and SIB can be equated through the saturated model and there are not significant differences when I equated means and variances between DZ and Sib
I would like to know if I can use the siblings in my analysis (correlations are a little bit different) or should I do something different to introduce siblings in my analysis?
I have read this thread "https://openmx.ssri.psu.edu/thread/4086" but I am not sure if I need to check the "twin effects" and what should I change in the script
Thank you so much
Here the script:
# Select Variables for Analysis vars <- 'APN' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis covVars <- c('age1',"age2", "Sex1" , "Sex2") nc <- 4 # Select Data for Analysis ordData <- twinData mzData <- subset(ordData, Zyg==1| Zyg==3, c(selVars, covVars)) dzData <- subset(ordData, Zyg==2 | Zyg==4| Zyg==5 , c(selVars, covVars)) sibData <- subset(ordData,Zyg==6| Zyg==7 | Zyg==8 , c(selVars, covVars)) mzDataF <- mzData dzDataF <- dzData mzDataF$APN1 <- mxFactor(mzDataF$APN1, levels =c(0,1)) mzDataF$APN2 <- mxFactor(mzDataF$APN2, levels =c(0,1)) dzDataF$APN1 <- mxFactor(dzDataF$APN1, levels =c(0,1)) dzDataF$APN2 <- mxFactor(dzDataF$APN2, levels =c(0,1)) sibData$APN1 <- mxFactor(sibData$APN1, levels =c(0,1)) sibData$APN2 <- mxFactor(sibData$APN2, levels =c(0,1)) # 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 # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.age1", "data.age2"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=.01, label=c("b11", "b12"), name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.Sex1","data.Sex2"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=.01, label=c("b21", "b22"), name="b2" ) # Create Algebra for expected Mean & Threshold Matrices meanG <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="meanG" ) expMean <- mxAlgebra( expression= meanG + b1*Age + b2*Sex, name="expMean" ) threG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svTh, labels="tob", name="threG" ) # Create Matrices for Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="a11", lbound=lbPa, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="c11", lbound=lbPa, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="e11", lbound=lbPa, name="e" ) # Create Algebra for 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" ) # 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" ) covSIB <- mxAlgebra( expression= 0.5%x%A+ C, name="cSIB" ) 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" ) expCovSIB <- mxAlgebra( expression= rbind( cbind(V, cSIB), cbind(t(cSIB), V)), name="expCovSIB" ) # Constrain Variance of Binary Variables var1 <- mxConstraint( expression=diag2vec(V)==1, name="Var1" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzDataF, type="raw" ) dataDZ <- mxData( observed=dzDataF, type="raw" ) dataSIB <- mxData( observed=sibData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG" ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG" ) expSIB <- mxExpectationNormal( covariance="expCovSIB", means="expMean", dimnames=selVars, thresholds="threG" ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(pathB1, pathB2, meanG, threG, pathA, pathC, pathE, covA, covC, covE, covP) 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 ) modelSIB <- mxModel( name="SIB", pars, defs, expMean, covSIB, expCovSIB, dataSIB, expSIB, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ","SIB") ) # 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" ) # Build Model with Confidence Intervals modelACE <- mxModel( "oneACEba", pars, var1, modelMZ, modelDZ,modelSIB, multi, estVC, ciACE ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=T ) sumACE <- summary( fitACE )
Printer-friendly version
At minimum, change the appropriate lines to this instead:
If you were to fix parameter "tw11" to zero, siblings and DZ twins would have equivalent model-expectations.
Thank you so much Rob.
I have used that in a univariate model and it works perfect!
I am also working with a tetravariate model (3 continuous variables and 1 binary).
I have made the changes and I have compared the ADE model (with TW effect) and ADE model (without TW effect) the comparison is non-significant p=0.96 and TW is less than 1% in all the variables
So, can I assume that there are not differences between sib and DZ and use Sib in my analysis?
I am also a little bit confuse about “Cholesky decomposition” and “correlated factor model” my results are from the Cholesky decomposition, right? How can I check the correlated factor model? Do I need a different script? I have searched it in Boulder materials but I am not completely sure.
Finally, can I know with this script the genetic/environmental influences shared by the four variables at the same time?
Thank you so much in advance!!!!
It would be more accurate to say that you find no evidence of twin-specific environmental effects, and therefore, you chose to treat full sibs like DZ twins in subsequent analyses.
Why do you ask? In most cases, the Cholesky and correlated-factor "models" are equivalent parameterizations of the same model. The only exception I can think of offhand is in the case of "non-scalar" or "qualitative" sex-limitation.
Yes, in the sense that it will estimate for you a 4x4 A, D, and E covariance matrix .
I thought they were different models (I am confused about it, sorry). Where can I find more information about that? more suggestions apart from Boulder’s materials and Neale’s book?
Here my results from the 4-trait script.
b11 b12 b13 b14 b21 b22 b23 -0.0249 0.0092 -0.0353 -0.0185 0.2984 0.4055 -0.0904 b24 mean_DE_Ln mean_AN_Ln mean_EXT_Ln mean_APN a_1_1 a_2_1 -0.0638 2.4840 3.6142 3.3896 -0.1613 0.7674 0.4689 a_3_1 a_4_1 a_2_2 a_3_2 a_4_2 a_3_3 a_4_3 0.7036 0.2898 0.3788 -0.1192 -0.1943 0.0115 -0.4787 a_4_4 e_1_1 e_2_1 e_3_1 e_4_1 e_2_2 e_3_2 0.0001 0.9260 0.3311 -0.1696 0.0531 0.5751 -0.0437 e_4_2 e_3_3 e_4_3 e_4_4 0.1980 0.8366 0.0591 0.7769 A A A A D D D D E E E E TW TW TW TW SA VC 0.5890 0.3598 0.5400 0.2224 0 0 0 0 0.8575 0.3066 -0.1571 0.0492 0 0 0 0 0.4072 VC 0.3598 0.3633 0.2848 0.0623 0 0 0 0 0.3066 0.4404 -0.0813 0.1315 0 0 0 0 0.5399 VC 0.5400 0.2848 0.5094 0.2215 0 0 0 0 -0.1571 -0.0813 0.7306 0.0318 0 0 0 0 1.4101 VC 0.2224 0.0623 0.2215 0.3508 0 0 0 0 0.0492 0.1315 0.0318 0.6492 0 0 0 0 0.8189 SA SA SA SD SD SD SD SE SE SE SE STW STW STW STW VC 0.5399 1.4101 0.8189 0 0 0 0 0.5928 0.4601 -0.4101 0.1811 0 0 0 0 VC 0.4520 1.3994 0.3214 0 0 0 0 0.4601 0.5480 -0.3994 0.6786 0 0 0 0 VC 1.3994 0.4108 0.8744 0 0 0 0 -0.4101 -0.3994 0.5892 0.1256 0 0 0 0 VC 0.3214 0.8744 0.3508 0 0 0 0 0.1811 0.6786 0.1256 0.6492 0 0 0 0 lbound estimate ubound twoAENOTWja.VC[1,1] 0.4639 0.5890 0.7224 MZ.rA[1,1] NA 1.0000 NA MZ.rA[2,1] 0.6836 0.7779 0.8612 MZ.rA[3,1] 0.9156 0.9858 0.9999 MZ.rA[4,1] 0.1929 0.4892 0.9797 MZ.rA[1,2] 0.6836 0.7779 0.8612 MZ.rA[2,2] NA 1.0000 NA MZ.rA[3,2] 0.5106 0.6620 0.8055 MZ.rA[4,2] -0.1524 0.1744 0.6370 MZ.rA[1,3] 0.9156 0.9858 0.9999 MZ.rA[2,3] 0.5106 0.6620 0.8055 MZ.rA[3,3] NA 1.0000 NA MZ.rA[4,3] 0.2093 0.5240 0.9994 MZ.rA[1,4] 0.1929 0.4892 0.9787 MZ.rA[2,4] -0.1524 0.1744 0.6370 MZ.rA[3,4] 0.2093 0.5240 0.9994 MZ.rA[4,4] NA 1.0000 NA MZ.rD[1,1] NA NaN NA MZ.rD[2,1] NA NaN NA MZ.rD[3,1] NA NaN NA MZ.rD[4,1] NA NaN NA MZ.rD[1,2] NA NaN NA MZ.rD[2,2] NA 0.0000 NA MZ.rD[3,2] NA 0.0000 NA MZ.rD[4,2] NA 0.0000 NA MZ.rD[1,3] NA NaN NA MZ.rD[2,3] NA 0.0000 NA MZ.rD[3,3] NA 0.0000 NA MZ.rD[4,3] NA 0.0000 NA MZ.rD[1,4] NA NaN NA MZ.rD[2,4] NA 0.0000 NA MZ.rD[3,4] NA 0.0000 NA MZ.rD[4,4] NA 0.0000 NA MZ.rE[1,1] NA 1.0000 NA MZ.rE[2,1] 0.4199 0.4989 0.5704 MZ.rE[3,1] -0.2885 -0.1984 -0.1036 MZ.rE[4,1] -0.1117 0.0659 0.2268 MZ.rE[1,2] 0.4199 0.4989 0.5704 MZ.rE[2,2] NA 1.0000 NA MZ.rE[3,2] -0.2442 -0.1433 -0.0373 MZ.rE[4,2] 0.0401 0.2459 0.4559 MZ.rE[1,3] -0.2885 -0.1984 -0.1036 MZ.rE[2,3] -0.2442 -0.1433 -0.0373 MZ.rE[3,3] NA 1.0000 NA MZ.rE[4,3] -0.1298 0.0462 0.2431 MZ.rE[1,4] -0.1117 0.0659 0.2267 MZ.rE[2,4] 0.0401 0.2459 0.4576 MZ.rE[3,4] -0.1298 0.0462 0.2431 MZ.rE[4,4] NA 1.0000 NA MZ.A[1,1] 0.4639 0.5890 0.7224 MZ.A[2,1] 0.2771 0.3598 0.4463 MZ.A[3,1] 0.4568 0.5400 0.6259 MZ.A[4,1] 0.0894 0.2224 0.3512 MZ.A[1,2] 0.2771 0.3598 0.4463 MZ.A[2,2] 0.2850 0.3633 0.4439 MZ.A[3,2] 0.2153 0.2848 0.3543 MZ.A[4,2] -0.0471 0.0623 0.1745 MZ.A[1,3] 0.4568 0.5400 0.6259 MZ.A[2,3] 0.2153 0.2848 0.3543 MZ.A[3,3] 0.3947 0.5094 0.6334 MZ.A[4,3] 0.0938 0.2215 0.3365 MZ.A[1,4] 0.0894 0.2224 0.3512 MZ.A[2,4] -0.0471 0.0623 0.1745 MZ.A[3,4] 0.0938 0.2215 0.3365 MZ.A[4,4] 0.0626 0.3508 0.6413 MZ.D[1,1] NA 0.0000 NA MZ.D[2,1] NA 0.0000 NA MZ.D[3,1] NA 0.0000 NA MZ.D[4,1] NA 0.0000 NA MZ.D[1,2] NA 0.0000 NA MZ.D[2,2] NA 0.0000 NA MZ.D[3,2] NA 0.0000 NA MZ.D[4,2] NA 0.0000 NA MZ.D[1,3] NA 0.0000 NA MZ.D[2,3] NA 0.0000 NA MZ.D[3,3] NA 0.0000 NA MZ.D[4,3] NA 0.0000 NA MZ.D[1,4] NA 0.0000 NA MZ.D[2,4] NA 0.0000 NA MZ.D[3,4] NA 0.0000 NA MZ.D[4,4] NA 0.0000 NA MZ.E[1,1] 0.7506 0.8575 0.9762 MZ.E[2,1] 0.2394 0.3066 0.3812 MZ.E[3,1] -0.2214 -0.1571 -0.0860 MZ.E[4,1] -0.0798 0.0492 0.1704 MZ.E[1,2] 0.2394 0.3066 0.3812 MZ.E[2,2] 0.3776 0.4404 0.5133 MZ.E[3,2] -0.1381 -0.0813 -0.0216 MZ.E[4,2] 0.0215 0.1315 0.2399 MZ.E[1,3] -0.2214 -0.1571 -0.0860 MZ.E[2,3] -0.1381 -0.0813 -0.0216 MZ.E[3,3] 0.6327 0.7306 0.8392 MZ.E[4,3] -0.0944 0.0318 0.1627 MZ.E[1,4] -0.0798 0.0492 0.1704 MZ.E[2,4] 0.0215 0.1315 0.2399 MZ.E[3,4] -0.0944 0.0318 0.1627 MZ.E[4,4] 0.3587 0.6492 NA twoAENOTWja.VC[2,1] 0.2771 0.3598 0.4463 twoAENOTWja.VC[3,1] 0.4568 0.5400 0.6259 twoAENOTWja.VC[4,1] 0.0894 0.2224 0.3512 twoAENOTWja.VC[1,2] 0.2771 0.3598 0.4463 twoAENOTWja.VC[2,2] 0.2850 0.3633 0.4439 twoAENOTWja.VC[3,2] 0.2153 0.2848 0.3543 twoAENOTWja.VC[4,2] -0.0471 0.0623 0.1745 twoAENOTWja.VC[1,3] 0.4568 0.5400 0.6259 twoAENOTWja.VC[2,3] 0.2153 0.2848 0.3543 twoAENOTWja.VC[3,3] 0.3947 0.5094 0.6334 twoAENOTWja.VC[4,3] 0.0938 0.2215 0.3365 twoAENOTWja.VC[1,4] 0.0894 0.2224 0.3512 twoAENOTWja.VC[2,4] -0.0471 0.0623 0.1745 twoAENOTWja.VC[3,4] 0.0938 0.2215 0.3365 twoAENOTWja.VC[4,4] 0.0626 0.3508 0.6413 twoAENOTWja.VC[1,5] NA 0.0000 NA twoAENOTWja.VC[2,5] NA 0.0000 NA twoAENOTWja.VC[3,5] NA 0.0000 NA twoAENOTWja.VC[4,5] NA 0.0000 NA twoAENOTWja.VC[1,6] NA 0.0000 NA twoAENOTWja.VC[2,6] NA 0.0000 NA twoAENOTWja.VC[3,6] NA 0.0000 NA twoAENOTWja.VC[4,6] NA 0.0000 NA twoAENOTWja.VC[1,7] NA 0.0000 NA twoAENOTWja.VC[2,7] NA 0.0000 NA twoAENOTWja.VC[3,7] NA 0.0000 NA twoAENOTWja.VC[4,7] NA 0.0000 NA twoAENOTWja.VC[1,8] NA 0.0000 NA twoAENOTWja.VC[2,8] NA 0.0000 NA twoAENOTWja.VC[3,8] NA 0.0000 NA twoAENOTWja.VC[4,8] NA 0.0000 NA twoAENOTWja.VC[1,9] 0.7506 0.8575 0.9762 twoAENOTWja.VC[2,9] 0.2394 0.3066 0.3812 twoAENOTWja.VC[3,9] -0.2214 -0.1571 -0.0860 twoAENOTWja.VC[4,9] -0.0798 0.0492 0.1704 twoAENOTWja.VC[1,10] 0.2394 0.3066 0.3812 twoAENOTWja.VC[2,10] 0.3776 0.4404 0.5133 twoAENOTWja.VC[3,10] -0.1381 -0.0813 -0.0216 twoAENOTWja.VC[4,10] 0.0215 0.1315 0.2399 twoAENOTWja.VC[1,11] -0.2214 -0.1571 -0.0860 twoAENOTWja.VC[2,11] -0.1381 -0.0813 -0.0216 twoAENOTWja.VC[3,11] 0.6327 0.7306 0.8392 twoAENOTWja.VC[4,11] -0.0944 0.0318 0.1627 twoAENOTWja.VC[1,12] -0.0798 0.0492 0.1704 twoAENOTWja.VC[2,12] 0.0215 0.1315 0.2399 twoAENOTWja.VC[3,12] -0.0944 0.0318 0.1627 twoAENOTWja.VC[4,12] 0.3587 0.6492 NA twoAENOTWja.VC[1,13] NA 0.0000 NA twoAENOTWja.VC[2,13] NA 0.0000 NA twoAENOTWja.VC[3,13] NA 0.0000 NA twoAENOTWja.VC[4,13] NA 0.0000 NA twoAENOTWja.VC[1,14] NA 0.0000 NA twoAENOTWja.VC[2,14] NA 0.0000 NA twoAENOTWja.VC[3,14] NA 0.0000 NA twoAENOTWja.VC[4,14] NA 0.0000 NA twoAENOTWja.VC[1,15] NA 0.0000 NA twoAENOTWja.VC[2,15] NA 0.0000 NA twoAENOTWja.VC[3,15] NA 0.0000 NA twoAENOTWja.VC[4,15] NA 0.0000 NA twoAENOTWja.VC[1,16] NA 0.0000 NA twoAENOTWja.VC[2,16] NA 0.0000 NA twoAENOTWja.VC[3,16] NA 0.0000 NA twoAENOTWja.VC[4,16] NA 0.0000 NA twoAENOTWja.VC[1,17] 0.3290 0.4072 0.4820 twoAENOTWja.VC[2,17] 0.4307 0.5399 0.6405 twoAENOTWja.VC[3,17] 1.2078 1.4101 1.6595 twoAENOTWja.VC[4,17] 0.4146 0.8189 1.3582 twoAENOTWja.VC[1,18] 0.4307 0.5399 0.6405 twoAENOTWja.VC[2,18] 0.3634 0.4520 0.5325 twoAENOTWja.VC[3,18] 1.0977 1.3994 1.8025 twoAENOTWja.VC[4,18] -0.2987 0.3214 0.8795 twoAENOTWja.VC[1,19] 1.2078 1.4101 1.6595 twoAENOTWja.VC[2,19] 1.0977 1.3994 1.8025 twoAENOTWja.VC[3,19] 0.3270 0.4108 0.4917 twoAENOTWja.VC[4,19] 0.3981 0.8744 1.4805 twoAENOTWja.VC[1,20] 0.4146 0.8189 1.3582 twoAENOTWja.VC[2,20] -0.2987 0.3214 0.8795 twoAENOTWja.VC[3,20] 0.3981 0.8744 1.4805 twoAENOTWja.VC[4,20] NA 0.3508 0.6413 twoAENOTWja.VC[1,21] NA 0.0000 NA twoAENOTWja.VC[2,21] NA 0.0000 NA twoAENOTWja.VC[3,21] NA 0.0000 NA twoAENOTWja.VC[4,21] NA 0.0000 NA twoAENOTWja.VC[1,22] NA 0.0000 NA twoAENOTWja.VC[2,22] NA 0.0000 NA twoAENOTWja.VC[3,22] NA 0.0000 NA twoAENOTWja.VC[4,22] NA 0.0000 NA twoAENOTWja.VC[1,23] NA 0.0000 NA twoAENOTWja.VC[2,23] NA 0.0000 NA twoAENOTWja.VC[3,23] NA 0.0000 NA twoAENOTWja.VC[4,23] NA 0.0000 NA twoAENOTWja.VC[1,24] NA 0.0000 NA twoAENOTWja.VC[2,24] NA 0.0000 NA twoAENOTWja.VC[3,24] NA 0.0000 NA twoAENOTWja.VC[4,24] NA 0.0000 NA twoAENOTWja.VC[1,25] 0.5180 0.5928 0.6710 twoAENOTWja.VC[2,25] 0.3595 0.4601 0.5693 twoAENOTWja.VC[3,25] -0.6595 -0.4101 -0.2078 twoAENOTWja.VC[4,25] -0.3582 0.1811 0.6434 twoAENOTWja.VC[1,26] 0.3595 0.4601 0.5693 twoAENOTWja.VC[2,26] 0.4675 0.5480 0.6367 twoAENOTWja.VC[3,26] -0.8025 -0.3994 -0.0977 twoAENOTWja.VC[4,26] 0.1205 0.6786 1.2987 twoAENOTWja.VC[1,27] -0.6595 -0.4101 -0.2078 twoAENOTWja.VC[2,27] -0.8025 -0.3994 -0.0977 twoAENOTWja.VC[3,27] 0.5083 0.5892 0.6730 twoAENOTWja.VC[4,27] -0.4805 0.1256 0.6019 twoAENOTWja.VC[1,28] -0.3582 0.1811 0.6434 twoAENOTWja.VC[2,28] 0.1205 0.6786 1.2987 twoAENOTWja.VC[3,28] -0.4805 0.1256 0.6019 twoAENOTWja.VC[4,28] 0.3587 0.6492 NA twoAENOTWja.VC[1,29] NA 0.0000 NA twoAENOTWja.VC[2,29] NA 0.0000 NA twoAENOTWja.VC[3,29] NA 0.0000 NA twoAENOTWja.VC[4,29] NA 0.0000 NA twoAENOTWja.VC[1,30] NA 0.0000 NA twoAENOTWja.VC[2,30] NA 0.0000 NA twoAENOTWja.VC[3,30] NA 0.0000 NA twoAENOTWja.VC[4,30] NA 0.0000 NA twoAENOTWja.VC[1,31] NA 0.0000 NA twoAENOTWja.VC[2,31] NA 0.0000 NA twoAENOTWja.VC[3,31] NA 0.0000 NA twoAENOTWja.VC[4,31] NA 0.0000 NA twoAENOTWja.VC[1,32] NA 0.0000 NA twoAENOTWja.VC[2,32] NA 0.0000 NA twoAENOTWja.VC[3,32] NA 0.0000 NA twoAENOTWja.VC[4,32] NA 0.0000 NABut I would like to know the genetic influences (also environmental influences) shared between the four variables. I mean I can know it for v1-v2, v1-v3,…… but can I know the shared genetic influences between v1, v2, v3 and v4 together?
By the way since I have a negative value for example in MZ.E[3,1] I can not know the % of the phenotypic correlation explained by genetic factors right?
Thank you so much again!!
Hmm. This is the kind of thing one typically learns from one's dissertation adviser and/or at behavior-genetics workshops (like the one in Boulder). Offhand, I can think of two articles that are somewhat relevant: Loehlin (1996) and Neale, Røysamb, & Jacobson (2006).
Oh, I see. Then, it sounds like you want to fit a common-pathway and/or independent-pathway model.
You can still estimate its contribution to the phenotypic correlation, but interpreting it as a proportion doesn't really make sense because of the negative value.
Thank you so much Rob!
I am going to work on the common-pathway and independent-pathway models.
Hello,
I am new to OpenMX, R, and ACTE modeling and am seeking feedback on my attempt had performing something similar to the OP.
In my instance, I am trying to do a univariate analysis with a continuous variable. My sample contains MZ, DZ, and siblings, and my siblings contain more than 2; in some cases they have 4 - 6 siblings. I followed the information below and incorporated siblings into my syntax. However, I received an error and I am unable to detect the source.
Any help would be greatly appreciated. Thank you!
> # Run ACE Model
> fitACE <- mxRun( modelACE, intervals=T )
Error in value[rows[[i]], cols[[i]]] <- startValue :
incorrect number of subscripts on matrix
> sumACE <- summary( fitACE )
Error in summary(fitACE) : object 'fitACE' not found
Please run
traceback()after you see this error. What doestraceback()give you?> traceback()
20: populateDefVarMatrix(matrix, model, defvar.row)
19: computeMatrix(matrix, model, labelsData, defvar.row, env, cache)
18: evaluateMatrix(lookup, model, labelsData, env, show, compute,
defvar.row, cache)
17: evaluateSymbol(formula, contextString, model, labelsData, env,
compute, show, outsideAlgebra, defvar.row, cache)
16: evaluateExpression(formula[[i]], contextString, model, labelsData,
env, compute, show, outsideAlgebra, defvar.row, cache)
15: evaluateExpression(formula[[i]], contextString, model, labelsData,
env, compute, show, outsideAlgebra, defvar.row, cache)
14: evaluateExpression(formula[[i]], contextString, model, labelsData,
env, compute, show, outsideAlgebra, defvar.row, cache)
13: evaluateExpression(algebra@formula, contextString, model, labelsData,
env, compute = TRUE, show, outsideAlgebra = FALSE, defvar.row,
cache)
12: computeAlgebra(algebra, model, labelsData, show = FALSE, defvar.row,
env, cache)
11: evaluateAlgebra(lookup, model, labelsData, env, show, compute,
defvar.row, cache)
10: evaluateSymbol(as.symbol(objname), objname, flatModel, labelsData,
globalenv(), compute = TRUE, show = FALSE, outsideAlgebra = FALSE,
defvar.row = 1, cache = cache)
9: eval(substitute(evaluateSymbol(x, objname, flatModel, labelsData,
globalenv(), compute = TRUE, show = FALSE, outsideAlgebra = FALSE,
defvar.row = 1, cache = cache), list(x = quote(as.symbol(objname)))))
8: eval(substitute(evaluateSymbol(x, objname, flatModel, labelsData,
globalenv(), compute = TRUE, show = FALSE, outsideAlgebra = FALSE,
defvar.row = 1, cache = cache), list(x = quote(as.symbol(objname)))))
7: evaluateMxObject(meansName, flatModel, labelsData, new.env(parent = emptyenv()))
6: updateExpectationDimnames(.Object, flatModel, labelsData)
5: genericExpConvertEntities(expectations[[i]], flatModel, namespace,
labelsData)
4: genericExpConvertEntities(expectations[[i]], flatModel, namespace,
labelsData)
3: expectationFunctionConvertEntities(flatModel, namespace, labelsData)
2: runHelper(model, frontendStart, intervals, silent, suppressWarnings,
unsafe, checkpoint, useSocket, onlyFrontend, useOptimizer,
beginMessage)
1: mxRun(modelACE, intervals = T)
OK. What do you get from
mxVersion()?Also, does your dataset really contain columns named "age1" thru "age5" and "sex1" thru "sex5"?
Here you go.
OpenMx version: 2.14.11 [GIT v2.14.11-dirty]
R version: R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32
Default optimizer: CSOLNP
NPSOL-enabled?: Yes
OpenMP-enabled?: No
And yes, I have my covariances sex and age labeled for each child. Since the maximum number of children is 6, i created a label for each one, including sex1...sex6 and age1...age6. However, as per my coauthor, we will be moderating by sex, so I will only be controlling for age now.
Would it just be easier to standardize values by age, and avoid this?
Scratch that. We wont be moderating by sex, and will still control for sex and age still.
Your OpenMx version is 2 releases out-of-date. I think the first error message you're seeing is due to a bug, which has hopefully been fixed since version 2.14. First, try updating OpenMx, and see if the problem persists.
OK, I understand now. Just remember that if you use sex as a moderator, you need to also adjust for its main effect.
No, I think your current approach is good.
Hello,
I was able to get everything up and running. Turns out it was how I was loading my data; it should be uploaded as .DAT and as a matrix, not as excel.
Right off the bat, I tested twin environment effects by dropping the T, and while it was not significantly different, the AIC slightly improved, so I can treat DZ and siblings the same.
Question regarding that. Does that mean in the beginning, when we are creating distinct subsets of data for MZ, DZ, and siblings, do we code siblings as DZ?
conData <- twindata
mzData <- subset(conData, zyg==1, c(selVars, covVars))
dzData <- subset(conData, zyg==2 | zyg==3, c(selVars, covVars))
sibData <- subset(conData, zyg==4 | zyg==5, c(selVars, covVars))
becomes
conData <- twindata
mzData <- subset(conData, zyg==1, c(selVars, covVars))
dzData <- subset(conData, zyg==2 | zyg==3 | zyg==4 | zyg==5, c(selVars, covVars))
If this is the case, and we have instances of more than 2 siblings, would this matter?
I don't see why you'd need to change your datasets. Once you drop the T source of variance, then your model should in effect model full sibs the same way it models DZ twins, except that the full sibs in some families will be more like DZ triplets, quadruplets, quintuplets, or sextuplets. Granted, the full sibs in a family won't all have the same age the way both twins in a pair do, but that's why we adjust for the main effects of age.
Wonderful! I figured this was the case, but I wanted to make sure. I will just move forward with ACE rather than ACTE. Thank you!!