Hello,
I am trying to fit a 5-variate model. But I am not sure about the best way to do it.
I have got these results from the univariate analyses
VAR1
ADE 0.34 0.10 0.57
AE 0.42 * 0.58
VAR2
ADE 0.26 0.16 0.58
AE 0.41 * 0.59
VAR3
ACE 0.29 0.13 0.58
AE 0.45 * 0.55
CE * 0.33 0.67
VAR4
ADE 0.23 0.17 0.59
AE 0.39 * 0.61
VAR5
ACE 0.35 0.05 0.60
AE 0.41 * 0.59
CE * 0.29 0.71
Therefore, I have:
AE model for Variable 1
AE model for Variable 2
ACE model for Variable 3
ADE model for Variable 4
AE model for Variable 5
I have fitted several models but I am not sure what model should I select.
5-ACE model: BIC= -310035.618; AIC=-10692.932
5-ADE model BIC= -310065.74; AIC=--10723.10
5-AE model BIC= -310144.380; AIC= -10703.607
*The comparison between ACE and AE Is non-significant therefore I can select the AE model but the comparison between ADE and AE is significant therefore I cannot select the nested model
I have also tried to fit a 5variables-ACE-ADE model (BIC=-310104.987; AIC=-10696.909) with (C for variables 3 and 5 and D for variables 1,2 and 4).
However, I am not sure if I adapted it correctly since I am getting much lower values for C (near 0) and D in comparison with the univariate results or 5-ACE or ADE models.
*The comparison between ADCE and AE is non-significant.
Here is my adapted script.
# Select Variables for Analysis vars <- c('VAR1','VAR2','VAR3','VAR4', "VAR5" ) # list of variables names nv <- 5 # 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('age', "Sex1" , "Sex2") nc <- 3 # number of covariates ordData <- twinData # 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)) # Set Starting Values svMe <- c(1.2,1.1,1.1, 1.2,1.2) # start value for means svPa <- .4 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .8 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix lbPa <- .0001 # start value for lower bounds lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements # Create Labels labMe <- paste("mean",vars,sep="_") labBe <- labFull("beta",nc,1) # ------------------------------------------------------------------------------ # PREPARE MODEL # ADCE Model defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=5, ncol=1, free=TRUE, values=.01, label=c("b11","b12", "b13","b14","b15"), name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=10, free=FALSE, labels=c("data.Sex1", "data.Sex1", "data.Sex1", "data.Sex1","data.Sex1", "data.Sex2","data.Sex2","data.Sex2", "data.Sex2", "data.Sex2"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=10, free=TRUE, values=.01, label=c("b21", "b22","b23","b24","b25", "b21", "b22", "b23","b24","b25"), 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 + cbind(t(b1%*%Age),t(b1%*%Age))+b2*Sex, name="expMean" ) # Create Matrices for Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("a",nv), lbound=lbPaD, name="a" ) pathC <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=c(F,F,T,F,T), values=c(0,0,0.3,0,0.3), label=c(NA,NA,"c11",NA,"c22"), name="c" ) pathD <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=c(T,T,F,T,F), values=c(0.3,0.3,0,0.3,0), label=c("d11","d22",NA,"d33",NA), name="d" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" ) # Create Algebra for Variance Comptwonts covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covD <- mxAlgebra( expression=d %*% t(d), name="D" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression=A+C+D+E, name="V" ) covMZ <- mxAlgebra( expression= A+C+D, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%A+C+0.25%x%D, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) # Create Algebra for Standardization matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Calculate genetic and environmental correlations corA <- mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="rA" ) #cov2cor() corC <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" ) corD <- mxAlgebra( expression=solve(sqrt(I*D))%&%D, name ="rD" ) corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(pathB1, pathB2, meanG, matI, invSD, pathA, pathC,pathD, pathE, covA, covC,covD, covE, covP, corA, corC,corD, corE) defs <- list(defAge, defSex) modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Variance Components rowVC <- rep('VC',nv) colVC <- rep(c('A','C','D','E','SA','SC','SD','SE'),each=nv) estVC <- mxAlgebra( expression=cbind(A,C,D,E,A/V,C/V,D/V,E/V), name="VC", dimnames=list(rowVC,colVC)) # Create Confidence Interval Objects ciADCE <- mxCI(c( "MZ.rA", "MZ.rC", "MZ.rE","MZ.A","MZ.C","MZ.E","VC")) # Build Model with Confidence Intervals modelADCE <- mxModel( "twoACDEca", pars, modelMZ, modelDZ, multi, estVC, ciADCE ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE Model fitADCE <- mxTryHard( modelADCE, intervals=F ) sumADCE <- summary( fitADCE ) # Compare with Saturated Model mxCompare( fit, fitADCE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitADCE) fitEsts(fitADCE) # ------------------------------------------------------------------------------ # RUN SUBMODELS modelACE <- mxModel( fitADCE, name="twoACEmodca" ) modelACE <- omxSetParameters(modelACE, labels=c("d11","d22","d33"),free=FALSE, values = 0) fitACE <- mxTryHard( modelACE, intervals=F ) mxCompare( fitADCE, fitACE ) fitGofs(fitACE) fitEsts(fitACE) # ------------------------------------------------------------------------------ # RUN SUBMODELS modelADE <- mxModel( fitADCE, name="twoADEmodca" ) modelADE <- omxSetParameters(modelADE, labels=c("c11","c22"),free=FALSE, values = 0) fitADE <- mxTryHard( modelADE, intervals=F ) mxCompare( fitADCE, fitADE ) fitGofs(fitADE) fitEsts(fitADE) # Run AE model modelAE <- mxModel( fitADCE, name="twoAEca" ) modelAE <- omxSetParameters(modelAE, labels=c("d11","d22","d33","c11","c22"),free=FALSE, values = 0) fitAE <- mxRun( modelAE, intervals=F ) mxCompare( fitADCE, fitAE ) fitGofs(fitAE) fitEsts(fitAE)