############################################################################################################## # FULL BIVARIATE MODERATION MODEL # ############################################################################################################## library(MASS) library(OpenMx) datmz <- read.csv("H:/datmz.csv") datdz <- read.csv("H:/datdz.csv") varnames=colnames(datmz)[3:6] # The dependent phenotypes m1 t1 m2 t2 # Big model MZ and DZ modmodelMZ=mxModel("modmodMZ", # Definition variables which we need to create the moderated paths (from M->T) mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm1"), name="m1" ) , mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm2"), name="m2" ) , # The means .. Not starting values are good. mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=0, labels=c("meanm","meant","meanm","meant"), name="expMean" ), # ------------------------------------------------------------------------ # ACE model moderator - cholesky elements. now defined as scalars (1x1), called lower for no good reason. # moderator parameters mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0m", name="A0M" ) , # Path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0m", name="C0M" ) , # Path coefficient C mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0m", name="E0M" ), # Path coefficient E # ACE model trait 0 trait parameters: main effects mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0t", name="A0T" ) , # Path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0t", name="C0T" ) , # Path coefficient C mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0t", name="E0T" ), # Path coefficient E # ACE model trait 1 moderator effects mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1t", name="A1T" ) , # Path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1t", name="C1T" ) , # Path coefficient C mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1t", name="E1T" ), # Path coefficient E # ACE model moderator-trait 0 moderated covariances; main effects mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0mt", name="A0MT" ) , # Path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0mt", name="C0MT" ) , # Path coefficient C mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0mt", name="E0MT" ), # Path coefficient E # ACE model moderator-trait 1 moderated covariances; moderation mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1mt", name="A1MT" ) , # Path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1mt", name="C1MT" ) , # Path coefficient C mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1mt", name="E1MT" ), # Path coefficient E # Assemble parameters into a Cholesky decomposition for twin 1 and twin 2 # Do this in three parts : A, C, and E. mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values= c(1,0,1,0,1,0,1,1,0,1),name="PsAmz"), mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values= c(1,0,1,0,1,0,1,1,0,1),name="PsC"), # Matrices generated to hold A and E computed Variance Components # This is a Cholesky decomposition of A for twin 1 and twin 2 (mz and dz) # Note that m1 appear in the top left part (twin 1) and m2 in the bottom right part (twin 2) # A chol mxAlgebra( rbind( cbind(A0M,0,0,0), cbind(A0MT+m1%x%A1MT,A0T+m1%x%A1T,0,0), cbind(0,0,A0M,0), cbind(0,0,A0MT+m2%x%A1MT,A0T+m2%x%A1T) ), name="chA"), # C chol mxAlgebra( rbind( cbind(C0M,0,0,0), cbind(C0MT+m1%x%C1MT,C0T+m1%x%C1T,0,0), cbind(0,0,C0M,0), cbind(0,0,C0MT+m2%x%C1MT,C0T+m2%x%C1T) ), name="chC"), # E chol mxAlgebra( rbind( cbind(E0M,0,0,0), cbind(E0MT+m1%x%E1MT,E0T+m1%x%E1T,0,0), cbind(0,0,E0M,0), cbind(0,0,E0MT+m2%x%E1MT,E0T+m2%x%E1T) ), name="chE"), # Get the expected covariance matrices mxAlgebra(chA%*%PsAmz%*%t(chA),name="Amz" ), # Variance component A mxAlgebra(chC%*%PsC%*%t(chC),name="C" ), # Variance component C mxAlgebra(chE%*%t(chE),name="E" ), # Variance component E # Assemble into phenotypic 4x4 covariance matrix of MZ and DZ mxAlgebra(Amz+C+E, name="expCovMZ" ), mxData( observed=datmz, type="raw" ), mxExpectationNormal(covariance="expCovMZ", means="expMean", dimnames=varnames), mxFitFunctionML(vector=FALSE) ) # Same as above modmodelDZ=mxModel("modmodDZ", # Definition variables mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm1"), name="m1" ) , mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm2"), name="m2" ) , # Means mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=0, labels=c("meanm","meant","meanm","meant"), name="expMean" ), # ACE model moderator mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0m", name="A0M" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0m", name="C0M" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0m", name="E0M" ), # ACE model trait 0 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0t", name="A0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0t", name="C0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0t", name="E0T" ), # ACE model trait 1 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1t", name="A1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1t", name="C1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1t", name="E1T" ), # ACE model moderator-trait 0 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0mt", name="A0MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0mt", name="C0MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0mt", name="E0MT" ), # ACE model moderator-trait 1 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1mt", name="A1MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1mt", name="C1MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1mt", name="E1MT" ), mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values= c(1,0,.5,0,1,0,.5,1,0,1),name="PsAdz"), mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values= c(1,0,1,0,1,0,1,1,0,1),name="PsC"), # Matrices generated to hold A and E computed Variance components # A chol mxAlgebra( rbind( cbind(A0M,0,0,0), cbind(A0MT+m1%x%A1MT,A0T+m1%x%A1T,0,0), cbind(0,0,A0M,0), cbind(0,0,A0MT+m2%x%A1MT,A0T+m2%x%A1T) ), name ="chA"), # C chol mxAlgebra( rbind( cbind(C0M,0,0,0), cbind(C0MT+m1%x%C1MT,C0T+m1%x%C1T,0,0), cbind(0,0,C0M,0), cbind(0,0,C0MT+m2%x%C1MT,C0T+m2%x%C1T) ), name ="chC"), # E chol mxAlgebra( rbind( cbind(E0M,0,0,0), cbind(E0MT+m1%x%E1MT,E0T+m1%x%E1T,0,0), cbind(0,0,E0M,0), cbind(0,0,E0MT+m2%x%E1MT,E0T+m2%x%E1T) ), name ="chE"), mxAlgebra(chA%*%PsAdz%*%t(chA),name="Adz" ), # Variance component A mxAlgebra(chC%*%PsC%*%t(chC),name="C" ), # Variance component C mxAlgebra(chE%*%t(chE),name="E" ), # Variance component E # Algebra for expected variance/covariance matrices in MZ & DZ twins # This should look familiar! Expected MZ covariance matrix as function of A and E mxAlgebra(Adz+C+E, name="expCovDZ" ), mxData( observed=datdz, type="raw" ), mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=varnames), mxFitFunctionML(vector=FALSE) ) modmodelA=mxModel("modmodel",modmodelMZ, modmodelDZ, mxFitFunctionMultigroup(c("modmodMZ","modmodDZ"))) modmodelAfit=mxRun(modmodelA) ### Constrain moderation parameters to 0 ### # Fix common A, C, E, and unique A and E moderation paths to 0 modmodelAf1=omxSetParameters(modmodelA, label=c("a1mt","c1mt","e1mt","a1t","e1t"), free=F, values=0) modmodelAfitf1=mxRun(modmodelAf1) mxCompare(modmodelAfit,modmodelAfitf1) # Fix all moderation paths to 0 modmodelAf2=omxSetParameters(modmodelA, label=c("a1mt","c1mt","e1mt","a1t","c1t","e1t"), free=F, values=0) modmodelAfitf2=mxRun(modmodelAf2) mxCompare(modmodelAfit,modmodelAfitf2)