# ----------------------------------------------------------------------- # Program: NTF.ASDE.R # Author: Matthew Keller & Sarah Medland # Date: Sept 16 2009 # # Nuclear Twin Family model to estimate causes of variation for # continuous data # Matrix style model input - Raw data input # # Revision History # Matthew Keller -- 02 17 2010 updated & reformatted # # NOTE: NTF design does not allow Vs & Vf to be estimated simultaneously for identifiability; # This script only estimates Vs # # ----------------------------------------------------------------------- require(OpenMx) #Set working directory setwd("~/Documents/RESEARCH/MxR/OpenMx") #Get nice functions written by Maes & York source("GenEpiHelperFunctions.R") #Get the raw data we'll be working with con <- url("http://matthewckeller.com/public/ntf.datasets.RData") showConnections(all=TRUE) load(con) close(con) # Look at Data # ----------------------------------------------------------------------- summary(mz) summary(dz) nv <- 1 # Print Descriptive Statistics # ----------------------------------------------------------------------- summary(mzData) colMeans(mzData,na.rm=T) cov(mzData,use="complete") summary(dzData) colMeans(dzData,na.rm=T) cov(dzData,use="complete") #Fit NTF Model with **Sibling Env** #----------------------------------------------------------------------- # Paths going to latent variables a <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.5,label="AddGenPath",name="a") d <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.1,label="DomPath",name="d") e <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.7,label="EnvPath",name="e") mu <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.2,label="AMCopath",name="mu") s <- mxMatrix(type="Lower", nrow=nv, ncol=nv, free=T, values=.2, label="SibPath", name="s") # Matrices for latent variances Vp1 <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=T,values=1.5,label="VarPhen",name="Vp1") delta1 <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=T,values=.5,label="CovPhenGen",name="delta1") q1 <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=T,values=1,label="LatentVarAddGen",name="q1") # mxAlgebra - nonlinear constraints Vp2 <- mxAlgebra(a %&% q1 + e %*% t(e) + d %*% t(d) + s %*% t(s),name="Vp2") delta2 <- mxAlgebra(q1 %*% a,name="delta2") q2 <- mxAlgebra(1 + delta1 %&% mu,name="q2") # Equating nonlinear constraints and parameters VpCon <- mxConstraint('Vp1',"=",'Vp2',name='VpCon') deltaCon <- mxConstraint('delta1',"=",'delta2',name='deltaCon') qCon <- mxConstraint('q1',"=",'q2',name='qCon') # mxAlgebra - relative covariances CvMz <- mxAlgebra(a %&% q1 + d %*% t(d) + s %*% t(s),name="CvMz") CvDz <- mxAlgebra(a %&% (q1-.5) + .25 %x% d %*% t(d) + s %*% t(s),name="CvDz") ParChild <- mxAlgebra(.5 %x% a %*% delta1 + .5 %x% a %*% delta1 %*% mu %*% Vp1,name="ParChild") CvSps <- mxAlgebra(Vp1 %&% mu,name="CvSps") # Put all these pathways, variances, and constraints together ntf.pths <- mxModel(model="NTFpaths",a,d,e,mu,s,Vp1,delta1,q1,Vp2,delta2,q2,VpCon,deltaCon,qCon,CvMz,CvDz,ParChild,CvSps) #make the expected mean and covariance matrices for MZ data and pull the raw MZ data in mzModel <- mxModel(model=ntf.pths,name = "MZNTF", mxMatrix(type="Full", nrow=1, ncol=4*nv, free=TRUE, values= .05, label="mean", dimnames=list(NULL, colnames(mz)), name="expMeanMz"), # Algebra for expected variance/covariance matrix in MZ mxAlgebra(rbind( cbind(Vp1, CvMz, ParChild, ParChild), cbind(CvMz, Vp1, ParChild, ParChild), cbind(ParChild, ParChild, Vp1, CvSps), cbind(ParChild, ParChild, CvSps, Vp1)), dimnames=list(colnames(mz),colnames(mz)),name="expCovMz"), mxData(observed=mz, type="raw"), mxFIMLObjective(covariance="expCovMz",means="expMeanMz")) #make the expected mean and covariance matrices for DZ data and pull the raw DZ data in dzModel <- mxModel(model=ntf.pths,name = "DZNTF", mxMatrix(type="Full", nrow=1, ncol=4*nv, free=TRUE, values= .05, label="mean", dimnames=list(NULL, colnames(dz)), name="expMeanDz"), # Algebra for expected variance/covariance matrix in DZ mxAlgebra(rbind( cbind(Vp1, CvDz, ParChild, ParChild), cbind(CvDz, Vp1, ParChild, ParChild), cbind(ParChild, ParChild, Vp1, CvSps), cbind(ParChild, ParChild, CvSps, Vp1)), dimnames=list(colnames(dz),colnames(dz)),name="expCovDz"), mxData(observed=dz, type="raw"), mxFIMLObjective(covariance="expCovDz",means="expMeanDz")) #create the mxModel that defines the objective function ASDE <- mxModel(model="ASDE", mzModel, dzModel, #MZNTF.objective is the automatic name for the -2LL of mzModel mxAlgebra(MZNTF.objective + DZNTF.objective, name="-2LL"), mxAlgebraObjective("-2LL")) #Run MX ASDE.Fit <- mxRun(ASDE) summary(ASDE.Fit) #Look at results res <- ASDE.Fit@output$estimate round(res,3) #drop a at 0 (you have to do it in 3 places because a exists in 3 places) SDE <- mxModel(ASDE.Fit,name="SDE", mxModel(ntf.pths,mxMatrix(type="Lower",nrow=nv,ncol=nv,free=F,values=0,label="AddGenPath",name="a")), mxModel(dzModel,mxMatrix(type="Lower",nrow=nv,ncol=nv,free=F,values=0,label="AddGenPath",name="a")), mxModel(mzModel,mxMatrix(type="Lower",nrow=nv,ncol=nv,free=F,values=0,label="AddGenPath",name="a"))) SDE.Fit <- mxRun(SDE) summary(SDE.Fit)$Minus2LogLikelihood - summary(ASDE.Fit)$Minus2LogLikelihood ##OLD STUFF DOWN HERE #Fit NTF Model with **Familial Env** #----------------------------------------------------------------------- # Paths common in all models a <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=TRUE,values=.5,label="AddGenPath",name="a"), d <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=TRUE,values=.1,label="DomPath",name="d"), e <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=TRUE,values=.7,label="EnvPath",name="e"), mu <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=TRUE,values=.2,label="AMCopath",name="mu"), # Paths for familial env m <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=FALSE,values=0,label="FamPath",name="m"), x1 <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=FALSE,values=0,label="LatentVarF",name="x1"), w1 <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=FALSE,values=0,label="CovFA",name="w1"), #fix s=0 for Vs=0 mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=.2, label="SibPath", name="s"), #Matrices for latent variances mxMatrix(type="Full",nrow=nv,ncol=nv,free=TRUE,values=1.5,label="VarPhen",name="Vp1"), mxMatrix(type="Full",nrow=nv,ncol=nv,free=TRUE,values=.5,label="CovPhenGen",name="delta1"), mxMatrix(type="Full",nrow=nv,ncol=nv,free=TRUE,values=1,label="LatentVarAddGen",name="q1"), #this parameter should always be fixed @1 regardless mxMatrix(type="Full",nrow=nv,ncol=nv,free=FALSE,values=1,label="LatentFamilialVar",name="f"), #mxAlgebra section - nonlinear constraints mxAlgebra(a %&% q1 + f %&% x1 + 2 %x% a %*% w1 %*% t(f) + e %*% t(e) + d %*% t(d) + s %*% t(s),name="Vp2"), mxAlgebra(2 %x% m %&% Vp1 + 2 %x% m %*% Vp1 %*% mu %*% t(Vp1) %*% t(m),name="x2"), mxAlgebra(q1 %*% a + w1 %*% f,name="delta2"), mxAlgebra(1 + delta1 %&% mu,name="q2"), mxAlgebra(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(a %&% q1 + f %&% x1 + 2 %x% a %*% w1 %*% t(f) + d %*% t(d) + s %*% t(s),name="CvMz"), mxAlgebra(a %&% (q1-.5) + .25 %x% d %*% t(d) + f %&% x1 + 2 %x% a %*% w1 %*% t(f)+ s %*% t(s),name="CvDz"), mxAlgebra(.5 %x% a %*% delta1 + .5 %x% a %*% delta1 %*% mu %*% Vp1 + m %*% Vp1 +m %*% Vp1 %&% mu,name="ParChild"), mxAlgebra(Vp1 %&% mu,name="CvSps")) #Look at results summary(NTFModelFit) res <- NTFModelFit@output$estimate round(res,3) RESULTS <- round(c(c(res[c(1:3,5)]^2,1) * c(res[9],rep(1,3),res[7]),res[4]*res[6],2*res[1]*res[10],res[6]),3) names(RESULTS) <- c('var.A','var.D','var.E','var.S','Var.F','Cor.Sps','CovAF','Var.P') #look at implied vs observed Covariances round(cov(mz,use='pairwise'),3) #observed round(NTFModelFit@output$algebras$MZNTF.expCovMz,3) #implied round(cov(dz,use='pairwise'),3) #observed round(NTFModelFit@output$algebras$DZNTF.expCovDz,3) #implied #compare to old Mx results Old.mx <- c(.442,.099,.294,.149,0,.206,0,.986) res.mat <- rbind(RESULTS,Old.mx) res.mat #drop A # var.A var.D var.E var.S Var.F Cor.Sps CovAF Var.P #OpenMx-Estimated 0.443 0.100 0.294 0.149 0 0.209 0 0.986 #Simulated 0.433 0.100 0.304 0.197 0 0.200 0 1.032 #Old.mx 0.442 0.099 0.294 0.149 0 0.206 0 0.986