# ----------------------------------------------------------------------- # 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 # Matt Keller -- March 4, 2014, changed OpenMx syntax to make it more 'objective' # and also included additional assumed values of C # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 1 - you do NOT need to change anything to make this work # ----------------------------------------------------------------------- # Starting things up # ----------------------------------------------------------------------- # Where is the working directory? getwd() #If it's not where you want it to be, change it (Under File) # or just use setwd("/home/xxx/xxx") function where the '/home/xxx/xxx' is the full path # 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) #get summary of each column in ntf.data summary(ntf.data) #YIKES - look at those big -999's in there - those are NAs #change -999 to NA (could have done this using na.strings argument in read.table) ntf.data[ntf.data==-999] <- NA #make sure changing the NA's worked: summary(ntf.data) #define variables nv <- 1 #number of variables; allows us to make these models multivariate with ease 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 & look at data # ----------------------------------------------------------------------- #Look at your data ! x11() op <- par(mfrow=c(2,1)) plot(mzData$tw1,mzData$tw2) plot(dzData$tw1,dzData$tw2) par(op) #Look at covariances cov(mzData,use="complete") #MZ Var/Covar 2x2 matrix cov(dzData,use="complete") #DZ Var/Covar 2x2 matrix #create object of MZ vs DZ cov's (cov.MZ <- cov(mzData,use="complete") [1,2]) (cov.DZ <- cov(dzData,use="complete") [1,2]) # ----------------------------------------------------------------------- # Get ****RANDOM**** start values # For identified models, start values should not affect the final best fit # ----------------------------------------------------------------------- (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 #Ask your neighbor what they got! # ----------------------------------------------------------------------- # Fit ADCE Model with RawData and Matrices Input # ----------------------------------------------------------------------- # 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! Go ahead and run this non-id'd model! a.mat <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,values=a.start, label="a11", name="a" ) d.mat <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,values=d.start, label="d11", name="d" ) c.mat <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,values=c.start, label="c11", name="c" ) e.mat <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,values=0, label="e11", name="e" ) # Matrices A, C, D, and E compute variance components A.mat <- mxAlgebra( expression=a %*% t(a), name="A" ) D.mat <- mxAlgebra( expression=d %*% t(d), name="D" ) C.mat <- mxAlgebra( expression=c %*% t(c), name="C" ) E.mat <- mxAlgebra( expression=e %*% t(e), name="E" ) # Matrix & Algebra for expected ***means*** vector Mn.mat <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE,values= 20, label="mean", name="Mean" ) MN <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean") # Algebra for expected variance/covariance matrix in ***MZ*** MZ.Cov <- 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*** DZ.Cov <- 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" ) # Put in the data dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) #Objective objects for two groups objMZ <- mxFIMLObjective(covariance="expCovMZ",means="expMean", dimnames=c('tw1','tw2')) objDZ <- mxFIMLObjective(covariance="expCovDZ",means="expMean", dimnames=c('tw1','tw2')) #Combine groups params <- list(a.mat,d.mat,c.mat,e.mat,A.mat,D.mat,C.mat,E.mat) modelMZ <- mxModel("MZ",params,Mn.mat,MN,MZ.Cov,dataMZ,objMZ) modelDZ <- mxModel("DZ",params,Mn.mat,MN,DZ.Cov,dataDZ,objDZ) minus2ll <- mxAlgebra(MZ.objective+DZ.objective,name='minus2sumLL') obj <- mxAlgebraObjective('minus2sumLL') ADCE.Model <- mxModel("ADCE",modelMZ,modelDZ,minus2ll,obj) #Run the model ADCE.Fit <- mxRun(ADCE.Model,intervals=FALSE) (ADCE.Sum <- summary(ADCE.Fit)) #YOUR start values (compare with neighbor) a.start d.start c.start #Get the -2LL - ask your neighbor if theirs is the same! ADCE.Sum$Minus2LogLikelihood #Estimate of VA - ask your neighbors what they got! ADCE.Fit$MZ$a@values^2 #Estimate of VD - ask your neighbors what they got! ADCE.Fit$MZ$d@values^2 #Estimate of VC - ask your neighbors what they got! ADCE.Fit$MZ$c@values^2 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # END PRACTICAL 1 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 2 - You DO need to change the ???'s to make this script work # ----------------------------------------------------------------------- # Fit ***ADE*** model by dropping C # ----------------------------------------------------------------------- NEW.c.start <- 0 # Go back and pickup the model so that we can run significance tests ADE.Model <- ADCE.Fit ADE.Model <- omxSetParameters(ADE.Model, label="c11", free=???, values=???) ADE.Fit <- mxRun(ADE.Model, intervals=FALSE) (ADE.Sum <- summary(ADE.Fit)) #Get the -2LL - Has it changed from your last model? Why or why not? ADE.Sum$Minus2LogLikelihood ADCE.Sum$Minus2LogLikelihood #your previous -2LL #Estimate of VA - ask your neighbors what they got! ADE.Fit$MZ$a@values^2 #Estimate of VD - ask your neighbors what they got! ADE.Fit$MZ$d@values^2 #Estimate of VC - zero by definition ADE.Fit$MZ$c@values^2 # ----------------------------------------------------------------------- # Fit ***ACE*** model # ***NOTE: we normally wouldn't run an ACE model bc CovMz > 2CovDz # Do you notice a difference in the -2LL? Why?? # ----------------------------------------------------------------------- NEW.d.start <- 0 # Go back and pickup the model so that we can run significance tests ACE.Model <- ADCE.Fit ACE.Model <- omxSetParameters(ACE.Model, label="???", free=FALSE, values=NEW.d.start ) ACE.Fit <- mxRun(ACE.Model, intervals=FALSE) (ACE.Sum <- summary(ACE.Fit)) #Get the -2LL - Has it changed from your last model? Why or why not? ACE.Sum$Minus2LogLikelihood #Previous -2LL's ADE.Sum$Minus2LogLikelihood ADCE.Sum$Minus2LogLikelihood #Estimate of VA - ask your neighbors what they got! ACE.Fit$MZ$a@values^2 #Estimate of VD - zero by definition ACE.Fit$MZ$d@values^2 #Estimate of VC - ask your neighbors what they got! ACE.Fit$MZ$c@values^2 # ----------------------------------------------------------------------- # Fit ***DCE*** model # ***NOTE: I am NOT recommeding this model. A phenotype with only D and # no A is VERY unlikely (for reasons Dr. Eaves discussed; interactions almost always have # main (or additive) effects associated with them). # Rather, the point is simply to but show that it is mathematically possible # in order to demonstrate the concept of parameter indeterminacy # ----------------------------------------------------------------------- NEW.a.start <- ??? # Go back and pickup the model so that we can run significance tests DCE.Model <- ADCE.Fit DCE.Model <- omxSetParameters(DCE.Model, label="a11", free=FALSE, values=NEW.a.start ) DCE.Fit <- mxRun(DCE.Model, intervals=FALSE) (DCE.Sum <- summary(DCE.Fit)) #Get the -2LL - Has it changed from your other models? Why or why not? DCE.Sum$Minus2LogLikelihood #Previous -2LL's ADCE.Sum$Minus2LogLikelihood ADE.Sum$Minus2LogLikelihood ACE.Sum$Minus2LogLikelihood #Estimate of VA - zero by definition DCE.Fit$MZ$a@values^2 #Estimate of VD - ask your neighbors what they got! DCE.Fit$MZ$d@values^2 #Estimate of VC - ask your neighbors what they got! DCE.Fit$MZ$c@values^2 # ----------------------------------------------------------------------- #compare the fits of the four models run RESULTS <- data.frame(name=c("rand.ADCE","ADE","ACE","DCE"),neg2ll=c(ADCE.Sum$Minus2LogLikelihood,ADE.Sum$Minus2LogLikelihood,ACE.Sum$Minus2LogLikelihood,DCE.Sum$Minus2LogLikelihood), VA=c(ADCE.Fit$MZ$a@values^2,ADE.Fit$MZ$a@values^2,ACE.Fit$MZ$a@values^2,DCE.Fit$MZ$a@values^2),VD=c(ADCE.Fit$MZ$d@values^2,ADE.Fit$MZ$d@values^2,ACE.Fit$MZ$d@values^2,DCE.Fit$MZ$d@values^2),VC=c(ADCE.Fit$MZ$c@values^2,ADE.Fit$MZ$c@values^2,ACE.Fit$MZ$c@values^2,DCE.Fit$MZ$c@values^2)) RESULTS # ----------------------------------------------------------------------- # END PRACTICAL 2 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 3 - You do NOT need to change anything to make this work # ----------------------------------------------------------------------- # What would have happened if, in the real world, C = .05, .10, .15, or .3 rather # than being its assumed value of 0?? # Fit ADE model where C is assumed to be ***.05*** # NOTE: YOU need to figure out what to change! # HINT: the path is the sqrt(x) of the variance. So if you want to simulate # variance of .05, the path should be sqrt(.05). # ----------------------------------------------------------------------- NEW.c.start <- sqrt(.05) # Go back and pickup the model so that we can run significance tests ADE.fixC05.Model <- ADCE.Fit ADE.fixC05.Model <- omxSetParameters(ADE.fixC05.Model, label="c11", free=FALSE, values=NEW.c.start ) ADE.fixC05.Fit <- mxRun(ADE.fixC05.Model, intervals=FALSE) (ADE.fixC05.Sum <- summary(ADE.fixC05.Fit)) #Get the -2LL - Has it changed from your last ADCE model? Why or why not? ADE.fixC05.Sum$Minus2LogLikelihood #Previous -2LL's DCE.Sum$Minus2LogLikelihood ADCE.Sum$Minus2LogLikelihood ADE.Sum$Minus2LogLikelihood ACE.Sum$Minus2LogLikelihood #Estimate of VA - ask your neighbors what they got! ADE.fixC05.Fit$MZ$a@values^2 #Estimate of VD - ask your neighbors what they got! ADE.fixC05.Fit$MZ$d@values^2 #Estimate of VC - .05 by definition ADE.fixC05.Fit$MZ$c@values^2 # ----------------------------------------------------------------------- # Fit ADE model where C is assumed to be ***.10*** # NOTE: YOU need to figure out what to change! # ----------------------------------------------------------------------- NEW.c.start <- sqrt(.10) # Go back and pickup the model so that we can run significance tests ADE.fixC10.Model <- ADCE.Fit ADE.fixC10.Model <- omxSetParameters(ADE.fixC10.Model, label="c11", free=FALSE, values=NEW.c.start ) ADE.fixC10.Fit <- mxRun(ADE.fixC10.Model, intervals=FALSE) (ADE.fixC10.Sum <- summary(ADE.fixC10.Fit)) #Get the -2LL - Has it changed from your last ADCE model? Why or why not? ADE.fixC10.Sum$Minus2LogLikelihood #Previous -2LL's for other values of C ADE.fixC05.Sum$Minus2LogLikelihood #Estimate of VA - ask your neighbors what they got! ADE.fixC10.Fit$MZ$a@values^2 #Estimate of VD - ask your neighbors what they got! ADE.fixC10.Fit$MZ$d@values^2 #Estimate of VC - .05 by definition ADE.fixC10.Fit$MZ$c@values^2 # ----------------------------------------------------------------------- # Fit ADE model where C is assumed to be ***.15*** # NOTE: YOU need to figure out what to change! # ----------------------------------------------------------------------- NEW.c.start <- sqrt(.15) # Go back and pickup the model so that we can run significance tests ADE.fixC15.Model <- ADCE.Fit ADE.fixC15.Model <- omxSetParameters(ADE.fixC15.Model, label="c11", free=FALSE, values=NEW.c.start ) ADE.fixC15.Fit <- mxRun(ADE.fixC15.Model, intervals=FALSE) (ADE.fixC15.Sum <- summary(ADE.fixC15.Fit)) #Get the -2LL - Has it changed from your last ADCE model? Why or why not? ADE.fixC15.Sum$Minus2LogLikelihood #Previous -2LL's for other values of C ADE.fixC10.Sum$Minus2LogLikelihood ADE.fixC05.Sum$Minus2LogLikelihood #Estimate of VA - ask your neighbors what they got! ADE.fixC15.Fit$MZ$a@values^2 #Estimate of VD - ask your neighbors what they got! ADE.fixC15.Fit$MZ$d@values^2 #Estimate of VC - .05 by definition ADE.fixC15.Fit$MZ$c@values^2 # ----------------------------------------------------------------------- # Fit ADE model where C is assumed to be ***.30*** # NOTE: YOU need to figure out what to change! # ----------------------------------------------------------------------- NEW.c.start <- sqrt(.30) # Go back and pickup the model so that we can run significance tests ADE.fixC30.Model <- ADCE.Fit ADE.fixC30.Model <- omxSetParameters(ADE.fixC30.Model, label="c11", free=FALSE, values=NEW.c.start ) ADE.fixC30.Fit <- mxRun(ADE.fixC30.Model, intervals=FALSE) (ADE.fixC30.Sum <- summary(ADE.fixC30.Fit)) #Get the -2LL - Has it changed from your last ADCE model? Why or why not? ADE.fixC30.Sum$Minus2LogLikelihood #Previous -2LL's for other values of C ADE.fixC15.Sum$Minus2LogLikelihood ADE.fixC10.Sum$Minus2LogLikelihood ADE.fixC05.Sum$Minus2LogLikelihood #Estimate of VA - ask your neighbors what they got! ADE.fixC30.Fit$MZ$a@values^2 #Estimate of VD - ask your neighbors what they got! ADE.fixC30.Fit$MZ$d@values^2 #Estimate of VC - .05 by definition ADE.fixC30.Fit$MZ$c@values^2 # ----------------------------------------------------------------------- #compare the fits of the four models run RESULTS3 <- data.frame(name=c("ADE","ADE.C05","ADE.C10","ADE.C15","ADE.C30"),neg2ll=c(ADE.Sum$Minus2LogLikelihood,ADE.fixC05.Sum$Minus2LogLikelihood,ADE.fixC10.Sum$Minus2LogLikelihood,ADE.fixC15.Sum$Minus2LogLikelihood,ADE.fixC30.Sum$Minus2LogLikelihood),VA=c(ADE.Fit$MZ$a@values^2,ADE.fixC05.Fit$MZ$a@values^2,ADE.fixC10.Fit$MZ$a@values^2,ADE.fixC15.Fit$MZ$a@values^2,ADE.fixC30.Fit$MZ$a@values^2),VD=c(ADE.Fit$MZ$d@values^2,ADE.fixC05.Fit$MZ$d@values^2,ADE.fixC10.Fit$MZ$d@values^2,ADE.fixC15.Fit$MZ$d@values^2,ADE.fixC30.Fit$MZ$d@values^2),VC=c(ADE.Fit$MZ$c@values^2,ADE.fixC05.Fit$MZ$c@values^2,ADE.fixC10.Fit$MZ$c@values^2,ADE.fixC15.Fit$MZ$c@values^2,ADE.fixC30.Fit$MZ$c@values^2)) RESULTS3 #VISUALIZING NON-IDENTIFIED MODELS! #Look at a 3d plot of 2LL on y as function of fixed values of VA, VC, and VD (VD as a color here) #The R code for setting up this data is below, just for your perusal # ----------------------------------------------------------------------- load("CTD.ACDE-param.indet.plot_2014.RData") #save(SOLS3,d.colors,D.scaled,file="CTD.ACDE-param.indet.plot_2014.RData") #Load in library to allow for 3d interactive plots require(rgl) #Look at VA vs. VC (color coding VD) plot3d(SOLS3$VC,SOLS3$VA,SOLS3$two.ll,col=d.colors[D.scaled],type='s',size=.5,xlab='VC.est',ylab='VA.est',zlab='2LL') # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # END PRACTICAL 3 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 4 (you do NOT need to do anything to make this work) # ----------------------------------------------------------------------- #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] #Look at your MZ data - run just the next 6 lines and look (and think!) x11() op <- par(mfrow=c(2,2)) plot(mz.ntf$tw1,mz.ntf$tw2,main=paste("r = ",round(cor(mz.ntf$tw1,mz.ntf$tw2,use='pairwise'),2))) plot(mz.ntf$tw1,mz.ntf$father,main=paste("r = ",round(cor(mz.ntf$tw1,mz.ntf$father,use='pairwise'),2))) plot(mz.ntf$tw1,mz.ntf$mother,main=paste("r = ",round(cor(mz.ntf$tw1,mz.ntf$mother,use='pairwise'),2))) plot(mz.ntf$father,mz.ntf$mother,main=paste("r = ",round(cor(mz.ntf$father,mz.ntf$mother,use='pairwise'),2))) par(op) #Look at your DZ data - run just the next 6 lines and look (and think!) x11() op <- par(mfrow=c(2,2)) plot(dz.ntf$tw1,dz.ntf$tw2,main=paste("r = ",round(cor(dz.ntf$tw1,dz.ntf$tw2,use='pairwise'),2))) plot(dz.ntf$tw1,dz.ntf$father,main=paste("r = ",round(cor(dz.ntf$tw1,dz.ntf$father,use='pairwise'),2))) plot(dz.ntf$tw1,dz.ntf$mother,main=paste("r = ",round(cor(dz.ntf$tw1,dz.ntf$mother,use='pairwise'),2))) plot(dz.ntf$father,dz.ntf$mother,main=paste("r = ",round(cor(dz.ntf$father,dz.ntf$mother,use='pairwise'),2))) par(op) #Fit NTF Model with **Sibling Env** (instead of "C" - both ~ the same) #----------------------------------------------------------------------- # 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 these relative covariances together into MZ relatives and DZ relatives matrices MZ.rel.cv <- 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="expCovMZRels") DZ.rel.cv <- 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="expCovDZRels") #Put the relatives means together into means matrix (same for MZ and DZ rels) Both.means <- mxMatrix(type="Full", nrow=1, ncol=4*nv, free=TRUE, values= .05, label="mean",dimnames=list(NULL, colnames(mz.ntf)), name="expMeanBoth") # 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) # Put in the data dataMZ <- mxData( observed=mz.ntf, type="raw" ) dataDZ <- mxData( observed=dz.ntf, type="raw" ) #Objective objects for two groups objMZRels <- mxFIMLObjective(covariance="expCovMZRels",means="expMeanBoth", dimnames=c('tw1','tw2','father','mother')) objDZRels <- mxFIMLObjective(covariance="expCovDZRels",means="expMeanBoth", dimnames=c('tw1','tw2','father','mother')) #Combine groups params <- list(a,d,e,mu,s,Vp1,delta1,q1,Vp2,delta2,q2,VpCon,deltaCon,qCon,CvMz,CvDz,ParChild,CvSps,MZ.rel.cv,Both.means,DZ.rel.cv) modelMZ <- mxModel("MZntf",params,dataMZ,objMZRels) modelDZ <- mxModel("DZntf",params,dataDZ,objDZRels) minus2ll <- mxAlgebra(MZntf.objective+DZntf.objective,name='minus2sumLL') obj <- mxAlgebraObjective('minus2sumLL') ADCE.NTF.Model <- mxModel("ntfADCE",modelMZ,modelDZ,minus2ll,obj) #Run the model ADCE.NTF.Fit <- mxRun(ADCE.NTF.Model,intervals=FALSE) (ADCE.NTF.Fit.Sum <- summary(ADCE.NTF.Fit)) # Do we notice a bias in these? # (Note, estimates never exactly equal true values due to sampling variability) # A (true simulated value was .30) ADCE.NTF.Fit$MZntf.a@values %*% ADCE.NTF.Fit$MZntf.q1@values %*% ADCE.NTF.Fit$MZntf.a@values # D (true simulated value was .30) ADCE.NTF.Fit$MZntf.d@values %*% ADCE.NTF.Fit$MZntf.d@values # S (true simulated value was .10) ADCE.NTF.Fit$MZntf.s@values %*% ADCE.NTF.Fit$MZntf.s@values # E (true simulated value was .30) ADCE.NTF.Fit$MZntf.e@values %*% ADCE.NTF.Fit$MZntf.e@values # ----------------------------------------------------------------------- # END PRACTICAL 4 # ----------------------------------------------------------------------- # ***NOTE*** this is commented out - the code is here just in case you're interested # This code builds the data that we used in the plot above (end practical 3) # ----------------------------------------------------------------------- #C.fix <- seq(0,.5,.025) #A.est <- D.est <- min2ll <- rep(NA,length(C.fix)) #SOLS <- cbind.data.frame(C.fix,A.est,D.est,min2ll) #dim(SOLS) #SOLS #Suppress NPSOL warnings #options(warn=-1) #for (IT in 1:nrow(SOLS)){ #New.Mod <- omxSetParameters(ADCE.Fit, label="c11", free=FALSE, values=sqrt(C.fix[IT])) #New.Fit <- mxRun(New.Mod, intervals=FALSE) #New.Fit.Sum <- summary(New.Fit) #SOLS$A.est[IT] <-New.Fit$MZ$a@values^2 #SOLS$D.est[IT] <-New.Fit$MZ$d@values^2 #SOLS$min2ll[IT] <- New.Fit.Sum$Minus2LogLikelihood #} #put warnings back on #options(warn=1) #SOLS$DminA <- round(SOLS$D.est - SOLS$A.est,2) # ----------------------------------------------------------------------- # Fit ADE model under multiple different fixed C's #This time, use results above but sample from a whole grid of C's, D's and A's # ----------------------------------------------------------------------- #SOLS2 <- expand.grid(SOLS$C.fix,SOLS$A.est) #ns <- expand.grid(SOLS$C.fix,SOLS$D.est) #SOLS3 <- cbind(SOLS2,ns[,2]) #names(SOLS3) <- c('VC','VA','VD') #SOLS3$min2ll <- NA #X <- matrix(NA,nrow=1,ncol=2) #write.table(X,file="min2ll") #Suppress warnings #options(warn=-1) #for (IT in 1:nrow(SOLS3)){ #print(IT) #New.Mod <- omxSetParameters(ADCE.Fit, label=c("c11","a11","d11"), free=FALSE, values=c(sqrt(SOLS3$VC[IT]),sqrt(SOLS3$VA[IT]),sqrt(SOLS3$VD[IT]))) #New.Fit <- mxRun(New.Mod, intervals=FALSE) #New.Fit.Sum <- summary(New.Fit) #SOLS3$min2ll[IT] <- round(New.Fit.Sum$Minus2LogLikelihood,4) #write.table(cbind(IT,round(New.Fit.Sum$Minus2LogLikelihood,4)),file='min2ll',append=TRUE,col.names=FALSE,row.names=FALSE) #} #put warnings back on #options(warn=1) #Create some variables for graphing #SOLS3$two.ll <- SOLS3$min2ll*-1 #SOLS3$AminD <- round(SOLS3$VA - SOLS3$VD,3) #d.colors <- rev(heat.colors(20)) #Color code by estimate of VD #SOLS3$D.scaled <- round((SOLS3$VD - min(SOLS3$VD))/(max(SOLS3$VD-min(SOLS3$VD))/20),0) # ----------------------------------------------------------------------- #save.image("CTD.ACDE-param.indet_2014.RData")