copied scripts from passing
[openmx:openmx.git] / demo / UnivariateTwinAnalysis_Matrix.R
1 require(OpenMx)
2
3 #Prepare Data
4 twinData <- read.table("myTwinData.txt", header=T, na.strings=".")
5 twinVars <- c('fam','age','zyg','part','wt1','wt2','ht1','ht2','htwt1','htwt2','bmi1','bmi2')
6 #dimnames(twinData) <- list(NULL, twinVars)
7 summary(twinData)
8 selVars <- c('bmi1','bmi2')
9 mzfData <- as.matrix(subset(twinData, zyg==1, c(bmi1,bmi2)))
10 dzfData <- as.matrix(subset(twinData, zyg==3, c(bmi1,bmi2)))
11 cov(mzfData)
12 cov(dzfData)
13
14 #Fit ACE Model with RawData and Matrices Input
15 twinACEModel <- mxModel("twinACE", 
16     mxMatrix("Full", 1, 2, T, c(20,20), labels= c("mean","mean"), dimnames=list(NULL, selVars), name="expMeanMZ"), 
17     mxMatrix("Full", 1, 2, T, c(20,20), labels= c("mean","mean"), dimnames=list(NULL, selVars), name="expMeanDZ"), 
18     mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="a", name="X"),
19     mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="c", name="Y"),
20     mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="e", name="Z"),
21     mxMatrix("Full", nrow=1, ncol=1, free=FALSE, values=.5, name="h"),
22     mxAlgebra(X * t(X), name="A"),
23     mxAlgebra(Y * t(Y), name="C"),
24     mxAlgebra(Z * t(Z), name="E"), 
25     mxAlgebra(rbind (cbind(A + C + E, A + C),
26                      cbind(A + C    , A + C + E)), dimnames = list(selVars, selVars), name="expCovMZ"),
27     mxAlgebra(rbind (cbind(A + C + E  , h %x% A + C),
28                      cbind(h %x% A + C, A + C + E)), dimnames = list(selVars, selVars), name="expCovDZ"),
29     mxModel("MZ",
30         mxData(mzfData, type="raw"), 
31         mxFIMLObjective("twinACE.expCovMZ", "twinACE.expMeanMZ")),
32     mxModel("DZ", 
33         mxData(dzfData, type="raw"), 
34         mxFIMLObjective("twinACE.expCovDZ", "twinACE.expMeanDZ")),
35     mxAlgebra(MZ.objective + DZ.objective, name="twin"), 
36     mxAlgebraObjective("twin"))
37
38 #Run ACE model
39 twinACEFit <- mxRun(twinACEModel)
40
41 MZc <- mxEvaluate(expCovMZ, twinACEFit)
42 DZc <- mxEvaluate(expCovDZ, twinACEFit)
43 M <- mxEvaluate(expMeanMZ, twinACEFit)
44 A <- mxEvaluate(A, twinACEFit)
45 C <- mxEvaluate(C, twinACEFit)
46 E <- mxEvaluate(E, twinACEFit)
47 V <- (A+C+E)
48 a2 <- A/V
49 c2 <- C/V
50 e2 <- E/V
51 ACEest <- rbind(cbind(A,C,E),cbind(a2,c2,e2))
52 LL_ACE <- mxEvaluate(objective, twinACEFit)
53
54
55 #Run Mx scripts
56 #mymatrices1 <- omxOriginalMx("mx-scripts/univACE.mx", "temp-files")                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
57 #Mx.A <- mymatrices1$A1.1
58 #Mx.C <- mymatrices1$C1.1
59 #Mx.E <- mymatrices1$E1.1
60 #Mx.M <- mymatrices1$M2.1
61 #Mx.LL_ACE <- mymatrices1$F4.1;
62
63 #Mx answers hard-coded
64 #1: Heterogeneity Model
65 Mx.A <- 0.6173023
66 Mx.C <- 5.595822e-14
67 Mx.E <- 0.1730462
68 Mx.M <- matrix(c(21.39293, 21.39293),1,2)
69 Mx.LL_ACE <- 4067.663
70
71
72 #Compare OpenMx results to Mx results (LL: likelihood; EC: expected covariance, EM: expected means)
73 omxCheckCloseEnough(LL_ACE,Mx.LL_ACE,.001)
74 omxCheckCloseEnough(A,Mx.A,.001)
75 omxCheckCloseEnough(C,Mx.C,.001)
76 omxCheckCloseEnough(E,Mx.E,.001)
77 omxCheckCloseEnough(M,Mx.M,.001)
78
79
80 #Run AE model
81 twinAEModel <- mxModel(twinACEModel, 
82     mxMatrix("Full", nrow=1, ncol=1, free=F, values=0, label="c", name="Y")
83     )
84 twinAEFit <- mxRun(twinAEModel)
85
86 MZc <- mxEvaluate(expCovMZ, twinAEFit)
87 DZc <- mxEvaluate(expCovDZ, twinAEFit)
88 A <- mxEvaluate(A, twinAEFit)
89 C <- mxEvaluate(C, twinAEFit)
90 E <- mxEvaluate(E, twinAEFit)
91 V <- (A + C + E)
92 a2 <- A / V
93 c2 <- C / V
94 e2 <- E / V
95 AEest <- rbind(cbind(A, C, E),cbind(a2, c2, e2))
96 LL_AE <- mxEvaluate(objective, twinAEFit)
97
98 LRT_ACE_AE <- LL_AE - LL_ACE
99
100 #Print relevant output
101 ACEest
102 AEest
103 LRT_ACE_AE