You are here

Siblings

17 posts / 0 new
Last post
JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Siblings

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 )
AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
twin effects

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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much!

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 )
AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
4-trait script
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 .

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you!

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

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
Where can I find more
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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much!!

Thank you so much Rob!

I am going to work on the common-pathway and independent-pathway models.

inooradd's picture
Offline
Joined: 03/11/2020 - 11:42
Siblings 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

File attachments: 
AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
traceback()
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?

inooradd's picture
Offline
Joined: 03/11/2020 - 11:42
re: traceback()

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

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
OK. What do you get from

OK. What do you get from mxVersion()?

Also, does your dataset really contain columns named "age1" thru "age5" and "sex1" thru "sex5"?

inooradd's picture
Offline
Joined: 03/11/2020 - 11:42
Here you go.

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?

inooradd's picture
Offline
Joined: 03/11/2020 - 11:42
Controlling for sex and age

Scratch that. We wont be moderating by sex, and will still control for sex and age still.

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
update OpenMx

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.

inooradd's picture
Offline
Joined: 03/11/2020 - 11:42
Treating DZ as siblings

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?

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
why do that?

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.

inooradd's picture
Offline
Joined: 03/11/2020 - 11:42
Wonderful! I figured this was

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