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 )