# ----------------------------------------------------------------------- # Program: UnivariateTwinAnalysis_MatrixRawCon.R # Author: Matt Keller & Hermine Maes # 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 # Matt Keller - March 6, 2018 - updated script to fit the variances directly and change # how we simulate data # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # PRACTICAL 3 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # begin practical 3A # ----------------------------------------------------------------------- #------------------ You might change this part of 3A later #Simulate TRUE VA, VC, VD, and VE PARAMETERS to make simulated twins VA <- .4 VD <- .2 VC <- .05 VE <- .35 #For simplicity, it's useful to have var's sum to 1 N.MZ <- 5e3 #5k MZ twins N.DZ <- 5e3 #5k DZ twins #------------------ Should not need to change anything in rest of 3A section # 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) require(MASS) #and make sure we're using the NPSOL optimizer mxOption(NULL,"Default optimizer","NPSOL") #Expected MZ and DZ covariances (SIM.MZ.COV <- VA+VD+VC) (SIM.DZ.COV <- .5*VA + .25*VD + VC) (MZ.VAR.COV.MAT <- matrix(c(1,SIM.MZ.COV,SIM.MZ.COV,1),nrow=2,byrow=T)) (DZ.VAR.COV.MAT <- matrix(c(1,SIM.DZ.COV,SIM.DZ.COV,1),nrow=2,byrow=T)) #Create the data - note "empirical=TRUE" means the var/cov matrix is #EXACTLY the values we gave it above. In real simulations, I wouldn't do this #(i.e., I'd say empirical=FALSE) to allow for sampling variability in the #variances and covariances, but for teaching, I make it exact here mzData <- mvrnorm(N.MZ,c(0,0),MZ.VAR.COV.MAT,empirical=TRUE) dzData <- mvrnorm(N.DZ,c(0,0),DZ.VAR.COV.MAT,empirical=TRUE) colnames(mzData) <- c('tw1','tw2') colnames(dzData) <- c('tw1','tw2') #Make sure it worked: cov(mzData) cov(dzData) # ----------------------------------------------------------------------- # begin practical 3A # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # begin practical 3B - Fit identified ADE model # ----------------------------------------------------------------------- # For identified models, start values should not affect the final best fit #We are fitting an ADE model (assuming VC=0), thus we should have an identified model # ----------------------------------------------------------------------- (VA.start <- runif(1,0,1)) (VD.start <- runif(1,0,1)) (VC.start <- 0) (VD.estimated <- TRUE) (VC.estimated <- FALSE) #look at c.mat below;"FALSE" means we fix VC to it's start value (0) #Look at your own personal start values and ask your neighbor what they got! # ----------------------------------------------------------------------- # Fit ADE 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-identified model! nv <- 1 a.mat <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=TRUE,values=VA.start, labels="VA", name="A" ) d.mat <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=VD.estimated,values=VD.start, labels="VD", name="D" ) c.mat <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=VC.estimated,values=VC.start, labels="VC", name="C" ) e.mat <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=TRUE,values=1, label="VE", name="E" ) # Matrix & Algebra for expected ***means*** vector Mn.mat <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE,values= .25, 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 <- mxExpectationNormal(covariance="expCovMZ",means="expMean", dimnames=c('tw1','tw2')) objDZ <- mxExpectationNormal(covariance="expCovDZ",means="expMean", dimnames=c('tw1','tw2')) #Combine groups myfitfun <- mxFitFunctionML() params <- list(a.mat,d.mat,c.mat,e.mat,myfitfun) modelMZ <- mxModel("MZ",params,Mn.mat,MN,MZ.Cov,dataMZ,objMZ) modelDZ <- mxModel("DZ",params,Mn.mat,MN,DZ.Cov,dataDZ,objDZ) obj <- mxFitFunctionMultigroup(c("MZ","DZ")) ADE.Model <- mxModel("ADE",modelMZ,modelDZ,obj) #Run the model ADE.Fit <- mxRun(ADE.Model,intervals=FALSE) (ADE.Sum <- summary(ADE.Fit)) #Get the -2LL - ask your neighbor if theirs is the same! ADE.Sum$Minus2LogLikelihood #Compare our Mx estimates to the algebriac expected ones: mxEval(A,model=ADE.Fit$MZ,T) (EXP.VAest.ADE <- 4*SIM.DZ.COV - SIM.MZ.COV) mxEval(D,model=ADE.Fit$MZ,T) (EXP.VDest.ADE <- 2*SIM.MZ.COV - 4*SIM.DZ.COV) mxEval(C,model=ADE.Fit$MZ,T) 0 #Obviously, we fixed VC to 0 so it's "estimated" at 0!! # ----------------------------------------------------------------------- # end practical 3B # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # begin practical 3C - Fit identified ACE model # ----------------------------------------------------------------------- # NOTE: we would NOT typically do this because covMZ > 2covDZ,but mathematically, it's fine #We are fitting an ACE model (assuming VD=0), thus we should have an identified model # ----------------------------------------------------------------------- (VA.start <- runif(1,0,1)) (VD.start <- 0) (VC.start <- runif(1,0,1)) (VD.estimated <- FALSE) #now we fix VD to 0 (VC.estimated <- TRUE) #and estimate VC #Look at your own personal start values and ask your neighbor what they got! # ----------------------------------------------------------------------- # Fit ACE 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-identified model! nv <- 1 a.mat <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE,values=VA.start, labels="VA", name="A" ) d.mat <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=VD.estimated,values=VD.start, labels="VD", name="D" ) c.mat <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=VC.estimated,values=VC.start, labels="VC", name="C" ) e.mat <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE,values=1, label="VE", name="E" ) # Matrix & Algebra for expected ***means*** vector Mn.mat <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE,values= .25, 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 <- mxExpectationNormal(covariance="expCovMZ",means="expMean", dimnames=c('tw1','tw2')) objDZ <- mxExpectationNormal(covariance="expCovDZ",means="expMean", dimnames=c('tw1','tw2')) #Combine groups myfitfun <- mxFitFunctionML() params <- list(a.mat,d.mat,c.mat,e.mat,myfitfun) modelMZ <- mxModel("MZ",params,Mn.mat,MN,MZ.Cov,dataMZ,objMZ) modelDZ <- mxModel("DZ",params,Mn.mat,MN,DZ.Cov,dataDZ,objDZ) obj <- mxFitFunctionMultigroup(c("MZ","DZ")) ACE.Model <- mxModel("ACE",modelMZ,modelDZ,obj) #Run the model ACE.Fit <- mxRun(ACE.Model,intervals=FALSE) (ACE.Sum <- summary(ACE.Fit)) #Compare LL of the ADE from b4 and ACE - they are the same when we allow negative variances! ADE.Sum$Minus2LogLikelihood ACE.Sum$Minus2LogLikelihood #Compare our Mx estimates to the algebriac expected ones: mxEval(A,model=ACE.Fit$MZ,T) (EXP.VAest.ACE <- 2*(SIM.MZ.COV-SIM.DZ.COV)) mxEval(C,model=ACE.Fit$MZ,T) (EXP.VCest.ACE <- 2*SIM.DZ.COV - SIM.MZ.COV) mxEval(D,model=ACE.Fit$MZ,T) 0 #Obviously, we fixed VD to 0 so it's "estimated" at 0 (even though in truth, it's .20) # ----------------------------------------------------------------------- # end practical 3C # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # begin practical 3D - Fit NON-identified ADCE model # ----------------------------------------------------------------------- # NOTE: we have more estimates than informative statistics. #Our model will NOT be identified. What happens to -2LL? What about the fits? # ----------------------------------------------------------------------- (VA.start <- runif(1,0,1)) (VD.start <- runif(1,0,1)) (VC.start <- runif(1,0,1)) (VD.estimated <- TRUE) #now we estimate VD (VC.estimated <- TRUE) #and estimate VC TOO!! #Look at your own personal start values and 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-identified model! nv <- 1 a.mat <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE,values=VA.start, labels="VA", name="A" ) d.mat <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=VD.estimated,values=VD.start, labels="VD", name="D" ) c.mat <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=VC.estimated,values=VC.start, labels="VC", name="C" ) e.mat <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE,values=1, label="VE", name="E" ) # Matrix & Algebra for expected ***means*** vector Mn.mat <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE,values= .25, 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 <- mxExpectationNormal(covariance="expCovMZ",means="expMean", dimnames=c('tw1','tw2')) objDZ <- mxExpectationNormal(covariance="expCovDZ",means="expMean", dimnames=c('tw1','tw2')) #Combine groups myfitfun <- mxFitFunctionML() params <- list(a.mat,d.mat,c.mat,e.mat,myfitfun) modelMZ <- mxModel("MZ",params,Mn.mat,MN,MZ.Cov,dataMZ,objMZ) modelDZ <- mxModel("DZ",params,Mn.mat,MN,DZ.Cov,dataDZ,objDZ) obj <- mxFitFunctionMultigroup(c("MZ","DZ")) ADCE.Model <- mxModel("ADCE",modelMZ,modelDZ,obj) #Run the model ADCE.Fit <- mxRun(ADCE.Model,intervals=FALSE) (ADCE.Sum <- summary(ADCE.Fit)) #Compare LL of the ADE, ACE, and ADCE models - are they different? ADE.Sum$Minus2LogLikelihood ADCE.Sum$Minus2LogLikelihood ADCE.Sum$Minus2LogLikelihood #Is the -2LL from YOUR ADCE model the same as your neighbors' -2LL?? #Compare our Mx estimates to the other ones #ALL estimates are equally likely, but the ADCE ones can be wacky, and depend on start values! mxEval(A,model=ADCE.Fit$MZ,T) mxEval(C,model=ADCE.Fit$MZ,T) mxEval(D,model=ADCE.Fit$MZ,T) #This is parameter indeterminacy in action!! # ----------------------------------------------------------------------- # end practical 3D # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # END PRACTICAL 3 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 4 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # 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*** ADE.fixC05.Model <- omxSetParameters(ADE.Fit, labels="VC", free=FALSE, values=.05) ADE.fixC05.Fit <- mxRun(ADE.fixC05.Model, intervals=FALSE) (ADE.fixC05.Sum <- summary(ADE.fixC05.Fit)) #Compare this model's -2LL to the other 3 - why are they all the same? ADE.fixC05.Sum$Minus2LogLikelihood ADCE.Sum$Minus2LogLikelihood ADE.Sum$Minus2LogLikelihood ACE.Sum$Minus2LogLikelihood #Get the 2 estimates and the fixed value of VC mxEval(A,ADE.fixC05.Fit$MZ,T) #estimated assuming VC=.05 mxEval(D,ADE.fixC05.Fit$MZ,T) #estimated assuming VC=.05 mxEval(C,ADE.fixC05.Fit$MZ,T) #fixed to .05 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # Fit ADE model where C is assumed to be ***.1*** ADE.fixC10.Model <- omxSetParameters(ADE.Fit, labels="VC", free=FALSE, values=.1) ADE.fixC10.Fit <- mxRun(ADE.fixC10.Model, intervals=FALSE) (ADE.fixC10.Sum <- summary(ADE.fixC10.Fit)) #Compare this model's -2LL to the others ADE.fixC10.Sum$Minus2LogLikelihood #Get the 2 estimates and the fixed value of VC mxEval(A,ADE.fixC10.Fit$MZ,T) #estimated assuming VC=.10 mxEval(D,ADE.fixC10.Fit$MZ,T) #estimated assuming VC=.10 mxEval(C,ADE.fixC10.Fit$MZ,T) #fixed to .10 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # Fit ADE model where C is assumed to be ***.15*** ADE.fixC15.Model <- omxSetParameters(ADE.Fit, labels="VC", free=FALSE, values=.15) ADE.fixC15.Fit <- mxRun(ADE.fixC15.Model, intervals=FALSE) (ADE.fixC15.Sum <- summary(ADE.fixC15.Fit)) #Compare this model's -2LL to the others ADE.fixC15.Sum$Minus2LogLikelihood #Get the 2 estimates and the fixed value of VC mxEval(A,ADE.fixC15.Fit$MZ,T) #estimated assuming VC=.15 mxEval(D,ADE.fixC15.Fit$MZ,T) #estimated assuming VC=.15 mxEval(C,ADE.fixC15.Fit$MZ,T) #fixed to .15 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # Fit ADE model where C is assumed to be ***.30*** ADE.fixC30.Model <- omxSetParameters(ADE.Fit, labels="VC", free=FALSE, values=.30) ADE.fixC30.Fit <- mxRun(ADE.fixC30.Model, intervals=FALSE) (ADE.fixC30.Sum <- summary(ADE.fixC30.Fit)) #Compare this model's -2LL to the others ADE.fixC30.Sum$Minus2LogLikelihood #Get the 2 estimates and the fixed value of VC mxEval(A,ADE.fixC30.Fit$MZ,T) #estimated assuming VC=.30 mxEval(D,ADE.fixC30.Fit$MZ,T) #estimated assuming VC=.30 mxEval(C,ADE.fixC30.Fit$MZ,T) #fixed to .30 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- #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(mxEval(A,ADE.Fit$MZ,T),mxEval(A,ADE.fixC05.Fit$MZ,T),mxEval(A,ADE.fixC10.Fit$MZ,T),mxEval(A,ADE.fixC15.Fit$MZ,T),mxEval(A,ADE.fixC30.Fit$MZ,T)), VD=c(mxEval(D,ADE.Fit$MZ,T),mxEval(D,ADE.fixC05.Fit$MZ,T),mxEval(D,ADE.fixC10.Fit$MZ,T),mxEval(D,ADE.fixC15.Fit$MZ,T),mxEval(D,ADE.fixC30.Fit$MZ,T)), VC=c(mxEval(C,ADE.Fit$MZ,T),mxEval(C,ADE.fixC05.Fit$MZ,T),mxEval(C,ADE.fixC10.Fit$MZ,T),mxEval(C,ADE.fixC15.Fit$MZ,T),mxEval(C,ADE.fixC30.Fit$MZ,T))) RESULTS3 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # END PRACTICAL 4 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # START PRACTICAL 5 - fitting Nuclear Twin Family Design # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # start practical 5A - Prep the data # ----------------------------------------------------------------------- #Get our dataset from the web. #This is a cool R trick if your file exist on a web server con <- "http://www.matthewckeller.com/public/ntf.data" ntf.data <- read.table(con,header=TRUE) # Look at 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) #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] # ----------------------------------------------------------------------- # end practical 5A # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # start practical 5B # ----------------------------------------------------------------------- #Fit NTF Model with **Sibling Env** (instead of "C" - both ~ the same) # Paths going to latent variables # Note that we now have to use path coefficient model approach for ETFDs 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 in the data dataMZ <- mxData( observed=mz.ntf, type="raw" ) dataDZ <- mxData( observed=dz.ntf, type="raw" ) #Objective objects for two groups objMZRels <- mxExpectationNormal(covariance="expCovMZRels",means="expMeanBoth", dimnames=c('tw1','tw2','father','mother')) objDZRels <- mxExpectationNormal(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,myfitfun) modelMZ <- mxModel("MZntf",params,dataMZ,objMZRels) modelDZ <- mxModel("DZntf",params,dataDZ,objDZRels) obj <- mxFitFunctionMultigroup(c("MZntf","DZntf")) ADCE.NTF.Model <- mxModel("ntfADCE",modelMZ,modelDZ,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 will not exactly equal true values due to sampling variability # which wasn't an issue above bc we forced the CVmz and CVdz be EXACTLY what we expected) # A (true simulated value was .30) mxEval(a%*%q1%*%a, ADCE.NTF.Fit$MZntf, T) # D (true simulated value was .30) mxEval(d%*%d, ADCE.NTF.Fit$MZntf, T) # S (true simulated value was .10) mxEval(s%*%s, ADCE.NTF.Fit$MZntf, T) # E (true simulated value was .30) mxEval(e%*%e, ADCE.NTF.Fit$MZntf, T) # ----------------------------------------------------------------------- # end practical 5B # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # END PRACTICAL 5 # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # -----------------------------------------------------------------------