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