Update demos directory to reflect new style header and new commenting style.
[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: 2009.08.01 
20 #
21 # ModelType: ACE
22 # DataType: Simulated Twin
23 # Field: Human Behavior Genetics
24 #
25 # Purpose: 
26 #      Univariate Twin Analysis Example in OpenMx:
27 #      From fitting saturated models to testing model assumptions 
28 #      To fitting the ACE model and a submodel
29 #
30 # RevisionHistory:
31 #      Hermine Maes -- 10 08 2009 updated & reformatted
32 #      Ross Gore -- 06 06 2011  added Model, Data & Field metadata
33 # -----------------------------------------------------------------------
34
35 require(OpenMx)
36 require(MASS)
37 # Load Library
38 # -----------------------------------------------------------------------
39
40 set.seed(200)
41 a2<-0.5         #Additive genetic variance component (a squared)
42 c2<-0.3         #Common environment variance component (c squared)
43 e2<-0.2         #Specific environment variance component (e squared)
44 rMZ <- a2+c2
45 rDZ <- .5*a2+c2
46 DataMZ <- mvrnorm (1000, c(0,0), matrix(c(1,rMZ,rMZ,1),2,2))
47 DataDZ <- mvrnorm (1000, c(0,0), matrix(c(1,rDZ,rDZ,1),2,2))
48
49 selVars <- c('t1','t2')
50 dimnames(DataMZ) <- list(NULL,selVars)
51 dimnames(DataDZ) <- list(NULL,selVars)
52 summary(DataMZ)
53 summary(DataDZ)
54 colMeans(DataMZ,na.rm=TRUE)
55 colMeans(DataDZ,na.rm=TRUE)
56 cov(DataMZ,use="complete")
57 cov(DataDZ,use="complete")
58 # Simulate Data: two standardized variables t1 & t2 for MZ's & DZ's
59 # -----------------------------------------------------------------------
60
61 twinSatModel <- mxModel("twinSat",
62         mxModel("MZ",
63                 mxMatrix("Full", 1, 2, T, c(0,0), name="expMeanMZ"), 
64                 mxMatrix("Lower", 2, 2, T, .5, name="CholMZ"), 
65                 mxAlgebra(CholMZ %*% t(CholMZ), name="expCovMZ"), 
66                 mxData(DataMZ, type="raw"),
67                 mxFIMLObjective("expCovMZ", "expMeanMZ", selVars)),  
68         mxModel("DZ",
69                 mxMatrix("Full", 1, 2, T, c(0,0), name="expMeanDZ"), 
70                 mxMatrix("Lower", 2, 2, T, .5, name="CholDZ"), 
71                 mxAlgebra(CholDZ %*% t(CholDZ), name="expCovDZ"), 
72                 mxData(DataDZ, type="raw"), 
73                 mxFIMLObjective("expCovDZ", "expMeanDZ", selVars)),
74         mxAlgebra(MZ.objective + DZ.objective, name="twin"), 
75         mxAlgebraObjective("twin"))
76 twinSatFit <- mxRun(twinSatModel, suppressWarnings=TRUE)
77 # Specify and Run Saturated Model with RawData and Matrix-style Input
78 # -----------------------------------------------------------------------
79
80
81 ExpMeanMZ <- mxEval(MZ.expMeanMZ, twinSatFit)
82 ExpCovMZ <- mxEval(MZ.expCovMZ, twinSatFit)
83 ExpMeanDZ <- mxEval(DZ.expMeanDZ, twinSatFit)
84 ExpCovDZ <- mxEval(DZ.expCovDZ, twinSatFit)
85 LL_Sat <- mxEval(objective, twinSatFit)
86 # Generate Saturated Model Output
87 # -----------------------------------------------------------------------
88
89 twinSatModelSub1 <- twinSatModel
90 twinSatModelSub1$MZ$expMeanMZ <- mxMatrix("Full", 1, 2, T, 0, "mMZ", name="expMeanMZ")
91 twinSatModelSub1$DZ$expMeanDZ <- mxMatrix("Full", 1, 2, T, 0, "mDZ", name="expMeanDZ")
92 twinSatFitSub1 <- mxRun(twinSatModelSub1, suppressWarnings=TRUE)
93 # Specify and Run Saturated SubModel 1 equating means across twin order
94 # -----------------------------------------------------------------------
95
96 twinSatModelSub2 <- twinSatModelSub1
97 twinSatModelSub2$MZ$expMeanMZ <- mxMatrix("Full", 1, 2, T, 0, "mean", name="expMeanMZ")
98 twinSatModelSub2$DZ$expMeanDZ <- mxMatrix("Full", 1, 2, T, 0, "mean", name="expMeanDZ")
99 twinSatFitSub2 <- mxRun(twinSatModelSub2, suppressWarnings=TRUE)
100 # Specify and Run Saturated SubModel 2 equating means across zygosity
101 # -----------------------------------------------------------------------
102
103 LL_Sat <- mxEval(objective, twinSatFit)
104 LL_Sub1 <- mxEval(objective, twinSatFitSub1)
105 LRT1 <- LL_Sub1 - LL_Sat
106 LL_Sub2 <- mxEval(objective, twinSatFitSub2)
107 LRT2 <- LL_Sub2 - LL_Sat
108 # Generate Saturated Model Comparison Output
109 # -----------------------------------------------------------------------
110
111 twinACEModel <- mxModel("twinACE", 
112         mxMatrix("Full", 1, 2, T, 20, "mean", name="expMean"), 
113         # Matrix expMean for expected mean 
114         # vector for MZ and DZ twins    
115         # -------------------------------------
116         
117         mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="a", name="X"),
118         mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="c", name="Y"),
119         mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="e", name="Z"),
120         # Matrices X, Y, and Z to store the 
121         # a, c, and e path coefficients
122         # -------------------------------------
123         
124         mxAlgebra(X * t(X), name="A"),
125         mxAlgebra(Y * t(Y), name="C"),
126         mxAlgebra(Z * t(Z), name="E"),  
127         # Matrixes A, C, and E to compute 
128         # A, C, and E variance components
129         # -------------------------------------
130         
131         mxAlgebra(rbind(cbind(A+C+E   , A+C),
132                                         cbind(A+C     , A+C+E)), name="expCovMZ"),
133         # Matrix expCOVMZ for expected 
134         # covariance matrix for MZ twins
135         # -------------------------------------
136         
137         mxModel("MZ",
138                 mxData(DataMZ, type="raw"), 
139                 mxFIMLObjective("twinACE.expCovMZ", "twinACE.expMean",selVars)),
140
141         mxAlgebra(rbind(cbind(A+C+E   , .5%x%A+C),
142                                         cbind(.5%x%A+C , A+C+E)), name="expCovDZ"),
143         # Matrix expCOVMZ for expected 
144         # covariance matrix for DZ twins
145         # -------------------------------------
146         
147         mxModel("DZ", 
148                 mxData(DataDZ, type="raw"), 
149                 mxFIMLObjective("twinACE.expCovDZ", "twinACE.expMean",selVars)),
150
151         mxAlgebra(MZ.objective + DZ.objective, name="twin"), 
152         mxAlgebraObjective("twin"))
153 twinACEFit <- mxRun(twinACEModel, suppressWarnings=TRUE)
154 # Specify and Run ACE Model with RawData and Matrix-style Input
155 # -----------------------------------------------------------------------
156
157 LL_ACE <- mxEval(objective, twinACEFit)
158 LRT_ACE= LL_ACE - LL_Sat
159         MZc <- mxEval(expCovMZ, twinACEFit)
160         DZc <- mxEval(expCovDZ, twinACEFit)
161         M   <- mxEval(expMean, twinACEFit)
162         # Retrieve expected mean vector and 
163         # expected covariance matrices
164         # -------------------------------------
165         
166         A <- mxEval(A, twinACEFit)
167         C <- mxEval(C, twinACEFit)
168         E <- mxEval(E, twinACEFit)
169         # Retrieve the A, C, and E 
170         # variance components
171         # -------------------------------------
172         
173         V <- (A+C+E)
174         a2 <- A/V
175         c2 <- C/V
176         e2 <- E/V
177         # Calculate standardized variance 
178         # components
179         # -------------------------------------
180         
181         ACEest <- rbind(cbind(A,C,E),cbind(a2,c2,e2)) 
182         ACEest <- data.frame(ACEest, row.names=c("Variance Components","Standardized VC"))
183         names(ACEest)<-c("A", "C", "E")
184         ACEest; LL_ACE; LRT_ACE
185         #Build and print reporting table with 
186         # row and column names
187         # -------------------------------------
188         
189 # Generate ACE Model Output
190 # -----------------------------------------------------------------------
191         
192 twinAEModel <- twinACEModel
193 twinAEModel$twinACE$Y <- mxMatrix("Full", 1, 1, F, 0, "c", name="Y")  # drop c
194 twinAEFit <- mxRun(twinAEModel, suppressWarnings=TRUE)
195 # Specify and Fit Reduced AE Model (Drop c @0)
196 # ----------------------------------------------------------------------
197
198 LL_AE <- mxEval(objective, twinAEFit)
199
200         MZc <- mxEval(expCovMZ, twinAEFit)
201         DZc <- mxEval(expCovDZ, twinAEFit)
202         M   <- mxEval(expMean, twinAEFit)
203         # Retrieve expected mean vector and 
204         # expected covariance matrices
205         # -------------------------------------
206         A <- mxEval(A, twinAEFit)
207         C <- mxEval(C, twinAEFit)
208         E <- mxEval(E, twinAEFit)
209         # Retrieve the A, C and E variance 
210         # components
211         # -------------------------------------
212         
213         V <- (A+C+E)
214         a2 <- A/V
215         c2 <- C/V
216         e2 <- E/V
217         # Calculate standardized variance 
218         # components
219         # -------------------------------------
220
221         AEest <- rbind(cbind(A,C,E),cbind(a2,c2,e2)) 
222         AEest <- data.frame(ACEest, row.names=c("Variance Components","Standardized VC"))
223         names(ACEest)<-c("A", "C", "E")
224         AEest; LL_AE; 
225         # Build and print reporting table with 
226         # row and column names
227         # -------------------------------------
228
229         LRT_ACE_AE <- LL_AE - LL_ACE
230         LRT_ACE_AE
231         # Calculate and print likelihood ratio 
232         # test
233         # -------------------------------------
234 # Generate ACE Model Output
235 # -----------------------------------------------------------------------