ifa: Store latent mean & cov in regular MxMatrix
[openmx:openmx.git] / models / failing / NTF_design.R
1 # SCRIPT: NTF_design.R - NTF design in OpenMx 
2 # Author: Matt Keller & Sarah Medland; edited by Tim Bates
3 # History:  Thu Sep 24 17:23:33 BST 2009
4 # For description of model, see Keller, Medland, Duncan, Hatemi, Neale, Maes, & Eaves (2009) TRHG, 21, p.8 - 18.
5 # 2009-09-26: (tb) read data from DEMO folder   
6 # OpenMx: http://www.openmx.virginia.com
7 # TODO (tb): Add closeenough calls
8 # TODO (tb): Resolve question about fixing m...
9 ##########################################
10
11 require(OpenMx);
12 require(OpenMx);
13 # Get Data
14 data(nuclear_twin_design_data)
15 selVars <- names(nuclear_twin_design_data)[1:4]
16 mzData  <- nuclear_twin_design_data[nuclear_twin_design_data$zyg=='mz',selVars]
17 dzData  <- nuclear_twin_design_data[nuclear_twin_design_data$zyg=='dz',selVars]
18
19 #Fit NTF Model with RawData and Matrices Input
20 ntf <- mxModel(model="NucTwFam",
21         # Matrices
22         # NOTE: NTF design does not allow Vs & Vf to be estimated simultaneously; for identifiability, choose either m or s to be free and the other to be fixed at 0             
23         # mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=0,  label="FamilialPath", name="m"),  # fix m=0 if you want Vf=0
24         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.7,  label="Env",          name="e"),
25         mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1,  label="FamilialVar",  name="f"),
26         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=0.1, label="FamilialPath", name="m"),  # fix m=0 if you want Vf=0  (nope... non pos def)
27         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.6,  label="AddGen",       name="a"),
28         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.3,  label="Sib",          name="s"),  # fix s=0 if you want Vs=0
29         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2,  label="Dominance",    name="d"),
30         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.1,  label="AMCopath",     name="mu"),
31         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=1.2, label="VarPhen",      name="Vp1"),
32         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2,  label="VarF",         name="x1"),  # keep this parameter free, even if Vf is fixed to 0
33         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.25, label="CovPhenGen",   name="delta1"),
34         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=1,   label="VarAddGen",    name="q1"),
35         mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.15, label="CovFA",        name="w1"),
36   #mxAlgebra section - nonlinear constraints
37   mxAlgebra(expression= a %*% q1 %*% t(a) + f %*% x1 %*% t(f) + 2 %x% a %*% w1 %*% t(f) + e %*% t(e) + d %*% t(d) + s %*% t(s), name="Vp2"),
38   mxAlgebra(expression= 2 %x% m %*% Vp1 %*% t(m) + 2 %x% m %*% Vp1 %*% mu %*% t(Vp1) %*% t(m), name="x2"),
39   mxAlgebra(expression= q1 %*% a + w1 %*% f, name="delta2"),
40   mxAlgebra(expression= 1 + delta1 %*% mu %*% t(delta1), name="q2"),
41   mxAlgebra(expression= delta1 %*% m + delta1 %*% mu %*% Vp1 %*% t(m), name="w2"),
42   #constraints - equating nonlinear constraints and parameters
43   mxConstraint('Vp1',"=",'Vp2',name='VpCon'),
44   mxConstraint('x1',"=",'x2',name='xCon'),
45   mxConstraint('delta1',"=",'delta2',name='deltaCon'),
46   mxConstraint('q1',"=",'q2',name='qCon'),
47   mxConstraint('w1',"=",'w2',name='wCon'),
48   #mxAlgebra section - relative covariances
49   mxAlgebra(expression= a %*% q1 %*% t(a) + f %*% x1 %*% t(f) + 2 %x% a %*% w1 %*% t(f) + d %*% t(d) + s %*% t(s), name="CvMz"),
50   mxAlgebra(expression= a %*% (q1-.5) %*% t(a) + .25 %x% d %*% t(d) + f %*% x1 %*% t(f) + 2 %x% a %*% w1 %*% t(f) + s %*% t(s), name="CvDz"),
51   mxAlgebra(expression= .5 %x% a %*% (q1 %*% a + w1 %*% f) + .5 %x% a %*% (q1 %*% a + w1 %*% f) %*% mu %*% Vp1 + m %*% Vp1 + m %*% Vp1 %*% mu %*% t(Vp1), name="ParChild"),
52   mxAlgebra(expression= Vp1 %*% mu %*% t(Vp1), name="CvSps")
53 )
54
55 mzModel <- mxModel(ntf, name = "MZNTF",
56   mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(mzData)), name="expMeanMz"),
57   # Algebra for expected variance/covariance matrix in MZF
58   mxAlgebra(expression=rbind(
59           cbind(Vp1,       CvMz,       ParChild, ParChild),
60           cbind(CvMz,      Vp1,        ParChild, ParChild),
61           cbind(ParChild,  ParChild,   Vp1,      CvSps),
62           cbind(ParChild,  ParChild,   CvSps,     Vp1)
63   ), dimnames=list(colnames(mzData),colnames(mzData)),name="expCovMz"),
64   mxData(observed=mzData, type="raw"),
65   mxFIMLObjective(covariance="expCovMz",means="expMeanMz")
66 )
67
68 dzModel <- mxModel(ntf, name = "DZNTF",
69   mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(dzData)), name="expMeanDz"),
70   # Algebra for expected variance/covariance matrix in MZF
71   mxAlgebra(expression=rbind(
72           cbind(Vp1,      CvDz,     ParChild, ParChild),
73           cbind(CvDz,     Vp1,      ParChild, ParChild),
74           cbind(ParChild, ParChild, Vp1,      CvSps),
75           cbind(ParChild, ParChild, CvSps,    Vp1)
76   ), dimnames=list(colnames(dzData),colnames(dzData)),name="expCovDz"),
77   mxData(observed=dzData, type="raw"),
78   mxFIMLObjective(covariance="expCovDz",means="expMeanDz")
79 )
80
81 model <- mxModel(model="NTF", mzModel, dzModel,
82   mxAlgebra(expression=MZNTF.objective + DZNTF.objective, name="ntffit"), #MZ.objective is the automatic name for the -2LL of mzModel
83   mxFitFunctionAlgebra("ntffit")
84 )
85
86 #Run MX
87 start <- proc.time()[1]; # record start time
88 fit <- mxRun(model)
89 endTime <- proc.time()[1]; (endTime -start)
90
91 #Look at results
92 estimate <- fit@output$estimate
93 round(estimate,3)
94 # TODO include simulation output
95 # compare to simulation
96 # estimate.mat <- rbind(round(c(estimate[1:5]^2,res[6:7]),3),round(ALL$track.changes[c('var.U','var.F','var.A','var.S','var.D','cor.spouses','var.cur.phenotype'),'data.t1'],3))
97 # dimnames(estimate.mat) <- list(c('OpenMx-Estimated','Simulated'),c('Var.E','Var.F','Var.A','Var.S','Var.D','Cor.Sps','Var.Phen'))
98 # estimate == estimate.mat
99
100
101 #NOTE: the minor difference between simulated & estimated parameters are to be expected & are due to sampling error
102 # > res.mat
103 #                  Var.E   Var.F  Var.A  Var.S   Var.D  Cor.Sps Var.Phen
104 #OpenMx-Estimated 0.443   0.100   0.294   0.149     0   0.209     0.986
105 #Old.mx           0.442   0.099   0.294   0.149     0   0.206     0.986
106 #Simulated        0.433   0.100   0.304   0.197     0   0.200     1.032