new 'suppressWarnings=TRUE' argument to mxRun()
[openmx:openmx.git] / demo / UnivariateTwinAnalysis.R
1 #
2 #   Copyright 2007-2010 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #   Unless required by applicable law or agreed to in writing, software
11 #   distributed under the License is distributed on an "AS IS" BASIS,
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16 # -----------------------------------------------------------------------
17 # Program: UnivariateTwinAnalysis.R  
18 #  Author: Hermine Maes
19 #    Date: 08 01 2009 
20 #
21 # Univariate Twin Analysis Example in OpenMx:
22 # From fitting saturated models to testing model assumptions 
23 # To fitting the ACE model and a submodel
24 #
25 # Revision History
26 #   Hermine Maes -- 10 08 2009 updated & reformatted
27 # -----------------------------------------------------------------------
28
29 require(OpenMx)
30
31 # Simulate Data: two standardized variables t1 & t2 for MZ's & DZ's
32 # -----------------------------------------------------------------------
33 require(MASS)
34 set.seed(200)
35 a2<-0.5         #Additive genetic variance component (a squared)
36 c2<-0.3         #Common environment variance component (c squared)
37 e2<-0.2         #Specific environment variance component (e squared)
38 rMZ <- a2+c2
39 rDZ <- .5*a2+c2
40 DataMZ <- mvrnorm (1000, c(0,0), matrix(c(1,rMZ,rMZ,1),2,2))
41 DataDZ <- mvrnorm (1000, c(0,0), matrix(c(1,rDZ,rDZ,1),2,2))
42
43 selVars <- c('t1','t2')
44 dimnames(DataMZ) <- list(NULL,selVars)
45 dimnames(DataDZ) <- list(NULL,selVars)
46 summary(DataMZ)
47 summary(DataDZ)
48 colMeans(DataMZ,na.rm=TRUE)
49 colMeans(DataDZ,na.rm=TRUE)
50 cov(DataMZ,use="complete")
51 cov(DataDZ,use="complete")
52
53 # Specify and Run Saturated Model with RawData and Matrix-style Input
54 # -----------------------------------------------------------------------
55 twinSatModel <- mxModel("twinSat",
56         mxModel("MZ",
57                 mxMatrix("Full", 1, 2, T, c(0,0), name="expMeanMZ"), 
58                 mxMatrix("Lower", 2, 2, T, .5, name="CholMZ"), 
59                 mxAlgebra(CholMZ %*% t(CholMZ), name="expCovMZ"), 
60                 mxData(DataMZ, type="raw"),
61                 mxFIMLObjective("expCovMZ", "expMeanMZ", selVars)),  
62         mxModel("DZ",
63                 mxMatrix("Full", 1, 2, T, c(0,0), name="expMeanDZ"), 
64                 mxMatrix("Lower", 2, 2, T, .5, name="CholDZ"), 
65                 mxAlgebra(CholDZ %*% t(CholDZ), name="expCovDZ"), 
66                 mxData(DataDZ, type="raw"), 
67                 mxFIMLObjective("expCovDZ", "expMeanDZ", selVars)),
68         mxAlgebra(MZ.objective + DZ.objective, name="twin"), 
69         mxAlgebraObjective("twin"))
70 twinSatFit <- mxRun(twinSatModel, suppressWarnings=TRUE)
71
72 # Generate Saturated Model Output
73 # -----------------------------------------------------------------------
74 ExpMeanMZ <- mxEval(MZ.expMeanMZ, twinSatFit)
75 ExpCovMZ <- mxEval(MZ.expCovMZ, twinSatFit)
76 ExpMeanDZ <- mxEval(DZ.expMeanDZ, twinSatFit)
77 ExpCovDZ <- mxEval(DZ.expCovDZ, twinSatFit)
78 LL_Sat <- mxEval(objective, twinSatFit)
79
80
81 # Specify and Run Saturated SubModel 1 equating means across twin order
82 # -----------------------------------------------------------------------
83 twinSatModelSub1 <- twinSatModel
84 twinSatModelSub1$MZ$expMeanMZ <- mxMatrix("Full", 1, 2, T, 0, "mMZ", name="expMeanMZ")
85 twinSatModelSub1$DZ$expMeanDZ <- mxMatrix("Full", 1, 2, T, 0, "mDZ", name="expMeanDZ")
86 twinSatFitSub1 <- mxRun(twinSatModelSub1, suppressWarnings=TRUE)
87
88 # Specify and Run Saturated SubModel 2 equating means across zygosity
89 # -----------------------------------------------------------------------
90 twinSatModelSub2 <- twinSatModelSub1
91 twinSatModelSub2$MZ$expMeanMZ <- mxMatrix("Full", 1, 2, T, 0, "mean", name="expMeanMZ")
92 twinSatModelSub2$DZ$expMeanDZ <- mxMatrix("Full", 1, 2, T, 0, "mean", name="expMeanDZ")
93 twinSatFitSub2 <- mxRun(twinSatModelSub2, suppressWarnings=TRUE)
94
95 # Generate Saturated Model Comparison Output
96 # -----------------------------------------------------------------------
97 LL_Sat <- mxEval(objective, twinSatFit)
98 LL_Sub1 <- mxEval(objective, twinSatFitSub1)
99 LRT1 <- LL_Sub1 - LL_Sat
100 LL_Sub2 <- mxEval(objective, twinSatFitSub2)
101 LRT2 <- LL_Sub2 - LL_Sat
102
103
104 # Specify and Run ACE Model with RawData and Matrix-style Input
105 # -----------------------------------------------------------------------
106 twinACEModel <- mxModel("twinACE", 
107         mxMatrix("Full", 1, 2, T, 20, "mean", name="expMean"), 
108                 # Matrix expMean for expected mean vector for MZ and DZ twins    
109
110         mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="a", name="X"),
111         mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="c", name="Y"),
112         mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="e", name="Z"),
113                 # Matrices X, Y, and Z to store the a, c, and e path coefficients
114
115         mxAlgebra(X * t(X), name="A"),
116         mxAlgebra(Y * t(Y), name="C"),
117         mxAlgebra(Z * t(Z), name="E"),  
118                 # Matrixes A, C, and E to compute A, C, and E variance components
119
120         mxAlgebra(rbind(cbind(A+C+E   , A+C),
121                                         cbind(A+C     , A+C+E)), name="expCovMZ"),
122                 # Matrix expCOVMZ for expected covariance matrix for MZ twins
123         mxModel("MZ",
124                 mxData(DataMZ, type="raw"), 
125                 mxFIMLObjective("twinACE.expCovMZ", "twinACE.expMean",selVars)),
126
127         mxAlgebra(rbind(cbind(A+C+E   , .5%x%A+C),
128                                         cbind(.5%x%A+C , A+C+E)), name="expCovDZ"),
129                 # Matrix expCOVMZ for expected covariance matrix for DZ twins
130         mxModel("DZ", 
131                 mxData(DataDZ, type="raw"), 
132                 mxFIMLObjective("twinACE.expCovDZ", "twinACE.expMean",selVars)),
133
134         mxAlgebra(MZ.objective + DZ.objective, name="twin"), 
135         mxAlgebraObjective("twin"))
136 twinACEFit <- mxRun(twinACEModel, suppressWarnings=TRUE)
137
138 # Generate ACE Model Output
139 # -----------------------------------------------------------------------
140 LL_ACE <- mxEval(objective, twinACEFit)
141 LRT_ACE= LL_ACE - LL_Sat
142 #Retrieve expected mean vector and expected covariance matrices
143         MZc <- mxEval(expCovMZ, twinACEFit)
144         DZc <- mxEval(expCovDZ, twinACEFit)
145         M   <- mxEval(expMean, twinACEFit)
146 #Retrieve the A, C, and E variance components
147         A <- mxEval(A, twinACEFit)
148         C <- mxEval(C, twinACEFit)
149         E <- mxEval(E, twinACEFit)
150 #Calculate standardized variance components
151         V <- (A+C+E)
152         a2 <- A/V
153         c2 <- C/V
154         e2 <- E/V
155 #Build and print reporting table with row and column names
156         ACEest <- rbind(cbind(A,C,E),cbind(a2,c2,e2)) 
157         ACEest <- data.frame(ACEest, row.names=c("Variance Components","Standardized VC"))
158         names(ACEest)<-c("A", "C", "E")
159         ACEest; LL_ACE; LRT_ACE
160
161 # Specify and Fit Reduced AE Model (Drop c @0)
162 # ----------------------------------------------------------------------
163 twinAEModel <- twinACEModel
164 twinAEModel$twinACE$Y <- mxMatrix("Full", 1, 1, F, 0, "c", name="Y")  # drop c
165 twinAEFit <- mxRun(twinAEModel, suppressWarnings=TRUE)
166
167 # Generate ACE Model Output
168 # -----------------------------------------------------------------------
169 LL_AE <- mxEval(objective, twinAEFit)
170 #Retrieve expected mean vector and expected covariance matrices
171         MZc <- mxEval(expCovMZ, twinAEFit)
172         DZc <- mxEval(expCovDZ, twinAEFit)
173         M   <- mxEval(expMean, twinAEFit)
174 #Retrieve the A, C and E variance components
175         A <- mxEval(A, twinAEFit)
176         C <- mxEval(C, twinAEFit)
177         E <- mxEval(E, twinAEFit)
178 #Calculate standardized variance components
179         V <- (A+C+E)
180         a2 <- A/V
181         c2 <- C/V
182         e2 <- E/V
183 #Build and print reporting table with row and column names
184         AEest <- rbind(cbind(A,C,E),cbind(a2,c2,e2)) 
185         AEest <- data.frame(ACEest, row.names=c("Variance Components","Standardized VC"))
186         names(ACEest)<-c("A", "C", "E")
187         AEest; LL_AE; 
188 #Calculate and print likelihood ratio test
189         LRT_ACE_AE <- LL_AE - LL_ACE
190         LRT_ACE_AE