# ----------------------------------------------------------------------- # Program: UnivariateTwinAnalysis_MatrixRawCon.R # Author: Hermine Maes & Matt Keller # Date: March 3, 2010 # # Univariate Twin Analysis model to estimate causes of variation (ADCE) # for continuous data # Matrix style model input - Raw data input # # Revision History # Matt Keller -- March 3, 2010, added D to an ACE script to demonstrate # identification issues # ----------------------------------------------------------------------- # Starting things up # ----------------------------------------------------------------------- # Where is the working directory? getwd() #If it's not where you want it to be, change it (Under File) # Get the two packages we'll use require(OpenMx) # Get some additional functions written by Hermine & Tim # **NOTE**: you can read in an R source file from the web!! Cool! source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R") #Get our dataset from the web. #This is a cool R trick: your file might exist on your local machine #OR it might exist on the web! con <- "http://www.matthewckeller.com/public/ntf.data" ntf.data <- read.table(con,header=TRUE) # ----------------------------------------------------------------------- # Prepare Data # ----------------------------------------------------------------------- head(ntf.data) str(ntf.data) #change -999 to NAs ntf.data[ntf.data==-999] <- NA #define variables nv <- 1 ntv <- nv*2 mzData <- ntf.data[ntf.data$twin.type=='MZ',3:4] dzData <- ntf.data[ntf.data$twin.type=='DZ',3:4] dim(mzData) dim(dzData) # ----------------------------------------------------------------------- # Print Descriptive Statistics # ----------------------------------------------------------------------- cov(mzData,use="complete") #MZ Var/Covar matrix cov(dzData,use="complete") #DZ Var/Covar matrix # ----------------------------------------------------------------------- # Get ****RANDOM**** start values # ----------------------------------------------------------------------- a.start <- runif(1,-1,1) d.start <- runif(1,-1,1) c.start <- runif(1,-1,1) #Look at your own personal start values a.start d.start c.start # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 1 # ----------------------------------------------------------------------- # Fit ADCE Model with RawData and Matrices Input # ----------------------------------------------------------------------- univADCEModel <- mxModel("univADCE", mxModel("ADCE", # Matrices a, c, and e to store a, c, and e path coefficients # ***NOTE*** IF a, d, and c are all free=TRUE, # your model will not be identified!! mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=a.start, label="a11", name="a" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=d.start, label="d11", name="d" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c.start, label="c11", name="c" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=0, label="e11", name="e" ), # Matrices A, C, and E compute variance components mxAlgebra( expression=a %*% t(a), name="A" ), mxAlgebra( expression=d %*% t(d), name="D" ), mxAlgebra( expression=c %*% t(c), name="C" ), mxAlgebra( expression=e %*% t(e), name="E" ), # Algebra to compute ***total variance*** mxAlgebra( expression=A+D+C+E, name="TotVar" ), # Matrix & Algebra for expected ***means*** vector mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 20, label="mean", name="Mean" ), mxAlgebra( expression= cbind(Mean,Mean), name="expMean"), # Algebra for expected variance/covariance matrix in ***MZ*** mxAlgebra( expression= rbind( cbind(A+D+C+E , A+D+C), cbind(A+D+C , A+D+C+E)),name="expCovMZ" ), # Algebra for expected variance/covariance matrix in ***DZ*** mxAlgebra( expression= rbind( cbind(A+D+C+E , 0.5%x%A+0.25%x%D+C), cbind(0.5%x%A+0.25%x%D+C , A+D+C+E)), name="expCovDZ" ) ), mxModel("MZ", mxData( observed=mzData, type="raw" ), mxFIMLObjective( covariance="ADCE.expCovMZ", means="ADCE.expMean", dimnames=c('tw1','tw2') ) ), mxModel("DZ", mxData( observed=dzData, type="raw" ), mxFIMLObjective( covariance="ADCE.expCovDZ", means="ADCE.expMean", dimnames=c('tw1','tw2') ) ), mxAlgebra( expression=MZ.objective + DZ.objective, name="-2sumll" ), mxAlgebraObjective("-2sumll") ) univADCEFit <- mxRun(univADCEModel) summary(univADCEFit) # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # END PRACTICAL 1 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 2 - YOU need to do something to make this script work # ----------------------------------------------------------------------- # Fit ***ADE*** model # ----------------------------------------------------------------------- univADEModel <- mxModel(univADCEFit, name="univADE", mxModel(univADCEFit$ADCE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) # drop c at 0 ) ) univADEFit <- mxRun(univADEModel) summary(univADEFit) # ----------------------------------------------------------------------- # Fit ***ACE*** model # ***NOTE: we normally wouldn't run an ACE model bc CovMz > 2CovDz # Can anyone notice a difference in the -2LL??? # ----------------------------------------------------------------------- univACEModel <- mxModel(univADCEFit, name="univACE", mxModel(univADCEFit$ADCE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="d11", name="d" ) # drop d at 0 ) ) univACEFit <- mxRun(univACEModel) summary(univACEFit) # ----------------------------------------------------------------------- # Fit ***DCE*** model # ***NOTE: We are NOT recommeding this model, but showing that it is possible # in order to demonstrate the concept of parameter indeterminacy # ----------------------------------------------------------------------- univDCEModel <- #YOU DO IT! Figure out how to run a DCE Model univDCEFit <- mxRun(univDCEModel) summary(univDCEFit) # ----------------------------------------------------------------------- tableFitStatistics(univADCEFit,list(univADEFit,univACEFit,univDCEFit)) # ----------------------------------------------------------------------- # END PRACTICAL 2 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 3 - YOU need to do something to make this work # ----------------------------------------------------------------------- # What would have happened if, in the real world, C = .05, .10, or .15 rather # than being its assumed value of 0?? # NOTE: YOU need to figure out what to change! # HINT: the path is the sqrt(x) of the variance. So if you want to fix the # variance of C = .05, the path coefficient should be sqrt(.05). # Fit ADE model where C is assumed to be ***.05*** # ----------------------------------------------------------------------- univADEModel.05 <- mxModel(univADCEFit, name="univADE", mxModel(univADCEFit$ADCE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) # drop c at .05 ) ) univADEFit.05 <- mxRun(univADEModel.05) summary(univADEFit.05) # ----------------------------------------------------------------------- # Fit ADE model where C is assumed to be ***.10*** # NOTE: YOU need to figure out what to change! # ----------------------------------------------------------------------- univADEModel.10 <- mxModel(univADCEFit, name="univADE", mxModel(univADCEFit$ADCE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) # drop c at .10 ) ) univADEFit.10 <- mxRun(univADEModel.10) summary(univADEFit.10) # ----------------------------------------------------------------------- # Fit ADE model where C is assumed to be ***.15*** # NOTE: YOU need to figure out what to change! # ----------------------------------------------------------------------- univADEModel.15 <- mxModel(univADCEFit, name="univADE", mxModel(univADCEFit$ADCE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) # drop c at .15 ) ) univADEFit.15 <- mxRun(univADEModel.15) summary(univADEFit.15) # ----------------------------------------------------------------------- # Compare the results of ***A-hat*** across assumed values of C: univADEFit@submodels$ADCE$a@values^2 univADEFit.05@submodels$ADCE$a@values^2 univADEFit.10@submodels$ADCE$a@values^2 univADEFit.15@submodels$ADCE$a@values^2 # Compare the results of ***D-hat*** across assumed values of C: univADEFit@submodels$ADCE$d@values^2 univADEFit.05@submodels$ADCE$d@values^2 univADEFit.10@submodels$ADCE$d@values^2 univADEFit.15@submodels$ADCE$d@values^2 # Compare the -2LLs of all these models univADEFit@algebras$'-2sumll'@result univADEFit.05@algebras$'-2sumll'@result univADEFit.10@algebras$'-2sumll'@result univADEFit.15@algebras$'-2sumll'@result # ----------------------------------------------------------------------- # END PRACTICAL 3 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 4 # ----------------------------------------------------------------------- #Create the data mz.ntf <- ntf.data[ntf.data$twin.type=='MZ',3:6] dz.ntf <- ntf.data[ntf.data$twin.type=='DZ',3:6] #Fit NTF Model with **Sibling Env** #NOTE: This is a different way of writing MX models than you're used to #----------------------------------------------------------------------- # Paths going to latent variables a <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.6, label="AddGenPath",name="a") d <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.6, label="DomPath",name="d") s <- mxMatrix(type="Lower", nrow=nv, ncol=nv, free=T, values=.4, label="SibPath", name="s") e <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.6, label="EnvPath",name="e") mu <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.0, label="AMCopath",name="mu") # 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.ntf)), 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.ntf),colnames(mz.ntf)),name="expCovMz"), mxData(observed=mz.ntf, 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.ntf)), 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.ntf),colnames(dz.ntf)),name="expCovDz"), mxData(observed=dz.ntf, 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) # ----------------------------------------------------------------------- # END PRACTICAL 4 # -----------------------------------------------------------------------