Hello,
I am trying to fit a trivariate model with 3 continuous variables. I am also trying to fit the independent and common pathways models but I do not know if my script is correct.
My models have:
ADE: 27 parameters
Independent pathway model: 27 parameters
Common pathway model: 24 parameters
I think that number of parameters is correct but I am not sure. Also, I have negative values (for example dl_1_1, fl_2_1) and I do not know if it is correct.
Finally, I am not sure how to get the standardized results for this model. would be correct the square of each value? For al_1_1, dl_1_1 and el_1_1, would be square root?
Thank you so much!
Here the script:
# Select Variables for Analysis vars <- c('IN_Ln','CA_Ln','MF_Ln') # list of variables names nv <- 3 # 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)) # Generate Descriptive Statistics colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") # Set Starting Values svMe <- c(1.2,1.9,1.1) # 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 # ADE Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12", "b13"), name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.Sex1", "data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23", "b21", "b22", "b23"), 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" ) 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" ) # Create Algebra for Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) 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+D+E, name="V" ) covMZ <- mxAlgebra( expression= A+D, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%A+ 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() 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, pathD, pathE, covA, covD, covE, covP, corA, 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','D','E','SA','SD','SE'),each=nv) estVC <- mxAlgebra( expression=cbind(A,D,E,A/V,D/V,E/V), name="VC", dimnames=list(rowVC,colVC)) # Create Confidence Interval Objects ciADE <- mxCI(c( "MZ.rA", "MZ.rD", "MZ.rE","MZ.A","MZ.D","MZ.E","VC")) # Build Model with Confidence Intervals modelADE <- mxModel( "twoADEca", pars, modelMZ, modelDZ, multi, estVC, ciADE ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ADE Model fitADE <- mxTryHard( modelADE, intervals=F ) sumADE <- summary( fitADE ) # Compare with Saturated Model mxCompare( fit, fitADE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitADE) fitEsts(fitADE) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Test Significance of Covariates modelCov <- mxModel( fitADE, name="testCov" ) modelCov <- omxSetParameters( modelCov, label=c("b11","b12","b13","b21", "b22","b23"), free=FALSE, values=0 ) fitCov <- mxRun( modelCov ) mxCompare( fitADE, fitCov ) # Run AE model modelAE <- mxModel( fitADE, name="twoAEca" ) modelAE <- omxSetParameters( modelAE, labels=labLower("d",nv), free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=F ) mxCompare( fitADE, fitAE ) fitGofs(fitAE) fitEsts(fitAE) # Run E model modelE <- mxModel( fitAE, name="twoEca" ) modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) mxCompare( fitAE, fitE ) fitGofs(fitE) fitEsts(fitE) # Print Comparative Fit Statistics mxCompare( fitADE, nested <- list(fitCov, fitAE, fitE) ) #round(rbind(fitADE$VC$result,fitAE$VC$result,fitE$VC$result),4) fitADE$algebras mxEval(cov2cor(V), fitADE, T) fitAE$algebras mxEval(cov2cor(V), fitAE, T) #----------------------------------------------- #COMMON AND INDEPENDENT MODELS # Fit Independent Pathway ADE Model # ------------------------------------------------------------------------------ nf <- 1 # number of factors # Matrices ac, cc, and ec to store a, d, and e path coefficients for common factors pathAc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=labFull("ac",nv,nf), name="ac" ) pathDc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=labFull("dc",nv,nf), name="dc" ) pathEc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=labFull("ec",nv,nf), name="ec" ) # Matrices as, ds, and es to store a, d, and e path coefficients for specific factors pathAs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=4, labels=labDiag("as",nv), lbound=.00001, name="as" ) pathDs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=4, labels=labDiag("ds",nv), lbound=.00001, name="ds" ) pathEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=5, labels=labDiag("es",nv), lbound=.00001, name="es" ) # Matrices A, D, and E compute variance components covA <- mxAlgebra( expression=ac %*% t(ac) + as %*% t(as), name="A" ) covD <- mxAlgebra( expression=dc %*% t(dc) + ds %*% t(ds), name="D" ) covE <- mxAlgebra( expression=ec %*% t(ec) + es %*% t(es), name="E" ) # Create Model Objects for Multiple Groups pars <- list(pathB1, pathB2,meanG, matI, invSD, pathAc, pathDc, pathEc, pathAs, pathDs, pathEs, covA, covD, covE, covP, corA, corD, corE) 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") ) # Build & Run Model modelIP <- mxModel( "mulIPc", pars, modelMZ, modelDZ, multi ) fitIP <- mxRun( modelIP, intervals=F ) sumIP <- summary( fitIP ) mxCompare( fitADE, fitIP ) fitGofs(fitIP) fitEsts(fitIP) parameterSpecifications(fitIP) # Fit Common Pathway ACE Model # ------------------------------------------------------------------------------ nl <- 1 # Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s) pathAl <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("al",nl), lbound=.00001, name="al" ) pathDl <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("dl",nl), lbound=.00001, name="dl" ) pathEl <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("el",nl), lbound=.00001, name="el" ) # Matrix and Algebra for constraint on variance of latent phenotype covarLP <- mxAlgebra( expression=al %*% t(al) + dl %*% t(dl) + el %*% t(el), name="CovarLP" ) varLP <- mxAlgebra( expression= diag2vec(CovarLP), name="VarLP" ) unit <- mxMatrix( type="Unit", nrow=nl, ncol=1, name="Unit") varLP1 <- mxConstraint( expression=VarLP == Unit, name="varLP1") # Matrix f for factor loadings on latent phenotype pathFl <- mxMatrix( type="Full", nrow=nv, ncol=nl, free=TRUE, values=.2, labels=labFull("fl",nv,nl), name="fl" ) # Matrices A, C, and E compute variance components covA <- mxAlgebra( expression=fl %&% (al %*% t(al)) + as %*% t(as), name="A" ) covD <- mxAlgebra( expression=fl %&% (dl %*% t(dl)) + ds %*% t(ds), name="D" ) covE <- mxAlgebra( expression=fl %&% (el %*% t(el)) + es %*% t(es), name="E" ) # Create Model Objects for Multiple Groups pars <- list(pathB1, pathB2,meanG, matI, invSD, pathAl, pathDl, pathEl, covarLP, varLP, unit, pathFl, pathAs, pathDs, pathEs, covA, covD, covE, covP) 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") ) # Build & Run Model modelCP <- mxModel( "mulCPc", pars, varLP1, modelMZ, modelDZ, multi ) fitCP <- mxRun(modelCP, intervals=F ) sumCP <- summary( fitCP ) mxCompare( fitADE, fitCP ) parameterSpecifications(fitCP) # ------------------------------------------------------------------------------ sink() save.image(paste(filename,".Ri",sep=""))