Siblings

Posted on
No user picture. JuanJMV Joined: 07/20/2016
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 )

Replied on Fri, 05/11/2018 - 11:40
Picture of user. AdminRobK Joined: 01/24/2014

At minimum, change the appropriate lines to this instead:

# 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" )
pathTw <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="tw11", lbound=lbPa, name="tw" )

# 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" )
covTw <- mxAlgebra( expression=tw %*% t(tw), name="Tw" )

# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= A+C+E+Tw, name="V" )
covMZ <- mxAlgebra( expression= A+C+Tw, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%A+ C+Tw, 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" )

If you were to fix parameter "tw11" to zero, siblings and DZ twins would have equivalent model-expectations.

Replied on Sat, 05/12/2018 - 10:53
No user picture. JuanJMV Joined: 07/20/2016

In reply to by AdminRobK

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

# Select Continuous Variables
varsc <- c('DE_Ln',"AN_Ln", 'EXT_Ln') # list of continuous variables names
nvc <- 3 # 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('APN') # 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

# Select Variables for Analysis
vars <- c('DE_Ln',"AN_Ln", 'EXT_Ln','APN') # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,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

covVars <- c('age1',"age2", "Sex1" , "Sex2")
nc <- 4 # number of covariates # number of covariates

# Select Data for Analysis
mzData <- subset(ordData, Zyg==1| Zyg==3, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2 | Zyg==4| Zyg==5, c(selVars, covVars))
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))
# Generate Descriptive Statistics

# 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(2,4,2,0) # start value for means
svPa <- .5 # start value for path coefficient
svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix
svPe <- .5 # start value for path coefficient for e
svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix
lbPa <- .00001 # 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

# ADE Model
# Create Matrices for Covariates and linear Regression Coefficients
defAge <- mxMatrix( type="Full", nrow=1, ncol=8, free=FALSE, labels=c("data.age1", "data.age1", "data.age1","data.age1","data.age2", "data.age2", "data.age2","data.age2"), name="Age" )
pathB1 <- mxMatrix( type="Full", nrow=1, ncol=8, free=TRUE, values=.01, label=c("b11", "b12","b13","b14", "b11", "b12", "b13","b14"), name="b1" )

defSex <- mxMatrix( type="Full", nrow=1, ncol=8, free=FALSE, labels=c("data.Sex1", "data.Sex1","data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2","data.Sex2"), name="Sex" )
pathB2 <- mxMatrix( type="Full", nrow=1, ncol=8, free=TRUE, values=.01, label=c("b21", "b22","b23","b24", "b21", "b22", "b23","b24"), name="b2" )

# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMe, name="meanG" )
expMean <- mxAlgebra( expression= meanG + b1*Age + b2*Sex , name="expMean" )
thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, 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" )
pathD <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("d",nv), lbound=lbPaD, name="d" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" )
pathTw <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD,label=labLower("tw",nv), lbound=lbPaD, name="tw" )

# Create Algebra for Variance Comptwonts
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covD <- mxAlgebra( expression=d %*% t(d), name="D" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
covTw <- mxAlgebra( expression=tw %*% t(tw), name="TW" )

# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= A+D+E+TW, name="V" )
covMZ <- mxAlgebra( expression= A+D+TW, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%A+ 0.25%x%D+TW, name="cDZ" )
covSIB <- mxAlgebra( expression= 0.5%x%A+ 0.25%x%D, 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" )

# 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()
corD <- mxAlgebra( expression=solve(sqrt(I*D))%&%D, name ="rD" )
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" )
dataSIB <- mxData( observed=sibData, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars )
expSIB <- mxExpectationNormal( covariance="expCovSIB", means="expMean", dimnames=selVars, thresholds="thinG",threshnames=ordVars )
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars <- list(pathB1,pathB2, meanG, thinG, matI, invSD, matUnv, matOc,
pathA, pathD, pathE, covA, covD, covE, covP, corA, corD, corE, pathTw,covTw)
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','D','E',"TW",'SA','SD','SE',"STW"),each=nv)
estVC <- mxAlgebra( expression=cbind(A,D,E,TW,A/V,D/V,E/V,TW/V), name="VC", dimnames=list(rowVC,colVC))
# Create Confidence Interval Objects
ciADE <- mxCI(c("VC[1,1]", "MZ.rA", "MZ.rD", "MZ.rE","MZ.A","MZ.D","MZ.E","VC"))

# Build Model with Confidence Intervals
modelADE <- mxModel( "twoACEja", pars, var1, modelMZ, modelDZ,modelSIB, multi, estVC, ciADE )

# ------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitADE <- mxTryHardOrdinal( modelADE, intervals=F, extraTries = 31 )
sumADE <- summary( fitADE )

# Compare with Saturated Model
mxCompare( fit, fitADE )

# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitADE)
fitEsts(fitADE)

# ------------------------------------------------------------------------------
# RUN SUBMODELS
#Without TW

modelADENOTW <- mxModel( fitADE, name="twoADENOTWja" )
modelADENOTW <- omxSetParameters( modelADENOTW, labels=labLower("tw",nv), free=FALSE, values=0 )
fitADENOTW <- mxTryHardOrdinal( modelADENOTW, intervals=F, extraTries = 11 )
mxCompare( fitADE, fitADENOTW )

Replied on Mon, 05/14/2018 - 16:26
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by JuanJMV

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?

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.

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.

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.

Finally, can I know with this script the genetic/environmental influences shared by the four variables at the same time?

Yes, in the sense that it will estimate for you a 4x4 A, D, and E covariance matrix .

Replied on Tue, 05/15/2018 - 09:49
No user picture. JuanJMV Joined: 07/20/2016

In reply to by AdminRobK

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 NA

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

Replied on Tue, 05/15/2018 - 10:51
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by JuanJMV

Where can I find more information about that? more suggestions apart from Boulder’s materials and Neale’s book?

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

But 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?

Oh, I see. Then, it sounds like you want to fit a common-pathway and/or independent-pathway model.

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?

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.

Replied on Fri, 03/13/2020 - 16:45
No user picture. inooradd Joined: 03/11/2020

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

File attachments
Replied on Mon, 03/16/2020 - 13:55
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by inooradd

Error in value[rows[[i]], cols[[i]]] <- startValue :
incorrect number of subscripts on matrix

Please run `traceback()` after you see this error. What does `traceback()` give you?

Replied on Mon, 03/16/2020 - 19:48
No user picture. inooradd Joined: 03/11/2020

In reply to by AdminRobK

> 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)
Replied on Tue, 03/17/2020 - 20:57
No user picture. inooradd Joined: 03/11/2020

In reply to by AdminRobK

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?

Replied on Wed, 03/18/2020 - 10:45
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by inooradd

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.

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.

OK, I understand now. Just remember that if you use sex as a moderator, you need to also adjust for its main effect.

Would it just be easier to standardize values by age, and avoid this?

No, I think your current approach is good.

Replied on Thu, 05/14/2020 - 17:30
No user picture. inooradd Joined: 03/11/2020

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?

- Select Data for Analysis
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

- Select Data for Analysis
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?

Replied on Thu, 05/14/2020 - 17:54
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by inooradd

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.