require(OpenMx); require(foreign); require(MASS); require(ggplot2); summaryACEFit = function(fit, selVars, accuracy=2) { nVar = length(selVars)/2; A = mxEval(z.A, fit); C = mxEval(z.C, fit); E = mxEval(z.E, fit); logLikelihood = mxEval(objective, fit); # Calculate standardised variance components Vtot = (A+C+E); a2 = A/Vtot; a2[upper.tri(a2)]=NA; c2 = C/Vtot; c2[upper.tri(c2)]=NA; e2 = E/Vtot; e2[upper.tri(e2)]=NA; aceEst = data.frame(cbind(a2,c2,e2), row.names=selVars[1:nVar]); names(aceEst)= paste(rep(c("A", "C", "E"), each=nVar), rep(1:nVar), sep=""); print(round(aceEst,accuracy)); print(logLikelihood[1,1]); } # read in a file called twinData data(twinData) selVars = c("wt1", "ht1", "wt2", "ht2") nVar = length(selVars)/2; # number of dependent variables ** per INDIVIDUAL ( so times-2 for a family)** mzData = as.matrix(subset(twinData, zyg==1, selVars)) dzData = as.matrix(subset(twinData, zyg==3, selVars)) obsMZmeans = colMeans(mzData,na.rm=TRUE); meanStarts = rep(obsMZmeans[1:nVar],2); meanLabels = paste("Trait", rep(1:nVar,2), "mean", sep=""); # make labels for the means matrix: same label for each trait in each twin ##### Fit ACE Model ##### share = mxModel("z", mxMatrix("Lower", nrow=nVar, ncol=nVar, free=TRUE, values=.6, name="a", byrow=TRUE), # Additive genetic path coefficient mxMatrix("Lower", nrow=nVar, ncol=nVar, free=TRUE, values=.6, name="c", byrow=TRUE), # Common environmental path coefficient mxMatrix("Lower", nrow=nVar, ncol=nVar, free=TRUE, values=.6, name="e", byrow=TRUE), # Unique environmental path coefficient # Multiply by each path coefficient by its inverse to get variance component mxAlgebra(a %*% t(a), name="A"), # additive genetic variance mxAlgebra(c %*% t(c), name="C"), # common environmental variance mxAlgebra(e %*% t(e), name="E"), # unique environmental variance mxMatrix("Full", nrow=1, ncol=(nVar*2), free=TRUE, values=meanStarts, label=meanLabels, dimnames=list( "means",selVars), name="expMean") ) mzModel = mxModel(name="MZ", # MZ ACE algebra mxAlgebra(rbind (cbind(z.A+z.C+z.E, z.A+z.C), cbind(z.A+z.C, z.A+z.C+z.E) ), dimnames = list(selVars, selVars), name="expCov"), mxData(mzData, type="raw"), mxFIMLObjective("expCov", "z.expMean") ) dzModel = mxModel(name="DZ", mxAlgebra(rbind (cbind(z.A+z.C+z.E , .5%x%z.A+z.C), cbind(.5%x%z.A+z.C, z.A+z.C+z.E) ), dimnames = list(selVars, selVars), name="expCov"), mxData(dzData, type="raw"), mxFIMLObjective("expCov", "z.expMean") ) #Build and Run ACE model model = mxModel("twinACE",share, mzModel, dzModel, mxAlgebra(MZ.objective + DZ.objective, name="twin"), mxAlgebraObjective("twin") ) fit = mxRun(model); summaryACEFit(fit, selVars); # so far so good #Now try and drop e newShare = mxModel(share, mxMatrix("Lower", nrow=nVar, ncol=nVar, free=FALSE, value=0, name="e", byrow=TRUE)) # drop e at 0 model = mxModel("twinAC",newShare, mzModel, dzModel, mxAlgebra(MZ.objective + DZ.objective, name="twin"), mxAlgebraObjective("twin") ) fit = mxRun(model); # ERROR in paste("The job for model", omxQuotes(flatModel@name), "exited abnormally with the error message:", : # cannot coerce type 'char' to vector of type 'character' # Now just switch to drop "c" instead of "e" newShare <- mxModel(share, mxMatrix("Lower", nrow=nVar, ncol=nVar, free=FALSE, value=0, name="c", byrow=TRUE)) # drop e at 0 model = mxModel("twinAE",newShare, mzModel, dzModel, mxAlgebra(MZ.objective + DZ.objective, name="twin"), mxAlgebraObjective("twin") ) fit = mxRun(model); # Hmm.. that worked fine apparently... summaryACEFit(fit, selVars); # now for the other strange bit: let's just try and drop one path in the 'a' matrix...say the second Cholesky latent factor # Build and Run reduced ACE model drop a path in 'a' newShare = mxModel(share, mxMatrix("Lower", nrow=nVar, ncol=nVar, free=c(T,T,F), values=c(.7,.7,0), name="a", byrow=TRUE) ) # drop 2:2 at 0 model = mxModel("twinCE",newShare, mzModel, dzModel, mxAlgebra(MZ.objective + DZ.objective, name="twin"), mxAlgebraObjective("twin") ) fit = mxRun(model); summaryACEFit(fit, selVars); # runs fine, but the path was not dropped. Some values now more than 1 or less than 0 # A1 A2 C1 C2 E1 E2 # wt1 0.64 NA 0.14 NA 0.22 NA # ht1 1.39 0.67 -0.33 0.18 -0.06 0.15 # [1] 6515.86 # Try dropping the first path: newShare = mxModel(share, mxMatrix("Lower", nrow=nVar, ncol=nVar, free=c(F,T,T), values=c(0,.7,.7), name="a", byrow=TRUE) ) # drop 2:2 at 0 model = mxModel("twinCE",newShare, mzModel, dzModel, mxAlgebra(MZ.objective + DZ.objective, name="twin"), mxAlgebraObjective("twin") ) fit = mxRun(model); summaryACEFit(fit, selVars); # The whole first column has been dropped # A1 A2 C1 C2 E1 E2 # wt1 0 NA 0.64 NA 0.36 NA # ht1 0 0.66 0.81 0.21 0.19 0.14 # [1] 6617.06 # Have a peak a the relevant matrix in the fitted model: it is fixed as desired, and there is a value for r2c1 # as there should be... perhaps this is my misunderstanding about to get from 'a' to 'A'? ... # fit@submodels$z@matrices$a # LowerMatrix 'a' # @labels # [,1] [,2] # [1,] "v1a1" NA # [2,] "v2a1" "v2a2" # @values # [,1] [,2] # [1,] 0.00000000 0.0000000 # [2,] -0.02042795 -0.0507354 # @free # [,1] [,2] # [1,] FALSE FALSE # [2,] TRUE TRUE