# SCRIPT: NTF_design.R - NTF design in OpenMx
# Author: Matt Keller & Sarah Medland; edited by Tim Bates
# History: Thu Sep 24 17:23:33 BST 2009
# For description of model, see Keller, Medland, Duncan, Hatemi, Neale, Maes, & Eaves (2009) TRHG, 21, p.8 - 18.
# 2009-09-26: (tb) read data from DEMO folder
# OpenMx: http://www.openmx.virginia.com
# TODO (tb): Add closeenough calls
# TODO (tb): Resolve question about fixing m...
##########################################
require(OpenMx);
require(OpenMx);
# Get Data
data(nuclear_twin_design_data)
selVars <- names(nuclear_twin_design_data)[1:4]
mzData <- nuclear_twin_design_data[nuclear_twin_design_data$zyg=='mz',selVars]
dzData <- nuclear_twin_design_data[nuclear_twin_design_data$zyg=='dz',selVars]
#Fit NTF Model with RawData and Matrices Input
ntf <- mxModel(model="NucTwFam",
# Matrices
# 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
# mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="FamilialPath", name="m"), # fix m=0 if you want Vf=0
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.7, label="Env", name="e"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="FamilialVar", name="f"),
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)
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="AddGen", name="a"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.3, label="Sib", name="s"), # fix s=0 if you want Vs=0
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2, label="Dominance", name="d"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.1, label="AMCopath", name="mu"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=1.2, label="VarPhen", name="Vp1"),
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
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.25, label="CovPhenGen", name="delta1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=1, label="VarAddGen", name="q1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.15, label="CovFA", name="w1"),
#mxAlgebra section - nonlinear constraints
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"),
mxAlgebra(expression= 2 %x% m %*% Vp1 %*% t(m) + 2 %x% m %*% Vp1 %*% mu %*% t(Vp1) %*% t(m), name="x2"),
mxAlgebra(expression= q1 %*% a + w1 %*% f, name="delta2"),
mxAlgebra(expression= 1 + delta1 %*% mu %*% t(delta1), name="q2"),
mxAlgebra(expression= delta1 %*% m + delta1 %*% mu %*% Vp1 %*% t(m), name="w2"),
#constraints - equating nonlinear constraints and parameters
mxConstraint('Vp1',"=",'Vp2',name='VpCon'),
mxConstraint('x1',"=",'x2',name='xCon'),
mxConstraint('delta1',"=",'delta2',name='deltaCon'),
mxConstraint('q1',"=",'q2',name='qCon'),
mxConstraint('w1',"=",'w2',name='wCon'),
#mxAlgebra section - relative covariances
mxAlgebra(expression= a %*% q1 %*% t(a) + f %*% x1 %*% t(f) + 2 %x% a %*% w1 %*% t(f) + d %*% t(d) + s %*% t(s), name="CvMz"),
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"),
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"),
mxAlgebra(expression= Vp1 %*% mu %*% t(Vp1), name="CvSps")
)
mzModel <- mxModel(ntf, name = "MZNTF",
mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(mzData)), name="expMeanMz"),
# Algebra for expected variance/covariance matrix in MZF
mxAlgebra(expression=rbind(
cbind(Vp1, CvMz, ParChild, ParChild),
cbind(CvMz, Vp1, ParChild, ParChild),
cbind(ParChild, ParChild, Vp1, CvSps),
cbind(ParChild, ParChild, CvSps, Vp1)
), dimnames=list(colnames(mzData),colnames(mzData)),name="expCovMz"),
mxData(observed=mzData, type="raw"),
mxFIMLObjective(covariance="expCovMz",means="expMeanMz")
)
dzModel <- mxModel(ntf, name = "DZNTF",
mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(dzData)), name="expMeanDz"),
# Algebra for expected variance/covariance matrix in MZF
mxAlgebra(expression=rbind(
cbind(Vp1, CvDz, ParChild, ParChild),
cbind(CvDz, Vp1, ParChild, ParChild),
cbind(ParChild, ParChild, Vp1, CvSps),
cbind(ParChild, ParChild, CvSps, Vp1)
), dimnames=list(colnames(dzData),colnames(dzData)),name="expCovDz"),
mxData(observed=dzData, type="raw"),
mxFIMLObjective(covariance="expCovDz",means="expMeanDz")
)
model <- mxModel(model="NTF", mzModel, dzModel,
mxAlgebra(expression=MZNTF.objective + DZNTF.objective, name="ntffit"), #MZ.objective is the automatic name for the -2LL of mzModel
mxFitFunctionAlgebra("ntffit")
)
#Run MX
start <- proc.time()[1]; # record start time
fit <- mxRun(model)
endTime <- proc.time()[1]; (endTime -start)
#Look at results
estimate <- fit$output$estimate
round(estimate,3)
# TODO include simulation output
# compare to simulation
# 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))
# dimnames(estimate.mat) <- list(c('OpenMx-Estimated','Simulated'),c('Var.E','Var.F','Var.A','Var.S','Var.D','Cor.Sps','Var.Phen'))
# estimate == estimate.mat
#NOTE: the minor difference between simulated & estimated parameters are to be expected & are due to sampling error
# > res.mat
# Var.E Var.F Var.A Var.S Var.D Cor.Sps Var.Phen
#OpenMx-Estimated 0.443 0.100 0.294 0.149 0 0.209 0.986
#Old.mx 0.442 0.099 0.294 0.149 0 0.206 0.986
#Simulated 0.433 0.100 0.304 0.197 0 0.200 1.032
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.