copied scripts from passing
[openmx:openmx.git] / demo / UnivariateTwinAnalysis_Path.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 aceVars <- c("A1","C1","E1","A2","C2","E2")
10 mzfData <- as.matrix(subset(twinData, zyg==1, c(bmi1,bmi2)))
11 dzfData <- as.matrix(subset(twinData, zyg==3, c(bmi1,bmi2)))
12 cov(mzfData)
13 cov(dzfData)
14
15 #Fit ACE Model with RawData and Path-Style Input
16 share <- mxModel("share", 
17     type="RAM",
18     manifestVars=selVars,
19     latentVars=aceVars,
20     mxPath(from=aceVars, arrows=2, free=FALSE, values=1),
21     mxPath(from="one", to=aceVars, arrows=1, free=FALSE, values=0),
22     mxPath(from="one", to=selVars, arrows=1, free=TRUE, values=20, labels= c("mean","mean")),
23     mxPath(from=c("A1","C1","E1"), to="bmi1", arrows=1, free=TRUE, values=.6, label=c("a","c","e")),
24     mxPath(from=c("A2","C2","E2"), to="bmi2", arrows=1, free=TRUE, values=.6, label=c("a","c","e")),
25     mxPath(from="C1", to="C2", arrows=2, free=FALSE, values=1)
26     )
27 twinACEModel <- mxModel("twinACE", 
28     mxModel(share,
29         mxPath(from="A1", to="A2", arrows=2, free=FALSE, values=1),
30         mxData(mzfData, type="raw"), 
31         type="RAM", name="MZ"),
32     mxModel(share, 
33         mxPath(from="A1", to="A2", arrows=2, free=FALSE, values=.5),
34         mxData(dzfData, type="raw"), 
35         type="RAM", name="DZ"),
36     mxAlgebra(MZ.objective + DZ.objective, name="twin"), 
37     mxAlgebraObjective("twin"))
38
39 #Run ACE model
40 twinACEFit <- mxRun(twinACEModel)
41
42 MZc <- mxEvaluate(MZ.covariance, twinACEFit)
43 DZc <- mxEvaluate(DZ.covariance, twinACEFit)
44 M <- mxEvaluate(MZ.means, twinACEFit)
45 A <- mxEvaluate(a*a, twinACEFit)
46 C <- mxEvaluate(c*c, twinACEFit)
47 E <- mxEvaluate(e*e, twinACEFit)
48 V <- (A+C+E)
49 a2 <- A/V
50 c2 <- C/V
51 e2 <- E/V
52 ACEest <- rbind(cbind(A,C,E),cbind(a2,c2,e2))
53 LL_ACE <- mxEvaluate(objective, twinACEFit)
54
55
56 #Mx answers hard-coded
57 #1: Heterogeneity Model
58 Mx.A <- 0.6173023
59 Mx.C <- 5.595822e-14
60 Mx.E <- 0.1730462
61 Mx.M <- matrix(c(21.39293, 21.39293),1,2)
62 Mx.LL_ACE <- 4067.663
63
64
65 #Compare OpenMx results to Mx results (LL: likelihood; EC: expected covariance, EM: expected means)
66 omxCheckCloseEnough(LL_ACE,Mx.LL_ACE,.001)
67 omxCheckCloseEnough(A,Mx.A,.001)
68 omxCheckCloseEnough(C,Mx.C,.001)
69 omxCheckCloseEnough(E,Mx.E,.001)
70 omxCheckCloseEnough(M,t(Mx.M),.001)
71
72
73 #Run AE model
74 twinAEModel <- mxModel(twinACEModel, type="RAM",
75     manifestVars=selVars,
76     latentVars=aceVars,
77     mxPath(from=c("A1","C1","E1"), to="bmi1", arrows=1, free=c(T,F,T), values=c(.6,0,.6), label=c("a","cfixed","e")),
78     mxPath(from=c("A2","C2","E2"), to="bmi2", arrows=1, free=c(T,F,T), values=c(.6,0,.6), label=c("a","cfixed","e"))
79     )
80 twinAEFit <- mxRun(twinAEModel)
81
82 MZc <- mxEvaluate(MZ.covariance, twinAEFit)
83 DZc <- mxEvaluate(DZ.covariance, twinAEFit)
84 M <- mxEvaluate(MZ.means, twinAEFit)
85 A <- mxEvaluate(a*a, twinAEFit)
86 C <- mxEvaluate(c*c, twinAEFit)
87 E <- mxEvaluate(e*e, twinAEFit)
88 V <- (A + C + E)
89 a2 <- A / V
90 c2 <- C / V
91 e2 <- E / V
92 AEest <- rbind(cbind(A, C, E),cbind(a2, c2, e2))
93 LL_AE <- mxEvaluate(objective, twinAEFit)
94
95 LRT_ACE_AE <- LL_AE - LL_ACE
96
97 #Print relevant output
98 ACEest
99 AEest
100 LRT_ACE_AE